An introductory guide to fluid models with anisotropic temperatures Part 1 -- CGL description and collisionless fluid hierarchy
P. Hunana, A. Tenerani, G. P. Zank, E. Khomenko, M. L. Goldstein, G., M. Webb, P. S. Cally, M. Collados, M. Velli, L. Adhikari
TL;DR
This paper provides a comprehensive introduction to collisionless fluid models with anisotropic temperatures, focusing on the CGL model and non-Landau closures, relevant for astrophysical plasma turbulence modeling.
Contribution
It offers a detailed guide to non-Landau fluid closures, including the CGL model, Hall and FLR effects, and the impact of heat flux vectors on plasma instabilities.
Findings
Non-gyrotropic heat flux significantly affects firehose instability growth rates.
Higher-order moment fluid models are generally unstable without Landau closures.
Normal closure at the 4th moment is the last stable non-Landau fluid model.
Abstract
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 modeling 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 CGL fluid model in great detail, together with fluid models that contain dispersive effects introduced by the Hall term and by the…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
An introductory guide to fluid models with anisotropic temperatures
Part 1 - CGL description and collisionless fluid hierarchy
P. Hunana11affiliation: Instituto de Astrofísica de Canarias (IAC), La Laguna, Tenerife, 38205, Spain; [email protected] 22affiliation: Universidad de La Laguna, La Laguna, Tenerife, 38206, Spain , A. Tenerani33affiliation: Department of Physics, The University of Texas at Austin, TX 78712, USA , G. P. Zank44affiliation: Center for Space Plasma and Aeronomic Research (CSPAR), University of Alabama, Huntsville, AL 35805, USA 55affiliation: Department of Space Science, University of Alabama, Huntsville, AL 35899, USA , E. Khomenko11affiliationmark: 22affiliationmark: , M. L. Goldstein66affiliation: Space Science Institute, Boulder, CO 80301, USA ,
G. M. Webb44affiliationmark: , P. S. Cally77affiliation: School of Mathematics, Monash University, Clayton, Victoria 3800, Australia , M. Collados11affiliationmark: 22affiliationmark: , M. Velli88affiliation: Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles, CA 90095, USA , L. Adhikari44affiliationmark:
Abstract
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 modeling 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 CGL fluid model in great detail, together with fluid models that contain dispersive effects introduced by the Hall term and by the finite Larmor radius (FLR) 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 4th-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 4th-order moment, Landau fluid closures are required.
Contents
-
4.3 Hall-CGL dispersion relation for arbitrary propagation angle
-
5.8 FLR corrections with non-gyrotropic heat flux vectors (FLR3)
-
11 Collisionless damping in fluid models - Landau fluid models
-
12.2 Impossibility to go beyond CGL2 without Landau fluid closures
1 Introduction
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 modeling, 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 (1958) to predict the existence of the solar wind, which was surprisingly several decades later than 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 analyzed.
The correct modeling and understanding of collisionless plasmas in a fluid framework, does not concern 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 non-linear 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., . In contrast, in collisionless systems with 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., . 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 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 and , even if the mean pressure values are equal and no mean temperature anisotropy exists.
The simplest collisionless fluid description is the CGL fluid model - named after Chew, Goldberger and Low Chew et al. (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 ), since the pressure fluctuations still remain anisotropic and the evolution equations for and 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 et al., 1954), see for example Snyder et al. (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 1) is derived from collisionless Vlasov equation with zero right hand side; and 2) that has two different pressure equations. Our definition therefore differs from an alternative view of for example Zank et al. (2014); Zank (2014), where models with scattering operator that reflects charged particle scattering by electromagnetic fluctuations on the right hand side of Vlasov equation are also viewed as collisionless.
Fluid modeling 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 to be 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 researchers new to the subject to follow the algebra easily. Instead of spending years, an interested reader should be able to get 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, Section 2, we start with the detailed derivation of the pressure tensor equation. The pressure tensor is then decomposed to its gyrotropic part (also referred to as ), 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 , see e.g. Braginskii (1958, 1965), the quantity is called the stress tensor. The decomposition procedure yields rigorously exact (even though still not closed) evolution equations for and , see e.g. Passot & Sulem (2004); Goswami et al. (2005); Oraevskii et al. (1968). Importantly, at the leading order (by neglecting the and also the non-gyrotropic heat flux contributions ), 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 used in that paper are related to the usual by relations and . 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 eq. (98)-(101). 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 process that can be encountered, and that are responsible for the breaking of the adiabatic invariants; see eq. (108), (109). 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. (2017b, a). That formulation is obtained by considering the total “internal” energy for each particle species (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 , , and formulate the total conservation of energy with the possibility of anisotropic plasma heating, see eq. (121)-(124).
The CGL description is analyzed in great detail in Section 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 . Collisionless plasmas can not 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 scale invariant, even though valid only on the largest scales. The instability thresholds present in the CGL model therefore can be referred to as “hard” thresholds, i.e. thresholds that are obtained in the long-wavelength low-frequency limit. The firehose and mirror instabilities are believed to play a crucial role in solar wind dynamics (see e.g. Hellinger et al. (2006); Bale et al. (2009)), and in Figure 1, they are plotted in the usual plane with logarithmic scales. While the firehose threshold matches the one obtained from kinetic theory, the CGL mirror threshold contains the well-known factor of 6 error for large values. The factor of 6 error can be interpreted as inadequacy of the adiabatic CGL closure in the very slow-dynamics context, such as the mirror instability. This is further addressed in Section 9.4, where we discuss that the “static” closure, which can be viewed as a generalization of isothermal closure in presence of temperature anisotropy and variations of magnetic field strength, reproduces the correct mirror threshold (Constantinescu, 2002; Chust & Belmont, 2006; Passot et al., 2006).
We discuss the core differences between MHD and CGL, which can be nicely summarized with the concept of adiabatic indices , related to the number of degrees of freedom by . While MHD is fully 3D with , the CGL can be viewed as composed of 1D and 2D dynamics with and . We also address the velocity and magnetic field eigenvectors. Similarly to the velocity field eigenvector in MHD (see Figure 2), which shows a “singular” behavior for strictly parallel propagation with , the CGL eigenvector (Figure 3) shows similar singularity for . We also briefly discuss fluid models with empirical “free” polytropic indices , studied for example by Hau & Sonnerup (1993); Hau et al. (1993); Abraham-Shrauner (1973).
In Section 4, we introduce the simplest dispersive effects by including the Hall term into the induction equation, and study dispersion relations of the Hall-CGL model. We focus on the parallel firehose instability, and show that the instability is indeed associated with the whistler mode, see Figures 4 and 5. We show that negative real frequencies have to be handled carefully, and that non-causal analytic solutions (438)-(441) have to be modified to the causal form (449)-(452). We briefly discuss the simplest ion-cyclotron resonances, and compare solutions of the Hall-CGL model with solutions of linear kinetic theory, see Figure 6.
In Section 5, we evaluate the FLR tensor at several levels of approximation. The evaluation of the FLR tensor is cumbersome, because the tensor is described by the pressure tensor equation implicitly. We first reproduce the fully nonlinear “inversion” procedure on how to obtain from expression that can be found for example in Hsu et al. (1986); Passot & Sulem (2004); Ramos (2005) as a brief note. By applying this inversion procedure at the pressure tensor equation evaluates the FLR corrections correctly along the magnetic field lines, but leads to an equation for that is still implicit. Nevertheless, by evaluating the resulting equation at the leading-order (technically first order in frequency and wavenumber), leads to an explicit expression for . We first recover the nonlinear result of Schekochihin et al. (2010), derived in that paper slightly differently, without using the inversion procedure. We point out that the result can be slightly simplified, and by rearranging the expression yields two different useful forms for writing the nonlinear . Finally, by using the non-dispersive (MHD) induction equation, we obtain the nonlinear result of Ramos (2005) (see also Macmahon (1965)). For further evaluation of nonlinear FLR corrections to higher-orders, an advanced reader is referred to Ramos (2005). We continue with the evaluation of the FLR tensor in the linear approximation, i.e. when the magnetic field lines are not too distorted. For the 1st-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 ), where is time between two collisions ( 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 Hunana & Zank (2017). 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 2nd-order FLR corrections (FLR2) are here defined as containing the Hall-term and the time-derivative . We also consider FLR corrections with the non-gyrotropic heat flux vectors , , 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 et al. (2017), see eq. (645)-(647). 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); Sulem & Passot (2015) and references therein.
In Section 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 Hunana & Zank (2017). 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); Rosin et al. (2011) and references therein.
In Section 7, we derive the heat flux tensor equation. The heat flux tensor is then decomposed into its gyrotropic and non-gyrotropic parts, . The procedure yields evolution equations for the gyrotropic parallel and perpendicular heat fluxes, and , that contain the tensor of the 4th-order moment . 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 4th-order moment as well, since there are several possible cancellations even at the linear level. The non-gyrotropic heat flux can be further decomposed to the non-gyrotropic heat flux vectors , and tensor . The detailed algebra of the non-gyrotropic heat flux vectors , , i.e. how to express them through lower-order moments, is presented in the 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 3rd-rank tensor , and an advanced reader is referred to Ramos (2005).
In Section 8, we consider the 4th-order moment which is a tensor of 4th-rank, . The moment is again decomposed to its gyrotropic and non-gyrotropic parts, and the gyrotropic part has three scalar components, , and . We show that for a bi-Maxwellian distribution function, the gyrotropic components can indeed be evaluated as , , and . 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 non-gyrotropic bi-Maxwellian through a combination of and , see e.g. Oraevskii et al. (1968).
In Section 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 2nd-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 et al. (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 Section 10, we consider the bi-Kappa distribution function. We show that the closure at the 4th-order moment is constructed by , , and , where the coefficient , and the closure is valid for . We call this closure and the associated fluid model as “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 , . We do not calculate the for the bi-Kappa distribution, and therefore we do not provide the second order non-gyrotropic heat flux vectors.
In Section 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 the 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 how the Vlasov equation is integrated over the velocity space. We introduce preliminary ideas how the Landau fluid closures will be constructed. For example, for closures performed at the 4th-order moment instead of the “normal” closure, one needs to consider perturbations around this state, and prescribe , , and . The deviations , , will be calculated from linear kinetic theory in Part 2, by performing Landau fluid closures.
In Section 12, we derive evolution equation for a general n-th order fluid moment (a tensor with components). We then consider fluid models in the simplified 1D geometry that can be viewed as an electrostatic case (or as a propagation along the mean magnetic field), that are closed at a general n-th order level by a Maxwellian fluid closure. A dispersion relation is obtained, which for 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 , 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.
2 Pressure tensor equation
Collisionless plasmas are described by the Vlasov equation, which in CGS units reads
[TABLE]
and which describes how a distribution function evolves in time. The is the index of species and for protons, for electrons, etc. The is the particle charge, the particle mass, the speed of light, the electric field vector and the magnetic field vector. The species index 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 . It is important to realize that the distribution function just describes the probability to find a particle with velocity at the position and that the “kinetic” velocity entering the distribution function is a completely independent variable from , i.e.
[TABLE]
Also, the magnetic and electric fields are macroscopic quantities that do not depend on and can be moved outside of velocity integrals over , or in another words and . The definitions of the fluid moments are
[TABLE]
where we have omitted the tensor product notation that is sometimes written down explicitly as and in the index notation . If two vectors or tensors are next to each other without an operator between them, a tensor product is always assumed. The number density is a scalar, the fluid velocity is a vector, the pressure tensor is a tensor of second rank ( matrix), the heat flux tensor is a tensor of third rank (with components), the tensor is of fourth rank (with components), etc. Directly from the definitions, it is obvious that all fluid tensors must be symmetric in all of their indices, i.e. , The fluctuating velocity is defined as
[TABLE]
and should not be confused with the speed of light .
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
[TABLE]
The is a general n-th order tensor, is a distribution function, the left hand side is a 3D integral calculated over the entire velocity volume and the right hand side is a surface integral calculated over a boundary of that volume , where 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 to in each velocity component and the integration on the right hand side of (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
[TABLE]
Since the area , the sufficient condition for the integral to vanish, can be estimated more precisely as
[TABLE]
When calculating the fluid hierarchy, the encountered expressions are always , where is a positive integer. Then, for example, for a Maxwellian distribution the limit (11) is always zero for all . Even for a slower converging distribution , the limit is always zero. For distribution functions proportional to a power law, for example as for a kappa-distribution, the situation is more restrictive, with some minimum required values of . When we will calculate the n-th order fluid moment (Section 12), we will see that two surface integrals are encountered, one with and one with . Therefore, for a kappa-distribution, the strict condition for yields the requirement . For example, many fluid models discussed in this guide are closed at the 4th-order moment , so , which implies the requirement . At first sight, the required limit (11) can be considered only a technical 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 cut-off 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 , and usually a cut-off is implicitly assumed. However, once we calculate the 4th-order moments for a kappa distribution, we will see that for example , where , so the restriction returns, and has to be applied regardless. When the limit (11) is satisfied, the neglect of the surface integrals (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 (9) that are zero asymptotically, by using and .
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 or (another possibility is ). 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
[TABLE]
We write down the derivatives with respect to time and velocity explicitly, but we abbreviate the derivative with respect to spatial coordinates as . We will need
[TABLE]
where in the last identity we used the result since the Levi-Civita tensor is antisymmetric and the Kronecker is a symmetric tensor. Integrating the first term of the Vlasov equation yields
[TABLE]
The second term of the Vlasov equation yields
[TABLE]
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 a 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 , the divergence means . This convention is the reason why in the above expression we used that and also the ordering of and inside reflects that in the index notation we have .
The third term of the Vlasov equation calculates
[TABLE]
Finally, the fourth term of the Vlasov equation calculates
[TABLE]
In the expression above, one encounters a vector product of a vector with a matrix, . This operator might appear unusual at first, but it is just a generalization of a vector product of two vectors, it is defined as and the result is a matrix. Similarly, a vector product between a vector and a tensor of 3rd rank is defined as and for a tensor of n-th rank (i.e., the second index in is the vector index and the third index in is the first index of tensor ). However, a definition of a vector product where a matrix is applied on a vector reads , and in general , but since the is here symmetric, . The vector product (cross product) is addressed further in the Appendix B. The term (24) is sometimes rewritten as
[TABLE]
Combining all the results together , yields the entire pressure tensor equation as
[TABLE]
One notices that, for example, the pressure tensor equation of Kulsrud (1983), eq. 57, contains a minus sign typo in the last term, which should be written as (25). That it is indeed a typo and not a definition of is evident from his subsequent eq. 59. It is useful to introduce a symmetric operator that acts on a matrix according to
[TABLE]
which yields a more compact form of the pressure tensor equation
[TABLE]
2.1 Pressure tensor decomposition
By introducing a unit vector in the direction of the local magnetic field and the gyrofrequency , the last term in the pressure tensor equation (26) can be rewritten as
[TABLE]
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 , where is the scalar pressure and is the unit matrix. In index notation , and we can verify that
[TABLE]
There is however a much more general solution that makes the term (29) equal to zero, which is
[TABLE]
where the “g” stands for gyrotropic. For quick analytic calculations it is sometimes easier to use since we have to deal with only once instead of twice. We need to verify that the term (29) indeed disappears,
[TABLE]
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
[TABLE]
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 and yields a scalar. Sometimes, the double contraction is defined as , which for symmetric matrices is equivalent to the previous definition. As a reminder, the usual matrix product of two matrices is and yields a matrix. It is useful to write down the following identities that involve the double contraction,
[TABLE]
Also interestingly, the trace of a matrix can be expressed through the double contraction with , i.e. . The pressure decomposition (35) can now be verified easily. Even more revealing is to apply the decomposition (35) directly at the definition of the entire pressure tensor (5), which yields
[TABLE]
where we have used the fact that the magnitude of the 3D velocity vector can be decomposed as , or equivalently .
A distribution function is called isotropic when the direction of the velocity vector does not matter, and the distribution function depends only on magnitude , as for example the Maxwellian distribution . A distribution function is called gyrotropic when the direction of the perpendicular velocity vector does not matter, and the function depends only on , as for example the bi-Maxwellian distribution . 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 , . For the pressure tensor, there are only two possibilities, represented by the parallel and perpendicular pressure (41), (42).
The part of the pressure tensor extracted by the decomposition (35) is therefore called the gyrotropic pressure , sometimes also called the CGL pressure . 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 (FLR) corrections to the (gyrotropic) pressure tensor, or sometimes non-gyrotropic or gyroviscous stress tensor. The entire pressure tensor is thus decomposed according to
[TABLE]
and for clarity we write down the decomposition explicitly in matrix notation
[TABLE]
In the momentum equation, the pressure tensor enters through its divergence, which is useful to break down to its components. Since , and , the divergence of the pressure tensor calculates
[TABLE]
Importantly, we have seen in (41), (42) that the scalar parallel and perpendicular pressures are extracted by applying the double contractions at the entire pressure tensor (and not necessarily only at the gyrotropic part). The decomposition (43) can be rewritten as
[TABLE]
By applying the double contraction to (55) yields the first requirement for the FLR tensor
[TABLE]
Similarly, by applying to (55) yields the second requirement for the FLR tensor,
[TABLE]
By further using the first requirement, the second requirement can be rewritten as
[TABLE]
By using the pressure decomposition (43) in the pressure tensor equation (26), and using the fact that the gyrotropic part of pressure satisfies , the pressure tensor equation can be rewritten as
[TABLE]
We can now derive the time dependent equations for parallel and perpendicular pressures by applying the usual double contractions on this pressure tensor equation.
2.2 Parallel pressure equation
We calculate the double contraction with term by term. We will need the following identities
[TABLE]
The first term calculates
[TABLE]
The second term calculates
[TABLE]
The third term calculates
[TABLE]
The fourth term calculates similarly
[TABLE]
and the final fifth term becomes
[TABLE]
The entire equation for the parallel pressure therefore reads
[TABLE]
Using the symmetric operator (27) and the definition for convective derivative , the equation is written in a simple form as
[TABLE]
which corresponds to eq. 26 of Passot & Sulem (2004), eq. 9 of Goswami et al. (2005) (see also Passot & Sulem (2007); Passot et al. (2012); Oraevskii et al. (1968)).
2.3 Perpendicular pressure equation
It is possible to either apply the double contraction on the pressure tensor equation, and directly obtain the equation for , or to apply the Trace (the double contraction ), obtain the equation for , and subtract the previously obtained equation for . Both approaches are equivalent, since . In another words, directly calculating the trace of the pressure tensor yields
[TABLE]
where we have used that . Applying the trace operator at the pressure tensor equation (59) term by term yields for the first term
[TABLE]
for the second term
[TABLE]
the third term
[TABLE]
the fourth term is equivalent to the third term
[TABLE]
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
[TABLE]
Combining all the terms, using the parallel pressure equation (68) and dividing by 2 yields the time dependent equation for the perpendicular pressure that reads
[TABLE]
The result corresponds to for example eq. 25 of Passot & Sulem (2004), eq. 8 of Goswami et al. (2005). It is important to note that the parallel 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 slightly so that we can directly use the induction equation and eliminate the time derivative. The term can be rewritten as
[TABLE]
We will need the following identity for the time derivative of the unit vector (see later in the text),
[TABLE]
which, when used in the above expression together with , allows us to write
[TABLE]
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 and we will show that if only the gyrotropic heat flux components are considered, the heat flux terms entering the pressure equations read
[TABLE]
Therefore, by splitting the heat flux into gyrotropic and non-gyrotropic parts, , the pressure equations read
[TABLE]
Terms with the non-gyrotropic heat flux can be further split to non-gyrotropic heat flux vectors , and heat flux tensor , which is addressed in the Appendix D, see equations (D120), (D122). 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 Laveder et al. (2013), 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, Goldberger, Low 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 , in the notation of that paper are equal to the usual pressures , , this is not the case for the heat fluxes. The gyrotropic heat flux tensor in that paper is decomposed into components , according to (Chew et al. 1956, eq. 34)
[TABLE]
In contrast, the usual gyrotropic heat flux decomposition to , reads (see later in the text)
[TABLE]
i.e. the last term is not present in their decomposition. However, this does not mean that their pressure equations (Chew et al. (1956), eq. 31, 32) are incorrect. It only implies that their equations are written in terms of
[TABLE]
By neglecting the FLR corrections and , the parallel and perpendicular pressure equations (83), (84) simplify to
[TABLE]
By using (87), the heat flux contributions in the parallel pressure equation above can be rewritten as
[TABLE]
which agrees with the parallel pressure equation obtained by Chew et al. (1956), eq. 31. The pressure equations of Chew et al. 1956 are therefore correct, but it has to be remembered that is not the usual parallel heat flux .111Sometimes the parallel heat flux is defined with the factor of 1/2 as , but this definition should be avoided. Actually, the authors do not even call the quantities heat flux components, but components of a “pressure-transport” tensor.
To conclude this subsection, we will go slightly ahead and use equations (176), (177) (without the zero right hand sides) that are the usual CGL pressure equations and valid only at long spatial scales when the Hall term in the induction equation is neglected. Under this assumption, the pressure equations (88), (89) are rewritten as
[TABLE]
It is often forgotten that the paper by Chew, Goldberger, Low 1956 derived the correct pressure equations with gyrotropic heat flux contributions. The authors however did not write the equations in the form (91), (92), and left them in the form (88), (89), using the “conservative” form only after prescribing the closure (their eq. 35, 36):
[TABLE]
Over time, equations (93), (94) became known as the CGL description. Equations (93), (94) are often replaced by the equations of state and , see e.g. Mjølhus (2009). Additionally, similarly to the definition of (appropriately normalized) entropy in MHD, , it is useful to define parallel and perpendicular entropy in the CGL model, according to
[TABLE]
The total entropy, , is equivalent to the MHD entropy when , see for example eq. 9-11 in Abraham-Shrauner (1967).
2.5 Physical meaning of CGL equations
The exact derivation of the CGL equation (93), (94) will be completed in the next section. Here we briefly discuss the physical meaning of these equations. The CGL equations (94), (93) can be interpreted as the conservation of the 1st and 2nd adiabatic invariants, as is nicely discussed for example by Kulsrud (1983) and by Gurnett & Bhattacharjee (2005). Considering many particles with velocity components and with respect to the mean magnetic field, the parallel and perpendicular pressures can be estimated as
[TABLE]
where the brackets represent an average over all particles. To complete the estimate, we need expressions for and , which come from the conservation of the adiabatic invariants. The 1st adiabatic invariant is the conservation of magnetic moment of a particle that is gyrating periodically around a mean magnetic field. The 2nd adiabatic invariant, sometimes also called the longitudinal invariant, is associated with a particle bouncing periodically between two magnetic mirror points. The conserved quantity is the integral over the parallel momentum of a particle , 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 labeled as L. The conservation of the 1st and 2nd adiabatic invariants is therefore
[TABLE]
Since we are considering only non-relativistic particles, the particle mass “m” is also a constant. The conservation of magnetic moment therefore implies , and using this in (96) yields , recovering (94). The CGL equation (94) therefore corresponds to conservation of magnetic moment, i.e., the 1st 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 , which yields an estimate for the area . The second principle is the conservation of the total mass of particles that are completely trapped in that flux tube and that can not escape, , yielding . The use of yields the final estimate for the non-intuitive length . The conservation of the 2nd adiabatic invariant therefore implies . Using this result in (96) implies , recovering (93). The CGL equation (93) 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 (93), (94) have to be modified with non-zero right hand sides.
In Part 2 of our guide (see section 4.1), we show that the conservation of magnetic moment is very useful for understanding the form of the perturbed distribution function in the gyrotropic limit. By performing calculations in the laboratory reference frame, we show that to obtain the correct in the gyrotropic limit actually requires the usual complicated kinetic integration around the unperturbed orbit (see Appendix of Part 2) and only then the gyrotropic limit can be imposed. In contrast, performing calculations in the guiding-center reference frame allows one to prescribe the conservation of magnetic moment from the beginning, and obtain the correct 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 et al., 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 at all. By reintroducing the species index r, the equations that were directly obtained by integrating the collisionless Vlasov equation read
[TABLE]
where the convective derivative . In the above equations, the pressure tensor for each species was just decomposed to a gyrotropic contributions , and the rest of the pressure tensor (the FLR pressure tensor ), according to
[TABLE]
which is a rigorous decomposition not introducing any simplifications. The heat flux contributions are here left at the most general level without any simplifications either. The equations are accompanied by the Maxwell’s equations, that for now we write down here with no simplifications to emphasize that no Maxwell’s equations were used in deriving the above system,
[TABLE]
Equations (98)-(104) represent 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 according to . If one does not like the symmetric operator in the expressions above, the symmetric operator is actually not necessary, since
[TABLE]
where we also got rid of the double contractions. Note that and also .
Equations (98)-(101) contain two very important catches - the entire heat flux tensor and the nongyrotropic pressure tensor 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 and the FLR pressure tensor . 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 (171), (172) (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
[TABLE]
Again, no simplifications were introduced and no Maxwell’s 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 (108)-(109) 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 that couples various species together through the Maxwell’s equations, and that also introduces the Hall-term responsible for the simplest dispersive effects. The heat flux tensor (a tensor of 3rd rank), is decomposed to its gyrotropic and non-gyrotropic (FLR) part . At large scales, the gyrotropic heat flux can be viewed as a gateway for the simplest forms of collisionless damping mechanisms, 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 & mirror force (Landau damping & 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 introduces further dispersive effects and therefore also modifies the collisionless damping rates, as well as the heat flux modifies the FLR pressure corrections . The quantities and are therefore generally coupled. Importantly, the FLR pressure tensor 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 nonlinear numerical simulations ( 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 and . The heat flux tensor 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 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 , simplest first order FLR corrections expanded around the proton scales 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 , but the first order frequency restriction is there nevertheless.
2.7 Conservation of energy
The equations (108)-(109), that correctly split the entire plasma heating to parallel and perpendicular heating, are indeed quite complicated. To gain further insight, let’s briefly consider a case in which we are not interested in parallel and perpendicular heating (whose net effect can be possibly 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. (2017b, a). By summing together the pressure equations (100), (101), we are interested in the evolution equation , or in another words we are interested in . Summing the equations together (actually the trace of the entire pressure tensor equation was calculated few pages back, that is how the equation was derived) yields
[TABLE]
Note that , and also . By integrating over the entire spatial volume of the plasma, one can define the internal energy (or thermal energy) for each particle species as
[TABLE]
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 (9), now in a spatial domain, and for periodic boundary conditions all divergence operators in (110) vanish, yielding
[TABLE]
Equations (110), (112) are equivalent to equations 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 , 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
[TABLE]
For the kinetic energy, by first applying at the momentum equation (99)
[TABLE]
which further implies
[TABLE]
and for the electromagnetic energy
[TABLE]
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)
[TABLE]
which beautifully clarifies how the energy can be transfered. 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 few terms in the pressure equations (100), (101) vanish. Considering anisotropic heating, and splitting the internal (thermal) energy into parallel and perpendicular parts
[TABLE]
the conservation of total energy can be expressed as
[TABLE]
where
[TABLE]
The anisotropic plasma heating is obviously a very complicated process.
3 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 modeling heliospheric and astrophysical plasmas, we recommend the reviews by Goldstein et al. (1995); Tu & Marsch (1995); Zank (1999); Zhou et al. (2004); Bruno & Carbone (2013). The CGL model is obtained by using the parallel and perpendicular pressure equations (83), (84) and prescribing zero heat flux , , and by neglecting the FLR pressure tensor . The pressure equations greatly simplify to the form
[TABLE]
By using the notation for the convective derivative , the pressure equations can be rewritten as
[TABLE]
Note that the CGL pressure equations are fully nonlinear and that the nonlinearities are actually of 4th order, i.e. CGL pressure equations are more nonlinear than the usual MHD pressure equation , where the nonlinearity is only of 2nd order. At first sight, expressions might appear to be a little unusual when written without any brackets. From MHD, one is familiar with nonlinear expressions of the type , which are also usually written without any brackets and can be interpreted in two equivalent forms, either as (a scalar applied on a vector ; or (a vector multiplied by a matrix ). In a similar fashion, the expression is meant to be interpreted as or ; but importantly, it is not meant to be interpreted as , because the derivative is meant to act only on in this case. For direct numerical simulations, the 4th order nonlinearity explicitly reads
[TABLE]
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 , where for protons, for electrons etc., the n-fluid CGL-type equations become
[TABLE]
The above CGL pressure equations are valid regardless of the form of the induction equation 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, 1) Gauss’s law where the total charge density , 2) , 3) Faraday’s induction equation and 4) Ampère’s law , where the total current . 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 and by neglecting the second term in the Ampère’s law 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 . It is important to emphasize that the Gauss’s law is replaced by the quasi-neutrality condition and no requirement is imposed on (other than it is small). For further discussion concerning the charge neutrality and , see for example Braginskii (1965) page 264, Webb et al. (2007) eq. 7, and Passot & Sulem (2007) eq. A1. The Maxwell’s equations in this approximation simplify to
[TABLE]
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 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 and , but models with more species can of course be considered. Charge neutrality implies that the proton and electron number densities are equal, i.e. . By using the definition of current , the electron velocities are related to the proton velocities by
[TABLE]
which after using (137) yields
[TABLE]
By using the momentum equation (132) for electrons, the electric field can be expressed as
[TABLE]
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 instead of increasing to infinity. When considering only relatively low frequencies, it is very useful to neglect the electron inertia term.222 The magnitude of terms on the r.h.s. of (140) can be easily compared to the advection term by estimating
(141)
where are characteristic length-scale and time-scale, the electron gyroradius and the electron thermal velocity. 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 (139), the electric field can be expressed as
[TABLE]
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
[TABLE]
and is accompanied by the parallel and perpendicular pressure equations (133), (134) for the proton and electron species. In the induction equation (145), 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 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 nonzero, since
[TABLE]
In contrast, by prescribing the electrons to be isotropic and isothermal (a typical assumption used in hybrid simulations) with pressure , , so that , eliminates the pressure term since . By prescribing the electrons to be isotropic but not isothermal, with pressure , the term initially does not disappear and it is equal to
[TABLE]
The term only disappears after the equation of state is prescribed, for example for . The term (147) 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 . The length is normalized with respect to (for now) arbitrary and the time with respect to . The normalized quantities are then
[TABLE]
Obviously . After dropping the tilde, the normalization yields the same density and pressure equations, whereas the normalized momentum and induction equation become
[TABLE]
where the Alfvén speed . The momentum equation implies that it is beneficial to choose the (so far unspecified) normalizing velocity to be the Alfvén speed . If the 2nd and 3rd 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 , similarly to the MHD system. If the Hall term is considered, it is beneficial to choose (so far unspecified) normalizing length-scale to be the ion inertial length
[TABLE]
Therefore by choosing , the normalizations that we use in all fluid models in this paper read
[TABLE]
The normalization also has the advantage that, when transferred to Fourier space, one obtains the dispersion relations for normalized wavenumber and frequency as
[TABLE]
It is useful to define the parallel and perpendicular thermal speeds
[TABLE]
and the definition of parallel and perpendicular temperatures are and . We use the usual convention with the Boltzmann constant .333Without this convention the temperature is defined as and the thermal speeds , so the plasma , expressions are the same as above. The parallel and perpendicular plasma beta are defined according to
[TABLE]
and later we will often use the abbreviation for the proton temperature anisotropy ratio
[TABLE]
The normalized momentum and induction equations therefore read (tilde are dropped)
[TABLE]
Another useful quantity that we will use later in the text is the proton Larmor radius (the gyroradius) , defined according to
[TABLE]
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 (NI) models of Zank et al. (2017); Zank & Matthaeus (1993).
3.2 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 (), 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’s work in physical units for a moment. The classical CGL model with cold electrons therefore reads (written in physical units)
[TABLE]
Similarly to the usual MHD model, the CGL equations are scale-invariant and do not contain any information about the physical length scale. A length scale is introduced by considering the Hall-term in the induction equation, and we discuss the Hall-CGL model in the next section. Also, the absence of the Hall-term allows the CGL pressure equations to be rewritten into “conservative” form as
[TABLE]
where is the convective derivative. This formulation is the most often cited form of the CGL pressure equations (often the is abbreviated only as , which can lead to confusion on how to correctly calculate the derivatives). We want to show the equivalence of (166) and (163), (164). We need the following identity,
[TABLE]
and the result is of course a scalar. Similarly for the spatial and convective derivative
[TABLE]
Later we will need
[TABLE]
Now we can directly calculate
[TABLE]
and, by using the density equation and the identity (168), we obtain
[TABLE]
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
[TABLE]
Now we use the induction equation
[TABLE]
that allows us to calculate
[TABLE]
and which, when used in (172)-(171) yields the final result
[TABLE]
Equations (176)-(177) verify that the CGL pressure equations (166) written in the “conservative” form and the directly derived CGL equations (163)-(164) are equivalent. We will see that the right hand side of the usual CGL equations (166) is nonzero 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 (161)-(165) are linearized with respect to mean values where the magnetic field is assumed to be in the z-direction. We assume that the mean value of the velocity is zero (the zero mean velocity value should not be confused with normalizing velocity as in (148), which is later chosen to be ). The terms are linearized as and the linearized system reads
[TABLE]
The divergence of the pressure tensor is
[TABLE]
and the first step of linearization yields, by components,
[TABLE]
which can be conveniently written in matrix notation as
[TABLE]
Notice that since , the derivatives of unit vectors are calculated according to
[TABLE]
and in the linear approximation one has
[TABLE]
which by components reads
[TABLE]
where, importantly, . The divergence of the pressure tensor in the liner approximation is therefore expressed as
[TABLE]
where the second term contributes if the temperature anisotropy is present, i.e. when . The vectors and are straightforward to calculate and are
[TABLE]
The entire set of linearized CGL equations reads
[TABLE]
Additionally, the 3 components of the induction equation are not independent, but satisfy the 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 are zero; alternatively, one can consider propagation in the y-z plane, where all expressions containing 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 (211)-(217) written in the x-z plane with , it is apparent that the equations for and decouple from the entire system and by applying at the equation for yields
[TABLE]
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 . For isotropic mean pressures the equation is equivalent to the MHD description with the usual Alfvén speed . Since
[TABLE]
where is the temperature ratio (157), the Alfvén wave equation can be rewritten as
[TABLE]
Considering a wave propagating obliquely in the x-z plane with the wavevector at an angle with respect to the z-direction (the direction), the parallel and perpendicular wavenumbers are defined as
[TABLE]
and the transformation to Fourier space is performed according to (see Appendix A)
[TABLE]
The wave equation (220) is then
[TABLE]
which yields the dispersion relation for the Alfvén waves
[TABLE]
It is important to emphasize that similarly to MHD, the Alfvén mode propagates in all the oblique directions, with frequency that is proportional to the projection of the wavenumber to the direction of , i.e. . If the expression under the square root in (224) becomes negative, i.e. if , 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 is . The firehose instability criterion can be rewritten in many forms, the most common being
[TABLE]
The first expression directly implies that if , no firehose instability can exist. The firehose instability therefore exists only for a relatively high plasma beta and only if the parallel temperature is higher than the perpendicular temperature. Importantly, the firehose instability criterium in the CGL model is equivalent to the one found from linear kinetic theory in the long-wavelength limit (for example, Rosenbluth (1956); Chandrasekhar et al. (1958); Parker (1958); Vedenov & Sagdeev (1958)). Alfvén wave propagation only affects components and which are decoupled from the rest of the system, and yields the other (eigenvector) components as being zero. The 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 in Fourier space, that measures the relative ratio of magnetic energy in the parallel component versus the total magnetic energy
[TABLE]
Since the Alfvén mode produces only fluctuations, with , the magnetic compressibility for this mode is
[TABLE]
for all the propagation directions. Sometimes the magnetic compressibility is defined as a ratio of parallel and perpendicular energies , 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 (226), 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
[TABLE]
and the Alfvén mode has for all propagation directions, since only the perpendicular velocity component is nonzero.
3.5 Slow and fast modes
By applying at the equations for and and using the evolution equations for one obtains
[TABLE]
The easiest way how to solve this coupled system is to transform to Fourier space, which yields
[TABLE]
The system has nontrivial solutions only if the determinant of the matrix is zero, yielding the CGL dispersion relation of 4th order in frequency that contains the slow and fast modes
[TABLE]
One can proceed in several ways with the above dispersion relation.
In the specific cases of strictly parallel () and strictly perpendicular () propagation directions it is actually easier to work directly with the system (237), because the off-diagonal components in the matrix are zero and the dynamics in the and components decouples. For strictly parallel propagation one solution is in the component, with the dispersion relation identical to the Alfvén wave (224). Since the magnetic and velocity fluctuations are only in the , components, the compressibility and for this mode. The second solution for parallel propagation is in the component, and it is called the sound mode, or the ion-acoustic mode, with dispersion relation . The parallel ion-acoustic mode has and . Considering strictly perpendicular propagation, the waves in the component vanish with the solution . The waves in the component are fast magnetosonic waves with dispersion relation , and and . It is also possible to define the parallel and perpendicular sound speeds ; , which yields the parallel acoustic mode dispersion and the perpendicular fast mode dispersion . 444Sometimes the parallel and perpendicular sounds speeds are defined as ; , which yields the parallel acoustic mode dispersion and the perpendicular fast mode dispersion
In the general case of oblique propagation, we choose to use the normalized frequencies and wavenumbers (154), and rewrite the dispersion relation as
[TABLE]
It is sometimes useful to introduce the parallel Alfvén phase speed , and its normalization , so that
[TABLE]
since can be used to write the CGL dispersion relation in a shorter form, as done later in (388). The solution of the polynomial (239) is simply , and it is obvious that the coefficient is always positive, since both the plasma beta and the temperature anisotropy ratio must be always positive. It is also possible to show that the discriminant under the square root (proof shown 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, becomes negative (and complex) when , implying that the slow mode becomes unstable when
[TABLE]
where we suppress the tildes on the wavenumbers since the above condition is valid for both and forms of wavenumbers. The above condition implies that there are two competing instabilities. The highly parallel propagation with yields instability threshold equivalent to (225) and the associated instability is called the parallel firehose instability. The highly oblique propagation with yields the mirror instability and the threshold in the CGL description reads 555To obtain the mirror instability threshold, it is not possible to consider strictly perpendicular propagation with because the coefficient would be exactly zero. The mirror threshold is revealed only in the highly oblique limit .
[TABLE]
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
[TABLE]
The threshold (242) implies that the necessary condition for the mirror instability to exist in the CGL description is , whereas in kinetic theory the necessary condition is . 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); 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 . By using (221) the wavenumber can be pulled out of the coefficients and the dispersion relation can be rewritten in terms of a phase speed as
[TABLE]
and (after staring for a sufficiently long time) it can be shown that
[TABLE]
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
[TABLE]
which is equivalent to dispersion relation of Abraham-Shrauner (1967), eq. 30, here written in the convenient units of and (for a direct comparison with units used in that paper, the relations are and ). The magnetic compressibility of the slow and fast modes can be calculated easily directly from the induction equations (215). The slow and fast modes generate only components and with the component . The transformation of (215) to Fourier space yields , and , and by applying the magnitude operator yields the magnetic compressibility
[TABLE]
For quasi-parallel propagation angles the magnetic compressibility is therefore not a good indicator for mode recognition, since all three modes have close to zero. The magnetic compressibility is however an excellent tool for oblique propagation angles, because the Alfvén mode has regardless of the propagation direction, whereas for the slow and fast mode increases with towards . We will leave discussion about the CGL velocity eigenvector (the ratio ) 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 and using this expression with the dispersion relation (244), the critical condition when the slow mode speed matches the Alfvén speed is , and since , this further yields . 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
[TABLE]
The expression is naturally satisfied for cases of strictly parallel and strictly perpendicular propagation. For , both the slow mode and the Alfvén mode have zero frequency. For , 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 (250) 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
[TABLE]
which agrees with the result obtained by Abraham-Shrauner 1967, eq. 42. The expression (251) reveals that there is a condition that can be satisfied for all directions of propagation and that reads
[TABLE]
This is a quadratic equation in which can be solved easily for a given (or solved for 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 (). In this case the inequality (252) becomes and the critical plasma beta at which the slow mode speed becomes larger than the Alfvén mode speed is
[TABLE]
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 (nonzero) propagation angles, the slow and fast mode expressions always get the correct ordering with . It is sometimes confusing about how the dispersion relations hold in the limit of parallel propagation which we will address here. In this limit, the dispersion relations (248) for the slow and fast modes can be directly evaluated as
[TABLE]
Now, because , for further calculation one needs to determine the sign of the expression under the square root. For , i.e. when , the above equations yield
[TABLE]
so the slow mode is the sound/acoustic mode and the fast mode coincides with the Alfvén mode. However, for , i.e. when , the dispersion relations yield
[TABLE]
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 , which is equivalent to
[TABLE]
all three modes propagate with the same phase speed . The criterion (257) can be satisfied only for and in the limit the asymptote is , so the possible range of is . This curve is plotted in Figure 1 as a green dotted line. The curve separates two regions. In the first region and . In the second region and . The interesting result to remember is that for isotropic temperatures () and parallel propagation, the critical plasma beta at which the slow mode speed matches the Alfvén mode speed is
[TABLE]
3.8 CGL with is not equivalent to MHD
It is important to emphasize that even for isotropic temperatures , the CGL dispersion (248) is not equivalent to the MHD dispersion for the slow and fast waves. The MHD dispersion is usually written in the form (see (C8) in the Appendix C)
[TABLE]
or alternatively
[TABLE]
which nicely shows that the expression under the square root is always greater or equal to zero. The MHD sound speed is defined as , where . Rewritten with parameters we use here , the MHD dispersion for the slow and fast waves reads
[TABLE]
or alternatively
[TABLE]
In the CGL description with isotropic temperatures the dispersion relation (248) for slow and fast waves reads
[TABLE]
or alternatively
[TABLE]
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 can not match the CGL dispersion relation simultaneously for all propagation directions . One could even play with an abstract idea to make the MHD dependent on 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 decouples from the dynamics in the other directions and the velocity eigenvector matrix for the slow and fast waves in MHD is
[TABLE]
or alternatively
[TABLE]
This yields the MHD dispersion relations for the slow and fast mode in the form of (259) or alternatively (260). It possible to either express the ratio or . We want to calculate the velocity ratio and we choose the first possibility that allows us to express that for the slow and fast modes (since )
[TABLE]
The ratio can be found for example from the second row of the matrix as
[TABLE]
or equivalently from the first row of the matrix as
[TABLE]
The possible range of the ratio is from to and the possible range for is from [math] to . Solutions can be easily divided into two categories - depending on the behavior for parallel propagation - and are plotted in Figure 2. For the sake of clarity, we only consider propagation with positive wavenumbers and , so the possible range of angle is from [math] to and , . Consider first the case . For parallel propagation, the dispersion relation (260) directly yields that for the fast mode and for the slow mode . Using the expression (280) then yields that for the fast mode and for the slow mode . This implies that for , all the curves start at
[TABLE]
regardless of the actual value of . The second case is straightforward, since for parallel propagation everything is just turned around; the fast mode dispersion relation is and the slow mode dispersion relation is . Implying that for , all the curves start at
[TABLE]
regardless of the actual value of . The perpendicular propagation limit is even easier since the split into two categories is not required and for both cases the dispersion (259) yields and , implying that for the fast mode and that for the slow mode . All the curves therefore “end” at
[TABLE]
whether the ratio or . The behavior of the curves is plotted in Figure 2.
In MHD, the ratio is sometimes used as the definition of plasma beta, which we will call here . As can be seen in Figure 2, for approaching 1, the curves display a “singular” behavior, and it is difficult to guess what will happen when reaches the value of 1 exactly. This special case is addressed in the next subsection. There are also two special cases, one for or and one for or . Both limits and are easily obtained by rewriting the dispersion relation (259) to a form
[TABLE]
In both limits the second term under the square root is small and the expression can be expanded as , yielding the slow and fast mode dispersion relations
[TABLE]
which further yields dispersion relations (287), (289). Note that the limit can be easily obtained directly from the velocity matrix (278) by neglecting the off-diagonal terms that decouples the two modes. In contrast, the limit is not that obvious and one needs to calculate the determinant.
In the limit , the MHD dispersion relation simplify to
[TABLE]
The expression (280) yields for the fast mode and the expression (281) yields for the slow mode , implying
[TABLE]
In the opposite limit , the MHD dispersions simplify to
[TABLE]
In this limit, the expression (280) yields , that describes the behavior of the fast mode for large beta values and the complete limit is . For the slow mode the expression (280) yields and the complete limit is . The solutions for large beta values therefore converge to
[TABLE]
and these two curves are plotted in Figure 2 as dashed lines.
3.9.1 MHD special case
Consider the very peculiar case of parallel propagation with . The velocity eigenvector matrix (278), when evaluated exactly at these values reads
[TABLE]
It is not possible to determine the ratio and the eigenvector must be found by taking the limit . For the dispersion relation for the slow and fast modes (260) simplifies to . As before, for the sake of simplicity we consider only positive wavenumbers, so the propagation angle is in the range and , . Evaluating the velocity matrix (278) at only, it is possible to pull the out of the system so that
[TABLE]
which nicely demonstrates that for the system approaches (297). Nevertheless, we can now calculate the velocity ratio (for example from the second row) as
[TABLE]
which can be evaluated at . For the fast mode and for the slow mode . This implies that for the peculiar case of , the velocity ratio for both the fast and slow modes. The expression (305) can of course be directly obtained by using the result in the expression (280) and by canceling the in the nominator and denominator. The expression (305) is valid for any angle . However, for the slow mode it becomes singular at , and has to be rearranged by multiplying and dividing with (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
[TABLE]
which yields in this peculiar case of , that the ratio for the slow and fast modes can be expressed as
[TABLE]
These solutions are plotted in Figure 2 as dotted lines.
3.10 CGL velocity eigenvector
To find the ratio for velocity components in the CGL description, it is useful for a moment to work with abbreviated symbols, and we write the matrix (237) for the and components as
[TABLE]
where the coefficients are
[TABLE]
The determinant must be zero which implies , and has a solution
[TABLE]
The last step was actually used in obtaining (247) which is otherwise quite difficult to see. Using this result in the above matrix, the eigenvector matrix becomes
[TABLE]
where
[TABLE]
The ratio of can be expressed (for example from the second row) as
[TABLE]
or alternatively (from the first row) as
[TABLE]
The behavior of these solutions is plotted in Figure 3, where we plot solutions for isotropic temperatures . Similarly to MHD, two classes of solutions can be characterized based upon the behavior for parallel propagation , depending if the quantity (324) evaluated at is positive or negative. The quantity was discussed in the section for parallel propagation and this quantity represents the green dotted curve in Figure 1 that separates the system into two regions. For , i.e. in the region , for the slow mode (use eq. 326) the ratio and for the fast mode (use eq. 325) the ratio . In this region, all the curves start at
[TABLE]
For , i.e. in the region , for the slow mode (use eq. 325) the ratio and for the fast mode (use eq. 326) the ratio . In this region, all the curves start at
[TABLE]
For perpendicular propagation the quantity is always positive and the separation into two regions is not necessary. For the slow mode and for the fast mode . All the curves, regardless of the value of therefore end at
[TABLE]
This summarizes the general behavior of the solutions. However, as can be seen from Figure 3, there is a special case when the parameters approach , i.e. when , that will be discussed in the next subsection. One can again consider the two limiting cases when is small or large. The case yields dispersion relations 666One can directly work with the dispersion relation (316) and for and bounded and less than an extremely large values, the parameter is positive and , on the other hand the parameters and , so an expansion of the square root is possible as
\frac{1}{2}(a+c)\pm\frac{1}{2}(a-c)\sqrt{1+\frac{4b^{2}}{(a-c)^{2}}}\simeq\frac{1}{2}(a+c)\pm\frac{1}{2}(a-c)\Big{(}1+\frac{2b^{2}}{(a-c)^{2}}\Big{)}\simeq\frac{1}{2}(a+c)\pm\frac{1}{2}\Big{[}(a-c)+\frac{2b^{2}}{a}\Big{]}.
(330)
The fast mode , and the slow mode , which yields low dispersion relations. With “less confidence” the same result can be directly obtained from the velocity matrix (314) by neglecting the off-diagonal terms.
[TABLE]
Using expression (326) for the slow mode yields and expression (325) for the fast mode yields , implying that in low beta limit the solutions are
[TABLE]
and the constant solutions approaching this limit are visible in the plot of Figure 3. The opposite limit in the CGL description is not particularly revealing since there is no obvious way to simplify the expansion, even for isotropic temperatures. For , the expansion yields
[TABLE]
where we neglected a term proportional to . The expression results in the correct parallel limit , and the correct perpendicular limit , . If the perpendicular fast mode limit is sufficient to be (which is acceptable since ), only the terms proportional to can be retained and the fast mode dispersion relation can be simplified to
[TABLE]
The expression (325) yields the velocity ratio of the fast mode as
[TABLE]
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 nonzero everywhere except exactly at (and also ). For , the entire contribution of comes from terms in (335) that are not proportional to . Nevertheless, the analytic velocity ratio (337) for the fast mode is sufficiently simple and to obtain analytic 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 . Our final large beta results therefore are
[TABLE]
Solutions are plotted in Figure 3 as dashed lines and fit the solutions of the full dispersion relations with very nicely.
3.10.1 CGL special case
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 . This CGL special case can be encountered only for because the temperature anisotropy must be . We start with the velocity eigenvector matrix (323) and for the special case considered , the expression (324) simplifies to . 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 can be pulled out of the matrix or alternatively one can use (325), (326) and calculate the limit. For the special case considered here, the parameter . Considering again for simplicity only positive wavenumbers, the expression (325) yields in this particular case
[TABLE]
The expression can be evaluated at , which yields for the slow mode and for the fast mode , implying that both modes start at
[TABLE]
For the perpendicular propagation the discussion is unchanged from the previous general case and yields and .
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
[TABLE]
One can introduce the concept of polytropic indices and to the CGL pressure equations (216), (217) and write
[TABLE]
where for the CGL model and . 777One can then define the parallel and perpendicular sound speeds as and , since in this model the parallel acoustic mode will have a dispersion and the perpendicular fast mode , as can be easily verified. 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 and on the other side. The important observation is that in the parallel pressure equation the acts only on the parallel velocity component , and in the perpendicular pressure equation the acts only on the perpendicular velocity components . In MHD, the acts on all three components . The values of polytropic indices are directly related to the number of degrees of freedom , through the well known relation . In the MHD description the number of degrees of freedom is and so . 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 and and for the perpendicular direction and . The subtle appearing, but nevertheless fundamental difference between (342) and (343), (344) 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 and equal to the mean MHD pressure . 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 (343), (344). If these pressure equations are “un-linearized”, i.e. if one performs a “reverse engineering” procedure to obtain a nonlinear equations (a procedure that is possible to do because of the known guidance by the CGL model), the nonlinear equations read
[TABLE]
which can be rewritten as
[TABLE]
Empirical fluid models of this type were studied for example by Hau & Sonnerup (1993) and Hau et al. (1993) with the motivation that spacecraft observations in the Earth’s magnetosheath often show behavior that cannot be explained by the CGL values of or by the isothermal values of , 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 (347) can be rewritten to a form
[TABLE]
or to an interesting alternative form
[TABLE]
Obviously, for CGL values of 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 (91), (92) which contain the gyrotropic heat flux contributions. Therefore, the deviations from the CGL values of 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 Hau et al. (1993) and also in Abraham-Shrauner (1973), who studied even more general empirical model with four polytropic indices (essentially, the , acting on and in (347) can be further split to two independent parameters). Interestingly, complex polytropic indices have been suggested to model Landau damping (Belmont & Mazelle, 1992). By considering the linearized CGL equations (211)-(215) and linearized polytropic pressure equations (343)-(344) 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
[TABLE]
and the transformation to Fourier space yields the velocity matrix
[TABLE]
or alternatively
[TABLE]
Calculating the determinant yields a dispersion relation for the slow and fast modes in the form
[TABLE]
Solutions for parallel propagation are the usual CGL Alfvén mode and the ion-acoustic mode , and for perpendicular propagation the fast mode . Similarly to the CGL description it can be shown that because of the trick (316), implying that the slow mode becomes unstable when , i.e. when
[TABLE]
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
[TABLE]
which is consistent with the result of Hau & Sonnerup (1993), their eq. 14, written in the form 888Note that the equation in the cited paper has opposite sign, which is not a typo, authors just call the mirror threshold a condition when a mode is stable - opposite to a more usual condition when a mode becomes unstable. Similarly, authors write the condition for the firehose stable region as .
[TABLE]
For , the condition (370) is naturally equivalent to the CGL mirror threshold. However, the condition offers a temptingly simple result to keep the perpendicular polytropic index unchanged and modify the parallel polytropic index to . In this way, the mirror threshold matches the correct kinetic threshold (243). 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 , 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 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 and focus on fluid models that derive the correct heat flux contributions and .
3.12 CGL and radially expanding solar wind
Similarly to MHD, the CGL fluid model can be very useful for understanding evolutions of temperatures in the radially expanding solar wind. A good discussion can be found for example in Matteini et al. (2007, 2013) and references therein. Let’s first consider the MHD case. Combining the pressure equation with the density equation yields
[TABLE]
which in a steady-state (or co-moving with the convective derivative) implies . By considering simple radially expanding solar wind, one prescribes that the density evolves with the heliocentric distance as , and by using , the ideal MHD model yields (in the absence of any dissipation and associated turbulent heating)
[TABLE]
In contrast, the CGL equations (166) are rewritten as
[TABLE]
which in a steady-state implies
[TABLE]
For the magnetic field it is possible to consider profiles of roughly , or (one can be more precise and consider the Parker spiral profile). The first magnetic field profile yields
[TABLE]
The second magnetic field profile yields
[TABLE]
There is also a curious case when the magnetic field profile is prescribed to be , i.e. a profile between two cases of and , that yields
[TABLE]
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 (376) will lead the system to a firehose threshold, and evolution according to (377) will lead the system to a mirror threshold. In general, if we fix the density profile to , any magnetic field profile steeper than will evolve the system towards a firehose threshold, and a magnetic field profile shallower than will evolve the system towards a mirror threshold. The estimations do not reveal 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 ). To find an instability, one of course needs to consider evolution equations for fluctuating variables and analyze associated dispersion relations, since it is the fluctuations that become unstable. Further discussion can be found in Hunana & Zank (2017).
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 et al. (2017) has shown that the firehose instability criterion reads (see eq. 9 in that paper)
[TABLE]
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 fluctuations starts to grow, the firehose criterion is modified according to (379). Therefore, the amplitude of fluctuations can not grow without bounds, and beyond some critical amplitude, the system becomes again stabilized. The condition (379) 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); Tenerani & Velli (2018).
To finish this subsection, the heuristically generalized CGL model with polytropic indices , yields that in the steady state the mean temperatures evolve according to
[TABLE]
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 and substituting one equation to 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
[TABLE]
where we have introduced the normalized parallel Alfvén speed (tilde on is dropped too). The determinant can be easily calculated and factorized in Maple, yielding the CGL dispersion relation in normalized units (tilde are dropped)
[TABLE]
which of course agrees with the previously obtained CGL solutions for the Alfvén mode (224), and the slow & fast modes (239).
4 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 thoery 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 (140). 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
[TABLE]
The notation is perhaps slightly misleading since only the first term in (390) is usually called the Hall electric field. Note that the first term in still has to be used in the proton momentum equation even for the simplest CGL (or MHD) models. The induction equation then reads
[TABLE]
which is then rewritten with the use of convective derivative as
[TABLE]
The last expression is then used in the completely general equations (171), (172), which yields
[TABLE]
By canceling pressure equations (128), (129), that are also the pressure equations of the Hall-CGL model, one obtains that for a general Hall-CGL model the electric field modifies the conservation laws according to
[TABLE]
implying the Hall term breaks the 1st and 2nd adiabatic invariants.
4.1 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 . The full model is described by the usual CGL equations (161)-(164), but the induction equation now reads
[TABLE]
The induction equation is linearized as
[TABLE]
By using (210), direct calculation of the second term yields
[TABLE]
and by using , the Hall-term contributions are
[TABLE]
The entire linearized induction equation reads
[TABLE]
which when written in the x-z plane (all ) yields
[TABLE]
These equations replace the 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 , 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’s quickly consider the two special cases of parallel and perpendicular propagation that can be solved analytically.
4.2 Parallel propagation
For propagation parallel to the mean magnetic field, all and the familiar sound mode /ion acoustic mode decouples in the , components as
[TABLE]
with the familiar dispersion relation . The dispersion relation of this mode does not depend on the Hall-term. However, the other two modes are now coupled through
[TABLE]
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
[TABLE]
Applying at the B-field equations yields
[TABLE]
which when transformed to Fourier space reads
[TABLE]
and the zero determinant requirement implies
[TABLE]
The solutions of the dispersion relation (quadratic polynomial in ) is simply
[TABLE]
where importantly, the was used. The solution (434) is the whistler mode, and the solution (435) is the ion-cyclotron mode. For isotropic temperatures (), the result is equal to the Hall-MHD dispersion relation. At long wavelengths () 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 , 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 , the limits are 999Note that the limit for the ion-cyclotron mode has to be calculated with a trick since both go to .
[TABLE]
For isotropic temperatures the frequency of the ion-cyclotron mode (that should be actually called the proton-cyclotron mode) converges towards the proton cyclotron frequency . For anisotropic temperatures, the frequency of the mode starts to be dependent and will converge to frequencies that can be higher or lower than the . 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 (434), (435) 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 (433) can be decomposed as 101010One can double check the result by working with quantities , yielding \Big{(}\omega^{2}\pm\frac{k_{\parallel}^{2}V_{A}^{2}}{\Omega_{p}}\omega-k_{\parallel}^{2}v_{A\parallel}^{2}\Big{)}B^{\mp}=0.
[TABLE]
The first bracket yields two solutions
[TABLE]
and the second bracket yields two solutions
[TABLE]
It is straightforward to check that yields the whistler mode solution (434), and yields the ion-cyclotron mode solution (435). Often, solutions (438)-(441) are written in an alternative form where the absolute value on 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
[TABLE]
with another two solutions obtained by substituting with .
4.2.1 Firehose instability
Now we can discuss the firehose instability. If the quantity , then all the solutions (438)-(441) have frequencies that are purely real for the entire range of wavenumbers. For isotropic temperatures , 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 (438)-(441) with are plotted in Figure 4 left panel, 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 , the solutions are naturally dependent. To show solutions in a fire-hose unstable regime, one has to choose since the instability does not exist for lower values of . In Figure 4, right panel, we have chosen a value of and varied the temperature anisotropy from down to a critical value of . In this regime, all the modes have purely real frequencies.
Importantly, if the temperature anisotropy parameter is further decreased, the frequencies of all the modes will become complex numbers (with a real and imaginary part) in the region where
[TABLE]
which can be viewed as a modified Hall-CGL firehose instability threshold that is now length-scale dependent. At long spatial scales (), the above criterium is naturally equivalent to the CGL firehose criterium. The term proportional to in (443) is always positive and as increases, after some critical wavenumber the frequency 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 and the temperature anisotropy is varied from to the lowest possible value of . The left panel shows real frequencies and the right panel 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 (443) is satisfied, for some range of the second terms in solutions (438)-(441) become purely imaginary and the first terms stay purely real. Consider one of the four possible quadrants in the Figure 5 left, for example the quadrant with . Up to some critical where the Hall-CGL firehose criterium is satisfied, both the blue mode and the cyan mode propagate with the same real frequency . The differences are in the imaginary frequency, see Figure 5 right, where in this quadrant the blue mode is unstable and the cyan mode is stable (and damped). After a critical (for example for it is ), 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 . We can therefore conclude that for positive real frequencies, , 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 , the whistler mode is stable and the ion-cyclotron mode is unstable. This is indeed contradictory. Keeping the , a change from positive to negative represents a change of direction of propagation along . 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 and , and the ion-cyclotron mode is stable. It is important to emphasize that from a perspective of solutions (438)-(441), 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 Section 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 can be possibly both or ). One should consider all the possible solutions and verify, which solutions are physically plausible. Therefore, after the firehose threshold is crossed, the analytic solutions (438)-(441) 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.
4.2.2 Generalization to causal (correct) solutions
The mind boggling puzzle, why the Hall-CGL solutions (438)-(441) 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 (419)-(420), in the following form
[TABLE]
where is a dimensionless real small positive number, for example . The dissipation will be later completely removed by the limit . The is there so that the parameter is indeed dimensionless (since here we usually normalize with respect to and , but other normalizations can be of course considered). The factor of 2 is there for pure convenience, so that the final dispersion relations (449)-(452) are of the utmost beauty, and that we don’t have to later redefine . The momentum equations are accompanied by the induction equations (421)-(422). The reader is encouraged to verify the calculations by switching to normalized variables , , , we continue in physical units.
By adding the dissipation, the dispersion relation (437) is modified to the following form
[TABLE]
It is useful to define the Hall-CGL firehose threshold
[TABLE]
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 (438)-(441) are “corrected” by the tiny dissipation to the following causal form
[TABLE]
The reader is encouraged to plot these solutions, for example with , and verify that the right panel of Figure 5, is indeed corrected. I.e., that the whistler is unstable for both and . We want to further address the limits. For clarity, let’s first consider the following simplest case
[TABLE]
where one can further suppress the small real on the right hand side. The limit (453) provides an answer, when the is equal to or . The result (453) can be easily generalized, and limits in expressions (449)-(452) are calculated according to
[TABLE]
Obviously, the whistler solutions (450), (451) are unstable, and the ion-cyclotron solutions (449), (452) 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 (449)-(452) are recovered with . I.e., by adding a 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 ), or more precisely through the current , see Part 2.
4.2.3 Polarization
The polarization of the magnetic field can be easily checked, actually without plugging any dispersion relations into the matrix (432), since the matrix can be represented simply as
[TABLE]
where and . It is useful to define the polarization as the angle between the complex and components, which is calculated as . 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 (461) is , so . The case (which is directly equal to the first bracket in (437) with solutions (438), (439)) yields that , so that
[TABLE]
The second case (which is equal to the second bracket in (437) with solutions (440), (441)) yields that , so that
[TABLE]
Similarly, one can calculate the polarization of electric field , and (when not specified that the field is considered), the words left/right polarized are usually meant for the field. The ratio can be easily calculated, for example from equations (486), (487) (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
[TABLE]
When plotting many solutions, it is useful to additionally divide the angle by , and define the and field polarizations as
[TABLE]
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 kinetic description; and also by Camporeale & Burgess (2017) who compared various linear modes in hybrid, gyrokinetic and fully kinetic descriptions. 111111Similar definition was used by Sahraoui et al. (2012), who however plot a bit confusing , since in the WHAMP code the electric field is normalized so that . When the wave is called left polarized and when , the wave is called right polarized. For the first group (462) this yields and (left polarized wave), and for the second group (463) and (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 (465) is redefined to
[TABLE]
A similar definition of polarization with 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 (440, blue line) the real frequency remains , and for the whistler mode (439, magenta line) . The situation is more confusing for the ion-cyclotron mode, where the firehose-stable and firehose-unstable regimes have to be considered separately. 1) In the firehose-stable regime, both ion-cyclotron modes are left polarized, since for the mode (438, green line) , and for the mode (441, cyan line) . 2) In contrast, in the firehose-unstable regime, both ion-cyclotron modes are right polarized, since for the mode (438, green line) , and for the mode (441, cyan line) .
4.2.4 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 from the beginning and expresses the fluid model through a matrix multiplied by a vector , leads to dispersion relations that are not singular explicitly. In contrast, expressing the model through a matrix multiplied by a vector , yields equivalent dispersion relations where the resonances are shown explicitly. The dispersion equation (437) can be easily rewritten as
[TABLE]
which, in the cold plasma limit () or in this case also for isotropic temperatures , simplifies to
[TABLE]
making it more evident that we have the simplest ion-cyclotron resonances here too. Indeed, for the frequency (first case) and (second case). For a plasma with finite and nonisotropic temperatures, the limit yields more general ion-cyclotron resonances (436).
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 and the electrons are prescribed to be cold with (in the WHAMP code, the value is chosen to be ). In the WHAMP code it is necessary to prescribe the ratio of the parallel proton thermal speed to the speed of light, and we choose . The solid lines are kinetic solutions with different . The blue curve is the whistler mode and the green curve is the ion-cyclotron mode. Solutions of the Hall-CGL model are shown in the left panel of Figure 6 and, for isotropic temperatures the solutions are completely independent and are represented with two dashed curves. It is shown that the kinetic solutions for the whistler mode are almost independent and only the mode with deviates from the other whistler solutions, and only between . For comparison, in the right panel of Figure 6 we plot solutions of the Hall-CGL-FLR2 model (the model is discussed later in the text).
4.3 Hall-CGL dispersion relation for arbitrary propagation angle
To make everything as clear as possible, let’s write down one more time the normalized equations we are solving, including the linearized and Fourier transformed equations in the x-z plane. The nonlinear Hall-CGL equations with cold electrons are (dropping the tilde)
[TABLE]
Linearizing in normalized units is even easier, because (temporarily introducing back the tilde for clarity) the ; ; ; . Linearizing these equations, writing them in the x-z plane, and Fourier transforming yields (dropping also the proton index ’p’)
[TABLE]
where the normalized . Alternatively, it is sometimes illuminating to work with the general induction equation where the electric field is used
[TABLE]
and in this case with cold electrons the electric field components are
[TABLE]
The induction equation (483)-(485) with the above electric field components is of course equivalent to the induction equation (478)-(480). All equations are normalized, the tilde are dropped, and the electric field , where . The Hall-CGL dispersion relation reads
[TABLE]
where the wavenumber . Perhaps the nicest form is to keep the CGL contributions on the l.h.s. (so the terms can be factorized to form the CGL dispersion relation) and put the Hall term contributions on the r.h.s., which yields the Hall-CGL dispersion relation
[TABLE]
Solutions for the parallel propagation () are of course the whistler and the ion-cyclotron mode
[TABLE]
with other two solutions obtained by exchanging with . By bringing back the tilde to the normalized solution (491), the result (442) in physical units is easily recovered. The other parallel solution is the ion-acoustic mode . For perpendicular propagation (), the fast mode solution is , 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” (Hunana & Zank, 2017), 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.
5 FLR corrections to the pressure tensor
The finite Larmor radius (FLR) corrections to the pressure tensor can be derived in several levels of approximation. The pressure tensor equation (59) can be rewritten as
[TABLE]
This equation is exact and since , the FLR tensor is described implicitly.
5.1 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); Ramos (2005) as a brief note. By applying to the left hand side, the first term calculates
[TABLE]
where we have used that . The second term on the left hand side of (492) calculates quite differently and we will need an identity
[TABLE]
Let’s also remind ourselves that in the index notation and . The second term calculates
[TABLE]
The entire left hand side therefore reads
[TABLE]
By performing a usual matrix (single) dot product of this result with matrix yields
[TABLE]
This further yields that
[TABLE]
and adding this result together with its transpose implies
[TABLE]
where we have used that the . One therefore derives that the FLR tensor can be extracted from the left hand side of equation (492) by performing an “inversion” procedure
[TABLE]
We can introduce matrix that will represent the right hand side of equation (492) as
[TABLE]
By performing operations (500) at this equation yields the FLR tensor in the following form
[TABLE]
The expression (503) is completely general since no simplifications were introduced, the expression is exact. However, the equation is still implicit, since on the right hand side the tensor contains the full pressure tensor . Nevertheless, the equation is extremely useful, once the expansion of the r.h.s. is performed, as we will do below. Note that the brackets are not unique, since
[TABLE]
where the last choice leaves the brackets unspecified (since both options are allowed), similarly to . Alternatively, one can get rid off the operator (which acting on a matrix yields the same matrix), and split the result into two parts
[TABLE]
Another possible form can be obtained by subtracting the evolution equations for the gyrotropic pressure components and . The scalar pressure equations were derived by applying and on the pressure tensor equation. The subtraction can be formally represented by introducing “overbar” projection operator that projects any matrix to
[TABLE]
The definition is of course motivated by the pressure decomposition (55). By applying this operator on (501), the l.h.s. remains the same, yielding . The solution (503) for the FLR tensor therefore rewrites
[TABLE]
The subtraction of scalar pressure equations is motivated by the observation, that when working in the linear approximation directly with equation (501), i.e. without performing the inversion procedure, the scalar pressure equation 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 (503) instead of (507).
Therefore, let’s work with (503). The leading-order nonlinear FLR corrections are obtained by neglecting the non-gyrotropic contributions and on the right hand side of (503), i.e. by prescribing
[TABLE]
Therefore, the only non-gyrotropic quantities that were retained on the r.h.s. of (503), are the perpendicular components of the velocity . The slightly unappealing factors in all the expressions can be removed, by redefining the cyclotron frequency with instead of , so that . Let’s calculate (508) step by step. For the brevity of calculations, it is useful to introduce convective derivative and work with
[TABLE]
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’s first write each term in in the index notation. By using
[TABLE]
in the following form
[TABLE]
and grouping similar terms in together yields
[TABLE]
The first line in the above expression is left unchanged, since it is the final result, defined below as a matrix , that can be written as \bm{W}=p_{\perp}\nabla\bm{u}+\nabla(q_{\perp}\bm{\hat{b}})+\big{(}p_{\perp}\nabla\bm{u}+\nabla(q_{\perp}\bm{\hat{b}})\big{)}^{T}=\big{(}p_{\perp}\nabla\bm{u}+\nabla(q_{\perp}\bm{\hat{b}})\big{)}^{S}. However, calculating the rest of the expression according to (508) simplifies several terms, and for example all terms in the third line of (516) vanish. The terms in the 3rd line containing vanish after applying , since . The terms in the 3rd line containing vanish after applying the symmetric operator since
[TABLE]
We therefore need to concentrate only on the 2nd line of (516), and slightly rewrite that one. For example, the first term in the 2nd line of (516), the term , calculates
[TABLE]
The second term in the 2nd line of (516), the term , vanishes after applying , since
[TABLE]
The third term in the 2nd line of (516), the term however does not vanish. Note that b_{j}\bm{\hat{b}}\cdot\nabla u_{i}=\big{(}(\bm{\hat{b}}\cdot\nabla\bm{u})\bm{\hat{b}}\big{)}_{ij} and the third term calculates
[TABLE]
The fourth term in the 2nd line of (516), calculates in the same fashion than the second and third term since , yielding
[TABLE]
Collecting all the results together yields the FLR tensor
[TABLE]
where is a matrix and is a vector. The result is equivalent to equations (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 and vector , the symbols and are used, where and . We wanted to avoid the use of and , since these symbols are typically used when decomposing the heat flux tensor. For a reader who just jumped to this result, the symmetric operator is defined as . For clarity, for any vector and any matrix , the following identities hold
[TABLE]
The last identity follows from an identity for general matrix , that reads , and for a matrix the , implying .
The result (522)-(524) can be slightly simplified. The matrix in the index notation reads
[TABLE]
and by applying eliminates only the last term
[TABLE]
implying that the matrix could be possibly redefined to
[TABLE]
However, now the matrix is not symmetric. Note that the brackets in the expression (522) are not unique, and one can for example pull the out, so that the expression reads
[TABLE]
Let’s multiply the last term in (529) by , it calculates
[TABLE]
Obviously, it is beneficial to move this term to the vector . Therefore, redefining and , the solution reads
[TABLE]
Now the matrix is again symmetric, and is just a vector. This solution with matrix (533) is slightly nicer than the solution with matrix (523). Alternatively, one can write the solution in its full form
[TABLE]
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 . The solution (532) in the index notation reads
[TABLE]
and direct calculation yields
[TABLE]
where it was useful to define new vector , by adding the quantity
[TABLE]
to the vector . By separating the part, the solution (532) for is therefore re-arranged to the following form
[TABLE]
The solution (540) has a nice advantage that it is decomposed with respect to the magnetic field lines. I.e. projecting the along the field lines by performing cancels the first term, and directly yields vector , that could be defined as vector .
5.1.1 Employing non-dispersive induction equation
Considering the long-wavelength low-frequency limit, the solution (534) for vector and the solution (542) for vector can be further simplified, by using the usual non-dispersive CGL/MHD induction equation, that can be written in the following form (see later in the text)
[TABLE]
By further applying at the induction equation, the last term on the r.h.s. disappears
[TABLE]
and vectors (534), (542) simplify to
[TABLE]
The result (546) can be further rewritten with vorticity , by using identity , finally yielding
[TABLE]
Results (540), (541), (547) are equivalent to eq. 49-51 of Ramos (2005) (see also eq. 11-13 of Macmahon (1965)). Further evaluation of nonlinear FLR corrections to higher precision is out of scope of this simple guide, and the interested reader is referred to Ramos (2005). We will continue in a much simpler linear approximation.
5.2 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 to derive the simplest FLR corrections directly in the linear approximation. Let’s start again, and work with the general equation (492). To derive the simplest finite Larmor radius corrections to the pressure tensor, one has to cancel the contributions on the right hand side of (492) and also neglect the heat flux contributions, which yields
[TABLE]
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 . 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 that is prescribed to be in the z-direction, and therefore . 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 along . However, it is important to note that the procedure of evaluating some terms with , while 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 (492), 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’s evaluate (548) along . On the left hand side we calculate each component in the matrix separately. Since , it is obvious that non-zero contributions are obtained only if , so the index , and . Now, for example the component calculates and the only nonzero contribution is obtained for index , which yields , where in the last step we used that the FLR tensor is symmetric. The entire left hand side and it is straightforward to calculate each component, which yields
[TABLE]
The gyrotropic pressure tensor is of course
[TABLE]
The first term on the r.h.s. of (548) calculates
[TABLE]
and using the matrix notation
[TABLE]
The second term in (564) contributes because of the induction equation . For now, we will not use the induction equation and keep our FLR calculations in the above form. Let’s continue with the second term on the r.h.s of (548). 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 , so that in (564) is replaced by . The remaining terms on the r.h.s. of (548) are straightforward to evaluate since no derivatives of are encountered, and the third term \bm{p}^{\textrm{g}}\nabla\cdot\bm{u}\big{|}_{\bm{\hat{b}}=\bm{\hat{b}}_{0}} is equal to matrix (556) multiplied by . The last two terms of (548) are together evaluated as
[TABLE]
Now we need to subtract the scalar CGL pressure equations obtained previously. Since , the scalar pressure equations read
[TABLE]
By canceling terms coming from these two equations (affecting only the diagonal components), the entire equation (548) rewrites
[TABLE]
The identity evaluated for implies that and the identity further implies that . The components of the pressure tensor therefore read
[TABLE]
Finally, the time derivative of the unit vector
[TABLE]
which evaluated for , yields
[TABLE]
The general induction equation reads
[TABLE]
which evaluates in the approximation as
[TABLE]
These expressions are used in the FLR tensor components and , and the FLR tensor reads
[TABLE]
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 and which yields the leading-order FLR pressure tensor evaluated along the magnetic field in the form
[TABLE]
which is equivalent to equation 5 of Yajima (1966). Perhaps surprisingly, the result is also equivalent to the highly collisional eq. (2.20)-(2.24) of Braginskii (1965), after one prescribes and after terms containing the collisional time are “ignored” ( where is the usual collisional frequency). A proper collisionless limit can not be obtained from the expressions of Braginskii (1965), since the results are proportional to (and also ).
Furthermore, evaluation along 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
[TABLE]
Notably, the tensor is very different than eq. (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. (2009); Hunana et al. (2011); Kobayashi et al. (2017); Perrone et al. (2018), or for the isotropic case see for example Ghosh & Parashar (2015), and references therein.
5.3 Hall-CGL-FLR1 fluid model
The normalized, linearized, Fourier transformed system of equations written in x-z plane reads (dropping tilde everywhere)
[TABLE]
where , and the first-order FLR tensor is
[TABLE]
The dispersion relation for the Hall-CGL-FLR1 fluid model can be written as
[TABLE]
where the polynomial on the right hand side
[TABLE]
The solution for strictly parallel propagation is
[TABLE]
For strictly perpendicular propagation the solution is
[TABLE]
and the result is consistent with eq. 17 of Yajima (1966), after neglecting the electron inertia in that paper (by ). 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 . The result (606) is also equal to eq. 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 (606) actually has an opposite sign in front of the FLR term . For further refinement with FLR2, see eq. (643) and the subsequent discussion.
5.4 CGL-FLR1 fluid model (no Hall term)
If the Hall term is neglected, the general dispersion relation of the CGL-FLR1 fluid model reads
[TABLE]
where the CGL parameters are (600) and the FLR1 polynomial on the right hand side has the form
[TABLE]
where . For strictly parallel propagation, the solution for the whistler and ion-cyclotron waves reads
[TABLE]
For strictly perpendicular propagation the solution is equivalent to (606) since the Hall contributions vanish.
5.5 2nd-order FLR corrections (FLR2)
If a 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 at first order. The pressure tensor equation (492) is simplified to
[TABLE]
For clarity, let’s first consider the case with the heat flux . Compared to the first-order FLR corrections in the previous section, we just have one more matrix present, . As before, the FLR tensor must be symmetric and also satisfy conditions and . Following the previous derivation, the system (580) now includes also the time derivative of the FLR tensor on the right hand side. Let’s write down the system directly in the linear approximation. The system reads
[TABLE]
There are only four FLR components that have to be considered and which can be rewritten separately as
[TABLE]
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 direct numerical simulations perspective it is preferred not to introduce 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 et al. (2016); Ng et al. (2017)). The above equations are therefore typically simplified by expanding the FLR tensor to the first and second-order as and the time derivative of the second-order tensor is neglected, i.e. and the system reads
[TABLE]
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 to be equal to set (588) and move the Hall-term to , or it is possible to keep the Hall term in that is equal to linearized set (586). The first choice, i.e. if is set (588) yields the second-order components
[TABLE]
The second choice, i.e. when is equal to linearized set (586) yields
[TABLE]
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.
5.6 Hall-CGL-FLR2 fluid model
The second order FLR tensor (635) is normalized and written in Fourier space according to
[TABLE]
where for cold (and massless) electrons in the x-z plane
[TABLE]
We purposely wrote the equations with so that the generalization to a more elaborate electric field will be straightforward. Note that the Hall contributions (638) in the FLR tensor completely vanish for isotropic temperatures , regardless of the form of , i.e. regardless if the electrons are cold or not.
To write down the solution for parallel propagation it is useful for brevity to introduce the following quantity
[TABLE]
and the parallel propagating whistler and ion-cyclotron modes factorize as
[TABLE]
and the two solutions are
[TABLE]
with two further solutions corresponding to on the left hand side. It is perhaps useful to rewrite expression (640) to the “ion-cyclotron” resonance form
[TABLE]
that clearly shows that the ion-cyclotron resonances (in the parallel direction) become more and more complicated. For strictly perpendicular propagation , where notably the Hall contributions (638) vanish, the solution of the model is
[TABLE]
It is useful to compare this result with the solution obtained for FLR1 corrections, eq. (606). For direct comparison, expanding the solution (643) for small yields
[TABLE]
Importantly, in comparison to FLR1 solution (606), the FLR2 solution (644) has an opposite sign in front of the FLR term . The equation (644) is equivalent to eq. 32 of Del Sarto et al. (2017), obtained by solving the pressure tensor equation without the expansion to and , and by expanding for small wavenumber afterwards. The difference between the FLR1 solution (606) and the FLR2 solution (644) is perhaps even more clear, when written in physical units for and , and by using the gyroradius , in the following form
[TABLE]
Equations (645), (646), (647) are eq. 34, 32, 56 of Del Sarto et al. (2017). As discussed in that paper, the FLR1 model does not only fail to capture the correction to the Alvén velocity, it also introduces correction to the thermal velocity with incorrect sign. The kinetic result can be obtained for example from eq. 2.10 of Mikhailovskii & Smolyakov (1985), when written for cold electrons , and rewritten to our notation (thermal speed in that paper does not have a factor of 2, so ).
Coming back to the generally oblique propagation, the second approach to write the FLR2 tensor, i.e. when the time derivative of the Hall term is not neglected and kept in the first-order tensor with components
[TABLE]
yields the second-order tensor
[TABLE]
The resulting dispersion relation for strictly parallel propagation is slightly different
[TABLE]
with explicit solutions
[TABLE]
For isotropic temperatures the result is equivalent to (641) since the Hall term contributions for the FLR tensor disappear for all propagation directions. Also, for perpendicular propagation, the dispersion relation is naturally equal to (643).
5.7 FLR corrections with gyrotropic heat flux
The equation (612) contains the divergence of the gyrotropic heat flux, , which in the index notation calculates
[TABLE]
Evaluation along magnetic field yields and the matrix is evaluated as
[TABLE]
where in the last ’zz’ component we also used , but for brevity we left the divergence intact in all expressions. Note that , which is not a problem as will be described below. The pressure tensor equation is here approximated as
[TABLE]
and evaluated with respect to . 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 along read
[TABLE]
The requirement again is that and . Direct calculation yields
[TABLE]
Note that after subtraction of the perpendicular pressure equation, the and components are again anti-symmetric, i.e. , and the system indeed satisfies that . Also, considering the component, after subtraction of the parallel pressure equation all the terms cancel, i.e. and the system indeed satisfies .
Finally, expressing the above system in the linear approximation, with the assumption that mean heat flux values are zero, i.e. , , the only gyrotropic heat flux contributions that remain are and , in equations (664), (665). The FLR tensor is again separated to and similarly to the Hall-term contributions, it is again a matter of choice if the heat flux contributions are pushed to or kept in . We prefer the first choice, i.e. when is equivalent to (588), which yields the second order tensor
[TABLE]
Note that we could have derived the contributions and in a much quicker way, if the expression (652) 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 along is performed. Right now, we can not use these FLR corrections to obtain a dispersion relation, since we have no means to determine a closure for the gyrotropic heat flux . We would have to consider CGL2 fluid model or more complicated Landau fluid models, that contains evolution equations for and . 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.
5.8 FLR corrections with non-gyrotropic heat flux vectors (FLR3)
It is possible to further increase precision of the FLR corrections by considering non-gyrotropic heat flux contributions. Detailed calculations with the non-gyrotropic heat flux vectors and are presented in the Appendix D. At the linear level, these non-gyrotropic heat flux vectors contribute to the linearized pressure equations
[TABLE]
and also to the linearized equations for the FLR tensor
[TABLE]
These equations are equivalent to equations 18-21 of Goswami et al. (2005) (by noting that the vector , so at the linear level ). The non-gyrotropic heat flux vectors are separated to the first and second order
[TABLE]
By using the nonlinear results for the first-order vectors derived in the Appendix D, eq. (D36), (D58), yields that at the linear level
[TABLE]
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 the Appendix D, see equations (D24), (D25) and (D59), (D60). It is important to emphasize that the expressions (674), (675) were derived for a bi-Maxwellian distribution function. As discussed later in the text, this is achieved by prescribing closures for the gyrotropic 4th-order moment in the form , , and . In (674), (675), we additionally neglected the perturbations (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 , and expressions similar to (674), (675) can be derived.
The second-order vectors are also derived in the Appendix D, this time directly in the linear approximation, eq. (D86), (D87) and the vectors read
[TABLE]
As a reminder, . For clarity, we wrote the expressions above explicitly in the index notation. The non-gyrotropic heat flux vectors are perpendicular to , and the non-zero components are for index . The expressions can be written in a more elegant vector form, by defining vector , vector , and matrix
[TABLE]
so that , yielding
[TABLE]
which are equations 53, 54 of Passot et al. (2012).
Normalizing the equations (dropping the tilde), Fourier transforming, and writing them in the x-z plane (with ) yields
[TABLE]
where the terms with Hall electric field components are specified in (638). 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 , . It is useful to briefly explore what fluid models are obtained, if we decide to keep only the first-order vectors , , and ignore the second-order corrections , . By considering only the , contributions, we can either put them into to or . However, it can be shown that both options yield dispersion relation for the perpendicular fast mode (long-wavelength, in physical units)
[TABLE]
Rather amusingly, in comparison to the FLR2 result (646), 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 , . 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’s state the entire model in all of its beauty. The linearized, normalized, Fourier transformed equations written in the x-z plane read
[TABLE]
where each component of the (FLR) pressure tensor , and each component of the (FLR) heat flux vectors , , is separated to
[TABLE]
The components of the tensor are given by
[TABLE]
where the Hall electric field contributions for cold massless electrons read
[TABLE]
and the components of the non-gyrotropic (FLR) heat fluxes are given by
[TABLE]
Importantly, the expressions contain both the 1st and 2nd order contributions from the heat flux vectors , . To explicitly see where the closure for the gyrotropic heat fluxes was performed, we scratched the terms involving these quantities. Additionally, we also scratched contributions from , which are here zero at the linear level.
Let’s check the solution of the Hall-CGL-FLR3 fluid model for the perpendicular propagation (), which in this case reads
[TABLE]
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
[TABLE]
and in physical units
[TABLE]
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 (700), instead of the correct FLR correction to the thermal speed, one obtains . 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 () ion-cyclotron and whistler modes yields the following dispersion relation
[TABLE]
where the quantities and . 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 to use analytic software such as Maple or Mathematica.
Moving the Hall-term to
For completeness, just in case we want to investigate later the influence of the Hall-term, moving it to yields the following dispersion relation for the parallel propagation
[TABLE]
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 in the equations of the Hall-CGL-FLR3 fluid model, the model greatly simplifies. The parallel magnetic field , the ion-acoustic mode decouples, and both the first and second order contributions to , , . The normalized, Fourier transformed equations written in the x-z plane read
[TABLE]
where the components of the non-gyrotropic (FLR) pressure tensor are given by
[TABLE]
and the components of the non-gyrotropic (FLR) heat flux read
[TABLE]
When the entire is neglected, yields the Hall-CGL solution
[TABLE]
Neglecting the contributions yields the Hall-CGL-FLR1 solution
[TABLE]
and neglecting the heat flux yields the Hall-CGL-FLR2 solution
[TABLE]
Finally, the full model yields Hall-CGL-FLR3 solution
[TABLE]
The quantities and used in the solutions are defined as
[TABLE]
Importantly, the solutions are written here for (and not ), 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 section 4.2 where the Hall-CGL model is discussed with final equation (442). The solutions can be written in various forms, and one possibility is to keep in front of the square roots.
The (parallel and oblique) firehose instability was studied in detail by Hunana & Zank (2017), 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 usual fluid models and 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 (left figure, blue dashed lines) and the Hall-CGL-FLR3 fluid model (right figure, black dashed lines).
The plasma , so the long-wavelength firehose threshold, that we call “hard threshold”, is at . The temperature anisotropy is varied from to . For solutions with the Hall-CGL and Hall-CGL-FLR1 models see Figure 3 and 4 of Hunana & Zank (2017). In the WHAMP code the value of cannot be used, and 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 ).
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 analyzed in Figure 8, where instead of the imaginary phase speed, the growth rate is plotted.
The same conclusion is obtained. For the temperature anisotropy , all simpler fluid models (Hall-CGL, FLR1, FLR2) are fully stable for all range of wavenumbers. In contrast, the FLR3 model still develops a strong firehose instability. Therefore, it is the non-gyrotropic heat flux 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 , where the quantities and . At the hard firehose threshold, the expression under the square root of the FLR2 solution (715) can be factorized as
[TABLE]
implying that the FLR2 model is always stable for all values of and . Such a factorization can not be achieved for the FLR3 model, and the solution (LABEL:eq:firehose-FLR3) still becomes unstable at some range of wavenumbers, where
[TABLE]
which is the (strictly parallel) firehose instability criterium of the FLR3 model. We note that the model can become unstable also for the temperature anisotropy . 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 is shown in Figure 9. The left Figure is from Hunana & Zank (2017), 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. The right Figure 9 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’s consider the Hall-CGL model. Assuming firehose unstable regime, the frequency of the Hall-CGL solution (713) can be split into , where the imaginary part . Then by calculating , yields the wavenumber and the maximum growth rate in the following form
[TABLE]
Similarly, the Hall-CGL-FLR1 model yields
[TABLE]
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 plane are shown in Figure 10. The plasma and the temperature anisotropy , so all the fluid models are in the firehose unstable regime. We compare four different fluid models, Hall-CGL (top left), and enhancement with FLR1 (top right), FLR2 (middle left) and FLR3 (middle right).
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 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 bottom left figure is the FLR3 model, and the bottom right figure 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.
6.3 Hellinger’s contours for Hall-CGL-FLR3 model
Here we prescribe a fixed value for the maximum growth rate, , and plot solutions in the plane, which is shown in Figure 11. The left panels show the parallel firehose instability, and the right panels show the oblique firehose instability. The top panels are plotted with the usual logarithmic scales, and the bottom panels with linear scales. Additionally, only solutions with 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 (left panels) and solid red lines (right panels). 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 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 . This does not matter for the parallel propagation, since 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 can not 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 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 (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.
7 Heat flux tensor equation
Multiply the Vlasov equation by and integrate over the velocity space. Naturally, the other possibilities are to multiply by , or , but 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
[TABLE]
i.e. the operator represents all possible cyclic permutations. For the first term we will need the following identities
[TABLE]
and the entire first term of the integrated Vlasov equation calculates
[TABLE]
For the second term we will need identities
[TABLE]
and the second term calculates
[TABLE]
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. . In the expressions above and . The last term can also be rewritten as , where the expression is familiar from MHD. For the third term we will need identities
[TABLE]
and the entire third term calculates
[TABLE]
For the fourth term we will need identities
[TABLE]
and the entire fourth term calculates
[TABLE]
Combining all the results together , the entire heat flux tensor equation obtained by direct integration of the Vlasov equation reads
[TABLE]
Now we will need to use 3 different momentum equations that will cancel various terms. We will also need identity
[TABLE]
We need to multiply the momentum equation for by , equation for by , and equation for by . All the three momentum equations that we want to subtract from the heat flux equation can be written together as
[TABLE]
and because of identity (739), this is equivalent to
[TABLE]
Note that because of the symmetric operator, it does not matter if the tensor is applied on the momentum equation from the left or right, since all the expressions are symmetric in this regard, i.e. for example \big{[}(\nabla\cdot\bm{p})\bm{p}\big{]}^{S}=\big{[}\bm{p}(\nabla\cdot\bm{p})\big{]}^{S}. By subtracting (741) from (738), the final heat flux tensor equation reads
[TABLE]
By defining the cyclotron frequency vector , this equation identifies with equation (A5) in Chust & Belmont (2006). By using scalar cyclotron frequency defined with respect to as , the heat flux tensor equation is written explicitly in the index notation
[TABLE]
which is equivalent to equation (4) in Goswami et al. (2005) with the term
[TABLE]
7.1 Heat flux tensor decomposition
Considering only the gyrotropic part of the heat flux tensor cube, the heat flux is decomposed to the scalar parallel and perpendicular heat flux components according to
[TABLE]
and which in the index notation reads
[TABLE]
This part of the pressure tensor represents only the gyrotropic part, and the full heat flux decomposition can be written as . In the heat flux tensor equation, the term is proportional to the cyclotron frequency and the equation rewrites
[TABLE]
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
[TABLE]
At first look, it is not that obvious that the decomposition (745) satisfy this equation. It is however possible to verify that indeed
[TABLE]
Putting all terms together
[TABLE]
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
[TABLE]
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 . For the clarity of the presented material, here we separate the non-gyrotropic heat flux to a separate term , and write the heat flux tensor equation in the following form
[TABLE]
where
[TABLE]
We will address the non-gyrotropic heat flux contributions in the Appendix D.
In a similar fashion to the pressure decomposition (35), we are going to frequently apply double contractions with and . Applying these operators on the heat flux tensor , must yield quantities that are vectors. It is therefore logical to define parallel and perpendicular heat flux vectors
[TABLE]
The scalar parallel and perpendicular heat flux components (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 , so that
[TABLE]
Briefly considering only the gyrotropic heat flux (746) on the right hand side, it is easy to verify that the parallel decomposition indeed works
[TABLE]
For the perpendicular decomposition, it is useful to specifically calculate components
[TABLE]
so that
[TABLE]
which verifies that the perpendicular decomposition (757) is satisfied for . Now we apply the decomposition (757) directly at the definition of the entire heat flux tensor (6), which yields
[TABLE]
and the perpendicular decomposition
[TABLE]
The decomposition (757) obviously works for the entire heat flux , as well as for the gyrotropic part . Similarly to the pressure decomposition, this further yields required properties that the non-gyrotropic heat flux must satisfy. By using (757), the entire heat flux decomposition reads
[TABLE]
Obviously, it would be useful to introduce a triple-contraction operator and instead of to write something like , but we do not want to introduce new notations. Also, an alternative and perhaps prettier expression is to move the to the left hand side, as , but we will keep the first choice. By applying and at the above equation yields
[TABLE]
implying the first requirement for the non-gyrotropic heat flux
[TABLE]
Similarly, the second requirement is obtained by applying and , yielding
[TABLE]
By using the first requirement, the second requirement simplifies to
[TABLE]
where obviously the trace and operators commute. The two requirements in the index notation read
[TABLE]
Instead of decomposing , an alternative and very useful decomposition of the entire heat flux tensor reads
[TABLE]
with the requirement and (or equivalently ). The heat flux vectors (756) therefore satisfy and . 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 , , it is useful to introduce the following decomposition
[TABLE]
The vectors , are referred to as the non-gyrotropic heat flux vectors, and their algebra is addressed in the Appendix D. Here we only state that the can be decomposed to vectors , and tensor according to
[TABLE]
The entire heat flux tensor is of course non-gyrotropic.
Now we need to verify the heat flux contributions (81), (82) that we used in the pressure equations (83), (84). To calculate the heat flux contributions to the pressure equations, we will need
[TABLE]
and the contributions are
[TABLE]
The left hand side of the above equation can be also written naturally as , 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 .
7.2 Parallel heat flux equation
Now we apply the decomposition (757) on the heat flux tensor equation (754) step by step. We consider only gyrotropic heat flux components. Starting with the equation for the parallel heat flux, the first term calculates
[TABLE]
the second term calculates
[TABLE]
Note that in all expressions here, the operator can be naturally moved to the left as . The third term calculates
[TABLE]
For the final fourth term we will need
[TABLE]
and the fourth term calculates
[TABLE]
Combining all the terms yields equation for the scalar parallel heat flux
[TABLE]
where . The heat flux equations contain the fourth-order moment , which we will consider in the next section and at this stage it is unspecified.
7.3 Perpendicular heat flux equation
We will apply the Trace operator on (754) and also perform step by step. In the index notation, we are basically multiplying the entire equation by . The first term calculates
[TABLE]
The second term calculates
[TABLE]
The components are nothing else but the trace of the moment
[TABLE]
The third term in (754) calculates
[TABLE]
where we have used that and . The fourth term in (754) calculates
[TABLE]
and we have to calculate each term separately. Since and by using already calculated equation (790) for and identities and one obtains
[TABLE]
and the final result of (807) is
[TABLE]
Collecting all the terms, one obtains
[TABLE]
and subtracting the parallel heat flux equation (794) and dividing by two yields the perpendicular heat flux equation
[TABLE]
where . The term containing in the above equation can be rewritten to many possible forms, for example
[TABLE]
The last expression is of course consistent with an alternative way for deriving the perpendicular heat flux equation, where instead of doing Trace of the heat flux tensor equation and subtracting the scalar parallel heat flux equation, one can directly apply operators and .
7.4 Scalar heat flux equations continued
The parallel and perpendicular heat flux equations (794), (812) contain expressions for the 4th-order moment , 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 4th order moment can be decomposed to its gyrotropic and non-gyrotropic part, . The gyrotropic part is decomposed according to (838), and it can be shown (see later in the text), that direct calculation yields
[TABLE]
which yields the parallel heat flux equation
[TABLE]
and the perpendicular heat flux equation
[TABLE]
Exact nonlinear expressions for and that represent contributions from the non-gyrotropic heat flux are calculated in the Appendix D, see eq. (D106), (D116). The scalar heat flux equations (816) and (817) are completely general at this stage, since no distribution function was prescribed yet, and no simplification was introduced. These equations are exact.
Nevertheless, equations (816), (817) are very complicated, and at this stage, it is beneficial to simplify. One possibility, is to cancel all the non-gyrotropic contributions , and , 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 (817) that is proportional to . Also, by following derivations in the Appendix D, it is easy to show that the terms , (that represent the non-gyrotropic heat flux ), do not contribute at the linear level. Therefore, the heat flux equations simplify
[TABLE]
Still, no specific distribution function was assumed. However, to evaluate the non-gyrotropic , and correctly evaluate possible cancellations with terms containing , we need to use decomposition (879). Importantly, the decomposition (879) is valid only for perturbations around a bi-Maxwellian distribution function. By keeping only terms that have non-zero contribution at the linear level, see equations (884), (886) later in the text, the terms are evaluated as
[TABLE]
Importantly, both non-gyrotropic contributions in the parallel heat flux equation (818) completely cancel out! The cancellation demonstrates the importance of keeping the non-gyrotropic , if one wants to keep the non-gyrotropic . Also, there is a partial cancellation in the perpendicular heat flux equation (819). The heat flux equations therefore read
[TABLE]
The last term can be evaluated with respect to the mean magnetic field, and assuming that is in the z-direction, the term is equal to , since . Essentially, the term should be written in a fully linearized form , since other non-gyrotropic nonlinear terms were neglected. It is important to emphasize that it is because of this one term, that the equations (822), (823) are valid only for perturbations around a bi-Maxwellian distribution function. If the term is neglected, i.e. if one neglects from the beginning the non-gyrotropic contributions in (818), (819), the heat flux equations are valid for perturbations around any distribution function. We will consider bi-Kappa distribution function later.
7.5 Heat flux equations with “normal” closure
Detailed calculations with the 4th-order moment will be presented in the next section. Here we want to finish the derivation, and by using the bi-Maxwellian “normal” fluid closure
[TABLE]
terms entering the heat flux equations directly calculate
[TABLE]
The use of these expressions in (822), (823) cancels various terms, and the scalar heat flux equations read
[TABLE]
8 Fourth-order fluid moment
The algebra of the fourth order fluid moment can be quite intimidating at first, since the moment is a 4D 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 n-th order fluid moment , see eq. (1031), and therefore for , the equation reads
[TABLE]
The symmetric operator “” is here defined as
[TABLE]
Here we are not that interested in the evolution equation for , 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 and . For the 4th moment , there are obviously only 3 possibilities: , and , and these gyrotropic components will be called and . We have already seen that the double contractions with and , were very useful operators to extract the gyrotropic components for the lower order moments, and . To obtain any scalar quantity from the 4th moment, we obviously need to apply two double contractions. There are 3 possibilities: we can apply twice, we can apply and , or we can apply twice. Not surprisingly, these double contractions with indeed extract the 3 possible gyrotropic parts, as it is easy to verify
[TABLE]
We are now ready to guess how to write the decomposition of the gyrotropic 4th order moment. Motivated with the previous decompositions, it obviously has to be something in the form of
[TABLE]
where we still did not determine how the symmetric operator acts here. Importantly, the symmetric operator “Sym” is not equivalent to the symmetric operator “” that cycles all the indices around, eq. (830). The fluid hierarchy obviously needs two symmetric operators, one unique “” that is used to derive the evolution equation of a given fluid moment , 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 the . For the , one does two double contractions with . The possible choices are ; and , however, all of these choices are equivalent. To obtain the , we perform double contractions with and , where the last operator contains a function . How many different delta functions we can obtain from 4 indices ? There are different possibilities:
[TABLE]
and all of them are non-equivalent. The symmetric operator acting in the second term therefore has 6 components
[TABLE]
Lastly, to obtain the component, one applies two double contractions with that results in combinations of . 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 . However, 3 possibilities have an equivalent pair, ; and and there are only 3 non-equivalent combinations. The symmetric operator acting on the last term in (834) can be determined to be
[TABLE]
The factor is actually not that obvious, and one needs to verify that the decomposition indeed satisfies (833). The entire 4th-order gyrotropic moment is decomposed in the index notation according to
[TABLE]
It might be tempting to rearrange this decomposition to a more compact form
[TABLE]
nevertheless, in actual calculations we find the form (838) to be more useful. It is important to verity that the decomposition (838) really works. A straightforward calculation yields
[TABLE]
which is consistent with (831). The component calculates
[TABLE]
which is consistent with (832). The component calculates
[TABLE]
which is consistent with (833). The decomposition (838) indeed works.
Now we need to verify expressions (814), (815) that were used in the scalar heat flux equations. The first expression calculates
[TABLE]
and the second expression calculates similarly
[TABLE]
Furthermore, by exploring the equation (829), at frequencies that are much smaller than the gyrofrequency, the gyrotropic part of should satisfy
[TABLE]
in the same way that the gyrotropic parts of and satisfied this requirement. By using the gyrotropic (838), it is easy to calculate for example
[TABLE]
and by adding together all the 4 representations of the “” operator, one can indeed verify that all the terms cancel, yielding (853). By decomposing the entire 4th order moment to its gyrotropic and non-gyrotropic part
[TABLE]
the evolution equation (829) can therefore be rewritten as
[TABLE]
8.1 Non-gyrotropic
Basic properties of the non-gyrotropic tensor can be easily determined with a similar procedure as those we used for the non-gyrotropic pressure and the non-gyrotropic heat flux . The decomposition of the entire 4th-order moment is
[TABLE]
and the meaning of the “Sym” operators were specified by (836), (837). By using definitions (831), (832), (833) for , , in the expression above, the full decomposition reads
[TABLE]
Now, by applying twice on both sides of the equation yields \big{(}\bm{r}:\bm{\hat{b}}\bm{\hat{b}}\big{)}:(\bm{\hat{b}}\bm{\hat{b}})=\big{(}\bm{r}:\bm{\hat{b}}\bm{\hat{b}}\big{)}:(\bm{\hat{b}}\bm{\hat{b}})+\big{(}\bm{r}^{\textrm{ng}}:\bm{\hat{b}}\bm{\hat{b}}\big{)}:(\bm{\hat{b}}\bm{\hat{b}}), implying \big{(}\bm{r}^{\textrm{ng}}:\bm{\hat{b}}\bm{\hat{b}}\big{)}:(\bm{\hat{b}}\bm{\hat{b}})=0, which can be also rewritten as
[TABLE]
Similarly, the other two properties are derived by either applying and , or twice , yielding
[TABLE]
which by using the property (859) are further reduced to
[TABLE]
The last two properties can be also written as , . Finally, writing all 3 properties in the index notation
[TABLE]
One would assume that the non-gyrotropic can be evaluated by rewriting (856) as
[TABLE]
and then expand the r.h.s. similarly to we did for the non-gyrotropic (and also the non-gyrotropic heat flux vectors in the 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.
8.2 Bi-Maxwellian distribution
Here we consider a special case of bi-Maxwellian distribution function
[TABLE]
It is useful to use the fluctuating velocity , where and . For the brevity of calculations we introduce
[TABLE]
where here the thermal speeds are meant to be spatially dependent, i.e. they are written with and not with . The bi-Maxwellian distribution therefore reads
[TABLE]
It is useful to remind us the following one dimensional integrals
[TABLE]
and the general integral with where is an integer reads
[TABLE]
The double factorial , with . The identity (868) is easily obtained by performing differentiation of the first result in (867), and this technique is useful in the calculation of fluid moments.
To get more familiar with the bi-Maxwellian distribution, it is useful to verify if integrals over indeed yield the expected fluid moments. Since the macroscopic velocity is independent of , changing from variable to just yields (similarly to substitution y=x+constant that yields dy=dx). We use notation and . Also, because the integrals are evaluated from to , substitution from to does not change these bounds. Integration in velocity space therefore yields
[TABLE]
For simplicity we integrated in Cartesian coordinates (instead of perhaps 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 it is very easy to verify that
[TABLE]
Continuing with the heat flux components yields
[TABLE]
It is important to emphasize, that the zero heat flux values were obtained by assuming some prescribed equilibrium , here bi-Maxwellian. A crucial principle used in constructing fluid models with higher-order moments is that the specific is assumed only for the last retained moment. Indeed, if one closes the hierarchy by prescribing , the CGL fluid model is obtained. However, to go higher in the fluid hierarchy, evolution equations for and cannot be eliminated, even if the hierarchy will be eventually closed, for example at the 4th-order moment level, by assuming a bi-Maxwellian . Or in another words, the zero heat flux values were obtained for a distribution function that is strictly , however, fluctuations/perturbations around 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 4th-order moment. Only this last moment is calculated by assuming specific distribution function, that is strictly . For a bi-Maxwellian, the gyrotropic 4th-order moments are calculated according to
[TABLE]
Therefore, for bi-Maxwellian, the following closure can be constructed
[TABLE]
and this closure is known as the “normal” closure, a name suggested by Chust & Belmont (2006).
8.3 Bi-Maxwellian non-gyrotropic contributions
Let’s 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 first few terms, one can express the 4th-order moment through 2nd-order (pressure) moments. The procedure is written down in Oraevskii et al. (1968), and the entire 4th-order moment can be decomposed in the following way: see equation (2.31) in Oraevskii et al. (1968) (where a small typo on the l.h.s. is present, where instead of , there should be ), equation (24) in Goswami et al. (2005), equation (9) in Passot & Sulem (2007)
[TABLE]
where the non-gyrotropic contributions read
[TABLE]
and where terms were neglected. By using , it can be indeed shown by direct straightforward calculation that the gyrotropic part
[TABLE]
This agrees with the decomposition (838) valid for a general distribution function, after one specifies that for a bi-Maxwellian distribution , , and (last term in (838) has there ).
Now it is possible to calculate the non-gyrotropic contributions , by using equation (879). It is useful to pre-calculate the trace of (879) that reads
[TABLE]
The two terms that enter the parallel and perpendicular heat flux equations (816), (817) calculate
[TABLE]
At this stage, we need to further simplify, and we keep only those terms that contribute at the linear level, yielding
[TABLE]
so the term entering the equation (817) reads
[TABLE]
Note that (\bm{\hat{b}}\cdot\nabla\bm{\Pi}):\bm{\hat{b}}\bm{\hat{b}}=-\bm{\Pi}:\big{(}\bm{\hat{b}}(\bm{\hat{b}}\cdot\nabla)\bm{\hat{b}}\big{)}^{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 form
[TABLE]
The results show that if the contributions are kept in the heat flux equations, the contributions can not just be straightforwardly neglected, since the and contributions partially cancel out.
9 Bi-Maxwellian fluid model, 2nd-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
[TABLE]
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 contributions. The nonlinear model reads
[TABLE]
To develop a vocabulary, let’s first consider the above fluid model with all the non-gyrotropic fluctuations neglected, i.e. when all the terms with , , are neglected. Additionally, let’s neglect the Hall-term in the induction equation. Such fluid model is non-dispersive, and represents generalization of the non-dispersive CGL model. Even though such fluid model was not explicitly considered by Chew, Goldberger and Low, the pressure equations that include the gyrotropic heat flux contributions , were indeed given by Chew et al. (1956). As discussed previously, the authors just used a different notation with and , but their pressure equations contain the gyrotropic heat flux. The authors did not consider evolution equations for the gyrotropic heat fluxes, eq. (894), (895). 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 to call the non-dispersive fluid model that used the above gyrotropic heat flux equations, as the “the 2nd-order CGL”, abbreviated as “CGL2”. The name CGL2 has a nice advantage, 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, 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 , , yields the Hall-CGL2-FLR3 model. Fluid models with the Hall term neglected, can be easily abbreviated as the “CGL2-FLR1”, “CGL2-FLR2” etc. Even though, if the dispersive effects by the FLR corrections are considered, there is really no reason to neglect the simple Hall-term, and such fluid models are discouraged to use. Obviously, the name “CGL2” is extremely beneficial and 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 , so both the parallel and perpendicular heat flux are normalized the same way (similarly, both and are normalized with respect to ). By dropping the tilde, the scalar pressure equations (892), (893) remain unchanged, and the normalized nonlinear heat flux equations read
[TABLE]
Now we need to linearize the system, where one prescribes and . Once normalized, linearized, transformed to Fourier space and written in the x-z plane, the model reads (dropping tilde everywhere)
[TABLE]
As previously, the .
9.1 CGL2 dispersion relation
Neglecting all the dispersive effects, by prescribing , , , and eliminating the Hall term in the induction equation (i.e. terms proportional to 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 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)
[TABLE]
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
[TABLE]
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 8th-order polynomial in (4th-order polynomial in ). Let’s 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 & MHD descriptions. There is no standardized vocabulary on how the modes should be called, perhaps an expression of “higher-order modes” is the most appropriate. For the highly oblique propagation, considering the mirror instability, sometimes an expression of “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 CGL2 model and kinetic description is, that unless a firehose or mirror instability threshold is reached, the fluid modes are un-damped.
For strictly parallel propagation , in addition to the parallel Alfvén mode (909) that was already separated, the solution of the 8th-order polynomial contains another Alfvén mode , which was of course expected and is equivalent to the CGL model, since the components de-couple for parallel propagation. The remaining solutions are
[TABLE]
Solutions (915) and (916) can be obtained by considering pure 1D geometry (where and disappear since the kinetic velocity disappear) and considering only fluctuations in quantities . The heat flux fluctuations in the CGL2 fluid model are therefore responsible for “splitting” of the CGL ion-acoustic mode to two modes (915) and (916).
Considering solutions for strictly perpendicular propagation , the only non-zero mode is the fast mode
[TABLE]
and the dispersion relation is unchanged from the CGL model (for perpendicular propagation the contributions in the equations for 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, , since the CGL2 fluid model contains the correct mirror instability threshold. The correct mirror threshold can be already seen in the quantity , where in contrast to the usual CGL model, the factor disappeared in the last term .
In the highly oblique limit, the fast mode (918) can be separated from the coupled 8-th degree polynomial (actually 4-th degree polynomial in ) in the following way. Importantly, one can not prescribe a completely perpendicular propagation , and a highly oblique limit , or , has to be considered. The polynomial coefficients can be approximated as
[TABLE]
An alternative approximation is to use in the above expressions. The coefficient does not contain any and is very large compared to that all contain . 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 8th order polynomial in (910) with a neat trick
[TABLE]
The trick can be verified by multiplying both brackets term by term, dividing by , and observing that is negligible. The fast mode is always obtained correctly in this way (when the 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 is separated, and the rest of the dispersion is of 6th order in and contains 3 “slow” modes. By canceling out the that is common in all the remaining coefficients (or , an obvious substitution is offered that reads , and that yields a dispersion equation of 3rd order in in the form
[TABLE]
Now, at exactly the mirror threshold , the last term disappears and we have one solution , and two solutions
[TABLE]
that are both positive. Slightly below or beyond the mirror threshold, we can prescribe , where is a small parameter, meaning that for , the system is slightly below the expected mirror threshold, i.e. stable. For , the system is slightly beyond the expected mirror threshold, i.e. unstable. The polynomial transforms to
[TABLE]
It is possible to show that this polynomial has a discriminant121212Discriminant of a cubic polynomial is calculated according to , implying the existence of 3 distinct real roots. Now, when is small, the two solutions (925) will remain almost the same and change only slightly, importantly, the solutions will remain positive. Solutions of a cubic polynomial form the last coefficient in that polynomial, , here equal to , and here and are positive. Therefore, (slightly below the threshold) implies ; and (slightly beyond the threshold) implies . The negative solution (in ) represents the mirror instability, which finishes the analytic proof that the (highly oblique) mirror threshold in CGL2 model is indeed . 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 (926) for several values of , and it is indeed possible to conclude that for (slightly below the threshold) the signs of solutions are ; and that for (slightly beyond the mirror threshold), the signs of solutions are . 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 (Sulem & Passot, 2012). Considering Landau fluid models, the mirror instability was numerically investigated in detail for example by Passot & Sulem (2007); Passot et al. (2012); Sulem & Passot (2015). As a suggestion for possible future work, it would be beneficial to numerically calculate “Hellinger’s” contours (see subsection 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.
9.3 Hall-CGL2 dispersion relation
Considering the Hall term, the dispersion relation of the Hall-CGL2 fluid model reads
[TABLE]
where the l.h.s. represents the CGL2 dispersion relation (910) and the r.h.s. is the Hall-term contribution. For strictly parallel propagation, the solutions are the usual Hall-CGL whistler and ion-cyclotron waves , , that are accompanied by the CGL2 solutions (915), (917) and (916). For strictly perpendicular propagation, the solution is (918) 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 to work with the full system of linearized equations and solve it numerically or to use analytic software such as Maple.
9.4 “Static” closure - generalized isothermal closure
There exists even simpler fluid closure that recovers the correct mirror threshold, called “static” or “quasi-static” closure (Constantinescu, 2002; Chust & Belmont, 2006; Passot et al., 2006). By using heat flux equations (894), (895) and considering quasi-static regime with , , yields
[TABLE]
where the FLR pressure tensor was neglected in the second equation. It is useful to define which represents gradient along the magnetic field lines, and to use the identity
[TABLE]
further yielding
[TABLE]
Eq. (930) describes evolution of parallel and perpendicular temperature fluctuations along the magnetic field lines, and it can be verified that the solution is
[TABLE]
where . Solution (931) is referred to as the “static” closure and it is for example equivalent to eq. 19, 20 of Passot et al. (2006), and eq. 2 of Constantinescu (2002). For isotropic mean temperatures , eq. (931) yields and both parallel and perpendicular temperatures are isothermal. A similar result is obtained for general with , which yields . The closure (931) can be therefore viewed as a generalization of isothermal closure in the presence of temperature anisotropy and variations of magnetic field strength. The closure is used directly in the momentum equation (890), and the time-dependent pressure equations (892), (893) are disregarded. Therefore, even though the closure was derived by considering 4th-order moments, the resulting fluid model is actually simpler than CGL2 and CGL models. As discussed for example by Passot et al. (2006), the closure prescribes limitation for the minimal value of that is required to prevent temperature singularity. Or in other words, by separating the requirement reads , meaning that the fluctuating can not 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 can not be too deep. To easily analyze dispersion relations, it is useful to write the closure at the linear level, in the following form
[TABLE]
The result can be obtained by linearizing (931) or (930), and at the linear level . The normalized closure reads (tilde are dropped) and . Considering the simplest non-dispersive model (by neglecting the Hall term and ), the dispersion relation of this fluid model reads
[TABLE]
In comparison with the dispersion relation (388) of the adiabatic CGL model, the “static” closure eliminates the erroneous factor in the coefficient and yields the correct mirror threshold. The erroneous 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 () are and . For perpendicular propagation () the solution is . 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” closure (932) are indeed crucial in recovering the mirror instability, let’s neglect them and quickly consider strictly isothermal closure ; , which yields the following dispersion relation
[TABLE]
As can be seen from the 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 . As a double check, after prescribing , the dispersion relation (934) is equal to (933), since the contributions in the “static” closure disappear. Also, we considered models with polytropic indices in Section 3.11, and the dispersion relation (934) is consistent with (368), after prescribing , .
10 Bi-kappa fluid model (BiKappa)
10.1 Bi-kappa distribution function
Often in the solar wind and in the space plasma physics, distribution functions with elongated high-energy “suprathermal” tails are observed, that cannot be modeled with the bi-Maxwellian distribution. A slightly more general distribution function is the bi-kappa distribution, that contains one free parameter , see for example Summers & Thorne (1992); Pierrard & Lazar (2010); Livadiotis & McComas (2013) and references therein. The parameter can be different in the parallel and transverse directions (see e.g. Basu (2009); dos Santos et al. (2015)). Here we consider only one free parameter. The distribution has the following form
[TABLE]
where the generalized thermal speeds
[TABLE]
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
[TABLE]
Why the generalized thermal speeds are defined this way will become clear from later calculations of moments, and the choice of the power law index , instead of perhaps more logical , is attributed to a purely historical reasons. Calculating 2nd-order (pressure) integrals with this distribution requires for convergence , and this value is also the minimum value for the generalized thermal speeds (936) to have real values. Sometimes, in the definition of the isotropic kappa distribution, a value of is used and . Again, for the brevity of calculations, it is useful to introduce , , and rewrite the bi-kappa distribution (935) to the following form
[TABLE]
The important integral is
[TABLE]
For the critical value of , the integral , and since , the integral diverges. For an isotropic kappa distribution , the 3-dimensional integrals can be conveniently evaluated by transformation to spherical co-ordinates
[TABLE]
In this case, the important integral is
[TABLE]
and by a direct substitution
[TABLE]
For an isotropic kappa distribution therefore
[TABLE]
which explains the normalization factors. For an anisotropic bi-kappa distributions , the integration is slightly more complicated because in contrast to a bi-Maxwellian distribution, the integration cannot be fully separated to two independent integrations over parallel and perpendicular velocity components, and have to be performed successively. It is possible to use a cylindrical coordinate system and integrate
[TABLE]
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 . A small added benefit is that we do not have to follow if the integrals are or , and we can drop writing the boundaries.
The most important integral is a slight generalization of (939) in the form
[TABLE]
where for simplicity are assumed to be positive constants. So one can calculate
[TABLE]
where importantly, the “power law” changed from to . Performing a successive 3-dimensional integral therefore yields
[TABLE]
and integrating the bi-kappa distribution over the velocity space
[TABLE]
which verifies that the normalization constants in the bi-kappa definition (938) are indeed correct, and that . For later calculations, it is useful to write down two partial integrals when the integration is done over and that read
[TABLE]
and
[TABLE]
Continuing with the calculation of velocity moments, it is easy to show that , which verifies . 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 to a bi-Maxwellian distribution use a trick with the differentiation with respect to as
[TABLE]
where we have also used a property of Gamma function , and the last integral requires . A slightly more general integral is
[TABLE]
which can be used to evaluate the following integrals
[TABLE]
Since we like to double check everything, the above integral can be also calculated as
[TABLE]
For the calculation of the parallel pressure, we first integrate over so we can use result (954) and then over
[TABLE]
where we have used . Calculating the perpendicular pressure, we here first integrate over so we can use result (951) and then over
[TABLE]
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 by calculating the parallel and perpendicular pressure integrals, yields the required forms ; . Then the split to ; 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 bit obscure power law instead of more logical , which is attributed to purely historical reasons. If one wants to define an anisotropic power-law distribution, a first guess would be
[TABLE]
which then yields normalization constants
[TABLE]
Then, by calculating the parallel and perpendicular pressures yields the requirement . By further writing yields a distribution function
[TABLE]
with generalized thermal speeds
[TABLE]
The convergence of the pressure integrals requires , 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 (965) is really correct, one can of course substitute , which fully recovers the original bi-kappa distribution (935), since . From now on, lets continue calculations with the original bi-kappa distribution (935), that uses the 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 are anti-symmetric in , and by a direct calculation
[TABLE]
implying
[TABLE]
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 4-th order moment level. To calculate the 4-th order fluid moments, we will need the following integrals
[TABLE]
The moment then calculates
[TABLE]
where in the last step we used that with . The moment calculates
[TABLE]
and finally the moment calculates
[TABLE]
To summarize, the following closure can be constructed
[TABLE]
where for the bi-kappa distribution (935), the coefficient reads
[TABLE]
It is noteworthy that all three 4-th order moments are just multiplied by the same constant . Obviously, in the limit , the , 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
[TABLE]
where the 3 proportionality constants have to be determined, once a specific distribution function is prescribed. Such unspecified closure, can account for a very wide class of distribution functions. Also, a much more complicated distribution functions can be modeled by considering multi-fluid models. For example, one could consider 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. (2005, 2008); Pierrard et al. (2011) and references therein.
Here we consider closure (976), (977). 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’s completely neglect the non-gyrotropic and in the heat flux equations (818) and (819). We need to calculate
[TABLE]
and the nonlinear heat flux equations for the bi-kappa fluid model read
[TABLE]
Obviously for (i.e. for ), the bi-Maxwellian heat flux equations are obtained. The model is now closed, and it is accompanied by equations (889)-(893). 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 a bi-Maxwellian nonlinear simulations, the temperature anisotropy is free to develop, and complicated plasma heating process can be studied (we recommend to add the FLR pressure stress forces to the scalar pressure equations). We call the non-dispersive fluid model that uses the above heat flux equations with abbreviation “BiKappa”. If the Hall term is considered in the induction equation, “Hall-BiKappa” fluid model is obtained. Similarly, by considering FLR1 and FLR2 corrections to the pressure tensor, yields “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
[TABLE]
where we have used the same plasma beta definition as always . The heat flux equations are accompanied by the linearized system (898)-(906).
To consider FLR3 corrections, one needs to derive non-gyrotropic heat flux vectors , for the bi-kappa distribution. For the first-order vectors the derivation is straightforward, and by following calculations in the Appendix D, it is easy to show that the nonlinear
[TABLE]
Here we neglected perturbations , since we will not consider Landau fluid models for bi-kappa distribution. For , the results are of course equivalent to bi-Maxwellian expressions (D36), (D58). The second-order heat flux vectors are not addressed here, since we would have to derive for bi-kappa distribution.
10.3 BiKappa dispersion relation
The BiKappa dispersion relation can be written in the following form
[TABLE]
and for it is of course equivalent to the CGL2 dispersion relation (910). The oblique Alfvén mode is not effected, and is identical to (909), implying that the long-wavelength threshold of the oblique firehose instability is identical as well.
For strictly parallel propagation , the solutions of the remaining 8th-order polynomial are
[TABLE]
Note that in the last solution (995), the expression under the square root remains always positive, regardless how large is (i.e. when the -index approaches 5/2). For strictly perpendicular propagation , the fast mode is not effected, and is equivalent to (918).
10.4 Mirror instability
By considering the highly oblique limit , and by assuming that is not enormously large (i.e. that is not very close to ) so that can be neglected with respect to , the polynomial coefficients can be approximated as
[TABLE]
The fast mode can be quickly separated with the trick (923), and by using substitution , the resulting polynomial reads
[TABLE]
At exactly the (assumed) mirror threshold , the coefficient disappears, so we have one solution ; from the remaining polynomial the can be factored-out and the solutions are
[TABLE]
both always positive. Slightly below or beyond the mirror threshold, one can prescribe , where is a small parameter and the polynomial transforms to
[TABLE]
and which for is equivalent to the bi-Maxwellian polynomial (926). As previously, for small the solutions (1001) remain almost the same, and positive. Furthermore, since multiplying all 3 solutions together must yield , it implies that below the threshold the 3rd solution is positive, and beyond the threshold the 3rd solution is negative, which analytically proves that the (highly oblique) mirror instability threshold is , regardless of . The long-wavelength limit or “hard” mirror threshold for the bi-kappa distribution is therefore the same as for the bi-Maxwellian distribution.
Conclusion
We investigated solutions of the (non-dispersive) BiKappa fluid model. The calculations were done with closure (976) for a general coefficient, i.e. without using the specific value (977) for the bi-kappa distribution. The calculations are therefore more general and valid for any distribution function that can be closed as (976). We conclude that all 3 long-wavelength “hard” thresholds
- for the parallel and oblique firehose instability and the highly oblique mirror instability - are not effected by the value of .
10.5 Hall-BiKappa dispersion relation
The Hall-BiKappa dispersion relation can be written in the following form
[TABLE]
where we kept the BiKappa dispersion relation (987)-(991) on the l.h.s., and the Hall term contributions on the r.h.s.. For the bi-Maxwellian value , the result is equivalent to the Hall-CGL2 dispersion relation (927). For strictly parallel propagation, the solutions are the usual dispersive ion-cyclotron and whistler waves (that are not affected by the coefficient) and the non-dispersive solutions (993)-(995) of the BiKappa model. For strictly perpendicular propagation, the Hall term naturally disappears and the fast mode solution is (918).
11 Collisionless damping in fluid models - Landau fluid models
In the previous sections, we have calculated the fluid hierarchy up to the 4th-order fluid moment . We have seen that unless we crossed some instability threshold, the dispersion relations yielded frequencies that were purely real, and there was no collisionless damping present. Remarkably, collisionless damping, i.e. Landau damping, will be absent even if we continue to develop the fluid hierarchy to higher-order moments. In fact, we will see in the next section by considering parallel 1D propagation (where all the fluctuations are along ), that fluid models closed by a bi-Maxwellian fluid closure at higher-order moments than , are always unstable. Such fluid models are physically ill-posed. To go higher in the fluid hierarchy, one needs to incorporate Landau damping phenomenon into the fluid framework, and consider Landau fluid closures.
Nevertheless, to capture the Landau damping, it is not necessary to go beyond the 4th-order moment and closures performed at the 4th moment level describe the Landau damping with very good accuracy. Even though the entire Part 2 of this text is dedicated to Landau fluid closures, it is beneficial to briefly show, why the Landau damping was omitted in the previous calculations. Let’s consider the simplest example of 1D electrostatic propagation. The collisionless Vlasov equation reads
[TABLE]
Direct integration of the Vlasov equation (over ) for a general distribution function , of course yields the usual density equation, . The same result is obtained if one considers specific , as for example 1D Maxwellian , 1D kappa distribution, etc.. The last term of the Vlasov equation vanishes by the integration, since . Now let’s do the calculation differently, by expanding the distribution function , where does not have any time or spatial dependence, so . We can consider Maxwellian where . The electric field is also expanded , and one assumes that there is no large-scale electric field, so . By neglecting term with two small quantities , the Vlasov equation reads
[TABLE]
and transformation to Fourier space yields
[TABLE]
The equation (1006) is actually not that much different from (1004), the Vlasov equation is just fully linearized and written in Fourier space. However, in kinetic calculations, the (1006) 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 (1006) is first divided by , and an expression for is obtained. Only then the equation is integrated over the velocity space, yielding
[TABLE]
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 , where is the cyclotron frequency of the considered species, and is an integer. The resonance for is called the Landau resonance, and the other resonances for are called cyclotron resonances. It is the presence of these resonances in integrals such as (1007), 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 is simple Maxwellian. The integral is just expressed through the famous plasma dispersion function , where , since the plasma dispersion function was specifically developed to describe this integral. Skipping many technical (but very important) details such as the definition of , for Maxwellian the density integral is simply expressed as
[TABLE]
One can define function , that is called the plasma response function. Now for example considering proton-electron plasma with at wavelengths that are much longer than the Debye length, yields the following dispersion relation and numerical solution
[TABLE]
The negative imaginary part represents the Landau damping, and the effect comes from the Landau resonance in the integral (1007). Obviously, integrating the Vlasov equation directly in the form (1004) or (1) neglects the wave-particle resonances in velocity space, which clarifies why traditional fluid models do not contain collisionless damping mechanisms. It is the - the fluctuations/perturbations around - 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” . We find it useful to introduce the following vocabulary:
- •
fluid moments = integrals over a general (unspecified) distribution function , or over a specific . 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 .
- •
kinetic moments = integrals over (perturbations , where a specific equilibrium is assumed), where the wave-particle resonances in velocity space are written down explicitly. The resulting expressions are linear and contain , 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 , instead of the “normal” closure, the 4th-order moment is decomposed in the following form
[TABLE]
The first terms represent the “normal” closure, i.e. when an exact bi-Maxwellian is prescribed. The tilde components represent fluctuations/perturbations around , and a good notation would also be . The decomposition (1010) yields the following form of the nonlinear heat flux equations
[TABLE]
Similarly to decomposition (1010), “kinetic” corrections with wave-particle resonances can also be considered for the non-gyrotropic , which is not addressed in this guide. Neglecting those contributions, the heat flux equations are equivalent for example to eq. 14, 15 of Passot & Sulem (2007), and eq. 9, 10 of Passot et al. (2012).
In Part 2 of this guide, we will calculate the “kinetic moments” for , , , by integrating over . In the 3D electromagnetic geometry, we will use more sophisticated in the gyrotropic limit, which describes both Landau damping and its magnetic analogue, the transit-time damping. The expressions of the entire “kinetic” hierarchy will contain the , or actually the , which is not a fluid variable. Nevertheless, we will see that by introducing the concept of Padé approximation to , that sometimes a rare possibility arises, when the “kinetic” moment can be expressed through lower-order “kinetic” moments. Such a closure, where the last retained kinetic moment (in this case ) is expressed through lower-order moments, in a way that the kinetic 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 are going to be (written in Fourier space)
[TABLE]
The first closure was obtained in the slab geometry by Hammett & Perkins (1990) and the remaining closures were obtained in the 3D geometry by Snyder et al. (1997) (note that thermal speeds in those papers are defined without the factor of 2). In Part 2 of our guide, we are going to map all the Landau fluid closures that can be constructed at the 4th-moment level, and we will identify the most precise closures. A summary of these findings, can be found in Hunana et al. (2018).
12 Evolution equation for n-th order fluid moment
Now we are proficient with calculating the fluid moments for general , and before we start with kinetic calculations, it is useful to derive the time-dependent evolution equation for a general n-th order fluid moment. Before going through this part, we highly recommend to go through the parts of the paper where the pressure tensor and heat flux tensor equations are obtained. Then the algebra and the notation used here will look hopefully logical and natural. We define the n-th order fluid moment according to
[TABLE]
The symmetric operator is defined that it cycles around through all possibilities according to
[TABLE]
and it is compatible with the previously used symmetric operators for a 2nd-order tensor and for a 3rd-order tensor . We multiply the Vlasov equation by the mass m and by the velocities and integrate over . To integrate the first term of the Vlasov equation we will need
[TABLE]
and the first term of the Vlasov equation calculates
[TABLE]
For the second term we will need
[TABLE]
[TABLE]
The second term calculates
[TABLE]
The third term calculates
[TABLE]
For the fourth term we will need
[TABLE]
[TABLE]
The fourth term calculates
[TABLE]
Combining all the results together , the evolution equation for n-th order fluid moment reads (written in the vector notation without indices)
[TABLE]
The final result should have only time derivative of and similarly to the heat flux tensor, we need to use the momentum equation to eliminate . We subtract a total of n momentum equations combined with in the form
[TABLE]
Again, because of the symmetric operator, it does not matter if the is applied to the momentum equation from the left or right. The subtraction also conveniently eliminates the electric field, however it introduces . We also need an identity
[TABLE]
where in the last step we emphasized that it does not matter if the applies only to or to the entire tensor , since the vector product operates only through the first component of the tensor anyway. The final evolution equation reads
[TABLE]
and it is valid for . Evaluation at recovers the heat flux tensor equation (742) and evaluation at recovers the pressure tensor equation (28). Note that , , , however, and very importantly, according to our definition (1016), and it is not equal to the velocity , which is defined according to (4). Additionally, . The momentum equation can still be recovered by evaluating (1028) at , and by dividing by (where naturally the symmetric operator does not have any influence on a vector, and ). Similar equation to our (1028) was also obtained for example by Waelbroeck (2010), eq. 22, even though we did not verify if it is consistent with ours since in that work is defined with instead of our .
We want to always decompose the tensor to its gyrotropic and nongyrotropic part, and to have a space to write the “g” and “ng”, lets move the index down, so that
[TABLE]
By generalizing the results we have seen for the gyrotropic , the same result is obtained for the n-th moment, and at low frequencies the gyrotropic part must satisfy
[TABLE]
which can be indeed viewed as a definition of gyrotropy. The general (1031) therefore rewrites
[TABLE]
yielding an implicit equation for the non-gyrotropic part in the form
[TABLE]
Ideally, an “inversion procedure” should be found for the l.h.s., so that an equation for can be found, and then by expanding the r.h.s. (for example at the leading order by considering only gyrotropic ), a fully nonlinear expressions for can be obtained. Nevertheless, we have seen that the inversion procedure was already quite complicated for the matrix , see eq. (500). For a tensor of 3rd-rank, such as the , the full inversion procedure is actually not addressed in this guide and an advanced reader is referred to eq. 43 of Ramos (2005). The inversion procedure was obtained only for a part of the , that can be decomposed to the non-gyrotropic heat flux vectors , , which is done in detail in the Appendix D. For the rest of the that consists of heat flux tensor , the inversion procedure cannot be avoided. Additionally, we have seen that for 4th-order moment , the eq. (1035) does not seem to yield anything useful even at the linear level, and a different procedure was used, where the specific bi-Maxwellian distribution was expanded in Hermite polynomials. Therefore, for higher-order moments, the eq. (1035) is immensely complicated.
It is useful to briefly consider, how many scalar gyrotropic moments has an n-th order tensor . The scalar gyrotropic moments can be constructed only from combinations of and . It is beneficial to write down the following self-explanatory table
[TABLE]
Obviously, the number of gyrotropic moments for tensor is equal to . The function “int” means integer part. In mathematics, the integer function is often abbreviated with square brackets, as . However, this notation might be possibly confusing (we used square brackets in many other places), and we prefer to write explicitly “int”. Another possibility would be to use the “floor” function.
12.1 Hierarchy of moments in 1D geometry (electrostatic)
Let’s see if we can use the n-th order moment equation (1031) for something useful, since the equation is immensely complicated. Let’s consider the special case of strictly parallel propagation along (in the z-direction), where all the moments (including the velocity) are along , i.e. let’s consider the electrostatic case. Let’s change the index of a moment from to , so that the index is not confused with the number density. Let’s stop writing the subscript on all the variables starting from the pressure, so , , are simply and the l-th order parallel moment is just . We keep the -subscript for to remind us we are writing 1D equations. In 1D geometry, the nonlinear evolution equation (1031) simplifies to
[TABLE]
If one wants to be more general, instead of propagation along , it is possible to consider propagation along magnetic field line , and write the same equation with and . The equation is accompanied by the density and momentum equations
[TABLE]
where for the charge it is useful to keep the species index , so that it is not confused with the heat flux . The last term disappears anyway since in the case considered here . Evaluating equation (1036) for yields
[TABLE]
and so on. The nonlinear equation (1036) can be linearized, in the first step by assuming zero mean value, yielding
[TABLE]
and calculating it for and by prescribing , and in general for odd , yields
[TABLE]
Alternatively, we can split the linearized equations for even and odd according to
[TABLE]
When closing the system, one has to be careful how to deal with =even quantities, since we want to separate the “deviations”, similarly as we did for . For a Maxwellian and and even, we therefore write
[TABLE]
or equivalently
[TABLE]
so for example , , , and linearization of the l-th order moment yields
[TABLE]
Using (1056) in the first term of =even equation (1052), and using the pressure and density equations, yields linear evolution equation for deviations of even moments
[TABLE]
and using (1056) in the second term of =odd equation (1053), and yields linear evolution equation for deviations of odd moments
[TABLE]
Finishing the calculation with evaluations of values, the evolution equations of deviations are
[TABLE]
or by using , equivalently
[TABLE]
Continuing to higher orders, the hierarchy of deviations of Maxwellian moments therefore evaluates
[TABLE]
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 to go beyond CGL2 without Landau fluid closures
One can analyze the dispersion relations easily in physical units, but since in this Part 1 we used normalized units with , and in almost all dispersion relations, let’s rewrite the fluid hierarchy to normalized units, so that everything feels more familiar. The is normalized with , the with , and with . By dropping the normalization tilde as usual (but obviously not for even moments such as , that should be perhaps called ), the equations in Fourier space read
[TABLE]
Now, by using just first 3 equations (1070)-(1072), and closing the system with a closure , yields the CGL model, with solution
[TABLE]
the familiar CGL ion-acoustic mode. By going higher in the fluid hierarchy and using the first 4 equations (1070)-(1073), and closing the system with (which is equivalent to prescribing ), yields the CGL2 model. The dispersion relation is , with solutions
[TABLE]
the already obtained (915), (916). 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 , and use the first 5 equations (1070)-(1074) with the closure . We might call this model CGL3. The dispersion relation reads , and the solutions are
[TABLE]
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 fluctuations on the r.h.s yields
[TABLE]
Therefore, if the , so when (1081) is used, the CGL2 model is obtained. However, when fluctuations in are considered, so when (1082) is used, the term cancels on both sides, with the resulting dispersion relation
[TABLE]
and by prescribing closure , 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 (1070)-(1075), the dispersion relation reads
[TABLE]
and by prescribing closure (which is equivalent to prescribing ), yields a model that we can call CGL4, with numerical solutions
[TABLE]
again, two modes being imaginary, and one unstable.
The variable , written here in normalized units, is easily rewritten to physical units (reintroducing tilde for clarity)
[TABLE]
which is directly related to the “kinetic” variable , that we will use in Part 2 of the text, defined as
[TABLE]
where the absolute value does not make any difference now (since all fluid models here contain ), and will become important only when Landau fluid closures are considered. The hierarchy of CGL dispersion relations then can be written as
[TABLE]
All the models beyond CGL2 contain modes that are unstable.
Double-check
Let’s 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.131313An alternative approach that would describe the Landau damping very precisely, is to keep the fully nonlinear CGL2 equations, and supplement them with a hierarchy of time-dependent linear equations for higher-order moments, closed with an appropriate Landau fluid closure. The CGL and CGL2 systems were already discussed at great length. Let’s check the CGL3 dispersion relation, and let’s work directly in physical units. The system of equations that we have in mind reads
[TABLE]
and the closure is performed by setting . Importantly, the 4th-order moment is kept undetermined, and it is not separated to its “core” and , nor any specific form of a distribution function is prescribed yet (other than the mean values of odd moments ). Calculating dispersion relation of this system (a matrix multiplied by vector ) yields result , which for recovers the CGL3 dispersion relation. Our calculations were therefore done right. Similarly, double checking the CGL4 dispersion, the system (1095) is supplemented with equation (1048), and a closure is performed with . Again, the is left untouched, and dispersion relation of this system recovers the CGL4 result. Checking the CGL5 dispersion relation, the system is supplemented with (1049), and a closure . Again, the is kept untouched and generally undetermined. Calculating the dispersion relation yields , and prescribing Maxwellian recovers the CGL5 result. Going higher in the hierarchy, and using closure recovers the CGL6 results. And finally, keeping all the variables up to undetermined and using closure yields dispersion relation , which for Maxwellian 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
[TABLE]
and for higher orders therefore
[TABLE]
where the validity can be proven by induction and using (1052), (1053). CGLl model (or CGL of l-th order) is therefore defined as performing a closure on moment. In the equations (1096)-(1104) we used the word CGL, even though no specific distribution function (nor any closure) was prescribed yet.
The closure in each model is performed by specifying the r.h.s. of (1103), (1104), which yields dispersion relations
[TABLE]
The results are written in physical units. Alternatively, by using the variable, the results are
[TABLE]
which for specific values recovers results (1088)-(1094).
The results (1105)-(1106) are written for a closure performed on moment , which is dictated by the CGLl vocabulary, and which might be perhaps a bit confusing. Alternatively, one can reformulate the results for closures performed on moment (corresponding name for a Maxwellian is ) and that read
[TABLE]
or rewritten with the variable
[TABLE]
a result reported in Hunana et al. (2018), eq. (10). Importantly, solutions of the above dispersion relations for , will always yield results with some complex numbers, and some solutions will be unstable.
Therefore, the last physically meaningful fluid model (without considering Landau fluid closures) is the CGL2 model, and the last closure is the “normal” closure, . To go higher in the fluid hierarchy, we necessarily need to make connection to kinetic theory, and consider Landau fluid closures. These closures are described in Part 2 of our guide, and for direct comparison with the simplest Landau fluid dispersion relations, see Part 2, Section 3.12 “Parallel ion-acoustic (sound) mode, cold electrons”.
13 Conclusions
We offer a brief summary of the major results discussed throughout the text.
- •
Collisionless plasma is completely described by the kinetic Vlasov equation. By directly integrating the Vlasov equation over velocity space, yields evolution equations for an infinite hierarchy of fluid moments. In addition to the usual density and momentum equations, it is possible to derive an evolution equation for a general n-th order fluid moment (), see eq. (1031). Then, evolution equations for the pressure tensor and heat flux tensor can be obtained by simply evaluating (1031) at and .
- •
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: 1) “Classical” (non-Landau fluid) closures, that neglect wave-particle resonances, and 2) Landau fluid closures, that account for Landau wave-particle resonances and associated Landau damping & 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 4th-order moment with classical bi-Maxwellian fluid closures , where represents bi-Maxwellian value of and its perturbation , since all the fluid models contain higher-order modes that are unstable. This is perhaps the most surprising result discussed in Part 1. For a detailed proof, see section 12.2, where 1D geometry is considered and a bi-Maxwellian closure is used at the n-th order fluid moment, which yields dispersion relation (1111). The same dispersion relation will be valid in 3D geometry for the propagation along the magnetic field, where additional de-coupled modes will be present as well. For , all the fluid models contain unphysical instabilities. The same result is expected for other distribution functions.
- •
The classical bi-Maxwellian closure at the 4th-order moment reads (or equivalently ), and is called the “normal” closure. Therefore, the “normal” closure is the last classical fluid closure ) and beyond the 4th-order moment, Landau fluid closures ) are required.
- •
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 the pressure tensor , the heat flux tensor , the 4th-order moment . The non-gyrotropic parts are also referred to as the Finite Larmor Radius (FLR) corrections, since they represent deviations from gyrotropy at small spatial scales, when the Larmor radius (gyroradius) is not infinitely small.
- •
A general n-th order fluid moment contains scalar gyrotropic moments, where the function “int” means integer part.
- •
The pressure tensor is 2nd-order fluid moment and contains two gyrotropic moments, and . Therefore, any fluid description of collisionless (or weakly collisional) plasmas, requires two separate evolution equations for and . 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 (), since pressure fluctuations in the directions parallel and perpendicular to the local magnetic field remain anisotropic. We often use proton pressure (temperature) anisotropy coefficient .
- •
The simplest fluid model that describes 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 . 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 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 Section 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 and FLR pressure tensor unspecified, it is possible to derive rigorously exact evolution equations for and , see eq. (100)-(101).
- •
Rewriting the pressure equations to an alternative form, eq. (108)-(109), 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 eq. (121)-(124). Only when the anisotropic heating is neglected, i.e. when only total plasma heating studied, the conservation of energy significantly simplify, see eq. (117)-(119).
- •
By using the concept of polytropic indices, the pressure equations in the ideal CGL model can be interpreted as having and . Consequently, the CGL dispersion relation is in general not equal to the MHD dispersion relation, even when . The exception is the Alfvén mode, which for propagates with the same phase speed in both CGL and MHD models, for all the propagation directions .
- •
Polytropic indices can be further related to the number of degrees of freedom through , which for the CGL model yields and . The CGL pressure equations can be therefore viewed as being composed of strongly coupled 1D and 2D dynamics, whereas the dynamics in the MHD model with and can be viewed as isotropicaly 3-dimensional.
- •
In the CGL model, the Alfvén mode propagates with the phase speed , which matches the (collisionless) kinetic theory in the long-wavelength limit. Consequently, for , 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 (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 is present in the CGL model and exists even for , see Section 3.6. The effect is very important for the parallel firehose instability.
- •
Considering strictly perpendicular propagation ( or ), the fast mode in the CGL model propagates with the phase speed , which matches the kinetic theory in the long-wavelength limit. The expression can be rewritten as , and can be contrasted with the MHD result . The factor of 2 in the CGL result comes from , 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 is appropriate, since for example the parallel propagating ion-acoustic (sound) mode is strongly Landau damped in kinetic theory (unless 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.
- •
Nevertheless, when electrons are hot () and the ion-acoustic mode is Landau damped only weakly, the real phase speed obtained from kinetic theory in the long-wavelength limit reads . Thus, in this example the kinetic theory is matched by a fluid model where protons are described adiabatically (with the CGL value ), while electrons are isothermal (with ). Also, a good example for is the real frequency of the Langmuir mode, , which is addressed in Part 2.
- •
For anisotropic temperatures (), some modes can become unstable. A threshold of an instability that is obtained in the long-wavelength (non-dispersive) limit is called “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 matches the kinetic theory. The instability requires and (in the long-wavelength limit) . Since the CGL model is non-dispersive, here the instability growth rate only has a simple 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 (438)-(441), one arrives at a contradiction that for the whistler mode is unstable, and that for 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 looses sense. The problem is resolved by introducing a small causal dissipation into the momentum equations, which is later removed by the limit . The procedure yields causal dispersion relations (449)-(452), where the whistler mode is firehose unstable for both positive and negative , and the ion-cyclotron mode is stable.
- •
The parallel firehose instability occurs for quasi-parallel (small ) propagation directions, with the maximum growth rate at , see Figure 10. In the non-dispersive CGL model it corresponds to the instability of the slow mode. Its threshold obtained for is equivalent to the oblique firehose instability (since ), 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 , , yields that the fast mode is the ion-acoustic (sound) mode and , see Figure 1. After dispersive effects are introduced, which split the Alfvén and slow mode to the ion-cyclotron and whistler mode, it is clarified that the instability is associated with the whistler mode. For quasi-parallel propagation directions the situation is more obvious, since (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 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 , and with respect to kinetic theory the factor is erroneous. The erroneous 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 Section 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 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 uphysical 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 and modifying the parallel polytropic index to . 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 is even below the isothermal value , 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 , 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 , that has two gyrotropic components, and . The system is closed at the 4th-order moment , that have 3 gyrotropic components, . 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 is known as the “normal” closure (see above). We call this non-dispersive bi-Maxwellian fluid model the “2nd-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 of the CGL model (388) and the CGL2 model (910). The Hall term introduces simplest dispersive effects, see dispersion relations of the Hall-CGL model (490), and the Hall-CGL2 model (927).
- •
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 ). 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.
- •
One can construct closures for other distribution functions than bi-Maxwellians. For example, one can consider the bi-Kappa distribution function (935) with a free parameter . The closure then reads , where the coefficient . The closure is restricted to , which can be seen from the coefficient, and the same requirement is necessary for the convergence of the velocity integrals. We call this closure and resulting fluid model simply as “biKappa”.
- •
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 all models contain unphysical instabilities), the closure will be restricted to , which increases with . 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 (987). 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 (1003).
- •
The non-gyrotropic (FLR) corrections to the pressure tensor were studied in Section 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 (532)-(534), or equivalently (540)-(542). 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 as “off-diagonal”, since the diagonal components are non-zero. In the collisionless case considered here with the mean field in the z-direction, only , and . In the collisional case (which was not considered) also .
- •
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 (588). The FLR2 tensor also contains and the Hall-term from the induction equation, see (635). It can also contain gyrotropic heat flux contributions, see (666). The FLR3 contains non-gyrotropic heat flux vectors , , described in Section 5.8, with the detailed algebra presented in Appendix D. It can also contain the non-gyrotropic heat flux tensor , which we neglected. For complete clarity in analyzing the dispersion relations, the entire Hall-CGL-FLR3 model (linearized in the x-z plane, normalized and Fourier transformed) is written down in Section 5.9, see eq. (689)-(700).
- •
Considering the strictly perpendicular fast mode, kinetic theory yields the leading order FLR corrections (in the long-wavelength limit) to the phase speed in the following form . 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 (688). In such a model, the correction to the Alfvén speed is captured correctly, however, the correction to the thermal speed is instead of , i.e., the correction has a wrong sign. This is surprising, since the solution of the FLR2 model, eq. (646), has a correction to the thermal speed , i.e. at least the sign is correct. Finally, the FLR1 model yields solution (645), and does not capture any correction to the Alfvén speed, additionally, the correction to the thermal speed has a 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 colorbars).
- •
In general, when the maximum growth rates are sufficiently large (let’s say ) and not tiny (such as ), 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 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.
14 Acknowledgments
We acknowledge support of the NSF EPSCoR RII-Track-1 Cooperative Agreement No. OIA-1655280 “Connecting the Plasma Universe to Plasma Technology in Alabama”, led by Gary P. Zank. This work was supported by the European Research Council in the frame of the Consolidating Grant ERC-2017-CoG771310-PI2FA “Partial Ionisation: Two-Fluid Approach”, led by Elena Khomenko. Anna Tenerani acknowledges support of the NASA Heliophysics Supporting Research Grant #80NSSC18K1211. PH thanks Thierry Passot, Monica Laurenza, Nikola Vitas, Petr Hellinger and S. Peter Gary for many useful discussions. We are also very thankful to two (out of three) anonymous referees whose comments and suggestions had a great impact on this text. Significant effort has been made to eliminate all the misprints from the equations. However, we will amend possible misprints, if found, in a corrigendum.
Appendix A Fourier transformations
If we have an equation in Real space, we decompose each quantity as a superposition of waves according to
[TABLE]
Most of plasma physics books and kinetic papers use in the definition (A1). Decomposing the frequency to real and imaginary part , so , implies that a wave with grows (i.e. has a positive growth rate) and a wave with is damped. Of course, an alternative definition of (A1) with is allowed, one just has to remember that a wave with is damped and a wave with grows. However, in part 2 of the text, we will calculate the Landau damping and we will see that following all the correct minus signs in the Landau integral can be very confusing. We therefore recommend to use the first decomposition, where waves with are damped. This choice yields that Fourier transformations in the x-z plane (with in the z-direction) are performed according to a shortcut
[TABLE]
The transformation (A1) is technically the inverse/backward Fourier transform . The forward Fourier transform reads
[TABLE]
The location of normalization constants is an ad-hoc choice, one just needs to be consistent in using them, especially when calculating convolutions.
Appendix B Generalized vector (cross) product
In the collisionless fluid hierarchy one encounters a vector product between vectors and tensors, as for example in the pressure tensor equation (26). It is useful to clarify the definition of this generalized vector product. Consider matrix , which is written as a tensor product between vectors , so , where as everywhere in this text, we omit writing the tensor product . Or in the index notation . Now consider vector , and how can one define a generalized vector product , by using the usual vector product between two vectors. The natural definition is , or in the index notation Now let’s consider how to define the much more confusing . The natural definition is , and in the index notation Now it is straightforward to show that
[TABLE]
Generalizing the vector product to higher order tensors is easy and in the fluid hierarchy we always encounter a vector , and tensors , let’s therefore write
[TABLE]
The vector product is very useful for decomposing quantities to the directions parallel and perpendicular with respect to . For example, the usual decomposition of the velocity reads
[TABLE]
Nevertheless, it is possible to directly obtain , by applying twice on , so that
[TABLE]
or . Similar decompositions can be made for tensors. For example, a general matrix can be decomposed according to
[TABLE]
where the perpendicular component can be written as , i.e., where all the act on from the right. Alternatively, one can write a decomposition when all the act on the from the left, , where the perpendicular part .
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
[TABLE]
where . As in advanced fluid models, it is beneficial to normalize the speed with respect to the Alfvén speed , and normalize the length with respect to the ion-inertial length , i.e. by using normalizations (148)-(153), which in Fourier space yields , . Normalizing the MHD equations and dropping the tilde yields unchanged density, induction and pressure equations, and the momentum equation reads
[TABLE]
One can decide how to rewrite the , either by introducing the usual MHD sound speed , or by introducing the plasma beta
[TABLE]
where we use the Boltzmann constant , see the footnote after the definition of thermal speeds (155). By specifying the mean magnetic field to be in the z-direction, the normalized equations are linearized according to
[TABLE]
and without a loss of generality, we consider propagation in the x-z plane (with ). By considering a wave propagating with wavenumber , i.e. a wave propagating in the direction that makes an angle with respect to , the MHD system written in Fourier space reads
[TABLE]
By calculating determinant of this system yields MHD dispersion relation in normalized units (bringing back tildes for clarity) that reads
[TABLE]
or alternatively by using (C3) with the sound speed
[TABLE]
or in physical units
[TABLE]
yielding the MHD dispersion relation for the Alfvén mode (sometimes called the shear Alfvén mode), and the slow and fast modes (259).
In the limit that corresponds to a cold plasma () or a plasma with a very strong magnetic field (), the slow mode , and the fast mode (sometimes called the compressional Alfvén mode). In the limit that corresponds to the incompressible MHD () or a plasma with vanishing magnetic field (), the slow mode , and the fast mode . Note the difference in and in these two limits.
Appendix D Non-gyrotropic heat flux tensor
The heat flux tensor equation (753) is rewritten as
[TABLE]
The “inversion” procedure for the l.h.s. of this is equation is very complicated, and will not be addressed here. Nevertheless, the inversion procedure exist, and an interested reader can check eq. (43) of Ramos (2005). The leading-order (first order in frequency and wavenumber) can be obtained by making the quantities on the r.h.s. gyrotropic, and we want to solve
[TABLE]
The heat flux tensor decomposition is , or alternatively , where , . Since
[TABLE]
one can define perpendicular components of the heat flux vectors
[TABLE]
Therefore, to obtain these vectors, we want to apply , and on the entire equation (D2). A very useful expression also is
[TABLE]
and the trace of the entire heat flux tensor reads
[TABLE]
D.1 Non-gyrotropic heat flux vector
Applying on the l.h.s. of (D2) yields
[TABLE]
Terms on the r.h.s. of (D2) calculate
[TABLE]
For the last term on the r.h.s. of (D2)
[TABLE]
and so for the gyrotropic
[TABLE]
Collecting all the results together yields
[TABLE]
The equation appears very complicated, however, it is not overly so, since it is only a vector equation. The equation just describes perpendicular components of . By multiplying it with , the l.h.s. is zero, and on the r.h.s. one recovers the gyrotropic parallel heat flux equation (822). To understand the equation better, let’s simplify for a moment and evaluate it with respect to . The x-component , so that
[TABLE]
and similarly for the y-component , yielding
[TABLE]
Now, by using the leading order MHD induction equation and (which was for example used to get the FLR1 pressure tensor), and by further considering specific example of a bi-Maxwellian expressions for the 4th-order moments, i.e. etc., yields
[TABLE]
The above equations are very useful if we ever need linearized expressions for , which is beneficial to do right now, and the partially linearized expressions read
[TABLE]
Note that the mean values of perturbations . It is also noteworthy to point out, that the contributions from the induction equation are completely eliminated by the linearization process, and more elaborate forms of induction equation will not bring additional precision at the linear level. Coming back to the expressions (D22), (D23), it is also possible to introduce vorticity , with components and . By noting that and , both equations can be written together as
[TABLE]
Now that we understand that the equation (D19) just describes perpendicular components of with respect to , since the parallel components are zero (directly from the decomposition of )
[TABLE]
we actually do not have to evaluate the eq. (D19) with respect to . For a general vector , the following identity holds . So for a vector which does not have any parallel components to , so that , the identity is , and
[TABLE]
We can apply on the entire equation (D19), and get the fully nonlinear expression for . Let’s apply on each term separately, since we will need these results later, the terms calculate
[TABLE]
To summarize, by applying and on the equation (D2), yields the following nonlinear expression
[TABLE]
Now the calculations proceed in a similar way as before. The usual (MHD) induction equation can be written as
[TABLE]
and since , further yielding
[TABLE]
Finally, by prescribing bi-Maxwellian perturbations for the 4th-order moment yields
[TABLE]
We can report a sign typo in eq. (3.7) of Sulem & Passot (2015), whose term has an opposite sign.
D.2 Conversion of other notations
In the paper of Ramos (2005), the following notation is used
[TABLE]
and his eq. 54, 57 are indeed equivalent to (D36).
In the paper by Macmahon (1965), the following notation is used
[TABLE]
and his eq. 14 is equivalent to (D35), if in Macmahon (1965) we continue the calculation of
[TABLE]
D.3 Non-gyrotropic heat flux vector
Now we want to apply on the equation (D2). Additionally, we now know that has only perpendicular components, and that afterward we will be also applying . The l.h.s. calculates trivially
[TABLE]
and applying therefore yields
[TABLE]
Let’s continue with applying trace on the equation (D2) term by term
[TABLE]
and by further applying on these terms yields
[TABLE]
and direct subtraction with results (D29)-(D33) gives
[TABLE]
To summarize, by applying and on the equation (D2) yields the following nonlinear expression
[TABLE]
By using the simple MHD induction equation gives
[TABLE]
and finally prescribing bi-Maxwellian values for the 4-th order moments yields
[TABLE]
Additionally, evaluation from linear kinetic theory in the gyrotropic limit (at long wavelengths) yields , as we will see in Part 2. With this requirement, the equation (D58) is equivalent to eq. 3.6 of Sulem & Passot (2015). The result (D58) is also equivalent to expressions 55, 58 of Ramos (2005). However, eq. 15 of Macmahon (1965) appears to have a typo, where in comparison to our (D57), his definition of appears to be missing a factor of 2 in front of the term. Ramos (2005) however states that his results are equivalent to those of Macmahon (1965), so we might not correctly understand the notation in that paper.
Partial linearization of result (D58) yields that by components
[TABLE]
where we have kept the contributions, just in case we need them in the future. Note that similarly to the linearization of the vector , the induction equation is completely eliminated at the linear level, and more elaborate forms of induction equation will not have additional contributions.
D.4 2nd-order heat flux vectors
In a previous section we have seen that by applying , and on the equation (D1), the fully nonlinear (but implicit) expressions for the non-gyrotropic heat flux vectors read
[TABLE]
Here we are interested only in the linear level contributions, which we want to use in the FLR pressure tensor. Therefore, since we know that here we will linearize everything at the end, we can immediately get rid of the two heat flux terms and , and also replace the by , so partial linearization yields
[TABLE]
Additionally, the term can be easily precalculated as
[TABLE]
and by applying , terms proportional to disappear, yielding
[TABLE]
The first-order vectors were obtained by keeping only gyrotropic quantities on the right hand side. Linearizing results (D36), (D58) indeed yields that at the linear level
[TABLE]
where the gradients are meant to be further linearized, and where we also neglected perturbations , since right now we do not want to calculate these quantities from linear kinetic theory. To obtain the second-order heat flux vectors at the linear level, we need to calculate
[TABLE]
To calculate the contributions it is enough to work with partially linearized (valid only for bi-Maxwellian)
[TABLE]
which further yields
[TABLE]
and further evaluation with (e.g. , or , but keeping the intact for now) yields
[TABLE]
We purposely kept the for now, since multiplying the results by and further linearizing yields
[TABLE]
which verifies that the contributions (820), (821) to the gyrotropic heat flux equations were indeed calculated correctly. In contrast, applying eliminates terms proportional to , implying
[TABLE]
where we have used perhaps a bit strange mixture of vector and index notations, but which should be otherwise clear. For any vector , the notation just means . Similarly straightforward calculation yields
[TABLE]
The final form of the second-order non-gyrotropic heat flux vectors reads
[TABLE]
Alternatively, one can write . The second order heat flux vectors contribute to scalar (gyrotropic) pressure equations, and also to the FLR pressure tensor . The non-gyrotropic pressure tensor , and the non-gyrotropic heat flux vectors , are therefore generally coupled. To make the system easily solvable, terms with on the r.h.s. of (D86), (D87) are approximated by . The equations (D86), (D87) can be written in a more elegant vector form, see eq. (682), (683), and the result is equivalent to equations 53, 54 of Passot et al. (2012).
Instead of splitting of , and to the first and second order components, one might be interested in the future to check the behaviour of the coupled system. I.e., one might be interested in solving the pressure tensor equation, that is coupled to the heat flux vectors, by making the associated variables as independent quantities. For example, one might be interested to double check the solutions for the parallel firehose instability. The procedure how to obtain time dependent equations for and can be perhaps a bit blurry from the steps outlined above, so let’s write it one more time, only for . One can work directly with the heat flux tensor equation (753), and keeping only terms that contribute at the linear level, the partially linearized heat flux tensor equation reads
[TABLE]
By applying and , yields that at the linear level
[TABLE]
or in the index notation
[TABLE]
Importantly, bi-Maxwellian distribution function was assumed. Normalizing the equations, writing them in the x-z plane and Fourier transforming yields
[TABLE]
Considering parallel propagation (), these equations are coupled with equations for , eq. (686), (687), and also with equations for , system (706). Now it should be easy to derive the dispersion relation, and for example check, if the solutions for the parallel firehose instability are improved.
D.5 Contributions of to various equations
Here we calculate contributions of to the FLR pressure tensor, scalar heat flux equations, and scalar pressure equations. Contributions to the FLR pressure tensor are calculated at the linear level by neglecting . Contributions to scalar heat flux equations and scalar pressure equations are calculated nonlinearly. In these calculations, no specific form of a distribution function is assumed.
Contributions to FLR pressure tensor
For contributions to the FLR tensor , we need to calculate . To calculate it, we will neglect the heat flux tensor , and use the decomposition of the heat flux tensor , that can be figured out to be
[TABLE]
or in the index notation
[TABLE]
Here we want to calculate only contributions at the linear level, and linearizing the yields
[TABLE]
and straightforward evaluation gives
[TABLE]
Contributions to scalar heat flux equations
When we calculated the scalar heat flux equations, we separated the non-gyrotropic heat flux by defining a quantity
[TABLE]
and we need to calculate
[TABLE]
Here we want to obtain exact nonlinear expressions. Direct calculations yields
[TABLE]
and the final fully nonlinear result reads
[TABLE]
Similar calculations yield nonlinear expression for the , which calculates
[TABLE]
To continue the calculation, we need to use (an exact expression)
[TABLE]
which allows us to calculate
[TABLE]
and therefore
[TABLE]
The final fully nonlinear result reads
[TABLE]
However, it is easy to see that at the linear level both and .
Contributions to scalar pressure equations
The non-gyrotropic heat flux term entering the parallel pressure equation calculate
[TABLE]
and since
[TABLE]
which yields the final nonlinear result
[TABLE]
Similarly, the term entering the perpendicular pressure equation calculates
[TABLE]
The nonlinear results (D120), (D122) are exact, and here we do not address how to further decompose the tensor . Finally, considering contributions at the linear level (assuming that the mean values of the entire heat flux ), it is easy to see that
[TABLE]
The heat flux tensor therefore does not contribute to the gyrotropic pressure equations at the linear level.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abraham-Shrauner (1967) Abraham-Shrauner, B. 1967 Propagation of hydromagnetic waves through an anisotropic plasma. J. Plasma Physics 1 , 361–378.
- 2Abraham-Shrauner (1973) Abraham-Shrauner, B. W. 1973 Small amplitude hydromagnetic waves for a plasma with a generalized polytrope law. Plasma Physics 15 , 375.
- 3Bale et al. (2009) Bale, S. D., Kasper, J. C., Howes, G. G., Quataert, E., Salem, C. & Sundkvist, D. 2009 Magnetic Fluctuation Power Near Proton Temperature Anisotropy Instability Thresholds in the Solar Wind. Phys. Rev. Lett. 103 , 211101.
- 4Barakat & Schunk (1982) Barakat, A. R. & Schunk, R. W. 1982 Transport equations for multicomponent anisotropic space plasmas - a review. Plasma Physics 24 , 389–418.
- 5Basu (2009) Basu, B. 2009 Hydromagnetic waves and instabilities in kappa distribution plasma. Phys. Plasmas 16 (5), 052106.
- 6Basu & Coppi (1984) Basu, B. & Coppi, B. 1984 Theory of field-swelling instability in anisotropic plasmas. Physics of Fluids 27 , 1187–1193.
- 7Belmont & Mazelle (1992) Belmont, G. & Mazelle, C. 1992 Polytropic indices in collisionless plasmas - Theory and measurements. J. Geophys. Res. 97 , 8327–8336.
- 8Bhatnagar et al. (1954) Bhatnagar, P. L., Gross, E. P. & Krook, M. 1954 A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Physical Review 94 , 511–525.