In this paper, we investigate the kinetic stability of classical, collisional plasma – that is, plasma in which the mean-free-path $\lambda$ of constituent particles is short compared to the length scale $L$ over which fields and bulk motions in the plasma vary macroscopically, and the collision time is short compared to the
evolution time. Fluid equations are typically used to describe such plasmas, since their distribution functions are close to being Maxwellian. The small deviations from the Maxwellian distribution are calculated via the Chapman-Enskog (CE) expansion in $\lambda/L \ll 1$, and determine macroscopic momentum and heat fluxes in the plasma. Such a calculation is only valid if the underlying CE distribution function is stable at collisionless length scales and/or time scales. We find that at sufficiently high plasma $\beta$, the CE distribution function can be subject to numerous microinstabilities across a wide range of scales. For a particular form of the CE distribution function arising in strongly magnetised plasma (viz., plasma in which the Larmor periods of particles are much smaller than collision times), we provide a detailed analytic characterisation of all significant microinstabilities, including peak growth rates and their associated wavenumbers. Of specific note is the discovery of several new microinstabilities, including one at sub-electron-Larmor scales (the `whisper instability') whose growth rate in certain parameter regimes is large compared to other instabilities. Our approach enables us to construct the kinetic stability maps of classical, two-species collisional plasma in terms of $\lambda$, the electron inertial scale $d_e$ and the plasma $\beta$. This work is of general consequence in emphasising the fact that high-$\beta$ collisional plasmas can be kinetically unstable; for strongly magnetised CE plasmas, the condition for instability is $\beta \gtrsim L/\lambda$. In this situation, the determination of transport coefficients via the standard CE approach is not valid.

§ INTRODUCTION

Answering the question of when a plasma can be described adequately
by fluid equations is fundamental for a comprehensive understanding of plasma dynamics.
It is well known that some physical effects in plasmas – for example, Landau damping –
specifically require a fully kinetic description in terms of distribution functions of the plasma's constituent particles [Landau, 1946].
However, for many other plasma processes, a detailed description of the underlying
particle distribution provides little additional understanding of the essential physics governing that process.
Characterising such processes with fluid equations, which describe the evolution of macroscopic physical quantities
such as density, fluid velocity and temperature,
often simplifies the description and therefore aids understanding.
Fluid equations are also easier to solve numerically than kinetic equations:
the latter reside in six-dimensional phase space (and time), with three
additional dimensions – the velocity space – when compared to the former.
The underlying difficulty associated with determining when a plasma is a fluid
is finding a closed set of equations in the macroscopic plasma variables. The derivation of fluid equations from the Maxwell-Vlasov-Landau equations
governing the evolution of the plasma's distribution functions is carried out by taking moments
(that is, integrating the governing equations and their outer products with velocity $\boldsymbol{v}$ over velocity space). However, the resulting equations
are not closed: the evolution equation of the zeroth-order moment (density) requires
knowledge of the evolution of the first-order moment, the evolution equation for the first-order moment
needs the second-order moment, and so on. For plasma fluid equations to be
able to describe the evolution of a plasma without reference to that plasma's
underlying distribution functions, a closure hypothesis or an approximation relating higher-order
moments to lower ones is required.

For a collisional plasma – i.e., one in which the mean free paths $\lambda_{s}$
and collision times $\tau_s$ of the ions and electrons ($s = i, e$) are much smaller than
the typical length scale $L$ and time scale $\tau_L$ on which macroscopic properties of the plasma
change – there is a procedure for achieving such a closure: the Chapman-Enskog (CE)
expansion [Chapman & Cowling, 1970, Enskog, 1917, Cercignani, 1988]. It is assumed that in a collisional plasma,
the small perturbations of the distribution functions away from a Maxwellian equilibrium have typical size $\epsilon \sim \lambda_{s}/L \sim \tau_s/\tau_L \ll 1$ (assuming sonic motions, and $\lambda_i \sim \lambda_e$).
Since the perturbation is small, its form can be determined
explicitly by performing an asymptotic expansion of the Maxwell-Vlasov-Landau equations.
Once the underlying distribution is known, the relevant moments can be
calculated – in particular, the momentum and heat fluxes are the second- and third-order moments of the $\textit{O}(\epsilon)$ non-Maxwellian component of the distribution
The CE expansion applied to a two-species magnetised plasma was worked out by Braginskii, 1965.
Subsequent studies have refined and extended various aspects of his calculation [Epperlein, 1984, Mikhailovskii & Tsypin, 1984, Epperlein & Haines, 1986, Helander et al., 1994, Simakov & Catto, 2004].
In this paper, we will refer to the distribution functions associated with the CE
expansion as CE distribution functions, and plasmas
with particle distribution functions given by CE distribution functions as CE

However, the theory constructed as outlined above is incomplete. For the CE
expansion to provide an adequate fluid closure, the resulting distribution functions
must be stable to all kinetic instabilities with
length scales shorter than the longest mean free path, and timescales shorter than the
macroscopic plasma timescale $\tau_L$. Such instabilities (if present) are known as
microinstabilities. We emphasise that these microinstabilities should be distinguished
conceptually from instabilities describable by the closed set of plasma-fluid equations:
for example, Rayleigh-Taylor [Rayleigh, 1883, Taylor, 1950, Takabe et al., 1985, Kull, 1991], magnetorotational [Balbus & Hawley, 1991, Hawley & Balbus, 1991], magnetoviscous [Quataert et al., 2002, Balbus, 2004, Islam & Balbus, 2005], or magnetothermal/heat-flux-driven buoyancy instabilities [Balbus, 2000, Balbus, 2001, Quataert, 2008, Kunz, 2011]. Kinetic microinstabilities should also be distinguished
from the small-scale instabilities that arise in solving higher-order $(\textit{O}(\epsilon^2)$) fluid equations obtained
from the CE asymptotic expansion <cit.>.
Such instabilities are not physical because they arise at scales
where the equations themselves do not apply [Bobylev, 1982].
Fluid instabilities do not call into question the
validity of the fluid equations themselves; in contrast, if
microinstabilities occur, the plasma-fluid equations obtained through the closure
hypothesis are physically invalid, irrespective of their own stability.

Microinstabilities have been studied in depth for a
wide range of classical plasmas by many authors; see, for example, Davidson, 1983, Gary, 1993, and Hasegawa, 2012 for three different general perspectives on microinstability theory.  Although it can be shown that a Maxwellian distribution is always immune to such instabilities [Bernstein, 1958, Krall & Trivelpiece, 1973],
anisotropic distribution functions are often not [Furth, 1963, Kahn, 1962, Kalman et al., 1968]. A notable example is the Weibel instability,
which occurs in counter-streaming unmagnetised plasmas [Weibel, 1959, Fried, 1959]. The linear theory of such instabilities is generally well known <cit.>.  Microinstabilities in magnetised plasma have also been comprehensively studied.
The ion firehose and mirror instabilities
are known to occur in plasmas with sufficient ion-pressure anisotropy and large enough plasma $\beta$ [Chandrasekhar et al., 1958, Parker, 1958, Vedenov & Sagdeev, 1958, Hasegawa, 1969, Hall, 1981, Hellinger, 2007],
while electron-pressure anisotropy can also result in microinstabilities of various types [Kennel & Petschek, 1966, Hollweg & Völk, 1970, Gary & Madland, 1985].

A number of authors have noted that microinstabilities, if present, will have a significant effect on the macroscopic transport properties of plasmas [Kahn, 1964, Schekochihin et al., 2005, Schekochihin et al., 2008, Melville et al., 2016, Riquelme et al., 2016, Komarov et al., 2016, Komarov et al., 2018, Roberg-Clark et al., 2018, Drake et al., 2021].
Typically (although not always), once the small-scale magnetic and electric fields
associated with microinstabilities have grown, they will start to scatter particles,
which in turn will alter the plasma's distribution functions. This has micro- and macroscopic consequences for plasma behaviour.
From the microscopic perspective, it changes the course of the evolution of the microinstabilities
themselves – by, e.g., reducing the anisotropy of the underlying particle distribution functions [Hellinger et al., 2014, Riquelme et al., 2018].
From the macroscopic perspective, the changes to the distribution functions will alter both heat and
momentum fluxes in the plasma (which, as previously mentioned, are determined by non-Maxwellian terms in the distribution
function). In this picture, a plasma subject to microinstabilities in some sense generates its own
effective anomalous collisionality [Schekochihin et al., 2008, Mogavero & Schekochihin, 2014, Kunz et al., 2014, Squire et al., 2017, Kunz et al., 2020].
The typical values of the altered fluxes attained must depend on the saturated state of microinstabilities [Schekochihin et al., 2010].
Exploring the mechanisms leading to saturation of both unmagnetised, Weibel-type instabilities <cit.>
and magnetised instabilities <cit.>
continues to be an active research area. Simulation results [Hellinger et al., 2009, Kunz et al., 2014, Guo et al., 2014, Riquelme et al., 2016, Melville et al., 2016, Guo et al., 2018, Bott et al., 2021]
support the claim that the saturation amplitude of
such microinstabilities is typically such that the plasma maintains itself close to
marginality of the relevant instability.

Do these kinetic instabilities afflict the CE distribution function? Naively, it might
be assumed not, since it is `almost' Maxwellian. However, it turns out that, provided the plasma $\beta$ is sufficiently high, small
distortions from a Maxwellian can be sufficient to lead to instability.
Instabilities of a CE distribution function in an unmagnetised plasma were first explored
by Kahn, 1964, who considered a collisional electron plasma  (mean free path $\lambda_e$)
with macroscopic variations in density, temperature and velocity (scale $\sim$$L$). He showed that
the CE distribution function in such a plasma would have two non-Maxwellian terms of order $\lambda_e/L$ – an antisymmetric term associated
with heat flux, and another term associated with velocity shear – and
that the latter term would result in the so-called transverse instability. Kahn, 1964
also claimed that this instability would lead to a significant change in the
plasma viscosity, and other transport coefficients. Albright, 1970, Albright, 1970 further
developed the theory of the transverse instability, including a quasi-linear
theory resulting in isotropisation of the underlying electron distribution

The stability of the CE distribution function was later considered by Ramani & Laval, 1978. They found that in an initially unmagnetised two-species plasma supporting
a fluid-scale electron-temperature gradient (scale $L_T$, no flow shear), the second-order terms (in $\lambda/L_T$) in the electron distribution function
could result in the formation of unstable waves, with typical real frequencies $\varpi \propto \lambda_{e}/L_T$, and growth rates $\gamma_{\rm RL} \propto
\left(\lambda_{e}/L_T\right)^2$. Similarly to Kahn, 1964, they argued that the presence of such instabilities would
suppress the macroscopic heat flux in the plasma (which in a collisional plasma is carried predominantly by
electrons). This particular instability has also been proposed as an explanation for the origin of the cosmic magnetic field [Okabe & Hattori, 2003].
Subsequent authors have explored further the idea that non-Maxwellian components of the electron distribution function required to support a
macroscopic heat flux can lead to kinetic instability.  Levinson & Eichler, 1992 considered the
effect of introducing a uniform, macroscopic magnetic field into the same problem, and found
that a faster instability feeding off first-order heat-flux terms in the CE distribution function – the whistler instability – arose at the electron Larmor scale, with $\gamma_{\rm whistler,T} \propto \lambda_{e}/L_T$. A
quasi-linear theory of this instability was
subsequently constructed by Pistinner & Eichler, 1998. Both Levinson & Eichler, 1992 and Pistinner & Eichler, 1998 proposed that the instability at saturation would result in a suppressed heat flux <cit.>.
More recently, the whistler instability has been studied in simulations of high-$\beta$ plasma–
with two groups independently finding both the onset of instability at electron
scales, and evidence of a suppression of heat flux [Roberg-Clark et al., 2016, Roberg-Clark et al., 2018, Komarov et al., 2018, Roberg-Clark et al., 2018].
Drake et al., 2021 constructed a theoretical model for whistler-regulated heat transport based on a set of
reasonable assumptions that were motivated by these prior simulations.

The possibility of microinstabilities associated with the ion CE distribution
function was also considered by Schekochihin et al., 2005, who found that
weakly collisional, magnetised plasma undergoing subsonic, turbulent shearing motions
can be linearly unstable to firehose and mirror instabilities at
sufficiently high $\beta_i$ (where $\beta_i$ is the ion plasma beta). This is because the shearing motions give rise
to an ion pressure anisotropy $\Delta_i \sim \lambda_i^2/L_V^2$, where $L_V$ is the length scale associated with the shearing
For $|\Delta_i| \gtrsim \beta_i^{-1}$, the mirror and firehose instability thresholds can
be crossed (the mirror instability is trigged by sufficiently positive pressure anisotropy,
the firehose instability by negative pressure anisotropy).
Beyond its threshold, the maximum firehose instability growth rate $\gamma_{\rm fire}$ was found to
satisfy $\gamma_{\rm fire} \propto |\Delta_i+2/\beta_i|^{1/2}$, whilst for the mirror instability, the maximum
growth rate was $\gamma_{\rm mirr} \propto \Delta_i-1/\beta_i$. Such destabilisation of shearing motions was confirmed numerically by Kunz et al., 2014, followed by many others <cit.>.

In this paper, we examine the criteria for the CE distribution function to be stable to
microinstabilities at collisionless scales – i.e., at $k \lambda_s \gg 1$ (where $k$ is the microinstability wavenumber), and $\gamma \tau_L \gg 1$.
In a two-species plasma with a fixed mass ratio $\mu_e \equiv m_e/m_i$ and a charge $Z$ that is not very large, these criteria turn out to be relationships between
three dimensionless parameters: $\lambda/L$, $d_e/L$, and $\beta$, where $\lambda \equiv \lambda_e =
\lambda_i$ is the mean free path for both ions and electrons, and $d_e$ is the electron inertial scale.
The first criterion (which we refer to as the $\beta$-stabilisation condition)
is that the ratio $\lambda/L$ be much smaller than the reciprocal of the plasma
$\beta$, viz. $\lambda \beta/L \ll 1$. This condition arises because the
microinstabilities discussed in this paper are stabilised (usually by Lorentz forces) at sufficiently low
$\beta$. The second criterion (the collisional-stabilisation condition) is that
the characteristic wavenumber $k_{\rm{peak}}$ of the fastest-growing
microinstability in the absence of collisional effects be comparable to (or smaller than) the
reciprocal of the mean-free-path: $k_{\rm{peak}} \lambda \lesssim 1$.
Unlike the $\beta$-stabilisation condition, we do not justify this
condition rigorously, because our calculations are only valid for wavenumbers $k$ such that $k \lambda \gg
1$; thus, we cannot say anything definitive about the $k \lambda \lesssim 1$
regime. We do, however, show that another, more restrictive
stabilisation condition that one might naively expect to exist on account of collisions – that
microinstabilities cannot occur if their growth rate $\gamma$ is smaller than
the collision frequency (viz., $\gamma \tau_s \lesssim 1$) – does not, in fact, apply
to the most significant microinstabilities in CE plasma.
There are good physical reasons to believe that the CE distribution
function is stable against collisionless microinstabilities if the collisional-stabilisation
condition $k_{\rm{peak}} \lambda \lesssim 1$ is satisfied: not least that the typical growth time of the fastest microinstability in CE plasma (calculated neglecting collisional damping of microinstabilities) becomes comparable to the macroscopic evolution time scale
$\tau_L$.  We thus assume the validity of
the collisional-stabilisation condition throughout this paper.
How $k_{\rm{peak}}$ relates to the other
physical parameters is in general somewhat complicated; however, typically
the collisional-stabilisation condition can be written as a lower bound on the ratio $d_e/L$. For example,
in the limit of very high $\beta$, it is $d_e/L > (m_e/m_i)^{-1/6}
(\lambda/L)^{2/3}$ (see section <ref>).

If both the $\beta$-stabilisation and collisional-stabilisation conditions are violated, we demonstrate that
CE plasma will be subject to at least one microinstability, and quite possibly
multiple microinstabilities across a wide range of scales. Some of these microinstabilities are thresholdless – that is, without including collisional
effects, they will occur for CE distributions
departing from a Maxwellian distribution by an asymptotically small amount. Note that
all significant microinstabilities associated with the CE distribution function are `low frequency':
their growth rate $\gamma$ satisfies $\gamma \ll k v_{\mathrm{th}s}$, where $k$ is
the typical wavenumber of the instability, and $v_{\mathrm{th}s}$ the thermal velocity of
the particles of species $s$. This property enables a small anisotropy of the
distribution function to create forces capable of
driving microinstabilities (see section <ref>).

In this paper, we characterise all significant microinstabilities that
arise at different values of $\lambda/L$, $\beta$, and $d_e/L$ for a particular
form of the CE distribution function appropriate for a strongly magnetised plasma – that is, a plasma
where the Larmor radii of ions and electrons are much smaller than the
corresponding mean free paths of these particles. We treat this particular case because of its importance to astrophysical
systems, which almost always possess macroscopic magnetic fields of sufficient
strength to magnetise their constituent particles [Schekochihin & Cowley, 2006].
Our characterisation of microinstabilities
focuses on providing the maximum microinstability growth rates, as well as the wavenumbers
at which this growth occurs. We find that there exist two general classes of
microinstabilities: those driven by the non-Maxwellian component of the CE
distribution associated with temperature gradients, and those driven by the
non-Maxwellian component associated with bulk velocity gradients (`shear'). We refer to
these two non-Maxwellian terms (which exist for both the ion and electron CE distribution functions)
as the CE temperature-gradient terms and the CE shear terms
respectively. Microinstabilities driven by the CE temperature-gradient terms are
called the CE temperature-gradient-driven (CET) microinstabilities, while those driven by
the CE shear terms are the CE shear-driven (CES) microinstabilities.

As expected, within this general microinstability classification scheme,
we recover a number of previously identified microinstabilities,
including the (electron-shear-driven) transverse instability (which we discuss in sections <ref> and <ref>), the whistler instability (section <ref>), the electron mirror instability (section <ref>),
the electron firehose instability (sections <ref> and <ref>), the
ordinary-mode instability (section <ref>), the
(electron-temperature-gradient-driven) whistler heat-flux instability (sections <ref> and <ref>), and the (ion-shear-driven) mirror (section <ref>) and firehose (sections  <ref>, <ref>, <ref>, <ref>, and <ref>) instabilities.
We also find four microinstabilities that, to our knowledge, have not been previously discovered: two ion-temperature-gradient-driven ones at ion Larmor scales – the
slow-hydromagnetic-wave instability (section <ref>) and the long-wavelength kinetic-Alfvén wave instability (section <ref>)
– and two electron-shear-driven ones – the electron-scale-transition (EST) instability (section <ref>) and the whisper instability (section <ref>)
– at electron-Larmor and sub-electron-Larmor scales, respectively.
Of these microinstabilities, the whisper instability seems to be of particular
significance: it has an extremely large growth rate in certain
parameter regimes, and is associated with a new high-$\beta$ wave in a Maxwellian plasma, which also appears to have previously escaped attention.
For convenience, a complete index of microinstabilities discussed in this paper is given in table
<ref>, while the peak growth rates of these microinstabilities
and the scales at which they occur (for a hydrogen CE plasma) are given in table <ref>.

Microinstability name     Section(s)       Other names
occurring in literature       Driving CE term

Mirror instability     <ref>     –     Ion-velocity shear

Firehose instability        <ref>, <ref>,
<ref>, <ref>,
<ref>        Garden-hose
instability        Ion-velocity shear

wave instability*      <ref>     –       Ion-temperature

Long-wavelength kinetic Alfvén
wave (KAW) instability*      <ref>     –       Ion-temperature

CES whistler
instability       <ref>       Electron-cyclotron
(whistler) instability         Electron-velocity

mirror instability       <ref>       KAW, field-swelling
instability        Electron-/ion-
velocity shear

firehose instability       <ref>, <ref>     KAW instability       Electron-/ion-
velocity shear

(EST) instability*      <ref>     –       Electron-velocity

Whisper instability*     <ref>     –       Electron-velocity

Transverse instability     <ref>, <ref>        Small-anisotropy
Weibel instability         Electron-velocity

instability       <ref>     –       Electron-velocity

CET whistler
instability       <ref>, <ref>       Whistler heat
flux instability         Electron-temp.

Index of microinstabilities. The microinstabilities listed here are those discussed in the main text, with the relevant sections indicated.
We also indicate whether
these microinstabilities are driven by macroscopic electron/ion temperature gradients associated with the CE
distribution function, or by macroscopic electron/ion velocity gradients (shears): see section <ref> for a discussion of this classification. Newly identified microinstabilities are indicated with an asterisk.

name      Wavenumber scale       Growth rate
($\times \Omega_e$)      $\beta$-threshold

instability       $k_{\|} \rho_i \lesssim k_{\perp} \rho_i \sim 1 $     $\mu_e \epsilon$     $\epsilon \beta \sim  1$

Parallel firehose
instability      $k_{\|} \rho_i \sim k_{\perp} \rho_i \sim \epsilon^{1/2}$, $ k_{\bot} \ll \epsilon^{1/4} k_{\|}$     $\mu_e \epsilon$     $\epsilon \beta \sim 1 $

Oblique firehose
instability       $ k_{\|} \rho_i \sim \epsilon^{1/4}$,  $k_{\bot} \sim k_{\|}$     $\mu_e \epsilon^{3/4}$     $\epsilon \beta \sim  1 $

Critical-line firehose
instability ($\epsilon \gtrsim 10^{-6}$)      $k_{\|} \rho_i \approx \sqrt{3/2} \, k_{\perp} \rho_i < 1 $     $\mu_e \epsilon^{1/2}$     $\epsilon \beta \sim  1 $

Critical-line firehose
instability ($\epsilon \ll 10^{-6}$)      $k_{\|} \rho_i \approx \sqrt{3/2} \, k_{\perp} \rho_i \sim \epsilon^{1/12} $     $\mu_e \epsilon^{7/12}$     $\epsilon \beta \sim  1 $

instability       $k_{\perp} \rho_i \lesssim k_{\|} \rho_i \sim 1 $     $ \mu_e^{5/4} \epsilon$     $\epsilon \beta \sim \mu_e^{-1/4}$

Long wavelength
KAW instability       $k_{\|} \rho_i < k_{\bot} \rho_i  \sim 1 $     $\mu_e^{5/4} \epsilon k_{\|}/k_{\perp}$     $\epsilon \beta \sim  \mu_e^{-1/4}$

CES whistler
instability       $k_{\bot} \rho_e \lesssim k_{\|} \rho_e \sim 1 $     $\mu_e^{1/2} \epsilon$     $\epsilon \beta \sim  \mu_e^{-1/2}$

Electron mirror
instability       $k_{\|} \rho_e \lesssim k_{\perp} \rho_e \sim 1 $     $\mu_e^{1/2} \epsilon$     $\epsilon \beta \sim  \mu_e^{-1/2}$

Parallel electron
firehose instability       $k_{\bot} \rho_e \lesssim k_{\|} \rho_e \sim 1 $     $\mu_e \epsilon$     $\epsilon \beta \sim  \mu_e^{-1/2}$

Oblique electron
firehose instability       $k_{\|} \rho_e \lesssim k_{\perp} \rho_e \sim 1 $     $\mu_e^{1/2} \epsilon$     $\epsilon \beta \sim  \mu_e^{-1/2}$

EST instability
($\epsilon \beta^{5/7} \lesssim \mu_e^{-1/2}$)       $k_{\|} \rho_e < 1 \lesssim k_{\perp} \rho_e \sim \epsilon^{1/2} \beta^{1/2} \mu_e^{1/4}$     $\mu_e^{5/4} \epsilon^{5/2} \beta^{3/2} $     $\epsilon \beta \sim \mu_e^{-1/2}$

EST instability
($\epsilon \beta^{5/7} \gtrsim \mu_e^{-1/2}$)       $k_{\|} \rho_e < 1 \ll k_{\perp} \rho_e \sim \epsilon^{1/5} \mu_e^{1/10}$     $\mu_e^{1/5} \epsilon^{2/5} $     $\epsilon \beta \sim \mu_e^{-1/2}$

instability       $k_{\|} \rho_e < 1 \ll k_{\perp} \rho_e \sim \epsilon^{1/2} \beta^{1/2} \mu_e^{1/4}$     $\mu_e^{3/8} \epsilon^{3/4} \beta^{1/4} $     $\epsilon \beta \sim \beta^{2/7} \mu_e^{-1/2}$

Parallel transverse
instability      $k_{\perp} \rho_e \lesssim k_{\|} \rho_e \sim \epsilon^{1/2} \beta^{1/2} \mu_e^{1/4}$     $\mu_e^{3/4} \epsilon^{3/2} \beta^{1/2}$     $\epsilon \beta \sim \mu_e^{-1/2}$

Oblique transverse
instability        $1 \lesssim k_{\|} \rho_e \lesssim k_{\bot} \rho_e \sim \epsilon^{1/2} \beta^{1/2} \mu_e^{1/4}$     $\mu_e^{3/4} \epsilon^{3/2} \beta^{1/2}$     $\epsilon \beta \sim \mu_e^{-1/2}$

instability       $k_{\|} = 0$, $k_{\perp} \rho_e \sim \epsilon^{1/2} \beta^{1/2} \mu_e^{1/4}$      $\mu_e^{3/4} \epsilon^{3/2} \beta^{1/2}$      $\epsilon \beta \sim \beta^{2/3} \mu_e^{-1/2} $

CET whistler
instability       $k_{\perp} \rho_e \lesssim k_{\|} \rho_e \sim \epsilon^{1/5} \beta^{1/5} \mu_e^{1/20} $     $\mu_e^{1/4} \epsilon$     $\epsilon \beta \sim  \mu_e^{-1/4}$

Properties of microinstabilities. Typical wavenumbers and maximum growth rates of microinstabilities in strongly magnetised hydrogen CE plasma, and their $\beta$-stabilisation thresholds. Here, $\mu_e = m_e/m_i$.
We assume scalings (<ref>) to relate the magnitude of CE temperature-gradient-driven and CE shear-driven microinstabilities. These scalings lead to the non-Maxwellian component of the ion distribution function having magnitude $\sim$$\epsilon = \mathrm{Ma} \, \lambda/L_V$, where $\lambda = \lambda_e = \lambda_i$, $L_V$ is the length scale of the CE plasma's bulk fluid motions in the direction parallel to the guide magnetic field [see (<ref>d)], and $\mathrm{Ma}$ is the Mach number of those bulk motions.
The quoted wavenumbers and growth rates apply when the $\beta$-stabilisation
threshold is exceeded by an order-unity or much larger factor.

There do exist microinstabilities in CE plasma that are not represented in tables <ref>
and <ref>; however, we claim that the instabilities discussed in this paper are the most significant, on account of their large growth rates and/or low
$\beta$-stabilisation thresholds compared to the unrepresented ones.

Having systematically identified all significant microinstabilities, we can construct `stability maps' of
strongly magnetised CE plasma using “phase diagrams” over a two-dimensional ($\lambda/L$, $d_e/L$)
parameter space at a fixed $\beta$. An example of such a map (for a hydrogen plasma with equal ion and electron temperatures) is shown
in figure <ref>.

Stability map for the CE distribution function. Idealised illustration of the stability
of strongly magnetised, classical, collisional hydrogen plasma to microinstabilities
for different (non-dimensionalised) values of the mean free path $\lambda = \lambda_e = \lambda_i$ and the electron inertial scale $d_e$.
Here, the length scale $L_V$ to which $\lambda$ and $d_e$ are normalised is
the length scale of the CE plasma's bulk fluid motions in the direction parallel to the guide magnetic field [see (<ref>d)];
we assume scalings (<ref>)
to relate the magnitude of CE temperature-gradient-driven and CE shear-driven
microinstabilities, so the CE expansion parameter is $\epsilon = \mathrm{Ma} \, \lambda/L_V$ (see the caption of table <ref> for definitions).
The white region of the ($d_e/L_V$,$\mathrm{Ma} \, \lambda/L_V$) stability map is stable; the coloured
regions are not. In the unstable regions, the fastest-growing microinstability is
indicated by colour according to the figure's legend; in the regions where multiple microinstabilities could be
operating simulataneously, multiple colours have been employed. The plasma beta $\beta$
here was taken to be $\beta = 10^{4}$, and the Mach number $\mathrm{Ma} = 1$.

The entire region of the ($\lambda/L$, $d_e/L$) space depicted in figure <ref> could naively be
characterised as pertaining to classical, collisional plasma, and thus
describable by fluid equations, with transport coefficients given by standard
CE theory. However, there is a significant region of the
parameter space (which is demarcated by boundaries corresponding to the $\beta$-stabilisation and collisional-stabilisation conditions)
that is unstable to microinstabilities.  In fact, in strongly magnetised plasma,
the collisional-stabilisation condition is never satisfied, because there exist
microinstabilities whose characteristic length scales are the ion and
electron Larmor radii, respectively; this being the case, only the $\beta$-stabilisation
condition guarantees kinetic stability.

The effect of microinstabilities being present in CE plasma would be to change the non-Maxwellian components of the
distribution function, and therefore to alter the CE-prescribed resistivity, thermal conductivity
and/or viscosity. Identifying the dominant microinstability or microinstabilities in such plasmas (as is done in figure <ref> for a hydrogen plasma) is
then necessary for calculating the true transport coefficients, which are
likely determined by the effective collisionality associated with the saturated state of the dominant microinstability rather than
by Coulomb collisions. Although such calculations are not undertaken in this paper, it seems
possible that the modified transport coefficients could be determined
self-consistently in terms of macroscopic plasma properties such as temperature gradients or velocity shears.
We note that the calculation presented here assumes that the CE distribution function is determined without
the microinstabilities and thus is only correct when the plasma is stable. Therefore, strictly speaking, the only conclusion one can make when
the CE plasma is unstable is that the naive CE values of transport coefficients should not be taken as

We emphasise that kinetic instability of CE plasmas is a phenomenon of practical importance as well as
academic interest. We illustrate this in tables <ref> and
<ref>, where the possibility of microinstabilities is considered
for a selection of physical systems composed of classical,
collisional plasma.

Environment     $T_{e}$ (eV)     $T_{i}$ (eV)     $n_{e}$ (cm$^{-3}$)    $B$ (G)     $L$ (cm)

Warm intergalatic medium (WIGM)     $10^{2}$     $10^{2}$     $10^{-5}$     $10^{-8}$     $3 \times 10^{24}$

Intracluster medium (ICM)     $10^{4}$     $10^{4}$     $10^{-2}$     $10^{-5}$     $3 \times 10^{23}$

IGM post reionisation      $1$     $1$     $10^{-6}$     $10^{-19}$     $3 \times 10^{24}$

Solar photosphere     $1$     $1$     $10^{17}$     $500$     $10^{7}$

Solar chromosphere     $1$     $1$     $10^{12}$     $10$     $10^{7}$

ICF hot spot (NIF)     $5 \times 10^3$     $5 \times 10^3$     $10^{25}$     $10^7$     $2 \times 10^{-3}$

Laser-ablated plasma (long pulse)     $10^3$     $5 \times 10^2$     $4 \times 10^{21}$     $10^6$     $10^{-2}$

NIF `TDYNO' laser-plasma     $10^3$     $10^3$     $5 \times 10^{20}$     $10^6$     $10^{-2}$

Plasma parameters for some physical systems composed of classical, collisional plasma.
The values of temperature and density of the WIGM given here are from Nicastro et al., 2008, while those of the ICM come from Fabian, 1994. The estimates of the typical magnetic-field strengths and scale lengths for both the WIGM and the ICM are from Ryu et al., 2008. For simplicity, we have assumed equal ion and
electron temperatures; however, we acknowledge that there is some uncertainty as to the validity of this assumption <cit.>. Barkana & Loeb, 2001 is the source of estimates for the IGM post reionisation. Estimates for typical solar parameters are from Wiegelmann et al., 2014
and Stix, 2012. The values of ion temperature and electron density
for ICF hot spots are from Hurricane et al., 2014, who reported DT experiments carried out on the National Ignition
Facility (NIF); the estimates of magnetic-field strength, electron temperature and scale
length come from numerical simulations of the same
experiment [Walsh et al., 2017]. The parameters for the laser-ablated CH plasma are from an
experiment on the OMEGA laser facility, with a 1 ns, 500 J pulse with a $0.351 \, \mu \mathrm{m}$ wavelength [Li et al., 2007];
we assume that the measured fields are found in front of the critical-density surface when
estimating the density. The `TDYNO' laser-plasma is a turbulent CH plasma that was
produced as part of a recent laboratory astrophysics experiment on the NIF which found evidence of suppressed
heat conduction <cit.>.
Naturally, the systems described here often support a range of density, temperatures
and magnetic fields, so the values provided should be
understood as representative, but negotiable.

Environment     $\lambda_e/L$     $\lambda_i/L$     $d_e/L$     $\beta$     $\beta \lambda_e/L$     $\rho_e/\lambda_e$     $\rho_i/\lambda_i$     $k_{\rm{peak}} \lambda_e$

WIGM     $2 \times 10^{-3}$     $2 \times 10^{-3}$     $ 2 \times 10^{-17}$     $10^{4}$     20     $ 10^{-12}$     $10^{-11}$     $10^{12}$

ICM     $ 10^{-2}$      $ 10^{-2}$     $10^{-17}$     $10^{2}$     $1$     $10^{-14}$
$10^{-13}$     $10^{14}$

Reion. IGM     $ 10^{-7}$      $ 10^{-7}$     $10^{-16}$     $10^{22}$     $10^{15}$     $0.5$
$20$     $10^{5}$

Photosphere     $ 6 \times 10^{-12} $      $6 \times 10^{-12}$      $2 \times 10^{-10}$     $30$     $10^{-10}$     $110$
$4 \times 10^{3}$     $10^{-4}$

Chromosphere     $ 2 \times 10^{-7}$      $2 \times 10^{-7}$     $5 \times 10^{-8}$     $1$     $10^{-7}$     $0.2$
$6$     $0.2$

ICF hot spot     $ 0.3$      $ 0.2$     $4 \times 10^{-5}$     $4 \times 10^{6}$     $10^6$     $0.1$
$10$     $10^3$

Laser-abl. pl.     $ 7 \times 10^{-3}$      $ 3 \times 10^{-4}$     $8 \times 10^{-4}$     $200$     $1$     $0.4$
$800$     $2.5$

NIF TDYNO     $ 6 \times 10^{-2}$      $2 \times 10^{-2}$     $2 \times 10^{-3}$     $45$     $2.5$     $0.1$
$200$     $10$

Derived plasma parameters for systems composed of classical, collisional plasma.
All parameters are
calculated using Huba, 1994, except for $k_{\mathrm{peak}} \lambda_e$. This is calculated by
considering all possible instabilities, and then finding the magnitude of $k_{\mathrm{peak}} \lambda_e$
for the fastest-growing instability satisfying $k_{\mathrm{peak}} \lambda_e \gtrsim 1$.  Depending on the values of other
parameters, the fastest-growing instability varies between systems; in the WIGM, ICM, laser-ablation and TDYNO plasmas, the whistler heat-flux
instability is the fastest-growing one, while in the reionised IGM or ICF hot spots, the transverse instability is.

We find that, while there exist some systems where CE plasmas are immune to
microinstabilities – for example, the photosphere and chromosphere – there
are many other astrophysical plasma systems that are likely susceptible to them. Similar considerations apply to a
range of laser plasmas, including plasmas generated in inertial-confinement-fusion
and laboratory-astrophysics
experiments. Indeed, a recent experiment carried out on the National Ignition Facility (NIF) – part of a wider
programme of work exploring magnetic-field amplification in turbulent
laser-plasmas [Tzeferacos et al., 2018, Bott et al., 2021, Bott et al., 2021, Bott et al., 2022] – found evidence for the existence of large-amplitude local temperature fluctuations over a range of scales, a finding that was inconsistent with Spitzer thermal conduction [Meinecke et al., 2022]. This claim was corroborated by MHD simulations (with the code FLASH) of the experiment that modelled thermal conduction either using the Spitzer model, or no explicit thermal conduction model: the latter simulations were found to be much closer to the actual experimental data. Because the plasma created in the NIF experiment is also anticipated by our theory to be susceptible to CE microinstabilities, observations of a discrepancy with CE-derived transport coefficients are tantalising. We note that the idea of microinstabilities emerging in both collisional
astrophysical plasmas and laser plasmas is not a new one: see, e.g. Schekochihin et al., 2005 or Hellinger & Trávníček, 2015
in the former context; in the latter, Epperlein & Bell, 1987 or Bell et al., 2020.
However, to our knowledge there does not exist
a systematic treatment of the general kinetic stability of
CE plasmas. This is the gap that this paper attempts to fill.

This paper has the following structure. In section <ref>, we discuss kinetic and fluid descriptions of
classical plasma. We then describe the CE expansion in collisional plasma:
we work out the CE distribution function arising in a two-species strongly magnetised plasma,
evaluate the friction forces, heat and momentum fluxes necessary to
construct a closed set of plasma-fluid equations, and systematically estimate the size of the non-Maxwellian components of this
distribution. Next, we discuss qualitatively
the existence and nature of microinstabilities potentially arising in CE plasma, before
presenting the methodology that we later use to perform the full linear, kinetic stability
calculation. We provide an overview of this methodology in section <ref>,
and then a much more detailed exposition of it in section <ref>: in particular,
we describe in the latter how a simple form of the dispersion relation for the fastest microinstabilities can be
obtained by considering the low-frequency limit $\gamma \ll k v_{\mathrm{th}s}$ of the hot-plasma dispersion relation,
and how this simplified dispersion relation can be solved analytically.
Readers who are uninterested in the technical details of this calculation are
encouraged to pass over section <ref>; knowledge of its
contents is not a pre-requisite for subsequent sections.
In sections <ref> and <ref>, we construct stability maps (analogous to figure <ref>) showing the parameter ranges in which the
CE distribution function is stable, to CET and CES microinstabilities, respectively. The parameters are $\beta$ and $\lambda/L$, and we construct separate
stability maps for CET and CES microinstabilities in order to take into account the fact that $L$ is in general not the same in the situations where these two types of microinstabilities occur.
In section <ref>, we also discuss the
significant CET microinstabilities that can occur (or not) at different values $\lambda/L$ and
$\beta$, and provide simple analytic characterisations of them; in section <ref>, we do the same for significant CES
microinstabilities. Finally, in section <ref>, we discuss the general implications of these
instabilities for classical, collisional plasmas, and consider future research directions.
Throughout this paper, most lengthy calculations are
exiled to appendices; a glossary of mathematical notation is given in appendix

§ PROBLEM SETUP

§.§ Kinetic versus fluid description of classical plasma

The evolution of classical plasma is most generally described by kinetic theory, via the solution
of Maxwell-Vlasov-Landau equations for the distribution functions
of constituent particles. More specifically, in a kinetic description of a plasma, the distribution function $f_s(\boldsymbol{r},\boldsymbol{v},t)$ of the
particle of species $s$ satisfies

\begin{equation}
\frac{\partial f_{s}}{\partial t} + \boldsymbol{v} \bcdot \bnabla f_{s} + \frac{Z_s e}{m_s} \left(\boldsymbol{E}+ \frac{\boldsymbol{v} \times \boldsymbol{B}}{c}\right) \bcdot \frac{\partial f_s}{\partial \boldsymbol{v}}=
\sum_{s'} \mathfrak{C}(f_s,f_{s'}) , \label{MaxVlasLan}
\end{equation}

where $t$ is time, $\boldsymbol{r}$ spatial position, $\boldsymbol{v}$ the velocity, $e$ the elementary
charge, $Z_s e$ the charge and $m_s$ the mass of species $s$, $\boldsymbol{E}$ the
electric field, $\boldsymbol{B}$ the magnetic field, $c$ the speed of
light, and $\mathfrak{C}(f_s,f_{s'})$ the collision operator for interactions between species $s$ and $s'$. Equation (<ref>) is coupled to
Maxwell's equations:

E     =
4 ∑_s Z_s e ∫d^3 v   f_s,
B     =     0,

×E     =
- 1/c B/t,

×B     =
1/c E/t + 4 /c ∑_s Z_s e ∫d^3 v   v   f_s

Together, (<ref>) and (<ref>) form a closed set of governing equations.

The density $n_s$, bulk fluid velocity $\boldsymbol{V}_{s}$ and
temperature $T_s$ of species $s$ can be formally defined in terms of moments of the
distribution function:

n_s     ≡    ∫d^3 v   f_s,

V_s     ≡    1/n_s ∫d^3 v   v   f_s,

T_s     ≡    1/n_s  ∫d^3 v   1/3 m_s |v-V_s|^2

Governing “fluid” equations are then derived by integrating (<ref>) or outer products of (<ref>)
and the velocity variable $\boldsymbol{v}$ with respect to $\boldsymbol{v}$:

D n_s/D t|_s + n_s ·V_s  = 0 ,

m_s n_s D V_s D t |_s = -p_s - ·π_s +Z_s e n_s (E + V_s ×B/c) + R_s ,

3/2 n_s D T_s D t |_s + p_s ·V_s = - ·q_s - π_s:V_s + 𝒬_s ,

\begin{equation}
\frac{\mathrm{D}}{\mathrm{D} t}\bigg{|}_s \equiv \frac{\p }{\p  t} + \boldsymbol{V}_s \bcdot \bnabla
\end{equation}

is the convective derivative with respect to the fluid motions of species $s$, $p_s$ the pressure, $\boldsymbol{\pi}_s$ the viscosity tensor, and $\boldsymbol{q}_s$ the heat flux of species $s$, $\boldsymbol{R}_s$
the friction force on this species due to collisional interactions with
other species, and $\mathcal{Q}_s$ the
heating rate due to inter-species collisions. The latter quantities are formally defined in terms of
the distribution function as follows:

p_s     ≡    ∫d^3 v   1/3 m_s |v-V_s|^2   f_s = n_s T_s ,

π_s     ≡    - p_s I + ∫d^3 v   m_s (v-V_s) (v-V_s)   f_s ,

q_s     ≡    ∫d^3 v   1/2 m_s |v-V_s|^2 (v-V_s)    f_s,

R_s     ≡    ∑_s' ∫d^3 v   m_s v   ℭ(f_s,f_s'),

𝒬_s     ≡    -R_s V_s + ∑_s' ∫d^3 v   1/2 m_s |v|^2

The distribution function only appears in Maxwell's equations via
its zeroth and first moments; namely, Gauss' law (<ref>a) and the Maxwell-Ampère law (<ref>d) can be
written as

E     =
4 ∑_s Z_s e n_s,
×B     =
1/c E/t + 4 /c ∑_s Z_s e n_s V_s.

Unlike the kinetic description, the fluid equations (<ref>) combined with Maxwell's equations (<ref>b), (<ref>d), (<ref>a), and (<ref>b) are not
a closed system: knowledge of the distribution function, not just of $n_s$, $\boldsymbol{V}_s$ or
$T_s$, is required to calculate momentum and heat fluxes, as well as the friction force or

As discussed in the Introduction, solving fluid equations
as opposed to kinetic equations is advantageous in many cases of interest. Since the dimensionality of the kinetic system
is greater (a six-dimensional phase space vs. three-dimensional position
space), solving the kinetic system introduces both significant
numerical and conceptual complexity. However, the system of fluid equations (<ref>) is only usable
if some type of closure can be introduced to calculate
$ \boldsymbol{\pi}_s$, $\boldsymbol{q}_{s}$, $\boldsymbol{R}_s$ and $\mathcal{Q}_s$
in terms of $n_s$, $\boldsymbol{V}_s$ and $T_s$. For classical plasmas, such a
closure is generally not possible, except in the case of strongly
collisional plasmas.

§.§ The Chapman-Enskog (CE) expansion

§.§.§ The CE distribution functions

For a classical, collisional plasma – i.e., a plasma where the mean free path $\lambda_s$ of particles of species $s$
satisfies $\lambda_s/L \ll 1$ for all $s$, $L$ being the length scale over which the macroscopic properties of the plasma vary – a formal procedure exists for deriving a closed system of fluid equations from a
kinetic description of the plasma.
This procedure is the Chapman-Enskog (CE) expansion, which gives distribution
functions that are close to, but not exactly, Maxwellian.
We call them Chapman-Enskog (CE) distribution functions.
The non-Maxwellian components of the CE distribution functions of particle species $s$ are proportional to $\lambda_s/L$, and
must be present in order to support gradients of $n_s$, $\boldsymbol{V}_s$ and
$T_s$ on $\textit{O}(L)$ length scales, because (<ref>b-e) are all zero for a Maxwellian plasma.

We consider a collisional electron-ion plasma (in which, by definition, $\mu_e \equiv m_e/m_i \ll 1$) with the property that
all constituent particle species are strongly magnetised by the macroscopically varying magnetic field $\boldsymbol{B}$:
that is, the Larmor radius $\rho_s \equiv m_s v_{\mathrm{th}s} c/|Z_s| e |\boldsymbol{B}|$ satisfies $\rho_s \ll \lambda_{s}$ both for the ions and for the electrons (here $v_{\mathrm{th}s} \equiv \sqrt{2 T_s/m_s}$
is the thermal speed of species $s$). Equivalently, a strongly magnetised plasma is one in which the Larmor frequency $\Omega_s \equiv e |Z_s|/m_s c$ satisifies $\Omega_s \tau_s \gg 1$, where $\tau_s$ is the collision time of species $s$. In such a plasma, the macroscopic variation of the fluid moments is locally anisotropic with respect to $\boldsymbol{B}$;
$L$ is the typical length scale of variation in the direction
locally parallel to $\boldsymbol{B}$. It can then be shown that, to first order of the Chapman-Enskog expansion in $\lambda_{s}/L \ll 1$,
and to zeroth order in $\rho_s/\lambda_{s} \ll 1$, the CE distribution
functions of the electrons and ions

f_e(ṽ_e,ṽ_e)     =    n_e/v_the^3 ^3/2 exp(-ṽ_e^2)

×{1+[η_e^T A_e^T(ṽ_e) + η_e^R A_e^R(ṽ_e) + η_e^u A_e^u(ṽ_e) ] ṽ_e

+ ϵ_e C_e(ṽ_e) (ṽ_e^2- ṽ_e^2/2)} ,

f_i(ṽ_i,ṽ_i)     =    n_i/v_thi^3 ^3/2 exp(-ṽ_i^2)

×{1+η_i A_i(ṽ_i) ṽ_i + ϵ_i C_i(ṽ_i) (ṽ_i^2- ṽ_i^2/2 )} .

Let us define the various symbols employed in (<ref>),
before discussing the origin of these expressions and their significance for formulating fluid equations (see section <ref>).

The particle velocity $\boldsymbol{v}$ (with the corresponding speed $v = |\boldsymbol{v}|$) is split into components parallel and perpendicular to the macroscopic magnetic field $\boldsymbol{B} = B \hat{\boldsymbol{z}}$ as $\boldsymbol{v} = v_{\|} \hat{\boldsymbol{z}} + \boldsymbol{v}_{\bot}$, and the perpendicular plane is in turn characterised by two vectors $\hat{\boldsymbol{x}}$ and $\hat{\boldsymbol{y}}$ chosen so that $\left\{\hat{\boldsymbol{x}},\hat{\boldsymbol{y}},\hat{\boldsymbol{z}}\right\}$ is an orthonormal basis. The perpendicular velocity is related to these basis vectors by the gyrophase angle $\phi$:

\begin{equation}
\boldsymbol{v}_{\bot} = v_{\bot} \left(\cos{\phi} \, \hat{\boldsymbol{x}} - \sin{\phi} \, \hat{\boldsymbol{y}} \right) .
\end{equation}

The non-dimensionalised peculiar velocity $\tilde{\boldsymbol{v}}_s$ in the rest frame of the ion fluid is defined by $\tilde{\boldsymbol{v}}_s \equiv (\boldsymbol{v}-\boldsymbol{V}_i)/v_{\mathrm{th}s}$, $\tilde{v}_s \equiv
|\tilde{\boldsymbol{v}}_s|$, $\tilde{v}_{s\|} \equiv \hat{\boldsymbol{z}} \bcdot
\tilde{\boldsymbol{v}}_s$, and $\tilde{v}_{s\bot} \equiv |\tilde{\boldsymbol{v}}_s-\tilde{v}_{s\|}
\hat{\boldsymbol{z}}|$. The number densities satisfy the quasi-neutrality condition

\begin{equation}
Z n_i = n_e \, ,
\end{equation}

where we have utilised $Z_{e} = -1$, and defined $Z \equiv Z_i$. We emphasise that $n_{s}$,
and $v_{\mathrm{th}s}$ all vary over length scales $L$ in the plasma, but
not on shorter scales (at least not in the direction locally parallel to $\boldsymbol{B}$). The functions $A_e^T(\tilde{v}_e)$, $A_e^R(\tilde{v}_e)$, $A_e^u(\tilde{v}_e)$, $C_e(\tilde{v}_e)$, $A_i(\tilde{v}_i)$ and
$C_i(\tilde{v}_i)$ are isotropic functions. Their magnitude is $\textit{O}(1)$ when $\tilde{v}_e \sim
1$ or $\tilde{v}_i \sim 1$, for electrons and ions respectively. Finally, the parameters $\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and $\epsilon_i$
are defined as follows:

η_e^T     =     λ_e ∇_ logT_e = sgn(∇_ logT_e) λ_e/L_T  ,

η_e^R     =      λ_e R_e/p_e ,

η_e^u     =      λ_e m_e u_ei/T_e τ_e ,

η_i     =     λ_i ∇_ logT_i = sgn(∇_ logT_i) λ_i/L_T_i,

ϵ_e     =
λ_e/v_the (ẑ ẑ - 1/3I )  :  W_e = sgn[(ẑ ẑ - 1/3I )  :  W_e] V_e/v_the λ_e/L_V_e ,

ϵ_i     =
λ_i/v_thi (ẑ ẑ - 1/3I )  :  W_iv_thi = sgn[(ẑ ẑ - 1/3I )  :  W_i]V_i/v_thi λ_i/L_V ,

where $\lambda_{e}$ is the electron mean free path, $\lambda_{i}$ the ion mean-free-path, $\tau_e$ the electron collision time,  $R_{e\|} \equiv \hat{\boldsymbol{z}} \bcdot \boldsymbol{R}_{e}$ the parallel electron friction force, $\boldsymbol{u}_{ei} \equiv \boldsymbol{V}_e - \boldsymbol{V}_i$
the relative electron-ion velocity, $u_{ei\|} \equiv \hat{\boldsymbol{z}} \bcdot \boldsymbol{u}_{ei}$,

\begin{equation}
\mathsfbi{W}_s =
\bnabla \boldsymbol{V}_s + \left(\bnabla \boldsymbol{V}_s\right)^{T} -\frac{2}{3}  \left(\bnabla \bcdot \boldsymbol{V}_s\right) \mathsfbi{I}
\label{ChapEnsTerms}
\end{equation}

the traceless rate-of-strain tensor of species $s$, $V_e$ ($V_i$) the bulk electron-(ion-)fluid speed, and

L_T     ≡    |∇_ logT_e |^-1   ,

L_T_i     ≡    |∇_ logT_i |^-1   ,

L_V_e     ≡    1/V_e|(ẑ ẑ - 1/3I )  :  W_e |^-1   ,

L_V     ≡    1/V_i |(ẑ ẑ - 1/3I )  :  W_i |^-1   ,

are, respectively, the electron- and ion-temperature and the electron- and ion-flow length scales parallel to the
background magnetic field.
The mean free paths are formally defined
for a two-species plasma by

λ_e     ≡    v_theτ_e ,

λ_i     ≡    v_thiτ_i ,

and the collision times $\tau_e$ and $\tau_i$ are given in terms of macroscopic plasma parameters by

τ_e     ≡    3 m_e^1/2 T_e^3/2/4 √(2 ) Z_i^2 e^4 n_i logΛ_CL ,

τ_i     ≡     3 m_i^1/2 T_i^3/2/4 √(2 ) Z_i^4 e^4 n_i logΛ_CL

where $\log{\Lambda_{\mathrm{CL}}}$ is the Coulomb logarithm [Braginskii, 1965][Braginskii defined his ion collision time as equal to (<ref>b) multiplied by a factor of $\sqrt{2}$; for the sake of species equality, we remove this factor.]. In a
collisional plasma, $\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and $\epsilon_i$
are assumed small. We note that all these parameters can be either positive or
negative, depending on the orientation of temperature and velocity

It is clear from their definitions that each of the non-Maxwellian terms associated with
the parameters $\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and $\epsilon_i$
is linked to a different macroscopic physical quantity. Thus, $\eta_e^T$ and $\eta_i$
are proportional to the electron- and ion-temperature gradients, respectively; we will therefore
refer to the associated non-Maxwellian terms as the CE electron-temperature-gradient term
and the CE ion-temperature-gradient term. We refer to the
non-Maxwellian term proportional to $\eta_e^R$
as the CE electron-friction term, to the non-Maxwellian term proportional to $\eta_e^u$
as the CE electron-ion-drift term, and the non-Maxwellian terms
proportional to $\epsilon_e$ and $\epsilon_i$ as the CE electron-shear term
and the CE ion-shear term. We note that the friction and electron-ion-drift terms appear in the electron CE distribution function
but not the ion CE distribution function because of our choice to define all velocities in the ion-fluid rest frame.

The derivation of the CE distribution functions (<ref>) for a two-species strongly magnetised plasma undergoing sonic motions (that is, $V_i \sim v_{\mathrm{th}i}$) from the kinetic equation (<ref>) was first
completed by Braginskii, 1965 for arbitrary values of $\rho_s/\lambda_{s}$. We do
not reproduce the full derivation in the main text, but, for the reader's convenience, we provide a derivation of $(\ref{ChapEnskogFunc})$ in appendix <ref>.
The gist of the full derivation is to assume that the distribution function is close to a Maxwellian,
with parameters that only evolve on a slow time scale $t' \sim t L/\lambda_e \sim t
L/\lambda_i \gg t$. The kinetic equation (<ref>) is then expanded
and solved order by order in $ \lambda_e/L \sim \lambda_i/L \ll 1$, allowing for
the calculation of the (small) non-Maxwellian components of the distribution
function. The small parameters $\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and
$\epsilon_i$, as well as the isotropic functions $A_e^T(\tilde{v}_e)$, $A_e^R(\tilde{v}_e)$, $A_e^u(\tilde{v}_e)$, $C_e(\tilde{v}_e)$, $A_i(\tilde{v}_i)$ and
$C_i(\tilde{v}_i)$ emerge during this calculation. The precise forms of these functions depend only on the collision operator assumed in the original
Maxwell-Vlasov-Landau system; in appendix <ref>, we provide a simple illustration
of this, by calculating these isotropic functions explicitly for
Krook [Bhatnagar et al., 1954] and Lorentz collision operators (Appendices <ref> and <ref>, respectively).
For the full Landau collision operator,
the equivalent calculation is more complicated, but can be performed (for example)
by expanding the isotropic functions in Sonine polynomials <cit.>.

§.§.§ Closure of fluid equations (<ref>)

Once the CE distribution function has been calculated, the desired fluid closure can be obtained by evaluating
the heat fluxes, the friction forces, and the momentum fluxes (<ref>) associated with the non-Maxwellian
components of the CE distribution functions. These
calculations were carried out in full for arbitrary values of $\rho_s/\lambda_{s}$ by Braginskii, 1965.

We do not reproduce the full fluid
closure relations here; instead, we illustrate how the non-Maxwellian terms
in the CE distribution functions (<ref>) give rise to the friction force and heat fluxes
parallel to the macroscopic magnetic field, as well as to the viscosity tensor. In a strongly
magnetised two-species plasma (where $\rho_s \ll \lambda_s$), parallel friction forces and heat fluxes
are typically much larger than their perpendicular or diamagnetic counterparts.

* Heat fluxes. Recalling (<ref>c), the parallel heat flux $q_{s\|} \equiv \hat{\boldsymbol{z}} \bcdot \boldsymbol{q}_{s}$ associated with
species $s$ is given by

\begin{equation}
q_{s\|} = \frac{1}{2} \int \mathrm{d}^3\boldsymbol{v}_s' \, m_s \left|\boldsymbol{v}_s'\right|^2 v_{s\|}' \,  f_{s}
, \label{par_heat_flux}
\end{equation}

where $\boldsymbol{v}_s' \equiv \boldsymbol{v}-\boldsymbol{V}_s$. Noting that the electron distribution function (<ref>a) is specified in the rest frame of the ions, not electrons,
it is necessary first to calculate the electron distribution function in the electron rest frame
before calculating the parallel electron heat flux. An expression for this quantity is given by (<ref>) in appendix
<ref> as part of our derivation of (<ref>a):

\begin{eqnarray}
f_{e}(v_{e\|}',v_{e\bot}') & = & \frac{n_{e}}{v_{\mathrm{th}e}^3 \upi^{3/2}} \exp \left(-\frac{|\boldsymbol{v}_{e}'|^2}{v_{\mathrm{th}e}^2}\right) \nonumber \\
&& \times \Bigg\{1 + \bigg[\eta_e^T A_e^{T}\!\left(\frac{|\boldsymbol{v}_e'|}{v_{\mathrm{th}e}}\right) +\eta_e^R A_e^{R}\!\left(\frac{|\boldsymbol{v}_e'|}{v_{\mathrm{th}e}}\right) + \eta_e^u \left(A_e^{u}\!\left(\frac{|\boldsymbol{v}_e'|}{v_{\mathrm{th}e}}\right) -1 \right) \bigg] \frac{v_{e\|}'}{v_{\mathrm{th}e}} \nonumber \\
&& \qquad \qquad + \epsilon_e C_e\!\left(\frac{|\boldsymbol{v}_e'|}{v_{\mathrm{th}e}}\right) \left(\frac{v_{e\|}'^2}{v_{\mathrm{th}e}^2} - \frac{v_{e\bot}'^2}{2 v_{\mathrm{th}e}^2} \right) \Bigg\} \, ,\label{CE_dist_electronframe_maintext}
\end{eqnarray}

Now substituting (<ref>)
into (<ref>) (with $s = e$), we find that the
parallel electron heat flux is

\begin{equation}
q_{e\|} = - n_{e} T_e v_{\mathrm{th}e} \left[\mathcal{A}_e^T \eta_e^T + \mathcal{A}_e^R \eta_e^R + \left(\mathcal{A}_e^u - \frac{1}{2}\right) \eta_e^u\right], \label{heatflux_expression}
\end{equation}

\begin{equation}
\mathcal{A}_e^{T,R,u} = -\frac{4}{3\sqrt{\upi}} \int_0^{\infty} \mathrm{d}\tilde{v}_e \, \tilde{v}_e^6 A_e^{T,R,u}(\tilde{v}_e)  \exp \left(-\tilde{v}_e^2\right)
\end{equation}

The minus signs in the definitions of $\mathcal{A}_e^{T,R,u}$
have been introduced so that $\mathcal{A}_e^{T,R,u} \geq 0$ for a typical
collision operator (determining that these constants are indeed positive for any given
collision operator is non-trivial, but it is a simple exercise to show this for a
Krook collision operator, using the expressions for $A_e^T(\tilde{v}_e)$, $A_e^R(\tilde{v}_e)$, and $A_e^u(\tilde{v}_e)$ given in appendix <ref>).
Expression (<ref>) for the electron heat flux can be rewritten as

\begin{equation}
q_{e\|} =  - \kappa_{e}^{\|} \nabla_{\|} T_e - \left[ \mathcal{A}_e^u - \frac{1}{2} - \frac{ \mathcal{A}_e^R}{\tilde{\mathcal{A}}_e^R}\left(\tilde{\mathcal{A}}_e^u-\frac{1}{2}\right)\right]
n_e T_e u_{ei\|} \, , \label{brag_heatflux_elec}
\end{equation}

where the parallel electron heat conductivity is defined by

\begin{equation}
\kappa_{e}^{\|} = 2 \left(\mathcal{A}_e^T - \frac{\mathcal{A}_e^R}{\tilde{\mathcal{A}}_e^R} \tilde{\mathcal{A}}_e^T\right) \frac{n_e T_e \tau_e}{m_e}
\, , \label{parallel_cond_electrons}
\end{equation}

\begin{equation}
\tilde{\mathcal{A}}_e^{T,R,u} = -\frac{4}{3\sqrt{\upi}} \int_0^{\infty} \mathrm{d}\tilde{v}_e \, \tilde{v}_e^4 A_e^{T,R,u}(\tilde{v}_e)  \exp \left(-\tilde{v}_e^2\right) \,
\end{equation}

Numerical evaluation of the coefficients $\mathcal{A}_e^{T,R,u}$ and $\tilde{\mathcal{A}}_e^{T,R,u}$ for the Landau collision operator
gives [Braginskii, 1965]

\begin{equation}
q_{e\|} \simeq  - 3.16 \frac{n_e T_e \tau_e}{m_e}  \nabla_{\|} T_e + 0.71 n_e T_e u_{ei\|} \,
\end{equation}

The ion heat flux can be calculated directly from (<ref>) ($s = i$) using (<ref>b):

\begin{equation}
q_{i\|} =  - n_{i} T_i v_{\mathrm{th}i} \mathcal{A}_i \eta_i , \label{ion_par_heat_flux}
\end{equation}

\begin{equation}
\mathcal{A}_i = -\frac{4}{3\sqrt{\upi}} \int_0^{\infty} \mathrm{d}\tilde{v}_i \, \tilde{v}_i^6 A_i(\tilde{v}_i)  \exp \left(-\tilde{v}_i^2\right) \, .
\end{equation}

This becomes

\begin{equation}
q_{i\|} =   - \kappa_{i}^{\|} \nabla_{\|} T_i , \label{brag_heatflux_ion}
\end{equation}

where the parallel ion heat conductivity is

\begin{equation}
\kappa_{i}^{\|} = -2 \mathcal{A}_i \frac{n_i T_i \tau_i}{m_i} \simeq - 3.9 \frac{n_i T_i \tau_i}{m_i}
\, . \label{parallel_cond_ions}
\end{equation}

The last equality is for the Landau collision operator [Braginskii, 1965].
Note that the absence of a term proportional to the
electron-ion-drift in the ion heat flux (<ref>) is physically due to the smallness of
the ion-electron collision operator [Helander & Sigmar, 2005].

* Friction force.
We evaluate the friction force by considering the electron-ion-drift associated with electron CE distribution function.
Namely, noting that

\begin{equation}
u_{ei\|} = \frac{v_{\mathrm{th}e}^4}{n_e} \int \mathrm{d}^3 \tilde{\boldsymbol{v}}_e \, \tilde{v}_{e\|} f_{e} ,
\end{equation}

it follows from (<ref>a) that

\begin{equation}
u_{ei\|} = v_{\mathrm{th}e} \left(\tilde{\mathcal{A}}_e^T \eta_e^T + \tilde{\mathcal{A}}_e^R \eta_e^R +  \tilde{\mathcal{A}}_e^u \eta_e^u \right) .
\label{electronparticleflux}
\end{equation}

This expression can in turn be used to relate the
parallel electron-friction force $R_{e\|}$, defined in (<ref>d), to electron flows and temperature

\begin{equation}
R_{e\|} = -\left(\frac{2 \tilde{\mathcal{A}}_e^u + 1}{2 \tilde{\mathcal{A}}_e^R}\right) \frac{n_e m_e u_{ei\|}}{\tau_e}
- \frac{\tilde{\mathcal{A}}_e^T}{\tilde{\mathcal{A}}_e^R} n_e \nabla_{\|} T_e \, .
\label{Friction_force}
\end{equation}

Evaluating the coefficients $\tilde{\mathcal{A}}_e^T$, $\tilde{\mathcal{A}}_e^R$ and $\tilde{\mathcal{A}}_e^u$ for the full Landau collision operator, one
finds [Braginskii, 1965]

\begin{equation}
R_{e\|} \simeq -0.51 \frac{n_e m_e u_{ei\|}}{\tau_e}
- 0.71 n_e \nabla_{\|} T_e \, .
\end{equation}

* Viscosity tensor.
For gyrotropic distributions such as the CE distribution functions (<ref>), the viscosity tensor $\boldsymbol{\pi}_s$
of species $s$ defined by (<ref>b)
– which is the momentum flux excluding the convective terms and isotropic pressure  – is given by

\begin{equation}
\boldsymbol{\pi}_s = \left(p_{s\|} - p_{s\perp}\right) \left(\hat{\boldsymbol{z}} \hat{\boldsymbol{z}}-\frac{1}{3}\mathsfbi{I}\right)  ,
\label{viscosity_tens_gyro}
\end{equation}

where the parallel pressure $p_{s\|}$ and the perpendicular pressure $p_{s\bot}$ are
defined by

p_s     ≡    ∫d^3 v_s'   m_s |v_s'|^2 f_s = n_s T_s (1-2/3ϵ_s 𝒞_s)   ,

p_s     ≡    1/2∫d^3 v_s'    m_s |v_s'|^2 f_s = n_s T_s (1+1/3ϵ_s 𝒞_s)

with the last expressions having being obtained on substitution of the CE distribution function (<ref>),

\begin{equation}
\mathcal{C}_s = -\frac{8}{5\sqrt{\upi}} \int_0^{\infty} \mathrm{d}\tilde{v}_s \, \tilde{v}_s^6 C_s(\tilde{v}_s)  \exp \left(-\tilde{v}_s^2\right) \, .
\label{pressureanisop_const}
\end{equation}

The sign of the constant $\mathcal{C}_s$ is again chosen so that $\mathcal{C}_s > 0$ for typical
collision operators; for the Landau collision operator, $\mathcal{C}_e \simeq 1.1$ and $\mathcal{C}_i \simeq 1.44$ [Braginskii, 1965].
We note for reference that the parameter $\epsilon_s$ [see (<ref>e-f)]
has a simple relationship to the pressure anisotropy of species $s$: utilising (<ref>),
one finds

\begin{equation}
\Delta_s \equiv \frac{p_{s\bot}-p_{s\|}}{p_{s}} = \mathcal{C}_s \epsilon_s \, . \label{pressureanisop}
\end{equation}

Using (<ref>), the viscosity tensor (<ref>) can
be written

\begin{equation}
\boldsymbol{\pi}_s = - \frac{\mu_{\mathrm{v}s}}{2} \left(\hat{\boldsymbol{z}} \hat{\boldsymbol{z}} - \frac{1}{3}\mathsfbi{I} \right)  {\left(\hat{\boldsymbol{z}} \hat{\boldsymbol{z}} - \frac{1}{3}\mathsfbi{I} \right) \! \boldsymbol{:} \! \mathsfbi{W}_s}
, \label{brag_viscosity}
\end{equation}

where the dynamic viscosity of species $s$ is

\begin{equation}
\mu_{\mathrm{v}s} \equiv 2 \mathcal{C}_s n_s T_s \tau_s \, . \label{dynam_viscosity}
\end{equation}

* Thermal energy transfer between species. It can be shown that for
the CE distribution functions $(\ref{ChapEnskogFunc})$, the rate of thermal
energy transfer from electrons to ions $\mathcal{Q}_e$ is simply

\begin{equation}
\mathcal{Q}_e = - \boldsymbol{R}_e \bcdot \boldsymbol{u}_{ei}  ,
\end{equation}

while the rate of thermal energy transfer from ions to electrons vanishes: $\mathcal{Q}_i \approx
0$. This is because the ion-electron collision rate is assumed small (by a factor of the mass ratio)
compared to the ion-ion collision rate when deriving
(<ref>b), and is thus neglected. Braginskii, 1965 shows that, in
there is a non-zero (but small) rate of transfer:

\begin{equation}
\mathcal{Q}_i = -\mathcal{Q}_e - \boldsymbol{R}_e \bcdot \boldsymbol{u}_{ei}
= \frac{3 n_e m_e}{m_i \tau_{e}} \left(T_e-T_i\right) \, .
\end{equation}

The time scale on which the ion and electron temperatures equilibrate is the ion-electron temperature equilibration time

\begin{equation}
\tau_{ie}^{\rm eq} \equiv \frac{1}{2}  \mu_e^{-1/2} \tau_{i}
. \label{ion_elec_equil_time}
\end{equation}

In summary, the non-Maxwellian components of the CE distribution function are
essential for a collisional plasma to be able to support fluxes of heat
and momentum. More specifically,  (<ref>) demonstrates that the
electron heat fluxes in a CE plasma are proportional to both temperature gradients and
electron-ion drifts, and are carried by the electron-temperature-gradient,
friction and electron-ion-drift terms of the CE distribution function. In contrast, the ion
heat fluxes (<ref>) are proportional only to ion temperature gradients (and carried by
the CE ion-temperature-gradient term). Momentum fluxes (<ref>) for electrons and ions
are carried by the CE electron- and ion-shear terms, respectively, and are
proportional to components of the rate-of-strain tensor.

§.§.§ Relative size of non-Maxwellian terms in the CE distribution function

In the case of magnetised, two-species plasma satisfying $T_i \sim
T_e$, (<ref>) can be used to estimate the size of
the small parameters $\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and
$\epsilon_i$. Although these parameters are a priori proportional to $\lambda_s/L$ for both ions and electrons, their
precise magnitudes are, in fact, subtly different. Namely, the terms associated
with $\eta_e^T$, $\eta_e^R$, $\eta_e^u$ and $\eta_i$ are gradients of the electron and ion temperatures and electron-ion relative parallel drift velocities, whereas terms associated with $\epsilon_e$ and $\epsilon_i$
involve gradients of the bulk flows [cf. (<ref>)] – and these gradients do not necessarily occur on the same
length scale. Recalling that the (electron) temperature
and the (ion) flow length scales parallel to the macroscopic magnetic field are defined by [cf. (<ref>)]

L_T     =     |∇_ logT_e |^-1   ,

L_V     =     1/V_i|(ẑ ẑ - 1/3I )  :  W_i |^-1   ,

where $\mathsfbi{W}_i$ is the ion rate-of-strain tensor (<ref>), and assuming that $L_{T_i} = \left(\nabla_{\|} \log T_i \right)^{-1} \sim L_{T}$ (an assumption we will check a posteriori),
it follows from (<ref>) that

η_e^T     ∼     λ_e/L_T,

η_e^R     ∼    λ_e R_e/p_e ∼λ_e/L_T ∼η_e^T,

η_e^u     ∼    u_ei/v_the ∼λ_e/L_T ∼η_e^T,

η_i     ∼    λ_i/L_T ∼1/Z^2 η_e^T ,

ϵ_e     ∼     V_i/v_the λ_e/L_V ∼Ma  μ_e^1/2 L_T/L_V η_e^T,

ϵ_i     ∼    V_i/v_thi λ_i/L_V ∼Ma   L_T/Z^2 L_V η_e^T,

where $\mathrm{Ma} \equiv V_i/v_{\mathrm{th}i}$ is the Mach number. Note that, to arrive at (<ref>b), we assumed that $R_{e\|} \sim p_e/L_{T}$ and $u_{ei\|} \sim v_{\mathrm{th}e} \lambda_e/L_T$, justified by (<ref>) and (<ref>), respectively. The relative magnitudes of
$\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and
$\epsilon_i$ therefore depend on the Mach number of the plasma, as well as on
the length scales $L_T$ and $L_V$.

In the work of Braginskii, 1965, who a priori presumes all “fluid” quantities in the plasma to vary on just a single scale $L \sim L_T \sim L_V$,
with sonic ordering $\mathrm{Ma} \lesssim 1$, determining the relative size of these parameters for a hydrogen plasma ($Z = 1$) is simple:

\begin{equation}
\epsilon_e \sim  \mu_e^{1/2} \epsilon_i \ll  \epsilon_i \sim  \eta_i \sim \eta_e^T \sim \eta_e^R \sim \eta_e^u
\, . \label{brag_onescale}
\end{equation}

However, in most interesting applications, this single-scale ordering is incorrect. In a plasma with $\lambda_s/L \ll 1$ under Braginskii's ordering,
motions on many scales naturally arise. The fluid Reynolds
number in such a plasma is given by

\begin{equation}
\mathrm{Re} \equiv \frac{V_0 L_0}{\nu} \, ,
\end{equation}

where $V_0$ is the typical fluid velocity at the scale $L_0$ of driving motions and $\nu \equiv \mu_{\mathrm{v}i}/m_i n_i \sim v_{
\mathrm{th}i} \lambda_i$
is the kinematic viscosity [see (<ref>)]. Typically, this number is large:

\begin{equation}
\mathrm{Re} \sim \frac{V_0}{v_{\mathrm{th}i}} \frac{L_0}{\lambda_{i}} \gtrsim \frac{1}{\epsilon_i} \gg 1 \, ,
\label{Reynolds_num}
\end{equation}

where we have assumed $\mathrm{Ma}_0 \equiv V_0/v_{\mathrm{th}i} \lesssim 1$, in line with Braginskii's sonic ordering. Therefore, such a plasma will naturally become turbulent and exhibit motions across a range of scales. As a consequence, velocity and temperature fluctuations
on the smallest (fluid) scales must be considered, since the associated shears and temperature gradients are the largest.
To estimate $\eta_e^T$, $\eta_e^R$, $\eta_e^u$, $\eta_i$, $\epsilon_e$ and
$\epsilon_i$ accurately, we must determine the magnitude of these gradients.

First, let $\ell_{\nu}$ be the smallest scale on which the velocity varies due to turbulent motions (the Kolmogorov scale), with velocity fluctuations on scales $\ell \ll \ell_{\nu}$ being suppressed by viscous diffusion.
Then it follows that $\mathrm{Re}_{\ell_{\nu}} \sim 1$, where $\mathrm{Re}_\ell \equiv V\!\left(\ell\right) \ell/\nu$
is the scale-dependent Reynolds number and $V\!\left(\ell\right)$ is the typical fluid velocity on scale
For Kolmogorov

\begin{equation}
\frac{V\!\left(\ell\right)}{V_0} \sim \left(\frac{\ell}{L_0}\right)^{1/3} \sim \left(\frac{\mathrm{Re}_\ell}{\mathrm{Re}} \right)^{1/4},
\end{equation}

and $\ell/L_0 \sim \left({\mathrm{Re}_\ell}/{\mathrm{Re}} \right)^{3/4}$, which gives $V\!\left(\ell\right)\!/\ell \sim (V_0/L_0) \left({\mathrm{Re}_\ell}/{\mathrm{Re}}
\right)^{-1/2}$, and thus, from (<ref>),

\begin{equation}
\frac{V\!\left(\ell_\nu\right)}{\ell_\nu} \sim  \frac{V_0}{L_0} \left(\frac{\mathrm{Re}_{\ell_\nu}}{\mathrm{Re}}
\right)^{-1/2} \sim \mathrm{Ma}_0^{1/2} \left(\frac{\lambda_{i}}{L_0}\right)^{-1/2}
\frac{V_0}{L_0} \, .
\end{equation}

We therefore conclude that

\begin{equation}
L_V \sim \ell_{\nu} \frac{V_0}{V\!\left(\ell_\nu\right)} \sim L_0 \mathrm{Ma}_0^{-1/2} \left(\frac{\lambda_{i}}{L_0}\right)^{1/2}
\end{equation}

Next, the smallest scale on which the electron temperature varies, $\ell_{\chi}$,
is the scale below which temperature fluctuations are suppressed by
thermal diffusion; it satisfies $\mathrm{Pe}_{\ell_{\chi}} \sim 1$, where $\mathrm{Pe}_\ell \equiv V\!\left(\ell\right) L/\chi$
is the scale-dependent Péclet number and $\chi \equiv 2 \kappa_{e}^{\|}/3 n_e \sim v_{\mathrm{th}e} \lambda_{e}$ is the (parallel) thermal diffusivity [see (<ref>)].
Because temperature is passively advected by the flow, the temperature fluctuation $T(\ell)$ at any scale $\ell > \ell_\chi$
obeys the same scaling as the bulk velocity:

\begin{equation}
\frac{T\!\left(\ell\right)}{T(L_0)} \sim \frac{V\!\left(\ell\right)}{V_0} \sim \left(\frac{\mathrm{Pe}_\ell}{\mathrm{Pe}} \right)^{1/4}  \, .
\end{equation}

In addition, the magnitude of temperature fluctuations at the driving scale is
related to the mean temperature by the Mach number of the driving-scale motions, $T(L_0) \sim T_0 \mathrm{Ma}_0$, which then gives

\begin{equation}
\frac{T\!\left(\ell\right)}{T_0} \sim \mathrm{Ma}_0 \left(\frac{\mathrm{Pe}_\ell}{\mathrm{Pe}} \right)^{1/4}  \, ,
\end{equation}

where $\mathrm{Pe} \equiv \mathrm{Pe}_{L_0}$. It follows from an analogous argument to that just given for the velocity fluctuations that

\begin{equation}
\frac{T\!\left(\ell_\chi\right)}{\ell_\chi} \sim  \frac{T_0}{L_0} \mathrm{Ma}_0 \mathrm{Pe}^{1/2}
\end{equation}

Under Braginskii's ordering, the Prandtl number of CE
plasma is

\begin{equation}
\mathrm{Pr} \equiv \frac{\nu}{\chi} = \frac{\mathrm{Pe}}{\mathrm{Re}} \sim \frac{v_{\mathrm{th}i} \lambda_{i}}{v_{\mathrm{th}e} \lambda_{e} }
\sim \mu_e^{1/2} \ll 1 \, ,
\end{equation}

and, therefore,

\begin{equation}
L_T \sim \ell_{\chi} \frac{T_0}{T\!\left(\ell_\chi\right)} \sim L_0 \mu_e^{-1/4} \mathrm{Ma}_0^{-3/2} \left(\frac{\lambda_{i}}{L_0}\right)^{1/2}
\end{equation}

Thus, $L_V \sim \mathrm{Ma}_0 \mu_e^{1/4} L_T \ll L_T$ under the assumed ordering.

Finally, we consider whether our a priori assumption that $L_{T_i} \sim L_T$ is, in fact, justified. A sufficient condition for ion-temperature gradients
to be the same as electron-temperature gradients is for the evolution time $\tau_{L}$ of all macroscopic motions to be much longer than the ion-electron temperature equilibration time
$\tau_{ie}^{\rm eq}$ defined by (<ref>). Since $\tau_{L} \gtrsim {\ell_\nu}/{V\!\left(\ell_\nu\right)}$, it follows that

\begin{equation}
\frac{\tau_{L}}{\tau_{ie}^{\rm eq}} \sim \left(\frac{m_i}{m_e}\right)^{1/2}\mathrm{Ma}_0^{3/2} \left(\frac{\lambda_i}{L_0}\right)^{1/2} \sim \epsilon_i \left(\frac{m_i}{m_e}\right)^{1/2} \, . \label{macroieequilrat}
\end{equation}

Thus, if $\epsilon_i \gg \mu_e^{1/2}$, we conclude that collisional equilibration of ion and electron temperatures might be too inefficient to regulate small-scale ion-temperature fluctuations, in which case it would follow that $L_{T_i} < L_{T}$. However, it has been previously demonstrated
via numerical solution of the Vlasov-Fokker-Planck equation that the CE expansion procedure breaks down due to
nonlocal transport effects if $\lambda_e/L$ is only moderately small [Bell et al., 1981]; thus, the only regime in which there is not ion-electron equilibration over all scales is one where
the CE expansion is not valid anyway. In short, we conclude that assuming $L_{T_i} \sim L_T$ is reasonable.

Bringing these considerations together with (<ref>), we find that

η_e^T     ∼    μ_e^1/4 Ma_0 λ_i/L_V ∼Ma_0^3/2 μ_e^1/4 (λ_i/L_0)^1/2 ∼η_e^R ∼η_e^u ∼η_i   ,

ϵ_e     ∼    μ_e^1/2 Ma_0 λ_i/L_V∼μ_e^1/2 Ma_0^3/2 (λ_i/L_0)^1/2   ,

ϵ_i     ∼    Ma_0 λ_i/L_V ∼Ma_0^3/2 (λ_i/L_0)^1/2

Thus, we conclude that the largest distortions of the ion CE distribution are due to
flow gradients, while temperature gradients cause the greatest distortions of the
electron CE distribution function.

§.§ Kinetic stability of classical, collisional plasma

§.§.§ Overview

We have seen that the CE expansion provides a procedure for the calculation of
the distribution functions arising in a classical, collisional plasma in terms of gradients of
temperature, electron-ion drifts and bulk fluid velocities; these calculations in turn allow for the
closure of the system (<ref>) of fluid equations. However, these same gradients are
sources of free energy in the plasma, so they can lead to instabilities. Some
of these instabilities will be `fluid', i.e., they are captured within the CE
description and are features of the fluid dynamics of plasmas; others
are kinetic (`microinstabilities'), and their existence implies that the CE expansion is, in fact, illegitimate.
Our primary purpose in this paper is to determine when such microinstabilities do not occur in a strongly magnetised two-species plasma.
If, however, they do occur, we wish to determine their growth rates.
We begin by making a few general qualitative comments concerning the existence and nature of these
microinstabilities, before presenting the technical details of their derivation.

§.§.§ Existence of microinstabilities in classical, collisional plasma

It might
naively be assumed that a classical, collisional plasma is kinetically stable, on two
grounds. The first of these is that the distribution function of such a plasma is `almost' Maxwellian,
and thus stable. While it is certainly the case that a plasma whose constituent particles have
Maxwellian distribution functions is kinetically stable [Bernstein, 1958, Krall & Trivelpiece, 1973], it is also known that a plasma with anisotropic particle distribution functions is typically
not [Furth, 1963, Kalman et al., 1968, Davidson, 1983, Gary, 1993].
The (small) non-Maxwellian component of the CE distribution
function is anisotropic (as, e.g., was explicitly demonstrated by the calculation of pressure anisotropy in section <ref>), and thus we cannot a priori rule out
microinstabilities associated with this anisotropy.

The second naive reason for dismissing the possibility of microinstabilities in
classical, collisional plasma is the potentially stabilising effect of collisional damping on microinstability growth
rates. If collisional processes are sufficiently dominant to
be responsible for the mediation of macroscopic momentum and heat
fluxes in the plasma, it might be naively inferred that they would also suppress microinstabilities. This is, in fact, far from guaranteed, for the following reason. The
characteristic scales of the microinstabilities are not fluid scales, but are rather
intrinsic plasma length scales related to quantities such as the Larmor radius $\rho_s$
or the inertial scale $d_s$ of species $s$, or the Debye length $\lambda_{\mathrm{D}}$
– quantities given in terms of macroscopic physical properties of plasma by

ρ_s     =     m_s v_ths c/Z_s e |B| ,

d_s     ≡    (4 Z_s^2 e^2 n_s/m_s c^2)^-1/2 = ρ_s β_s^-1/2 ,

λ_D     ≡    (∑_s 4 Z_s^2 e^2 n_s/T_s)^-1/2 = (∑_s 2 c/d_s^2 v_ths )^-1/2,

\begin{equation}
\beta_s \equiv  \frac{8 \upi n_s T_s}{B^2} \label{plasma_beta_def}
\end{equation}

is the plasma beta of species $s$. The crucial observation is then that the dynamics on characteristic
microinstability scales may
be collisionless. For a classical, collisional hydrogen plasma
(where $\lambda \equiv \lambda_e \sim \lambda_i$ for $T_e \sim T_i$), the mean free path is
much larger than the Debye length: $\lambda/\lambda_{\mathrm{D}} \sim n_e \lambda_{\mathrm{D}}^3 \gg
1$; so there exists a range of wavenumbers $k$ on which microinstabilities
are both possible ($k \lambda_{\mathrm{D}} \lesssim 1$) and collisionless
($k \lambda \gg 1$). For a strongly magnetised collisional
plasma, $\lambda_s \gg \rho_s$ for all species by definition; thus, any microinstability with a characteristic scale
comparable to the Larmor radius of any constituent particle will be
effectively collisionless. We note that such a range of collisionless wavenumbers only
exists in classical (viz., weakly coupled) plasmas; in strongly coupled plasmas, for which $\lambda \lesssim
\lambda_{\mathrm{D}}$, all hypothetically possible microinstability wavenumber scales are collisional.
Thus the phenomenon of microinstabilities in collisional plasmas is solely a concern for the classical

§.§.§ A simple example: the firehose instability in CE plasmas

Perhaps the simplest example of a microinstability that can occur in CE
plasma is the firehose instability. This example was previously discussed
by Schekochihin et al., 2005, but we nonetheless outline it here to illustrate the
central concept of our paper.

Consider bulk fluid motions of the plasma on length scales $L_V$ that are much smaller
than the mean free path $\lambda_i$, but much larger than the ion Larmor radius
$\rho_i$; the characteristic frequencies associated with these motions are assumed to be much smaller that the ion Larmor
frequency $\Omega_i$, but much larger than the inverse of the ion collision time $\tau_{i}^{-1}$.  Under these assumptions,
the following four statements can be shown to be true [Schekochihin et al., 2005]:

* The bulk velocities of the electron and ion species are approximately
equal: $\boldsymbol{V}_e \approx \boldsymbol{V}_i$.

* The electric field in a frame co-moving with the ion fluid vanishes; transforming to the stationary frame of the system,
this gives

\begin{equation}
\boldsymbol{E} = -\frac{\boldsymbol{V}_i \times \boldsymbol{B}}{c} \, .  \label{idealMHDohmslaw}
\end{equation}

* The contribution of the displacement current to the Maxwell-Ampère law
(<ref>d) is negligible, and so

\begin{equation}
e n_e \left(\boldsymbol{V}_i - \boldsymbol{V}_e \right) \approx \frac{c}{4 \upi}
\bnabla \times \boldsymbol{B} \, .
\end{equation}

* The electron and ion viscosity tensors both take the form
(<ref>), and the electron pressure anisotropy, defined by (<ref>), is small compared to the ion pressure anisotropy: $\Delta_e \ll \Delta_i$.

It then follows directly from (<ref>b), summed over both ion and electron species, that

\begin{equation}
m_i n_i {\mathrm{D} \boldsymbol{V}_i \over \mathrm{D} t} \bigg{|}_i = -\bnabla \left(\frac{B^2}{8 \upi} + p_{e\perp} + p_{i\perp}\right) -  \bnabla \bcdot \left[ \hat{\boldsymbol{z}} \hat{\boldsymbol{z}} \left(p_{i\perp}-p_{i\|}\right)
\right] + \frac{\boldsymbol{B} \bcdot \bnabla \boldsymbol{B}}{4 \upi}
\, . \label{momequation_ionscales}
\end{equation}

We remind the reader that $\hat{\boldsymbol{z}} = \boldsymbol{B}/B$, and emphasize that we have neglected
the electron inertial term on the grounds that it is small compared to the ion inertial term:

\begin{equation}
m_e n_e {\mathrm{D} \boldsymbol{V}_e \over \mathrm{D} t} \bigg{|}_e \ll m_i n_i {\mathrm{D} \boldsymbol{V}_i \over \mathrm{D} t} \bigg{|}_i \, .
\end{equation}

evolution of the magnetic field is described by the induction equation,

\begin{equation}
{\mathrm{D} \boldsymbol{B} \over \mathrm{D} t} \bigg{|}_i = \boldsymbol{B} \bcdot
\bnabla \boldsymbol{V}_{i} -  \boldsymbol{B} \bnabla \bcdot \boldsymbol{V}_{i}
\, , \label{inductioneqn}
\end{equation}

which is derived by substituting (<ref>) into Faraday's law

Now consider small-amplitude perturbations with respect to a particular macroscale state of the

δV_i     =
δV_i⊥ exp{i(k r - ωt)}  ,

δB     =
δB_⊥ exp{i(k r - ωt)}

whose characteristic frequency $\omega$ is much greater than that of the plasma's bulk fluid motions (but is still much smaller than $\Omega_i$),
whose wavevector $\boldsymbol{k} = k_{\|} \hat{\boldsymbol{z}}$ is parallel to
$\boldsymbol{B}$, and assume also that the velocity and magnetic-field perturbations are
perpendicular to $\boldsymbol{B}$. It is then easy to show that (<ref>) and (<ref>)

- i m_i n_i ωδV_i⊥      =
i (B_0^2/4  + p_i⊥ - p_i) k_  δB_⊥ /B

- i ωδB_⊥     =
i B k_ δV_i⊥   ,

where $p_{i\perp}$ and $p_{i\|}$ are the perpendicular and parallel ion pressures
associated with the macroscale state (which, on account of its comparatively slow evolution compared to the perturbation, can be
regarded a quasi-equilibrium). Physically, the macroscale flow gives rise
to different values of $p_{i\perp}$ and $p_{i\|}$, and thereby an ion pressure
anisotropy $\Delta_i$, because it changes the strength $B$ of the macroscale magnetic
field; thanks to the effective conservation of the first and second adiabatic
moments of the ions on the evolution timescale of the macroscale flow [Chew et al., 1956],
an increase (decrease) in $B$ results in an increase (decrease) in $p_{i\perp}$,
and a decrease (increase) in $p_{i\|}$.
The dispersion relation for the perturbation is then

\begin{equation}
\omega^2 = k_{\|}^2 v_{\mathrm{th}i}^2 \left(\frac{1}{\beta_i} + \frac{\Delta_i}{2}\right) \label{firehosegrowth_Sec2}
\, ,
\end{equation}

where $\beta_i $, defined by (<ref>), is the ion plasma beta. For a sufficiently negative ion pressure anisotropy, viz., $\Delta_i < -2/\beta_i$, the perturbation is unstable. This instability is known as the (parallel) firehose instability.

The underlying physics of the parallel firehose instability has been discussed extensively
elsewhere <cit.>.
Here, we simply note that the firehose instability arises in a magnetised plasma with sufficiently negative
pressure anisotropy as compared to the inverse of the ion plasma beta; because the ion CE distribution
function has a small, non-zero pressure anisotropy, this statement applies to CE plasma at
large $\beta_i$. We also observe that the product of the growth rate (<ref>) of the firehose instability with the ion-ion collision time

\begin{equation}
\omega \tau_i \sim k_{\|} \lambda_i \left|\frac{1}{\beta_i} + \Delta_i \right|^{1/2} \sim
\frac{1}{\beta_i} \frac{\lambda_i}{\rho_i}
\, ,
\end{equation}

where we have assumed that $\Delta_i \lesssim 2 \beta_i^{-1}$, and employed the (non-trivial) result that the peak growth of the parallel firehose instability occurs at wavenumbers satisfying $k_{\|} \rho_i \sim \beta_i^{-1/2}$ (see sections <ref> and
Thus, if $\beta_i \ll \lambda_i/\rho_i$
– a condition easily satisifed in weakly collisional astrophysical environments such as the ICM (see table <ref>)
– it follows that $\omega \tau_i \gg 1$, and so collisional damping is unable to inhibit the parallel
firehose in a CE plasma[In fact, the naive condition $\gamma \tau_i \lesssim 1$
is not sufficient to ensure collisional stabilisation of the firehose
instability; the true stabilisation condition is instead $k_{\|} \lambda_i \lesssim 1$
(see section <ref> for a discussion of this claim).]. This
failure is directly attributable to its characteristic wavelength being at
collisionless scales: the parallel wavenumber satisfies $k_{\|} \lambda_i \sim \beta_i^{-1/2} \lambda_i/\rho_i \gg

This simple example clearly illustrates
that microinstabilities are indeed possible in a classical, collisional plasma,
for precisely the reasons given in section <ref>.

§.§.§ Which microinstabilities are relevant

Although the naive arguments described in section <ref> do not imply kinetic
stability of CE plasma, these same arguments do lead to significant restrictions
on the type of microinstabilities that can arise. Namely, for some plasma modes,
the small anisotropy of CE distribution functions is an insufficient free-energy
source for overcoming the competing collisionless damping mechanisms that ensure stability for pure Maxwellian distribution functions – e.g., Landau damping or
cyclotron damping. For other plasma modes, the characteristic length scales
are so large that collisional damping does suppress growth. In magnetised plasmas, there also
exist cyclotron harmonic oscillations that, despite minimal damping, can only become unstable
for sufficiently large anisotropy of the particle distribution function: e.g., the electrostatic Harris instability [Harris, 1959, Hall et al., 1964]. Since the anisotropy threshold for such microinstabilities is typically $\Delta_s \gtrsim 1$ [Shima & Hall, 1965], they cannot operate in a CE plasma.

We claim that there are only two classes of microinstabilities that can be triggered in a CE plasma.
The first are quasi-cold plasma modes: these are modes whose frequency is
so large that resonant wave-particle interactions (Landau or cyclotron resonances)
only occur with electrons whose speed greatly exceeds the electron thermal speed
$v_{\mathrm{th}e}$. Collisionless damping of such modes is typically very weak,
and thus small anisotropies of particle distribution functions can be sufficient to
drive an instability. Well-known examples of a small non-Maxwellian part of
the distribution function giving rise to microinstabilities include the bump-on-tail instability associated with
a fast beam of electrons <cit.>,
or the whistler instability for small temperature anisotropies <cit.>.
The existence of such instabilities for the CE distribution can be
demonstrated explicitly: e.g., the peak growth rate of the bump-on-tail instability associated with the CE distribution function (`the CE bump-on-tail instability') is calculated
in appendix <ref>. However, the growth rates $\gamma$ of such instabilities are
exponentially small in $\lambda_e/L \ll 1$. This claim, which is explicitly
proven for the CE bump-on-tail instability in appendix
<ref>, applies to all electrostatic instabilities (see appendix <ref>), and it can be argued that it also applies to all quasi-cold plasma modes (see
appendix <ref>).
When combined with the constraint that the resonant wave-particle interactions
required for such instabilities cannot occur if $\gamma \tau_{r} \lesssim 1$, where $\tau_{r}$ is the collision time of the resonant particles,
the exponential smallness of the growth rate suggests that such microinstabilities will not be
significant provided $\lambda_e/L$ really is small. As discussed in section <ref>, plasmas in which $\lambda_e/L$ is only moderately small are not well modelled as CE plasmas anyway, and thus, for the rest of
this paper, we will not study quasi-cold-plasma-mode instabilities.

The second class of allowed microinstabilities comprises modes that are electromagnetic and low-frequency in the sense that
the complex frequency $\omega$ of the microinstability satisfies, for at least one particle species $s$,

\begin{equation}
\frac{\omega}{k v_{\mathrm{th}s}} \sim \left(\frac{\lambda_s}{L}\right)^{\iota} \label{omegascale_A}
\ll 1 ,
\end{equation}

where $\iota$ is some order-unity number. Low-frequency electromagnetic modes are
in general only subject to weak Landau and cyclotron damping (of order $\omega/k v_{\mathrm{th}s} \ll 1 $ or less), and thus can become unstable for small distribution-function anisotropies. By contrast, electromagnetic
modes satisfying $\omega \sim k v_{\mathrm{th}s}$ would typically generate strong
inductive electric fields, which would in turn be subject to significant Landau or cyclotron
damping, overwhelming any unstable tendency.
The firehose instability introduced in section <ref> is one example of this
type of microinstability: it satisfies (<ref>) with $\iota = 1/2$, provided its $\beta$-stabilisation
threshold is surpassed.

In this paper, we will focus on microinstabilities in this second class.
Whilst small compared
to the streaming rate $k v_{\mathrm{th}s}$ of species $s$, the growth
rates satisfying (<ref>) can still be significantly larger than the rate at which the plasma evolves
on macroscopic scales, and thus invalidate the CE expansion.
We do not in this paper present a rigorous proof that there are no
microinstabilities of the CE distribution function which do not fall into either of the two classes considered
above. However, there do exist more precise arguments supporting the latter claim than
those based on physical intuition just presented; these are discussed further in sections <ref> and <ref>.

The microinstabilities satisfying
(<ref>) fall into two sub-classes. The first sub-class consists of microinstabilities driven by
the CE temperature-gradient, CE electron-friction and CE electron-ion-drift terms in the CE distribution functions (<ref>); we refer to
these collectively as CE temperature-gradient-driven microinstabilities, or CET microinstabilities,
on account of the parameters $\eta_e^R$ and $\eta_e^u$ scaling with temperature gradients (see section <ref>).
The second sub-class is microinstabilities driven by the CE shear terms, or
CE shear-driven microinstabilities (CES microinstabilities). This sub-classification is
necessary for two reasons. First, the velocity-space anisotropy associated with the CE shear
terms is different from other non-Maxwellian terms, and thus different types of
microinstabilities can emerge for the two sub-classes. Secondly, as was discussed in section <ref> for the case
of CE plasma, the typical size of small parameters $\eta_e^T$, $\eta_e^R$, $\eta_e^u$ and $\eta_i$
is different from that of $\epsilon_e$ and $\epsilon_i$. In our initial
overview of our calculations (section <ref>) and in the more detailed discussion of our method (section <ref>),
we will consider all microinstabilities driven by the
non-Maxwellian terms of the CE distribution together; however, when it comes to
presenting detailed results, we will consider CET and CES microinstabilities
separately (sections <ref> and <ref>,

§.§ Linear stability calculation: overview

§.§.§ General dispersion relation

Our linear kinetic stability calculation proceeds as follows: we consider an electromagnetic perturbation with wavevector $\boldsymbol{k}$ and (complex) frequency $\omega$ of the form

δE     =
δE exp{i(k r - ωt)}  ,

δB     =
δB exp{i(k r - ωt)}

in a plasma with the equilibrium electron and ion distribution functions
given by (<ref>a) and (<ref>b), respectively. We
assume that all macroscopic parameters in the CE distribution function are
effectively constant on the time scales and length scales associated with
microinstabilities: this is equivalent to assuming that $k \lambda_e, k \lambda_i \gg 1$
(where $k \equiv |\boldsymbol{k}|$ is the wavenumber of the perturbation), and $|\omega| \tau_L \gg 1$.
To minimise confusion between quantities evolving on short, collisionless time scales, and those on long, fluid time scales,
we relabel the equilibrium number density of species $s$ as $n_{s0}$, and the macroscopic magnetic field
as $\boldsymbol{B}_0$ in subsequent calculations. For notational convenience, we define

\begin{equation}
\eta_e \equiv \eta_e^T ,
\end{equation}

\begin{equation}
A_e(\tilde{v}_e) \equiv A_e^T(\tilde{v}_e) + \frac{\eta_e^R}{\eta_e^T} A_e^R(\tilde{v}_e)
+ \frac{\eta_e^u}{\eta_e^T} A_e^u(\tilde{v}_e)
\, ,
\end{equation}

which in turn allows for the equilibrium distribution function of species $s$ to be
written as

\begin{equation}
f_{s0}(\tilde{v}_{s\|},\tilde{v}_{s\bot}) = \frac{n_{s0}}{v_{\mathrm{th}s}^3 \upi^{3/2}} \exp \left(-\tilde{v}_{s}^2\right) \left[1+\eta_s A_s(\tilde{v}_s) \tilde{v}_{s\|} + \epsilon_s C_s(\tilde{v}_s) \left(\tilde{v}_{s\|}^2- \frac{\tilde{v}_{s\bot}^2}{2} \right)\right] . \quad
\label{ChapEnskogFunc_s}
\qquad \;
\end{equation}

Finally, without loss of generality, we can set $\boldsymbol{V}_{i} = 0$
by choosing to perform the kinetic calculation in the frame of the ions; thus, $\tilde{\boldsymbol{v}}_s =
\boldsymbol{v}/v_{\mathrm{th}s}$.

It is well known [Stix, 1962, Parra, 2017] that the electric field of all linear electromagnetic perturbations in a collisionless, magnetised
plasma with equilibrium distribution function $f_{s0}$
must satisfy

\begin{equation}
\left[\frac{c^2 k^2}{\omega^2} \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right)+ \boldsymbol{\mathfrak{E}} \right]
\bcdot \widehat{\delta \boldsymbol{E}} = 0 \, , \label{singulareigenvaleqnFull}
\end{equation}

where $\hat{\boldsymbol{k}} \equiv \boldsymbol{k}/k$
is the direction of the perturbation,

\begin{equation}
\boldsymbol{\mathfrak{E}} \equiv \mathsfbi{I} + \frac{4 \upi \mathrm{i}}{\omega}
\boldsymbol{\sigma} \label{dielecttensfull}
\end{equation}

the plasma dielectric tensor, and $\boldsymbol{\sigma}$ the plasma
conductivity tensor.
The hot-plasma dispersion relation is then given by

\begin{equation}
\mbox{det}\left[\frac{c^2 k^2}{\omega^2} \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right)+ \boldsymbol{\mathfrak{E}} \right]=0
\label{hotplasmadisprel}
\end{equation}

The conductivity tensor in a hot, magnetised plasma is best displayed in an orthogonal coordinate system with
basis vectors $\left\{\hat{\boldsymbol{x}},\hat{\boldsymbol{y}},\hat{\boldsymbol{z}}\right\}$ defined in terms of $\boldsymbol{B}_0$ and

\begin{equation}
\hat{\boldsymbol{z}} \equiv \frac{\boldsymbol{B}_0}{B_0}, \quad
\hat{\boldsymbol{x}}  \equiv \frac{\boldsymbol{k}_{\perp}}{k_{\perp}} \equiv \frac{\boldsymbol{k}- k_{\|} \hat{\boldsymbol{z}}}{k_{\bot}} ,
\quad
\hat{\boldsymbol{y}}  \equiv \hat{\boldsymbol{z}} \times \hat{\boldsymbol{x}} \label{basisdef}
, \label{xyzcoordinatebasis}
\end{equation}

where $B_0 \equiv |\boldsymbol{B}_0|$, $k_{\|} \equiv \boldsymbol{k}\bcdot \hat{\boldsymbol{z}}$, and $k_{\bot} \equiv \left|\boldsymbol{k}_{\perp}\right|$.
In this notation, $\boldsymbol{k} = k_{\|} \hat{\boldsymbol{z}} + k_{\bot}
\hat{\boldsymbol{x}}$.
The conductivity tensor is then given by

\begin{eqnarray}
\boldsymbol{\sigma} = \sum_s \boldsymbol{\sigma} _s & = & - \frac{\mathrm{i}}{4 \upi \omega} \sum_s \omega_{\mathrm{p}s}^2 \bigg[ \frac{2}{\sqrt{\upi}} \frac{k_{\|}}{|k_{\|}|}  \int_{-\infty}^{\infty} \mathrm{d} \tilde{w}_{s\|} \, \tilde{w}_{s\|} \int_0^{\infty} \mathrm{d} \tilde{v}_{s\bot} \Lambda_s(\tilde{w}_{s\|},\tilde{v}_{s\bot}) \hat{\boldsymbol{z}} \hat{\boldsymbol{z}} \nonumber \\
&& \mbox{} + \tilde{\omega}_{s\|} \frac{2}{\sqrt{\upi}} \int_{C_L} \mathrm{d} \tilde{w}_{s\|} \int_0^{\infty} \mathrm{d} \tilde{v}_{s\bot} \tilde{v}_{s\bot}^2 \Xi_s(\tilde{w}_{s\|},\tilde{v}_{s\bot}) \sum_{n = -\infty}^{\infty} \frac{\mathsfbi{R}_{sn}}{\zeta_{sn} -\tilde{w}_{s\|}}  \bigg] \, , \qquad \label{conductivity}
\end{eqnarray}

\begin{equation}
\omega_{\mathrm{p}s} \equiv \sqrt{\frac{4 \upi Z_s^2 e^2 n_{s0}}{m_s}} , \label{plasmafrequency_def}
\end{equation}

\begin{equation}
\tilde{w}_{s\|} \equiv \frac{k_{\|} \tilde{v}_{s\|}}{|k_{\|}|},   \label{w_pl_def}
\end{equation}

\begin{equation}
\tilde{\rho}_s \equiv \frac{m_s c v_{\mathrm{th}s}}{Z_{s} e B_0} =  \frac{|Z_{s}|}{Z_{s}} \rho_s,
\end{equation}

\begin{equation}
\tilde{\omega}_{s\|} \equiv \frac{\omega}{|k_{\|}| v_{\mathrm{th}s}} , \label{normfrequency}
\end{equation}

\begin{equation}
\zeta_{sn} \equiv \tilde{\omega}_{s\|} - \frac{n}{|k_{\|}| \tilde{\rho}_s} , \label{zeta_def}
\end{equation}

\begin{equation}
\tilde{f}_{s0}(\tilde{v}_{s\|},\tilde{v}_{s\bot}) \equiv \frac{\upi^{3/2} v_{\mathrm{th}s}^3}{n_{s0}} f_{s0}\left(\frac{k_{\|}}{|k_{\|}|} v_{\mathrm{th}s} \tilde{w}_{s\|},v_{\mathrm{th}s} \tilde{v}_{s\bot}\right)
, \label{renorm_distfunc_def}
\end{equation}

\begin{equation}
\Lambda_s(\tilde{w}_{s\|},\tilde{v}_{s\bot}) \equiv \tilde{v}_{s\bot} \frac{\p \tilde{f}_{s0}}{\p \tilde{w}_{s\|}}-\tilde{w}_{s\|} \frac{\p \tilde{f}_{s0}}{\p \tilde{v}_{s\bot}}
, \label{anisotropyfunc}
\end{equation}

\begin{equation}
\Xi_s(\tilde{w}_{s\|},\tilde{v}_{s\bot}) \equiv \frac{\p \tilde{f}_{s0}}{\p \tilde{v}_{s\bot}}
+ \frac{\Lambda_s(\tilde{w}_{s\|},\tilde{v}_{s\bot})}{\tilde{\omega}_{s\|}} , \label{IntgradCond}
\end{equation}

(R_sn )_xx     ≡    n^2 J_n(k_ ρ̃_s ṽ_s)^2/k_^2 ρ̃_s^2 ṽ_s^2 ,

(R_sn )_xy     ≡    i n J_n(k_ ρ̃_s ṽ_s) J_n'(k_ ρ̃_s ṽ_s)/k_ ρ̃_s ṽ_s ,

(R_sn )_xz     ≡    n J_n(k_ ρ̃_s ṽ_s)^2/k_ ρ̃_s ṽ_s k_ w̃_s/|k_| ṽ_s  ,

(R_sn )_yx     ≡    - (R_sn )_xy ,

(R_sn )_yy     ≡    J_n'(k_ ρ̃_s ṽ_s)^2 ,

(R_sn )_yz     ≡    -i n J_n(k_ ρ̃_s ṽ_s) J_n'(k_ ρ̃_s ṽ_s) k_ w̃_s/|k_| ṽ_s ,

(R_sn )_zx     ≡    (R_sn )_xz ,

(R_sn )_zy     ≡    -(R_sn )_yz ,

(R_sn )_zz     ≡    w̃_s^2/ṽ_s^2 J_n(k_ ρ̃_s ṽ_s)^2 .

Here $(\mathsfbi{R}_{sn})_{xy} = \hat{\boldsymbol{x}} \bcdot \mathsfbi{R}_{sn} \bcdot \hat{\boldsymbol{y}}$, and similarly
for other components of $\mathsfbi{R}_{sn}$. For the reader's convenience, a summary of the derivation of the hot-plasma
dispersion relation is given in appendix <ref>.

We note that the dielectric and conductivity tensors have the following symmetries:

\begin{equation}
\mathfrak{E}_{yx} = - \mathfrak{E}_{xy} \, , \quad  \mathfrak{E}_{zx} = \mathfrak{E}_{xz} \, ,
\quad
\mathfrak{E}_{zy} = - \mathfrak{E}_{yz} \, , \label{dielectricsimsgen}
\end{equation}

\begin{equation}
\sigma_{yx} = - \sigma_{xy} \, , \quad \sigma_{zx} = \sigma_{xz} \, , \quad \sigma_{zy} = - \sigma_{yz} \, , \label{conductsimsgen}
\end{equation}

where, for tensors with no species subscript, we use the notation $ \mathfrak{E}_{xy} \equiv \hat{\boldsymbol{x}} \bcdot  \boldsymbol{\mathfrak{E}} \bcdot
\hat{\boldsymbol{y}}$. We also observe that if $f_{s0}(v_{\|},v_{\bot})$
is an even function with respect to $v_{\|}$, then, for $k_{\|} > 0$,

σ_xx(-k_)     =     σ_xx(k_)   ,

σ_xy(-k_)     =     σ_xy(k_)   ,

σ_xz(-k_)     =     -σ_xz(k_)   ,

σ_yy(-k_)     =     σ_yy(k_)   ,

σ_yz(-k_)     =     -σ_yz(k_)   ,

σ_zz(-k_)     =     σ_zz(k_)   ,

with the remaining components of the conductivity tensor given by equations
(<ref>). If $f_{s0}(v_{\|},v_{\bot})$
is an odd function with respect to $v_{\|}$, then

σ_xx(-k_)     =     -σ_xx(k_)   ,

σ_xy(-k_)     =     -σ_xy(k_)   ,

σ_xz(-k_)     =    σ_xz(k_)   ,

σ_yy(-k_)     =    -σ_yy(k_)   ,

σ_yz(-k_)     =     σ_yz(k_)   ,

σ_zz(-k_)     =    -σ_zz(k_)   .

These symmetries can be used to determine completely the behaviour of perturbations with $k_{\|} < 0$
directly from perturbations with  $k_{\|} > 0$, without any additional
calculations. Thus, unless stated otherwise, from this point on, we assume $k_{\|} >
0$, and thus $\tilde{w}_{s\|} = \tilde{v}_{s\|}$ [see (<ref>)].

§.§.§ Simplifications of dispersion relation: overview of our approach

The full hot-plasma dispersion relation (<ref>) is a transcendental equation, and
thus, for general distribution functions, the growth rates of perturbations can only be determined numerically;
this hinders the systematic investigation of stability over wide-ranging parameter regimes.
However, adopting a few simplifications both to the form of the CE distribution
functions (<ref>) and to the type of microinstabilities being
considered (see section <ref>) turns out to be advantageous when attempting a systematic study. It enables us to obtain simple analytical results
for microinstability growth rates and characteristic wavenumbers,
as well as greatly reducing the numerical cost of evaluating these quantities.
The former allows us to make straightforward comparisons between
microinstabilities, while the latter facilitates the calculation of stability plots over a wide range of parameters
without requiring intensive computational resources.

First, we choose a Krook collision operator, with constant collision time $\tau_s$ for each species $s$ [Bhatnagar et al., 1954], when evaluating the isotropic
functions $A_e^T(\tilde{v}_e)$, $A_e^R(\tilde{v}_e)$, $A_e^u(\tilde{v}_e)$,
$A_i(\tilde{v}_i)$, $C_e(\tilde{v}_e)$, and $C_i(\tilde{v}_i)$ in (<ref>).
As was explained in section <ref>, these functions are determined
by the collision operator.
While the full Landau collision operator might seem to be the most appropriate choice, the conductivity tensor $\boldsymbol{\sigma}$ defined
by (<ref>) cannot be written
in terms of standard mathematical functions if this choice is made. Instead, the relevant integrals must be done numerically.
If a simplified collision operator is assumed, $\boldsymbol{\sigma}$
can be evaluated analytically with only a moderate amount of algebra. In
appendix <ref>, we show that for the Krook collision operator,

A_e^T(ṽ_e)     =     -(ṽ_e^2-5/2)   ,

A_e^R(ṽ_e)     =     -1   ,

A_e^u(ṽ_e)     =     0   ,

A_i(ṽ_e)     =     -(ṽ_i^2-5/2)   ,

C_e(ṽ_e)     =     -1   ,

C_i(ṽ_i)     =     -1   ,

where it is assumed that $\tilde{v}_e, \tilde{v}_i \ll \eta_e^{-1/3}, \epsilon_i^{-1/2}$ in order that the
CE distribution functions retain positive signs (the vanishing of the CE electron-ion-drift term is discussed in appendix <ref>).
Adopting the Krook collision operator has the additional advantage of allowing
a simple prescription for collisional damping of microinstabilities to be introduced self-consistently
into our stability calculation (see section <ref> for further discussion of this).

Secondly, as discussed in section <ref>, the most important microinstabilities associated with the CE
distribution function are low-frequency, i.e., they
satisfy (<ref>). Therefore, instead of solving the full hot-plasma dispersion
relation, we can obtain
a less complicated algebraic dispersion relation. We also always consider
electromagnetic rather than electrostatic perturbations. This is because it can
be shown for a CE plasma that purely electrostatic microinstabilities are limited
to the quasi-cold plasma modes (see appendix
<ref>). Describing how the simplified dispersion relation for
low-frequency, electromagnetic perturbations is obtained from the full hot-plasma
dispersion relation requires a rather lengthy exposition, and necessitates the introduction of a substantial amount of additional mathematical
notation. In addition to this, certain shortcomings of this approach warrant an extended discussion.
Readers who are interested these details will find them in the next section (section <ref>).
Readers who are instead keen to see the results of the stability calculations as soon as possible are encouraged to jump to
sections <ref> and <ref>.

§.§ Linear stability calculation: detailed methodology

§.§.§ Low-frequency condition in a magnetised plasma

Before applying to the hot-plasma dispersion relation
(<ref>) the simplifications discussed in section <ref>, we refine the low-frequency condition (<ref>)
based on the specific form (<ref>) of the conductivity tensor for a magnetised plasma.
It is clear that the equilibrium distribution function only affects the conductivity tensor via the functions $\Lambda_s(\tilde{v}_{s\|},\tilde{v}_{s\bot})$
and $\Xi_s(\tilde{v}_{s\|},\tilde{v}_{s\bot})$ [see (<ref>) and (<ref>)]. For a distribution function of the form (<ref>),
it can be shown that

\begin{equation}
\Lambda_s(\tilde{v}_{s\|},\tilde{v}_{s\bot})= - \tilde{v}_{s\bot} \exp \left(-\tilde{v}_{s}^2\right) \left[\eta_s A_s(\tilde{v}_s) - 3 \epsilon_s C_s(\tilde{v}_s)  \tilde{v}_{s\|} \right]
, \label{anisopparameter}
\end{equation}

\begin{eqnarray}
\Xi_s(\tilde{v}_{s\|},\tilde{v}_{s\bot}) & = &  - \tilde{v}_{s\bot} \exp \left(-\tilde{v}_{s}^2\right) \Bigg[ 2 + 2 \tilde{v}_{s\|} \eta_s A_s(\tilde{v}_s) - \frac{\tilde{v}_{s\|}}{\tilde{v}_s} \eta_s A_s'(\tilde{v}_s)
\nonumber \\
&& \quad + 2 \epsilon_s C_s(\tilde{v}_s) \left(\tilde{v}_{s\|}^2-\frac{\tilde{v}_{s\bot}^2}{2}+\frac{1}{2}\right)
- \frac{1}{\tilde{v}_s} \left(\tilde{v}_{s\|}^2-\frac{\tilde{v}_{s\bot}^2}{2}\right)
\epsilon_s C_s'(\tilde{v}_s) \nonumber \\
&& \quad \quad + \frac{\eta_s}{ \tilde{\omega}_{s\|}} A_s(\tilde{v}_s) - 3 \frac{\epsilon_s}{ \tilde{\omega}_{s\|}} C_s(\tilde{v}_s)  \tilde{v}_{s\|} \Bigg]
\, ,\label{Xi_CE}
\end{eqnarray}

where the first term in the square brackets in (<ref>) originates from the Maxwellian part of the distribution function.
A comparison of the size of the second, third, fourth, and fifth terms with the first
indicates that for $\tilde{v}_s \sim 1$ – for which $\Xi_s$ attains its largest characteristic values – the non-Maxwellian terms of the CE distribution function only provide a small, $\textit{O}(\eta_e, \epsilon_e)$
contribution, and thus the conductivity is only altered slightly.
However, considering the sixth and seventh terms in the square brackets in (<ref>)
(which are only present thanks to the anisotropy of the CE distribution function), it
is clear that the non-Maxwellian contribution to the conductivity tensor can be
significant for $\tilde{v}_s \sim 1$ provided the frequency (<ref>)
satisfies one of

\begin{equation}
\tilde{\omega}_{s\|} \sim \eta_s \ll 1  \quad \mathrm{or} \quad \tilde{\omega}_{s\|} \sim \epsilon_s \ll 1 \, . \label{omegascale}
\end{equation}

Thus, the relevant low-frequency condition in a magnetised plasma involves the
parallel particle streaming rate $k_{\|} v_{\mathrm{th}s}$.

There do exist certain
caveats to the claim that it is necessary for microinstabilities of CE
plasma to satisfy (<ref>); we defer detailed statement and discussion of these caveats – as well as of other potential shortcomings of our approach – to sections
<ref>, <ref> and

§.§.§ Simplification I: non-relativistic electromagnetic fluctuations

The requirement that the mode be electromagnetic, combined with the fact we are
interested in non-relativistic fluctuations ($\omega \ll k c$) enables our first simplification. We see from (<ref>) that for any
perturbation of interest, the dielectric tensor must satisfy $\|\boldsymbol{\mathfrak{E}}\| \gtrsim k^2
c^2/\omega^2 \gg 1$ (where $\| \cdot \|$ is the Euclidean tensor norm); therefore, it simplifies to

\begin{equation}
\boldsymbol{\mathfrak{E}} \approx \frac{4 \upi \mathrm{i}}{\omega}
\boldsymbol{\sigma} \, . \label{dielectric_nonrelsimp}
\end{equation}

This amounts to ignoring the displacement current in the
Ampère-Maxwell law, leaving Ampère's original equation. For convenience of exposition, we
denote the contribution of each species $s$ to (<ref>) by

\begin{equation}
\boldsymbol{\mathfrak{E}}_s \equiv \frac{4 \upi \mathrm{i}}{\omega}
\boldsymbol{\sigma}_s \, . \label{dielectric_species_s}
\end{equation}

§.§.§ Simplification II: expansion of dielectric tensor in $\omega \ll k_{\|} v_{\mathrm{th}s}$

The next simplification involves an expansion of the matrices $\boldsymbol{\mathfrak{E}}_s$ in
the small parameters $\tilde{\omega}_{s\|} \sim \eta_s \sim \epsilon_s \ll 1$. The general principle of the expansion is as follows. We first divide the matrix $\boldsymbol{\mathfrak{E}}_s$ [see (<ref>), (<ref>), and (<ref>)] into the Maxwellian contribution $\mathsfbi{M}_s$ and the non-Maxwellian one $\mathsfbi{P}_s$:

\begin{equation}
\boldsymbol{\mathfrak{E}}_s = \frac{\omega_{\mathrm{p}s}^2}{\omega^2} \left(\mathsfbi{M}_s + \mathsfbi{P}_s \right) , \label{Maxnonmaxsep_s}
\end{equation}

where the $\omega_{\mathrm{p}s}^2/\omega^2$ factor is introduced for later convenience. Next, we note
that for a Maxwellian distribution, $\Lambda_s(\tilde{v}_{s\|},\tilde{v}_{s\bot}) =
0$ [see (<ref>)], whereas $\Lambda_s \sim \epsilon_{s}, \eta_s$ for the non-Maxwellian component of the CE distribution function. Thus, from (<ref>)
considered under the ordering $k \rho_s \sim 1$, $\mathsfbi{M}_s = \textit{O}(\tilde{\omega}_{s\|})$
as $\tilde{\omega}_{s\|} \rightarrow 0$, while $\mathsfbi{P}_s = \textit{O}(\eta_s,
\epsilon_s)$. The expansion of $\mathsfbi{M}_s$ and $\mathsfbi{P}_s$ in $\tilde{\omega}_{s\|}$
is, therefore,

M_s(ω̃_s,k)     ≡    ω̃_s M_s^(0)(k  )  + ω̃_s^2 M_s^(1)(k ) + ...

P_s(ω̃_s,k)      ≡    P_s^(0)(k) + ω̃_s P_s^(1)(k) + ...

where the matrices $\mathsfbi{M}_s^{(0)}$ and
$\mathsfbi{M}_s^{(1)}$ are $\textit{O}(1)$ functions of $\boldsymbol{k}$ only, and $\mathsfbi{P}_s^{(0)}$
and $\mathsfbi{P}_s^{(1)}$ are $\textit{O}(\eta_s,
\epsilon_s)$. We then expand $\boldsymbol{\mathfrak{E}}_s$ as follows:

\begin{equation}
\boldsymbol{\mathfrak{E}}_s = \tilde{\omega}_{s\|} \boldsymbol{\mathfrak{E}}_s^{(0)} + \tilde{\omega}_{s\|}^2 \boldsymbol{\mathfrak{E}}_s^{(1)}
+ ... \, , \label{dielectric_expand}
\end{equation}

𝔈_s^(0)     ≡    ω_ps^2/ω^2 [ M_s^(0)(k )  + 1/ω̃_s P_s^(0)(k) ]

𝔈_s^(1)      ≡    ω_ps^2/ω^2 [ M_s^(1)(k ) + 1/ω̃_s  P_s^(1)(k) ]

§.§.§ Additional symmetries of low-frequency dielectric tensor $ \boldsymbol{\mathfrak{E}}_s^{(0)} $

The tensor $\boldsymbol{\mathfrak{E}}_s^{(0)}$ defined by (<ref>a) has some rather convenient additional symmetries, which lead to significant
simplification of the dispersion relation.
In appendix <ref> we show that in combination with the general symmetries
(<ref>), which apply to $\boldsymbol{\mathfrak{E}}_s^{(0)}$ in addition to $\boldsymbol{\mathfrak{E}}$,
for any distribution function of particle species $s$ with a
small anisotropy,

(𝔈_s^(0))_xz     =     - k_/k_ (𝔈_s^(0))_xx   ,

(𝔈_s^(0))_yz     =     k_/k_ (𝔈_s^(0))_xy   ,

(𝔈_s^(0))_zz     =     k_^2/k_^2 (𝔈_s^(0))_xx   .

These symmetries have the consequence that

\begin{equation}
\hat{\boldsymbol{k}} \bcdot \boldsymbol{\mathfrak{E}}_s^{(0)} =
\boldsymbol{\mathfrak{E}}_s^{(0)}
\bcdot \hat{\boldsymbol{k}}
= 0 \, . \label{dielectric_zeroth_order}
\end{equation}

As a result of this identity, it is convenient to calculate the components of $\boldsymbol{\mathfrak{E}}_s^{(0)}$ (and $\boldsymbol{\mathfrak{E}}_s$)
in the coordinate basis $\{\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3\}$
defined by

\begin{equation}
\boldsymbol{e}_1 \equiv \hat{\boldsymbol{y}} \times \hat{\boldsymbol{k}}
\, , \quad
\boldsymbol{e}_2 \equiv \hat{\boldsymbol{y}} \, , \quad
\boldsymbol{e}_3 \equiv \hat{\boldsymbol{k}} \, . \label{e_coordinatebasis_def}
\end{equation}

Carrying out this calculation (see appendix <ref>), we find

(𝔈_s^(0))_11     =      k^2/k_^2 (𝔈_s^(0))_xx    ,

(𝔈_s^(0))_12     =     -  (𝔈_s^(0))_21 = k/k_  (𝔈_s^(0))_xy    ,

(𝔈_s^(0))_22     =      (𝔈_s^(0))_yy   ,

(𝔈_s^(0))_13     =      (𝔈_s^(0))_31  =  (𝔈_s^(0))_23 =  (𝔈_s^(0))_32 =  (𝔈_s^(0))_33 = 0   ,

where $(\boldsymbol{\mathfrak{E}}_s^{(0)})_{ij}$ is the $(i,j)$-th component of $\boldsymbol{\mathfrak{E}}_s^{(0)}$ in the basis $\{\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3
\}$. We conclude that, if $k \rho_s \sim 1$ and $\tilde{\omega}_{s\|}  \ll 1$, the components of $\boldsymbol{\mathfrak{E}}_s$

\begin{equation}
(\boldsymbol{\mathfrak{E}}_s)_{13} \sim (\boldsymbol{\mathfrak{E}}_s)_{23} \sim (\boldsymbol{\mathfrak{E}}_s)_{33} \sim \tilde{\omega}_{s\|} (\boldsymbol{\mathfrak{E}}_s)_{11} \sim \tilde{\omega}_{s\|}  (\boldsymbol{\mathfrak{E}}_s)_{12} \sim \tilde{\omega}_{s\|}  (\boldsymbol{\mathfrak{E}}_s)_{22} \, .\label{dielectric3_smallterm}
\end{equation}

These components can be written in terms of the components of
$\boldsymbol{\mathfrak{E}}_s$ in the
$\{\hat{\boldsymbol{x}},\hat{\boldsymbol{y}},\hat{\boldsymbol{z}}\}$ coordinate
frame [see (<ref>)] via a coordinate transformation;
the resulting expressions are rather bulky, so we do not reproduce them here – they are detailed in appendix

§.§.§ Consequences for dispersion relation

On account of the additional symmetries described in the previous section, a
simplified dispersion relation for low-frequency modes can be derived in place of the full
hot-plasma dispersion relation (<ref>). However, depending on the frequency and characteristic wavelengths of modes, this derivation
has a subtlety because of the large discrepancy between ion and electron masses.
In, e.g., a two-species plasma with $\mu_e = m_e/m_i \ll 1$ (and ion charge
$Z$), we have

\begin{equation}
\frac{\tilde{\omega}_{e\|}}{\tilde{\omega}_{i\|}} = \sqrt{\mu_e \tau} \,
, \label{ratio_tildeomegas}
\end{equation}

where $\tau = T_i/T_e$. If $\tau \sim 1$ [as would be expected in a collisional plasma on macroscopic evolution time scales $\tau_L$ greater than the ion-electron temperature equilibration time $\tau_{ie}^{\rm eq}$ – cf. (<ref>)], then $\tilde{\omega}_{i\|} \sim {\mu_e}^{-1/2} \tilde{\omega}_{e\|} \gg
\tilde{\omega}_{e\|}$. Thus, in general, $\tilde{\omega}_{i\|} \not \sim \tilde{\omega}_{e\|}$, and any dispersion relation will in principle depend on an additional (small) dimensionless parameter $\mu_e$. This introduces various complications to the simplified dispersion relation's derivation, most significant of which being that, since $\rho_e = Z \mu_e^{1/2} \tau^{-1/2} \rho_i \ll \rho_i$ (for $Z \gtrsim 1$), to assume the ordering $k \rho_s \sim 1$
for both ions and electrons is inconsistent (see section <ref>).

To avoid the description of our approach being obscured by these complications, we consider a special case at first: we adopt the ordering $k \rho_e \sim 1$ in a two-species plasma and assume that $\tilde{\omega}_{i\|} \sim \mu_e^{-1/2} \tilde{\omega}_{e\|} \ll 1$. In this case, $\tilde{\omega}_{i\|} \|\boldsymbol{\mathfrak{E}}_i^{(0)}\| \sim \mu_e^{1/2} Z \tau^{-1/2} \tilde{\omega}_{e\|} \|\boldsymbol{\mathfrak{E}}_e^{(0)}\| \ll \tilde{\omega}_{e\|} \|\boldsymbol{\mathfrak{E}}_e^{(0)}\|$, and so the dielectric tensor $\boldsymbol{\mathfrak{E}}$ is given by

\begin{equation}
\boldsymbol{\mathfrak{E}} = \tilde{\omega}_{e\|} \boldsymbol{\mathfrak{E}}^{(0)} + \tilde{\omega}_{e\|}^2 \boldsymbol{\mathfrak{E}}^{(1)}
+ ... \, , \label{dielectric_expand_multispecies}
\end{equation}

𝔈^(0)     ≡    𝔈_e^(0) + ω̃_i/ω̃_e 𝔈_i^(0) ≈𝔈_e^(0)

𝔈_s^(1)      ≡    𝔈_e^(1) + ω̃_i^2/ω̃_e^2 𝔈_i^(1)

Thus, to leading order in the $\tilde{\omega}_{e\|} \ll 1$ expansion, only the electron species contributes to the dielectric tensor for electron-Larmor-scale modes. We revisit the derivation of simplified dispersion relations for CE microinstabilities more generally in section <ref>.

To derive the simplified dispersion relation for electron-Larmor-scale modes, we start by considering the
component of (<ref>) for the electric field that is parallel to the wavevector $\hat{\boldsymbol{k}}$,

\begin{equation}
\hat{\boldsymbol{k}} \bcdot \boldsymbol{\mathfrak{E}} \bcdot \widehat{\delta \boldsymbol{E}}
= 0 \, , \label{singulareigenvaleqnFull_33}
\end{equation}

and then substitute the expanded form (<ref>) of the
dielectric tensor (with $s = e$). The orthogonality of $\boldsymbol{\mathfrak{E}}_e^{(0)}$ to $\hat{\boldsymbol{k}}$ – viz., (<ref>)
– implies that (<ref>) becomes

\begin{equation}
\hat{\boldsymbol{k}} \bcdot \boldsymbol{\mathfrak{E}}^{(1)} \bcdot \widehat{\delta \boldsymbol{E}}
=\mathfrak{E}_{33}^{(1)} \hat{\boldsymbol{k}} \bcdot \widehat{\delta \boldsymbol{E}} + \hat{\boldsymbol{k}} \bcdot \boldsymbol{\mathfrak{E}}^{(1)} \bcdot \widehat{\delta
\boldsymbol{E}}_T =\textit{O}(\tilde{\omega}_{e\|} |\widehat{\delta
\boldsymbol{E}}|) \, , \label{consistencycond}
\end{equation}

where the transverse electric field is
defined by $\widehat{\delta \boldsymbol{E}}_{T} \equiv \widehat{\delta \boldsymbol{E}} \bcdot \left(\mathsfbi{I} -  \hat{\boldsymbol{k}} \hat{\boldsymbol{k}}
\right)$.
In appendix <ref>, we show that for $\tilde{\omega}_{e\|}, \tilde{\omega}_{i\|} \ll 1$,

\begin{equation}
\mathfrak{E}_{33}^{(1)} \approx \frac{\omega_{\mathrm{p}e}^2}{\omega^2} \frac{2 k_{\|}^2}{k^2} (1+Z\tau^{-1}) \left[1+ \textit{O}(\eta_e, \epsilon_e)\right]
\, .
\end{equation}

Since this is strictly positive, we can rewrite (<ref>) to
give the electrostatic field in terms of the transverse electric field:

\begin{equation}
\hat{\boldsymbol{k}} \bcdot \widehat{\delta \boldsymbol{E}} = -\left(\mathfrak{E}_{33}^{(1)} \right)^{-1} \left(\hat{\boldsymbol{k}} \bcdot \boldsymbol{\mathfrak{E}}^{(1)} \bcdot \widehat{\delta
\boldsymbol{E}}_T \right)
\, . \label{electrostaticcomp}
\end{equation}

We conclude that $|\hat{\boldsymbol{k}} \bcdot \widehat{\delta \boldsymbol{E}}| \sim |\widehat{\delta
\boldsymbol{E}}_T|$ for all low-frequency perturbations with $k_{\|} \sim k$; a corollary of this result is that there can be no low-frequency
purely electrostatic perturbations (see appendix
<ref> for an alternative demonstration of this).

We can now derive the dispersion relation from the other two components of (<ref>),

\begin{equation}
\left[\frac{c^2 k^2}{\omega^2} \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right)+  \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right) \bcdot \boldsymbol{\mathfrak{E}} \right]
\bcdot \widehat{\delta \boldsymbol{E}} = 0 \, , \label{singulareigenvaleqnFull_transverse}
\end{equation}

by (again) substituting the expanded dielectric tensor (<ref>) into (<ref>):

\begin{equation}
\left[\tilde{\omega}_{e\|} \boldsymbol{\mathfrak{E}}^{(0)} + \frac{c^2 k^2}{\omega^2} \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right) \right]
\bcdot \widehat{\delta \boldsymbol{E}}_T = - \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right) \bcdot \left( \boldsymbol{\mathfrak{E}} - \tilde{\omega}_{e\|} \boldsymbol{\mathfrak{E}}^{(0)} \right) \bcdot \widehat{\delta \boldsymbol{E}} \, , \label{perpdieletricA}
\end{equation}

where we have used the identity

\begin{equation}
\boldsymbol{\mathfrak{E}}^{(0)} =  \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right) \bcdot  \boldsymbol{\mathfrak{E}}^{(0)}
\bcdot  \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right) \, ,
\end{equation}

and ordered $k^2 c^2/\omega^2 \sim \tilde{\omega}_{e\|}
\|\boldsymbol{\mathfrak{E}}^{(0)}\|$. The ratio of the right-hand side of (<ref>)
to the left-hand side is $\textit{O}(\tilde{\omega}_{e\|})$; we thus conclude
that, to leading order in the $\tilde{\omega}_{e\|} \ll 1$ expansion,

\begin{equation}
\left[\tilde{\omega}_{e\|} \boldsymbol{\mathfrak{E}}_e^{(0)} + \frac{c^2 k^2}{\omega^2} \left(\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}-\mathsfbi{I}\right) \right]
\bcdot \widehat{\delta \boldsymbol{E}}_T = 0 \, , \label{perpdieletric}
\end{equation}

and the dispersion relation is approximately

\begin{equation}
\left[\tilde{\omega}_{e\|} (\boldsymbol{\mathfrak{E}}_e^{(0)})_{11}-\frac{k^2 c^2}{\omega^2}\right]\left[\tilde{\omega}_{e\|} (\boldsymbol{\mathfrak{E}}_e^{(0)})_{22}-\frac{k^2 c^2}{\omega^2}\right]+\left[\tilde{\omega}_{e\|} (\boldsymbol{\mathfrak{E}}_e^{(0)})_{12}\right]^2
= 0 \, . \label{disprel_simp_1}
\end{equation}

Finally, writing the dielectric tensor in terms of $\mathsfbi{M}_{e}$ and
$\mathsfbi{P}_{e}$ as defined by (<ref>a), we find

\begin{eqnarray}
\left[\tilde{\omega}_{e\|} (\mathsfbi{M}_{e}^{(0)})_{11} + (\mathsfbi{P}_{e}^{(0)})_{11}- k^2 d_e^2\right]\left[\tilde{\omega}_{e\|} (\mathsfbi{M}_{e}^{(0)})_{22} + (\mathsfbi{P}_{e}^{(0)})_{22} - k^2 d_e^2\right]
\nonumber \\
+ \left[\tilde{\omega}_{e\|} (\mathsfbi{M}_{e}^{(0)})_{12} + (\mathsfbi{P}_{e}^{(0)})_{12}\right]^2 = 0 \,
\qquad \label{disprel_simp_2}
\end{eqnarray}

where $d_e = c/\omega_{\mathrm{p}e}$ is the electron inertial scale [see (<ref>b)]. This can be re-written as a quadratic equation in $\omega$ – and thus, expressions for the complex
frequency of any low-frequency perturbation can be found for any given positive wavenumber.
We note that the electron inertial scale is related to the electron Larmor
radius by $d_e = \rho_e \beta_e^{-1/2}$; therefore, our expansion scheme is only
consistent with the low-frequency assumption (<ref>) under our assumed ordering,
$\tilde{\omega}_{e\|} \sim \beta_e^{-1}$, when $\beta_e \gg 1$.

We note that one only needs to know $\boldsymbol{\mathfrak{E}}_e^{(0)}$ in order
to obtain the dispersion relation of low-frequency perturbations and the transverse component of the
electric field, whereas to determine the electrostatic component of the electric
field (and other quantities, such as the density perturbation – see appendix <ref>), one must go to higher order in the $\tilde{\omega}_{e\|} \ll 1$
expansion. Since we are primarily interested in microinstability growth rates and
wavenumber scales, we will not explicitly calculate the
electrostatic fields associated with perturbations using
(<ref>), and thus can avoid the rather laborious calculation of
$\boldsymbol{\mathfrak{E}}^{(1)}$ for CE distribution functions.
We do, however, in appendix <ref> derive an explicit expression for
$\boldsymbol{\mathfrak{E}}^{(1)}$ for a plasma with Maxwellian distribution
functions for all particle species;
this in turn allows us to relate the electrostatic
electric field to the transverse field
for such a plasma (see appendix <ref>).

For the sake of completeness, we also observe that if the non-Maxwellian part
of the CE distribution function is even with respect to
$v_{\|}$, the transformation rules (<ref>) combined with (<ref>) imply
that a perturbation with a negative parallel wavenumber $k_{\|}$
will obey exactly the same dispersion relation as a perturbation for a
positive parallel wavenumber, viz., for $k_{\|} > 0$,

\begin{equation}
\mathsfbi{P}_e^{(0)}\!\left(-k_{\|}, k_{\bot}\right) = \mathsfbi{P}_e^{(0)}\!\left(k_{\|}, k_{\bot} \right) \, .
\end{equation}

If instead the non-Maxwellian part
is odd, then, for $k_{\|} > 0$,

\begin{equation}
\mathsfbi{P}_e^{(0)}\!\left(-k_{\|}, k_{\bot}\right) = -\mathsfbi{P}_e^{(0)}\!\left(k_{\|}, k_{\bot}
\right) \, .
\end{equation}

The dispersion relation for perturbations with $k_{\|} < 0$ can, therefore, be recovered by considering
perturbations with $k_{\|} > 0$, but under the substitution $\mathsfbi{P}_{e}^{(0)}
\rightarrow -\mathsfbi{P}_{e}^{(0)}$. Thus, we can
characterise all unstable perturbations under the assumption that $k_{\|} > 0$.

In all subsequent calculations, we require the Maxwellian part
$\mathsfbi{M}_e^{(0)}$ of the dielectric tensor. The elements of the matrix $\mathsfbi{M}_s^{(0)}$ of species $s$ are as follows:

(M_s^(0))_11     =     i k^2/k_^2  F(k_ ρ̃_s,k_ ρ̃_s)   ,

(M_s^(0))_12     =     -i k/k_ G(k_ ρ̃_s,k_ ρ̃_s)   ,

(M_s^(0))_21     =     i k/k_ G(k_ ρ̃_s,k_ ρ̃_s)   ,

(M_s^(0))_22     =     i H(k_ ρ̃_s,k_ ρ̃_s)   ,

where the functions $F\!\left(x,y\right)$, $G\!\left(x,y\right)$ and $H\!\left(x,y\right)$ are

F(x,y)     ≡    4 √()/y^2 exp(-y^2/2) ∑_m=1^∞ m^2 I_m(y^2/2)exp(-m^2/x^2)   ,

G(x,y)     ≡    exp(-y^2/2) ∑_m = -∞^∞ m     Z(m/x) [ I_m'(y^2/2)- I_m(y^2/2)]   ,

H(x,y)     ≡     F(x,y) + √() y^2 exp(-y^2/2) ∑_m = -∞^∞ [ I_m(y^2/2)- I_m'(y^2/2)] exp(-m^2/x^2)   ,

$I_m\!\left(\alpha\right)$ is the $m$-th modified Bessel function, and

\begin{equation}
Z\!\left(z\right) = \frac{1}{\sqrt{\upi}} \int_{C_L} \frac{\mathrm{d}u \exp{\left(-u^2\right)}}{u-z} \label{Plasma_disp_func_def}
\end{equation}

is the plasma dispersion function ($C_L$ is the Landau contour) [Fried & Conte, 1961].  The derivation of
these results from the full dielectric tensor (which is calculated in appendix <ref>) for a plasma whose constituent particles all have Maxwellian distributions is presented in
Appendices <ref>  (expansion in the $\{\hat{\boldsymbol{x}},\hat{\boldsymbol{y}},\hat{\boldsymbol{z}}\}$ basis) and <ref> (expansion in the $\{\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3
\}$ basis).

§.§.§ Effect of multiple species on dispersion-relation derivations

We now relax the assumptions adopted in section <ref> that the low-frequency modes of interest are on electron Larmor scales, and discuss how we derive simplified dispersion relations for (low-frequency) CE microinstabilities more generally.

First, it is unnecessarily restrictive to assume that, for all CE microinstabilities, $\tilde{\omega}_{s\|} \ll 1$
for all particle species. There are some instabilities for which $\tilde{\omega}_{e\|} \sim \eta_e \sim \epsilon_e \ll 1$
while $\tilde{\omega}_{i\|} \gtrsim 1$. Recalling the orderings $\tilde{\omega}_{e\|} \sim \beta_e^{-1}$ and $k \rho_e \sim 1$ that were adopted for the electron-Larmor-scale instabilities described in section <ref>, it follows that $\tilde{\omega}_{i\|} \gtrsim 1$ whenever $\beta_e \lesssim \tau^{-1/2} \mu_e^{-1/2}$; in other words, electron-Larmor-scale CE microinstabilities in plasmas with $\beta_e$ that is not too large will satisfy $\tilde{\omega}_{i\|} \gtrsim 1$. Therefore, we cannot naively apply our low-frequency approximation to both $\boldsymbol{\mathfrak{E}}_e$ and $\boldsymbol{\mathfrak{E}}_i$
in all cases of interest. We will remain cognisant of this in the calculations
that follow – a concrete example of $\tilde{\omega}_{i\|} \gtrsim 1$ will
be considered in section <ref>.

Secondly, because of the large separation between electron and ion Larmor scales, it is
necessary to consider whether the approximation $\mathsfbi{M}_s\!\left(\tilde{\omega}_{s\|},\boldsymbol{k}\right) \approx \tilde{\omega}_{s\|} \mathsfbi{M}_s^{(0)}\!\left(\boldsymbol{k}\right)$
remains valid for parallel or perpendicular wavenumbers much larger or smaller than the inverse Larmor radii of
each species. We show in appendix <ref> that the leading-order term in the $\tilde{\omega}_{s\|} \ll 1$ expansion remains larger than higher-order terms for all $k_{\|} \rho_s \gtrsim 1$ (as, indeed, was implicitly assumed in section <ref>).
However, for $k_{\|} \rho_s$ sufficiently small, the same statement does not hold for all components of
$\mathsfbi{M}_s$. More specifically, it is shown in the same appendix that the dominant contribution to $\mathsfbi{M}_s\!\left(\boldsymbol{k}\right)$
when $k_{\|} \rho_s \ll 1$ instead comes from the quadratic term $\tilde{\omega}_{s\|}^2 \mathsfbi{M}_s^{(1)}\!\left(\boldsymbol{k}\right)$ (rather than any higher-order
term). Thus, in general, our simplified dispersion relation for low-frequency modes in a two-species plasma has the form of
a quartic in $\omega$, rather than a quadratic, if $k_{\|} \rho_s \ll 1$ for at least
the electron species. Physically, the reason why a
quadratic dispersion relation is no longer a reasonable approximation
is the existence of more than two low-frequency modes in a two-species Maxwellian plasma in certain
wavenumber regimes.
For example, for quasi-parallel modes with characteristic parallel wavenumbers satisfying $k_{\|} \rho_i \ll 1$,
there are four low-frequency modes (see section <ref>).
Nevertheless, in other situations, the components of $\mathsfbi{M}_s\!\left(\boldsymbol{k}\right)$
for which the $\mathsfbi{M}_s\!\left(\tilde{\omega}_{s\|},\boldsymbol{k}\right) \approx \tilde{\omega}_{s\|} \mathsfbi{M}_s^{(0)}\!\left(\boldsymbol{k}\right)$
approximation breaks down are not important, on account of their small size
compared with terms in the dispersion relation associated with other Maxwellian
components. In this case, the original quadratic dispersion relation is
sufficient. An explicit wavenumber regime in which this is realised is $k_{\|} \rho_e \sim k_{\perp} \rho_e \ll
1$ but $k \rho_i \gg 1$ – see sections <ref> and <ref>.

Taking these multiple-species effects into account, the reasons behind the decision made in section <ref> to
consider the CES microinstabilities separately from the CET microinstabilities
come into plain focus. First, the characteristic sizes of the CE
electron-temperature-gradient and ion-temperature-gradient terms are comparable ($\eta_i \sim