$\displaystyle+\frac{75Z_{i}^{2}(5729Z_{i}^{2}+6711Z_{i}\sqrt{2}+4728)}{D_{2}^{2}}\frac{n_{e}}{\nu_{ei}}\frac{\partial\nabla
T_{e}}{\partial t},$

and for $Z_{i}=1$ the full momentum exchange rates become

$\displaystyle\boldsymbol{R}_{e}$ $\displaystyle=$
$\displaystyle-0.5129\nu_{ei}\rho_{e}\delta\boldsymbol{u}-0.7110n_{e}\nabla
T_{e}+0.3880\frac{p_{e}}{\nu_{ei}}\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}-0.2068\rho_{e}\frac{\partial(\delta\boldsymbol{u})}{\partial
t}+0.8098\frac{n_{e}}{\nu_{ei}}\frac{\partial\nabla T_{e}}{\partial t},$ (127)

where $\delta\boldsymbol{u}=\boldsymbol{u}_{e}-\boldsymbol{u}_{i}$. Only the
first two terms of (127) were considered by Braginskii (1965) and Spitzer &
Härm (1953) (the latter having slightly different proportionality constants;
see Appendix I). A further generalization by keeping the full
$\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}$ &
$\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(4)}$ viscosity tensors in the last terms
of (104) & (105) brings another 3 terms to $\boldsymbol{R}_{e}$ (not shown).
Naturally, in a highly-collisional regime ($\nu_{ei}\gg\omega$) all additional
terms are small in comparison to the first two terms of (127). Nevertheless,
at higher frequencies (shorter length-scales) these additional contributions
might become significant. Interestingly, the rate-of-strain tensor
$\bar{\bar{\boldsymbol{W}}}_{e}$ enters the momentum exchange rates (even at
the linear level), with contribution
$\nabla\cdot\bar{\bar{\boldsymbol{W}}}_{e}=\nabla^{2}\boldsymbol{u}_{e}+(1/3)\nabla(\nabla\cdot\boldsymbol{u}_{e})$.
Note that some terms are proportional to $1/\nu_{ei}$ and become unbounded
(divergent) in a regime of low-collisionality, which is a consequence of the
expansion procedure (i.e. a quasi-static approximation). Evolution equations
(97)-(100) are of course well-defined in the regime of low collisionality.

##  7 Multi-fluid 22-moment model

Here we consider a natural generalization of the 21-moment model, by
accounting for a fully contracted perturbation of the 4th-order fluid moment
$X_{ijkl}^{a(4)}=m_{a}\int c_{i}^{a}c_{j}^{a}c_{k}^{a}c_{l}^{a}f_{a}d^{3}v$.
The fully contracted (scalar) moment is decomposed into its Maxwellian core
and a perturbation $\widetilde{X}^{(4)}_{a}$ (denoted with tilde), according
to

$X^{(4)}_{a}=m_{a}\int|\boldsymbol{c}_{a}|^{4}f_{a}d^{3}v=15\frac{p^{2}_{a}}{\rho_{a}}+\widetilde{X}^{(4)}_{a},$
(128)

meaning a definition
$\widetilde{X}^{(4)}_{a}=m_{a}\int|\boldsymbol{c}_{a}|^{4}(f_{a}-f_{a}^{(0)})d^{3}v$,
where $f_{a}^{(0)}$ is Maxwellian. The scalar perturbation
$\widetilde{X}^{(4)}_{a}$ enters the decomposition of the 4th-order moment

$\displaystyle X_{ijkl}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{1}{15}\Big{(}15\frac{p_{a}^{2}}{\rho_{a}}+\widetilde{X}^{(4)}_{a}\Big{)}\big{(}\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\big{)}$
(129)
$\displaystyle+\frac{1}{7}\Big{[}\Pi_{ij}^{a(4)}\delta_{kl}+\Pi_{ik}^{a(4)}\delta_{jl}+\Pi_{il}^{a(4)}\delta_{jk}+\Pi_{jk}^{a(4)}\delta_{il}+\Pi_{jl}^{a(4)}\delta_{ik}+\Pi_{kl}^{a(4)}\delta_{ij}\Big{]}+\sigma_{ijkl}^{a(4)}\,{}^{\prime},$

where we neglect the traceless tensor $\sigma_{ijkl}^{a(4)}\,{}^{\prime}$, and
the entire model now represents the 22-moment model. The fully non-linear
model is given by evolution equations (6)-(9) which are unchanged, together
with

$\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}{6}\nabla\widetilde{X}^{(4)}_{a}+\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};$
(130)
$\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})-\,\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+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}\big{)}^{S}+\frac{7}{15}\big{(}15\frac{p_{a}^{2}}{\rho_{a}}+\widetilde{X}^{(4)}_{a}\big{)}\bar{\bar{\boldsymbol{W}}}_{a}$
$\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{]};$
(131)
$\displaystyle\frac{d_{a}}{dt}\widetilde{X}^{(4)}_{a}+\nabla\cdot\vec{\boldsymbol{X}}^{(5)}_{a}-20\frac{p_{a}}{\rho_{a}}\nabla\cdot\vec{\boldsymbol{q}}_{a}+\frac{7}{3}\widetilde{X}^{(4)}_{a}(\nabla\cdot\boldsymbol{u}_{a})+4\big{(}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-5\frac{p_{a}}{\rho_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}\big{)}:\nabla\boldsymbol{u}_{a}$
$\displaystyle-\frac{8}{\rho_{a}}(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})\cdot\vec{\boldsymbol{q}}_{a}=\widetilde{Q}^{(4)}_{a}\,{}^{\prime}\equiv\textrm{Tr}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}-20\frac{p_{a}}{\rho_{a}}Q_{a}-\frac{8}{\rho_{a}}\boldsymbol{R}_{a}\cdot\vec{\boldsymbol{q}}_{a};$
(132)
$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(5)}_{a}+\frac{1}{3}\nabla\widetilde{X}^{(6)}_{a}+\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(6)}_{a}$
$\displaystyle+\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{7}{3\rho_{a}}\big{(}\nabla\cdot\bar{\bar{\boldsymbol{p}}}^{a}\big{)}\widetilde{X}^{(4)}_{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\vec{\boldsymbol{Q}}^{(5)}_{a}-35\frac{p_{a}^{2}}{\rho_{a}^{2}}\boldsymbol{R}_{a}-\frac{7}{3\rho_{a}}\boldsymbol{R}_{a}\widetilde{X}^{(4)}_{a}-\frac{4}{\rho_{a}}\boldsymbol{R}_{a}\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}.$
(133)

The last equation (133) is closed with closure (13) for the stress-tensor
$\bar{\bar{\boldsymbol{\Pi}}}^{(6)}_{a}$, together with a closure for the
scalar perturbation (derived from a Hermite closure)

$\widetilde{X}^{(6)}_{a}=m_{a}\int|\boldsymbol{c}_{a}|^{6}(f_{a}-f_{a}^{(0)})d^{3}v=21\frac{p_{a}}{\rho_{a}}\widetilde{X}^{(4)}_{a}.$
(134)

In the semi-linear approximations the 22-moment model reads

$\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};$ (135)
$\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+\frac{1}{6}\nabla\widetilde{X}^{(4)}_{a}=\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime};$
(136)
$\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};$ (137)
$\displaystyle\frac{d_{a}}{dt}\widetilde{X}^{(4)}_{a}+\nabla\cdot\vec{\boldsymbol{X}}^{(5)}_{a}-20\frac{p_{a}}{\rho_{a}}\nabla\cdot\vec{\boldsymbol{q}}_{a}=\widetilde{Q}^{(4)}_{a}\,{}^{\prime}=\textrm{Tr}\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}-20\frac{p_{a}}{\rho_{a}}Q_{a};$
(138)
$\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+7\frac{p_{a}}{\rho_{a}}\nabla\widetilde{X}^{(4)}_{a}=\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime}.$
(139)

As discussed in Section 2.4, in the semi-linear approximation we are
neglecting terms such as $(\nabla p_{a})\vec{\boldsymbol{q}}_{a}$ which might
become significant in the presence of large-scale gradients, together with
other terms that are neglected. In comparison to the 21-moment model given by
(97)-(100), evolution equations (135) & (137) for stress-tensors
$\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$ &
$\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}$ remain unchanged. Importantly,
collisional contributions $\boldsymbol{R}_{a}$;
$\bar{\bar{\boldsymbol{Q}}}_{a}^{(2)}\,{}^{\prime}$;
$\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}$;
$\bar{\bar{\boldsymbol{Q}}}_{a}^{(4)}\,{}^{\prime}$;
$\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime}$ given in Section 2.1 remain
unchanged as well. The only differences are: 1) scalar perturbations
$\widetilde{X}^{(4)}_{a}$ now enter the left hand sides of evolution equations
(136) & (139) for heat fluxes $\vec{\boldsymbol{q}}_{a}$ &
$\vec{\boldsymbol{X}}^{(5)}_{a}$; 2) a new evolution equation (138) for scalar
$\widetilde{X}^{(4)}_{a}$ is present, with collisional contributions
$\widetilde{Q}^{(4)}_{a}\,{}^{\prime}$ that need to be specified; 3) the
energy exchange rates $Q_{a}$ entering the scalar pressure equation (8) are
modified, and given below.

###  7.1 Collisional contributions (arbitrary temperatures)

The energy exchange rates entering equation (8) are now given by

$\displaystyle Q_{a}=\sum_{b\neq a}Q_{ab}=\sum_{b\neq
a}\frac{\rho_{a}\nu_{ab}}{(m_{a}+m_{b})}\Big{\\{}3(T_{b}-T_{a})+\hat{P}_{ab(1)}\frac{\rho_{a}}{n_{a}p_{a}}\widetilde{X}^{(4)}_{a}-\hat{P}_{ab(2)}\frac{\rho_{b}}{n_{b}p_{b}}\widetilde{X}^{(4)}_{b}\Big{\\}},$
(140)

with mass-ratio coefficients

$\displaystyle\hat{P}_{ab(1)}=\frac{3T_{a}m_{b}(5T_{b}m_{b}+4T_{b}m_{a}-T_{a}m_{b})}{40(T_{a}m_{b}+T_{b}m_{a})^{2}};\qquad\hat{P}_{ab(2)}=\frac{3T_{b}m_{a}(5T_{a}m_{a}+4T_{a}m_{b}-T_{b}m_{a})}{40(T_{a}m_{b}+T_{b}m_{a})^{2}}.$
(141)

Interestingly, scalar perturbations $\widetilde{X}^{(4)}_{a}$ thus enter the
energy exchange rates. For self-collisions all the contributions naturally
disappear. As discussed also later in Section 8.1, for multi-fluid models the
conservation of energy
$Q_{ab}+Q_{ba}=(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\cdot\boldsymbol{R}_{ab}$
is satisfied only approximately, because in the semi-linear approximation the
differences in drifts $\boldsymbol{u}_{b}-\boldsymbol{u}_{a}$ are assumed to
be small, meaning $Q_{ab}+Q_{ba}=0$ holds. To satisfy the energy conservation
exactly, the collisional integrals would have to be calculated non-linearly
with unrestriced drifts (i.e. with the runaway effect). Nevertheless, for a
plasma consisting of only two species (such as a one ion-electron plasma) the
conservation of energy can be imposed by hand, by calculating $Q_{ab}$
according to (140), (141) and prescribing
$Q_{ba}=-Q_{ab}+(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\cdot\boldsymbol{R}_{ab}$.

Collisional exchange rates entering evolution equation (138) are given by

$\displaystyle\widetilde{Q}^{(4)}_{a}\,{}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{4}{5}\nu_{aa}\widetilde{X}^{(4)}_{a}+\sum_{b\neq
a}\nu_{ab}\Big{\\{}-\frac{p_{a}^{2}}{\rho_{a}}\frac{(T_{b}-T_{a})}{T_{a}}\hat{S}_{ab(0)}-\widetilde{X}^{(4)}_{a}\hat{S}_{ab(1)}+\frac{p_{a}^{2}\rho_{b}}{p_{b}^{2}\rho_{a}}\widetilde{X}^{(4)}_{b}\hat{S}_{ab(2)}\Big{\\}},$
(142)

with mass-ratio coefficients

$\displaystyle\hat{S}_{ab(0)}$ $\displaystyle=$
$\displaystyle\frac{36T_{a}m_{a}m_{b}}{(T_{a}m_{b}+T_{b}m_{a})(m_{b}+m_{a})};$
$\displaystyle\hat{S}_{ab(1)}$ $\displaystyle=$
$\displaystyle-\,\Big{\\{}m_{a}\big{(}17T_{a}^{3}m_{b}^{3}-36T_{a}^{2}T_{b}m_{a}m_{b}^{2}-69T_{a}^{2}T_{b}m_{b}^{3}+12T_{a}T_{b}^{2}m_{a}^{2}m_{b}-48T_{a}T_{b}^{2}m_{a}m_{b}^{2}$
$\displaystyle\quad-40T_{b}^{3}m_{a}^{3}-84T_{b}^{3}m_{a}^{2}m_{b}\big{)}\Big{\\}}\Big{[}10(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{b}+m_{a})\Big{]}^{-1};$
$\displaystyle\hat{S}_{ab(2)}$ $\displaystyle=$
$\displaystyle\frac{3T_{b}^{2}m_{a}^{2}m_{b}(7T_{a}m_{a}+4T_{a}m_{b}-3T_{b}m_{a})}{2(T_{a}m_{b}+T_{b}m_{a})^{3}(m_{b}+m_{a})},$
(143)

where the self-collisional contributions are represented by the first term of
(142).

####  7.1.1 Small temperature differences

For small temperature differences the mass-ratio coefficients become

$\displaystyle\hat{P}_{ab(1)}$ $\displaystyle=$
$\displaystyle\frac{3m_{b}}{10(m_{b}+m_{a})};\qquad\hat{P}_{ab(2)}=\frac{3m_{a}}{10(m_{b}+m_{a})};$
$\displaystyle\hat{S}_{ab(0)}$ $\displaystyle=$
$\displaystyle\frac{36m_{a}m_{b}}{(m_{b}+m_{a})^{2}};\qquad\hat{S}_{ab(1)}=\frac{2m_{a}(10m_{a}^{2}+8m_{a}m_{b}+13m_{b}^{2})}{5(m_{b}+m_{a})^{3}};\qquad\hat{S}_{ab(2)}=\frac{6m_{a}^{2}m_{b}}{(m_{b}+m_{a})^{3}},$
(144)

and for example for self-collisions $\hat{S}_{aa(1)}=31/20$ and
$\hat{S}_{aa(2)}=3/4$. We further consider a one ion-electron plasma.

###  7.2 Ion species (self-collisions)

In a quasi-static approximation the solution of equation (138) becomes

$\widetilde{X}^{(4)}_{a}=-\,\frac{5}{4\nu_{aa}}\Big{[}\nabla\cdot\vec{\boldsymbol{X}}^{(5)}_{a}-20\frac{p_{a}}{\rho_{a}}\nabla\cdot\vec{\boldsymbol{q}}_{a}\Big{]}.$
(145)

The quasi-static solution is thus completely determined by the heat fluxes
$\vec{\boldsymbol{q}}_{a}$ & $\vec{\boldsymbol{X}}^{(5)}_{a}$ and for a
magnetized plasma it has the following form

$\displaystyle\widetilde{X}^{(4)}_{a}$ $\displaystyle=$
$\displaystyle-\,\frac{5}{4\nu_{aa}}\Big{\\{}\nabla\cdot\Big{[}\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{)}\Big{]}$ (146)
$\displaystyle\quad-20\frac{p_{a}}{\rho_{a}}\nabla\cdot\Big{(}-\kappa_{\parallel}^{a}\nabla_{\parallel}T_{a}-\kappa_{\perp}^{a}\nabla_{\perp}T_{a}+\kappa_{\times}^{a}\boldsymbol{\hat{b}}\times\nabla
T_{a}\Big{)}\Big{\\}},$

where the thermal conductivities are given by (43), (46).

It feels natural to define thermal conductivities (of the moment
$\widetilde{X}^{(4)}_{a}$)

$\kappa_{\parallel}^{a(4)}=\frac{5}{4}\big{(}\kappa_{\parallel}^{a(5)}-20\kappa_{\parallel}^{a}\big{)};\qquad\kappa_{\perp}^{a(4)}=\frac{5}{4}\big{(}\kappa_{\perp}^{a(5)}-20\kappa_{\perp}^{a}\big{)};\qquad\kappa_{\times}^{a(4)}=\frac{5}{4}\big{(}\kappa_{\times}^{a(5)}-20\kappa_{\times}^{a}\big{)},$
(147)

and result (146) then transforms into

$\displaystyle\widetilde{X}^{(4)}_{a}$ $\displaystyle=$
$\displaystyle-\,\frac{p_{a}}{\nu_{aa}\rho_{a}}\nabla\cdot\Big{[}-\kappa_{\parallel}^{a(4)}\nabla_{\parallel}T_{a}-\kappa_{\perp}^{a(4)}\nabla_{\perp}T_{a}+\kappa_{\times}^{a(4)}\boldsymbol{\hat{b}}\times\nabla
T_{a}\Big{]}$ (148)
$\displaystyle-\frac{5}{4\nu_{aa}}\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{)}\cdot\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)},$

with thermal conductivities

$\displaystyle{\kappa}_{\parallel}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{1375}{24}\frac{p_{a}}{\nu_{aa}m_{a}};$
$\displaystyle{\kappa}_{\perp}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\frac{5x^{2}+(9504/245)}{x^{4}+(3313/1225)x^{2}+(20736/30625)};$
$\displaystyle{\kappa}_{\times}^{a(4)}$ $\displaystyle=$
$\displaystyle\frac{p_{a}}{\nu_{aa}m_{a}}\frac{25x^{3}+(3810/49)x}{x^{4}+(3313/1225)x^{2}+(20736/30625)}.$
(149)

The second term of (148) is strictly non-linear and may be neglected for
simplicity. The solution for $\widetilde{X}^{(4)}_{a}$ thus can be written as
a divergence of a heat flux vector defined by the expression in the square
brackets of (148). We have used 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$.

The result (148) can be further simplified in the semi-linear approximation,
where one may use $\nabla\cdot(\boldsymbol{\hat{b}}\times\nabla T_{a})=0$ and
so

$\widetilde{X}^{(4)}_{a}=+\frac{p_{a}}{\nu_{aa}\rho_{a}}\Big{[}\kappa_{\parallel}^{a(4)}\nabla_{\parallel}^{2}T_{a}+\kappa_{\perp}^{a(4)}\nabla_{\perp}^{2}T_{a}\Big{]},$
(150)

with scalars
$\nabla_{\parallel}^{2}=\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}:\nabla\nabla$
and $\nabla_{\perp}^{2}=\nabla^{2}-\nabla_{\parallel}^{2}$, and for zero
magnetic field

$\widetilde{X}^{(4)}_{a}=+\underbrace{\frac{1375}{24}}_{57.292}\frac{p_{a}^{2}}{\nu_{aa}^{2}\rho_{a}m_{a}}\nabla^{2}T_{a}.$
(151)

Note that the result is proportional to $1/\nu_{aa}^{2}$ and thus small in a
highly-collisional regime.

###  7.3 Electron species (one ion-electron plasma)

Here we consider a one ion-electron plasma with small temperature differences.
Similar to Braginskii, an exact energy conservation can be imposed by hand,
according to

$\displaystyle
Q_{ie}=\frac{\rho_{i}\nu_{ie}}{m_{i}}\Big{[}3(T_{e}-T_{i})+\frac{3}{10}m_{e}\Big{(}\frac{\widetilde{X}^{(4)}_{i}}{p_{i}}-\frac{\widetilde{X}^{(4)}_{e}}{p_{e}}\Big{)}\Big{]};\qquad
Q_{ei}=-Q_{ie}-(\boldsymbol{u}_{e}-\boldsymbol{u}_{i})\cdot\boldsymbol{R}_{ei}.$
(152)

The electron coefficients (144) become $\hat{S}_{ei(1)}=(26/5)(m_{e}/m_{i})$
and $\hat{S}_{ei(2)}=6(m_{e}/m_{i})^{2}$, and collisional contributions (142)
have a simple form

$\widetilde{Q}^{(4)}_{e}\,{}^{\prime}=-\frac{4}{5}\nu_{ee}\widetilde{X}^{(4)}_{e},$
(153)

determined solely by the electron-electron collisions. A quasi-static solution
of equation (138) then becomes

$\widetilde{X}^{(4)}_{e}=-\,\frac{5\sqrt{2}Z_{i}}{4\nu_{ei}}\Big{[}\nabla\cdot\vec{\boldsymbol{X}}^{(5)}_{e}-20\frac{p_{e}}{\rho_{e}}\nabla\cdot\vec{\boldsymbol{q}}_{e}\Big{]},$
(154)

where we have used $\nu_{ee}=\nu_{ei}/(Z_{i}\sqrt{2})$. The electron heat
fluxes are given by (53) and (61) and are of course determined by both
electron-electron and electron-ion collisions. The full solution thus consists
of thermal and frictional parts
$\widetilde{X}^{(4)}_{e}=\widetilde{X}^{(4)T}_{e}+\widetilde{X}^{(4)u}_{e}$,
where

$\displaystyle\widetilde{X}^{(4)T}_{e}$ $\displaystyle=$
$\displaystyle-\,\frac{5\sqrt{2}Z_{i}}{4\nu_{ei}}\Big{\\{}\nabla\cdot\Big{[}\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{)}\Big{]}$ (155)
$\displaystyle\quad-20\frac{p_{e}}{\rho_{e}}\nabla\cdot\Big{(}-\kappa_{\parallel}^{e}\nabla_{\parallel}T_{e}-\kappa_{\perp}^{e}\nabla_{\perp}T_{e}+\kappa_{\times}^{e}\boldsymbol{\hat{b}}\times\nabla
T_{e}\Big{)}\Big{\\}};$ $\displaystyle\widetilde{X}^{(4)u}_{e}$
$\displaystyle=$
$\displaystyle-\,\frac{5Z_{i}\sqrt{2}}{4\nu_{ei}}\Big{\\{}\nabla\cdot\Big{[}\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{)}\Big{]}$
(156)
$\displaystyle\quad-20\frac{p_{e}}{\rho_{e}}\nabla\cdot\Big{(}\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}\Big{)}\Big{\\}},$

with $\delta\boldsymbol{u}=\boldsymbol{u}_{e}-\boldsymbol{u}_{i}$. It is again
natural to define electron thermal conductivities (of the moment
$\widetilde{X}^{(4)}_{e}$)

$\kappa_{\parallel}^{e(4)}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\kappa_{\parallel}^{e(5)}-20\kappa_{\parallel}^{e}\big{)};\qquad\kappa_{\perp}^{e(4)}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\kappa_{\perp}^{e(5)}-20\kappa_{\perp}^{e}\big{)};\qquad\kappa_{\times}^{e(4)}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\kappa_{\times}^{e(5)}-20\kappa_{\times}^{e}\big{)},$
(157)

together with transport coefficients

$\displaystyle\beta_{0}^{(4)}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\beta_{0}^{(5)}-20\beta_{0}\big{)};\qquad\beta_{1}^{(4)^{\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\beta_{1}^{(5)^{\prime}}-20\beta_{1}^{\prime}\big{)};\qquad\beta_{0}^{(4)^{\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\beta_{0}^{(5)^{\prime}}-20\beta_{0}^{\prime}\big{)};$
$\displaystyle\beta_{1}^{(4)^{\prime\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\beta_{1}^{(5)^{\prime\prime}}-20\beta_{1}^{\prime\prime}\big{)};\qquad\beta_{0}^{(4)^{\prime\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\beta_{0}^{(5)^{\prime\prime}}-20\beta_{0}^{\prime\prime}\big{)};$
(158)
$\displaystyle\gamma_{0}^{(4)}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\gamma_{0}^{(5)}-20\gamma_{0}\big{)};\qquad\gamma_{1}^{(4)^{\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\gamma_{1}^{(5)^{\prime}}-20\gamma_{1}^{\prime}\big{)};\qquad\gamma_{0}^{(4)^{\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\gamma_{0}^{(5)^{\prime}}-20\gamma_{0}^{\prime}\big{)};$
$\displaystyle\gamma_{1}^{(4)^{\prime\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\gamma_{1}^{(5)^{\prime\prime}}-20\gamma_{1}^{\prime\prime}\big{)};\qquad\gamma_{0}^{(4)^{\prime\prime}}=\frac{5\sqrt{2}Z_{i}}{4}\big{(}\gamma_{0}^{(5)^{\prime\prime}}-20\gamma_{0}^{\prime\prime}\big{)}.$
(159)

The thermal and frictional parts then become

$\displaystyle\widetilde{X}^{(4)T}_{e}$ $\displaystyle=$
$\displaystyle-\,\frac{p_{e}}{\nu_{ei}\rho_{e}}\nabla\cdot\Big{(}-\kappa_{\parallel}^{e(4)}\nabla_{\parallel}T_{e}-\kappa_{\perp}^{e(4)}\nabla_{\perp}T_{e}+\kappa_{\times}^{e(4)}\boldsymbol{\hat{b}}\times\nabla
T_{e}\Big{)}$ (160)
$\displaystyle-\,\frac{5\sqrt{2}Z_{i}}{4\nu_{ei}}\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{)}\cdot\nabla\Big{(}\frac{p_{e}}{\rho_{e}}\Big{)};$
$\displaystyle\widetilde{X}^{(4)u}_{e}$ $\displaystyle=$
$\displaystyle-\,\frac{p_{e}}{\nu_{ei}\rho_{e}}\nabla\cdot\Big{(}\beta_{0}^{(4)}p_{e}\delta\boldsymbol{u}_{\parallel}+\frac{\beta_{1}^{(4)^{\prime}}x^{2}+\beta_{0}^{(4)^{\prime}}}{\triangle}p_{e}\delta\boldsymbol{u}_{\perp}-\frac{\beta_{1}^{(4)^{\prime\prime}}x^{3}+\beta_{0}^{(4)^{\prime\prime}}x}{\triangle}p_{e}\boldsymbol{\hat{b}}\times\delta\boldsymbol{u}\Big{)}$
(161)
$\displaystyle-\,\frac{5\sqrt{2}Z_{i}}{4\nu_{ei}}\nabla\cdot\Big{(}\beta_{0}^{(5)}p_{e}\delta\boldsymbol{u}_{\parallel}+\frac{\beta_{1}^{(5)^{\prime}}x^{2}+\beta_{0}^{(5)^{\prime}}}{\triangle}p_{e}\delta\boldsymbol{u}_{\perp}-\frac{\beta_{1}^{(5)^{\prime\prime}}x^{3}+\beta_{0}^{(5)^{\prime\prime}}x}{\triangle}p_{e}\boldsymbol{\hat{b}}\times\delta\boldsymbol{u}\Big{)}\cdot\nabla\Big{(}\frac{p_{e}}{\rho_{e}}\Big{)},$

where the second terms of (160) & (161) are purely non-linear and may be
neglected for simplicity. The thermal conductivities are

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

and the transport coefficients become

$\displaystyle\beta_{0}^{(4)}=\frac{150Z_{i}^{2}\sqrt{2}(16\sqrt{2}+29Z_{i})}{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288};\qquad\beta_{1}^{(4)^{\prime}}=-\,\frac{3\sqrt{2}(548\sqrt{2}+1261Z_{i})}{224};$
$\displaystyle\beta_{0}^{(4)^{\prime}}=\frac{3\sqrt{2}(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(16\sqrt{2}+29Z_{i})}{9800Z_{i}^{2}};$
$\displaystyle\beta_{1}^{(4)^{\prime\prime}}=-\,\frac{15Z_{i}\sqrt{2}}{4};\qquad\beta_{0}^{(4)^{\prime\prime}}=\frac{3\sqrt{2}(3079Z_{i}^{2}+3181Z_{i}\sqrt{2}+1420)}{490Z_{i}};$
(163)
$\displaystyle\gamma_{0}^{(4)}=\frac{250Z_{i}^{2}\sqrt{2}(66\sqrt{2}+229Z_{i})}{217Z_{i}^{2}+604Z_{i}\sqrt{2}+288};\qquad\gamma_{1}^{(4)^{\prime}}=\frac{5\sqrt{2}(4\sqrt{2}-17Z_{i})}{8};$
$\displaystyle\gamma_{0}^{(4)^{\prime}}=\frac{\sqrt{2}(217Z_{i}^{2}+604Z_{i}\sqrt{2}+288)(66\sqrt{2}+229Z_{i})}{1960Z_{i}^{2}};$
$\displaystyle\gamma_{1}^{(4)^{\prime\prime}}=25Z_{i}\sqrt{2};\qquad\gamma_{0}^{(4)^{\prime\prime}}=\frac{\sqrt{2}(176437Z_{i}^{2}+102558Z_{i}\sqrt{2}+30480)}{784Z_{i}};$
(164) $\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};\qquad\delta_{1}=\frac{586601Z_{i}^{2}+330152Z_{i}\sqrt{2}+106016}{78400Z_{i}^{2}},$
(165)

and with numerical values for $Z_{i}=1$

$\displaystyle\beta_{0}^{(4)}$ $\displaystyle=$ $\displaystyle
8.0576;\qquad\beta_{1}^{(4)^{\prime}}=-38.5624;\qquad\beta_{0}^{(4)^{\prime}}=30.3787;\qquad\beta_{1}^{(4)^{\prime\prime}}=-5.3033;\qquad\beta_{0}^{(4)^{\prime\prime}}=77.9054;$
$\displaystyle\gamma_{0}^{(4)}$ $\displaystyle=$ $\displaystyle
83.8471;\qquad\gamma_{1}^{(4)^{\prime}}=-10.0260;\qquad\gamma_{0}^{(4)^{\prime}}=316.1179;\qquad\gamma_{1}^{(4)^{\prime\prime}}=35.3553;\qquad\gamma_{0}^{(4)^{\prime\prime}}=634.8735;$
$\displaystyle\delta_{0}$ $\displaystyle=$ $\displaystyle
3.7702;\qquad\delta_{1}=14.7898.$ (166)

At the semi-linear level the solution becomes

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

and for zero magnetic field

$\displaystyle\widetilde{X}^{(4)}_{e}$ $\displaystyle=$
$\displaystyle\gamma_{0}^{(4)}\frac{p_{e}^{2}}{\nu_{ei}^{2}\rho_{e}m_{e}}\nabla^{2}T_{e}-\beta_{0}^{(4)}\frac{p_{e}^{2}}{\nu_{ei}\rho_{e}}\nabla\cdot\delta\boldsymbol{u}.$
(169)

##  8 Discussion & Conclusions

Here we discuss various topics that we find of importance.

###  8.1 Energy conservation

Collisional integrals were calculated in a semi-linear approximation, where
all quantities such as
$\vec{\boldsymbol{q}}_{a}\cdot(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})$ or
$|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}$ are neglected and considered
small. This approach is typically used for calculations with Landau or
Boltzmann collisional operators, and is for example used in the models of
Burgers (1969) and Schunk (1977). Importantly, an exact energy conservation
$Q_{ab}+Q_{ba}=(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\cdot\boldsymbol{R}_{ab}$
cannot be achieved, because the collisional integrals would have to be
calculated non-linearly. An exact conservation of energy can be achieved only
in two particular cases, the first being a one ion-electron plasma (or a two-
species plasma) where the conservation of energy can be imposed by hand,
according to

$Q_{ie}=3n_{e}\nu_{ei}(T_{e}-T_{i})\frac{m_{e}}{m_{i}};\qquad
Q_{ei}=-Q_{ie}-(\boldsymbol{u}_{e}-\boldsymbol{u}_{i})\cdot\boldsymbol{R}_{ei},$
(170)

which is the choice of Braginskii (1965); see his equation (2.18). Such a
construction cannot be done in general for multi-species plasmas, and
conservation of energy is thus satisfied only approximately.

The second particular case is by neglecting all heat fluxes and stress-
tensors, and considering only a 5-moment model with perturbation $\chi_{a}=0$.
In this specific example of collisions between strict Maxwellians, multi-fluid
calculations can be done for unrestricted drifts (see Burgers (1969), Schunk
(1977), and our Appendix G.3), yielding

$\displaystyle\boldsymbol{R}_{ab}$ $\displaystyle=$
$\displaystyle\rho_{a}\nu_{ab}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})\Phi_{ab};$
(171) $\displaystyle Q_{ab}$ $\displaystyle=$
$\displaystyle\rho_{a}\nu_{ab}\Big{[}3\frac{T_{b}-T_{a}}{m_{a}+m_{b}}\Psi_{ab}+\frac{m_{b}}{m_{a}+m_{b}}|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}\Phi_{ab}\Big{]},$
(172)

where one defines functions

$\displaystyle\Psi_{ab}$ $\displaystyle=$ $\displaystyle
e^{-\epsilon^{2}};\qquad\Phi_{ab}=\Big{(}\frac{3}{4}\sqrt{\pi}\frac{\textrm{erf}(\epsilon)}{\epsilon^{3}}-\frac{3}{2}\frac{e^{-\epsilon^{2}}}{\epsilon^{2}}\Big{)};\qquad\epsilon=\frac{|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|}{\sqrt{v_{\textrm{th}a}^{2}+v_{\textrm{th}b}^{2}}},$
(173)

thermal velocities $v_{\textrm{th}a}^{2}=2T_{a}/m_{a}$, and collisional
frequencies (178). Because $\rho_{a}\nu_{ab}=\rho_{b}\nu_{ba}$ holds, both
momentum and energy are conserved. Collisional exchange rates (171), (172)
represent the “runaway” effect, and the function $\Phi_{ab}$ is directly
related to the Chandrasekhar function; for further details see Appendix G.3
and Figure 6.

For a particular case when differences in drift velocities
$|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|$ become much smaller than thermal
velocities so that $\epsilon\ll 1$, functions $\Phi_{ab}\to 1$ and
$\Psi_{ab}\to 1$ and
$\boldsymbol{R}_{ab}=\rho_{a}\nu_{ab}(\boldsymbol{u}_{b}-\boldsymbol{u}_{a})$.
To correctly account for small $|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}$
contributions in $Q_{ab}$ while keeping the differences in temperatures
unrestricted is achieved by $\Psi_{ab}=1-\epsilon^{2}$, yielding the following
equivalent forms

$\displaystyle Q_{ab}$ $\displaystyle=$
$\displaystyle\rho_{a}\nu_{ab}\Big{[}3\frac{T_{b}-T_{a}}{m_{a}+m_{b}}\Big{(}1-\frac{|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}}{\frac{2T_{a}}{m_{a}}+\frac{2T_{b}}{m_{b}}}\Big{)}+\frac{m_{b}}{m_{b}+m_{a}}|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}\Big{]};$
(174) $\displaystyle Q_{ab}$ $\displaystyle=$
$\displaystyle\rho_{a}\nu_{ab}\Big{[}3\frac{T_{b}-T_{a}}{m_{a}+m_{b}}+\frac{3}{2}\Big{(}\frac{T_{a}m_{b}}{T_{a}m_{b}+T_{b}m_{a}}-\frac{1}{3}\frac{m_{b}}{m_{b}+m_{a}}\Big{)}|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}\Big{]};$
(175)

see also (G32). Energy is still conserved. When additionally the differences
in temperatures are small as well (with respect to their mean temperature),
the frictional part simplifies into

$Q_{ab}=\rho_{a}\nu_{ab}\Big{[}3\frac{T_{b}-T_{a}}{m_{a}+m_{b}}+\frac{m_{b}}{m_{a}+m_{b}}|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}\Big{]}.$
(176)

One can of course neglect the runaway effect from the beginning, and account
for small $|\boldsymbol{u}_{b}-\boldsymbol{u}_{a}|^{2}$ contributions either
through the center-of-mass velocity transformation, as is for example done in
the Appendix of Braginskii (1965), or by using the Rosenbluth potentials; see
Appendix G.1, G.2.

Note that considering the 22-moment model, the fully contracted scalar
perturbations $\widetilde{X}^{(4)}$ modify the energy conservation, according
to

$\displaystyle
Q_{ab}=\frac{\rho_{a}\nu_{ab}}{(m_{a}+m_{b})}\Big{[}3(T_{b}-T_{a})+\hat{P}_{ab(1)}\frac{\rho_{a}}{n_{a}p_{a}}\widetilde{X}^{(4)}_{a}-\hat{P}_{ab(2)}\frac{\rho_{b}}{n_{b}p_{b}}\widetilde{X}^{(4)}_{b}\Big{]};$
$\displaystyle\hat{P}_{ab(1)}=\frac{3T_{a}m_{b}(5T_{b}m_{b}+4T_{b}m_{a}-T_{a}m_{b})}{40(T_{a}m_{b}+T_{b}m_{a})^{2}};\qquad\hat{P}_{ab(2)}=\frac{3T_{b}m_{a}(5T_{a}m_{a}+4T_{a}m_{b}-T_{b}m_{a})}{40(T_{a}m_{b}+T_{b}m_{a})^{2}},$
(177)

and for only two species one can again impose an exact energy conservation by
hand; see e.g. (152).

###  8.2 Collisional frequencies for ion-electron plasma

The Landau collisional operator yields the following collisional frequencies
(see for example Hinton (1983) or our Appendix G.1)

$\nu_{ab}=\tau^{-1}_{ab}=\frac{16}{3}\sqrt{\pi}\frac{n_{b}e^{4}Z_{a}^{2}Z_{b}^{2}\ln\Lambda}{m_{a}^{2}(v_{\textrm{th}a}^{2}+v_{\textrm{th}b}^{2})^{3/2}}\Big{(}1+\frac{m_{a}}{m_{b}}\Big{)},$
(178)

where $v_{\textrm{th}a}^{2}=2T_{a}/m_{a}$, and
$\rho_{a}\nu_{ab}=\rho_{b}\nu_{ba}$ holds. Equivalently, in the form of
Burgers (1969) & Schunk (1977)

$\nu_{ab}=\frac{16}{3}\sqrt{\pi}\Big{(}\frac{\mu_{ab}}{2T_{ab}}\Big{)}^{3/2}\frac{m_{b}n_{b}}{m_{a}+m_{b}}\frac{e^{4}Z_{a}^{2}Z_{b}^{2}\ln\Lambda}{\mu_{ab}^{2}},$
(179)

where the reduced mass $\mu_{ab}$ and reduced temperature $T_{ab}$ are defined
in (14). For a particular case of self-collisions

$\nu_{aa}=\frac{4}{3}\sqrt{\pi}\frac{n_{a}e^{4}Z_{a}^{4}\ln\Lambda}{T_{a}^{3/2}\sqrt{m_{a}}}.$
(180)

For a particular case of $T_{a}=T_{b}=T$,

$\nu_{ab}=\frac{4}{3}\sqrt{2\pi}\frac{n_{b}e^{4}Z_{a}^{2}Z_{b}^{2}\ln\Lambda}{T^{3/2}}\frac{\sqrt{\mu_{ab}}}{m_{a}},$
(181)

which identifies with equation (7.6) of Braginskii (1965) (after one uses
$\nu_{ab}=n_{b}\mu_{ab}\alpha_{ab}^{\prime}/m_{a}$). For a particular case of
a one ion-electron plasma, collisional frequencies simplify into

$\displaystyle\nu_{ii}=\frac{4}{3}\sqrt{\pi}\frac{n_{i}e^{4}Z_{i}^{4}\ln\Lambda}{T^{3/2}_{i}\sqrt{m_{i}}};\qquad\nu_{ie}=\frac{4}{3}\sqrt{2\pi}\frac{n_{e}e^{4}Z_{i}^{2}\ln\Lambda}{T^{3/2}_{e}\sqrt{m_{i}}}\sqrt{\frac{m_{e}}{m_{i}}};$
$\displaystyle\nu_{ee}=\frac{4}{3}\sqrt{\pi}\frac{n_{e}e^{4}\ln\Lambda}{T^{3/2}_{e}\sqrt{m_{e}}};\qquad\nu_{ei}=\frac{4}{3}\sqrt{2\pi}\frac{n_{i}e^{4}Z_{i}^{2}\ln\Lambda}{T^{3/2}_{e}\sqrt{m_{e}}},$
(182)

where one assumes $T_{i}/m_{i}\ll T_{e}/m_{e}$, so the ions cannot be
extremely hot. Obviously, $\nu_{ii}\gg\nu_{ie}$ (by a factor of
$\sqrt{m_{i}/m_{e}}$ for equal temperatures and $Z_{i}=1$), however
$\nu_{ee}\sim\nu_{ei}$, with exact relation $\nu_{ei}=Z_{i}\sqrt{2}\nu_{ee}$
after one uses $n_{e}=Z_{i}n_{i}$. The relation
$\rho_{i}\nu_{ie}=\rho_{e}\nu_{ei}$ holds exactly in (182). Note the important
difference that while $\nu_{ei}$ contains a factor of $\sqrt{2}$, the
$\nu_{ii}$ does not. Thus, comparing Braginskii (1965) expressions (2.5i) and
(2.5e) with definitions (182), Braginskii clearly uses

$\tau_{i}=\tau_{ii};\qquad\tau_{e}=\tau_{ei},$ (183)

which also agrees with his definition (7.6), equivalent to our (181).

However, very often when considering ion-electron plasma, a different
definition of $\nu_{ab}$ is used without the reduced mass, in the following
form

$\displaystyle m_{a}\ll m_{b}:\qquad\nu_{ab}=\tau_{ab}^{-1}$ $\displaystyle=$
$\displaystyle\frac{16}{3}\sqrt{\pi}\frac{n_{b}e^{4}Z_{a}^{2}Z_{b}^{2}\ln\Lambda}{m_{a}^{2}v_{\textrm{th}a}^{3}}=\frac{4}{3}\sqrt{2\pi}\frac{n_{b}e^{4}Z_{a}^{2}Z_{b}^{2}\ln\Lambda}{T_{a}^{3/2}\sqrt{m_{a}}},$
(184)

which for example agrees with the Appendix of Helander & Sigmar (2002) (page
277, after using cgs units $\epsilon_{0}\to 1/(4\pi)$). We have added the
$m_{a}\ll m_{b}$ designation even though it is not present in Helander &
Sigmar (2002), because obviously it is the only way how to obtain (184) from
the general (178). Importantly, $\rho_{a}\nu_{ab}\neq\rho_{b}\nu_{ba}$, and if
one would use (184) to calculate $\nu_{ie}$, the result would be erroneous.
Instead, the $\nu_{ie}$ must be calculated from $\nu_{ei}$ so that the
momentum is conserved. Technically, (184) should not be used for self-
collisions either. Nevertheless, using (184) yields the following collisional
frequencies

$\displaystyle\nu_{ii}=\frac{4}{3}\sqrt{2\pi}\frac{n_{i}e^{4}Z_{i}^{4}\ln\Lambda}{T^{3/2}_{i}\sqrt{m_{i}}};\qquad\nu_{ie}=\frac{m_{e}n_{e}}{m_{i}n_{i}}\nu_{ei}=\frac{4}{3}\sqrt{2\pi}\frac{n_{e}e^{4}Z_{i}^{2}\ln\Lambda}{T^{3/2}_{e}\sqrt{m_{i}}}\sqrt{\frac{m_{e}}{m_{i}}};$
$\displaystyle\nu_{ee}=\frac{4}{3}\sqrt{2\pi}\frac{n_{e}e^{4}\ln\Lambda}{T^{3/2}_{e}\sqrt{m_{e}}};\qquad\nu_{ei}=\frac{4}{3}\sqrt{2\pi}\frac{n_{i}e^{4}Z_{i}^{2}\ln\Lambda}{T^{3/2}_{e}\sqrt{m_{e}}}.$
(185)

Now $\nu_{ii}$ contains a factor of $\sqrt{2}$, leading to an interpretation
that Braginskii uses

$\tau_{i}=\sqrt{2}\tau_{ii};\qquad\tau_{e}=\tau_{ei}.$ (186)

Also, now for $Z_{i}=1$ relation $\nu_{ee}=\nu_{ei}$ holds. These definitions
of collisional frequencies are used in a majority of the modern plasma
literature, where one argues that it seems unnatural to introduce asymmetry
between $\nu_{ii}$ and $\nu_{ei}$ (see e.g. Part 1 of Balescu (1988), p.192,
p.274). Obviously, for multi-species plasmas collisional frequencies (178)
have to be used, and we thus find it much more natural to use the original
Braginskii (1965) definitions (182), (183) also for an ion-electron plasma. Of
course, for the Landau operator both approaches yield the same results because
the collisional integrals are properly calculated. However, a difference
arises for the phenomenological operators such as the BGK or the Dougherty
(Lenard-Bernstein) operators, where for example one needs to add
$\nu_{ee}+\nu_{ei}$. Calculating this addition according to (185) would be
incorrect, and one has to use (182) instead. Comparison of Braginskii
viscosities and heat conductivities with the BGK operator can be found in
Appendix E.3.

###  8.3 Fluid hierarchy

Even though we do not calculate the collisional integrals for general n-th
order moments, we find it useful to discuss the fluid hierarchy and formulate
it for a general collisional operator $C(f_{a})$. One defines heat flux
vectors, stress-tensors and fully contracted moments according to

$\displaystyle\vec{\boldsymbol{X}}^{(2n+1)}_{a}$ $\displaystyle=$
$\displaystyle
m_{a}\int\boldsymbol{c}_{a}|\boldsymbol{c}_{a}|^{2n}f_{a}d^{3}v;$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2n)}_{a}$ $\displaystyle=$
$\displaystyle
m_{a}\int\big{(}\boldsymbol{c}_{a}\boldsymbol{c}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}|\boldsymbol{c}_{a}|^{2}\big{)}|\boldsymbol{c}_{a}|^{2n-2}f_{a}d^{3}v;$
$\displaystyle X^{(2n)}_{a}$ $\displaystyle=$ $\displaystyle
m_{a}\int|\boldsymbol{c}_{a}|^{2n}f_{a}d^{3}v=(2n+1)!!\frac{p^{n}_{a}}{\rho^{n-1}_{a}}+\widetilde{X}^{(2n)}_{a},$
(187)

together with collisional contributions

$\displaystyle\vec{\boldsymbol{\mathcal{Q}}}^{(2n+1)}_{a}$ $\displaystyle=$
$\displaystyle
m_{a}\int|\boldsymbol{c}_{a}|^{2n}\boldsymbol{c}_{a}C(f_{a})d^{3}v;$
$\displaystyle\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(2n)}_{a}$
$\displaystyle=$ $\displaystyle
m_{a}\int|\boldsymbol{c}_{a}|^{2n-2}\boldsymbol{c}_{a}\boldsymbol{c}_{a}C(f_{a})d^{3}v;$
$\displaystyle Q^{(2n)}_{a}$ $\displaystyle=$ $\displaystyle
m_{a}\int|\boldsymbol{c}_{a}|^{2n}C(f_{a})d^{3}v;\qquad
Q_{a}=\frac{m_{a}}{2}\int|\boldsymbol{c}_{a}|^{2}C(f_{a})d^{3}v,$ (188)

where to prevent incompatibility with the previous notation, for vectors and
matrices we use $\boldsymbol{\mathcal{Q}}$ (mathcal Q) instead of
$\boldsymbol{Q}$. The new notation fixes the problem that for example
$\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}$ was used for the right hand side
of evolution equation of the heat flux $\vec{\boldsymbol{q}}_{a}$, and not for
$\vec{\boldsymbol{X}}^{(3)}_{a}$. It also clarifies that in the vector
notation the matrix
$\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(2n)}_{a}=\textrm{Tr}\textrm{Tr}\ldots\textrm{Tr}\bar{\bar{\boldsymbol{Q}}}_{a}^{(2n)}$.
Note that $Q^{(2)}_{a}=2Q_{a}$; $X^{(2)}_{a}=3p_{a}$;
$\widetilde{X}^{(2)}_{a}=0$;
$\vec{\boldsymbol{X}}^{(3)}_{a}=2\vec{\boldsymbol{q}}_{a}$ and
$\vec{\boldsymbol{X}}^{(1)}_{a}=0$.

Fully non-linear evolution equations are given in Appendix D; see (D13)-(D15).
In the semi-linear approximation, these simplify into evolution equations for
vectors valid for $n\geq 1$

$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(2n+1)}_{a}+\frac{1}{3}\nabla\widetilde{X}^{(2n+2)}_{a}+\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2n+2)}_{a}-\frac{(2n+3)!!}{3}\frac{p_{a}^{n}}{\rho_{a}^{n}}\nabla\cdot\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(2n+1)}_{a}$
$\displaystyle+(2n+3)!!\frac{(n)}{3}\frac{p_{a}^{n}}{\rho_{a}^{n-1}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}=\vec{\boldsymbol{\mathcal{Q}}}^{(2n+1)}_{a}\,{}^{\prime}=\vec{\boldsymbol{\mathcal{Q}}}^{(2n+1)}_{a}-\frac{(2n+3)!!}{3}\frac{p_{a}^{n}}{\rho_{a}^{n}}\boldsymbol{R}_{a},$
(189)

stress-tensors valid for $n\geq 1$

$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2n)}_{a}+\frac{1}{5}\Big{[}\big{(}\nabla\vec{\boldsymbol{X}}^{(2n+1)}_{a}\big{)}^{S}-\frac{2}{3}\bar{\bar{\boldsymbol{I}}}\nabla\cdot\vec{\boldsymbol{X}}^{(2n+1)}_{a}\Big{]}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2n)}_{a}\big{)}^{S}$
$\displaystyle+\frac{(2n+3)!!}{15}\frac{p^{n}_{a}}{\rho_{a}^{n-1}}\bar{\bar{\boldsymbol{W}}}_{a}=\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(2n)}_{a}\,{}^{\prime}=\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(2n)}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}Q^{(2n)}_{a},$
(190)

and scalar perturbations valid for $n\geq 2$

$\displaystyle\frac{d_{a}}{dt}\widetilde{X}^{(2n)}_{a}+\nabla\cdot\vec{\boldsymbol{X}}^{(2n+1)}_{a}-(2n+1)!!\frac{(2n)}{3}\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}^{n-1}\nabla\cdot\vec{\boldsymbol{q}}_{a}$
$\displaystyle=\widetilde{Q}^{(2n)}_{a}\,{}^{\prime}=Q^{(2n)}_{a}-(2n+1)!!\frac{(2n)}{3}\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}^{n-1}Q_{a},$
(191)

where (n) without species index should not be confused with the number
density. Equation (191) is also valid for $n=1$, but it is identically zero.
In comparison to previous notation
$\vec{\boldsymbol{\mathcal{Q}}}^{(3)}_{a}\,{}^{\prime}=2\vec{\boldsymbol{Q}}^{(3)}_{a}\,{}^{\prime}$;
$\vec{\boldsymbol{\mathcal{Q}}}^{(5)}_{a}\,{}^{\prime}=\vec{\boldsymbol{Q}}^{(5)}_{a}\,{}^{\prime}$;
$\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(2)}_{a}\,{}^{\prime}=\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a}\,{}^{\prime}$;
and
$\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(4)}_{a}\,{}^{\prime}=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}\,{}^{\prime}$.

###  8.4 Reducible & irreducible Hermite polynomials

The _irreducible_ Hermite polynomials $H(\tilde{\boldsymbol{c}})$ (notation
without tilde) are usually defined through Laguerre-Sonine polynomials
$L(\tilde{\boldsymbol{c}})$ (see for example equation (G1.4.4), page 326 of
Balescu (1988))

$\displaystyle H^{(2n)}(\tilde{c})$ $\displaystyle=$ $\displaystyle
L_{n}^{(1/2)}(\frac{\tilde{c}^{2}}{2});\qquad
H^{(2n+1)}_{i}(\tilde{\boldsymbol{c}})=\sqrt{\frac{3}{2}}\tilde{c}_{i}L_{n}^{(3/2)}(\frac{\tilde{c}^{2}}{2});$
$\displaystyle H^{(2n)}_{ij}(\tilde{\boldsymbol{c}})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{15}{8}}(\tilde{c}_{i}\tilde{c}_{j}-\frac{\tilde{c}^{2}}{3}\delta_{ij})L_{n-1}^{(5/2)}(\frac{\tilde{c}^{2}}{2}),$
(192)

where we use tilde for the normalized fluctuating velocity
$\tilde{\boldsymbol{c}}=\sqrt{m_{a}/T_{a}}\boldsymbol{c}_{a}$ with species
indices dropped. In our calculations we find it more natural to use the
_reducible_ Hermite polynomials $\tilde{H}(\tilde{\boldsymbol{c}})$ (notation
with tilde) of Grad defined according to

$\tilde{H}^{(m)}_{r_{1}r_{2}\ldots
r_{m}}(\tilde{\boldsymbol{c}})=(-1)^{m}e^{\frac{\tilde{c}^{2}}{2}}\frac{\partial}{\partial\tilde{c}_{r_{1}}}\frac{\partial}{\partial\tilde{c}_{r_{2}}}\cdots\frac{\partial}{\partial\tilde{c}_{r_{m}}}e^{-\frac{\tilde{c}^{2}}{2}}.$
(193)

By applying a sufficient number of contractions then yields definitions of
fully contracted scalars, vectors and matrices

$\tilde{H}^{(2n)}=\tilde{H}^{(2n)}_{r_{1}r_{1}\ldots
r_{n}r_{n}};\qquad\tilde{H}^{(2n+1)}_{i}=\tilde{H}^{(2n+1)}_{ir_{1}r_{1}\ldots
r_{n}r_{n}};\qquad\tilde{H}^{(2n)}_{ij}=\tilde{H}^{(2n)}_{ijr_{1}r_{1}\ldots
r_{n-1}r_{n-1}},$ (194)

together with conveniently defined traceless matrices (notation with hat)

$\hat{H}_{ij}^{(2n)}=\tilde{H}_{ij}^{(2n)}-\frac{1}{3}\delta_{ij}\tilde{H}^{(2n)}.$
(195)

The relation between irreducible and reducible Hermite polynomials then can be
shown to be

$\displaystyle H^{(2n)}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{1}{2^{n}n!(2n+1)!!}\Big{)}^{1/2}\tilde{H}^{(2n)};\qquad
H^{(2n+1)}_{i}=\Big{(}\frac{3}{2^{n}n!(2n+3)!!}\Big{)}^{1/2}\tilde{H}_{i}^{(2n+1)};$
$\displaystyle H^{(2n)}_{ij}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{15}{2^{n}(n-1)!(2n+3)!!}\Big{)}^{1/2}\hat{H}_{ij}^{(2n)},$
(196)

and both approaches use essentially the same polynomials, the only difference
is the location of normalization factors. The reducible Hermite polynomials
are used to define Hermite moments

$\displaystyle\tilde{h}^{(2n)}=\frac{1}{n_{a}}\int
f_{a}\tilde{H}^{(2n)}d^{3}c;\qquad\tilde{h}^{(2n+1)}_{i}=\frac{1}{n_{a}}\int
f_{a}\tilde{H}^{(2n+1)}_{i}d^{3}c;\qquad\hat{h}^{(2n)}_{ij}=\frac{1}{n_{a}}\int
f_{a}\hat{H}^{(2n)}_{ij}d^{3}c,$ (197)

and analogously for the irreducible ones. Note that the scalar
$\tilde{h}^{(2)}=0$, and we thus often use
$\hat{h}^{(2)}_{ij}=\tilde{h}^{(2)}_{ij}=(1/n_{a})\int
f_{a}\tilde{H}^{(2)}_{ij}d^{3}c$. Finally, by using orthogonality relations
one obtains perturbation $\chi_{a}$ of the distribution function
$f_{a}=f_{a}^{(0)}(1+\chi_{a})$ around Maxwellian $f_{a}^{(0)}$, in the
following form

$\displaystyle\chi_{a}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\Big{[}\frac{15}{2^{n}(n-1)!(2n+3)!!}\hat{h}^{(2n)}_{ij}\hat{H}^{(2n)}_{ij}+\frac{1}{2^{n}n!(2n+1)!!}\tilde{h}^{(2n)}\tilde{H}^{(2n)}$
(198)
$\displaystyle+\frac{3}{2^{n}n!(2n+3)!!}\tilde{h}^{(2n+1)}_{i}\tilde{H}^{(2n+1)}_{i}\Big{]};$
$\displaystyle\chi_{a}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\Big{[}{h}^{(2n)}_{ij}{H}^{(2n)}_{ij}+{h}^{(2n)}{H}^{(2n)}+{h}^{(2n+1)}_{i}{H}^{(2n+1)}_{i}\Big{]},$
(199)

and the two approaches are equivalent. Alternatively, because
$\hat{h}^{(2n)}_{ij}$ are traceless, it is possible to use
$\hat{h}^{(2n)}_{ij}\hat{H}^{(2n)}_{ij}=\hat{h}^{(2n)}_{ij}\tilde{H}^{(2n)}_{ij}$.
Note that $2^{n}n!(2n+1)!!=(2n+1)!$. The 13-moment model of Burgers-Schunk is
obtained by $N=1$. Prescribing $N=2$ yields the 22-moment model

$\chi_{a}=\frac{1}{2}\hat{h}_{ij}^{(2)}\hat{H}_{ij}^{(2)}+\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)}+\frac{1}{28}\hat{h}_{ij}^{(4)}\hat{H}_{ij}^{(4)}+\frac{1}{120}\tilde{h}^{(4)}\tilde{H}^{(4)}+\frac{1}{280}\tilde{h}_{i}^{(5)}\tilde{H}_{i}^{(5)},$
(200)

with Hermite polynomials

$\displaystyle\tilde{H}^{(3)}_{i}$ $\displaystyle=$
$\displaystyle\tilde{c}_{i}(\tilde{c}^{2}-5);\qquad\tilde{H}^{(5)}_{i}=\tilde{c}_{i}(\tilde{c}^{4}-14\tilde{c}^{2}+35);$
$\displaystyle\hat{H}^{(2)}_{ij}$ $\displaystyle=$
$\displaystyle\big{(}\tilde{c}_{i}\tilde{c}_{j}-\frac{\delta_{ij}}{3}\tilde{c}^{2}\big{)};\qquad\hat{H}^{(4)}_{ij}=\big{(}\tilde{c}_{i}\tilde{c}_{j}-\frac{\delta_{ij}}{3}\tilde{c}^{2}\big{)}(\tilde{c}^{2}-7);\qquad\tilde{H}^{(4)}=\tilde{c}^{4}-10\tilde{c}^{2}+15,$
(201)

and neglecting $\tilde{h}^{(4)}=0$ (meaning $\widetilde{X}^{(4)}=0$) yields
the 21-moment model.

Transformation from Hermite to fluid moments is done according to

$\displaystyle\tilde{\boldsymbol{h}}^{(3)}_{a}=\frac{2}{p_{a}}\sqrt{\frac{m_{a}}{T_{a}}}\vec{\boldsymbol{q}}_{a};\qquad\tilde{\boldsymbol{h}}^{(5)}_{a}=\frac{1}{p_{a}}\sqrt{\frac{m_{a}}{T_{a}}}\Big{(}\frac{m_{a}}{T_{a}}\vec{\boldsymbol{X}}^{(5)}_{a}-28\vec{\boldsymbol{q}}_{a}\Big{)};$
$\displaystyle\hat{\bar{\bar{\boldsymbol{h}}}}^{(2)}_{a}=\tilde{\bar{\bar{\boldsymbol{h}}}}^{(2)}_{a}=\frac{1}{p_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a};\qquad\hat{\bar{\bar{\boldsymbol{h}}}}^{(4)}_{a}=\frac{\rho_{a}}{p_{a}^{2}}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}-\frac{7}{p_{a}}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a};\qquad\tilde{h}^{(4)}_{a}=\frac{\rho_{a}}{p_{a}^{2}}\widetilde{X}^{(4)}_{a}.$
(202)

Various models are summarized in Tables 8.1 & 8.2. In Table 8.1 the
perturbation $\chi_{a}$ is given in reducible Hermite moments, and in Table
8.2 the perturbation is given in fluid moments.

Model name | Corresponding perturbation of $f_{a}=f_{a}^{(0)}(1+\chi_{a})$ in Hermite moments
---|---
5-moment; |  $\chi_{a}=0$;
8-moment; | $\chi_{a}=\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)};$
10-moment; | $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)};$
13-moment; | $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)}+\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)};$
20-moment; | $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)}+\frac{1}{6}\tilde{h}_{ijk}^{(3)}\tilde{H}_{ijk}^{(3)};$
21-moment; |  $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)}+\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)}+\frac{1}{28}\hat{h}_{ij}^{(4)}\tilde{H}_{ij}^{(4)}+\frac{1}{280}\tilde{h}_{i}^{(5)}\tilde{H}_{i}^{(5)}$;
22-moment; | $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)}+\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)}+\frac{1}{28}\hat{h}_{ij}^{(4)}\tilde{H}_{ij}^{(4)}+\frac{1}{120}\tilde{h}^{(4)}\tilde{H}^{(4)}+\frac{1}{280}\tilde{h}_{i}^{(5)}\tilde{H}_{i}^{(5)};$
9-moment; | $\chi_{a}=\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)}+\frac{1}{120}\tilde{h}^{(4)}\tilde{H}^{(4)};$
11-moment; | $\chi_{a}=\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)}+\frac{1}{280}\tilde{h}_{i}^{(5)}\tilde{H}_{i}^{(5)};$
12-moment; | $\chi_{a}=\frac{1}{10}\tilde{h}_{i}^{(3)}\tilde{H}_{i}^{(3)}+\frac{1}{120}\tilde{h}^{(4)}\tilde{H}^{(4)}+\frac{1}{280}\tilde{h}_{i}^{(5)}\tilde{H}_{i}^{(5)};$
15-moment; | $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)}+\frac{1}{28}\hat{h}_{ij}^{(4)}\tilde{H}_{ij}^{(4)};$
16-moment; | $\chi_{a}=\frac{1}{2}\tilde{h}_{ij}^{(2)}\tilde{H}_{ij}^{(2)}+\frac{1}{28}\hat{h}_{ij}^{(4)}\tilde{H}_{ij}^{(4)}+\frac{1}{120}\tilde{h}^{(4)}\tilde{H}^{(4)};$
Table 8.1: Summary of various models with the perturbation $\chi_{a}$ given in reducible Hermite moments. Species indices “a” are dropped. The upper half of the table contains “major” models, and the lower half contains other possibilities. Note that the 16-moment model should not be confused with the anisotropic (bi-Maxwellian based) 16-moment model described in Section 8.9. Model name | Corresponding perturbation of $f_{a}=f_{a}^{(0)}(1+\chi_{a})$ in fluid moments
---|---
5-moment; | $\chi_{a}=0;$
8-moment; | $\chi_{a}=-\,\frac{m_{a}}{p_{a}T_{a}}\big{(}\vec{\boldsymbol{q}}_{a}\cdot\boldsymbol{c}_{a}\big{)}\Big{(}1-\frac{m_{a}}{5T_{a}}c_{a}^{2}\Big{)};$
10-moment; | $\chi_{a}=\frac{m_{a}}{2p_{a}T_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\boldsymbol{c}_{a}\boldsymbol{c}_{a}\big{)};$
13-moment; | $\chi_{a}=\frac{m_{a}}{2p_{a}T_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\boldsymbol{c}_{a}\boldsymbol{c}_{a}\big{)}-\frac{m_{a}}{p_{a}T_{a}}\big{(}\vec{\boldsymbol{q}}_{a}\cdot\boldsymbol{c}_{a}\big{)}\Big{(}1-\frac{m_{a}}{5T_{a}}c_{a}^{2}\Big{)};$
20-moment; | $\chi_{a}=\frac{m_{a}}{2p_{a}T_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\boldsymbol{c}_{a}\boldsymbol{c}_{a}\big{)}+\frac{m_{a}^{2}}{6p_{a}T_{a}^{2}}\big{(}\boldsymbol{c}_{a}\cdot\bar{\bar{\boldsymbol{q}}}_{a}:\boldsymbol{c}_{a}\boldsymbol{c}_{a}\big{)}-\frac{m_{a}}{p_{a}T_{a}}\big{(}\vec{\boldsymbol{q}}_{a}\cdot\boldsymbol{c}_{a}\big{)};$
21-moment; | $\chi_{a}=\frac{1}{2p_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\tilde{\boldsymbol{c}}_{a}\tilde{\boldsymbol{c}}_{a}\big{)}+\frac{1}{28}\Big{[}\frac{\rho_{a}}{p_{a}^{2}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}:\tilde{\boldsymbol{c}}_{a}\tilde{\boldsymbol{c}}_{a}\big{)}-\frac{7}{p_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\tilde{\boldsymbol{c}}_{a}\tilde{\boldsymbol{c}}_{a}\big{)}\Big{]}(\tilde{c}_{a}^{2}-7)$
|
$\qquad+\frac{1}{5p_{a}}\sqrt{\frac{m_{a}}{T_{a}}}(\vec{\boldsymbol{q}}_{a}\cdot\tilde{\boldsymbol{c}}_{a})(\tilde{c}_{a}^{2}-5)+\frac{1}{280p_{a}}\sqrt{\frac{m_{a}}{T_{a}}}\Big{[}\frac{\rho_{a}}{p_{a}}(\vec{\boldsymbol{X}}^{(5)}_{a}\cdot\tilde{\boldsymbol{c}}_{a})-28(\vec{\boldsymbol{q}}_{a}\cdot\tilde{\boldsymbol{c}}_{a})\Big{]}(\tilde{c}_{a}^{4}-14\tilde{c}_{a}^{2}+35);$
22-moment; | $\chi_{a}=\frac{1}{2p_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\tilde{\boldsymbol{c}}_{a}\tilde{\boldsymbol{c}}_{a}\big{)}+\frac{1}{28}\Big{[}\frac{\rho_{a}}{p_{a}^{2}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{a}:\tilde{\boldsymbol{c}}_{a}\tilde{\boldsymbol{c}}_{a}\big{)}-\frac{7}{p_{a}}\big{(}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}:\tilde{\boldsymbol{c}}_{a}\tilde{\boldsymbol{c}}_{a}\big{)}\Big{]}(\tilde{c}_{a}^{2}-7)$
|
$\qquad+\frac{1}{5p_{a}}\sqrt{\frac{m_{a}}{T_{a}}}(\vec{\boldsymbol{q}}_{a}\cdot\tilde{\boldsymbol{c}}_{a})(\tilde{c}_{a}^{2}-5)+\frac{1}{280p_{a}}\sqrt{\frac{m_{a}}{T_{a}}}\Big{[}\frac{\rho_{a}}{p_{a}}(\vec{\boldsymbol{X}}^{(5)}_{a}\cdot\tilde{\boldsymbol{c}}_{a})-28(\vec{\boldsymbol{q}}_{a}\cdot\tilde{\boldsymbol{c}}_{a})\Big{]}(\tilde{c}_{a}^{4}-14\tilde{c}_{a}^{2}+35)$
|
$\qquad+\frac{1}{120}\frac{\rho_{a}}{p_{a}^{2}}\widetilde{X}^{(4)}_{a}(\tilde{c}_{a}^{4}-10\tilde{c}^{2}+15);$
Table 8.2: Summary of various models with the perturbation $\chi_{a}$ given in
fluid moments. Results for the 21 & 22-moment models are written with
normalized $\tilde{\boldsymbol{c}}_{a}=\sqrt{m_{a}/T_{a}}\boldsymbol{c}_{a}$.

###  8.5 Rosenbluth potentials (22-moment model)

Here we summarize the Rosenbluth potentials, defined according to

$\displaystyle H_{b}(\boldsymbol{v})$ $\displaystyle=$
$\displaystyle\int\frac{f_{b}(\boldsymbol{v}^{\prime})}{|\boldsymbol{v}^{\prime}-\boldsymbol{v}|}d^{3}v^{\prime};\qquad
G_{b}(\boldsymbol{v})=\int|\boldsymbol{v}^{\prime}-\boldsymbol{v}|f_{b}(\boldsymbol{v}^{\prime})d^{3}v^{\prime},$
(203)

where the first potential should not be confused with the irreducible Hermite
polynomials. For the 22-moment model the fully non-linear results read

$\displaystyle H_{b}(\boldsymbol{v})$ $\displaystyle=$ $\displaystyle
n_{b}\sqrt{\frac{m_{b}}{T_{b}}}\Big{\\{}\frac{1}{\tilde{y}}\textrm{erf}\Big{(}\frac{\tilde{y}}{\sqrt{2}}\Big{)}-\sqrt{\frac{2}{\pi}}\frac{e^{-\tilde{y}^{2}/2}}{10}\Big{(}\tilde{\boldsymbol{y}}\cdot\tilde{\boldsymbol{h}}^{(3)}_{b}+(\tilde{y}^{2}-5)\frac{\tilde{\boldsymbol{y}}\cdot\tilde{\boldsymbol{h}}^{(5)}_{b}}{28}\Big{)}$
(204)
$\displaystyle+\frac{1}{2}(\hat{\bar{\bar{\boldsymbol{h}}}}^{(2)}_{b}:\tilde{\boldsymbol{y}}\tilde{\boldsymbol{y}})\Big{[}\frac{3}{\tilde{y}^{5}}\textrm{erf}\Big{(}\frac{\tilde{y}}{\sqrt{2}}\Big{)}-\sqrt{\frac{2}{\pi}}\big{(}\frac{1}{\tilde{y}^{2}}+\frac{3}{\tilde{y}^{4}}\big{)}e^{-\tilde{y}^{2}/2}\Big{]}$
$\displaystyle\qquad-\frac{1}{28}(\hat{\bar{\bar{\boldsymbol{h}}}}^{(4)}_{b}:\tilde{\boldsymbol{y}}\tilde{\boldsymbol{y}})\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^{2}/2}+\frac{1}{120}\tilde{h}^{(4)}_{b}(3-\tilde{y}^{2})\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^{2}/2}\Big{\\}};$
$\displaystyle G_{b}(\boldsymbol{v})$ $\displaystyle=$ $\displaystyle
n_{b}\sqrt{\frac{T_{b}}{m_{b}}}\Big{\\{}\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^{2}/2}+\big{(}\tilde{y}+\frac{1}{\tilde{y}}\big{)}\textrm{erf}\big{(}\frac{\tilde{y}}{\sqrt{2}}\big{)}$
(205)
$\displaystyle+\Big{(}\frac{\textrm{erf}(\tilde{y}/\sqrt{2})}{5\tilde{y}^{3}}-\sqrt{\frac{2}{\pi}}\frac{e^{-\tilde{y}^{2}/2}}{5\tilde{y}^{2}}\Big{)}\tilde{\boldsymbol{y}}\cdot\tilde{\boldsymbol{h}}^{(3)}_{b}-\sqrt{\frac{2}{\pi}}\frac{e^{-\tilde{y}^{2}/2}}{140}\tilde{\boldsymbol{y}}\cdot\tilde{\boldsymbol{h}}^{(5)}_{b}$
$\displaystyle\qquad-\frac{1}{2}(\hat{\bar{\bar{\boldsymbol{h}}}}^{(2)}_{b}:\tilde{\boldsymbol{y}}\tilde{\boldsymbol{y}})\Big{[}\frac{3}{\tilde{y}^{4}}\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^{2}/2}+\big{(}\frac{1}{\tilde{y}^{3}}-\frac{3}{\tilde{y}^{5}}\big{)}\textrm{erf}\Big{(}\frac{\tilde{y}}{\sqrt{2}}\Big{)}\Big{]}$
$\displaystyle\qquad-\frac{1}{14}(\hat{\bar{\bar{\boldsymbol{h}}}}^{(4)}_{b}:\tilde{\boldsymbol{y}}\tilde{\boldsymbol{y}})\Big{[}\sqrt{\frac{2}{\pi}}\big{(}\frac{1}{\tilde{y}^{2}}+\frac{3}{\tilde{y}^{4}}\big{)}e^{-\tilde{y}^{2}/2}-\frac{3}{\tilde{y}^{5}}\textrm{erf}\Big{(}\frac{\tilde{y}}{\sqrt{2}}\Big{)}\Big{]}-\frac{1}{60}\tilde{h}^{(4)}_{b}\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^{2}/2}\Big{\\}},$

where we use the variable

$\tilde{\boldsymbol{y}}=\sqrt{\frac{m_{b}}{T_{b}}}(\boldsymbol{v}-\boldsymbol{u}_{b}).$
(206)

These Rosenbluth potentials are used to calculate the dynamical friction
vector $\boldsymbol{A}_{ab}$ and the diffusion tensor
$\bar{\bar{\boldsymbol{D}}}_{ab}$, which then form the Landau collisional
operator, according to

$\displaystyle\boldsymbol{A}_{ab}(\boldsymbol{v})=2\frac{c_{ab}}{m_{a}^{2}}\big{(}1+\frac{m_{a}}{m_{b}}\big{)}\frac{\partial
H_{b}(\boldsymbol{v})}{\partial\boldsymbol{v}};\qquad\bar{\bar{\boldsymbol{D}}}_{ab}(\boldsymbol{v})=2\frac{c_{ab}}{m_{a}^{2}}\frac{\partial^{2}G_{b}(\boldsymbol{v})}{\partial\boldsymbol{v}\partial\boldsymbol{v}};\qquad
c_{ab}=2\pi e^{4}Z_{a}^{2}Z_{b}^{2}\ln\Lambda;$ $\displaystyle
C_{ab}(f_{a},f_{b})=-\,\frac{\partial}{\partial\boldsymbol{v}}\cdot\Big{[}\boldsymbol{A}_{ab}f_{a}-\frac{1}{2}\frac{\partial}{\partial\boldsymbol{v}}\cdot(\bar{\bar{\boldsymbol{D}}}_{ab}f_{a})\Big{]}.$
(207)

The dynamical friction vectors and diffusion tensors can be found in the
Appendix; see equations (K15)-(K16); (L13)-(L14) and (M4)-(M5). For clarity,
we split the calculations into heat fluxes (Appendix K), viscosities (Appendix
L) and scalar perturbations (Appendix M). These results are fully non-linear
and could be potentially useful to construct more sophisticated models that
could capture collisional effects beyond the semi-linear approximation, or
perhaps to explore the runaway effect numerically. All the equations can be
transformed from Hermite moments to fluid moments by (202).

###  8.6 Hermite closures

The general hierarchy of evolution equations (189)-(191) needs to be closed
with appropriate closures at the last retained fluid moment. A correct form of
a fluid closure is obtained in the Hermite space, by cutting the perturbation
$\chi_{a}$ given by (198) at an appropriate $N$. For example, the 22-moment
model is obtained with Hermite closures $\tilde{h}^{(6)}_{a}=0$ and
$\tilde{h}^{a(6)}_{ij}=0$, which translate into fluid closures (134) & (13).

It is useful to summarize closures for higher-order moments, with details
given in Appendix B. It can be shown that for vectors and scalars, fluid
closures derived from Hermite closures read

$\displaystyle\vec{\boldsymbol{X}}^{(2n+1)}_{a}$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{n-1}(-1)^{m+n+1}\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}^{n-m}\frac{n!}{m!(n-m)!}\frac{(2n+3)!!}{(2m+3)!!}\vec{\boldsymbol{X}}^{(2m+1)}_{a};$
$\displaystyle\widetilde{X}^{(2n)}_{a}$ $\displaystyle=$
$\displaystyle\sum_{m=2}^{n-1}(-1)^{m+n+1}\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}^{n-m}\frac{n!}{m!(n-m)!}\frac{(2n+1)!!}{(2m+1)!!}\widetilde{X}^{(2m)}_{a},$
(208)

together with closures for stress-tensors

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2n)}$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{n-2}(-1)^{m+n}\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}^{n-m-1}\frac{(n-1)!}{m!(n-m-1)!}\frac{(2n+3)!!}{(2m+5)!!}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2m+2)},$
(209)

where the result is zero if the upper summation index is less than the lower
summation index, yielding closures $\vec{\boldsymbol{X}}^{(3)}_{a}=0$;
$\widetilde{X}^{(4)}_{a}=0$ and $\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}=0$.
The closures are summarized bellow in Tables 8.3 and 8.4.

Hermite closures | Fluid closures
---|---
$\tilde{h}_{i}^{(3)}=0$; |  $X^{(3)}_{i}=0$;
$\tilde{h}^{(4)}=0$; |  $\widetilde{X}^{(4)}=0$;
$\tilde{h}_{i}^{(5)}=0$; |  $X^{(5)}_{i}=14\frac{p}{\rho}X^{(3)}_{i}$;
$\tilde{h}^{(6)}=0$; |  $\widetilde{X}^{(6)}=21\frac{p}{\rho}\widetilde{X}^{(4)}$;
$\tilde{h}_{i}^{(7)}=0$; |  $X^{(7)}_{i}=27\frac{p}{\rho}X^{(5)}_{i}-189\frac{p^{2}}{\rho^{2}}X_{i}^{(3)}$;
$\tilde{h}^{(8)}=0$; |  $\widetilde{X}^{(8)}=36\frac{p}{\rho}\widetilde{X}^{(6)}-378\frac{p^{2}}{\rho^{2}}\widetilde{X}^{(4)}$;
$\tilde{h}_{i}^{(9)}=0$; |  $X^{(9)}_{i}=44\frac{p}{\rho}X^{(7)}_{i}-594\frac{p^{2}}{\rho^{2}}X^{(5)}_{i}+2772\frac{p^{3}}{\rho^{3}}X^{(3)}_{i}$.
Table 8.3: Summary of (MHD) Hermite closures, together with corresponding fluid closures. Species indices “a” are dropped. The usual heat flux $q_{i}=X_{i}^{(3)}/2$. Note that beyond the 4th-order moment both classes start to differ. It can be shown that erroneously prescribing closures at the last retained moment such as $X^{(5)}_{i}=0$ or $\widetilde{X}^{(6)}=0$ leads to unphysical instabilities (unless one prescribes $X^{(3)}_{i}$ or $\widetilde{X}^{(4)}=0$ as well), which is later demonstrated in Appendix B.8, Table B.1. A general form for closures corresponding to $\tilde{h}^{(2n+1)}_{i}=0$ and $\tilde{h}^{(2n)}=0$ is given by (208). An analogous table can be constructed for CGL parallel closures; see Appendix B.9, Table B.2. Hermite closures | Fluid closures
---|---
$\hat{h}^{(2)}_{ij}=0$; |  $\Pi_{ij}^{(2)}=0$;
$\hat{h}^{(4)}_{ij}=0$; |  $\Pi_{ij}^{(4)}=7\frac{p}{\rho}\Pi_{ij}^{(2)}$;
$\hat{h}^{(6)}_{ij}=0$; |  $\Pi_{ij}^{(6)}=18\frac{p}{\rho}\Pi_{ij}^{(4)}-63\frac{p^{2}}{\rho^{2}}\Pi_{ij}^{(2)}$;
$\hat{h}^{(8)}_{ij}=0$; |  $\Pi_{ij}^{(8)}=33\frac{p}{\rho}\Pi_{ij}^{(6)}-297\frac{p^{2}}{\rho^{2}}\Pi_{ij}^{(4)}+693\frac{p^{3}}{\rho^{3}}\Pi_{ij}^{(2)}$;
Table 8.4: Similar to Table 8.3, but for Hermite closures
$\hat{h}^{(2n)}_{ij}=0$. A general form for closures corresponding to
$\hat{h}^{(2n)}_{ij}=0$ is given by (209).

Here we need to address one incorrect interpretation that we used in some of
our previous papers. In the last paragraph of Hunana et al. (2018), and also
in Hunana et al. (2019a, b) it is claimed that Landau fluid closures are
necessary to go beyond the 4th-order moment in the fluid hierarchy. This
interpretation was obtained in the CGL framework for parallel moments by
considering closures at the last retained moment $\widetilde{X}^{(2n)}_{a}=0$
and $\vec{\boldsymbol{X}}^{(2n+1)}_{\parallel a}=0$. It was shown (see
detailed proof in Section 12.2 in Hunana et al. (2019b)) that beyond the 4th-
order moment, all fluid models become unstable if these closures are used. The
proof is constructed correctly. What is incorrect is the interpretation, that
the proof implies that Landau fluid closures are required to overcome this
issue. The much simpler Hermite closures overcome this difficulty as well.

In another words, beyond the 4th-order moment it is not possible to cut the
fluid hierarchy by simply neglecting the next order moment with closures such
as $\vec{\boldsymbol{X}}^{(5)}_{a}=0$ or $\widetilde{X}^{(6)}_{a}=0$, and such
closures should be viewed as erroneous. For the CGL model the closures have
different coefficients than for the MHD model because the moments are defined
differently (a brief summary is given in Appendix B.9, Table B.2). The CGL
closures will be addressed in detail in a separate publication.

Importantly, the problem also disappears when one decouples the fluid
hierarchy. For example, higher-order Laguerre (Hermite) schemes that are
typically used to obtain more precise transport coefficients for
$\vec{\boldsymbol{q}}_{a}$ and $\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$,
neglect all the scalar perturbations
$\widetilde{X}^{(4)}_{a}=\cdots=\widetilde{X}^{(2n)}_{a}=0$, together with
neglecting coupling between heat fluxes and stress-tensors. In our formulation
this yields a system

$\displaystyle\frac{d_{a}}{dt}\vec{\boldsymbol{X}}^{(2n+1)}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\vec{\boldsymbol{X}}^{(2n+1)}_{a}+(2n+3)!!\frac{(n)}{3}\frac{p_{a}^{n}}{\rho_{a}^{n-1}}\nabla\Big{(}\frac{p_{a}}{\rho_{a}}\Big{)}$
$\displaystyle\qquad=\vec{\boldsymbol{\mathcal{Q}}}^{(2n+1)}_{a}-\frac{(2n+3)!!}{3}\frac{p_{a}^{n}}{\rho_{a}^{n}}\boldsymbol{R}_{a};$
(210)
$\displaystyle\frac{d_{a}}{dt}\bar{\bar{\boldsymbol{\Pi}}}^{(2n)}_{a}+\Omega_{a}\big{(}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{\Pi}}}^{(2n)}_{a}\big{)}^{S}+\frac{(2n+3)!!}{15}\frac{p^{n}_{a}}{\rho_{a}^{n-1}}\bar{\bar{\boldsymbol{W}}}_{a}$
$\displaystyle\qquad=\bar{\bar{\boldsymbol{\mathcal{Q}}}}^{(2n)}_{a}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}Q^{(2n)}_{a}.$
(211)

Closures (208), (209) are not required, because the equations are de-coupled.
We did not calculate collisional contributions for higher-order moments, but
in the semi-linear approximation equations (210)-(211) remain de-coupled and
represent two independent hierarchies. An essential feature of the Landau (or
the Boltzmann) collisional operator is that the operator couples all the heat
fluxes together, and it also couples all the stress-tensors together. Thus by
going higher and higher in the fluid hierarchy does not create new
contributions in a quasi-static approximation, but yields increasingly precise
transport coefficients for $\vec{\boldsymbol{q}}_{a}$ and
$\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}$. Also, because the momentum exchange
rates $\boldsymbol{R}_{a}$ contain contributions from all the heat fluxes
$\vec{\boldsymbol{X}}^{(3)}_{a}\ldots\vec{\boldsymbol{X}}^{(2n+1)}_{a}$, they
become increasingly precise as well. System (210)-(211) nicely clarifies how
higher-order schemes can be viewed. Reinstating the coupling between heat
fluxes and viscosity-tensors introduces additional contributions but does not
change the transport coefficients of the de-coupled system. A brief comparison
of various models is presented in Appendix I.

###  8.7 Inclusion of gravity

We have not explicitly considered the force of gravity during our calculations
in the Appendix, nevertheless, its inclusion is trivial. With the
gravitational acceleration $\boldsymbol{G}$ included, the Boltzmann equation
reads

$\frac{\partial f_{a}}{\partial t}+\boldsymbol{v}\cdot\nabla
f_{a}+\Big{[}\boldsymbol{G}+\frac{eZ_{a}}{m_{a}}(\boldsymbol{E}+\frac{1}{c}\boldsymbol{v}\times\boldsymbol{B})\Big{]}\cdot\nabla_{v}f_{a}=C(f_{a}).$
(212)

We use big $\boldsymbol{G}$ instead of small $\boldsymbol{g}$ to clearly
distinguish it from the heat flux $\boldsymbol{q}$. Gravity does not enter the
collisional operator, and collisional integrals with the right hand side are
not effected. Gravity enters the left hand side, and when the Boltzmann
equation is integrated gravity of course enters the fluid hierarchy of
moments. With the two exceptions of the density equation and the pressure
tensor equation, gravity enters evolution equations for all other moments,
analogously as the electric field does. An explicit collisionless equation for
a general n-th order moment with the electric field present is for example
equation (12.13) of Hunana et al. (2019b). Because no Maxwell’s equations are
used in deriving the fluid hierarchy, the presence of gravity can be accounted
for by simply replacing

$\frac{eZ_{a}}{m_{a}}\boldsymbol{E}\to\boldsymbol{G}+\frac{eZ_{a}}{m_{a}}\boldsymbol{E}.$
(213)

Furthermore, such a hierarchy is not very useful because the evolution
equation for an n-th order moment is coupled with “n” momentum equations.
Subtracting these momentum equations yields final equation (12.16) in Hunana
et al. (2019b), where the electric field is not present, meaning that gravity
is not present either. In other words, the collisionless equation (12.16) of
Hunana et al., as well as our new collisional equation (A12) remain valid in
the presence of gravitational force. The inclusion of gravity to the entire
model is thus achieved trivially by adding $-\boldsymbol{G}$ into the left
hand side of the momentum equation (7) (which we have done), and no additional
calculations are required. In the main text, the only other equation which
contains gravity is the electric field equation (96).

###  8.8 Precision of $m_{e}/m_{i}$ expansions (unmagnetized proton-electron
plasma)

The multi-fluid formulation is also an excellent tool to double-check the
precision of $m_{e}/m_{i}$ expansions. It is possible to again consider a one
ion-electron plasma, but this time calculate the transport coefficients
precisely, without any expansions in the smallness of $m_{e}/m_{i}$. As an
example we consider an unmagnetized proton-electron plasma ($Z_{p}=1$,
$m_{p}/m_{e}=1836.15267$) with similar temperatures $T_{e}=T_{p}=T_{ep}$.
Charge neutrality implies $n_{e}=n_{p}$ and so $p_{e}=p_{p}$. We however
maintain $\nabla T_{e}\neq\nabla T_{p}$, because the gradients can be
different. We first calculate heat fluxes. For clarity, we are solving 4
coupled evolution equations, which are explicitly given in Appendix N; see
equations (N1)-(N4).

Precise calculation should not use simplified collisional times (182) where
expansions in $m_{e}/m_{i}$ have been made, but exact collisional times (178)
with numerical values $\nu_{ee}=0.707299\nu_{ep}$ and
$\nu_{pp}=0.0165063\nu_{ep}$ (we take $\ln\Lambda$ to be constant). The quasi-
static approximation then yields heat fluxes

$\displaystyle\vec{\boldsymbol{q}}_{e}$ $\displaystyle=$
$\displaystyle\Big{[}-3.159370\nabla T_{e}+8.301\times 10^{-6}\nabla
T_{p}\Big{]}\frac{p_{e}}{m_{e}\nu_{ep}}+0.711046p_{e}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{X}}^{(5)}_{e}$ $\displaystyle=$
$\displaystyle\Big{[}-110.5793\nabla T_{e}+1.376\times 10^{-3}\nabla
T_{p}\Big{]}\frac{p_{e}^{2}}{\rho_{e}m_{e}\nu_{ep}}+18.78249\frac{p_{e}^{2}}{\rho_{e}}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{q}}_{p}$ $\displaystyle=$
$\displaystyle\Big{[}-3.302411\nabla T_{p}+0.2516\times 10^{-3}\nabla
T_{e}\Big{]}\frac{p_{p}}{m_{p}\nu_{pp}}+0.206535\times
10^{-4}p_{p}\delta\boldsymbol{u};$
$\displaystyle\vec{\boldsymbol{X}}^{(5)}_{p}$ $\displaystyle=$
$\displaystyle\Big{[}-103.3984\nabla T_{p}+0.7863\times 10^{-2}\nabla
T_{e}\Big{]}\frac{p_{p}^{2}}{\rho_{p}m_{p}\nu_{pp}}+0.646475\times
10^{-3}\frac{p_{p}^{2}}{\rho_{p}}\delta\boldsymbol{u},$ (214)

where $\delta\boldsymbol{u}=(\boldsymbol{u}_{e}-\boldsymbol{u}_{p})$. For the
electron heat flux $\vec{\boldsymbol{q}}_{e}$, note the difference of the
thermal conductivity 3.1594 from the Braginskii value 3.1616. The difference
is caused by calculating the mass-ratio coefficients (27), (28) exactly
without $m_{e}/m_{p}$ expansions, together with slightly different ratios of
frequencies (a less-precise calculation with neglecting proton-proton
collisions by $\nu_{pp}=0$ and using simplified $\nu_{ee}=\nu_{ep}/\sqrt{2}$
yields 3.1600).

For the proton heat flux $\vec{\boldsymbol{q}}_{p}$, the relatively large
difference between the thermal conductivity 3.302 and Braginskii self-
collisional value $125/32=3.906$ is caused by the proton-electron collisions.
Similarly for the $\vec{\boldsymbol{X}}^{(5)}_{p}$, where the self-collisional
value is $2975/24=123.96$. Calculating the coupled system exactly has a nice
advantage that one can calculate the momentum exchange rates in two different
ways

$\displaystyle\boldsymbol{R}_{e}$ $\displaystyle=$
$\displaystyle\nu_{ep}\Big{\\{}-\rho_{e}\delta\boldsymbol{u}+\frac{\mu_{ep}}{T_{ep}}\Big{[}V_{ep(1)}\vec{\boldsymbol{q}}_{e}-V_{ep(2)}\frac{\rho_{e}}{\rho_{p}}\vec{\boldsymbol{q}}_{p}\Big{]}-\frac{3}{56}\Big{(}\frac{\mu_{ep}}{T_{ep}}\Big{)}^{2}\Big{[}\vec{\boldsymbol{X}}^{(5)}_{e}-\frac{\rho_{e}}{\rho_{p}}\vec{\boldsymbol{X}}^{(5)}_{p}\Big{]}\Big{\\}};$
$\displaystyle\boldsymbol{R}_{p}$ $\displaystyle=$
$\displaystyle\nu_{pe}\Big{\\{}+\rho_{p}\delta\boldsymbol{u}+\frac{\mu_{ep}}{T_{ep}}\Big{[}V_{pe(1)}\vec{\boldsymbol{q}}_{p}-V_{pe(2)}\frac{\rho_{p}}{\rho_{e}}\vec{\boldsymbol{q}}_{e}\Big{]}-\frac{3}{56}\Big{(}\frac{\mu_{ep}}{T_{ep}}\Big{)}^{2}\Big{[}\vec{\boldsymbol{X}}^{(5)}_{p}-\frac{\rho_{p}}{\rho_{e}}\vec{\boldsymbol{X}}^{(5)}_{e}\Big{]}\Big{\\}},$
(215)

and both options yield the same result

$\boldsymbol{R}_{e}=-\boldsymbol{R}_{p}=-0.711046n_{e}\nabla
T_{e}-0.2065\times 10^{-4}n_{e}\nabla
T_{p}-0.513306\rho_{e}\nu_{ep}\delta\boldsymbol{u}.$ (216)

Viscosities of proton-electron plasma are (for clarity, we are solving 4
equations in 4 unknowns, explicitly given by (N5)-(N8))

$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{e}$ $\displaystyle=$
$\displaystyle\Big{[}-0.730622\bar{\bar{\boldsymbol{W}}}_{e}-0.2800\times
10^{-2}\bar{\bar{\boldsymbol{W}}}_{p}\Big{]}\frac{p_{e}}{\nu_{ei}};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{e}$ $\displaystyle=$
$\displaystyle\Big{[}-6.542519\bar{\bar{\boldsymbol{W}}}_{e}+3.1509\times
10^{-2}\bar{\bar{\boldsymbol{W}}}_{p}\Big{]}\frac{p_{e}^{2}}{\rho_{e}\nu_{ei}};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{p}$ $\displaystyle=$
$\displaystyle\Big{[}-0.892105\bar{\bar{\boldsymbol{W}}}_{p}-0.4621\times
10^{-4}\bar{\bar{\boldsymbol{W}}}_{e}\Big{]}\frac{p_{p}}{\nu_{pp}};$
$\displaystyle\bar{\bar{\boldsymbol{\Pi}}}^{(4)}_{p}$ $\displaystyle=$
$\displaystyle\Big{[}-7.250870\bar{\bar{\boldsymbol{W}}}_{p}-0.3759\times
10^{-3}\bar{\bar{\boldsymbol{W}}}_{e}\Big{]}\frac{p_{p}^{2}}{\rho_{p}\nu_{pp}},$
(217)

and for proton species the relatively large differences from self-collisional
values $1025/1068=0.960$ and $8435/1068=7.898$ are again caused by proton-
electron collisions. In Appendix N, we consider another examples of coupling
between two species, and we calculate heat fluxes and viscosities for protons
& alpha particles (fully ionized Helium), and for the deuterium-tritium plasma
used in plasma fusion.

###  8.9 Limitations of our approach

It is important to clarify the limitations of our model. In the highly-
collisional regime, our limitations are the same as for the model of
Braginskii (1965). For example, we describe only Coulomb collisions and we do
not take into account ionization & recombination and radiative transfer.
Additionally, our approach shows that coupling of stress-tensors and heat
fluxes should be ideally investigated with the 22-moment model. Even though
this model is fully formulated in Section 7, including its collisional
contributions calculated with the Landau operator, we did not use this model
to further explore the resulting coupling.

####  8.9.1 Weakly collisional regime: expansions around bi-Maxwellians

The situation becomes more complicated in the weakly-collisional regime where
there might not be enough collisions to keep the distribution function
sufficiently close to the equilibrium Maxwellian $f_{a}^{(0)}$. The
distribution function might evolve to such an extent that the core assumptions
in the entire derivation break down, i.e. equation (1) looses its validity. A
better approach is then to consider expansions similar to equation (1), but
around a bi-Maxwellian $f_{a}^{(0)}$ (see e.g. Oraevskii et al. (1968);
Chodura & Pohl (1971); Demars & Schunk (1979); Barakat & Schunk (1982), and
references therein), which can handle much larger departures from the highly-
collisional Maxwellian distribution. In order to point out the differences and
difficulties associated with this approach, it is of interest to briefly
describe how expansions around an anisotropic bi-Maxwellian would look like.
The simplest anisotropic model is known as the CGL, after the pioneering work
of Chew, Goldberger and Low (Chew et al., 1956). The difference with our
current approach starts with the decomposition of the pressure tensor
$p_{ij}^{a}$ defined in (A2), and the decomposition reads

$\displaystyle\textrm{isotropic:}\qquad\bar{\bar{\boldsymbol{p}}}_{a}$
$\displaystyle=$ $\displaystyle
p_{a}\bar{\bar{\boldsymbol{I}}}+\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)};$ (218)
$\displaystyle\textrm{anisotropic:}\qquad\bar{\bar{\boldsymbol{p}}}_{a}$
$\displaystyle=$ $\displaystyle p_{\parallel
a}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}+p_{\perp
a}(\bar{\bar{\boldsymbol{I}}}-\boldsymbol{\hat{b}}\boldsymbol{\hat{b}})+\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)\textrm{CGL}}$
$\displaystyle=$ $\displaystyle p_{a}\bar{\bar{\boldsymbol{I}}}+(p_{\parallel
a}-p_{\perp
a})\Big{(}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}-\frac{\bar{\bar{\boldsymbol{I}}}}{3}\Big{)}+\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)\textrm{CGL}},$

with scalar pressures

$\displaystyle p_{\parallel a}$ $\displaystyle=$
$\displaystyle\bar{\bar{\boldsymbol{p}}}_{a}:\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}=m_{a}\int
c_{\parallel a}^{2}f_{a}d^{3}v;\qquad p_{\perp
a}=\bar{\bar{\boldsymbol{p}}}_{a}:\bar{\bar{\boldsymbol{I}}}_{\perp}/2=\frac{m_{a}}{2}\int|\boldsymbol{c}_{\perp
a}|^{2}f_{a}d^{3}v.$ (220)

Directly from the above definitions, the stress-tensors satisfy

$\textrm{Tr}\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}=\textrm{Tr}\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)\textrm{CGL}}=0;\qquad\bar{\bar{\boldsymbol{\Pi}}}^{(2)}_{a}:\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\neq
0;\qquad\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)\textrm{CGL}}:\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}=0,$
(221)

and while $\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)}$ has 5 independent
components, $\bar{\bar{\boldsymbol{\Pi}}}_{a}^{(2)\textrm{CGL}}$ has only 4.

The decomposition of the heat flux tensor $q_{ijk}^{a}$ defined by equation
(A2) is slightly more complicated. In an arbitrary-collisional regime one
needs to define two heat flux vectors

$\boldsymbol{S}^{\parallel}_{a}=\bar{\bar{\boldsymbol{q}}}_{a}:\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}=m_{a}\int
c_{\parallel
a}^{2}\boldsymbol{c}_{a}f_{a}d^{3}v;\qquad\boldsymbol{S}^{\perp}_{a}=\bar{\bar{\boldsymbol{q}}}_{a}:\bar{\bar{\boldsymbol{I}}}_{\perp}/2=\frac{m_{a}}{2}\int|\boldsymbol{c}_{\perp
a}|^{2}\boldsymbol{c}_{a}f_{a}d^{3}v.$ (222)

These heat flux vectors are further split by projecting them along the
$\boldsymbol{\hat{b}}$, which defines the _gyrotropic_ (scalar) heat fluxes
$q_{\parallel a}$ & $q_{\perp a}$, and the perpendicular projection defines
the _non-gyrotropic heat flux vectors_ $\boldsymbol{S}^{\parallel}_{\perp a}$
& $\boldsymbol{S}^{\perp}_{\perp a}$, according to

$\displaystyle q_{\parallel a}$ $\displaystyle=$
$\displaystyle\boldsymbol{\hat{b}}\cdot\boldsymbol{S}^{\parallel}_{a}=m_{a}\int
c_{\parallel a}^{2}c_{\parallel a}f_{a}d^{3}v;\qquad q_{\perp
a}=\boldsymbol{\hat{b}}\cdot\boldsymbol{S}^{\perp}_{a}=\frac{m_{a}}{2}\int|\boldsymbol{c}_{\perp
a}|^{2}c_{\parallel a}f_{a}d^{3}v;$
$\displaystyle\boldsymbol{S}^{\parallel}_{\perp a}$ $\displaystyle=$
$\displaystyle\bar{\bar{\boldsymbol{I}}}_{\perp}\cdot\boldsymbol{S}^{\parallel}_{a}=m_{a}\int
c_{\parallel a}^{2}\boldsymbol{c}_{\perp
a}f_{a}d^{3}v;\qquad\boldsymbol{S}^{\perp}_{\perp
a}=\bar{\bar{\boldsymbol{I}}}_{\perp}\cdot\boldsymbol{S}^{\perp}_{a}=\frac{m_{a}}{2}\int|\boldsymbol{c}_{\perp
a}|^{2}\boldsymbol{c}_{\perp a}f_{a}d^{3}v.$ (223)

The following relations then hold $\boldsymbol{S}^{\parallel}_{a}=q_{\parallel
a}\boldsymbol{\hat{b}}+\boldsymbol{S}^{\parallel}_{\perp a}$; &
$\boldsymbol{S}^{\perp}_{a}=q_{\perp
a}\boldsymbol{\hat{b}}+\boldsymbol{S}^{\perp}_{\perp a}$; together with
$\boldsymbol{\hat{b}}\cdot\boldsymbol{S}^{\parallel}_{\perp a}=0$; &
$\boldsymbol{\hat{b}}\cdot\boldsymbol{S}^{\perp}_{\perp a}=0$. The two
different decompositions of the entire heat flux tensor then read

$\displaystyle\textrm{isotropic:}\qquad\bar{\bar{\boldsymbol{q}}}_{a}$
$\displaystyle=$
$\displaystyle\frac{2}{5}\Big{[}\vec{\boldsymbol{q}}_{a}\bar{\bar{\boldsymbol{I}}}\Big{]}^{S}+\bar{\bar{\boldsymbol{\sigma}}}_{a}^{\prime};$
(224) $\displaystyle\textrm{anisotropic:}\qquad\bar{\bar{\boldsymbol{q}}}_{a}$
$\displaystyle=$ $\displaystyle q_{\parallel
a}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}+q_{\perp
a}\Big{[}\boldsymbol{\hat{b}}\bar{\bar{\boldsymbol{I}}}_{\perp}\Big{]}^{S}+\Big{[}\boldsymbol{S}^{\parallel}_{\perp
a}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\Big{]}^{S}+\frac{1}{2}\Big{[}\boldsymbol{S}^{\perp}_{\perp
a}\bar{\bar{\boldsymbol{I}}}_{\perp}\Big{]}^{S}+\bar{\bar{\boldsymbol{\sigma}}}_{a};$
(225) $\displaystyle=$ $\displaystyle q_{\perp
a}\big{[}\boldsymbol{\hat{b}}\bar{\bar{\boldsymbol{I}}}\big{]}^{S}+(q_{\parallel
a}-3q_{\perp
a})\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}+\frac{1}{2}\Big{[}\boldsymbol{S}^{\perp}_{\perp
a}\bar{\bar{\boldsymbol{I}}}\Big{]}^{S}+\Big{[}\big{(}\boldsymbol{S}^{\parallel}_{\perp
a}-\frac{\boldsymbol{S}^{\perp}_{\perp
a}}{2}\big{)}\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}\Big{]}^{S}+\bar{\bar{\boldsymbol{\sigma}}}_{a},$

where both $\bar{\bar{\boldsymbol{\sigma}}}_{a}^{\prime}$ and
$\bar{\bar{\boldsymbol{\sigma}}}_{a}$ are traceless. Neglecting these
traceless contributions, the isotropic approach accounts for 3 (out of 10)
scalar components of $\bar{\bar{\boldsymbol{q}}}_{a}$ and represents a
13-moment model (1 density, 3 velocity, 1 scalar pressure, 5 stress-tensor
components, 3 heat flux $\vec{\boldsymbol{q}}_{a}$ components). The
anisotropic approach accounts for 6 scalar components of
$\bar{\bar{\boldsymbol{q}}}_{a}$ and represents a 16-moment model, described
by sixteen scalar evolution equations (1 density, 3 velocity, 2 scalar
pressures, 4 stress-tensor components, 3 for each heat flux vector
$\boldsymbol{S}^{\parallel}_{a}$ & $\boldsymbol{S}^{\perp}_{a}$.)
Unfortunately, such a complicated decomposition of the heat flux is necessary
in an arbitrary-collisional regime, and we only used decomposition (224). For
clarity, direct relation with the usual heat flux vector
$\vec{\boldsymbol{q}}_{a}$ reads

$\vec{\boldsymbol{q}}_{a}=\frac{1}{2}\boldsymbol{S}^{\parallel}_{a}+\boldsymbol{S}^{\perp}_{a}=\Big{(}\frac{1}{2}q_{\parallel
a}+q_{\perp
a}\Big{)}\boldsymbol{\hat{b}}+\frac{1}{2}\boldsymbol{S}^{\parallel}_{\perp
a}+\boldsymbol{S}^{\perp}_{\perp a}.$ (226)

Note that both $q_{\parallel a}$ & $q_{\perp a}$ denote components along the
$\boldsymbol{\hat{b}}$. The highly-collisional limit is achieved by
$q_{\parallel a}=3q_{\perp a}$ and $\boldsymbol{S}^{\parallel}_{\perp
a}=\boldsymbol{S}^{\perp}_{\perp a}/2$, in which case
$\vec{\boldsymbol{q}}_{a}=(5/2)q_{\perp
a}\boldsymbol{\hat{b}}+(5/4)\boldsymbol{S}^{\perp}_{\perp a}$ or equivalently
$\vec{\boldsymbol{q}}_{a}=(5/6)q_{\parallel
a}\boldsymbol{\hat{b}}+(5/2)\boldsymbol{S}^{\parallel}_{\perp a}$. We used the
same notation as for example collisionless papers by Passot & Sulem (2007);
Sulem & Passot (2015); Hunana et al. (2019a, b).

These anisotropic decompositions must be retained in an arbitrary-collisional
regime. However, calculations with the Landau (Boltzmann) collisional
operators then become very complicated. Notably, Chodura & Pohl (1971); Demars
& Schunk (1979); Barakat & Schunk (1982) used the anisotropic 16-moment model
as described above and calculated the collisional contributions for several
interaction potentials. Judging from the papers above, maintaining the
precision of our current model (where the 4th & 5th order moments are
considered) and extending it to an anisotropic (bi-Maxwellian) regime seems to
be so complicated, that it might not be worth the effort. Curiously, in a
simplified spherically symmetric radial geometry Cuperman et al. (1980, 1981);
Cuperman & Dryer (1985) considered what seems like a mixture of anisotropic
and isotropic moments, with anisotropic temperatures, isotropic heat flux
vector, and the parallel (anisotropic) perturbation of the 4th-order moment
(which we call $\widetilde{r}_{\parallel\parallel a}$).

####  8.9.2 Landau fluid closures for the collisionless case

In contrast to the free-streaming formula of Hollweg (1974, 1976), in plasma
physics the collisionless heat flux is typically associated with the
phenomenon of Landau damping. For example, collisionless linear kinetic theory
expanded around bi-Maxwellian with mean _zero drifts_ in gyrotropic limit
yields in Fourier space perturbation of the distribution function
$f_{a}=f_{a}^{(0)}(1+\chi_{a})$ in the following form

$\chi_{a}=\frac{B_{\parallel}^{(1)}}{B_{0}}\frac{m_{a}}{2T_{\perp
a}^{(0)}}\bigg{[}v_{\perp}^{2}+\frac{T_{\perp a}^{(0)}}{T_{\parallel
a}^{(0)}}\frac{k_{\parallel}v_{\parallel}v_{\perp}^{2}}{(\omega-
k_{\parallel}v_{\parallel})}\bigg{]}+\Phi\frac{eZ_{a}}{T_{\parallel
a}^{(0)}}\frac{k_{\parallel}v_{\parallel}}{(\omega-
k_{\parallel}v_{\parallel})},$ (227)

with the electrostatic potential $\Phi=iE_{\parallel}^{(1)}/k_{\parallel}$.
Integrating (227) then yields a parallel collisionless heat flux

$q_{\parallel a}^{(1)}=-v_{\textrm{th}\parallel a}n_{a}^{(0)}T_{\parallel
a}^{(0)}\textrm{sign}(k_{\parallel})\Big{(}\zeta_{a}+2\zeta_{a}^{3}R(\zeta_{a})-3\zeta_{a}R(\zeta_{a})\Big{)}\Big{[}\frac{B_{\parallel}^{(1)}}{B_{0}}\frac{T_{\perp
a}^{(0)}}{T_{\parallel a}^{(0)}}+\Phi\frac{eZ_{a}}{T_{\parallel
a}^{(0)}}\Big{]},$ (228)

with variable $\zeta_{a}=\omega/(|k_{\parallel}|v_{\textrm{th}\parallel a})$;
parallel thermal speed $v_{\textrm{th}\parallel a}=\sqrt{2T_{\parallel
a}/m_{a}}$, plasma response function $R(\zeta_{a})=1+\zeta_{a}Z(\zeta_{a})$
and plasma dispersion function
$Z(\zeta_{a})=i\sqrt{\pi}\exp(-\zeta_{a}^{2})[1+\textrm{erf}(i\zeta_{a})]$.
Such a kinetic answer can be expressed in fluid variables by searching for
Landau fluid closures, for example by replacing the $R(\zeta_{a})$ with its
3-pole Padé approximants

$\displaystyle R_{3,2}(\zeta_{a})$ $\displaystyle=$
$\displaystyle\frac{1-i\frac{\sqrt{\pi}}{2}\zeta_{a}}{1-i\frac{3\sqrt{\pi}}{2}\zeta_{a}-2\zeta_{a}^{2}+i\sqrt{\pi}\zeta_{a}^{3}};\qquad
R_{3,1}(\zeta_{a})=\frac{1-i\frac{(4-\pi)}{\sqrt{\pi}}\zeta_{a}}{1-i\frac{4}{\sqrt{\pi}}\zeta_{a}-2\zeta_{a}^{2}+2i\frac{(4-\pi)}{\sqrt{\pi}}\zeta_{a}^{3}}.$
(229)

The use of Padé approximants allows one to express (228) through lower-order
moments and eliminate the explicit dependence on $\zeta_{a}$, yielding
collisionless heat fluxes in Fourier space

$\displaystyle R_{3,2}(\zeta_{a}):\qquad q_{\parallel a}^{(1)}$
$\displaystyle=$
$\displaystyle-i\frac{2}{\sqrt{\pi}}n_{a}^{(0)}v_{\textrm{th}\parallel
a}\textrm{sign}(k_{\parallel})T_{\parallel a}^{(1)};$ (230) $\displaystyle
R_{3,1}(\zeta_{a}):\qquad q_{\parallel a}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{3\pi-8}{4-\pi}\,p^{(0)}_{\parallel a}u^{(1)}_{\parallel
a}-i\frac{\sqrt{\pi}}{4-\pi}n^{(0)}_{a}v_{\textrm{th}\parallel
a}\textrm{sign}(k_{\parallel})T^{(1)}_{\parallel a},$ (231)

where $T^{(1)}_{\parallel a}$ is perturbed temperature, and
$u^{(1)}_{\parallel a}$ is perturbed fluid velocity (mean value
$u^{(0)}_{\parallel a}=0$ is assumed). The heat flux closure (230) was
obtained by Hammett & Perkins (1990) & Snyder et al. (1997) and closure (231)
is equation (2) in Hunana et al. (2018) (or equation (3.211) in Hunana et al.
(2019a)). In real space these collisionless heat fluxes become

$\displaystyle R_{3,2}(\zeta_{a}):\qquad q_{\parallel a}(z)$ $\displaystyle=$
$\displaystyle-\frac{2}{\pi^{3/2}}n_{a}^{(0)}v_{\textrm{th}\parallel
a}\textrm{V.P.}\int_{0}^{\infty}\frac{T^{(1)}_{\parallel
a}(z+z^{\prime})-T^{(1)}_{\parallel a}(z-z^{\prime})}{z^{\prime}}dz^{\prime};$
(232) $\displaystyle R_{3,1}(\zeta_{a}):\qquad q_{\parallel a}(z)$
$\displaystyle=$ $\displaystyle\frac{3\pi-8}{4-\pi}p^{(0)}_{\parallel
a}u^{(1)}_{\parallel a}-\frac{n^{(0)}_{a}v_{\textrm{th}\parallel
a}}{\sqrt{\pi}(4-\pi)}\textrm{V.P.}\int_{0}^{\infty}\frac{T^{(1)}_{\parallel
a}(z+z^{\prime})-T^{(1)}_{\parallel a}(z-z^{\prime})}{z^{\prime}}dz^{\prime},$
(233)

where the non-locality presents itself as an integral over the entire magnetic
field line, where temperatures everywhere along that field line matter to
determine the heat flux at a specific spatial point. Note that the thermal
part of (233) is almost two times larger then (232). The Cauchy principal
value can be replaced by $\lim_{\epsilon\to+0}\int_{\epsilon}^{\infty}$. This
approach thus indeed allows one to have expressions for collisionless heat
fluxes in a quasi-static approximation. However, as is well-known these
expressions are not very precise with respect to kinetic theory. For example,
the precision can be easily compared by plotting normalized heat fluxes
$\hat{q}_{\parallel
a}=\zeta_{a}+2\zeta_{a}^{3}R(\zeta_{a})-3\zeta_{a}R(\zeta_{a})$, which is
shown in Figure 1. Weakly-damped regime with real valued $\zeta_{a}$ is
considered. The left panel shows the imaginary part of $\hat{q}_{\parallel a}$
and the right panel shows the real part of $\hat{q}_{\parallel a}$. Exact
kinetic heat flux is solid black line, heat flux $R_{3,2}$ is dashed magenta
line and heat flux $R_{3,1}$ is dashed cyan line. For comparison, higher-order
fluid models with approximants $R_{5,3}$ (dotted blue line) and $R_{7,5}$
(dashed red line) are shown as well (see equation (A11) & (A38) in Hunana et
al. (2019a)). The $R_{5,3}$ model represents a dynamic closure at the 4th-
order moment and the $R_{7,5}$ model represents a dynamic closure at the 6th-
order moment, given by equation (5) & (8) of Hunana et al. 2018. The heat
fluxes in these higher-order models are thus described by their usual
evolution equations, nevertheless, their precision can be compared with the
same technique. Which quasi-static heat flux is a better choice depends on the
value of $\zeta_{a}$, because the $R_{3,1}$ has a higher power-series
precision (for small $\zeta_{a}$) and the $R_{3,2}$ has a higher asymptotic-
series precision (for large $\zeta_{a}$). Regime $\zeta_{a}\ll 1$ can be
viewed as isothermal and regime $\zeta_{a}\gg 1$ can be viewed as adiabatic.
In the left panel of Figure 1 the $R_{3,1}$ is more precise up to roughly
$\zeta_{a}=2.3$ and in the right panel up to $\zeta_{a}=1.6$. For larger
$\zeta_{a}$ values than shown, the $R_{3,1}$ heat flux converges much slower
to the correct zero values than the $R_{3,2}$, especially for the real part.

Figure 1: Comparison of normalized collisionless heat fluxes
$\hat{q}_{\parallel a}=\zeta_{a}+(2\zeta_{a}^{3}-3\zeta_{a})R(\zeta_{a})$ in a
weakly-damped regime with real valued $\zeta_{a}$. Left panel: imaginary part
of $\hat{q}_{\parallel a}$. Right panel: real part of $\hat{q}_{\parallel a}$.
Colors are described in the text. Our Braginskii-type models do not contain
these collisionless heat fluxes.

The major obstacle in precision for the quasi-static heat fluxes of Landau
fluid models actually comes from the perpendicular heat flux $q_{\perp a}$
(which is along the $\boldsymbol{\hat{b}}$), because only a closure of Snyder
et al. (1997) with a crude Padé approximant
$R_{1}(\zeta_{a})=1/(1-i\sqrt{\pi}\zeta_{a})$ is available. As a consequence,
for large $\zeta_{a}$ values the quasi-static heat flux $q_{\perp a}$ fails to
disappear and instead converges to a constant value. To recover the adiabatic
behavior for $q_{\perp a}$, one has to abandon the idea of quasi-static
$q_{\perp a}$ and consider its evolution equation, with a closure at the 4th-
order moment. There is a vast amount of literature about Landau fluids with
various approaches; see e. g. Hammett & Perkins (1990); Hammett et al. (1992);
Snyder et al. (1997); Snyder & Hammett (2001); Goswami et al. (2005); Passot &
Sulem (2007); Passot et al. (2012); Sulem & Passot (2015); Joseph & Dimits
(2016); Hunana et al. (2018); Ji & Joseph (2018); Chen et al. (2019); Wang et
al. (2019) and references therein, where some authors also include collisional
effects. For a simple introductory guide to collisionless Landau fluids see
Hunana et al. (2019a). As a side note, Landau fluid closures are not
constructed with any specific mode in mind (as incorrectly criticized for
example by Scudder (2021)). The closures are constructed universally for all
the modes, so that numerical simulations can be performed; see e.g. Perrone et
al. (2018). Interestingly, as discussed by Meyrand et al. (2019), from a non-
linear perspective the effect of Landau damping might be canceled out by the
effect of plasma echo. From a linear perspective, the presence of drifts also
modifies the Landau damping, because the variable $\zeta_{a}$ which enters the
plasma response function $R(\zeta_{a})$ then contains the drift velocity
$u_{\parallel a}$. For sufficiently large drifts the sound mode can be
generated by the current-driven ion-acoustic instability; see e.g. Gurnett &
Bhattacharjee (2005) p. 368, or Fitzpatrick (2015) p. 258; and for a
3-component plasma which allows the net current to be zero by the ion-ion (or
the electron-ion & electron-electron) acoustic instability, see Gary (1993) p.
44-55.

####  8.9.3 Ion-sound wave damping in homogeneous media: comparison of various
models

To further clarify our limitations, it is useful to explore the linear
properties of waves propagating along the ambient magnetic field (assumed to
be straight and aligned with the $z$-coordinate) in a homogeneous medium, in
regimes that range from the highly-collisional to the weakly-collisional ones.
Let us in particular consider the damping of a monochromatic ion-sound wave of
parallel wavenumber $k_{\|}$ in a proton-electron plasma where the electrons
are assumed to be cold. The latter assumption is not physically appropriate
because kinetic theory is not well-defined for cold electrons (see e.g.
discussion in Hunana et al. (2019a) p. 73) but it allows one to simplify the
presentation with a goal to describe the general behavior and not to provide
precise values of the damping rates. Four different models are compared in
Figure 2, all using the heuristic BGK collisional operator, which leads to
much simpler calculations for models with a distribution function expanded
around a bi-Maxwellian.

Figure 2: Normalized damping rate
$\zeta_{i}=\omega_{i}/(|k_{\parallel}|v_{\textrm{th}\parallel})$ for a
parallel propagating ion-sound wave as a function of
$k_{\parallel}\lambda_{\textrm{mfp}}$, where $\lambda_{\textrm{mfp}}$ is a
mean-free-path. Red line: Braginskii-type (isotropic) 13-moment model with
quasi-static stress-tensor and heat flux; Green line: Braginskii-type
13-moment model with evolution equations for the stress-tensor and heat flux;
Blue line: bi-Maxwellian 16-moment model with evolution equations for parallel
and perpendicular pressures and (gyrotropic) heat fluxes; Magenta line: Landau
fluid model with quasi-static heat fluxes of Snyder et al. (1997).

The x-axis shows $k_{\|}\lambda_{\textrm{mfp}}$ where
$\lambda_{\textrm{mfp}}=v_{{\textrm{th}}\|}/\nu$ is the ion mean-free-path and
$\nu$ is the collisional frequency, so that $k_{\|}\lambda_{\textrm{mfp}}\ll
1$ represents a highly-collisional regime and $k_{\|}\lambda_{\textrm{mfp}}\gg
1$ represents a weakly-collisional regime. The y-axis shows a damping rate as
an imaginary part of $\zeta=\omega/(|k_{\parallel}|v_{\textrm{th}\parallel})$.
The usual isotropic 13-moment model (green line) and the anisotropic 16-moment
model (blue line) with all the moments described by their time-dependent
(dynamical) evolution equations were discussed after equation (225). For the
parallel sound mode at the linear level considered here, the 13-moment model
is reduced to evolution equations for $\rho,u_{z},p,\Pi_{zz},q_{z}$ (we
consider the case where $\Pi_{zz}$ and $q_{z}$ are coupled) and the 16-moment
model reduces to evolution equations for
$\rho,u_{z},p_{\parallel},p_{\perp},q_{\parallel},q_{\perp}$ (we consider mean
equal pressures $p_{\parallel}^{(0)}=p_{\perp}^{(0)}$). Figure 2 shows that
these two models behave in a similar way: both reach a maximum damping rate
around $k_{\|}\lambda_{\textrm{mfp}}\sim 0.5-1$ and converge towards zero
damping rate in the collisionless regime (with only a small shift in
$k_{\|}\lambda_{\textrm{mfp}}$ between them). In contrast, the red line,
corresponding to the 13-moment model with the $\Pi_{zz}$ & $q_{z}$ taken in
the quasi-static approximation, shows that the damping rate does not reach a
maximum and instead continues to increase in a weakly-collisional regime, and
around $k_{\|}\lambda_{\textrm{mfp}}\sim 6.3$ the sound mode stops existing
(it becomes non-propagating with zero real frequency). This is consequence of
the quasi-static approximation for the stress-tensor $\Pi_{zz}\sim 1/\nu$,
which in the collisionless regime becomes unbounded (the parallel heat flux
$q_{z}\sim 1/\nu$ becomes unbounded as well, but this simply reflects an
isothermal behavior with no damping present). While a vanishing damping is
preferred against a quantity that blows up in a weakly-collisional regime, all
three models are technically incorrect, because Landau damping provides a
significant contribution for the damping rate as the plasma becomes weakly
collisional. To illustrate the importance of Landau damping, we display by a
magenta line the damping rate obtained with a Landau fluid model which
contains evolution equations for $\rho,u_{z},p_{\parallel},p_{\perp}$, but
where the quasi-static $q_{\parallel}$ & $q_{\perp}$ are given by the
collisionaly modified 3+1 closures of Snyder et al. (1997), i.e. their
equations (48)-(49), which for isotropic mean temperatures $T_{\parallel
a}^{(0)}=T_{\perp a}^{(0)}$ considered here, are equivalent to (our thermal
speed contains a factor of 2, which is not the case in that paper)

$\displaystyle q_{\parallel a}^{(1)}$ $\displaystyle=$
$\displaystyle-\,\frac{\frac{4}{3\pi-8}n_{a}^{(0)}v_{\textrm{th}\parallel
a}^{2}}{\bar{\nu}_{a}+\frac{2\sqrt{\pi}}{3\pi-8}v_{\textrm{th}\parallel
a}|k_{\parallel}|}ik_{\parallel}T_{\parallel a}^{(1)};\qquad q_{\perp
a}^{(1)}=-\,\frac{\frac{1}{2}n_{a}^{(0)}v_{\textrm{th}\parallel
a}^{2}}{\bar{\nu}_{a}+\frac{\sqrt{\pi}}{2}v_{\textrm{th}\parallel
a}|k_{\parallel}|}ik_{\parallel}T_{\perp a}^{(1)},$ (234)

where in general $\bar{\nu}_{a}=\sum_{b}\nu_{ab}$. Technically, closures (234)
are only applicable to a weakly-collisional regime because $q_{\parallel
a}\neq 3q_{\perp a}$ in the highly-collisional limit. In spite of this and the
additional difficulty associated with the cold electron limit considered here,
an interesting point is that the behavior of the damping rate is very close to
the prediction of the three other models in the highly-collisional regime,
while the damping rate converges to a constant value in the collisionless
case. This is in fact analogous to the case of the damping of a pure sound
wave in rarefied media, which was considered by Stubbe (1994) and Stubbe &
Sukhorukov (1999). In the former paper, the result of an experiment by Meyer &
Sessler (1957) (measuring the damping length of a sound wave of a given
frequency $\omega$, emitted at one end of the domain filled with a rarefied
neutral gas) are compared with various theoretical models. The results are
very similar to those presented here, and show in particular that the damping
is dominated by a nonlocal effect analogous to Landau damping when
$2\nu/\omega$ decreases below unity (see Fig. 6 and 7 of Stubbe (1994)). This
simple result on the damping of an ion-sound wave shows that, in a homogeneous
medium, a Braginskii-type model provides reasonable predictions as long as the
typical wavelength is larger than the mean free path, or equivalently, when
its frequency stays below the collision frequency. More sophisticated models
are needed in the weakly collisional case, that should retain new
contributions originating from a Landau-fluid closure.

####  8.9.4 Large gradients and large drifts

It is now of interest to consider inhomogeneous situations, where other
applicability conditions apply for the Braginskii-type models. In high energy-
density laser-produced plasmas, there are often situations relevant for
inertial confinement fusion experiments, where the typical electron mean-free-
path becomes of the order of the typical scale of electron temperature
gradients, or even larger. In this case the usual Braginskii formulas, used
for example for the Nernst effect (see e.g. Lancia et al. (2014)) become
invalid and have to be replaced by nonlocal expressions. In this context an
explicit nonlocal formula was proposed by Luciani et al. (1983) for the
electron thermal heat flux due to steep temperature gradients, offering an
improvement (in the one-dimensional case) for the Spitzer-Braginskii heat
flux, where one required proportionality constant is obtained by fitting from
Fokker-Planck simulations. Further extension to three dimensions was proposed
by Schurtz et al. (2000), but it is to be noted that this approach is not
appropriate in the very weakly collisional case, as for example in the Solar
corona when the density has significantly decreased.

Additional complications arise in a regime of weak collisionality. In space
physics, the collisionless heat flux is typically associated with the free-
streaming formula of Hollweg (1974, 1976)

$\vec{\boldsymbol{q}}_{e}^{\textrm{Hollweg}}=\frac{3}{2}p_{e}\boldsymbol{u}_{sw}\alpha,$
(235)

where one multiplies the thermal energy of one electron $(3/2)T_{e}$ (we take
$k_{b}=1$ in the entire paper) with the number density $n_{e}$ and the solar
wind speed $\boldsymbol{u}_{sw}$. The free “bugger factor” $\alpha$ as Hollweg
(1974) calls it, is dependent on a given form of an electron distribution
function where the tail had departed and run away. Note that the parallel
_frictional_ heat fluxes (i.e. due to _small_ differences in drifts
$\delta\boldsymbol{u}$) of Spitzer & Härm (1953) and Braginskii (1965) are
also independent of collisional frequencies, even though derived from
collisions, and up to the numerical values have the same form as (235). As a
side note, in the numerical model of Spitzer & Härm (1953) the frictional heat
flux is technically incorrect because it does not satisfy the Onsager
symmetry; see our Tables I.4 & I.2, which was already criticized by Balescu
(1988) (p. 268). Of course, in our usual fluid formalism a tail of a
distribution function can not suddenly depart. Even though our model contains
evolution equations for the perturbation of the 4th-order moment (i.e. a
“reduced kurtosis” which describes if a distribution is tail-heavy or tail-
light) and also for the 5th-order moment (sometimes called a hyper-skewness),
our distribution functions still have to remain sufficiently close to
Maxwellian. For the isotropic 5-moment model (i.e. strict Maxwellians), the
runaway effect is just represented through collisional contributions
$\boldsymbol{R}_{ab}$ and $Q_{ab}$ which decrease to zero for large drifts
(see equations (171)-(173) derived in Appendix G.3; see also Dreicer (1959);
Tanenbaum (1967); Burgers (1969); Schunk (1977); Balescu (1988)). We note that
for sufficiently large drifts between species various instabilities can
develop with a subsequent development of turbulence, which should restrict the
runaway effect long before relativistic effects. Importantly, it is unclear
how the heat flux collisional contributions
$\vec{\boldsymbol{Q}}^{(3)}_{ab}\,{}^{\prime}$ (and higher) would look like
for unrestricted drifts, because the collisional integrals seem exceedingly
complicated. Even if calculated, only the drifts between species would be
allowed to be unrestricted, the distribution of each species will have to be
restricted to remain close to Maxwellian. For the simplest CGL plasmas (i.e.
considering colliding strict bi-Maxwellians with no stress-tensors or heat
fluxes), the corresponding collisional integrals where numerically evaluated
for selected cases by Barakat & Schunk (1981). For a further particular case
of unrestricted drifts only _along_ the magnetic field and of Coulomb
collisions, Hellinger & Trávníček (2009) obtained exact analytic forms for the
collisional integrals (for $p_{\parallel}$ & $p_{\perp}$), which are however
expressed through a double hypergeometric functions. Judging from the two
papers above, a proper extension of our model to an anisotropic regime with
unrestricted drifts seems overly complicated. Other approach for the heat flux
modeling was presented by Canullo et al. (1996).

####  8.9.5 Comments on the positivity of the perturbed distribution function

An additional complication arises in a low-collisionality regime in the
presence of sufficiently strong large-scale gradients. Perturbations of the
distribution function considered in equation (1) might become so large, that
the corresponding model might become invalid. The distribution function around
which to expand is indeed not well-defined in this case. Strictly speaking, in
a weakly-collisional (or a collisionless) regime, one should abandon the
construction of fluid models derived from the Boltzmann equation, and perform
kinetic simulations by directly evolving the Boltzmann equation. Perhaps the
best example is a radially expanding flow, such as the solar corona with
emerging solar wind, where the spherical expansion creates strong large-scale
gradients and simultaneously drives the system towards a collisionless regime.
It seems that in this extreme case it might be indeed possible (but not with
certainty) that the underlying distribution function can even become negative,
$f_{a}<0$, which is of course unphysical. We expect that our 21 & 22-moment
models might fail in this particular situation, even if evolution equations
are retained, but as discussed below, we were unfortunately not able to reach
clear conclusion and further research is needed to clearly establish the areas
of validity.

The $f_{a}<0$ was criticized for example by Scudder (2021); Cranmer & Schiff
(2021) (and references therein), on an example of 8-moment models in a quasi-
static approximation. It is in fact questionable if the $f_{a}<0$ can be shown
in a quasi-static approximation. It is necessary to distinguish between two
different cases, depending if large-scale gradients are present or absent
during the transition into the low-collisionality regime. In the homogeneous
case, the situation is clear because one needs to describe the presence of
waves with frequencies $\omega$, and neglecting the time-derivative $d/dt$ in
the evolution equations automatically imposes requirement $\omega\ll\nu$, i.e.
the collisional frequencies $\nu$ must remain sufficiently large. In this
case, it is erroneous to simply take the quasi-static heat flux
$\vec{\boldsymbol{q}}_{a}\sim 1/\nu$, evaluate it for some arbitrarily small
$\nu$ and claim that $f_{a}<0$. Instead, it is necessary to retain the
evolution equations with $d\vec{\boldsymbol{q}}_{a}/dt$; see e.g. (41), (51)
or the coupled system (135)-(139), which precludes one from reaching direct
interpretation that $f_{a}<0$ (unless one calculates the eigenvector and shows
otherwise). The negativity of the distribution function may not take place,
and as a consequence, the procedure seems inadequate for disproving the moment
method of Grad in a homogeneous low-collisionality regime. The situation is
much less clear when large-scale gradients are present, as in the example of
the solar wind expansion. In that case, it is possible to argue that keeping
the evolution equations and solving an initial value problem might help only
temporarily, because the system eventually has to converge to some stationary
solution, which might show that $f_{a}<0$. Such a possibility seems to be
implied by the simple 1D radially expanding quasi-static models (see e.g.
Cranmer & Schiff (2021), and references therein). The quasi-static
approximation can however be questioned in this case as well, but from a
different perspective. Introducing a heat flux or a stress-tensor is analogous
to introducing a new degree of freedom into a system, and if this new degree
of freedom is not restricted in any way, it might of course yield an
unphysical system. In plasma physics, degrees of freedom are usually
restricted by associated instabilities that develop, which can not be revealed
in a quasi-static approximation (even if an instability is non-propagating).
Useful examples are the anisotropic CGL and 16-moment models described above.
Using a quasi-static approximation, one might erroneously conclude that the
temperature anisotropy can grow without bounds in these models, whereas
considering evolution equations reveals the firehose and mirror instabilities,
which can restrict the anisotropy. Similar situation might be applicable here,
where sufficiently large drifts (and possibly large heat fluxes and stress-
tensors) might cause various instabilities and also development of turbulence,
but further clarifications are needed if our fluid models contain some of
these instabilities, especially considering that our collisional contributions
are valid only when differences in drifts between species are much smaller
than their thermal velocities. In this regard, it is not clear if it is
appropriate to neglect the Alfvénic fluctuations in the radially expanding
models. Finally, it is also not clear if it is physically meaningful to show
$f_{a}<0$ by skipping the stress-tensor in the expansions of Grad (which is a
2nd-order moment before the 3rd-order heat flux moment), because its
contributions to a total $f_{a}$ might be significant. For a sufficient proof
that the $f_{a}$ can become negative, it might be necessary to consider at
least the 13-moment model, where both stress-tensors and heat fluxes are
retained.

###  8.10 Conclusions

We have discussed various generalizations of the 21-moment model of Braginskii
(1958, 1965): 1) We have presented a multi-fluid formulation for arbitrary
masses $m_{a}$ & $m_{b}$ and arbitrary temperatures $T_{a}$ & $T_{b}$. 2) All
the fluid moments are described by their evolution equations, whose left hand
sides are given in a fully non-linear form. 3) Formulation with evolution
equations has an important consequence that the model does not become
divergent (unbounded) if a regime of low-collisionality is encountered. 4) For
a one ion-electron plasma we have provided all the Braginskii transport
coefficients in a fully analytic form for a general ion charge $Z_{i}$ (and
arbitrary strength of magnetic field). 5) We have also provided fully analytic
higher-order transport coefficients (for $\bar{\bar{\boldsymbol{\Pi}}}^{(4)}$
and $\vec{\boldsymbol{X}}^{(5)}$), which are not typically given. 6) All the
electron coefficients were further generalized to multi-ion plasmas. 7) We
have considered coupling between viscosity-tensors and heat fluxes, where a
heat flux enters a viscosity-tensor and a viscosity-tensor enters a heat flux.
As a consequence, we have introduced new higher-order physical effects even
for the simplest case of unmagnetized one ion-electron plasma of Spitzer &
Härm (1953). For example, the electron rate of strain tensor
$\bar{\bar{\boldsymbol{W}}}_{e}$ enters the electron heat fluxes even
linearly, and thus, it subsequently linearly enters the momentum exchange
rates; see equation (127). 8) We have formulated the 22-moment model which is
a natural extension of the 21-moment model, where one takes into account fully
contracted scalar perturbations $\widetilde{X}^{(4)}_{a}$ entering the
decomposition of the 4-th order moment $X_{ijkl}^{a(4)}$; see equation (129).
Collisional contributions for this model with arbitrary masses and
temperatures are given in Section 7.1 and supplement those given in Section
2.1 for the 21-moment model. Interestingly, scalar perturbations
$\widetilde{X}^{(4)}_{a}$ modify the energy exchange rates, see equation (140)
or (177). In the quasi-static approximation, scalar perturbations
$\widetilde{X}^{(4)}_{a}$ can be written as divergence of heat flux vectors
with their own heat conductivities; see for example solutions for a one ion-
electron plasma with the ion heat conductivities (149) and the electron heat
conductivities (162). These corrections remain small in the highly-collisional
regime, but might become significant at small wavelengths and/or at large
frequencies.

Our model can be useful for direct numerical simulations, as well as for quick
calculation of transport coefficients in a quasi-static approximation. We
provide three examples for coupling between two species. Thermal
conductivities and viscosities for unmagnetized proton-electron plasma
(without $m_{e}/m_{p}$ expansions) were presented in Section 8.8, and two
examples for proton-alpha particles and deuterium-tritium were moved to
Appendix N. Our model can also be useful from an observational perspective.
For example, the parallel thermal heat flux $\vec{\boldsymbol{q}}_{e}$ of
Braginskii (1965) and Spitzer & Härm (1953) (they differ only by 3.16 vs 3.20
factors rounded as 3.2) is sometimes analyzed in observational studies; see
e.g. Salem et al. (2003); Bale et al. (2013); Halekas et al. (2021) and
Verscharen et al. (2019) (page 61). It is also measured in (exospheric)
kinetic numerical simulations (Landi et al., 2014). Our model suggests that it
would be beneficial to analyze both parallel heat fluxes, which for $Z_{i}=1$
read

$\vec{\boldsymbol{q}}_{e}=\frac{\vec{\boldsymbol{X}}^{(3)}_{e}}{2}=-3.2\frac{p_{e}}{m_{e}\nu_{ei}}\nabla
T_{e};\qquad\vec{\boldsymbol{X}}^{(5)}_{e}=-110.7\frac{p_{e}^{2}}{\rho_{e}m_{e}\nu_{ei}}\nabla
T_{e},$ (236)

and which can be analyzed with the same techniques. For long parallel mean-
free-paths (low collisionality regime), both heat fluxes naturally have to
become non-local and independent of the mean-free-path. Our limitations are
described in Section 8.9, and “flattening/saturation” of heat fluxes due to
the runaway effect and Landau damping is not captured in our model. Our model
is aimed at the highly-collisional regime and in the low-collisionality regime
our heat fluxes are just described by their evolution equations, where the
collisional right hand sides are small. Nevertheless, it would be interesting
to see if in observational studies or kinetic simulations the
$\vec{\boldsymbol{X}}^{(5)}_{e}$ could be described by a free-streaming
formula similar to the one of Hollweg (1974, 1976), in a form
$\vec{\boldsymbol{X}}^{(5)}_{e}=(3/2)(p_{e}^{2}/\rho_{e})\boldsymbol{u}_{sw}\alpha_{5}$,
where the “bugger factor” $\alpha_{5}$ has to be determined from a given form
of a distribution function, or if such a concept does not apply for
$\vec{\boldsymbol{X}}^{(5)}_{e}$. As a side note, concerning collisionless
heat fluxes for plasmas where spherical expansion and large drifts are not
present and Landau damping dominates, our model actually implies that a
correct interpretation should not be that Landau damping diminishes/saturates
the heat flux in a low-collisionality regime. The correct interpretation is
that Landau damping creates the collisionless heat flux. Collisionless Landau
fluid closures for quasi-static parallel scalar $X^{(5)}_{\parallel}$ can be
found in Hunana et al. (2019a) (p. 84). In addition to (236), it might be also
useful to analyze the scalar perturbation, which for $Z_{i}=1$ reads

$\widetilde{X}^{(4)}_{e}=+83.8\frac{p_{e}^{2}}{\nu_{ei}^{2}\rho_{e}m_{e}}\nabla^{2}T_{e}.$
(237)

Our multi-fluid model might be also useful for modeling of enrichment of minor
ion abundancies in stellar atmospheres, because of the very precise thermal
force (thermal diffusion). Let us summarize the thermal force description in
three major models: the model of Burgers (1969)-Schunk (1977), the model of
Killie et al. (2004), and our model. Of course, all three models are
formulated as a general multi-fluid, but for a simplicity of the discussion
let us simplify and compare only thermal forces given by

$\displaystyle\textrm{Burgers-Schunk:}\qquad\boldsymbol{R}_{e}^{T}$
$\displaystyle=$
$\displaystyle+\frac{3}{5}\frac{\rho_{e}}{p_{e}}\nu_{ei}\vec{\boldsymbol{q}}_{e};$
(238) $\displaystyle\textrm{Killie et al.:}\qquad\boldsymbol{R}_{e}^{T}$
$\displaystyle=$
$\displaystyle+\frac{6}{35}\frac{\rho_{e}}{p_{e}}\nu_{ei}\vec{\boldsymbol{q}}_{e};$
(239) $\displaystyle\textrm{present paper:}\qquad\boldsymbol{R}_{e}^{T}$
$\displaystyle=$
$\displaystyle+\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}.$
(240)

Note that the viscosity-tensors are not required to describe the thermal
force, and focusing only at the heat fluxes, instead of the 13-moment model of
Burgers-Schunk, one can consider only the 8-moment model. Similarly, instead
of our 21 & 22-moment models, one can consider only the 11-moment model. In
general, the parallel thermal heat flux is given by
$\vec{\boldsymbol{q}}_{e}=-\gamma_{0}p_{e}/(m_{e}\nu_{ei})\nabla T_{e}$ and
the resulting parallel thermal force by
$\boldsymbol{R}_{e}^{T}=-\beta_{0}n_{e}\nabla T_{e}$, with coefficients
$\gamma_{0}$ and $\beta_{0}$. From the work of Spitzer & Härm (1953), for
$Z_{i}=1$ the correct coefficient of thermal conductivity is
$\gamma_{0}=3.203$ and the correct coefficient of thermal force is
$\beta_{0}=0.703$. The model of Burgers-Schunk (238) has thermal conductivity
$\gamma_{0}=1.34$, and with that value it describes the thermal force actually
quite accurately, yielding $\beta_{0}=0.804$ (for other $Z_{i}$ values see
comparison in Table I.2 in Appendix I). However, a problem arises if one uses
the correct value of thermal conductivity $\gamma_{0}=3.2$ in the expression
(238), which overestimates the thermal force. Killie et al. (2004) developed a
different 8-moment model, where the expansion is done differently than in
equation (1), with the goal to improve the heat flux and the thermal force of
Burgers-Schunk. The model is described in Appendix I.2. For $Z_{i}=1$, its
heat flux value is $\gamma_{0}=3.92$, which greatly improves the model of
Burgers-Schunk and for that value it also improves the thermal force, yielding
$\beta_{0}=0.672$. Additionally, now one can use the correct $\gamma_{0}=3.2$
value in expression (239) and the thermal force will be roughly correct (and
7/2 times smaller than Burgers-Schunk). However, as we point out in Appendix I
(see Table I.4), the model of Killie et al. (2004) breaks the Onsager symmetry
between the frictional heat flux and the thermal force. The numerical model of
Spitzer & Härm (1953) also does not satisfy the Onsager symmetry and its
frictional heat flux is technically incorrect, even though in this case the
discrepancies are small. Our model satisfies the Onsager symmetry, it has
thermal conductivity $\gamma_{0}=3.1616$ and thermal force $\beta_{0}=0.711$
(the same as Braginskii). In summary, our multi-fluid model has a very precise
thermal force (240) with precision equal to Braginskii (1965), and we thus
offer an improvement to the multi-fluid models of Burgers (1969)-Schunk (1977)
and Killie et al. (2004).

##  9 Acknowledgments

This work was supported by the European Research Council in the frame of the
Consolidating Grant ERC-2017-CoG771310-PI2FA “Partial Ionisation: Two-Fluid
Approach”, led by Elena Khomenko. Anna Tenerani acknowledges support of the
NASA Heliophysics Supporting Research Grant #80NSSC18K1211. We acknowledge
support of the NSF EPSCoR RII-Track-1 Cooperative Agreement No. OIA-1655280
“Connecting the Plasma Universe to Plasma Technology in Alabama”, led by Gary
P. Zank. Gary M. Webb was funded in part by NASA grant 80NSSC19K0075.

##  Appendix A General evolution equations

We consider Boltzmann equation (in CGS units)

$\frac{\partial f_{a}}{\partial t}+\boldsymbol{v}\cdot\nabla
f_{a}+\frac{eZ_{a}}{m_{a}}(\boldsymbol{E}+\frac{1}{c}\boldsymbol{v}\times\boldsymbol{B})\cdot\nabla_{v}f_{a}=C(f_{a}),$
(A1)

where “a” is a species index and $C(f_{a})=\sum_{b}C_{ab}(f_{a},f_{b})$ is the
Landau collisional operator, so equation (A1) is called the Landau equation.
One defines the usual number density $n_{a}=\int f_{a}d^{3}v$; density
$\rho_{a}=m_{a}n_{a}$, fluid velocity
$\boldsymbol{u}_{a}=(1/n_{a})\int\boldsymbol{v}f_{a}d^{3}v$ and fluctuating
velocity $\boldsymbol{c}_{a}=\boldsymbol{v}-\boldsymbol{u}_{a}$, and further
defines the pressure tensor $\bar{\bar{\boldsymbol{p}}}_{a}$, heat flux tensor
$\bar{\bar{\boldsymbol{q}}}_{a}$, 4th-order moment
$\bar{\bar{\boldsymbol{r}}}_{a}$, and 5th-order & 6th-order moments
$\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}$, $\bar{\bar{\boldsymbol{X}}}^{(6)}_{a}$
according to

$\displaystyle\bar{\bar{\boldsymbol{p}}}_{a}$ $\displaystyle=$ $\displaystyle
m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}f_{a}d^{3}v;\qquad\bar{\bar{\boldsymbol{q}}}_{a}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}f_{a}d^{3}v;\qquad\bar{\bar{\boldsymbol{r}}}_{a}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}f_{a}d^{3}v;$
(A2) $\displaystyle\bar{\bar{\boldsymbol{X}}}_{a}^{(5)}$ $\displaystyle=$
$\displaystyle
m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}f_{a}d^{3}v;\qquad\bar{\bar{\boldsymbol{X}}}_{a}^{(6)}=m_{a}\int\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}\boldsymbol{c}_{a}f_{a}d^{3}v.$
(A3)

Writing of the tensor product $\otimes$ is suppressed everywhere and
$\boldsymbol{c}_{a}\boldsymbol{c}_{a}=\boldsymbol{c}_{a}\otimes\boldsymbol{c}_{a}$.
For complicated fluid models the species index ${}^{\prime}a^{\prime}$ often
blurs the clarity of the tensor algebra, and thus in the vector notation (A2)
we emphasize tensors of second-rank and above with the double overbar symbol.
Sometimes we move the index ${}^{\prime}a^{\prime}$ freely up and down (which
here does not represent any mathematical operation), and in the index notation
the index ${}^{\prime}a^{\prime}$ is often dropped completely, so for example
$p_{ij}^{a}=m_{a}\int c_{i}^{a}c_{j}^{a}f_{a}d^{3}v$ and $p_{ij}=m\int
c_{i}c_{j}fd^{3}v$ are equivalent. The Einstein summation convention does not
apply for the species index ${}^{\prime}a^{\prime}$, and summations over other
particle species are written down explicitly. The divergence is defined
through the first index
$(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})_{j}=\partial_{i}p^{a}_{ij}$.

Here we do not consider ionization and recombination processes and the Landau
collisional operator conserves the number of particles $\int C(f_{a})d^{3}v=0$
for each species. One defines a unit vector in the direction of magnetic field
$\boldsymbol{\hat{b}}=\boldsymbol{B}/|\boldsymbol{B}|$, cyclotron frequency
$\Omega_{a}=eZ_{a}|\boldsymbol{B}|/(m_{a}c)$ and convective derivative
$d_{a}/dt=\partial/\partial t+\boldsymbol{u}_{a}\cdot\nabla$. It is also
useful to define a symmetric operator ${}^{\prime}S^{\prime}$, which acts on a
matrix as $A_{ij}^{S}=A_{ij}+A_{ji}$ and on a tensor of 3rd-rank as
$A_{ijk}^{S}=A_{ijk}+A_{jki}+A_{kij}$, i.e. it cycles around all indices. We
often use operator trace Tr and unit matrix $\bar{\bar{\boldsymbol{I}}}$,
where
$\textrm{Tr}\bar{\bar{\boldsymbol{A}}}=\bar{\bar{\boldsymbol{I}}}:\bar{\bar{\boldsymbol{A}}}$,
and operator ’:’ represents double contraction. We also use
$\bar{\bar{\boldsymbol{I}}}_{\perp}=\bar{\bar{\boldsymbol{I}}}-\boldsymbol{\hat{b}}\boldsymbol{\hat{b}}$.

To derive the model of Braginskii (1965) with the moment method of Grad, it is
necessary to consider the evolution equation for the 5th-order moment
$\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}$ and perform a closure at
$\bar{\bar{\boldsymbol{X}}}^{(6)}_{a}$. Integrating (A1) over velocity space
yields the the following hierarchy of evolution equations

$\displaystyle\frac{\partial n_{a}}{\partial
t}+\nabla\cdot(n_{a}\boldsymbol{u}_{a})=0;$ (A4)
$\displaystyle\frac{\partial\boldsymbol{u}_{a}}{\partial
t}+\boldsymbol{u}_{a}\cdot\nabla\boldsymbol{u}_{a}+\frac{1}{\rho_{a}}\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a}-\frac{eZ_{a}}{m_{a}}\Big{(}\boldsymbol{E}+\frac{1}{c}\boldsymbol{u}_{a}\times\boldsymbol{B}\Big{)}=\frac{\boldsymbol{R}_{a}}{\rho_{a}};$
(A5) $\displaystyle\frac{\partial\bar{\bar{\boldsymbol{p}}}_{a}}{\partial
t}+\nabla\cdot\big{(}\bar{\bar{\boldsymbol{q}}}_{a}+\boldsymbol{u}_{a}\bar{\bar{\boldsymbol{p}}}_{a}\big{)}+\Big{[}\bar{\bar{\boldsymbol{p}}}_{a}\cdot\nabla\boldsymbol{u}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{p}}}_{a}\Big{]}^{S}=\bar{\bar{\boldsymbol{Q}}}^{(2)}_{a};$
(A6) $\displaystyle\frac{\partial\bar{\bar{\boldsymbol{q}}}_{a}}{\partial
t}+\nabla\cdot\big{(}\bar{\bar{\boldsymbol{r}}}_{a}+\boldsymbol{u}_{a}\bar{\bar{\boldsymbol{q}}}_{a}\big{)}+\Big{[}\bar{\bar{\boldsymbol{q}}}_{a}\cdot\nabla\boldsymbol{u}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{q}}}_{a}-\frac{1}{\rho_{a}}(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})\bar{\bar{\boldsymbol{p}}}_{a}\Big{]}^{S}=\bar{\bar{\boldsymbol{Q}}}^{(3)}_{a}-\frac{1}{\rho_{a}}\Big{[}\boldsymbol{R}_{a}\bar{\bar{\boldsymbol{p}}}_{a}\Big{]}^{S};$
(A7) $\displaystyle\frac{\partial}{\partial
t}\bar{\bar{\boldsymbol{r}}}_{a}+\nabla\cdot\big{(}\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}+\boldsymbol{u}_{a}\bar{\bar{\boldsymbol{r}}}_{a}\big{)}+\Big{[}\bar{\bar{\boldsymbol{r}}}_{a}\cdot\nabla\boldsymbol{u}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{r}}}_{a}-\frac{1}{\rho_{a}}\big{(}\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a}\big{)}\bar{\bar{\boldsymbol{q}}}_{a}\Big{]}^{S}=\bar{\bar{\boldsymbol{Q}}}^{(4)}_{a}-\frac{1}{\rho_{a}}\big{[}\boldsymbol{R}_{a}\bar{\bar{\boldsymbol{q}}}_{a}\big{]}^{S};$
(A8) $\displaystyle\frac{\partial}{\partial
t}\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}+\nabla\cdot\big{(}\bar{\bar{\boldsymbol{X}}}^{(6)}_{a}+\boldsymbol{u}_{a}\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}\big{)}+\Big{[}\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}\cdot\nabla\boldsymbol{u}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{X}}}^{(5)}_{a}-\frac{1}{\rho_{a}}\big{(}\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a}\big{)}\bar{\bar{\boldsymbol{r}}}_{a}\Big{]}^{S}=\bar{\bar{\boldsymbol{Q}}}^{(5)}_{a}-\frac{1}{\rho_{a}}\big{[}\boldsymbol{R}_{a}\bar{\bar{\boldsymbol{r}}}_{a}\big{]}^{S},$
(A9)

where the collisional contributions on the right hand sides are given by (5).
It is also possible to define a general n-th-order moment
$\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}$ and collisional contributions
$\bar{\bar{\boldsymbol{Q}}}^{(n)}_{a}$

$X^{(n)}_{r_{1}r_{2}\ldots r_{n}}=m\int c_{r_{1}}c_{r_{2}}\ldots
c_{r_{n}}fd^{3}v;\qquad Q^{(n)}_{r_{1}r_{2}\ldots r_{n}}=m\int
c_{r_{1}}c_{r_{2}}\ldots c_{r_{n}}C(f)d^{3}v,$ (A10)

together with symmetric operator ’S’ that cycles around all of its indices

$\big{[}X^{(n)}\big{]}^{S}_{r_{1}r_{2}r_{3}\ldots
r_{n}}=X^{(n)}_{r_{1}r_{2}r_{3}\ldots r_{n}}+X^{(n)}_{r_{2}r_{3}\ldots
r_{n}r_{1}}+X^{(n)}_{r_{3}\ldots
r_{n}r_{1}r_{2}}+\cdots\cdots+X^{(n)}_{r_{n}r_{1}r_{2}r_{3}\ldots r_{n-1}},$
(A11)

(so that it contains “n” terms) and derive the following evolution equation
for $\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}$

$\displaystyle\frac{\partial}{\partial
t}\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}+\nabla\cdot\big{(}\bar{\bar{\boldsymbol{X}}}^{(n+1)}_{a}+\boldsymbol{u}_{a}\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}\big{)}+\Big{[}\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}\cdot\nabla\boldsymbol{u}_{a}+\Omega_{a}\boldsymbol{\hat{b}}\times\bar{\bar{\boldsymbol{X}}}^{(n)}_{a}-\frac{1}{\rho_{a}}(\nabla\cdot\bar{\bar{\boldsymbol{p}}}_{a})\bar{\bar{\boldsymbol{X}}}^{(n-1)}_{a}\Big{]}^{S}$
$\displaystyle=\bar{\bar{\boldsymbol{Q}}}^{(n)}_{a}-\frac{1}{\rho_{a}}\Big{[}\boldsymbol{R}_{a}\bar{\bar{\boldsymbol{X}}}^{(n-1)}_{a}\Big{]}^{S},$
(A12)

valid for $n\geq 2$. The left hand side of (A12) is equal to the collisionless
equation (12.16) of Hunana et al. (2019b). Evolution equations (A6)-(A9) then
can be obtained easily by evaluation of (A12). Note that definition (A10)
yields $\bar{\bar{\boldsymbol{X}}}^{(2)}=\bar{\bar{\boldsymbol{p}}}$,
$\bar{\bar{\boldsymbol{X}}}^{(3)}=\bar{\bar{\boldsymbol{q}}}$,