99affiliationtext: Centre for mathematical Plasma Astrophysics, KU Leuven,
3001 Leuven, Belgium

# Generalized Fluid Models of the Braginskii Type

P. Hunana11affiliation: Instituto de Astrofísica de Canarias (IAC), La Laguna,
Tenerife, 38205, Spain<EMAIL_ADDRESS>22affiliation: Universidad de
La Laguna, La Laguna, Tenerife, 38206, Spain , T. Passot33affiliation:
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire
Lagrange, Bd de l’Observatoire, CS 34229, 06304 Nice cedex 4, France , E.
Khomenko11affiliationmark:  22affiliationmark:  , D. Martínez-
Gómez11affiliationmark:  22affiliationmark:  , M. Collados11affiliationmark:
22affiliationmark:  ,
A. Tenerani44affiliation: Department of Physics, University of Texas, Austin,
TX 78712, USA , G. P. Zank55affiliation: Center for Space Plasma and Aeronomic
Research (CSPAR), University of Alabama, Huntsville, AL 35805, USA
66affiliation: Department of Space Science, University of Alabama, Huntsville,
AL 35899, USA , Y. Maneva77affiliation: Royal Observatory of Belgium, 1180
Brussels, Belgium , M. L. Goldstein88affiliation: Space Science Institute,
Boulder, CO 80301, USA , G. M. Webb55affiliationmark:

###### Abstract

Several generalizations of the well-known fluid model of Braginskii (Rev. of
Plasma Phys., 1965) are considered. We use the Landau collisional operator and
the moment method of Grad. We focus on the 21-moment model that is analogous
to the Braginskii model, and we also consider a 22-moment model. Both models
are formulated for general multi-species plasmas with arbitrary masses and
temperatures, where all the fluid moments are described by their evolution
equations. The 21-moment model contains two “heat flux vectors” (3rd and 5th-
order moments) and two “viscosity-tensors” (2nd and 4th-order moments). The
Braginskii model is then obtained as a particular case of a one ion-electron
plasma with similar temperatures, with de-coupled heat fluxes and viscosity-
tensors expressed in a quasi-static approximation. We provide all the
numerical values of the Braginskii model in a fully analytic form (together
with the 4th and 5th-order moments). For multi-species plasmas, the model
makes calculation of transport coefficients straightforward. Formulation in
fluid moments (instead of Hermite moments) is also suitable for implementation
into existing numerical codes. It is emphasized that it is the quasi-static
approximation which makes some Braginskii coefficients divergent in a weakly-
collisional regime. Importantly, we show that the heat fluxes and viscosity-
tensors are coupled even in the linear approximation, and that the fully
contracted (scalar) perturbations of the 4th-order moment, which are accounted
for in the 22-moment model, modify the energy exchange rates. We also provide
several Appendices, which can be useful as a guide for deriving the Braginskii
model with the moment method of Grad.

###### Contents

  1. 1 Introduction
  2. 2 Multi-fluid generalization of Braginskii (21-moment model)
    1. 2.1 Collisional contributions (arbitrary masses and temperatures)
    2. 2.2 Collisional contributions for small temperature differences
    3. 2.3 Reduction to 13-moment model
    4. 2.4 Semi-linear approximation (de-coupled stress tensors and heat fluxes)
  3. 3 One ion-electron plasma
    1. 3.1 Ion heat flux $\vec{\boldsymbol{q}}_{a}$ of Braginskii (self-collisions)
    2. 3.2 Ion heat flux $\vec{\boldsymbol{X}}^{(5)}_{a}$ (self-collisions)
    3. 3.3 Electron heat flux $\vec{\boldsymbol{q}}_{e}$ of Braginskii
    4. 3.4 Electron heat flux $\vec{\boldsymbol{X}}_{e}^{(5)}$
    5. 3.5 Ion viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ of Braginskii (self-collisions)
    6. 3.6 Ion viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ (self-collisions)
    7. 3.7 Electron viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}$ of Braginskii
    8. 3.8 Electron viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}$
  4. 4 Generalized electron coefficients for multi-species plasmas
  5. 5 Generalization with coupling of stress-tensors and heat fluxes
  6. 6 Coupling for unmagnetized one ion-electron plasma.
    1. 6.1 Ion species (self-collisions)
    2. 6.2 Electron species
    3. 6.3 Momentum exchange rates
  7. 7 Multi-fluid 22-moment model
    1. 7.1 Collisional contributions (arbitrary temperatures)
      1. 7.1.1 Small temperature differences
    2. 7.2 Ion species (self-collisions)
    3. 7.3 Electron species (one ion-electron plasma)
  8. 8 Discussion & Conclusions
    1. 8.1 Energy conservation
    2. 8.2 Collisional frequencies for ion-electron plasma
    3. 8.3 Fluid hierarchy
    4. 8.4 Reducible & irreducible Hermite polynomials
    5. 8.5 Rosenbluth potentials (22-moment model)
    6. 8.6 Hermite closures
    7. 8.7 Inclusion of gravity
    8. 8.8 Precision of $m_{e}/m_{i}$ expansions (unmagnetized proton-electron plasma)
    9. 8.9 Limitations of our approach
      1. 8.9.1 Weakly collisional regime: expansions around bi-Maxwellians
      2. 8.9.2 Landau fluid closures for the collisionless case
      3. 8.9.3 Ion-sound wave damping in homogeneous media: comparison of various models
      4. 8.9.4 Large gradients and large drifts
      5. 8.9.5 Comments on the positivity of the perturbed distribution function
    10. 8.10 Conclusions
  9. 9 Acknowledgments
  10. A General evolution equations
  11. B Tensorial Hermite Decomposition
    1. B.1 Usual perturbations of Grad
      1. B.1.1 20-moment model
      2. B.1.2 13-moment model
      3. B.1.3 Double-checking the fluid moments
    2. B.2 Higher-order perturbations (full $\bar{\bar{\boldsymbol{X}}}^{(4)}$ and $\bar{\bar{\boldsymbol{X}}}^{(5)}$ moments)
      1. B.2.1 Viscosity $\Pi^{(4)}_{ij}$ of the 4th-order moment $X_{ij}^{(4)}$
      2. B.2.2 Simplified perturbations (21-moment model)
    3. B.3 Irreducible Hermite polynomials
      1. B.3.1 Higher-order tensorial “anisotropies”
    4. B.4 Decomposition of $X^{(4)}_{ijkl}$
    5. B.5 Decomposition of $X^{(5)}_{ijklm}$
    6. B.6 Table of useful integrals
    7. B.7 General n-th order perturbation
    8. B.8 Hierarchy of MHD Hermite closures
    9. B.9 Hierarchy of CGL (parallel) Hermite closures
  12. C Evolution equations for 22-moment model
    1. C.1 Decomposition of moments
    2. C.2 Evolution equation for scalar pressure $p_{a}$
    3. C.3 Evolution equation for viscosity tensor $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$
    4. C.4 Evolution equation for heat flux vector $\vec{\boldsymbol{q}}_{a}$
    5. C.5 Evolution equation for viscosity-tensor $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$
    6. C.6 Evolution equation for perturbation $\widetilde{X}^{(4)}_{a}$
    7. C.7 Evolution equation for heat flux vector $\vec{\boldsymbol{X}}^{(5)}_{a}$
  13. D Simplified general fluid hierarchy
  14. E BGK collisional operator
    1. E.1 Viscosity-tensor $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$
    2. E.2 Heat flux vector $\vec{\boldsymbol{q}}_{a}$
    3. E.3 BGK vs Braginskii comparison
    4. E.4 Nonlinear stress-tensor decomposition
      1. E.4.1 Decomposition of a matrix
      2. E.4.2 Symmetric and traceless matrices
    5. E.5 BGK operator and electric field
    6. E.6 Ambipolar diffusion of two ion species
      1. E.6.1 Damping of Alfvén waves
  15. F General Fokker-Planck collisional operator
  16. G Landau collisional operator (5-moment model)
    1. G.1 Momentum exchange rates $\boldsymbol{R}_{ab}$
    2. G.2 Energy exchange rates $Q_{ab}$
    3. G.3 $R_{ab}$ and $Q_{ab}$ for unrestricted drifts $\boldsymbol{u}_{b}-\boldsymbol{u}_{a}$ (runaway effect)
    4. G.4 Difficulties with Rosenbluth potentials
  17. H 8-moment model (heat flux and thermal force)
    1. H.1 Rosenbluth potentials
    2. H.2 Momentum exchange rates $\boldsymbol{R}_{ab}$
    3. H.3 Heat flux exchange rates
    4. H.4 One ion-electron plasma
  18. I Comparison of various models with Braginskii (electrons)
    1. I.1 Notation of Spitzer-Härm 1953
    2. I.2 Model of Killie et al. 2004
  19. J 10-moment model (viscosity)
    1. J.1 Rosenbluth potentials
    2. J.2 Viscosity calculation
    3. J.3 Small temperature differences
    4. J.4 One ion-electron plasma
    5. J.5 Strong magnetic field limit
    6. J.6 Table of integrals
      1. J.6.1 Spherical integration
  20. K Braginskii heat flux (11-moment model)
    1. K.1 Rosenbluth potentials
    2. K.2 Dynamical friction vector and diffusion tensor
    3. K.3 Distribution function for species “a”
    4. K.4 Momentum exchange rates $R_{ab}$
    5. K.5 Heat flux exchange rates
    6. K.6 5th-order moment exchange rates
  21. L Braginskii viscosity (15-moment model)
    1. L.1 Rosenbluth potentials
    2. L.2 Dynamical friction vector and diffusion tensor
    3. L.3 Distribution function for species “a”
    4. L.4 Pressure tensor exchange rates
    5. L.5 Viscosity-tensor exhange rates
    6. L.6 4th-order moment exchange rates
    7. L.7 Exchange rates $\bar{\bar{\boldsymbol{Q}}}_{a}^{(4)}\,{}^{\prime}$
  22. M Collisional contributions for scalar $\widetilde{X}^{(4)}$
    1. M.1 Pressure tensor exchange rates
    2. M.2 4th-order moment exchange rates
    3. M.3 Collisional contributions $\widetilde{Q}^{(4)}_{a}\,{}^{\prime}$
  23. N Coupling of two species
    1. N.1 Protons and alpha particles (unmagnetized)
    2. N.2 Deuterium & tritium plasma (unmagnetized)

##  1 Introduction

The fluid model of Braginskii (1958, 1965) represents a cornerstone of plasma
transport theory and it is used in many different areas, from solar physics to
laboratory plasmas. The Braginskii model and its generalizations can be
derived through two major classical routes 1) Chapman-Enskog expansions
(Chapman & Cowling, 1939) and 2) the moment method of Grad (1949a, b, 1958).
There also exists a more modern route with the projection operator (Krommes,
2018a, b). Both classical routes were originally developed for gases, where
the full Boltzmann collisional operator has to be used. As was shown by Landau
(1936, 1937), for charged particles interacting through Coulomb collisions the
Boltzmann operator can be partially simplified, and this collisional operator
is known as the Landau operator. It is now well-established that for Coulomb
collisions both Landau and Boltzmann operators yield the same results, if in
the Boltzmann operator one introduces integration cut-offs that remove the
divergences in the same way as the Coulomb logarithm does. With the Landau
operator, the Boltzmann equation is then typically called the Landau equation.
By introducing Rosenbluth potentials the Landau operator can be re-written
into a general Fokker-Planck form, and the name Fokker-Planck equation is
often used as well. Nevertheless, many authors use the Boltzmann operator
during calculations even when Coulomb collisions are considered, because the
simplification is not exceedingly large. Braginskii used the Landau operator.
Of course, both routes through Chapman-Enskog expansions and the moment method
of Grad have sub-variations on how the methods are implemented that were
developed over the years. For the Chapman-Enskog method where the distribution
function is expanded in Laguerre-Sonine polynomials, see for example
Braginskii (1958); Hinton (1983); Helander & Sigmar (2002); Kunz (2021).

Here we use the moment method of Grad, which consists of expanding the
distribution function in tensorial Hermite polynomials. Concerning only
viscosity-tensors and heat fluxes (and neglecting fully contracted scalar
perturbations and higher-order tensorial “anisotropies” as Balescu (1988)
calls them), the method of Grad consists of approximating the distribution
function as a series

$f_{a}=f_{a}^{(0)}(1+\chi_{a});\qquad\chi_{a}=\sum_{n=1}^{N}\Big{[}{h}^{(2n)}_{ij}{H}^{(2n)}_{ij}+{h}^{(2n+1)}_{i}{H}^{(2n+1)}_{i}\Big{]},$
(1)

where $f_{a}^{(0)}$ is Maxwellian, “a” is species index, indices i and j run
from 1 to 3, $H$ are Hermite polynomials and $h$ are Hermite moments. Matrices
$h_{ij}^{(2n)}$ are traceless and can be viewed as viscosity-tensors (stress-
tensors), and vectors ${h}^{(2n+1)}_{i}$ can be viewed as heat fluxes. The
series is cut at some chosen $N$, and this distribution function is then used
in the Landau (or Boltzmann) equation, which is integrated to obtain a
corresponding fluid model. The usual quasi-static approximation does not have
to be applied, and one obtains evolution equations for all the considered
moments. For example, prescribing a strict Maxwellian with perturbation
$\chi_{a}=0$ (or equivalently $N=0$) represents the 5-moment model, with
evolution equations for density, fluid velocity and scalar pressure
(temperature), where stress-tensors and heat fluxes are zero. Prescribing
$N=1$ represents the 13-moment model, which contains evolution equation for
one traceless viscosity tensor (5 independent components) and an evolution
equation for one heat flux vector (3 independent components). This model thus
contains the main ingredients of the model of Braginskii, i.e. the usual
viscosity tensor and the usual heat flux vector are present. However,
prescribing a quasi-static approximation, yields for example the coefficient
of parallel electron heat conductivity (for a one-ion electron plasma with ion
charge $Z_{i}=1$) with value 1.34 instead of the Braginskii value 3.16,
meaning the model is not sufficiently precise. Prescribing $N=2$ represents
21-moment model, and this model can be viewed as containing evolution
equations for two viscosity tensors and two heat flux vectors. It can be shown
that expressing the viscosity tensors and heat fluxes in a quasi-static
approximation yields a model that is equivalent to Braginskii (1965). In fact,
as pointed out for example by Balescu (1988), the Hermite polynomials are
directly related to the Laguerre-Sonine polynomials; see equation (192), and
thus the Chapman-Enskog method and the moment method of Grad have to yield
equivalent results at the end. In general, if both heat fluxes and viscosities
are accounted for, an $N$-Laguerre model therefore represents a
$(5+8N)$-moment model. For a summary of various possible models, see Section
8.4 with Tables 8.1 & 8.2.

Of course, the model of Braginskii can be generalized in many different ways.
Naturally, one might focus at the case of one-ion electron plasma considered
by Braginskii, and increase the order of $N$ to study convergence of transport
coefficients with higher-order Laguerre (Hermite) schemes. Several studies of
this kind were done in the past (some numerically imprecise, some considering
only unmagnetized plasmas, and some only ion charge $Z_{i}=1$). For example
before Braginskii, Landshoff (1949, 1951) calculated several transport
coefficients with models from $N=1$ to $N=4$. Kaneko (1960) improved the
numerical accuracy of Landshoff and also considered $N=5$. Kaneko & Taguchi
(1978); Kaneko & Yamao (1980) performed large calculations with up to a
$N=49$. Perhaps the most comprehensive study to this date was done by Ji &
Held (2013), who studied convergence of all the transport coefficients with up
to $N=160$. Other useful references can be found in Epperlein & Haines (1986).
The last two studies emphasize that while the transport coefficients parallel
to the magnetic field (or for unmagnetized plasma) converge rapidly for $N\geq
2$, this is not the case for some perpendicular transport coefficients. For
clarity, in the famous work of Spitzer & Härm (1953) and the previous work of
Cohen et al. (1950) where only unmagnetized plasma is considered and
viscosity-tensors are neglected, the perturbation $\chi_{a}$ which satisfies
the Landau equation was found numerically, and the model thus technically
corresponds to $N=\infty$. Their work is criticized (even though a bit
unfairly) in the monograph of Balescu (1988) (Part 1, p. 266), who calculated
all the usual transport coefficients with the moment method of Grad for the
$N=2$ and $N=3$ cases (i.e. the 21-moment model and the 29-moment model). Note
that the 3-Laguerre calculations of Balescu (1988) were shown to be incorrect
by Ji & Held (2013), who were able to trace the problem to his analytic
collisional matrices (they have also corrected coefficients in collisional
matrices of Braginskii (1958), which were fortunately not used in his $N=2$
calculation). That there is a problem with the $N=3$ transport coefficients of
Balescu (1988) can be also seen for example by comparison with Kaneko (1960).
Here we focus at the 2-Laguerre approximation used by Braginskii (1965), i.e.
the 21-moment model, with the goal to extract more physical information from
that scheme.

For the 5-moment model and the 13-moment model, the method of Grad was
explored in great detail by Burgers (1969) and Schunk (1975, 1977) (see other
references therein). The Boltzmann operator was used and several interaction
potentials were considered, such as collisions between neutral particles (hard
sphere interaction), between charges (Coulomb interaction), or an induced
dipole interaction when an ion polarizes a colliding neutral (so called
Maxwell molecule interaction). These models have two important properties that
the Braginskii model does not have: 1) because the formulation is with
evolution equations for stress-tensors and heat fluxes rather than with quasi-
static approximation, these models do not become divergent if a regime of low
collisionality is encountered; 2) the formulation is as a general multi-fluid
description with arbitrary masses $m_{a}$, $m_{b}$ and temperatures $T_{a}$,
$T_{b}$. Note that the review paper of Braginskii (1965) also contains Section
7 about multi-component plasmas that is often implicitly cited in the solar
literature, but this section should be viewed as heuristic from a perspective
that no heat fluxes or stress-tensors were calculated. In plasma physics, the
work of Braginskii (1958, 1965) is celebrated for his results for a one ion-
electron plasma. Here we use the Landau operator and consider only Coulomb
collisions. Nevertheless, we will employ the 21-moment model, and we thus
improve the precision of the 13-moment model of Burgers (1969)-Schunk (1977)
for this interaction potential, so that the precision matches Braginskii. We
will use a restriction that the relative drift velocity between two colliding
species must be small in comparison to their thermal speeds. The same
restriction applies for the Braginskii model, for the Burgers-Schunk 13-moment
model (the exception is the Maxwell molecule interaction) and higher-order
schemes. For Coulomb collisions and hard-sphere collisions, only the simplest
5-moment model has been calculated fully analytically without this restriction
(Burgers, 1969; Tanenbaum, 1967; Schunk, 1977), yielding the runaway effect.

Several various multi-fluid descriptions with the level of precision of
Braginskii have been considered in the past; see for example Hinton (1983),
Zhdanov (2002) (orig. pub. 1982), Ji & Held (2006) (who actually consider
general $N$), Simakov & Molvig (2014, 2016a, 2016b), or for the case of
neoclassical theory (toroidal geometry applicable to tokamaks) see Hirshman &
Sigmar (1977, 1981). Our model seems to be very close to the model of Zhdanov
(2002), Chapter 8.1, who indeed uses the method of Grad and calculates the
21-moment model with it. We did not verify full equivalence because of his
puzzling notation. Even if equivalence is eventually shown for the case of
small temperature differences between ions, we consider a more general case
where temperatures of all the species are arbitrary. Our clear formulation
with fluid moments (instead of Hermite moments) might be also easier to
implement into existing numerical codes. Arbitrary temperatures were also
considered by Ji & Held (2006), but we did not verify equivalence with their
model either. We only verified equivalence with their model for the special
case of a one ion-electron plasma with small temperature differences of
Braginskii, by using collisional matrices from Ji & Held (2013).

Additionally, for all the considered moments we provide the left-hand-sides of
our evolution equations in a fully non-linear form, which is important for
direct numerical simulations and which are not typically given. An important
difference then arises even at the linear level, because calculations are
typically performed with de-coupled viscosity-tensors and heat fluxes, meaning
that the two viscosity-tensors interact only with each other, and the two heat
fluxes interact only with each other. We consider coupling between heat fluxes
and stress-tensors, where (even at the linear level in a quasi-static
approximation) a heat flux enters a stress-tensor and a stress-tensor enters a
heat flux. Such couplings are often considered in the collisionless regime;
see e.g. Macmahon (1965),Mikhailovskii & Smolyakov (1985), Ramos (2005),
Goswami et al. (2005), Passot et al. (2012), Hunana et al. (2019a, b), where
the effect is important for example for the perpendicular fast mode, or for
the growth rate of the firehose instability (see e.g. Figure 10 in Hunana et
al. (2019b)). The coupling might be important also in the highly-collisional
regime if sufficiently high frequencies (or short wavelengths) are considered.
The coupling was neglected by Braginskii (1958, 1965) and Spitzer & Härm
(1953); Spitzer (1962), and as an example we consider unmagnetized one ion-
electron plasma in detail, and we provide stress-tensors and heat fluxes where
this coupling is taken into account.

The coupling between viscosity-tensors and heat fluxes then inevitably leads
to the next step, by replacing equation (1) with

$f_{a}=f_{a}^{(0)}(1+\chi_{a});\qquad\chi_{a}=\sum_{n=1}^{N}\Big{[}{h}^{(2n)}_{ij}{H}^{(2n)}_{ij}+{h}^{(2n+1)}_{i}{H}^{(2n+1)}_{i}+h^{(2n)}H^{(2n)}\Big{]},$
(2)

where the scalar hermite moments $h^{(2n)}$ can be viewed as fully contracted
(scalar) perturbations of fluid moments. The lowest-order moment $h^{(2)}=0$
and all higher-order ones are generally non-zero. Thus, prescribing $N=1$
still yields the 13-moment model, however, precribing $N=2$ now yields the
22-moment model. This model is a natural extension of the Braginskii model,
because it takes into account fully contracted perturbations
$\widetilde{X}^{(4)}_{a}=m_{a}\int|\boldsymbol{c}_{a}|^{4}(f_{a}-f_{a}^{(0)})d^{3}v$
of the 4th-order fluid moment. Accounting for the scalar perturbations
according to (2), for $N\geq 1$ an $N$-Laguerre model then represents a
$(4+9N)$-moment model. Another possibility for writing the equation (2) is to
separate the matrices $\sum_{n=1}^{N}{h}^{(2n)}_{ij}{H}^{(2n)}_{ij}$, and for
the vectors and scalars to write the sum from $n=0$, with an imposed
requirement that $h^{(0)}=0$; $h^{(2)}=0$ & $h^{(1)}_{i}=0$ (where the first
one is non-trivial). This is for example the choice of Balescu (1988), p. 174,
his equations (3.11) & (3.16).

Finally, the main purpose of this work is to make the moment method of Grad
and the exciting work of Braginskii more understandable, which is reflected in
our relatively lengthy Appendix.

The entire paper is separated into eight Sections and fourteen Appendices. The
main paper summarizes the obtained results, while the Appendices provide the
detailed calculations.

In Section 2, we formulate the entire 21-moment model. We start with a
formulation valid for a general collisional operator $C(f_{a})$, where both
the left-hand-sides and the collisional right-hand-sides of evolution
equations are given in a fully non-linear form. We then provide collisional
contributions for arbitrary masses and temperatures calculated with the Landau
operator. Collisional contributions are calculated in the usual semi-linear
approximation where relative drifts between species are small in comparison to
their thermal speeds (i.e. the runaway effect is not considered), and the
product of $f_{a}f_{b}$ is approximated as
$f_{a}f_{b}=f_{a}^{(0)}f_{b}^{(0)}(1+\chi_{a}+\chi_{b})$, where the
“cross”-contributions $\chi_{a}\chi_{b}$ are neglected. We then provide a
simplified model where differences in temperatures between species are small.
For clarity, we also reduce our model to the 13-moment model and we provide a
formulation that is more compact than the one given by Burgers (1969)-Schunk
(1977) (because we only consider Coulomb collisions). We then simplify the
evolution equations of our 21-moment model into a semi-linear approximation
where viscosity-tensors and heat fluxes are de-coupled, and these are used in
Sections 3 & 4.

In Section 3, we compare our model to Braginskii (1965) by considering a one
ion electron plasma with similar temperatures, i.e. where the temperature
differences between species are small with respect to their mean values. We
provide all the transport coefficients in a fully analytic form, and we verify
the entire Table II of Braginskii (1965) (two of his coefficients are not
precise). Parallel electron coefficients (or equivalently for an unmagnetized
plasma), can be also found in Simakov & Molvig (2014). We also provide
analytic results for the viscosity of the 4th-order fluid moment and the heat
flux of the 5th-order fluid moment, which are not typically given.

In Section 4, we use the idea of for example Hinton (1983); Zhdanov (2002);
Simakov & Molvig (2014) that because of the smallness of electron/ion mass
ratios, the _electron_ coefficients of Braginskii can be straightforwardly
generalized to multiple ion species by introducing an effective ion charge and
effective ion velocity. All the electron analytic coefficients that are given
in Section 3 are thus generalized to multi-ion species with a simple
transformation.

In Section 5, we discuss the coupling between viscosity tensors and heat
fluxes. We provide evolution equations in the semi-linear approximation where
this coupling is retained, and we introduce a technique on how to split the
moments into their first and second orders.

In Section 6, we consider an example of an unmagnetized one ion-electron
plasma and explicitly calculate the coupling of stress-tensors and heat
fluxes. All the results are given in a fully analytic form, as well as with
numerical values for the ion charge $Z_{i}=1$.

In Section 7, we first formulate the fully non-linear 22-moment model for a
general collisional operator. We then provide the multi-fluid collisional
contributions calculated with the Landau operator in the semi-linear
approximation, and we show that the perturbations $\widetilde{X}^{(4)}$ modify
the energy exchange rates. We also provide quasi-static solutions for a one
ion-electron plasma, and we show that the perturbations $\widetilde{X}^{(4)}$
have their own heat conductivities.

In Section 8, we discuss various topics. 1) We discuss energy conservation. 2)
We clarify that from a multi-fluid perspective the Braginskii choice of ion
collisional time $\tau_{i}$ should be interpreted as $\tau_{i}=\tau_{ii}$, and
not as $\tau_{i}=\sqrt{2}\tau_{ii}$. 3) To clarify the higher-order schemes
and to double-check our evolution equations, we calculate the fluid hierarchy
for a general $N$, with an unspecified collisional operator. 4) We discuss
irreducible and reducible Hermite polynomials and show that both yield the
same results. 5) We provide fully non-linear Rosenbluth potentials for the
22-moment model, which might be useful in further studies of the runaway
effect with this scheme. 6) We discuss Hermite closures and their relation to
fluid closures, which are required to close the fluid hierarchy. We also
correct our previous erroneous interpretation that Landau fluid closures are
necessary to go beyond the 4th-order moment. 7) We discuss the inclusion of
gravity. 8) We use our multi-fluid formulation to double-check the precision
of $m_{e}/m_{i}$ expansions. We consider unmagnetized proton-electron plasma,
and calculate the transport coefficients exactly, without using the smallness
of $m_{e}/m_{i}$. 9) We discuss limitations of our approach. 10) We provide
conclusions, with examples where our model might be useful.

Appendix A introduces the general concept of tensorial fluid moments and
provides an evolution equation for an n-th-order fluid moment
$\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}$ in presence of a general (unspecified)
collisional operator, equation (A12). This evolution equation also remains
valid in presence of gravity; see the discussion in Section 8.7.

Appendix B introduces the tensorial Hermite polynomials of Grad (1949a, b,
1958), and discusses in detail the construction of perturbations around the
Maxwellian distribution function, i.e. equations (1) & (2), which are
summarized in Tables 8.1 and 8.2. The construction of Hermite closures is
addressed as well.

Appendix C derives evolution equations for the 22-moment model (for an
unspecified collisional operator), by applying contractions at the evolution
equations from Appendix A and by using decomposition of moments and Hermite
closures from Appendix B.

Appendix D uses a different technique and instead of applying contractions at
equations of Appendix A, simplified fluid hierarchy of a general n-th-order is
obtained directly, which only consists of evolution equations for scalars,
vectors and matrices. Evaluation of these equations for a specific “n”
recovers the 22-moment equations of Appendix C.

Appendix E introduces the BGK (relaxation-type) collisional operator of
Bhatnagar et al. (1954); Gross & Krook (1956), which greatly clarifies the
analytic forms of the Braginskii viscosity-tensors and heat fluxes.
Viscosities and heat conductivities of both models are directly compared in
Figures 3-5. The nonlinear solution for the viscosity-tensor (with respect to
a general direction of magnetic field $\boldsymbol{\hat{b}}$) is addressed in
Appendix E.4, and Appendix E.6 clarifies the ambipolar diffusion between two
ion species.

Appendix F introduces a general (unspecified) Fokker-Planck collisional
operator with its dynamical friction vector $\boldsymbol{A}_{ab}$ and
diffusion tensor $\bar{\bar{\boldsymbol{D}}}_{ab}$. General relations for the
collisional integrals (of n-th-order) are provided, which can be used once the
$\boldsymbol{A}_{ab}$ & $\bar{\bar{\boldsymbol{D}}}_{ab}$ are specified.

Appendix G introduces the Landau collisional operator, where the
$\boldsymbol{A}_{ab}$ & $\bar{\bar{\boldsymbol{D}}}_{ab}$ are expressed in the
usual form through the Rosenbluth potentials. The 5-moment model (strict
Maxwellians) is then considered, and the usual collisional momentum exchange
rates $\boldsymbol{R}_{ab}$ and energy exchange rates $Q_{ab}$ with the
assumption of small drifts between species are derived in detail in Appendices
G.1 & G.2. Both contributions are then re-calculated with unrestricted drifts
in Appendix G.3, where instead of the Rosenbluth potentials, the “center-of-
mass” transformation typically used with the Boltzmann collisional operator
has to be used, because the collisional integrals seem to be too complicated
to calculate directly. This is further discussed in Appendix G.4.

Appendix H considers the 8-moment model, where the simplest heat flux is
present, and the multi-fluid model of Burgers (1969)-Schunk (1977) is
calculated in detail. For a direct comparison with Braginskii, a one ion-
electron plasma is then considered and quasi-static heat fluxes, together with
the resulting momentum exchange rates are obtained as well. It is shown that
in the limit of strong magnetic field, the perpendicular and cross-
conductivities $\kappa_{\perp}$ & $\kappa_{\times}$ match the Braginskii model
exactly (for both the ion and electron species) and only the parallel
conductivities $\kappa_{\parallel}$ are different.

Appendix I compares the parallel heat fluxes and momentum exchange rates of
Braginskii (1965) with models of Burgers (1969)-Schunk (1977), Killie et al.
(2004), Landshoff (1949, 1951) and Spitzer & Härm (1953), see Tables I.1-I.5.
Useful conversion relations for the results of Kaneko (1960) and Balescu
(1988) are provided as well. The notation of Spitzer & Härm (1953) is
clarified in Appendix I.1 and it is shown that their model, as well as the
model of Killie et al. (2004), break the Onsager symmetry.

Appendix J calculates in detail the 10-moment multi-fluid model of Burgers
(1969)-Schunk (1977), where the simplest viscosity-tensor is present. It is
shown that in the limit of strong magnetic field, the perpendicular
viscosities and gyroviscosities $\eta_{1},\eta_{2},\eta_{3},\eta_{4}$ match
the Braginskii model exactly (for both the ion and electron species) and only
the parallel viscosities $\eta_{0}$ are different.

Appendix K calculates in detail the momentum exchange rates and collisional
contributions for the heat fluxes in our 21 & 22-moment multi-fluid models.
The calculations are shown on the 11-moment model where only the heat fluxes
are present (and viscosities and scalar perturbation are absent), because in
the semi-linear approximation the calculations can be split. Similarly,
collisional contributions for viscosity-tensors are calculated in Appendix L,
and contributions for the scalar perturbation of the 4th-order moment in
Appendix M.

Appendix N uses our 21-moment model and calculates heat conductivities and
viscosities for two examples of an unmagnetized plasma consisting of two ion
species (collisions with electrons are neglected). The first example (Appendix
N.1) is a plasma consisting of protons and alpha-particles (fully ionized
Helium), typical in astrophysical applications. The second example (Appendix
N.2) is a deuterium-tritium plasma used in plasma fusion.

##  2 Multi-fluid generalization of Braginskii (21-moment model)

Our model is formulated with heat flux vectors

$\displaystyle\vec{\boldsymbol{X}}^{(3)}_{a}=m_{a}\int\boldsymbol{c}_{a}|\boldsymbol{c}_{a}|^{2}f_{a}d^{3}v=2\vec{\boldsymbol{q}}_{a};\qquad\vec{\boldsymbol{X}}^{(5)}_{a}=m_{a}\int\boldsymbol{c}_{a}|\boldsymbol{c}_{a}|^{4}f_{a}d^{3}v,$
(3)

and traceless viscosity tensors

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ $\displaystyle=$
$\displaystyle
m_{a}\int\big{(}\boldsymbol{c}_{a}\boldsymbol{c}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}|\boldsymbol{c}_{a}|^{2}\big{)}f_{a}d^{3}v;\qquad\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}=m_{a}\int\big{(}\boldsymbol{c}_{a}\boldsymbol{c}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}|\boldsymbol{c}_{a}|^{2}\big{)}|\boldsymbol{c}_{a}|^{2}f_{a}d^{3}v,$
(4)

where the fluctuating velocity
$\boldsymbol{c}_{a}=\boldsymbol{v}-\boldsymbol{u}_{a}$, and “a” is the species
index. We are using free wording because $\vec{\boldsymbol{X}}^{(5)}_{a}$ is
not really a heat flux and $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ is not
really a viscosity-tensor. Also, we use the wording viscosity-tensor and
stress-tensor interchangeably in the entire text. The species indices are
moved freely up and down. We also define the usual rate-of-strain tensor
$\bar{\bar{\boldsymbol{W}}}_{a}=(\nabla\boldsymbol{u}_{a})^{S}-(2/3)\bar{\bar{\boldsymbol{I}}}\nabla\cdot\boldsymbol{u}_{a}$,
symmetric operator $A_{ij}^{S}=A_{ij}+A_{ji}$, and gravitational acceleration
$\boldsymbol{G}$. All other definitions are addressed in Appendix A. We note
that the definition of heat flux in equation (1.21) of Braginskii (1965)
contains two well-known misprints with prime symbols missing on his
fluctuating velocities $\boldsymbol{v}^{\prime}$. The heat flux is defined
correctly in Braginskii (1958).

We first present a formulation with a general (unspecified) collisional
operator $C(f_{a})$. We define (tensorial) collisional contributions

$\displaystyle\boldsymbol{R}_{a}=m_{a}\int\boldsymbol{v}C(f_{a})d^{3}v;\qquad
Q_{a}=\frac{m_{a}}{2}\int|\boldsymbol{c}_{a}|^{2}C(f_{a})d^{3}v;$
$\displaystyle\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}C(f_{a})d^{3}v;\qquad\bar{\bar{\boldsymbol{Q}}}^{(3)}_{a}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}C(f_{a})d^{3}v;$
$\displaystyle\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}C(f_{a})d^{3}v;\qquad\bar{\bar{\boldsymbol{Q}}}^{(5)}_{a}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}C(f_{a})d^{3}v,$
(5)

where $\boldsymbol{R}_{a}$ are the usual momentum exchange rates and $Q_{a}$
are the usual energy exchange rates. Then it can be shown that integration of
the Boltzmann equation yields the following nonlinear 21-moment model (see
details in Appendix C), where the basic evolution equations read

$\displaystyle\frac{d_{a}}{dt}n_{a}+n_{a}\nabla\cdot\boldsymbol{u}_{a}=0;$ (6)
$\displaystyle\frac{d_{a}}{dt}\boldsymbol{u}_{a}+\frac{1}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a}-\boldsymbol{G}-\frac{eZ_{a}}{m_{a}}\Big{(}\boldsymbol{E}+\frac{1}{c}\boldsymbol{u}_{a}\times\boldsymbol{B}\Big{)}=\frac{\boldsymbol{R}_{a}}{\rho_{a}};$
(7)
$\displaystyle\frac{d_{a}}{dt}p_{a}+\frac{5}{3}p_{a}\nabla\cdot\boldsymbol{u}_{a}+\frac{2}{3}\nabla\cdot\vec{\boldsymbol{q}}_{a}+\frac{2}{3}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}:(\nabla\boldsymbol{u}_{a})=\frac{2}{3}Q_{a},$
(8)

and are accompanied by evolution equations for the stress-tensors and heat
flux vectors

$\displaystyle\frac{d_{a}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}}{dt}+\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\nabla\cdot\boldsymbol{u}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}^{S}+\big{(}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\cdot\nabla\boldsymbol{u}_{a}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}(\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}:\nabla\boldsymbol{u}_{a})$
$\displaystyle\qquad+\frac{2}{5}\Big{[}(\nabla\vec{\boldsymbol{q}}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\vec{\boldsymbol{q}}_{a}\Big{]}+p_{a}\bar{\bar{\boldsymbol{W}}}_{a}=\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a}\,{}^{\prime}\equiv\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a};$
(9)
$\displaystyle\frac{d_{a}\vec{\boldsymbol{q}}_{a}}{dt}+\frac{7}{5}\vec{\boldsymbol{q}}_{a}\nabla\cdot\boldsymbol{u}_{a}+\frac{7}{5}\vec{\boldsymbol{q}}_{a}\cdot\nabla\boldsymbol{u}_{a}+\frac{2}{5}(\nabla\boldsymbol{u}_{a})\cdot\vec{\boldsymbol{q}}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{a}+\frac{5}{2}p_{a}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}$
$\displaystyle+\frac{1}{2}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-\frac{5}{2}\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}-\frac{1}{\rho_{a}}(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$
$\displaystyle\qquad=\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}\equiv\frac{1}{2}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(3)}_{a}-\frac{5}{2}\frac{p_{a}}{\rho_{a}}\boldsymbol{R}_{a}-\frac{1}{\rho_{a}}\boldsymbol{R}_{a}\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a};$
(10)
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\frac{1}{5}\Big{[}(\nabla\vec{\boldsymbol{X}}^{(5)}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}(\nabla\cdot\vec{\boldsymbol{X}}^{(5)}_{a})\Big{]}+\frac{9}{7}(\nabla\cdot\boldsymbol{u}_{a})\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\frac{9}{7}(\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\cdot\nabla\boldsymbol{u}_{a})^{S}$
$\displaystyle+\frac{2}{7}\big{(}(\nabla\boldsymbol{u}_{a})\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}-\frac{22}{21}\bar{\bar{\boldsymbol{I}}}(\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}:\nabla\boldsymbol{u}_{a})+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+7\frac{p_{a}^{2}}{\rho_{a}}\bar{\bar{\boldsymbol{W}}}_{a}$
$\displaystyle-\,\frac{14}{5\rho_{a}}\Big{[}\big{(}(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})\vec{\boldsymbol{q}}_{a}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})\cdot\vec{\boldsymbol{q}}_{a}\Big{]}$
$\displaystyle=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}\,{}^{\prime}\equiv\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}\textrm{Tr}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}-\frac{14}{5\rho_{a}}\Big{[}(\boldsymbol{R}_{a}\vec{\boldsymbol{q}}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}(\boldsymbol{R}_{a}\cdot\vec{\boldsymbol{q}}_{a})\Big{]};$
(11)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5)}_{a}+\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(6)}_{a}+\frac{9}{5}\vec{\boldsymbol{X}}^{(5)}_{a}(\nabla\cdot\boldsymbol{u}_{a})+\frac{9}{5}\vec{\boldsymbol{X}}^{(5)}_{a}\cdot\nabla\boldsymbol{u}_{a}+\frac{4}{5}(\nabla\boldsymbol{u}_{a})\cdot\vec{\boldsymbol{X}}^{(5)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{a}$
$\displaystyle+70\frac{p_{a}^{2}}{\rho_{a}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}-35\frac{p_{a}^{2}}{\rho_{a}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}-\frac{4}{\rho_{a}}\big{(}\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a}\big{)}\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$
$\displaystyle=\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime}\equiv\textrm{Tr}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(5)}_{a}-35\frac{p_{a}^{2}}{\rho_{a}^{2}}\boldsymbol{R}_{a}-\frac{4}{\rho_{a}}\boldsymbol{R}_{a}\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}.$
(12)

The last equation is closed with a fluid closure (derived from a Hermite
closure)

$\bar{\bar{\boldsymbol{\Pi}}}^{(6)}_{a}=m_{a}\int\big{(}\boldsymbol{c}_{a}\boldsymbol{c}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}|\boldsymbol{c}_{a}|^{2}\big{)}|\boldsymbol{c}_{a}|^{4}f_{a}d^{3}v=18\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-63\frac{p_{a}^{2}}{\rho_{a}^{2}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}.$
(13)

The system above thus represents a generalized model of Braginskii (1965),
where evolution equations for all the moments are fully non-linear and valid
for a general collisional operator $C(f_{a})$. It is a 21-moment model (1
density, 3 velocity, 1 scalar pressure, 3 for each heat flux vector, 5 for
each viscosity tensor).

###  2.1 Collisional contributions (arbitrary masses and temperatures)

We use the Landau collisional operator. All collisional contributions are
evaluated in a semi-linear approximation, with an assumption that differences
in drift velocities $\boldsymbol{u}_{b}-\boldsymbol{u}_{a}$ are small with
respect to thermal velocities. All the non-linear quantities such as
$\vec{\boldsymbol{q}}_{a}\cdot(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})$,
including $|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}$ are thus neglected in
the multi-fluid description, which is consistent with models of Burgers (1969)
and Schunk (1977). For energy conservation and a particular case of a one ion-
electron plasma, see Section 8.1. The wording semi-linear just means that
expressions containing pressures and densities such as
$(p_{a}/\rho_{a})\vec{\boldsymbol{q}}_{a}$ are retained and not fully
linearized with their mean pressure/density values. However, for example the
last terms of collisional contributions in the equations (10), (11), (12)
proportional to $\boldsymbol{R}_{a}\vec{\boldsymbol{q}}_{a}$,
$\boldsymbol{R}_{a}\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ and
$\boldsymbol{R}_{a}\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ are neglected
in the semi-linear approximation.

We introduce the usual reduced mass and reduced temperature

$\mu_{ab}=\frac{m_{a}m_{b}}{m_{a}+m_{b}};\qquad
T_{ab}=\frac{m_{a}T_{b}+m_{b}T_{a}}{m_{a}+m_{b}},$ (14)

together with collisional frequency (178). The momentum exchange rates are
given by

$\displaystyle\boldsymbol{R}_{a}$ $\displaystyle=$ $\displaystyle\sum_{b\neq
a}\nu_{ab}\Big{\\{}\rho_{a}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})+\frac{\mu_{ab}}{T_{ab}}\Big{[}V_{ab(1)}\vec{\boldsymbol{q}}_{a}-V_{ab(2)}\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{q}}_{b}\Big{]}$
(15)
$\displaystyle-\frac{3}{56}\Big{(}\frac{\mu_{ab}}{T_{ab}}\Big{)}^{2}\Big{[}\vec{\boldsymbol{X}}^{(5)}_{a}-\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{X}}^{(5)}_{b}\Big{]}\Big{\\}},$

with coefficients that include both masses and temperatures, but which we
simply call “mass-ratio coefficients”

$\displaystyle V_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{(21/10)T_{a}m_{b}+(3/5)T_{b}m_{a}}{T_{a}m_{b}+T_{b}m_{a}};\qquad
V_{ab(2)}=\frac{(3/5)T_{a}m_{b}+(21/10)T_{b}m_{a}}{T_{a}m_{b}+T_{b}m_{a}}.$
(16)

These and other mass-ratio coefficients given below come from the Landau
collisional operator introduced in Appendices F & G, where one uses perturbed
distribution functions of the 21-moment model; see Section 8.4 and Appendix B,
with calculations of the collisional integrals in Appendices K & L. Energy
exchange rates are given by

$Q_{a}=\sum_{b\neq a}3\rho_{a}\nu_{ab}\frac{T_{b}-T_{a}}{m_{a}+m_{b}},$ (17)

where $|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}$ are neglected as discussed
above. The heat flux exchange rates are given by

$\displaystyle\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}$ $\displaystyle=$
$\displaystyle-\Big{[}2\nu_{aa}+\sum_{b\neq
a}\nu_{ab}\hat{D}_{ab(1)}\Big{]}\vec{\boldsymbol{q}}_{a}+\sum_{b\neq
a}\nu_{ab}\hat{D}_{ab(2)}\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{q}}_{b}$
(18) $\displaystyle+\Big{[}\frac{3}{70}\nu_{aa}+\sum_{b\neq
a}\nu_{ab}\hat{E}_{ab(1)}\Big{]}\frac{\rho_{a}}{p_{a}}\vec{\boldsymbol{X}}^{(5)}_{a}-\sum_{b\neq
a}\nu_{ab}\hat{E}_{ab(2)}\frac{\rho_{b}}{p_{b}}\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{X}}^{(5)}_{b}$
$\displaystyle-p_{a}\sum_{b\neq
a}\nu_{ab}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\hat{U}_{ab(1)},$

with mass-ratio coefficients

$\displaystyle\hat{U}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{3m_{b}(3T_{a}m_{a}+T_{a}m_{b}-2T_{b}m_{a})}{2(T_{a}m_{b}+T_{b}m_{a})(m_{a}+m_{b})};$
$\displaystyle\hat{D}_{ab(1)}$ $\displaystyle=$
$\displaystyle\big{\\{}75T_{a}^{3}m_{a}m_{b}^{3}+95T_{a}^{3}m_{b}^{4}+174T_{a}^{2}T_{b}m_{a}m_{b}^{3}+300T_{a}T_{b}^{2}m_{a}^{3}m_{b}+498T_{a}T_{b}^{2}m_{a}^{2}m_{b}^{2}+60T_{b}^{3}m_{a}^{4}$
$\displaystyle+104T_{b}^{3}m_{a}^{3}m_{b}\big{\\}}\big{[}20(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{a}+m_{b})\big{]}^{-1};$
$\displaystyle\hat{D}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{9T_{a}m_{b}^{2}(10T_{a}^{2}m_{a}m_{b}+6T_{a}^{2}m_{b}^{2}+45T_{a}T_{b}m_{a}^{2}+27T_{a}T_{b}m_{a}m_{b}-14T_{b}^{2}m_{a}^{2})}{20(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{a}+m_{b})};$
$\displaystyle\hat{E}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{3T_{a}m_{b}(19T_{a}^{2}m_{a}m_{b}^{2}+23T_{a}^{2}m_{b}^{3}-2T_{a}T_{b}m_{a}^{2}m_{b}+36T_{a}T_{b}m_{a}m_{b}^{2}+84T_{b}^{2}m_{a}^{3}+118T_{b}^{2}m_{a}^{2}m_{b})}{560(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{a}+m_{b})};$
$\displaystyle\hat{E}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{9T_{a}T_{b}m_{a}m_{b}^{2}(7T_{a}m_{a}+5T_{a}m_{b}-2T_{b}m_{a})}{112(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{a}+m_{b})}.$
(19)

The 5th-order moment exchange rates are given by

$\displaystyle\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime}$ $\displaystyle=$
$\displaystyle-\Big{[}\frac{76}{5}\nu_{aa}+\sum_{b\neq
a}\nu_{ab}\hat{F}_{ab(1)}\Big{]}\frac{p_{a}}{\rho_{a}}\vec{\boldsymbol{q}}_{a}+\sum_{b\neq
a}\nu_{ab}\hat{F}_{ab(2)}\frac{p_{a}}{\rho_{a}}\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{q}}_{b}$
(20) $\displaystyle-\Big{[}\frac{3}{35}\nu_{aa}+\sum_{b\neq
a}\nu_{ab}\hat{G}_{ab(1)}\Big{]}\vec{\boldsymbol{X}}^{(5)}_{a}-\sum_{b\neq
a}\nu_{ab}\hat{G}_{ab(2)}\frac{p_{a}}{p_{b}}\vec{\boldsymbol{X}}^{(5)}_{b}$
$\displaystyle-\frac{p_{a}^{2}}{\rho_{a}}\sum_{b\neq
a}\nu_{ab}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\hat{U}_{ab(2)},$

with mass-ratio coefficients

$\displaystyle\hat{U}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{3m_{b}(17T_{a}^{2}m_{a}m_{b}+9T_{a}^{2}m_{b}^{2}+42T_{a}T_{b}m_{a}^{2}+6T_{a}T_{b}m_{a}m_{b}-28T_{b}^{2}m_{a}^{2})}{(T_{a}m_{b}+T_{b}m_{a})^{2}(m_{a}+m_{b})};$
$\displaystyle\hat{F}_{ab(1)}$ $\displaystyle=$
$\displaystyle\big{\\{}855T_{a}^{5}m_{a}m_{b}^{4}+759T_{a}^{5}m_{b}^{5}+2340T_{a}^{4}T_{b}m_{a}^{2}m_{b}^{3}+1972T_{a}^{4}T_{b}m_{a}m_{b}^{4}+2640T_{a}^{3}T_{b}^{2}m_{a}^{3}m_{b}^{2}$
$\displaystyle+2332T_{a}^{3}T_{b}^{2}m_{a}^{2}m_{b}^{3}+5880T_{a}^{2}T_{b}^{3}m_{a}^{4}m_{b}+3324T_{a}^{2}T_{b}^{3}m_{a}^{3}m_{b}^{2}-3080T_{a}T_{b}^{4}m_{a}^{4}m_{b}-560T_{b}^{5}m_{a}^{5}\big{\\}}$
$\displaystyle\times\big{[}10(T_{a}m_{b}+T_{b}m_{a})^{4}(m_{a}+m_{b})T_{a}\big{]}^{-1};$
$\displaystyle\hat{F}_{ab(2)}$ $\displaystyle=$ $\displaystyle
3T_{a}m_{b}^{2}\big{\\{}70T_{a}^{3}m_{a}m_{b}^{2}+102T_{a}^{3}m_{b}^{3}+385T_{a}^{2}T_{b}m_{a}^{2}m_{b}+561T_{a}^{2}T_{b}m_{a}m_{b}^{2}+1890T_{a}T_{b}^{2}m_{a}^{3}$
$\displaystyle+1446T_{a}T_{b}^{2}m_{a}^{2}m_{b}-588T_{b}^{3}m_{a}^{3}\big{\\}}\big{[}10(T_{a}m_{b}+T_{b}m_{a})^{4}(m_{a}+m_{b})\big{]}^{-1};$
$\displaystyle\hat{G}_{ab(1)}$ $\displaystyle=$
$\displaystyle-\big{\\{}565T_{a}^{4}m_{a}m_{b}^{4}+533T_{a}^{4}m_{b}^{5}+1270T_{a}^{3}T_{b}m_{a}^{2}m_{b}^{3}+1190T_{a}^{3}T_{b}m_{a}m_{b}^{4}+1020T_{a}^{2}T_{b}^{2}m_{a}^{3}m_{b}^{2}$
$\displaystyle+1152T_{a}^{2}T_{b}^{2}m_{a}^{2}m_{b}^{3}+3640T_{a}T_{b}^{3}m_{a}^{4}m_{b}+1916T_{a}T_{b}^{3}m_{a}^{3}m_{b}^{2}-1400T_{b}^{4}m_{a}^{5}-3304T_{b}^{4}m_{a}^{4}m_{b}\big{\\}}$
$\displaystyle\times\big{[}280(T_{a}m_{b}+T_{b}m_{a})^{4}(m_{a}+m_{b})\big{]}^{-1};$
$\displaystyle\hat{G}_{ab(2)}$ $\displaystyle=$
$\displaystyle-\,\frac{3T_{a}T_{b}m_{a}m_{b}^{2}(3T_{a}^{2}m_{a}m_{b}-5T_{a}^{2}m_{b}^{2}-42T_{a}T_{b}m_{a}^{2}-38T_{a}T_{b}m_{a}m_{b}+12T_{b}^{2}m_{a}^{2})}{8(T_{a}m_{b}+T_{b}m_{a})^{4}(m_{a}+m_{b})}.$
(21)

Exchange rates for the usual stress-tensor are given by

$\displaystyle\bar{\bar{\boldsymbol{Q}}}_{a}^{(2)}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\,\frac{21}{10}\nu_{aa}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}+\frac{9}{70}\nu_{aa}\frac{\rho_{a}}{p_{a}}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(4)}$
(22) $\displaystyle+\sum_{b\neq
a}\frac{\rho_{a}\nu_{ab}}{m_{a}+m_{b}}\Big{[}-\hat{K}_{ab(1)}\frac{1}{n_{a}}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}+\hat{K}_{ab(2)}\frac{1}{n_{b}}\bar{\bar{\boldsymbol{\Pi}}}_{b}^{(2)}$
$\displaystyle+L_{ab(1)}\frac{\rho_{a}}{n_{a}p_{a}}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(4)}-L_{ab(2)}\frac{\rho_{b}}{n_{b}p_{b}}\bar{\bar{\boldsymbol{\Pi}}}_{b}^{(4)}\Big{]},$

with mass-ratio coefficients

$\displaystyle\hat{K}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{10T_{a}^{2}m_{a}m_{b}^{2}+15T_{a}^{2}m_{b}^{3}+35T_{a}T_{b}m_{a}^{2}m_{b}+42T_{a}T_{b}m_{a}m_{b}^{2}+10T_{b}^{2}m_{a}^{3}+12T_{b}^{2}m_{a}^{2}m_{b}}{5(T_{a}m_{b}+T_{b}m_{a})^{2}m_{a}};$
$\displaystyle\hat{K}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{6T_{a}^{2}m_{a}m_{b}+4T_{a}^{2}m_{b}^{2}+21T_{a}T_{b}m_{a}^{2}+14T_{a}T_{b}m_{a}m_{b}-5T_{b}^{2}m_{a}^{2}}{5(T_{a}m_{b}+T_{b}m_{a})^{2}};$
$\displaystyle L_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{3T_{a}m_{b}(2T_{a}m_{a}m_{b}+3T_{a}m_{b}^{2}+7T_{b}m_{a}^{2}+8T_{b}m_{a}m_{b})}{35(T_{a}m_{b}+T_{b}m_{a})^{2}m_{a}};$
$\displaystyle L_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{3m_{a}T_{b}(5T_{a}m_{a}+4T_{a}m_{b}-T_{b}m_{a})}{35(T_{a}m_{b}+T_{b}m_{a})^{2}}.$
(23)

Finally, the 4th-order stress-tensor exchange rates are given by

$\displaystyle\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\,\frac{53}{20}\nu_{aa}\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}-\frac{79}{140}\nu_{aa}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\sum_{b\neq
a}\nu_{ab}\Big{[}-\hat{M}_{ab(1)}\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$
(24)
$\displaystyle+\hat{M}_{ab(2)}\frac{p_{a}^{2}}{\rho_{a}p_{b}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{b}-N_{ab(1)}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-N_{ab(2)}\frac{p_{a}^{2}\rho_{b}}{p_{b}^{2}\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{b}\Big{]},$

with mass-ratio coefficients

$\displaystyle\hat{M}_{ab(1)}$ $\displaystyle=$
$\displaystyle\Big{\\{}48T_{a}^{4}m_{a}m_{b}^{3}+36T_{a}^{4}m_{b}^{4}+216T_{a}^{3}T_{b}m_{a}^{2}m_{b}^{2}+107T_{a}^{3}T_{b}m_{a}m_{b}^{3}+378T_{a}^{2}T_{b}^{2}m_{a}^{3}m_{b}$
$\displaystyle+36T_{a}^{2}T_{b}^{2}m_{a}^{2}m_{b}^{2}-315T_{a}T_{b}^{3}m_{a}^{3}m_{b}-70T_{b}^{4}m_{a}^{4}\Big{\\}}\Big{[}5(T_{a}m_{b}+T_{b}m_{a})^{3}T_{a}(m_{b}+m_{a})\Big{]}^{-1};$
$\displaystyle\hat{M}_{ab(2)}$ $\displaystyle=$
$\displaystyle-\,\Big{\\{}T_{b}m_{a}\big{(}18T_{a}^{3}m_{a}m_{b}^{2}-4T_{a}^{3}m_{b}^{3}+81T_{a}^{2}T_{b}m_{a}^{2}m_{b}-18T_{a}^{2}T_{b}m_{a}m_{b}^{2}-147T_{a}T_{b}^{2}m_{a}^{3}$
$\displaystyle-189T_{a}T_{b}^{2}m_{a}^{2}m_{b}+35T_{b}^{3}m_{a}^{3}\big{)}\Big{\\}}\Big{[}5(T_{a}m_{b}+T_{b}m_{a})^{3}T_{a}(m_{b}+m_{a})\Big{]}^{-1};$
$\displaystyle N_{ab(1)}$ $\displaystyle=$
$\displaystyle-\,\Big{\\{}16T_{a}^{3}m_{a}m_{b}^{3}+12T_{a}^{3}m_{b}^{4}+72T_{a}^{2}T_{b}m_{a}^{2}m_{b}^{2}+21T_{a}^{2}T_{b}m_{a}m_{b}^{3}+126T_{a}T_{b}^{2}m_{a}^{3}m_{b}$
$\displaystyle\quad-54T_{a}T_{b}^{2}m_{a}^{2}m_{b}^{2}-140T_{b}^{3}m_{a}^{4}-273T_{b}^{3}m_{a}^{3}m_{b}\Big{\\}}\Big{[}35(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{b}+m_{a})\Big{]}^{-1};$
$\displaystyle N_{ab(2)}$ $\displaystyle=$
$\displaystyle-\,\frac{3T_{b}^{2}m_{a}^{2}(35T_{a}^{2}m_{a}m_{b}+12T_{a}^{2}m_{b}^{2}-35T_{a}T_{b}m_{a}^{2}-51T_{a}T_{b}m_{a}m_{b}+7T_{b}^{2}m_{a}^{2})}{35(T_{a}m_{b}+T_{b}m_{a})^{3}T_{a}(m_{b}+m_{a})}.$
(25)

The entire system is now fully specified, and represents a multi-fluid
generalization of the model of Braginskii (1965). Coupled with Maxwell’s
equations, it can be used in multi-fluid numerical simulations. Importantly,
when collisional frequencies become small, the right hand sides of evolution
equations just become small and no coefficients become divergent, which is in
contrast to the model of Braginskii, where the quasi-static approximation is
used for the stress-tensors and heat fluxes. For a detailed discussion on the
limitations of our model in a regime of low-collisionality, see Section 8.9.
The model of Braginskii is obtained as a particular case of a one ion-electron
plasma with similar temperatures, in a quasi-static and quasi-linear
approximation for the viscosity tensors and heat fluxes, where additionally,
the coupling between viscosity tensors and heat fluxes is neglected.

###  2.2 Collisional contributions for small temperature differences

In many instances, it might be satisfactory to consider a situation when the
temperature differences between species are small. The mass-ratio coefficients
(16) then become

$\displaystyle V_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{(21/10)m_{b}+(3/5)m_{a}}{m_{b}+m_{a}};\qquad
V_{ab(2)}=\frac{(3/5)m_{b}+(21/10)m_{a}}{m_{b}+m_{a}},$ (26)

mass-ratio coefficients (19) simplify into

$\displaystyle\hat{D}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{3m_{a}^{3}+(86/5)m_{a}^{2}m_{b}+(77/10)m_{a}m_{b}^{2}+(19/4)m_{b}^{3}}{(m_{a}+m_{b})^{3}};$
$\displaystyle\hat{D}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{(279/20)m_{a}m_{b}^{2}+(27/10)m_{b}^{3}}{(m_{a}+m_{b})^{3}};$
$\displaystyle\hat{E}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{(9/20)m_{a}^{2}m_{b}+(6/35)m_{a}m_{b}^{2}+(69/560)m_{b}^{3}}{(m_{a}+m_{b})^{3}};$
$\displaystyle\hat{E}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{(45/112)m_{a}m_{b}^{2}}{(m_{a}+m_{b})^{3}};\qquad\hat{U}_{ab(1)}=\frac{3}{2}\frac{m_{b}}{(m_{a}+m_{b})},$
(27)

mass-ratio coefficients (21) become

$\displaystyle\hat{F}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{(-56)m_{a}^{4}+336m_{a}^{3}m_{b}+(1302/5)m_{a}^{2}m_{b}^{2}+(1034/5)m_{a}m_{b}^{3}+(759/10)m_{b}^{4}}{(m_{a}+m_{b})^{4}};$
$\displaystyle\hat{F}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{(1953/5)m_{a}^{2}m_{b}^{2}+(1587/10)m_{a}m_{b}^{3}+(153/5)m_{b}^{4}}{(m_{a}+m_{b})^{4}};$
$\displaystyle\hat{G}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{5m_{a}^{4}-(31/5)m_{a}^{3}m_{b}-(30/7)m_{a}^{2}m_{b}^{2}-(611/140)m_{a}m_{b}^{3}-(533/280)m_{b}^{4}}{(m_{a}+m_{b})^{4}};$
$\displaystyle\hat{G}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{(45/4)m_{a}^{2}m_{b}^{2}+(15/8)m_{a}m_{b}^{3}}{(m_{a}+m_{b})^{4}};\qquad\hat{U}_{ab(2)}=\frac{42m_{a}m_{b}+27m_{b}^{2}}{(m_{a}+m_{b})^{2}},$
(28)

mass-ratio coefficients (23) become

$\displaystyle\hat{K}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{10m_{a}^{2}+37m_{a}m_{b}+15m_{b}^{2}}{5m_{a}(m_{b}+m_{a})};\qquad\hat{K}_{ab(2)}=\frac{4(4m_{a}+m_{b})}{5(m_{b}+m_{a})};$
$\displaystyle L_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{3(7m_{a}+3m_{b})m_{b}}{35m_{a}(m_{b}+m_{a})};\qquad
L_{ab(2)}=\frac{12m_{a}}{35(m_{a}+m_{b})},$ (29)

and mass-ratio coefficients (25) simplify into

$\displaystyle\hat{M}_{ab(1)}$ $\displaystyle=$
$\displaystyle-\,\frac{70m_{a}^{3}-133m_{a}^{2}m_{b}-119m_{a}m_{b}^{2}-36m_{b}^{3}}{5(m_{b}+m_{a})^{3}};\qquad\hat{M}_{ab(2)}=\frac{4m_{a}(28m_{a}^{2}-m_{a}m_{b}+m_{b}^{2})}{5(m_{b}+m_{a})^{3}};$
$\displaystyle N_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{140m_{a}^{3}+7m_{a}^{2}m_{b}-25m_{a}m_{b}^{2}-12m_{b}^{3}}{35(m_{b}+m_{a})^{3}};\qquad
N_{ab(2)}=\frac{12m_{a}^{2}(7m_{a}-3m_{b})}{35(m_{b}+m_{a})^{3}}.$ (30)

###  2.3 Reduction to 13-moment model

As a partial double-check of our calculations, neglecting the evolution
equations (11)-(12) for $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ &
$\vec{\boldsymbol{X}}^{(5)}_{a}$, and in the evolution equations (9)-(10) for
$\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ & $\vec{\boldsymbol{q}}_{a}$
prescribing closures (which are derived from Hermite closures)

$\vec{\boldsymbol{X}}^{(5)}_{a}=28\frac{p_{a}}{\rho_{a}}\vec{\boldsymbol{q}}_{a};\qquad\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}=7\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a},$
(31)

our 21-moment model simplifies into the 13-moment model, given by collisional
contributions

$\displaystyle\boldsymbol{R}_{a}$ $\displaystyle=$ $\displaystyle\sum_{b\neq
a}\nu_{ab}\Big{[}\rho_{a}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})+\frac{3}{5}\frac{\mu_{ab}}{T_{ab}}\big{(}\vec{\boldsymbol{q}}_{a}-\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{q}}_{b}\big{)}\Big{]};$
$\displaystyle\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{4}{5}\nu_{aa}\vec{\boldsymbol{q}}_{a}+\sum_{b\neq
a}\nu_{ab}\Big{[}-\hat{D}^{*}_{ab(1)}\vec{\boldsymbol{q}}_{a}+\hat{D}^{*}_{ab(2)}\frac{\rho_{a}}{\rho_{b}}\vec{\boldsymbol{q}}_{b}-p_{a}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\hat{U}_{ab(1)}\Big{]};$
$\displaystyle\bar{\bar{\boldsymbol{Q}}}_{a}^{(2)}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\,\frac{6}{5}\nu_{aa}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}+\sum_{b\neq
a}\frac{\rho_{a}\nu_{ab}}{m_{a}+m_{b}}\Big{[}-\hat{K}^{*}_{ab(1)}\frac{1}{n_{a}}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}+\hat{K}^{*}_{ab(2)}\frac{1}{n_{b}}\bar{\bar{\boldsymbol{\Pi}}}_{b}^{(2)}\Big{]},$
(32)

with mass-ratio coefficients

$\displaystyle\hat{D}^{*}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{9T_{a}^{2}m_{a}m_{b}^{2}+13T_{a}^{2}m_{b}^{3}-6T_{a}T_{b}m_{a}^{2}m_{b}+20T_{a}T_{b}m_{a}m_{b}^{2}+30T_{b}^{2}m_{a}^{3}+52T_{b}^{2}m_{a}^{2}m_{b}}{10(m_{a}+m_{b})(T_{a}m_{b}+T_{b}m_{a})^{2}};$
$\displaystyle\hat{D}^{*}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{9T_{a}m_{b}^{2}(5T_{a}m_{a}+3T_{a}m_{b}-2T_{b}m_{a})}{10(m_{a}+m_{b})(T_{a}m_{b}+T_{b}m_{a})^{2}};$
$\displaystyle\hat{K}^{*}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{2(2T_{a}m_{a}m_{b}+3T_{a}m_{b}^{2}+5T_{b}m_{a}^{2}+6T_{b}m_{a}m_{b})}{5m_{a}(T_{a}m_{b}+T_{b}m_{a})};\qquad\hat{K}^{*}_{ab(2)}=\frac{2(3T_{a}m_{a}+2T_{a}m_{b}-T_{b}m_{a})}{5(T_{a}m_{b}+T_{b}m_{a})},$
(33)

where $\hat{U}_{ab(1)}$ is unchanged from the 21-moment model. It can be shown
that for Coulomb collisions, this model is equivalent to equations (44)-(49)
of Schunk (1977), first calculated by Burgers (1969). For small temperature
differences the mass-ratio coefficients become

$\displaystyle\hat{D}^{*}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{30m_{a}^{2}+16m_{a}m_{b}+13m_{b}^{2}}{10(m_{a}+m_{b})^{2}};\qquad\hat{D}^{*}_{ab(2)}=\frac{27m_{b}^{2}}{10(m_{a}+m_{b})^{2}};$
$\displaystyle\hat{K}^{*}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{2m_{a}+(6/5)m_{b}}{m_{a}};\qquad\hat{K}^{*}_{ab(2)}=\frac{4}{5};\qquad\hat{U}_{ab(1)}=\frac{3}{2}\frac{m_{b}}{(m_{a}+m_{b})}.$
(34)

Our new 21-moment model thus can be viewed as a generalization of the multi-
fluid description of Burgers (1969) & Schunk (1977), where the heat fluxes and
stress-tensors are described more accurately, and with the same level of
precision as in Braginskii (1965). Nevertheless, we only use the Landau
collisional operator applicable for Coulomb collisions, whereas Burgers-Schunk
use the more general Boltzmann collisional operator and account for several
different interaction potentials.

###  2.4 Semi-linear approximation (de-coupled stress tensors and heat fluxes)

Here we consider the 21-moment model with evolution equations (9)-(12) in the
semi-linear approximation, where additionally viscosity-tensors and heat
fluxes are de-coupled. It will be shown later that the contributions
introduced by the coupling are smaller by a factor of $1/\nu_{aa}$. Within the
semi-linear approximation we also assume that there are no large-scale
gradients of considered fluid moments. For example, the de-coupling removes
the last terms at the left hand side of equations (10), (11), (12)
proportional to $(\nabla p_{a})\vec{\boldsymbol{q}}_{a}$, $(\nabla
p_{a})\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ and $(\nabla
p_{a})\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$. We neglect these terms
within the semi-linear approximation also when the coupling is considered (see
Sections 5 & 6). In the presence of large-scale gradients in
pressure/temperature these terms might become significant, together with many
other terms that are neglected in the semi-linear approximation. Evolution
equations for heat flux vectors simplify into

$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{q}}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{a}+\frac{5}{2}p_{a}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime};$
(35)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{a}+70\frac{p_{a}^{2}}{\rho_{a}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime},$
(36)

and evolution equations for viscosity-tensors become

$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}^{S}+p_{a}\bar{\bar{\boldsymbol{W}}}_{a}=\bar{\bar{\boldsymbol{Q}}}_{a}^{(2)}\,{}^{\prime};$
(37)
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+7\frac{p_{a}^{2}}{\rho_{a}}\bar{\bar{\boldsymbol{W}}}_{a}=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}\,{}^{\prime}.$
(38)

The above system will be used to recover the transport coefficients of
Braginskii (1965). In some instances, it might be actually advantageous to
suppress the non-linearities in numerical simulations, and perform multi-fluid
simulations with system (35)-(38) instead of the system (9)-(12).

##  3 One ion-electron plasma

###  3.1 Ion heat flux $\vec{\boldsymbol{q}}_{a}$ of Braginskii (self-
collisions)

Here we consider a one ion-electron plasma of similar temperatures, which is
the choice of Braginskii (1965). For the ion heat flux, Braginskii neglects
ion-electron collisions. Considering only self-collisions, evolution equations
for ion heat fluxes read

$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{q}}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{a}+\frac{5}{2}p_{a}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=-\frac{4}{5}\nu_{aa}\vec{\boldsymbol{q}}_{a}+\frac{3}{70}\nu_{aa}\Big{(}\frac{\rho_{a}}{p_{a}}\vec{\boldsymbol{X}}^{(5)}_{a}-28\vec{\boldsymbol{q}}_{a}\Big{)};$
(39)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{a}+70\frac{p_{a}^{2}}{\rho_{a}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=-\frac{88}{5}\nu_{aa}\frac{p_{a}}{\rho_{a}}\vec{\boldsymbol{q}}_{a}-\frac{3}{35}\nu_{aa}\Big{(}\vec{\boldsymbol{X}}^{(5)}_{a}-28\frac{p_{a}}{\rho_{a}}\vec{\boldsymbol{q}}_{a}\Big{)}.$
(40)

Neglecting the evolution equation (40) and prescribing closure (31) which
neglects the second term on the r.h.s. of (39), yields the ion heat flux model
of Burgers-Schunk, with the well-known $-4/5$ constant. However, now the
equations read

$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{q}}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{a}+\frac{5}{2}p_{a}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=-2\nu_{aa}\vec{\boldsymbol{q}}_{a}+\frac{3}{70}\nu_{aa}\frac{\rho_{a}}{p_{a}}\vec{\boldsymbol{X}}^{(5)}_{a};$
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{a}+70\frac{p_{a}^{2}}{\rho_{a}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=-\frac{76}{5}\nu_{aa}\frac{p_{a}}{\rho_{a}}\vec{\boldsymbol{q}}_{a}-\frac{3}{35}\nu_{aa}\vec{\boldsymbol{X}}^{(5)}_{a}.$
(41)

Prescribing the quasi-static approximation (by canceling the $d_{a}/dt$),
yields an analytic solution (see for example a general vector equation (E23)
with solution (E24))

$\vec{\boldsymbol{q}}_{a}=-\kappa_{\parallel}^{a}\nabla_{\parallel}T_{a}-\kappa_{\perp}^{a}\nabla_{\perp}T_{a}+\kappa_{\times}^{a}\boldsymbol{\hat{b}}\times\nabla
T_{a},$ (42)

and thermal conductivities

$\displaystyle\kappa_{\parallel}^{a}$ $\displaystyle=$
$\displaystyle\frac{125}{32}\frac{p_{a}}{\nu_{aa}m_{a}};$
$\displaystyle\kappa_{\perp}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\,\frac{2x^{2}+(648/245)}{x^{4}+(3313/1225)x^{2}+(20736/30625)};$
$\displaystyle\kappa_{\times}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\,\frac{(5/2)x^{3}+(2277/490)x}{x^{4}+(3313/1225)x^{2}+(20736/30625)},$
(43)

where $x=\Omega_{a}/\nu_{aa}$. Alternatively, by using numerical values

$\displaystyle\kappa_{\parallel}^{a}$ $\displaystyle=$ $\displaystyle
3.906\frac{p_{a}}{\nu_{aa}m_{a}};$ $\displaystyle\kappa_{\perp}^{a}$
$\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\,\frac{2x^{2}+2.645}{x^{4}+2.704x^{2}+0.6771};$
$\displaystyle\kappa_{\times}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\,\frac{(5/2)x^{3}+4.647x}{x^{4}+2.704x^{2}+0.6771},$
(44)

which recovers the ion heat flux of Braginskii (1965), his equation (4.40). We
use Braginskii notation with vectors
$\nabla_{\parallel}=\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\cdot\nabla$ and
$\nabla_{\perp}=\bar{\bar{\boldsymbol{I}}}_{\perp}\cdot\nabla=-\boldsymbol{\hat{b}}\times\boldsymbol{\hat{b}}\times\nabla$.

###  3.2 Ion heat flux $\vec{\boldsymbol{X}}^{(5)}_{a}$ (self-collisions)

The solution for the vector $\vec{\boldsymbol{X}}^{(5)}_{a}$ has a similar
form

$\vec{\boldsymbol{X}}^{(5)}_{a}=\frac{p_{a}}{\rho_{a}}\Big{[}-\kappa_{\parallel}^{a(5)}\nabla_{\parallel}T_{a}-\kappa_{\perp}^{a(5)}\nabla_{\perp}T_{a}+\kappa_{\times}^{a(5)}\boldsymbol{\hat{b}}\times\nabla
T_{a}\Big{]},$ (45)

with “thermal conductivities”

$\displaystyle\kappa_{\parallel}^{a(5)}$ $\displaystyle=$
$\displaystyle\underbrace{\frac{2975}{24}}_{123.96}\frac{p_{a}}{\nu_{aa}m_{a}};$
$\displaystyle\kappa_{\perp}^{a(5)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\frac{44x^{2}+(14688/175)}{x^{4}+(3313/1225)x^{2}+(20736/30625)};$
$\displaystyle\kappa_{\times}^{a(5)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\frac{70x^{3}+(1086/7)x}{x^{4}+(3313/1225)x^{2}+(20736/30625)}.$
(46)

###  3.3 Electron heat flux $\vec{\boldsymbol{q}}_{e}$ of Braginskii

Considering a one-ion electron plasma with similar temperatures, and keeping
only the dominant term in a $m_{e}/m_{i}$ expansion, the mass-ratio
coefficients (26), (27), (28) simplify into

$\displaystyle V_{ei(1)}=\frac{21}{10};\quad V_{ei(2)}=\frac{3}{5};$
$\displaystyle\hat{D}_{ei(1)}=\frac{19}{4};\quad\hat{D}_{ei(2)}=\frac{27}{10};\quad\hat{E}_{ei(1)}=\frac{69}{560};\quad\hat{E}_{ei(2)}=\frac{45}{112}\frac{m_{e}}{m_{i}};\quad\hat{U}_{ei(1)}=\frac{3}{2}$
$\displaystyle\hat{F}_{ei(1)}=\frac{759}{10};\quad\hat{F}_{ei(2)}=\frac{153}{5};\quad\hat{G}_{ei(1)}=-\,\frac{533}{280};\quad\hat{G}_{ei(2)}=\frac{15}{8}\frac{m_{e}}{m_{i}};\quad\hat{U}_{ei(2)}=27,$
(47)

collisional exchange rates become

$\displaystyle\boldsymbol{R}_{e}$ $\displaystyle=$
$\displaystyle-\rho_{e}\nu_{ei}\delta\boldsymbol{u}+\frac{21}{10}\frac{\rho_{e}}{p_{e}}\nu_{ei}\vec{\boldsymbol{q}}_{e}-\frac{3}{56}\frac{\rho_{e}^{2}}{p_{e}^{2}}\nu_{ei}\vec{\boldsymbol{X}}^{(5)}_{e};$
(48) $\displaystyle\vec{\boldsymbol{Q}}^{(3)}_{e}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle+\frac{3}{2}p_{e}\nu_{ei}\delta\boldsymbol{u}-\Big{[}2\nu_{ee}+\frac{19}{4}\nu_{ei}\Big{]}\vec{\boldsymbol{q}}_{e}+\Big{[}\frac{3}{70}\nu_{ee}+\frac{69}{560}\nu_{ei}\Big{]}\frac{\rho_{e}}{p_{e}}\vec{\boldsymbol{X}}^{(5)}_{e};$
(49) $\displaystyle\vec{\boldsymbol{Q}}^{(5)}_{e}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle+27\frac{p_{e}^{2}}{\rho_{e}}\nu_{ei}\delta\boldsymbol{u}-\Big{[}\frac{76}{5}\nu_{ee}+\frac{759}{10}\nu_{ei}\Big{]}\frac{p_{e}}{\rho_{e}}\vec{\boldsymbol{q}}_{e}-\Big{[}\frac{3}{35}\nu_{ee}-\frac{533}{280}\nu_{ei}\Big{]}\vec{\boldsymbol{X}}^{(5)}_{e},$
(50)

where $\delta\boldsymbol{u}=\boldsymbol{u}_{e}-\boldsymbol{u}_{i}$, and enter
the right hand side of the electron momentum equation, and evolution equations
for the electron heat flux vectors

$\displaystyle\frac{d_{e}}{dt}\vec{\boldsymbol{q}}_{e}+\Omega_{e}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{e}+\frac{5}{2}p_{e}\nabla\Big{(}\frac{p_{e}}{\rho_{e}}\Big{)}=\vec{\boldsymbol{Q}}^{(3)}_{e}\,{}^{\prime};$
$\displaystyle\frac{d_{e}}{dt}\vec{\boldsymbol{X}}^{(5)}_{e}+\Omega_{e}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{e}+70\frac{p_{e}^{2}}{\rho_{e}}\nabla\Big{(}\frac{p_{e}}{\rho_{e}}\Big{)}=\vec{\boldsymbol{Q}}^{(5)}_{e}\,{}^{\prime}.$
(51)

In Braginskii (1965), the results are expressed through the collisional
frequency $\nu_{ei}$, and conversion with $\nu_{ee}=\nu_{ei}/(Z_{i}\sqrt{2})$
yields

$\displaystyle\vec{\boldsymbol{Q}}^{(3)}_{e}\,{}^{\prime}$ $\displaystyle=$
$\displaystyle+\frac{3}{2}p_{e}\nu_{ei}\delta\boldsymbol{u}-\Big{[}\frac{\sqrt{2}}{Z_{i}}+\frac{19}{4}\Big{]}\nu_{ei}\vec{\boldsymbol{q}}_{e}+\Big{[}\frac{3}{70\sqrt{2}Z_{i}}+\frac{69}{560}\Big{]}\nu_{ei}\frac{\rho_{e}}{p_{e}}\vec{\boldsymbol{X}}^{(5)}_{e};$
$\displaystyle\vec{\boldsymbol{Q}}^{(5)}_{e}\,{}^{\prime}$ $\displaystyle=$
$\displaystyle+27\frac{p_{e}^{2}}{\rho_{e}}\nu_{ei}\delta\boldsymbol{u}-\Big{[}\frac{76}{5\sqrt{2}Z_{i}}+\frac{759}{10}\Big{]}\nu_{ei}\frac{p_{e}}{\rho_{e}}\vec{\boldsymbol{q}}_{e}-\Big{[}\frac{3}{35\sqrt{2}Z_{i}}-\frac{533}{280}\Big{]}\nu_{ei}\vec{\boldsymbol{X}}^{(5)}_{e}.$
(52)

In a quasi-static approximation, the solution of (51), (52) recovers the
famous electron heat flux of Braginskii (1965), together with vector
$\vec{\boldsymbol{X}}^{(5)}_{e}$ (which is of course not given by Braginskii).
Substituting these results into the momentum exchange rates (48), recovers the
$\boldsymbol{R}_{e}$ of Braginskii.

We use the same notation as Braginskii (1965) with $x=\Omega_{e}/\nu_{ei}$,
except (as is the norm in more recent papers) our $\Omega_{e}$ is formulated
as a general $\Omega_{a}$ and is thus negative, whereas in Braginskii
$\Omega_{e}$ is defined as positive. This yields a simple change of signs in
front of the “cross” ($\times$) terms with respect to Braginskii. In a quasi-
static approximation, the electron heat flux is split into a thermal and
frictional part
$\vec{\boldsymbol{q}}_{e}=\vec{\boldsymbol{q}}_{e}^{T}+\vec{\boldsymbol{q}}_{e}^{u}$,
where

$\displaystyle\vec{\boldsymbol{q}}_{e}^{T}$ $\displaystyle=$
$\displaystyle-\kappa_{\parallel}^{e}\nabla_{\parallel}T_{e}-\kappa_{\perp}^{e}\nabla_{\perp}T_{e}+\kappa_{\times}^{e}\boldsymbol{\hat{b}}\times\nabla
T_{e};$ $\displaystyle\vec{\boldsymbol{q}}_{e}^{u}$ $\displaystyle=$
$\displaystyle\beta_{0}p_{e}\delta\boldsymbol{u}_{\parallel}+p_{e}\delta\boldsymbol{u}_{\perp}\frac{\beta_{1}^{\prime}x^{2}+\beta_{0}^{\prime}}{\triangle}-p_{e}\boldsymbol{\hat{b}}\times\delta\boldsymbol{u}\frac{\beta_{1}^{\prime\prime}x^{3}+\beta_{0}^{\prime\prime}x}{\triangle},$
(53)

and the heat conductivities are given by

$\kappa_{\parallel}^{e}=\frac{p_{e}}{m_{e}\nu_{ei}}\gamma_{0};\qquad\kappa_{\perp}^{e}=\frac{p_{e}}{m_{e}\nu_{ei}}\frac{\gamma_{1}^{\prime}x^{2}+\gamma_{0}^{\prime}}{\triangle};\qquad\kappa_{\times}^{e}=\frac{p_{e}}{m_{e}\nu_{ei}}\frac{\gamma_{1}^{\prime\prime}x^{3}+\gamma_{0}^{\prime\prime}x}{\triangle}.$
(54)

The momentum exchange rates are also split into a thermal and frictional part
$\boldsymbol{R}_{e}=\boldsymbol{R}_{e}^{T}+\boldsymbol{R}_{e}^{u}$ (thermal
force and friction force), according to

$\displaystyle\boldsymbol{R}_{e}^{u}$ $\displaystyle=$
$\displaystyle-\alpha_{0}\rho_{e}\nu_{ei}\delta\boldsymbol{u}_{\parallel}-\rho_{e}\nu_{ei}\delta\boldsymbol{u}_{\perp}\Big{(}1-\frac{\alpha_{1}^{\prime}x^{2}+\alpha_{0}^{\prime}}{\triangle}\Big{)}-\rho_{e}\nu_{ei}\boldsymbol{\hat{b}}\times\delta\boldsymbol{u}\frac{\alpha_{1}^{\prime\prime}x^{3}+\alpha_{0}^{\prime\prime}x}{\triangle};$
$\displaystyle\boldsymbol{R}_{e}^{T}$ $\displaystyle=$
$\displaystyle-\beta_{0}n_{e}\nabla_{\parallel}T_{e}-n_{e}\nabla_{\perp}T_{e}\frac{\beta_{1}^{\prime}x^{2}+\beta_{0}^{\prime}}{\triangle}+n_{e}\boldsymbol{\hat{b}}\times\nabla
T_{e}\frac{\beta_{1}^{\prime\prime}x^{3}+\beta_{0}^{\prime\prime}x}{\triangle}.$
(55)

Instead of a numerical Table II on page 25 of Braginskii (1965), we provide
all the coefficients in a fully analytic form for a general ion charge
$Z_{i}$, which are given by

$\displaystyle\alpha_{0}$ $\displaystyle=$
$\displaystyle\frac{4(16Z_{i}^{2}+61Z_{i}\sqrt{2}+72)}{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288};\qquad\beta_{0}=\frac{30Z_{i}(11Z_{i}+15\sqrt{2})}{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288};$
$\displaystyle\gamma_{0}$ $\displaystyle=$
$\displaystyle\frac{25Z_{i}(433Z_{i}+180\sqrt{2})}{4(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)};$
(56)

$\displaystyle\triangle=x^{4}+\delta_{1}x^{2}+\delta_{0};$
$\displaystyle\delta_{0}=\Big{(}\frac{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288}{700Z_{i}^{2}}\Big{)}^{2};$
$\displaystyle\delta_{1}=\frac{586601Z_{i}^{2}+330152Z_{i}\sqrt{2}+106016}{78400Z_{i}^{2}};$
(57) $\displaystyle\alpha_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{9(40337Z_{i}+10996\sqrt{2})}{78400Z_{i}};$
$\displaystyle\alpha_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{9(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(17Z_{i}+40\sqrt{2})}{490000Z_{i}^{3}};$
$\displaystyle\alpha_{1}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{477}{280};\qquad\alpha_{0}^{\prime\prime}=\frac{9(64Z_{i}^{2}+151Z_{i}\sqrt{2}+253)}{6125Z_{i}^{2}};$
(58) $\displaystyle\beta_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{3(709Z_{i}+172\sqrt{2})}{560Z_{i}};$
$\displaystyle\beta_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{3(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(11Z_{i}+15\sqrt{2})}{49000Z_{i}^{3}};$
$\displaystyle\beta_{1}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{3}{2};\qquad\beta_{0}^{{}^{\prime\prime}}=\frac{3(5729Z_{i}^{2}+6711Z_{i}\sqrt{2}+4728)}{19600Z_{i}^{2}};$
(59) $\displaystyle\gamma_{1}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{13Z_{i}+4\sqrt{2}}{4Z_{i}};$
$\displaystyle\gamma_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(433Z_{i}+180\sqrt{2})}{78400Z_{i}^{3}};$
$\displaystyle\gamma_{1}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{5}{2};\qquad\gamma_{0}^{\prime\prime}=\frac{320797Z_{i}^{2}+202248Z_{i}\sqrt{2}+72864}{31360Z_{i}^{2}}.$
(60)

Numerical values for $Z_{i}=1$ are given in the first column of Table II of
Braginskii (1965) and for example the parallel coefficients are
$\alpha_{0}=0.5129$; $\beta_{0}=0.7110$; $\gamma_{0}=3.1616$, matching his
values exactly. We checked the entire Table II of Braginskii and his table is
very precise, except for two values. For the $\alpha_{0}$ coefficient, values
for $Z_{i}=2,3$ should be changed as $0.4408\to 0.4309$; $0.3965\to 0.3954$.
The rest of his table is calculated very accurately, with around handful of
irrelevant last digit rounding changes (such as $3.7703\to 3.7702$ in
$\delta_{0}(Z_{i}=1)$, $0.2400\to 0.2399$ in
$\alpha_{0}^{\prime\prime}(Z_{i}=3)$; and for $Z_{i}=4$ charge $0.3752\to
0.3751$ in $\alpha_{0}$; $9.055\to 9.056$ in $\delta_{0}$ $0.4478\to 0.4477$
in $\beta_{0}^{\prime}$ etc.).

Analytic results (56) for parallel coefficients
$\alpha_{0},\beta_{0},\gamma_{0}$ were also obtained by Simakov & Molvig
(2014); see later Section 4. To triple-check our other results, we re-
calculated our approach with analytic collisional matrices of Ji & Held
(2013), equations (28a)-(28f), together with their equations (40)-(44) and
other formulas, which yielded the same analytic expressions. Unfortunately,
the analytic results of Balescu (1988) are written in a such a complicated
form, (see his page 236, with collisional matrices on page 198 and the
required conversion equation (5.7.13) on page 270), that we were able to
verify only an analytic match with his parallel coefficients. The formulation
of Balescu (1988) is so different from Braginskii, that Balescu by himself
(page 275) only claims a match of below 1% for the 21-moment model, not
further analyzing possible discrepancies.

###  3.4 Electron heat flux $\vec{\boldsymbol{X}}_{e}^{(5)}$

Similarly to the usual electron heat flux $\vec{\boldsymbol{q}}_{e}$, a quasi-
static solution for the heat flux vector $\vec{\boldsymbol{X}}_{e}^{(5)}$ has
to be split into a thermal and frictional part, according to

$\displaystyle\vec{\boldsymbol{X}}_{e}^{(5)T}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\rho_{e}}\Big{[}-\kappa_{\parallel}^{e(5)}\nabla_{\parallel}T_{e}-\kappa_{\perp}^{e(5)}\nabla_{\perp}T_{e}+\kappa_{\times}^{e(5)}\boldsymbol{\hat{b}}\times\nabla
T_{e}\Big{]};$ $\displaystyle\vec{\boldsymbol{X}}_{e}^{(5)u}$ $\displaystyle=$
$\displaystyle\frac{p_{e}^{2}}{\rho_{e}}\Big{[}\beta_{0}^{(5)}\delta\boldsymbol{u}_{\parallel}+\frac{\beta_{1}^{(5)^{\prime}}x^{2}+\beta_{0}^{(5)^{\prime}}}{\triangle}\delta\boldsymbol{u}_{\perp}-\frac{\beta_{1}^{(5)^{\prime\prime}}x^{3}+\beta_{0}^{(5)^{\prime\prime}}x}{\triangle}\boldsymbol{\hat{b}}\times\delta\boldsymbol{u}\Big{]},$
(61)

with thermal conductivities

$\kappa_{\parallel}^{e(5)}=\frac{p_{e}}{m_{e}\nu_{ei}}\gamma_{0}^{(5)};\qquad\kappa_{\perp}^{e(5)}=\frac{p_{e}}{m_{e}\nu_{ei}}\frac{\gamma_{1}^{(5)^{\prime}}x^{2}+\gamma_{0}^{(5)^{\prime}}}{\triangle};\qquad\kappa_{\times}^{e(5)}=\frac{p_{e}}{m_{e}\nu_{ei}}\frac{\gamma_{1}^{(5)^{\prime\prime}}x^{3}+\gamma_{0}^{(5)^{\prime\prime}}x}{\triangle}.$
(62)

The analytic coefficients are given by

$\displaystyle\beta_{0}^{(5)}$ $\displaystyle=$
$\displaystyle\frac{840Z_{i}(13\sqrt{2}+12Z_{i})}{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288};$
$\displaystyle\beta_{1}^{(5)^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{3(5829Z_{i}+1172\sqrt{2})}{280Z_{i}};$
$\displaystyle\beta_{0}^{(5)^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{3(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(12Z_{i}+13\sqrt{2})}{1750Z_{i}^{3}};$
$\displaystyle\beta_{1}^{(5)^{\prime\prime}}$ $\displaystyle=$ $\displaystyle
27;\qquad\beta_{0}^{(5)^{\prime\prime}}=\frac{3(7611Z_{i}^{2}+8429Z_{i}\sqrt{2}+5000)}{700Z_{i}^{2}},$
(63)

and

$\displaystyle\gamma_{0}^{(5)}$ $\displaystyle=$
$\displaystyle\frac{175Z_{i}(204\sqrt{2}+571Z_{i})}{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288};$
$\displaystyle\gamma_{1}^{(5)^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{113Z_{i}+44\sqrt{2}}{2Z_{i}};$
$\displaystyle\gamma_{0}^{(5)^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(571Z_{i}+204\sqrt{2})}{2800Z_{i}^{3}};$
$\displaystyle\gamma_{1}^{(5)^{\prime\prime}}$ $\displaystyle=$ $\displaystyle
70;\qquad\gamma_{0}^{(5)^{\prime\prime}}=\frac{430783Z_{i}^{2}+261672Z_{i}\sqrt{2}+86880}{1120Z_{i}^{2}},$
(64)

with $\triangle$ unchanged and given by (57). These results were substituted
into the momentum exchange rates $\boldsymbol{R}_{e}$, equation (48), to
obtain the final expression for the friction force and thermal force. Useful
relations are

$\displaystyle\alpha_{0}$ $\displaystyle=$ $\displaystyle
1-\frac{21}{10}\beta_{0}+\frac{3}{56}\beta_{0}^{(5)};\qquad\alpha_{1}^{\prime}=\frac{21}{10}\beta_{1}^{\prime}-\frac{3}{56}\beta_{1}^{(5)^{\prime}};$
$\displaystyle\alpha_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{21}{10}\beta_{0}^{\prime}-\frac{3}{56}\beta_{0}^{(5)^{\prime}};\qquad\alpha_{1}^{\prime\prime}=\frac{21}{10}\beta_{1}^{\prime\prime}-\frac{3}{56}\beta_{1}^{(5)^{\prime\prime}};\qquad\alpha_{0}^{\prime\prime}=\frac{21}{10}\beta_{0}^{\prime\prime}-\frac{3}{56}\beta_{0}^{(5)^{\prime\prime}};$
$\displaystyle\beta_{0}$ $\displaystyle=$
$\displaystyle\frac{21}{10}\gamma_{0}-\frac{3}{56}\gamma_{0}^{(5)};\qquad\beta_{1}^{\prime}=\frac{21}{10}\gamma_{1}^{\prime}-\frac{3}{56}\gamma_{1}^{(5)^{\prime}};$
$\displaystyle\beta_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{21}{10}\gamma_{0}^{\prime}-\frac{3}{56}\gamma_{0}^{(5)^{\prime}};\qquad\beta_{1}^{\prime\prime}=\frac{21}{10}\gamma_{1}^{\prime\prime}-\frac{3}{56}\gamma_{1}^{(5)^{\prime\prime}};\qquad\beta_{0}^{\prime\prime}=\frac{21}{10}\gamma_{0}^{\prime\prime}-\frac{3}{56}\gamma_{0}^{(5)^{\prime\prime}}.$
(65)

For $Z_{i}=1$, transport coefficients (63), (64) have numerical values

$\displaystyle\beta_{0}^{(5)}$ $\displaystyle=$ $\displaystyle
18.778;\qquad\beta_{1}^{(5)^{\prime}}=80.212;\qquad\beta_{0}^{(5)^{\prime}}=70.797;\qquad\beta_{1}^{(5)^{\prime\prime}}=27;\qquad\beta_{0}^{(5)^{\prime\prime}}=105.135;$
$\displaystyle\gamma_{0}^{(5)}$ $\displaystyle=$ $\displaystyle
110.664;\qquad\gamma_{1}^{(5)^{\prime}}=87.613;\qquad\gamma_{0}^{(5)^{\prime}}=417.221;\qquad\gamma_{1}^{(5)^{\prime\prime}}=70;\qquad\gamma_{0}^{(5)^{\prime\prime}}=792.610.$
(66)

###  3.5 Ion viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ of Braginskii
(self-collisions)

Considering self-collisions, evolution equations for the ion viscosity-tensors
read

$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}^{S}+p_{a}\bar{\bar{\boldsymbol{W}}}_{a}=-\,\frac{6}{5}\nu_{aa}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}+\frac{9}{70}\nu_{aa}\Big{(}\frac{\rho_{a}}{p_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-7\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\Big{)};$
(67)
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+7\frac{p_{a}^{2}}{\rho_{a}}\bar{\bar{\boldsymbol{W}}}_{a}=-\,\frac{33}{5}\nu_{aa}\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}-\frac{79}{140}\nu_{aa}\Big{(}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-7\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\Big{)}.$
(68)

Neglecting (68) and prescribing closure (31) which neglects the second term on
the r.h.s. of (67), yields the ion-viscosity model of Burgers-Schunk, with the
well-known $-6/5$ constant. However, now the equations read

$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}^{S}+p_{a}\bar{\bar{\boldsymbol{W}}}_{a}=-\,\frac{21}{10}\nu_{aa}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}+\frac{9}{70}\nu_{aa}\frac{\rho_{a}}{p_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a};$
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+7\frac{p_{a}^{2}}{\rho_{a}}\bar{\bar{\boldsymbol{W}}}_{a}=-\,\frac{53}{20}\nu_{aa}\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}-\frac{79}{140}\nu_{aa}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}.$
(69)

In a quasi-static approximation, solution of (69) yields
$\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ in the following form (see for
example Appendix E.4)

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ $\displaystyle=$
$\displaystyle-\eta_{0}^{a}\bar{\bar{\boldsymbol{W}}}_{0}-\eta_{1}^{a}\bar{\bar{\boldsymbol{W}}}_{1}-\eta_{2}^{a}\bar{\bar{\boldsymbol{W}}}_{2}+\eta_{3}^{a}\bar{\bar{\boldsymbol{W}}}_{3}+\eta_{4}^{a}\bar{\bar{\boldsymbol{W}}}_{4};$
(70) $\displaystyle\bar{\bar{\boldsymbol{W}}}_{0}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\big{(}\bar{\bar{\boldsymbol{W}}}_{a}:\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\big{)}\Big{(}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}\Big{)};$
$\displaystyle\bar{\bar{\boldsymbol{W}}}_{1}$ $\displaystyle=$
$\displaystyle\bar{\bar{\boldsymbol{I}}}_{\perp}\cdot\bar{\bar{\boldsymbol{W}}}_{a}\cdot\bar{\bar{\boldsymbol{I}}}_{\perp}+\frac{1}{2}\big{(}\bar{\bar{\boldsymbol{W}}}_{a}:\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\big{)}\bar{\bar{\boldsymbol{I}}}_{\perp};$
$\displaystyle\bar{\bar{\boldsymbol{W}}}_{2}$ $\displaystyle=$
$\displaystyle\big{(}\bar{\bar{\boldsymbol{I}}}_{\perp}\cdot\bar{\bar{\boldsymbol{W}}}_{a}\cdot\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\big{)}^{S};$
$\displaystyle\bar{\bar{\boldsymbol{W}}}_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{W}}}_{a}\cdot\bar{\bar{\boldsymbol{I}}}_{\perp}\big{)}^{S};$
$\displaystyle\bar{\bar{\boldsymbol{W}}}_{4}$ $\displaystyle=$
$\displaystyle\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{W}}}_{a}\cdot\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\big{)}^{S},$
(71)

which is equivalent to equations (4.41) & (4.42) of Braginskii (1965), after
one prescribes in his $\bar{\bar{\boldsymbol{W}}}_{0}$ that the matrix
$\bar{\bar{\boldsymbol{W}}}_{a}$ is traceless. Alternatively, with respect to
$\boldsymbol{\hat{b}}=(0,0,1)$ (straight magnetic field applied in the
z-direction)

$\displaystyle\Pi_{xx}^{a(2)}=-\frac{\eta_{0}^{a}}{2}(W_{xx}^{a}+W_{yy}^{a})-\frac{\eta_{1}^{a}}{2}(W_{xx}^{a}-W_{yy}^{a})-\eta_{3}^{a}W_{xy}^{a};$
$\displaystyle\Pi_{xy}^{a(2)}=\frac{\eta_{3}^{a}}{2}(W_{xx}^{a}-W_{yy}^{a})-\eta_{1}^{a}W_{xy}^{a};$
$\displaystyle\Pi_{xz}^{a(2)}=-\eta_{4}^{a}W_{yz}^{a}-\eta_{2}^{a}W_{xz}^{a};$
$\displaystyle\Pi_{yy}^{a(2)}=-\frac{\eta_{0}^{a}}{2}(W_{xx}^{a}+W_{yy}^{a})+\frac{\eta_{1}^{a}}{2}(W_{xx}^{a}-W_{yy}^{a})+\eta_{3}^{a}W_{xy}^{a};$
$\displaystyle\Pi_{yz}^{a(2)}=\eta_{4}^{a}W_{xz}^{a}-\eta_{2}^{a}W_{yz}^{a};$
$\displaystyle\Pi_{zz}^{a(2)}=-\eta_{0}^{a}W_{zz}^{a},$ (72)

which is equation (2.21) of Braginskii (1965). The ion viscosities are

$\displaystyle\eta_{0}^{a}$ $\displaystyle=$
$\displaystyle\frac{1025}{1068}\frac{p_{a}}{\nu_{aa}};$
$\displaystyle\eta_{2}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{(6/5)x^{2}+(10947/4900)}{x^{4}+(79321/19600)x^{2}+(71289/30625)};$
$\displaystyle\eta_{4}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{x^{3}+(46561/19600)x}{x^{4}+(79321/19600)x^{2}+(71289/30625)},$
(73)

where $x=\Omega_{a}/\nu_{aa}$, and $\eta_{1}^{a}(x)=\eta_{2}^{a}(2x)$;
$\eta_{3}^{a}(x)=\eta_{4}^{a}(2x)$ (The solution is easily obtained for the
parallel “zz” direction with $\Omega_{a}=0$, and for perpendicular directions
for example by choosing coupled “xz” and “yz” directions, and solving 4
equations in 4 unknowns). Alternatively, using numerical values

$\displaystyle\eta_{0}^{a}$ $\displaystyle=$ $\displaystyle
0.960\frac{p_{a}}{\nu_{aa}};$ $\displaystyle\eta_{2}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{(6/5)x^{2}+2.234}{x^{4}+4.047x^{2}+2.328};$
$\displaystyle\eta_{4}^{a}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{x^{3}+2.376x}{x^{4}+4.047x^{2}+2.328},$
(74)

recovering ion viscosities of Braginskii (1965), his equation (4.44).
Numerical values in Braginskii are evaluated precisely, with the sole
exception of one value in the denominator, where his rounded 4.03 should be
replaced by 4.05.

###  3.6 Ion viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ (self-
collisions)

The ion viscosity tensor $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ is given by

$\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}=\frac{p_{a}}{\rho_{a}}\Big{[}-\eta_{0}^{a(4)}\bar{\bar{\boldsymbol{W}}}_{0}-\eta_{1}^{a(4)}\bar{\bar{\boldsymbol{W}}}_{1}-\eta_{2}^{a(4)}\bar{\bar{\boldsymbol{W}}}_{2}+\eta_{3}^{a(4)}\bar{\bar{\boldsymbol{W}}}_{3}+\eta_{4}^{a(4)}\bar{\bar{\boldsymbol{W}}}_{4}\Big{]},$
(75)

with matrices $\bar{\bar{\boldsymbol{W}}}_{0}-\bar{\bar{\boldsymbol{W}}}_{4}$
(71) unchanged, and viscosities

$\displaystyle\eta_{0}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{8435}{1068}\frac{p_{a}}{\nu_{aa}};$
$\displaystyle\eta_{2}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{(33/5)x^{2}+(64347/3500)}{x^{4}+(79321/19600)x^{2}+(71289/30625)};$
$\displaystyle\eta_{4}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{7x^{3}+(59989/2800)x}{x^{4}+(79321/19600)x^{2}+(71289/30625)},$
(76)

where $\eta_{1}^{a(4)}(x)=\eta_{2}^{a(4)}(2x)$,
$\eta_{3}^{a(4)}(x)=\eta_{4}^{a(4)}(2x)$ holds, or with numerical values

$\displaystyle\eta_{0}^{a(4)}$ $\displaystyle=$ $\displaystyle
7.898\frac{p_{a}}{\nu_{aa}};$ $\displaystyle\eta_{2}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{6.600x^{2}+18.385}{x^{4}+4.047x^{2}+2.328};$
$\displaystyle\eta_{4}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}}\,\frac{7x^{3}+21.425x}{x^{4}+4.047x^{2}+2.328}.$
(77)

###  3.7 Electron viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}$ of
Braginskii

For a one ion-electron plasma with similar temperatures, the mass-ratio
coefficients (29), (30) simplify into

$\displaystyle\hat{K}_{ei(1)}$ $\displaystyle=$ $\displaystyle
3\frac{m_{i}}{m_{e}};\qquad\hat{K}_{ei(2)}=\frac{4}{5};\qquad
L_{ei(1)}=\frac{9}{35}\frac{m_{i}}{m_{e}};\qquad
L_{ei(2)}=\frac{12}{35}\frac{m_{e}}{m_{i}};$ $\displaystyle\hat{M}_{ei(1)}$
$\displaystyle=$
$\displaystyle\frac{36}{5};\qquad\hat{M}_{ei(2)}=\frac{4}{5}\frac{m_{e}}{m_{i}};\qquad
N_{ei(1)}=-\,\frac{12}{35};\qquad
N_{ei(2)}=-\,\frac{36}{35}\frac{m_{e}^{2}}{m_{i}^{2}},$ (78)

and collisional exchange rates for the viscosity-tensors become

$\displaystyle\bar{\bar{\boldsymbol{Q}}}_{e}^{(2)}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\Big{(}\frac{21}{10}\nu_{ee}+3\nu_{ei}\Big{)}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}+\Big{(}\frac{9}{70}\nu_{ee}+\frac{9}{35}\nu_{ei}\Big{)}\frac{\rho_{e}}{p_{e}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e};$
$\displaystyle\bar{\bar{\boldsymbol{Q}}}_{e}^{(4)}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\Big{(}\frac{53}{20}\nu_{ee}+\frac{36}{5}\nu_{ei}\Big{)}\frac{p_{e}}{\rho_{e}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}+\Big{(}-\frac{79}{140}\nu_{ee}+\frac{12}{35}\nu_{ei}\Big{)}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}.$
(79)

Converting everything to $\nu_{ei}$ with $\nu_{ee}=\nu_{ei}/(Z_{i}\sqrt{2})$
yields

$\displaystyle\bar{\bar{\boldsymbol{Q}}}_{e}^{(2)}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\Big{(}\frac{21}{10Z_{i}\sqrt{2}}+3\Big{)}\nu_{ei}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}+\Big{(}\frac{9}{70Z_{i}\sqrt{2}}+\frac{9}{35}\Big{)}\nu_{ei}\frac{\rho_{e}}{p_{e}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e};$
$\displaystyle\bar{\bar{\boldsymbol{Q}}}_{e}^{(4)}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle-\Big{(}\frac{53}{20Z_{i}\sqrt{2}}+\frac{36}{5}\Big{)}\nu_{ei}\frac{p_{e}}{\rho_{e}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}+\Big{(}-\frac{79}{140Z_{i}\sqrt{2}}+\frac{12}{35}\Big{)}\nu_{ei}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e},$
(80)

and these contributions enter the r.h.s. of evolution equations

$\displaystyle\frac{d_{e}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}+\Omega_{e}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}\big{)}^{S}+p_{e}\bar{\bar{\boldsymbol{W}}}_{e}=\bar{\bar{\boldsymbol{Q}}}_{e}^{(2)}\,{}^{\prime};$
$\displaystyle\frac{d_{e}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}+\Omega_{e}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}\big{)}^{S}+7\frac{p_{e}^{2}}{\rho_{e}}\bar{\bar{\boldsymbol{W}}}_{e}=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{e}\,{}^{\prime}.$
(81)

In a quasi-static approximation, solution of (80), (81) yields the electron
viscosity tensor $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}$ in form (70), (71),
with electron viscosities

$\displaystyle\eta_{0}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\,\frac{5Z_{i}(408Z_{i}+205\sqrt{2})}{6(192Z_{i}^{2}+301Z_{i}\sqrt{2}+178)};$
$\displaystyle\eta_{2}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\Big{[}\frac{3\sqrt{2}+6Z_{i}}{5Z_{i}}x^{2}+\frac{3(192Z_{i}^{2}+301Z_{i}\sqrt{2}+178)(408Z_{i}+205\sqrt{2})}{196000Z_{i}^{3}}\Big{]}/\triangle;$
$\displaystyle\eta_{4}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}x\Big{[}x^{2}+\frac{119520Z_{i}^{2}+101784\sqrt{2}Z_{i}+46561}{39200Z_{i}^{2}}\Big{]}/\triangle;$
$\displaystyle\triangle$ $\displaystyle=$ $\displaystyle
x^{4}+\frac{212256Z_{i}^{2}+176376\sqrt{2}Z_{i}+79321}{39200Z_{i}^{2}}x^{2}+\Big{(}\frac{3(192Z_{i}^{2}+301Z_{i}\sqrt{2}+178)}{700Z^{2}}\Big{)}^{2},$
(82)

where $x=\Omega_{e}/\nu_{ei}$, and relations
$\eta_{1}^{e}(x)=\eta_{2}^{e}(2x)$, $\eta_{3}^{e}(x)=\eta_{4}^{e}(2x)$. For
the particular case of $Z_{i}=1$ these electron viscosities become

$\displaystyle\eta_{0}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\,\frac{2040+1025\sqrt{2}}{2220+1806\sqrt{2}};$
$\displaystyle\eta_{2}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\Big{[}\frac{3\sqrt{2}+6}{5}x^{2}+\frac{297987}{98000}\sqrt{2}+\frac{82311}{19600}\Big{]}/\triangle;$
$\displaystyle\eta_{4}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}x\Big{[}x^{2}+\frac{12723}{4900}\sqrt{2}+\frac{166081}{39200}\Big{]}/\triangle;$
$\displaystyle\triangle$ $\displaystyle=$ $\displaystyle
x^{4}+\Big{(}\frac{22047}{4900}\sqrt{2}+\frac{291577}{39200}\Big{)}x^{2}+\Big{(}\frac{1431459}{245000}+\frac{14319}{3500}\sqrt{2}\Big{)},$
(83)

or with numerical values

$\displaystyle\eta_{0}^{e}$ $\displaystyle=$ $\displaystyle
0.73094\frac{p_{e}}{\nu_{ei}};$ $\displaystyle\eta_{2}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\big{(}2.049x^{2}+8.500\big{)}/\triangle;$
$\displaystyle\eta_{4}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}x\big{(}x^{2}+7.909\big{)}/\triangle;$
$\displaystyle\triangle$ $\displaystyle=$ $\displaystyle
x^{4}+13.801x^{2}+11.628,$ (84)

recovering the electron viscosity of Braginskii (1965), his equation (4.45).
It appears that the Braginskii parallel viscosity value of 0.733 is slightly
imprecise and should be 0.731 instead. The analytic result for parallel
viscosity $\eta^{e}_{0}$ agrees with Simakov & Molvig (2014), and the value
0.73094 agrees with Ji & Held (2013); see the inset of their Figure 3
(curiously, in a more precise 3-Laguerre approximation the coefficient changes
to 0.733). Note that for $x\to 0$ viscosity $\eta_{2}^{e}\to\eta_{0}^{e}$. As
discussed previously, our $\Omega_{e}$ is negative and in Braginskii it is
positive, yielding an opposite sign in front of $\eta_{4}^{e}$. Braginskii
offers electron viscosities only for $Z_{i}=1$. The analytic result (82) is
useful to quickly calculate electron viscosities for any $Z_{i}$. Ji & Held
(2013, 2015) also provide useful fitting formulas.

###  3.8 Electron viscosity $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}$

The solution for electron viscosity tensor
$\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}$ has form (75) with viscosities

$\displaystyle\eta_{0}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\,\frac{35Z_{i}(552Z_{i}+241\sqrt{2})}{6(192Z_{i}^{2}+301Z_{i}\sqrt{2}+178)};$
$\displaystyle\eta_{2}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\Big{[}\frac{33\sqrt{2}+48Z_{i}}{10Z_{i}}x^{2}+\frac{3(192Z_{i}^{2}+301Z_{i}\sqrt{2}+178)(552Z_{i}+241\sqrt{2})}{28000Z_{i}^{3}}\Big{]}/\triangle;$
$\displaystyle\eta_{4}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}x\Big{[}7x^{2}+\frac{173088Z_{i}^{2}+142032Z_{i}\sqrt{2}+59989}{5600Z_{i}^{2}}\Big{]}/\triangle,$
(85)

where the denominator $\triangle$ is equivalent to (82). For the particular
case of $Z_{i}=1$ these electron viscosities become

$\displaystyle\eta_{0}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\,\frac{35(241\sqrt{2}+552)}{6(301\sqrt{2}+370)};$
$\displaystyle\eta_{2}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\Big{[}\frac{33\sqrt{2}+48}{10}x^{2}+\frac{382983}{14000}\sqrt{2}+\frac{523983}{14000}\Big{]}/\triangle;$
$\displaystyle\eta_{4}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}x\Big{[}7x^{2}+\frac{8877}{350}\sqrt{2}+\frac{233077}{5600}\Big{]}/\triangle,$
(86)

with $\triangle$ equal to (83), and with numerical values

$\displaystyle\eta_{0}^{e(4)}$ $\displaystyle=$ $\displaystyle
6.546\frac{p_{e}}{\nu_{ei}};$ $\displaystyle\eta_{2}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}\big{(}9.467x^{2}+76.114\big{)}/\triangle;$
$\displaystyle\eta_{4}^{e(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{\nu_{ei}}x\big{(}7x^{2}+77.489\big{)}/\triangle;$
$\displaystyle\triangle$ $\displaystyle=$ $\displaystyle
x^{4}+13.801x^{2}+11.628.$ (87)

##  4 Generalized electron coefficients for multi-species plasmas

Here we use the idea of Simakov & Molvig (2014), and before that for example
by Zhdanov (2002) (orig publ. 1982) and Hinton (1983), who pointed out that
because of the smallness of mass ratios $m_{e}/m_{i}$, the electron
coefficients of Braginskii (1965) can be straightforwardly generalized for
multi-species plasmas. Simakov & Molvig (2014) considered unmagnetized plasmas
and provide analytic parallel coefficients $\alpha_{0}$, $\beta_{0}$,
$\gamma_{0}$ together with the parallel electron viscosity $\eta_{0}^{e}$.
Here we show that the same construction applies when a magnetic field is
present, and that all the electron coefficients provided in the previous
section can be easily generalized in the same way. One starts by considering
the general multi-species description with collisional contributions given in
Section 2.2. Because of the smallness of $m_{e}/m_{i}$, mass-ratio
coefficients for each ion species simplify into (47). One introduces an
effective ion charge together with an effective ion velocity

$Z_{\textrm{eff}}=\frac{\sum_{i}\nu_{ei}}{\sqrt{2}\nu_{ee}}=\frac{\sum_{i}n_{i}Z_{i}^{2}}{n_{e}};\qquad\langle\boldsymbol{u}_{i}\rangle_{\textrm{eff}}=\frac{\sum_{i}\nu_{ei}\boldsymbol{u}_{i}}{\sum_{i}\nu_{ei}},$
(88)

and it is straightforward to show that the collisional contributions for a one
ion-electron plasma (48), (52) are then replaced by

$\displaystyle\boldsymbol{R}_{e}$ $\displaystyle=$
$\displaystyle-\rho_{e}(\sum_{i}\nu_{ei})(\boldsymbol{u}_{e}-\langle\boldsymbol{u}_{i}\rangle_{\textrm{eff}})+\frac{21}{10}\frac{\rho_{e}}{p_{e}}(\sum_{i}\nu_{ei})\vec{\boldsymbol{q}}_{e}-\frac{3}{56}\frac{\rho_{e}^{2}}{p_{e}^{2}}(\sum_{i}\nu_{ei})\vec{\boldsymbol{X}}^{(5)}_{e};$
(89) $\displaystyle\vec{\boldsymbol{Q}}^{(3)}_{e}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle+\frac{3}{2}p_{e}(\sum_{i}\nu_{ei})(\boldsymbol{u}_{e}-\langle\boldsymbol{u}_{i}\rangle_{\textrm{eff}})-\Big{[}\frac{\sqrt{2}}{Z_{\textrm{eff}}}+\frac{19}{4}\Big{]}(\sum_{i}\nu_{ei})\vec{\boldsymbol{q}}_{e}+\Big{[}\frac{3}{70\sqrt{2}Z_{\textrm{eff}}}+\frac{69}{560}\Big{]}(\sum_{i}\nu_{ei})\frac{\rho_{e}}{p_{e}}\vec{\boldsymbol{X}}^{(5)}_{e};$
(90) $\displaystyle\vec{\boldsymbol{Q}}^{(5)}_{e}\,{}^{\prime}$
$\displaystyle=$
$\displaystyle+27\frac{p_{e}^{2}}{\rho_{e}}(\sum_{i}\nu_{ei})(\boldsymbol{u}_{e}-\langle\boldsymbol{u}_{i}\rangle_{\textrm{eff}})-\Big{[}\frac{76}{5\sqrt{2}Z_{\textrm{eff}}}+\frac{759}{10}\Big{]}(\sum_{i}\nu_{ei})\frac{p_{e}}{\rho_{e}}\vec{\boldsymbol{q}}_{e}-\Big{[}\frac{3}{35\sqrt{2}Z_{\textrm{eff}}}-\frac{533}{280}\Big{]}(\sum_{i}\nu_{ei})\vec{\boldsymbol{X}}^{(5)}_{e}.$
(91)

Contributions (90), (91) enter the right-hand-sides of electron evolution
equations (51). The system is completely the same as for the one ion-electron
plasma, if in (48), (52) the following replacement is applied

$Z_{i}\to
Z_{\textrm{eff}};\qquad\nu_{ei}\to\sum_{i}\nu_{ei};\qquad\delta\boldsymbol{u}=\boldsymbol{u}_{e}-\boldsymbol{u}_{i}\to\boldsymbol{u}_{e}-\langle\boldsymbol{u}_{i}\rangle_{\textrm{eff}}.$
(92)

If evolution equations can be obtained with the transformation (92), of course
their solution can be obtained with the same transformation as well. The same
transformation applies for the viscous evolution equations (80), (81) and
their solutions. As an example, the generalized (thermal) electron heat of
Braginskii (1965) for multi-species plasmas reads

$\displaystyle\vec{\boldsymbol{q}}_{e}^{T}$ $\displaystyle=$
$\displaystyle-\kappa_{\parallel}^{e}\nabla_{\parallel}T_{e}-\kappa_{\perp}^{e}\nabla_{\perp}T_{e}+\kappa_{\times}^{e}\boldsymbol{\hat{b}}\times\nabla
T_{e};$ (93) $\displaystyle\kappa_{\parallel}^{e}$ $\displaystyle=$
$\displaystyle\frac{p_{e}}{m_{e}(\sum_{i}\nu_{ei})}\gamma_{0};\qquad\kappa_{\perp}^{e}=\frac{p_{e}}{m_{e}(\sum_{i}\nu_{ei})}\frac{\gamma_{1}^{\prime}x^{2}+\gamma_{0}^{\prime}}{\triangle};\qquad\kappa_{\times}^{e}=\frac{p_{e}}{m_{e}(\sum_{i}\nu_{ei})}\frac{\gamma_{1}^{\prime\prime}x^{3}+\gamma_{0}^{\prime\prime}x}{\triangle};$
(94) $\displaystyle\gamma_{0}$ $\displaystyle=$
$\displaystyle\frac{25Z_{\textrm{eff}}(433Z_{\textrm{eff}}+180\sqrt{2})}{4(217Z_{\textrm{eff}}^{2}+604Z_{\textrm{eff}}\sqrt{2}+288)};\qquad\gamma_{1}^{\prime}=\frac{13Z_{\textrm{eff}}+4\sqrt{2}}{4Z_{\textrm{eff}}};\qquad\gamma_{1}^{\prime\prime}=\frac{5}{2};$
$\displaystyle\gamma_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{(217Z_{\textrm{eff}}^{2}+604Z_{\textrm{eff}}\sqrt{2}+288)(433Z_{\textrm{eff}}+180\sqrt{2})}{78400Z_{\textrm{eff}}^{3}};\qquad\gamma_{0}^{\prime\prime}=\frac{320797Z_{\textrm{eff}}^{2}+202248Z_{\textrm{eff}}\sqrt{2}+72864}{31360Z_{\textrm{eff}}^{2}};$
$\displaystyle\triangle$ $\displaystyle=$ $\displaystyle
x^{4}+\delta_{1}x^{2}+\delta_{0};\qquad\delta_{0}=\Big{(}\frac{217Z_{\textrm{eff}}^{2}+604Z_{\textrm{eff}}\sqrt{2}+288}{700Z_{\textrm{eff}}^{2}}\Big{)}^{2};$
$\displaystyle\delta_{1}$ $\displaystyle=$
$\displaystyle\frac{586601Z_{\textrm{eff}}^{2}+330152Z_{\textrm{eff}}\sqrt{2}+106016}{78400Z_{\textrm{eff}}^{2}},$
(95)

where $x=\Omega_{e}/(\sum_{i}\nu_{ei})$. With recipe (92) one obtains
generalized solutions for the frictional electron heat flux
$\vec{\boldsymbol{q}}_{e}^{u}$, together with solutions for
$\vec{\boldsymbol{X}}^{(5)}_{e}$ and viscosity-tensors
$\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}$,
$\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}$ which are not repeated here.

From the electron momentum equation, the electric field then becomes

$\displaystyle\boldsymbol{E}$ $\displaystyle=$
$\displaystyle-\frac{1}{c}\boldsymbol{u}_{e}\times\boldsymbol{B}-\frac{1}{en_{e}}\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{e}+\frac{m_{e}}{e}\boldsymbol{G}$
(96)
$\displaystyle+(\sum_{i}\nu_{ei})\Big{[}+\frac{m_{e}}{e}\big{(}\langle\boldsymbol{u}_{i}\rangle_{\textrm{eff}}-\boldsymbol{u}_{e}\big{)}+\frac{21}{10}\frac{m_{e}}{ep_{e}}\vec{\boldsymbol{q}}_{e}-\frac{3}{56}\frac{\rho_{e}^{2}}{en_{e}p_{e}^{2}}\vec{\boldsymbol{X}}^{(5)}_{e}\Big{]}-\frac{m_{e}}{e}\frac{d_{e}\boldsymbol{u}_{e}}{dt},$

and expressions for heat fluxes $\vec{\boldsymbol{q}}_{e}$ and
$\vec{\boldsymbol{X}}^{(5)}_{e}$ enter the electric field.

##  5 Generalization with coupling of stress-tensors and heat fluxes

Here we consider the coupling between viscosity-tensors and heat fluxes. Using
the semi-linear approximation and retaining the coupling, the 21-moment model
(9)-(12) simplifies into

$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}^{S}+p_{a}\bar{\bar{\boldsymbol{W}}}_{a}+\frac{2}{5}\Big{(}(\nabla\vec{\boldsymbol{q}}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\vec{\boldsymbol{q}}_{a}\Big{)}$
$\displaystyle\qquad=\bar{\bar{\boldsymbol{Q}}}_{a}^{(2)}\,{}^{\prime}=\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a};$
(97)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{q}}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{a}+\frac{5}{2}p_{a}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}+\frac{1}{2}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-\frac{5}{2}\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$
$\displaystyle\qquad=\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}=\frac{1}{2}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(3)}_{a}-\frac{5}{2}\frac{p_{a}}{\rho_{a}}\boldsymbol{R}_{a};$
(98)
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+7\frac{p_{a}^{2}}{\rho_{a}}\bar{\bar{\boldsymbol{W}}}_{a}+\frac{1}{5}\Big{[}(\nabla\vec{\boldsymbol{X}}^{(5)}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}(\nabla\cdot\vec{\boldsymbol{X}}^{(5)}_{a})\Big{]}$
$\displaystyle\qquad=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}\,{}^{\prime}=\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}\textrm{Tr}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a};$
(99)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{a}+70\frac{p_{a}^{2}}{\rho_{a}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}+18\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-98\frac{p_{a}^{2}}{\rho_{a}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$
$\displaystyle\qquad=\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime}=\textrm{Tr}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(5)}_{a}-35\frac{p_{a}^{2}}{\rho_{a}^{2}}\boldsymbol{R}_{a}.$
(100)

Terms such as $(\nabla p_{a})\vec{\boldsymbol{q}}_{a}$ were neglected and
large-scale gradients are assumed to be small (see Section 2.4). The right-
hand-sides were given in Sections 2.1 & 2.2, and for one ion-electron plasmas
in Section 3. The system now represents a generalization of Braginskii (1965),
where heat fluxes and stress-tensors are coupled. For the highest-level of
precision, one should solve dispersion relations directly with the above
system, where all the heat fluxes and stress-tensors are _independent_
variables. At the lowest-level of precision, one prescribes the quasi-static
approximation and cancels the time-derivatives $d/dt$. Nevertheless, for
sufficiently low frequencies there exists a “middle-route” procedure known
from the algebra of collisionless models, by decomposing each moment into its
first and second orders

$\displaystyle\vec{\boldsymbol{q}}_{a}$ $\displaystyle=$
$\displaystyle\vec{\boldsymbol{q}}_{a}^{(1)}+\vec{\boldsymbol{q}}_{a}^{(2)};\qquad\vec{\boldsymbol{X}}^{(5)}_{a}=\vec{\boldsymbol{X}}^{(5,1)}_{a}+\vec{\boldsymbol{X}}^{(5,2)}_{a};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ $\displaystyle=$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}+\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{a};\qquad\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}=\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}+\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{a},$
(101)

and by neglecting the time derivative of the second-order moments. One can
consider

$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}^{S}+p_{a}\bar{\bar{\boldsymbol{W}}}_{a}+\frac{2}{5}\Big{(}(\nabla\vec{\boldsymbol{q}}_{a}^{(1)})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\vec{\boldsymbol{q}}_{a}^{(1)}\Big{)}=\bar{\bar{\boldsymbol{Q}}}_{a}^{(2)}\,{}^{\prime};$
(102)
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+7\frac{p_{a}^{2}}{\rho_{a}}\bar{\bar{\boldsymbol{W}}}_{a}+\frac{1}{5}\Big{[}(\nabla\vec{\boldsymbol{X}}^{(5,1)}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}(\nabla\cdot\vec{\boldsymbol{X}}^{(5,1)}_{a})\Big{]}=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}\,{}^{\prime};$
(103)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{q}}_{a}^{(1)}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{q}}_{a}+\frac{5}{2}p_{a}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}+\frac{1}{2}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}-\frac{5}{2}\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}=\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime};$
(104)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5,1)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(5)}_{a}+70\frac{p_{a}^{2}}{\rho_{a}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}+18\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}-98\frac{p_{a}^{2}}{\rho_{a}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}=\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime},$
(105)

where the collisional contributions on the right hand sides contain full
moments $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$,
$\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$, $\vec{\boldsymbol{q}}_{a}$,
$\vec{\boldsymbol{X}}^{(5)}_{a}$. In the collisionless regime, a similar
procedure was used for example by Macmahon (1965), Mikhailovskii & Smolyakov
(1985), Ramos (2005), Goswami et al. (2005), Passot et al. (2012), and it is
well-known that retaining the time-derivatives $d/dt$ is crucial for
recovering the dispersion relation of perpendicular fast mode with respect to
kinetic theory (its wavenumber dependence in the long-wavelength limit). It is
straightforward to further increase the precision, by for example retaining
full $\vec{\boldsymbol{q}}_{a}$ & $\vec{\boldsymbol{X}}^{(5)}_{a}$ in the last
terms of (102) & (103), or by retaining full
$\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}$ &
$\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(4)}$ in the last terms of (104) & (105)
(which we do not show). The procedure and its application is described in
detail in Hunana et al. (2019b) (see Sections 5.8 & 5.9), and the coupling of
stress-tensors and heat fluxes is also crucial for the firehose instability
(see Figures 7 & 10 there; see also figures with simpler models in Hunana &
Zank (2017)).

##  6 Coupling for unmagnetized one ion-electron plasma.

We further focus on the particular case of a one ion-electron plasma with
similar temperatures. It is of course possible to algebraically solve the
entire system (102)-(105) with a magnetic field present, which will be
presented elsewhere. Here for clarity and to demonstrate our point we find it
sufficient to focus on an unmagnetized plasma. Equivalently, we thus only
consider solutions for parallel moments along the magnetic field, similarly to
the heat flux model of Spitzer & Härm (1953). For the heat flux equations
(102), (103) it is beneficial to introduce matrices

$\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{a}=(\nabla\vec{\boldsymbol{q}}_{a}^{(1)})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\vec{\boldsymbol{q}}_{a}^{(1)};\qquad\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{a}=(\nabla\vec{\boldsymbol{X}}^{(5,1)}_{a})^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\vec{\boldsymbol{X}}^{(5,1)}_{a},$
(106)

which are symmetric and traceless, analogously to matrix
$\bar{\bar{\boldsymbol{W}}}_{a}$.

###  6.1 Ion species (self-collisions)

For the ion species, the viscosity-tensors have the following form

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}$ $\displaystyle=$
$\displaystyle-\frac{1025}{1068}\frac{p_{a}}{\nu_{aa}}\bar{\bar{\boldsymbol{W}}}_{a};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{a}$ $\displaystyle=$
$\displaystyle-\frac{1}{\nu_{aa}}\Big{[}\frac{79}{534}\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{a}+\frac{3}{178}\frac{\rho_{a}}{p_{a}}\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{a}+\frac{395}{1068}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}}{\partial
t}+\frac{15}{178}\frac{\rho_{a}}{p_{a}}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}}{\partial
t}\Big{]};$ $\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}$
$\displaystyle=$
$\displaystyle-\frac{8435}{1068}\frac{p_{a}^{2}}{\rho_{a}\nu_{aa}}\bar{\bar{\boldsymbol{W}}}_{a};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{a}$ $\displaystyle=$
$\displaystyle+\frac{1}{\nu_{aa}}\Big{[}+\frac{371}{534}\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{a}-\frac{49}{178}\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{a}+\frac{1855}{1068}\frac{p_{a}}{\rho_{a}}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}}{\partial
t}-\frac{245}{178}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}}{\partial
t}\Big{]},$ (107)

and heat fluxes become

$\displaystyle\vec{\boldsymbol{q}}_{a}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{125}{32}\frac{p_{a}}{m_{a}\nu_{aa}}\nabla T_{a};$
$\displaystyle\vec{\boldsymbol{q}}_{a}^{(2)}$ $\displaystyle=$
$\displaystyle+\frac{1}{\nu_{aa}}\Big{[}+\frac{515}{96}\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}-\frac{95}{96}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}-\frac{5}{48}\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{a}}{\partial
t}-\frac{5}{96}\frac{\rho_{a}}{p_{a}}\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{a}}{\partial
t}\Big{]};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,1)}_{a}$ $\displaystyle=$
$\displaystyle-\frac{2975}{24}\frac{p_{a}^{2}}{\rho_{a}m_{a}\nu_{aa}}\nabla
T_{a};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,2)}_{a}$ $\displaystyle=$
$\displaystyle+\frac{1}{\nu_{aa}}\Big{[}+\frac{p^{2}_{a}}{\rho^{2}_{a}}\frac{13825}{72}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{a}-\frac{2485}{72}\frac{p_{a}}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{a}+\frac{665}{36}\frac{p_{a}}{\rho_{a}}\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{a}}{\partial
t}-\frac{175}{72}\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{a}}{\partial
t}\Big{]}.$ (108)

The model is fully specified and closed, and can be used in the given form.
Nevertheless, it is possible to further apply the semi-linear approximation,
in which case the viscosity corrections simplify into

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{a}$ $\displaystyle=$
$\displaystyle+\underbrace{\frac{45575}{17088}}_{2.6671}\frac{p_{a}}{m_{a}\nu_{aa}^{2}}\Big{[}2\nabla\nabla
T_{a}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla^{2}T_{a}\Big{]}+\underbrace{\frac{1164025}{1140624}}_{1.0205}\frac{p_{a}}{\nu_{aa}^{2}}\frac{\partial\bar{\bar{\boldsymbol{W}}}_{a}}{\partial
t};$ $\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{a}$ $\displaystyle=$
$\displaystyle+\underbrace{\frac{536725}{17088}}_{31.4095}\frac{p_{a}^{2}}{\rho_{a}m_{a}\nu_{aa}^{2}}\Big{[}2\nabla\nabla
T_{a}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla^{2}T_{a}\Big{]}+\underbrace{\frac{10498075}{1140624}}_{9.2038}\frac{p_{a}^{2}}{\rho_{a}\nu_{aa}^{2}}\frac{\partial\bar{\bar{\boldsymbol{W}}}_{a}}{\partial
t},$ (109)

and the heat flux corrections become

$\displaystyle\vec{\boldsymbol{q}}_{a}^{(2)}$ $\displaystyle=$
$\displaystyle+\underbrace{\frac{45575}{17088}}_{2.6671}\frac{p_{a}^{2}}{\rho_{a}\nu_{aa}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{a}+\underbrace{\frac{31625}{4608}}_{6.8631}\frac{p_{a}}{m_{a}\nu_{aa}^{2}}\frac{\partial\nabla
T_{a}}{\partial t};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,2)}_{a}$
$\displaystyle=$
$\displaystyle+\underbrace{\frac{1131725}{12816}}_{88.3056}\frac{p_{a}^{3}}{\rho_{a}^{2}\nu_{aa}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{a}+\underbrace{\frac{791875}{3456}}_{229.1305}\frac{p_{a}^{2}}{\rho_{a}m_{a}\nu_{aa}^{2}}\frac{\partial\nabla
T_{a}}{\partial t}.$ (110)

###  6.2 Electron species

For the electron species, it is useful to introduce denominator

$D_{1}=192Z_{i}^{2}+301\sqrt{2}Z_{i}+178,$ (111)

and solutions for the stress-tensors are

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}$ $\displaystyle=$
$\displaystyle-\frac{5Z_{i}(205\sqrt{2}+408Z_{i})}{6D_{1}}\frac{p_{e}}{\nu_{ei}}\bar{\bar{\boldsymbol{W}}}_{e};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{e}$ $\displaystyle=$
$\displaystyle-\frac{1}{D_{1}\nu_{ei}}\Big{[}\frac{Z_{i}}{3}(79\sqrt{2}-96Z_{i})\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{e}+3Z_{i}(\sqrt{2}+4Z_{i})\frac{\rho_{e}}{p_{e}}\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{e}$
$\displaystyle\qquad+\frac{5}{6}Z_{i}(79\sqrt{2}-96Z_{i})\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}}{\partial
t}+15Z_{i}(\sqrt{2}+4Z_{i})\frac{\rho_{e}}{p_{e}}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}}{\partial
t}\Big{]};$ $\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}$
$\displaystyle=$
$\displaystyle-\frac{35Z_{i}(241\sqrt{2}+552Z_{i})}{6D_{1}}\frac{p_{e}^{2}}{\rho_{e}\nu_{ei}}\bar{\bar{\boldsymbol{W}}}_{e};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{e}$ $\displaystyle=$
$\displaystyle+\frac{1}{D_{1}\nu_{ei}}\Big{[}\frac{7}{3}Z_{i}(53\sqrt{2}+288Z_{i})\frac{p_{e}}{\rho_{e}}\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{e}-7Z_{i}(7\sqrt{2}+20Z_{i})\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{e}$
(112)
$\displaystyle\qquad+\frac{35}{6}Z_{i}(53\sqrt{2}+288Z_{i})\frac{p_{e}}{\rho_{e}}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}}{\partial
t}-35Z_{i}(7\sqrt{2}+20Z_{i})\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}}{\partial
t}\Big{]},$

with matrices $\bar{\bar{\boldsymbol{Y}}}_{e}$ defined by (106). For the heat
fluxes it is useful to define denominator

$D_{2}=217Z_{i}^{2}+604Z_{i}\sqrt{2}+288,$ (113)

together with $\delta\boldsymbol{u}=\boldsymbol{u}_{e}-\boldsymbol{u}_{i}$,
and the results read

$\displaystyle\vec{\boldsymbol{q}}_{e}^{(1)}$ $\displaystyle=$
$\displaystyle-\frac{25Z_{i}(180\sqrt{2}+433Z_{i})}{4D_{2}}\frac{p_{e}}{m_{e}\nu_{ei}}\nabla
T_{e}+\frac{30Z_{i}(15\sqrt{2}+11Z_{i})}{D_{2}}p_{e}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{q}}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle+\frac{1}{D_{2}\nu_{ei}}\Big{[}\frac{5}{4}Z_{i}(1236\sqrt{2}+4097Z_{i})\frac{p_{e}}{\rho_{e}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}-\frac{5}{4}Z_{i}(228\sqrt{2}+709Z_{i})\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}$
$\displaystyle\qquad-\frac{5}{2}Z_{i}(12\sqrt{2}-533Z_{i})\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{e}}{\partial
t}-\frac{15}{4}Z_{i}(4\sqrt{2}+23Z_{i})\frac{\rho_{e}}{p_{e}}\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{e}}{\partial
t}\Big{]};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,1)}_{e}$ $\displaystyle=$
$\displaystyle-\frac{175Z_{i}(204\sqrt{2}+571Z_{i})}{D_{2}}\frac{p_{e}^{2}}{\rho_{e}m_{e}\nu_{ei}}\nabla
T_{e}+\frac{840Z_{i}(13\sqrt{2}+12Z_{i})}{D_{2}}\frac{p^{2}_{e}}{\rho_{e}}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{X}}^{(5,2)}_{e}$ $\displaystyle=$
$\displaystyle\frac{1}{D_{2}\nu_{ei}}\Big{[}+175Z_{i}(316\sqrt{2}+1103Z_{i})\frac{p_{e}^{2}}{\rho_{e}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}-35Z_{i}(284\sqrt{2}+951Z_{i})\frac{p_{e}}{\rho_{e}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}$
(114)
$\displaystyle\qquad+70Z_{i}(76\sqrt{2}+759Z_{i})\frac{p_{e}}{\rho_{e}}\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{e}}{\partial
t}-175Z_{i}(4\sqrt{2}+19Z_{i})\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{e}}{\partial
t}\Big{]}.$

The system is now fully specified and can be used in this form. For the
particular case of $Z_{i}=1$ numerical values become

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}$ $\displaystyle=$
$\displaystyle-0.7309\frac{p_{e}}{\nu_{ei}}\bar{\bar{\boldsymbol{W}}}_{e};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{e}$ $\displaystyle=$
$\displaystyle-\frac{1}{\nu_{ei}}\Big{[}0.006587\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{e}+0.02041\frac{\rho_{e}}{p_{e}}\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{e}+0.01647\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}}{\partial
t}+0.1021\frac{\rho_{e}}{p_{e}}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}}{\partial
t}\Big{]};$ $\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}$
$\displaystyle=$
$\displaystyle-6.5455\frac{p_{e}^{2}}{\rho_{e}\nu_{ei}}\bar{\bar{\boldsymbol{W}}}_{e};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{e}$ $\displaystyle=$
$\displaystyle+\frac{1}{\nu_{ei}}\Big{[}1.0644\frac{p_{e}}{\rho_{e}}\bar{\bar{\boldsymbol{Y}}}^{(3,1)}_{e}-0.2630\bar{\bar{\boldsymbol{Y}}}^{(5,1)}_{e}+2.6609\frac{p_{e}}{\rho_{e}}\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}}{\partial
t}-1.3152\frac{\partial\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}}{\partial
t}\Big{]};$ (115) $\displaystyle\vec{\boldsymbol{q}}_{e}^{(1)}$
$\displaystyle=$ $\displaystyle-3.1616\frac{p_{e}}{m_{e}\nu_{ei}}\nabla
T_{e}+0.7110p_{e}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{q}}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle+\frac{1}{\nu_{ei}}\Big{[}5.3754\frac{p_{e}}{\rho_{e}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}-0.9486\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}+0.9492\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{e}}{\partial
t}-0.07906\frac{\rho_{e}}{p_{e}}\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{e}}{\partial
t}\Big{]};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,1)}_{e}$ $\displaystyle=$
$\displaystyle-110.664\frac{p_{e}^{2}}{\rho_{e}m_{e}\nu_{ei}}\nabla
T_{e}+18.7783\frac{p^{2}_{e}}{\rho_{e}}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{X}}^{(5,2)}_{e}$ $\displaystyle=$
$\displaystyle\frac{1}{\nu_{ei}}\Big{[}199.554\frac{p_{e}^{2}}{\rho_{e}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}-34.831\frac{p_{e}}{\rho_{e}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}+44.625\frac{p_{e}}{\rho_{e}}\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{e}}{\partial
t}-3.1747\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{e}}{\partial t}\Big{]}.$
(116)

By further applying the quasi-linear approximation, corrections to the
electron viscosities become

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{e}$ $\displaystyle=$
$\displaystyle+\frac{25Z_{i}^{2}(119520Z_{i}^{2}+101784Z_{i}\sqrt{2}+46561)}{18D_{1}^{2}}\frac{p_{e}}{\nu_{ei}^{2}}\frac{\partial\bar{\bar{\boldsymbol{W}}}_{e}}{\partial
t}$ (117)
$\displaystyle-\frac{10Z_{i}^{2}(11040Z_{i}^{2}+15557Z_{i}\sqrt{2}+8922)}{D_{2}D_{1}}\frac{p_{e}}{\nu_{ei}}\Big{[}\big{(}\nabla\delta\boldsymbol{u}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\delta\boldsymbol{u}\Big{]}$
$\displaystyle+\frac{25Z_{i}^{2}(534000Z_{i}^{2}+366451Z_{i}\sqrt{2}+131256)}{12D_{2}D_{1}}\frac{p_{e}}{\nu_{ei}^{2}m_{e}}\Big{[}2\nabla\nabla
T_{e}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla^{2}T_{e}\Big{]};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{e}$ $\displaystyle=$
$\displaystyle+\frac{175Z_{i}^{2}(173088Z_{i}^{2}+142032Z_{i}\sqrt{2}+59989)}{18D_{1}^{2}}\frac{p_{e}^{2}}{\nu_{ei}^{2}\rho_{e}}\frac{\partial\bar{\bar{\boldsymbol{W}}}_{e}}{\partial
t}$ (118)
$\displaystyle-\frac{70Z_{i}^{2}(16992Z_{i}^{2}+23993Z_{i}\sqrt{2}+13698)}{D_{2}D_{1}}\frac{p_{e}^{2}}{\nu_{ei}\rho_{e}}\Big{[}\big{(}\nabla\delta\boldsymbol{u}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\delta\boldsymbol{u}\Big{]}$
$\displaystyle+\frac{175Z_{i}^{2}(834576Z_{i}^{2}+603679Z_{i}\sqrt{2}+220824)}{12D_{2}D_{1}}\frac{p_{e}^{2}}{\nu_{ei}^{2}m_{e}\rho_{e}}\Big{[}2\nabla\nabla
T_{e}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla^{2}T_{e}\Big{]};$

together with corrections for the heat fluxes

$\displaystyle\vec{\boldsymbol{q}}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{25Z_{i}^{2}(534000Z_{i}^{2}+366451Z_{i}\sqrt{2}+131256)}{12D_{1}D_{2}}\frac{p_{e}^{2}}{\rho_{e}\nu_{ei}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}$
$\displaystyle-\frac{75Z_{i}^{2}(5729Z_{i}^{2}+6711Z_{i}\sqrt{2}+4728)}{D_{2}^{2}}\frac{p_{e}}{\nu_{ei}}\frac{\partial\delta\boldsymbol{u}}{\partial
t}$
$\displaystyle+\frac{125Z_{i}^{2}(320797Z_{i}^{2}+202248Z_{i}\sqrt{2}+72864)}{8D_{2}^{2}}\frac{p_{e}}{\nu_{ei}^{2}m_{e}}\frac{\partial\nabla
T_{e}}{\partial t};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,2)}_{e}$
$\displaystyle=$
$\displaystyle\frac{175Z_{i}^{2}(712272Z_{i}^{2}+463249Z_{i}\sqrt{2}+155208)}{3D_{1}D_{2}}\frac{p_{e}^{3}}{\rho_{e}^{2}\nu_{ei}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}$
(119)
$\displaystyle-\frac{2100Z_{i}^{2}(7611Z_{i}^{2}+8429Z_{i}\sqrt{2}+5000)}{D_{2}^{2}}\frac{p_{e}^{2}}{\nu_{ei}\rho_{e}}\frac{\partial\delta\boldsymbol{u}}{\partial
t}$
$\displaystyle+\frac{875Z_{i}^{2}(430783Z_{i}^{2}+261672Z_{i}\sqrt{2}+86880)}{2D_{2}^{2}}\frac{p_{e}^{2}}{\nu_{ei}^{2}\rho_{e}m_{e}}\frac{\partial\nabla
T_{e}}{\partial t}.$

For the ion charge $Z_{i}=1$, the numerical values read

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2,2)}_{e}$ $\displaystyle=$
$\displaystyle+0.6801\frac{p_{e}}{\nu_{ei}^{2}}\frac{\partial\bar{\bar{\boldsymbol{W}}}_{e}}{\partial
t}-0.3880\frac{p_{e}}{\nu_{ei}}\Big{[}\big{(}\nabla\delta\boldsymbol{u}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\delta\boldsymbol{u}\Big{]}$
(120) $\displaystyle+2.2799\frac{p_{e}}{\nu_{ei}^{2}m_{e}}\Big{[}2\nabla\nabla
T_{e}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla^{2}T_{e}\Big{]};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4,2)}_{e}$ $\displaystyle=$
$\displaystyle+6.6638\frac{p_{e}^{2}}{\nu_{ei}^{2}\rho_{e}}\frac{\partial\bar{\bar{\boldsymbol{W}}}_{e}}{\partial
t}-4.1827\frac{p_{e}^{2}}{\nu_{ei}\rho_{e}}\Big{[}\big{(}\nabla\delta\boldsymbol{u}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\delta\boldsymbol{u}\Big{]}$
(121)
$\displaystyle+25.7440\frac{p_{e}^{2}}{\nu_{ei}^{2}m_{e}\rho_{e}}\Big{[}2\nabla\nabla
T_{e}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla^{2}T_{e}\Big{]};$

together with

$\displaystyle\vec{\boldsymbol{q}}_{e}^{(2)}$ $\displaystyle=$ $\displaystyle
2.2799\frac{p_{e}^{2}}{\rho_{e}\nu_{ei}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}-0.8098\frac{p_{e}}{\nu_{ei}}\frac{\partial\delta\boldsymbol{u}}{\partial
t}+5.7487\frac{p_{e}}{\nu_{ei}^{2}m_{e}}\frac{\partial\nabla T_{e}}{\partial
t};$ $\displaystyle\vec{\boldsymbol{X}}^{(5,2)}_{e}$ $\displaystyle=$
$\displaystyle
82.1278\frac{p_{e}^{3}}{\rho_{e}^{2}\nu_{ei}^{2}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}-27.8859\frac{p_{e}^{2}}{\nu_{ei}\rho_{e}}\frac{\partial\delta\boldsymbol{u}}{\partial
t}+210.2318\frac{p_{e}^{2}}{\nu_{ei}^{2}\rho_{e}m_{e}}\frac{\partial\nabla
T_{e}}{\partial t}.$ (122)

The rate-of-strain tensor $\bar{\bar{\boldsymbol{W}}}_{e}$ obviously enters
the electron heat fluxes, even in a quasi-static approximation.

###  6.3 Momentum exchange rates

Collisional momentum exchange rates $\boldsymbol{R}_{e}=-\boldsymbol{R}_{i}$
given by (48) can also be split into the first and second order
$\boldsymbol{R}_{e}=\boldsymbol{R}_{e}^{(1)}+\boldsymbol{R}_{e}^{(2)}$,
according to

$\displaystyle\boldsymbol{R}_{e}^{(1)}$ $\displaystyle=$
$\displaystyle-\rho_{e}\nu_{ei}\delta\boldsymbol{u}+\frac{21}{10}\frac{\rho_{e}}{p_{e}}\nu_{ei}\vec{\boldsymbol{q}}_{e}^{(1)}-\frac{3}{56}\frac{\rho_{e}^{2}}{p_{e}^{2}}\nu_{ei}\vec{\boldsymbol{X}}^{(5,1)}_{e};$
$\displaystyle\boldsymbol{R}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle+\frac{21}{10}\frac{\rho_{e}}{p_{e}}\nu_{ei}\vec{\boldsymbol{q}}_{e}^{(2)}-\frac{3}{56}\frac{\rho_{e}^{2}}{p_{e}^{2}}\nu_{ei}\vec{\boldsymbol{X}}^{(5,2)}_{e}.$
(123)

Then by using results given in the previous section

$\displaystyle\boldsymbol{R}_{e}^{(1)}$ $\displaystyle=$
$\displaystyle-\nu_{ei}\rho_{e}\frac{(D_{2}-153Z_{i}^{2}-360Z_{i}\sqrt{2})}{D_{2}}\delta\boldsymbol{u}-\frac{30Z_{i}(15\sqrt{2}+11Z_{i})}{D_{2}}n_{e}\nabla
T_{e};$ $\displaystyle\boldsymbol{R}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle+\frac{6Z_{i}(47\sqrt{2}+69Z_{i})}{D_{2}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}-\frac{6Z_{i}(11\sqrt{2}+13Z_{i})}{D_{2}}\frac{\rho_{e}}{p_{e}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}$
(124)
$\displaystyle-\frac{12Z_{i}(29\sqrt{2}+4Z_{i})}{D_{2}}\frac{\rho_{e}}{p_{e}}\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{e}}{\partial
t}+\frac{3Z_{i}(2\sqrt{2}-Z_{i})}{D_{2}}\frac{\rho_{e}^{2}}{p_{e}^{2}}\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{e}}{\partial
t},$

or for a particular case of $Z_{i}=1$

$\displaystyle\boldsymbol{R}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle+0.5980\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2,1)}_{e}-0.1261\frac{\rho_{e}}{p_{e}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4,1)}_{e}$
(125)
$\displaystyle-0.3974\frac{\rho_{e}}{p_{e}}\frac{\partial\vec{\boldsymbol{q}}^{(1)}_{e}}{\partial
t}+0.004036\frac{\rho_{e}^{2}}{p_{e}^{2}}\frac{\partial\vec{\boldsymbol{X}}^{(5,1)}_{e}}{\partial
t}.$

Finally, at a semi-linear level

$\displaystyle\boldsymbol{R}_{e}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{10Z_{i}^{2}(11040Z_{i}^{2}+15557Z_{i}\sqrt{2}+8922)}{D_{1}D_{2}}\frac{p_{e}}{\nu_{ei}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}$
(126)
$\displaystyle-\frac{720Z_{i}^{2}(64Z_{i}^{2}+151Z_{i}\sqrt{2}+253)}{D_{2}^{2}}\rho_{e}\frac{\partial(\delta\boldsymbol{u})}{\partial
t}$