# Generalized Fluid Models of the Braginskii Type

P. Hunana, T. Passot, E. Khomenko, D. Martínez-Gómez, M. Collados, A.

Tenerani, G. Zank, Y. Maneva, M. Goldstein, G. Webb

► To cite this version:

P. Hunana, T. Passot, E. Khomenko, D. Martínez-Gómez, M. Collados, et al.. Generalized Fluid Models of the Braginskii Type. The Astrophysical Journal Supplement Series, 2022, 260 (2), pp.26. 10.3847/1538-4365/ac5044 . hal-03698172

HAL Id: hal-03698172

https://hal.science/hal-03698172

Submitted on 18 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Generalized Fluid Models of the Braginskii Type

P. Hunana<fnref n="1" />,<fnref n="2" /> ID, T. Passot<fnref n="3" /> ID, E. Khomenko<fnref n="1" />,<fnref n="2" /> ID, D. Martínez-Gómez<fnref n="1" />,<fnref n="2" /> ID, M. Collados<fnref n="1" />,<fnref n="2" /> ID, A. Tenerani<fnref n="4" /> ID, G. P. Zank<fnref n="5" />,<fnref n="6" /> ID,
Y. Maneva<fnref n="7" /> ID, M. L. Goldstein<fnref n="8" /> ID, and G. M. Webb<fnref n="5" /> ID

Received 2021 October 18; revised 2022 January 27; accepted 2022 January 29; published 2022 June 17

## Abstract

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

Unified Astronomy Thesaurus concepts: Collision processes (2065); Plasma physics (2089); Space plasmas (1544)

## 1. Introduction

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

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

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Instituto de Astrofísica de Canarias (IAC), La Laguna, Tenerife, E-38205, Spain; peter.hunana@gmail.com</footnote>

<footnote n="2">² Universidad de La Laguna, La Laguna, Tenerife, E-38206, Spain</footnote>

<footnote n="3">³ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, F-06304 Nice cedex 4, France</footnote>

<footnote n="4">⁴ Department of Physics, University of Texas, Austin, TX 78712, USA</footnote>

<footnote n="5">⁵ Center for Space Plasma and Aeronomic Research (CSPAR), University of Alabama, Huntsville, AL 35805, USA</footnote>

<footnote n="6">⁶ Department of Space Science, University of Alabama, Huntsville, AL 35899, USA</footnote>

<footnote n="7">⁷ Royal Observatory of Belgium, B-1180 Brussels, Belgium</footnote>

<footnote n="8">⁸ Space Science Institute, Boulder, CO 80301, USA</footnote>
function as a series

$$f_a = f_a^{(0)}(1 + \chi_a); \quad \chi_a = \sum_{n=1}^{N} [h_{ij}^{(2n)} H_{ij}^{(2n)} + h_i^{(2n+1)} H_i^{(2n+1)}], \quad (1)$$

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

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

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

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

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

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

$$f_a = f_a^{(0)}(1 + \chi_a); \quad \chi_a = \sum_{n=1}^{N} [h_{ij}^{(2n)} H_{ij}^{(2n)} + h_i^{(2n+1)} H_i^{(2n+1)} + h^{(2n)} H^{(2n)}], \quad (2)$$

where the scalar hermite moments $h^{(2n)}$ can be viewed as fully contracted (scalar) perturbations of fluid moments. The lowest-order moment $h^{(2)} = 0$ and all of the higher-order ones are generally nonzero. Thus, prescribing $N=1$ still yields the 13-moment model, but prescribing $N=2$ now yields the 22-moment model. This model is a natural extension of the Braginskii model, because it takes into account the fully contracted perturbations $\tilde{X}_a^{(4)} = m_a \int |c_a|^4 (f_a - f_a^{(0)}) d^3v$ of the fourth-order fluid moment. Accounting for the scalar perturbations according to (2), for $N \ge 1$ an N-Laguerre model then represents a $(4+9N)$-moment model. Another possibility for writing Equation (2) is to separate the matrices $\sum_{n=1}^N h_{ij}^{(2n)} H_{ij}^{(2n)}$, and to write the sum for the vectors and scalars from $n=0$, with an imposed requirement that $h^{(0)} = 0$, $h^{(2)} = 0$, and $h_i^{(1)} = 0$ (where the first one is nontrivial). This is the choice of Balescu (1988, p. 174), for example, in his Equations (3.11) and (3.16).

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

The entire paper is separated into eight sections and 14 appendices. The main paper summarizes the obtained results, while the appendices provide the detailed calculations.

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

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

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

In Section 5, we discuss the coupling between viscosity tensors and heat fluxes. We provide evolution equations in the semilinear approximation, where this coupling is retained, and we introduce a technique for splitting the moments into their first and second orders.

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

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

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

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

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

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

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

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

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

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

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

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

Appendix J calculates in detail the 10-moment multifluid model of Burgers (1969)-Schunk (1977), where the simplest viscosity tensor is present. It is shown that in the limit of a strong magnetic field, the perpendicular viscosities and gyroviscosities
$\eta_1, \eta_2, \eta_3, \eta_4$ match the Braginskii model exactly (for both the ion and electron species), and only the parallel viscosities $\eta_0$ are different.

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

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

## 2. Multifluid Generalization of Braginskii (21-moment Model)

Our model is formulated with heat flux vectors

$$X_a^{(3)} = m_a \int \mathbf{c}_a |\mathbf{c}_a|^2 f_a d^3v = 2\mathbf{q}_a; \quad X_a^{(5)} = m_a \int \mathbf{c}_a |\mathbf{c}_a|^4 f_a d^3v, \quad (3)$$

and traceless viscosity tensors

$$\bar{\bar{\Pi}}_a^{(2)} = m_a \int \left( \mathbf{c}_a \mathbf{c}_a - \frac{\bar{\bar{I}}}{3} |\mathbf{c}_a|^2 \right) f_a d^3v; \quad \bar{\bar{\Pi}}_a^{(4)} = m_a \int \left( \mathbf{c}_a \mathbf{c}_a - \frac{\bar{\bar{I}}}{3} |\mathbf{c}_a|^2 \right) |\mathbf{c}_a|^2 f_a d^3v, \quad (4)$$

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

We first present a formulation with a general (unspecified) collisional operator $C(f_a)$. We define the (tensorial) collisional contributions

$$\begin{align}
\mathbf{R}_a &= m_a \int \mathbf{v} C(f_a) d^3v; & Q_a &= \frac{m_a}{2} \int |\mathbf{c}_a|^2 C(f_a) d^3v; \nonumber \\
\bar{\bar{Q}}_a^{(2)} &= m_a \int \mathbf{c}_a \mathbf{c}_a C(f_a) d^3v; & \bar{\bar{Q}}_a^{(3)} &= m_a \int \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a C(f_a) d^3v; \nonumber \\
\bar{\bar{Q}}_a^{(4)} &= m_a \int \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a C(f_a) d^3v; & \bar{\bar{Q}}_a^{(5)} &= m_a \int \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a C(f_a) d^3v, \tag{5}
\end{align}$$

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

$$\frac{d}{dt} n_a + n_a \nabla \cdot \mathbf{u}_a = 0; \quad (6)$$

$$\frac{d}{dt} \mathbf{u}_a + \frac{1}{\rho_a} \nabla \cdot \bar{\bar{\mathbf{p}}}_a - \mathbf{G} - \frac{eZ_a}{m_a} \left( \mathbf{E} + \frac{1}{c} \mathbf{u}_a \times \mathbf{B} \right) = \mathbf{R}_a; \quad (7)$$

$$\frac{d}{dt} p_a + \frac{5}{3} p_a \nabla \cdot \mathbf{u}_a + \frac{2}{3} \nabla \cdot \mathbf{q}_a + \frac{2}{3} \bar{\bar{\Pi}}_a^{(2)} : (\nabla \mathbf{u}_a) = \frac{2}{3} Q_a, \quad (8)$$
and are accompanied by evolution equations for the stress tensors and heat flux vectors:

$$ 
\begin{aligned}
& \frac{d_a \bar{\Pi}_a^{(2)}}{dt} + \bar{\Pi}_a^{(2)} \nabla \cdot \mathbf{u}_a + \Omega_a (\hat{\mathbf{b}} \times \bar{\Pi}_a^{(2)})^S + (\bar{\Pi}_a^{(2)} \cdot \nabla \mathbf{u}_a)^S - \frac{2}{3} \bar{I} (\bar{\Pi}_a^{(2)} : \nabla \mathbf{u}_a) \\
& \qquad + \frac{2}{5} \left[ (\nabla \mathbf{q}_a)^S - \frac{2}{3} \bar{I} \nabla \cdot \mathbf{q}_a \right] + p_a \bar{W}_a = \bar{Q}_a^{(2)\prime} \equiv \bar{Q}_a^{(2)} - \frac{\bar{I}}{3} \mathrm{Tr} \bar{Q}_a^{(2)}; \tag{9} \\
& \frac{d_a \mathbf{q}_a}{dt} + \frac{7}{5} \mathbf{q}_a \nabla \cdot \mathbf{u}_a + \frac{7}{5} \mathbf{q}_a \cdot \nabla \mathbf{u}_a + \frac{2}{5} (\nabla \mathbf{u}_a) \cdot \mathbf{q}_a + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) \\
& \qquad + \frac{1}{2} \nabla \cdot \bar{\Pi}_a^{(4)} - \frac{5}{2} \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(2)} - \frac{1}{\rho_a} (\nabla \cdot \bar{\rho}_a) \cdot \bar{\Pi}_a^{(2)} \\
& \qquad = \mathbf{Q}_a^{(3)\prime} \equiv \frac{1}{2} \mathrm{Tr} \bar{Q}_a^{(3)} - \frac{5}{2} \frac{p_a}{\rho_a} \mathbf{R}_a - \frac{1}{\rho_a} \mathbf{R}_a \cdot \bar{\Pi}_a^{(2)};
\end{aligned}
\tag{10}
$$

$$ 
\begin{aligned}
& \frac{d_a \bar{\Pi}_a^{(4)}}{dt} + \frac{1}{5} \left[ (\nabla \mathbf{X}_a^{(5)})^S - \frac{2}{3} \bar{I} (\nabla \cdot \mathbf{X}_a^{(5)}) \right] + \frac{9}{7} (\nabla \cdot \mathbf{u}_a) \bar{\Pi}_a^{(4)} + \frac{9}{7} (\bar{\Pi}_a^{(4)} \cdot \nabla \mathbf{u}_a)^S \\
& \qquad + \frac{2}{7} ((\nabla \mathbf{u}_a) \cdot \bar{\Pi}_a^{(4)})^S - \frac{22}{21} \bar{I} (\bar{\Pi}_a^{(4)} : \nabla \mathbf{u}_a) + \Omega_a (\hat{\mathbf{b}} \times \bar{\Pi}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{W}_a \\
& \qquad - \frac{14}{5 \rho_a} \left[ ((\nabla \cdot \bar{\rho}_a) \mathbf{q}_a)^S - \frac{2}{3} \bar{I} (\nabla \cdot \bar{\rho}_a) \cdot \mathbf{q}_a \right] \\
& = \bar{Q}_a^{(4)\prime} \equiv \mathrm{Tr} \bar{Q}_a^{(4)} - \frac{\bar{I}}{3} \mathrm{Tr} \bar{Q}_a^{(4)} - \frac{14}{5 \rho_a} \left[ (\mathbf{R}_a \mathbf{q}_a)^S - \frac{2}{3} \bar{I} (\mathbf{R}_a \cdot \mathbf{q}_a) \right]; \tag{11} \\
& \frac{d_a \mathbf{X}_a^{(5)}}{dt} + \nabla \cdot \bar{\Pi}_a^{(6)} + \frac{9}{5} \mathbf{X}_a^{(5)} (\nabla \cdot \mathbf{u}_a) + \frac{9}{5} \mathbf{X}_a^{(5)} \cdot \nabla \mathbf{u}_a + \frac{4}{5} (\nabla \mathbf{u}_a) \cdot \mathbf{X}_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} \\
& \qquad + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) - 35 \frac{p_a^2}{\rho_a^2} \nabla \cdot \bar{\Pi}_a^{(2)} - \frac{4}{\rho_a} (\nabla \cdot \bar{\rho}_a) \cdot \bar{\Pi}_a^{(4)} \\
& = \mathbf{Q}_a^{(5)\prime} \equiv \mathrm{Tr} \mathrm{Tr} \bar{Q}_a^{(5)} - 35 \frac{p_a^2}{\rho_a^2} \mathbf{R}_a - \frac{4}{\rho_a} \mathbf{R}_a \cdot \bar{\Pi}_a^{(4)}.
\end{aligned}
\tag{12}
$$

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

$$ \bar{\Pi}_a^{(6)} = m_a \int \left( c_a c_a - \frac{\bar{I}}{3} |c_a|^2 \right) |c_a|^4 f_a d^3 v = 18 \frac{p_a}{\rho_a} \bar{\Pi}_a^{(4)} - 63 \frac{p_a^2}{\rho_a^2} \bar{\Pi}_a^{(2)}. \tag{13} $$

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

## 2.1. Collisional Contributions (Arbitrary Masses and Temperatures)

We use the Landau collisional operator. All of the collisional contributions are evaluated in a semilinear approximation, with an assumption that the differences in the drift velocities $\mathbf{u}_b - \mathbf{u}_a$ are small with respect to thermal velocities. All of the nonlinear quantities, such as $\mathbf{q}_a \cdot (\mathbf{u}_b - \mathbf{u}_a)$, including $|\mathbf{u}_b - \mathbf{u}_a|^2$, are thus neglected in the multifluid description, which is consistent with models of Burgers (1969) and Schunk (1977). For energy conservation and a particular case of a one ion-electron plasma, see Section 8.1. The wording “semilinear” just means that expressions containing pressures and densities, such as $(p_a/\rho_a)\mathbf{q}_a$, are retained and not fully linearized with their mean pressure/density values. However, for example, the last terms of the collisional contributions in Equations (10)–(12) proportional to $\mathbf{R}_a \mathbf{q}_a$, $\mathbf{R}_a \cdot \bar{\Pi}_a^{(2)}$, and $\mathbf{R}_a \cdot \bar{\Pi}_a^{(4)}$ are neglected in the semilinear approximation.

We introduce the usual reduced mass and reduced temperature

$$ \mu_{ab} = \frac{m_a m_b}{m_a + m_b}; \qquad T_{ab} = \frac{m_a T_b + m_b T_a}{m_a + m_b}, \tag{14} $$

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

$$ \begin{aligned} \mathbf{R}_a = & \sum_{b \neq a} \nu_{ab} \left\{ \rho_a (\mathbf{u}_b - \mathbf{u}_a) + \frac{\mu_{ab}}{T_{ab}} \left[ V_{ab(1)} \mathbf{q}_a - V_{ab(2)} \frac{\rho_a}{\rho_b} \mathbf{q}_b \right] \right. \\ & \left. - \frac{3}{56} \left( \frac{\mu_{ab}}{T_{ab}} \right)^2 \left[ \mathbf{X}_a^{(5)} - \frac{\rho_a}{\rho_b} \mathbf{X}_b^{(5)} \right] \right\}, \end{aligned} \tag{15} $$

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

$$ V_{ab(1)} = \frac{(21/10)T_a m_b + (3/5)T_b m_a}{T_a m_b + T_b m_a}; \qquad V_{ab(2)} = \frac{(3/5)T_a m_b + (21/10)T_b m_a}{T_a m_b + T_b m_a}. \tag{16} $$
These and the other mass-ratio coefficients given below come from the Landau collisional operator introduced in Appendices F and G, where one uses perturbed distribution functions of the 21-moment model; see Section 8.4 and Appendix B, with the calculations of the collisional integrals in Appendices K and L. The energy exchange rates are given by

$$Q_a = \sum_{b \neq a} 3 \rho_a \nu_{ab} \frac{T_b - T_a}{m_a + m_b}, \qquad (17)$$

where $|\mathbf{u}_b - \mathbf{u}_a|^2$ are neglected, as discussed above. The heat flux exchange rates are given by

$$\begin{aligned} \mathbf{Q}_a^{(3)\prime} = & -[2\nu_{aa} + \sum_{b \neq a} \nu_{ab} \hat{D}_{ab(1)}] \mathbf{q}_a + \sum_{b \neq a} \nu_{ab} \hat{D}_{ab(2)} \frac{\rho_a}{\rho_b} \mathbf{q}_b \\ & + \left[ \frac{3}{70}\nu_{aa} + \sum_{b \neq a} \nu_{ab} \hat{E}_{ab(1)} \right] \frac{\rho_a}{\rho_a} \mathbf{X}_a^{(5)} - \sum_{b \neq a} \nu_{ab} \hat{E}_{ab(2)} \frac{\rho_b}{\rho_b} \mathbf{X}_b^{(5)} \\ & - p_a \sum_{b \neq a} \nu_{ab}(\mathbf{u}_b - \mathbf{u}_a) \hat{U}_{ab(1)}, \end{aligned} \qquad (18)$$

with mass-ratio coefficients

$$\hat{U}_{ab(1)} = \frac{3m_b(3T_a m_a + T_a m_b - 2T_b m_a)}{2(T_a m_b + T_b m_a)(m_a + m_b)};$$

$$\hat{D}_{ab(1)} = \{75T_a^3 m_a m_b^3 + 95T_a^3 m_b^4 + 174T_a^2 T_b m_a m_b^3 + 300T_a T_b^2 m_a^3 m_b + 498T_a T_b^2 m_a^2 m_b^2 + 60T_b^3 m_a^4 \\ + 104T_b^3 m_a^3 m_b\} [20(T_a m_b + T_b m_a)^3(m_a + m_b)]^{-1};$$

$$\hat{D}_{ab(2)} = \frac{9T_a m_b^2 (10T_a^2 m_a m_b + 6T_a^2 m_b^2 + 45T_a T_b m_a^2 + 27T_a T_b m_a m_b - 14T_b^2 m_a^2)}{20(T_a m_b + T_b m_a)^3 (m_a + m_b)};$$

$$\hat{E}_{ab(1)} = \frac{3T_a m_b (19T_a^2 m_a m_b^2 + 23T_a^2 m_b^3 - 2T_a T_b m_a^2 m_b + 36T_a T_b m_a m_b^2 + 84T_b^2 m_a^3 + 118T_b^2 m_a^2 m_b)}{560(T_a m_b + T_b m_a)^3 (m_a + m_b)};$$

$$\hat{E}_{ab(2)} = \frac{9T_a T_b m_a m_b^2 (7T_a m_a + 5T_a m_b - 2T_b m_a)}{112(T_a m_b + T_b m_a)^3 (m_a + m_b)}. \qquad (19)$$

The fifth-order-moment exchange rates are given by

$$\begin{aligned} \mathbf{Q}_a^{(5)\prime} = & -\left[\frac{76}{5}\nu_{aa} + \sum_{b \neq a} \nu_{ab} \hat{F}_{ab(1)}\right] \frac{p_a}{\rho_a} \mathbf{q}_a + \sum_{b \neq a} \nu_{ab} \hat{F}_{ab(2)} \frac{p_a}{\rho_a} \frac{\rho_a}{\rho_b} \mathbf{q}_b \\ & -\left[\frac{3}{35}\nu_{aa} + \sum_{b \neq a} \nu_{ab} \hat{G}_{ab(1)}\right] \mathbf{X}_a^{(5)} - \sum_{b \neq a} \nu_{ab} \hat{G}_{ab(2)} \frac{p_a}{\rho_b} \mathbf{X}_b^{(5)} \\ & - \frac{p_a^2}{\rho_a} \sum_{b \neq a} \nu_{ab}(\mathbf{u}_b - \mathbf{u}_a) \hat{U}_{ab(2)}, \end{aligned} \qquad (20)$$

with mass-ratio coefficients

$$\hat{U}_{ab(2)} = \frac{3m_b(17T_a^2 m_a m_b + 9T_a^2 m_b^2 + 42T_a T_b m_a^2 + 6T_a T_b m_a m_b - 28T_b^2 m_a^2)}{(T_a m_b + T_b m_a)^2 (m_a + m_b)};$$

$$\hat{F}_{ab(1)} = \{855T_a^5 m_a m_b^4 + 759T_a^5 m_b^5 + 2340T_a^4 T_b m_a^2 m_b^3 + 1972T_a^4 T_b m_a m_b^4 + 2640T_a^3 T_b^2 m_a^3 m_b^2 \\ +2332T_a^3 T_b^2 m_a^2 m_b^3 + 5880T_a^2 T_b^3 m_a^4 m_b + 3324T_a^2 T_b^3 m_a^3 m_b^2 - 3080T_a T_b^4 m_a^4 m_b - 560T_b^5 m_a^5\} \\ \times [10(T_a m_b + T_b m_a)^4 (m_a + m_b)T_a]^{-1};$$

$$\hat{F}_{ab(2)} = 3T_a m_b^2 \{70T_a^3 m_a m_b^2 + 102T_a^3 m_b^3 + 385T_a^2 T_b m_a^2 m_b + 561T_a^2 T_b m_a m_b^2 + 1890T_a T_b^2 m_a^3 \\ +1446T_a T_b^2 m_a^2 m_b - 588T_b^3 m_a^3\} [10(T_a m_b + T_b m_a)^4 (m_a + m_b)]^{-1};$$

$$\hat{G}_{ab(1)} = \{-565T_a^4 m_a m_b^4 + 533T_a^4 m_b^5 + 1270T_a^3 T_b m_a^2 m_b^3 + 1190T_a^3 T_b m_a m_b^4 + 1020T_a^2 T_b^2 m_a^3 m_b^2 \\ +1152T_a^2 T_b^2 m_a^2 m_b^3 + 3640T_a T_b^3 m_a^4 m_b + 1916T_a T_b^3 m_a^3 m_b^2 - 1400T_b^4 m_a^5 - 3304T_b^4 m_a^4 m_b\} \\ \times [280(T_a m_b + T_b m_a)^4 (m_a + m_b)]^{-1};$$

$$\hat{G}_{ab(2)} = - \frac{3T_a T_b m_a m_b^2 (3T_a^2 m_a m_b - 5T_a^2 m_b^2 - 42T_a T_b m_a^2 - 38T_a T_b m_a m_b + 12T_b^2 m_a^2)}{8(T_a m_b + T_b m_a)^4 (m_a + m_b)}. \qquad (21)$$
The exchange rates for the usual stress tensor are given by

$$ \bar{Q}_a^{(2)} = -\frac{21}{10}\nu_{aa}\bar{\Pi}_a^{(2)} + \frac{9}{70}\nu_{aa}\frac{\rho_a}{p_a}\bar{\Pi}_a^{(4)} + \sum_{b \neq a} \frac{\rho_a \nu_{ab}}{m_a + m_b} \left[ -\hat{K}_{ab(1)} \frac{1}{n_a}\bar{\Pi}_a^{(2)} + \hat{K}_{ab(2)} \frac{1}{n_b}\bar{\Pi}_b^{(2)} + L_{ab(1)} \frac{\rho_a}{n_a p_a}\bar{\Pi}_a^{(4)} - L_{ab(2)} \frac{\rho_b}{n_b p_b}\bar{\Pi}_b^{(4)} \right], \quad (22) $$

with mass-ratio coefficients

$$ \begin{aligned} \hat{K}_{ab(1)} &= \frac{10T_a^2 m_a m_b^2 + 15T_a^2 m_b^3 + 35T_a T_b m_a^2 m_b + 42T_a T_b m_a m_b^2 + 10T_b^2 m_a^3 + 12T_b^2 m_a^2 m_b}{5(T_a m_b + T_b m_a)^2 m_a}; \\ \hat{K}_{ab(2)} &= \frac{6T_a^2 m_a m_b + 4T_a^2 m_b^2 + 21T_a T_b m_a^2 + 14T_a T_b m_a m_b - 5T_b^2 m_a^2}{5(T_a m_b + T_b m_a)^2}; \\ L_{ab(1)} &= \frac{3T_a m_b (2T_a m_a m_b + 3T_a m_b^2 + 7T_b m_a^2 + 8T_b m_a m_b)}{35(T_a m_b + T_b m_a)^2 m_a}; \\ L_{ab(2)} &= \frac{3m_a T_b (5T_a m_a + 4T_a m_b - T_b m_a)}{35(T_a m_b + T_b m_a)^2}. \end{aligned} \quad (23) $$

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

$$ \bar{Q}_a^{(4)} = -\frac{53}{20}\nu_{aa}\frac{p_a}{\rho_a}\bar{\Pi}_a^{(2)} - \frac{79}{140}\nu_{aa}\bar{\Pi}_a^{(4)} + \sum_{b \neq a} \nu_{ab} \left[ -\hat{M}_{ab(1)} \frac{p_a}{\rho_a}\bar{\Pi}_a^{(2)} + \hat{M}_{ab(2)} \frac{p_a^2}{\rho_a p_b}\bar{\Pi}_b^{(2)} - N_{ab(1)}\bar{\Pi}_a^{(4)} - N_{ab(2)} \frac{p_a^2}{p_b^2}\frac{\rho_b}{\rho_a}\bar{\Pi}_b^{(4)} \right], \quad (24) $$

with mass-ratio coefficients

$$ \begin{aligned} \hat{M}_{ab(1)} &= \{48T_a^4 m_a m_b^3 + 36T_a^4 m_b^4 + 216T_a^3 T_b m_a^2 m_b^2 + 107T_a^3 T_b m_a m_b^3 + 378T_a^2 T_b^2 m_a^3 m_b \\ &\quad +36T_a^2 T_b^2 m_a^2 m_b^2 - 315T_a T_b^3 m_a^3 m_b - 70T_b^4 m_a^4\} [5(T_a m_b + T_b m_a)^3 T_a(m_b + m_a)]^{-1}; \\ \hat{M}_{ab(2)} &= -\{T_b m_a (18T_a^3 m_a m_b^2 - 4T_a^3 m_b^3 + 81T_a^2 T_b m_a^2 m_b - 18T_a^2 T_b m_a m_b^2 - 147T_a T_b^2 m_a^3 \\ &\quad -189T_a T_b^2 m_a^2 m_b + 35T_b^3 m_a^3\} [5(T_a m_b + T_b m_a)^3 T_a(m_b + m_a)]^{-1}; \\ N_{ab(1)} &= -\{16T_a^3 m_a m_b^3 + 12T_a^3 m_b^4 + 72T_a^2 T_b m_a^2 m_b^2 + 21T_a^2 T_b m_a m_b^3 + 126T_a T_b^2 m_a^3 m_b \\ &\quad -54T_a T_b^2 m_a^2 m_b^2 - 140T_a^3 m_a^4 - 273T_b^3 m_a^3 m_b\} [35(T_a m_b + T_b m_a)^3 (m_b + m_a)]^{-1}; \\ N_{ab(2)} &= -\frac{3T_b^2 m_a^2 (35T_a^2 m_a m_b + 12T_a^2 m_b^2 - 35T_a T_b m_a^2 - 51T_a T_b m_a m_b + 7T_b^2 m_a^2)}{35(T_a m_b + T_b m_a)^3 T_a(m_b + m_a)}. \end{aligned} \quad (25) $$

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

## 2.2. Collisional Contributions for Small Temperature Differences

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

$$ V_{ab(1)} = \frac{(21/10)m_b + (3/5)m_a}{m_b + m_a}; \qquad V_{ab(2)} = \frac{(3/5)m_b + (21/10)m_a}{m_b + m_a}, \quad (26) $$
the mass-ratio coefficients (19) simplify into

$$
\begin{align}
\hat{D}_{ab(1)} &= \frac{3m_a^3 + (86/5)m_a^2 m_b + (77/10)m_a m_b^2 + (19/4)m_b^3}{(m_a + m_b)^3}; \\
\hat{D}_{ab(2)} &= \frac{(279/20)m_a m_b^2 + (27/10)m_b^3}{(m_a + m_b)^3}; \\
\hat{E}_{ab(1)} &= \frac{(9/20)m_a^2 m_b + (6/35)m_a m_b^2 + (69/560)m_b^3}{(m_a + m_b)^3}; \\
\hat{E}_{ab(2)} &= \frac{(45/112)m_a m_b^2}{(m_a + m_b)^3}; \qquad \hat{U}_{ab(1)} = \frac{3}{2} \frac{m_b}{(m_a + m_b)}, \tag{27}
\end{align}
$$

the mass-ratio coefficients (21) become

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

the mass-ratio coefficients (23) become

$$
\begin{align}
\hat{K}_{ab(1)} &= \frac{10m_a^2 + 37m_a m_b + 15m_b^2}{5m_a(m_b + m_a)}; \qquad \hat{K}_{ab(2)} = \frac{4(4m_a + m_b)}{5(m_b + m_a)}; \\
L_{ab(1)} &= \frac{3(7m_a + 3m_b)m_b}{35m_a(m_b + m_a)}; \qquad L_{ab(2)} = \frac{12m_a}{35(m_a + m_b)}, \tag{29}
\end{align}
$$

and the mass-ratio coefficients (25) simplify into

$$
\begin{align}
\hat{M}_{ab(1)} &= - \frac{70m_a^3 - 133m_a^2 m_b - 119m_a m_b^2 - 36m_b^3}{5(m_b + m_a)^3}; & \hat{M}_{ab(2)} &= \frac{4m_a(28m_a^2 - m_a m_b + m_b^2)}{5(m_b + m_a)^3}; \\
N_{ab(1)} &= \frac{140m_a^3 + 7m_a^2 m_b - 25m_a m_b^2 - 12m_b^3}{35(m_b + m_a)^3}; & N_{ab(2)} &= \frac{12m_a^2(7m_a - 3m_b)}{35(m_b + m_a)^3}. \tag{30}
\end{align}
$$

2.3. Reduction to the 13-moment Model

As a partial double-check of our calculations, neglecting the evolution Equations (11)–(12) for $\bar{\Pi}_a^{(4)}$ and $X_a^{(5)}$, and, in the evolution Equations (9)–(10) for $\bar{\Pi}_a^{(2)}$ and $q_a$, prescribing the closures (which are derived from Hermite closures)

$$
\mathbf{X}_a^{(5)} = 28 \frac{p_a}{\rho_a} \mathbf{q}_a; \qquad \bar{\mathbf{\Pi}}_a^{(4)} = 7 \frac{p_a}{\rho_a} \bar{\mathbf{\Pi}}_a^{(2)}, \tag{31}
$$

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

$$
\begin{align}
\mathbf{R}_a &= \sum_{b \neq a} \nu_{ab} \left[ \rho_a (\mathbf{u}_b - \mathbf{u}_a) + \frac{3}{5} \frac{\mu_{ab}}{T_{ab}} \left( \mathbf{q}_a - \frac{\rho_a}{\rho_b} \mathbf{q}_b \right) \right]; \\
\mathbf{Q}_a^{(3)'} &= -\frac{4}{5} \nu_{aa} \mathbf{q}_a + \sum_{b \neq a} \nu_{ab} \left[ -\hat{D}_{ab(1)}^* \mathbf{q}_a + \hat{D}_{ab(2)}^* \frac{\rho_a}{\rho_b} \mathbf{q}_b - p_a (\mathbf{u}_b - \mathbf{u}_a) \hat{U}_{ab(1)} \right]; \\
\bar{\mathbf{Q}}_a^{(2)'} &= -\frac{6}{5} \nu_{aa} \bar{\mathbf{\Pi}}_a^{(2)} + \sum_{b \neq a} \frac{\rho_a \nu_{ab}}{m_a + m_b} \left[ -\hat{K}_{ab(1)}^* \frac{1}{n_a} \bar{\mathbf{\Pi}}_a^{(2)} + \hat{K}_{ab(2)}^* \frac{1}{n_b} \bar{\mathbf{\Pi}}_b^{(2)} \right], \tag{32}
\end{align}
$$
with the mass-ratio coefficients

$$
\begin{align}
\hat{D}_{ab(1)}^* &= \frac{9T_a^2 m_a m_b^2 + 13T_a^2 m_b^3 - 6T_a T_b m_a^2 m_b + 20T_a T_b m_a m_b^2 + 30T_b^2 m_a^3 + 52T_b^2 m_a^2 m_b}{10(m_a + m_b)(T_a m_b + T_b m_a)^2}; \\
\hat{D}_{ab(2)}^* &= \frac{9T_a m_b^2 (5T_a m_a + 3T_a m_b - 2T_b m_a)}{10(m_a + m_b)(T_a m_b + T_b m_a)^2}; \\
\hat{K}_{ab(1)}^* &= \frac{2(2T_a m_a m_b + 3T_a m_b^2 + 5T_b m_a^2 + 6T_b m_a m_b)}{5m_a(T_a m_b + T_b m_a)}, \quad \hat{K}_{ab(2)}^* = \frac{2(3T_a m_a + 2T_a m_b - T_b m_a)}{5(T_a m_b + T_b m_a)}, \tag{33}
\end{align}
$$

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

$$
\begin{align}
\hat{D}_{ab(1)}^* &= \frac{30m_a^2 + 16m_a m_b + 13m_b^2}{10(m_a + m_b)^2}; & \hat{D}_{ab(2)}^* &= \frac{27m_b^2}{10(m_a + m_b)^2}; \nonumber \\
\hat{K}_{ab(1)}^* &= \frac{2m_a + (6/5)m_b}{m_a}; & \hat{K}_{ab(2)}^* &= \frac{4}{5}; & \hat{U}_{ab(1)} &= \frac{3}{2} \frac{m_b}{m_a}. \tag{34}
\end{align}
$$

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

2.4. Semilinear Approximation (Decoupled Stress Tensors and Heat Fluxes)

Here we consider the 21-moment model with evolution Equations (9)–(12) in the semilinear approximation, where additionally viscosity tensors and heat fluxes are decoupled. It will be shown later that the contributions introduced by the coupling are smaller by a factor of 1/ν<sub>aa</sub>. Within the semilinear approximation, we also assume that there are no large-scale gradients of considered fluid moments. For example, the decoupling removes the last terms on the left-hand sides of Equations (10), (11), and (12) proportional to (∇p<sub>a</sub>)q<sub>a</sub>, (∇p<sub>a</sub>) · Π̄<sub>a</sub><sup>(2)</sup> and (∇p<sub>a</sub>) · Π̄<sub>a</sub><sup>(4)</sup>. We also neglect these terms within the semilinear approximation when the coupling is considered (see Sections 5 and 6). In the presence of large-scale gradients in pressure/temperature, these terms might become significant, together with many other terms that are neglected in the semilinear approximation. The evolution equations for the heat flux vectors simplify into

$$
\frac{d_a}{dt} \boldsymbol{q}_a + \Omega_a \hat{\boldsymbol{b}} \times \boldsymbol{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) = \boldsymbol{Q}_a^{(3)\prime}, \qquad (35)
$$

$$
\frac{d_a}{dt} \mathbf{X}_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) = \mathbf{Q}_a^{(5)\prime}, \quad (36)
$$

and the evolution equations for the viscosity tensors become

$$
\frac{d_a}{dt} \bar{\mathbf{\Pi}}_a^{(2)} + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(2)})^S + p_a \bar{\mathbf{W}}_a = \bar{\mathbf{Q}}_a^{(2)\prime}; \qquad (37)
$$

$$
\frac{d_a}{dt} \bar{\mathbf{\Pi}}_a^{(4)} + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{\mathbf{W}}_a = \bar{\mathbf{Q}}_a^{(4)\prime}. \qquad (38)
$$

The above system will be used to recover the transport coefficients of Braginskii (1965). In some instances, it might actually be
advantageous to suppress the nonlinearities in numerical simulations, and to perform multifluid simulations with the system (35)–
(38) instead of the system (9)–(12).
### 3. One Ion–Electron Plasma

#### 3.1. Ion Heat Flux $\mathbf{q}_a$ of Braginskii (Self-collisions)

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

$$ \frac{d_a}{dt} \mathbf{q}_a + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) = -\frac{4}{5} \nu_{aa} \mathbf{q}_a + \frac{3}{70} \nu_{aa} \left( \frac{\rho_a}{p_a} \mathbf{X}_a^{(5)} - 28 \mathbf{q}_a \right); \quad (39) $$

$$ \frac{d_a}{dt} \mathbf{X}_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) = -\frac{88}{5} \nu_{aa} \frac{p_a}{\rho_a} \mathbf{q}_a - \frac{3}{35} \nu_{aa} \left( \mathbf{X}_a^{(5)} - 28 \frac{p_a}{\rho_a} \mathbf{q}_a \right). \quad (40) $$

Neglecting the evolution Equation (40), and prescribing closure (31), which neglects the second term on the right-hand side of (39), yields the ion heat flux model of Burgers–Schunk, with the well-known $-4/5$ constant. However, now the equations read:

$$ \begin{aligned} \frac{d_a}{dt} \mathbf{q}_a + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) &= -2\nu_{aa}\mathbf{q}_a + \frac{3}{70}\nu_{aa}\frac{\rho_a}{p_a}\mathbf{X}_a^{(5)}; \\ \frac{d_a}{dt} \mathbf{X}_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) &= -\frac{76}{5}\nu_{aa}\frac{p_a}{\rho_a}\mathbf{q}_a - \frac{3}{35}\nu_{aa}\mathbf{X}_a^{(5)}. \end{aligned} \quad (41) $$

Prescribing the quasistatic approximation (by canceling the $d_a/dt$) yields an analytic solution (see, for example, the general vector Equation (E23) with solution (E24))

$$ \mathbf{q}_a = -\kappa_{\parallel}^a \nabla_{\parallel} T_a - \kappa_{\perp}^a \nabla_{\perp} T_a + \kappa_{\times}^a \hat{\mathbf{b}} \times \nabla T_a, \quad (42) $$

and thermal conductivities

$$ \begin{aligned} \kappa_{\parallel}^a &= \frac{125}{32} \frac{p_a}{\nu_{aa} m_a}; \\ \kappa_{\perp}^a &= \frac{p_a}{\nu_{aa} m_a} \frac{2x^2 + (648/245)}{x^4 + (3313/1225)x^2 + (20736/30625)}; \\ \kappa_{\times}^a &= \frac{p_a}{\nu_{aa} m_a} \frac{(5/2)x^3 + (2277/490)x}{x^4 + (3313/1225)x^2 + (20736/30625)}, \end{aligned} \quad (43) $$

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

$$ \begin{aligned} \kappa_{\parallel}^a &= 3.906 \frac{p_a}{\nu_{aa} m_a}; \\ \kappa_{\perp}^a &= \frac{p_a}{\nu_{aa} m_a} \frac{2x^2 + 2.645}{x^4 + 2.704x^2 + 0.6771}; \\ \kappa_{\times}^a &= \frac{p_a}{\nu_{aa} m_a} \frac{(5/2)x^3 + 4.647x}{x^4 + 2.704x^2 + 0.6771}, \end{aligned} \quad (44) $$

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

#### 3.2. Ion Heat Flux $\mathbf{X}_a^{(5)}$ (Self-collisions)

The solution for the vector $\mathbf{X}_a^{(5)}$ has a similar form:

$$ \mathbf{X}_a^{(5)} = \frac{p_a}{\rho_a} [-\kappa_{\parallel}^{a(5)} \nabla_{\parallel} T_a - \kappa_{\perp}^{a(5)} \nabla_{\perp} T_a + \kappa_{\times}^{a(5)} \hat{\mathbf{b}} \times \nabla T_a], \quad (45) $$
with “thermal conductivities”

$$ 
\begin{aligned}
\kappa_{\parallel}^{a(5)} &= \frac{2975}{24} \frac{p_a}{\nu_{aa}m_a}; \\
\kappa_{\perp}^{a(5)} &= \frac{P_a}{\nu_{aa}m_a} \frac{44x^2 + (14688/175)}{x^4 + (3313/1225)x^2 + (20736/30625)}; \\
\kappa_{\times}^{a(5)} &= \frac{P_a}{\nu_{aa}m_a} \frac{70x^3 + (1086/7)x}{x^4 + (3313/1225)x^2 + (20736/30625)}.
\end{aligned}
\quad (46) $$

### 3.3. Electron Heat Flux $\mathbf{q}_e$ of Braginskii

Considering a one ion-electron plasma with similar temperatures, and keeping only the dominant term in an $m_e/m_i$ expansion, the mass-ratio coefficients (26), (27), and (28) simplify into

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

and the collisional exchange rates become

$$ \mathbf{R}_e = - \rho_e \nu_{ei} \delta \mathbf{u} + \frac{21}{10} \frac{\rho_e}{p_e} \nu_{ei} \mathbf{q}_e - \frac{3}{56} \frac{\rho_e^2}{p_e^2} \nu_{ei} \mathbf{X}_e^{(5)}; \quad (48) $$

$$ \mathbf{Q}_e^{(3) \prime} = + \frac{3}{2} p_e \nu_{ei} \delta \mathbf{u} - \left[ 2\nu_{ee} + \frac{19}{4} \nu_{ei} \right] \mathbf{q}_e + \left[ \frac{3}{70} \nu_{ee} + \frac{69}{560} \nu_{ei} \right] \frac{\rho_e}{\rho_e} \mathbf{X}_e^{(5)}; \quad (49) $$

$$ \mathbf{Q}_e^{(5) \prime} = +27 \frac{p_e^2}{\rho_e} \nu_{ei} \delta \mathbf{u} - \left[ \frac{76}{5} \nu_{ee} + \frac{759}{10} \nu_{ei} \right] \frac{p_e}{\rho_e} \mathbf{q}_e - \left[ \frac{3}{35} \nu_{ee} - \frac{533}{280} \nu_{ei} \right] \mathbf{X}_e^{(5)}, \quad (50) $$

where $\delta \mathbf{u} = \mathbf{u}_e - \mathbf{u}_i$, and enter the right-hand sides of the electron momentum equation and the evolution equations for the electron heat flux vectors

$$ 
\begin{aligned}
\frac{d_e}{dt} \mathbf{q}_e + \Omega_e \hat{\mathbf{b}} \times \mathbf{q}_e + \frac{5}{2} p_e \nabla \left( \frac{p_e}{\rho_e} \right) &= \mathbf{Q}_e^{(3) \prime}; \\
\frac{d_e}{dt} \mathbf{X}_e^{(5)} + \Omega_e \hat{\mathbf{b}} \times \mathbf{X}_e^{(5)} + 70 \frac{p_e^2}{\rho_e} \nabla \left( \frac{p_e}{\rho_e} \right) &= \mathbf{Q}_e^{(5) \prime}.
\end{aligned}
\quad (51) $$

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

$$ 
\begin{aligned}
\mathbf{Q}_e^{(3) \prime} &= + \frac{3}{2} p_e \nu_{ei} \delta \mathbf{u} - \left[ \frac{\sqrt{2}}{Z_i} + \frac{19}{4} \right] \nu_{ei} \mathbf{q}_e + \left[ \frac{3}{70 \sqrt{2} Z_i} + \frac{69}{560} \right] \nu_{ei} \frac{\rho_e}{p_e} \mathbf{X}_e^{(5)}; \\
\mathbf{Q}_e^{(5) \prime} &= + 27 \frac{p_e^2}{\rho_e} \nu_{ei} \delta \mathbf{u} - \left[ \frac{76}{5 \sqrt{2} Z_i} + \frac{759}{10} \right] \nu_{ei} \frac{p_e}{\rho_e} \mathbf{q}_e - \left[ \frac{3}{35 \sqrt{2} Z_i} - \frac{533}{280} \right] \nu_{ei} \mathbf{X}_e^{(5)}.
\end{aligned}
\quad (52) $$

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

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

$$ 
\begin{aligned}
\mathbf{q}_e^T &= -\kappa_{\parallel}^e \nabla_{\parallel} T_e - \kappa_{\perp}^e \nabla_{\perp} T_e + \kappa_{\times}^e \hat{\mathbf{b}} \times \nabla T_e; \\
\mathbf{q}_e^u &= \beta_0 p_e \delta \mathbf{u}_{\parallel} + p_e \delta \mathbf{u}_{\perp} \frac{\beta_1' x^2 + \beta_0'}{\Delta} - p_e \hat{\mathbf{b}} \times \delta \mathbf{u} \frac{\beta_1'' x^3 + \beta_0'' x}{\Delta},
\end{aligned}
\quad (53) $$
and the heat conductivities are given by

$$ \kappa_{\parallel}^{e} = \frac{p_e}{m_e \nu_{ei}} \gamma_0; \quad \kappa_{\perp}^{e} = \frac{p_e}{m_e \nu_{ei}} \frac{\gamma'_1 x^2 + \gamma'_0}{\Delta}; \quad \kappa_{\times}^{e} = \frac{p_e}{m_e \nu_{ei}} \frac{\gamma''_1 x^3 + \gamma''_0 x}{\Delta}. \qquad (54) $$

The momentum exchange rates are also split into a thermal part and a frictional part $R_e = R_e^T + R_e^u$ (thermal force and friction force), according to

$$ \begin{aligned} \mathbf{R}_e^u &= -\alpha_0 \rho_e \nu_{ei} \delta \mathbf{u}_\parallel - \rho_e \nu_{ei} \delta \mathbf{u}_\perp \left(1 - \frac{\alpha'_1 x^2 + \alpha'_0}{\Delta}\right) - \rho_e \nu_{ei} \hat{\mathbf{b}} \times \delta \mathbf{u} \frac{\alpha''_1 x^3 + \alpha''_0 x}{\Delta}; \\ \mathbf{R}_e^T &= -\beta_0 n_e \nabla_{\parallel} T_e - n_e \nabla_{\perp} T_e \frac{\beta'_1 x^2 + \beta'_0}{\Delta} + n_e \hat{\mathbf{b}} \times \nabla T_e \frac{\beta''_1 x^3 + \beta''_0 x}{\Delta}. \end{aligned} \qquad (55) $$

Instead of the numerical Table II on p. 25 of Braginskii (1965), we provide all of the coefficients in a fully analytic form for a general ion charge $Z_i$, which are given by

$$ \begin{aligned} \alpha_0 &= \frac{4(16Z_i^2 + 61Z_i\sqrt{2} + 72)}{217Z_i^2 + 604Z_i\sqrt{2} + 288}; & \beta_0 &= \frac{30Z_i(11Z_i + 15\sqrt{2})}{217Z_i^2 + 604Z_i\sqrt{2} + 288}; \\ \gamma_0 &= \frac{25Z_i(433Z_i + 180\sqrt{2})}{4(217Z_i^2 + 604Z_i\sqrt{2} + 288)}; \end{aligned} \qquad (56) $$

$$ \begin{aligned} \Delta &= x^4 + \delta_1 x^2 + \delta_0; \\ \delta_0 &= \left( \frac{217Z_i^2 + 604Z_i\sqrt{2} + 288}{700Z_i^2} \right)^2; \\ \delta_1 &= \frac{586601Z_i^2 + 330152Z_i\sqrt{2} + 106016}{78400Z_i^2}; \end{aligned} \qquad (57) $$

$$ \begin{aligned} \alpha'_1 &= \frac{9(40337Z_i + 10996\sqrt{2})}{78400Z_i}; \\ \alpha'_0 &= \frac{9(217Z_i^2 + 604Z_i\sqrt{2} + 288)(17Z_i + 40\sqrt{2})}{490000Z_i^3}; \\ \alpha''_1 &= \frac{477}{280}; \quad \alpha''_0 = \frac{9(64Z_i^2 + 151Z_i\sqrt{2} + 253)}{6125Z_i^2}; \end{aligned} \qquad (58) $$

$$ \begin{aligned} \beta'_1 &= \frac{3(709Z_i + 172\sqrt{2})}{560Z_i}; \\ \beta'_0 &= \frac{3(217Z_i^2 + 604Z_i\sqrt{2} + 288)(11Z_i + 15\sqrt{2})}{49000Z_i^3}; \\ \beta''_1 &= \frac{3}{2}; \quad \beta''_0 = \frac{3(5729Z_i^2 + 6711Z_i\sqrt{2} + 4728)}{19600Z_i^2}; \end{aligned} \qquad (59) $$

$$ \begin{aligned} \gamma'_1 &= \frac{13Z_i + 4\sqrt{2}}{4Z_i}; \\ \gamma'_0 &= \frac{(217Z_i^2 + 604Z_i\sqrt{2} + 288)(433Z_i + 180\sqrt{2})}{78400Z_i^3}; \\ \gamma''_1 &= \frac{5}{2}; \quad \gamma''_0 = \frac{320797Z_i^2 + 202248Z_i\sqrt{2} + 72864}{31360Z_i^2}. \end{aligned} \qquad (60) $$

The numerical values for $Z_i = 1$ are given in the first column of Table II of Braginskii (1965), and the parallel coefficients are $\alpha_0 = 0.5129$, $\beta_0 = 0.7110$, and $\gamma_0 = 3.1616$, for example, matching his values exactly. We checked the entire Table II of Braginskii, and his table is very precise, except for two values. For the $\alpha_0$ coefficient, the values for $Z_i = 2, 3$ should be changed according to $0.4408 \rightarrow 0.4309, 0.3965 \rightarrow 0.3954$. The rest of his table is calculated very accurately, with a handful of irrelevant last digit rounding changes (such as $3.7703 \rightarrow 3.7702$ in $\delta_0$ ($Z_i = 1$) and $0.2400 \rightarrow 0.2399$ in $\alpha''_0$ ($Z_i = 3$); and for the $Z_i = 4$ charge, $0.3752 \rightarrow 0.3751$ in $\alpha_0, 9.055 \rightarrow 9.056$ in $\delta_0$, and $0.4478 \rightarrow 0.4477$ in $\beta'_0$, etc.).

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

### 3.4. Electron Heat Flux $X_e^{(5)}$

Similar to the usual electron heat flux $q_e$, a quasistatic solution for the heat flux vector $X_e^{(5)}$ has to be split into a thermal part and a frictional part, according to

$$ 
\begin{aligned}
X_e^{(5)T} &= \frac{P_e}{\rho_e} [-\kappa_{\parallel}^{e(5)} \nabla_{\parallel} T_e - \kappa_{\perp}^{e(5)} \nabla_{\perp} T_e + \kappa_{\times}^{e(5)} \hat{b} \times \nabla T_e]; \\
X_e^{(5)u} &= \frac{P_e^2}{\rho_e} \left[ \beta_0^{(5)} \delta u_{\parallel} + \frac{\beta_1^{(5)'} x^2 + \beta_0^{(5)'} \Delta}{\Delta} \delta u_{\perp} - \frac{\beta_1^{(5)''} x^3 + \beta_0^{(5)''} \Delta}{\Delta} \hat{b} \times \delta u \right],
\end{aligned} 
\quad (61) $$

with the thermal conductivities

$$ \kappa_{\parallel}^{e(5)} = \frac{P_e}{m_e \nu_{ei}} \gamma_0^{(5)}; \quad \kappa_{\perp}^{e(5)} = \frac{P_e}{m_e \nu_{ei}} \frac{\gamma_1^{(5)'} x^2 + \gamma_0^{(5)'} \Delta}{\Delta}; \quad \kappa_{\times}^{e(5)} = \frac{P_e}{m_e \nu_{ei}} \frac{\gamma_1^{(5)''} x^3 + \gamma_0^{(5)''} \Delta}{\Delta}. \quad (62) $$

The analytic coefficients are given by

$$ 
\begin{aligned}
\beta_0^{(5)} &= \frac{840Z_i(13\sqrt{2} + 12Z_i)}{217Z_i^2 + 604Z_i\sqrt{2} + 288}; \\
\beta_1^{(5)'} &= \frac{3(5829Z_i + 1172\sqrt{2})}{280Z_i}; \\
\beta_0^{(5)'} &= \frac{3(217Z_i^2 + 604Z_i\sqrt{2} + 288)(12Z_i + 13\sqrt{2})}{1750Z_i^3}; \\
\beta_1^{(5)''} &= 27; \quad \beta_0^{(5)''} = \frac{3(7611Z_i^2 + 8429Z_i\sqrt{2} + 5000)}{700Z_i^2},
\end{aligned} 
\quad (63) $$

and

$$ 
\begin{aligned}
\gamma_0^{(5)} &= \frac{175Z_i(204\sqrt{2} + 571Z_i)}{217Z_i^2 + 604Z_i\sqrt{2} + 288}; \\
\gamma_1^{(5)'} &= \frac{113Z_i + 44\sqrt{2}}{2Z_i}; \\
\gamma_0^{(5)'} &= \frac{(217Z_i^2 + 604Z_i\sqrt{2} + 288)(571Z_i + 204\sqrt{2})}{2800Z_i^3}; \\
\gamma_1^{(5)''} &= 70; \quad \gamma_0^{(5)''} = \frac{430783Z_i^2 + 261672Z_i\sqrt{2} + 86880}{1120Z_i^2},
\end{aligned} 
\quad (64) $$

with $\Delta$ unchanged and given by (57). These results were substituted into the momentum exchange rates $R_e$, Equation (48), to obtain the final expression for the friction force and the thermal force. The useful relations are

$$ 
\begin{aligned}
\alpha_0 &= 1 - \frac{21}{10}\beta_0 + \frac{3}{56}\beta_0^{(5)}; \quad \alpha_1' = \frac{21}{10}\beta_1' - \frac{3}{56}\beta_1^{(5)'}; \\
\alpha_0' &= \frac{21}{10}\beta_0' - \frac{3}{56}\beta_0^{(5)'}; \quad \alpha_1'' = \frac{21}{10}\beta_1'' - \frac{3}{56}\beta_1^{(5)''}; \quad \alpha_0'' = \frac{21}{10}\beta_0'' - \frac{3}{56}\beta_0^{(5)''}; \\
\beta_0 &= \frac{21}{10}\gamma_0 - \frac{3}{56}\gamma_0^{(5)}; \quad \beta_1' = \frac{21}{10}\gamma_1' - \frac{3}{56}\gamma_1^{(5)'}; \\
\beta_0' &= \frac{21}{10}\gamma_0' - \frac{3}{56}\gamma_0^{(5)'}; \quad \beta_1'' = \frac{21}{10}\gamma_1'' - \frac{3}{56}\gamma_1^{(5)''}; \quad \beta_0'' = \frac{21}{10}\gamma_0'' - \frac{3}{56}\gamma_0^{(5)''}.
\end{aligned} 
\quad (65) $$

For $Z_i = 1$, the transport coefficients (63), (64) have numerical values

$$ 
\begin{aligned}
\beta_0^{(5)} &= 18.778; \quad \beta_1^{(5)'} = 80.212; \quad \beta_0^{(5)'} = 70.797; \quad \beta_1^{(5)''} = 27; \quad \beta_0^{(5)''} = 105.135; \\
\gamma_0^{(5)} &= 110.664; \quad \gamma_1^{(5)'} = 87.613; \quad \gamma_0^{(5)'} = 417.221; \quad \gamma_1^{(5)''} = 70; \quad \gamma_0^{(5)''} = 792.610.
\end{aligned} 
\quad (66) $$
### 3.5. Ion Viscosity $\bar{\bar{\Pi}}_a^{(2)}$ of Braginskii (Self-collisions)

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

$$ \frac{d_a \bar{\bar{\Pi}}_a^{(2)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(2)})^S + p_a \bar{\bar{W}}_a = - \frac{6}{5} \nu_{aa} \bar{\bar{\Pi}}_a^{(2)} + \frac{9}{70} \nu_{aa} \left( \frac{\rho_a}{p_a} \bar{\bar{\Pi}}_a^{(4)} - 7 \bar{\bar{\Pi}}_a^{(2)} \right); \quad (67) $$

$$ \frac{d_a \bar{\bar{\Pi}}_a^{(4)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{\bar{W}}_a = - \frac{33}{5} \nu_{aa} \frac{p_a}{\rho_a} \bar{\bar{\Pi}}_a^{(2)} - \frac{79}{140} \nu_{aa} \left( \bar{\bar{\Pi}}_a^{(4)} - 7 \frac{p_a}{\rho_a} \bar{\bar{\Pi}}_a^{(2)} \right). \quad (68) $$

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

$$ \begin{aligned} \frac{d_a \bar{\bar{\Pi}}_a^{(2)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(2)})^S + p_a \bar{\bar{W}}_a &= - \frac{21}{10} \nu_{aa} \bar{\bar{\Pi}}_a^{(2)} + \frac{9}{70} \nu_{aa} \frac{\rho_a}{p_a} \bar{\bar{\Pi}}_a^{(4)}; \\ \frac{d_a \bar{\bar{\Pi}}_a^{(4)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{\bar{W}}_a &= - \frac{53}{20} \nu_{aa} \frac{p_a}{\rho_a} \bar{\bar{\Pi}}_a^{(2)} - \frac{79}{140} \nu_{aa} \bar{\bar{\Pi}}_a^{(4)}. \end{aligned} \quad (69) $$

In a quasistatic approximation, the solution of (69) yields $\bar{\bar{\Pi}}_a^{(2)}$ in the following form (see, for example, Appendix E.4):

$$ \begin{aligned} \bar{\bar{\Pi}}_a^{(2)} &= - \eta_0^a \bar{\bar{W}}_0 - \eta_1^a \bar{\bar{W}}_1 - \eta_2^a \bar{\bar{W}}_2 + \eta_3^a \bar{\bar{W}}_3 + \eta_4^a \bar{\bar{W}}_4; \\ \bar{\bar{W}}_0 &= \frac{3}{2} (\bar{\bar{W}}_a : \hat{\mathbf{b}}\hat{\mathbf{b}}) \left( \hat{\mathbf{b}}\hat{\mathbf{b}} - \frac{\bar{\bar{I}}}{3} \right); \\ \bar{\bar{W}}_1 &= \bar{\bar{I}}_\perp : \bar{\bar{W}}_a : \bar{\bar{I}}_\perp + \frac{1}{2} (\bar{\bar{W}}_a : \hat{\mathbf{b}}\hat{\mathbf{b}}) \bar{\bar{I}}_\perp; \\ \bar{\bar{W}}_2 &= (\bar{\bar{I}}_\perp : \bar{\bar{W}}_a : \hat{\mathbf{b}}\hat{\mathbf{b}})^S; \\ \bar{\bar{W}}_3 &= \frac{1}{2} (\hat{\mathbf{b}} \times \bar{\bar{W}}_a : \bar{\bar{I}}_\perp)^S; \\ \bar{\bar{W}}_4 &= (\hat{\mathbf{b}} \times \bar{\bar{W}}_a : \hat{\mathbf{b}}\hat{\mathbf{b}})^S, \end{aligned} \quad (71) $$

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

$$ \begin{aligned} \Pi_{xx}^{a(2)} &= -\frac{\eta_0^a}{2}(W_{xx}^a + W_{yy}^a) - \frac{\eta_1^a}{2}(W_{xx}^a - W_{yy}^a) - \eta_3^a W_{xy}^a; \\ \Pi_{xy}^{a(2)} &= \frac{\eta_3^a}{2}(W_{xx}^a - W_{yy}^a) - \eta_1^a W_{xy}^a; \\ \Pi_{xz}^{a(2)} &= -\eta_4^a W_{yz}^a - \eta_2^a W_{xz}^a; \\ \Pi_{yy}^{a(2)} &= -\frac{\eta_0^a}{2}(W_{xx}^a + W_{yy}^a) + \frac{\eta_1^a}{2}(W_{xx}^a - W_{yy}^a) + \eta_3^a W_{xy}^a; \\ \Pi_{yz}^{a(2)} &= \eta_4^a W_{xz}^a - \eta_2^a W_{yz}^a; \\ \Pi_{zz}^{a(2)} &= -\eta_0^a W_{zz}^a, \end{aligned} \quad (72) $$

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

$$ \begin{aligned} \eta_0^a &= \frac{1025}{1068} \frac{p_a}{\nu_{aa}}; \\ \eta_2^a &= \frac{p_a}{\nu_{aa}} \frac{(6/5)x^2 + (10947/4900)}{x^4 + (79321/19600)x^2 + (71289/30625)}; \\ \eta_4^a &= \frac{p_a}{\nu_{aa}} \frac{x^3 + (46561/19600)x}{x^4 + (79321/19600)x^2 + (71289/30625)}, \end{aligned} \quad (73) $$

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

$$
\begin{align}
\eta_0^a &= 0.960 \frac{p_a}{\nu_{aa}}; \\
\eta_2^a &= \frac{p_a}{\nu_{aa}} \frac{(6/5)x^2 + 2.234}{x^4 + 4.047x^2 + 2.328}; \\
\eta_4^a &= \frac{p_a}{\nu_{aa}} \frac{x^3 + 2.376x}{x^4 + 4.047x^2 + 2.328}, \tag{74}
\end{align}
$$

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

3.6. Ion Viscosity $\bar{\bar{\Pi}}_a^{(4)}$ (Self-collisions)

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

$$
\bar{\bar{\Pi}}_a^{(4)} = \frac{P_a}{\rho_a} [-\eta_0^{a(4)} \bar{\bar{W}}_0 - \eta_1^{a(4)} \bar{\bar{W}}_1 - \eta_2^{a(4)} \bar{\bar{W}}_2 + \eta_3^{a(4)} \bar{\bar{W}}_3 + \eta_4^{a(4)} \bar{\bar{W}}_4], \quad (75)
$$

with the matrices $\bar{\bar{W}}_0 - \bar{\bar{W}}_4$ (71) unchanged, and the viscosities

$$
\begin{align}
\eta_0^{a(4)} &= \frac{8435}{1068} \frac{P_a}{\nu_{aa}}; \\
\eta_2^{a(4)} &= \frac{P_a}{\nu_{aa}} \frac{(33/5)x^2 + (64347/3500)}{x^4 + (79321/19600)x^2 + (71289/30625)}; \\
\eta_4^{a(4)} &= \frac{P_a}{\nu_{aa}} \frac{7x^3 + (59989/2800)x}{x^4 + (79321/19600)x^2 + (71289/30625)}, \tag{76}
\end{align}
$$

where $\eta_1^{a(4)}(x) = \eta_2^{a(4)}(2x)$, $\eta_3^{a(4)}(x) = \eta_4^{a(4)}(2x)$ holds, or with numerical values,

$$
\begin{align}
\eta_0^{a(4)} &= 7.898 \frac{P_a}{\nu_{aa}}; \\
\eta_2^{a(4)} &= \frac{P_a}{\nu_{aa}} \frac{6.600x^2 + 18.385}{x^4 + 4.047x^2 + 2.328}; \\
\eta_4^{a(4)} &= \frac{P_a}{\nu_{aa}} \frac{7x^3 + 21.425x}{x^4 + 4.047x^2 + 2.328}. \tag{77}
\end{align}
$$

3.7. Electron Viscosity $\bar{\bar{\Pi}}_e^{(2)}$ of Braginskii

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

$$
\begin{align}
\hat{K}_{ei(1)} &= 3 \frac{m_i}{m_e}; & \hat{K}_{ei(2)} &= \frac{4}{5}; & L_{ei(1)} &= \frac{9}{35} \frac{m_i}{m_e}; & L_{ei(2)} &= \frac{12}{35} \frac{m_e}{m_i}; \nonumber \\
\hat{M}_{ei(1)} &= \frac{36}{5}; & \hat{M}_{ei(2)} &= \frac{4}{5} \frac{m_e}{m_i}; & N_{ei(1)} &= -\frac{12}{35}; & N_{ei(2)} &= -\frac{36}{35} \frac{m_e^2}{m_i^2}, \tag{78}
\end{align}
$$

and the collisional exchange rates for the viscosity tensors become

$$
\begin{align}
\bar{\bar{Q}}_e^{(2)\prime} &= -\left(\frac{21}{10}\nu_{ee} + 3\nu_{ei}\right)\bar{\bar{\Pi}}_e^{(2)} + \left(\frac{9}{70}\nu_{ee} + \frac{9}{35}\nu_{ei}\right)\frac{\rho_e}{\rho_e}\bar{\bar{\Pi}}_e^{(4)}; \nonumber \\
\bar{\bar{Q}}_e^{(4)\prime} &= -\left(\frac{53}{20}\nu_{ee} + \frac{36}{5}\nu_{ei}\right)\frac{p_e}{\rho_e}\bar{\bar{\Pi}}_e^{(2)} + \left(-\frac{79}{140}\nu_{ee} + \frac{12}{35}\nu_{ei}\right)\bar{\bar{\Pi}}_e^{(4)}. \tag{79}
\end{align}
$$

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

$$
\begin{align}
\bar{\bar{Q}}_e^{(2)\prime} &= -\left(\frac{21}{10Z_i\sqrt{2}} + 3\right)\nu_{ei}\bar{\bar{\Pi}}_e^{(2)} + \left(\frac{9}{70Z_i\sqrt{2}} + \frac{9}{35}\right)\nu_{ei}\frac{\rho_e}{\rho_e}\bar{\bar{\Pi}}_e^{(4)}; \\
\bar{\bar{Q}}_e^{(4)\prime} &= -\left(\frac{53}{20Z_i\sqrt{2}} + \frac{36}{5}\right)\nu_{ei}\frac{p_e}{\rho_e}\bar{\bar{\Pi}}_e^{(2)} + \left(-\frac{79}{140Z_i\sqrt{2}} + \frac{12}{35}\right)\nu_{ei}\bar{\bar{\Pi}}_e^{(4)}, \tag{80}
\end{align}
$$
and these contributions enter the right-hand sides of the evolution equations

$$ 
\begin{aligned}
& \frac{d_e \bar{\Pi}_e^{(2)}}{dt} + \Omega_e (\hat{\mathbf{b}} \times \bar{\Pi}_e^{(2)})^S + p_e \bar{W}_e = \bar{Q}_e^{(2)}' ; \\
& \frac{d_e \bar{\Pi}_e^{(4)}}{dt} + \Omega_e (\hat{\mathbf{b}} \times \bar{\Pi}_e^{(4)})^S + 7 \frac{p_e^2}{\rho_e} \bar{W}_e = \bar{Q}_e^{(4)}' .
\end{aligned}
\qquad (81)
$$

In a quasistatic approximation, the solution of (80), (81) yields the electron viscosity tensor $\bar{\Pi}_e^{(2)}$ in the form (70), (71), with the electron viscosities

$$ 
\begin{aligned}
\eta_0^e &= \frac{p_e}{\nu_{ei}} \frac{5Z_i(408Z_i + 205\sqrt{2})}{6(192Z_i^2 + 301Z_i\sqrt{2} + 178)} ; \\
\eta_2^e &= \frac{p_e}{\nu_{ei}} \left[ \frac{3\sqrt{2} + 6Z_i}{5Z_i} x^2 + \frac{3(192Z_i^2 + 301Z_i\sqrt{2} + 178)(408Z_i + 205\sqrt{2})}{196000Z_i^3} \right] / \Delta ; \\
\eta_4^e &= \frac{p_e}{\nu_{ei}} x \left[ x^2 + \frac{119520Z_i^2 + 101784\sqrt{2}Z_i + 46561}{39200Z_i^2} \right] / \Delta ; \\
\Delta &= x^4 + \frac{212256Z_i^2 + 176376\sqrt{2}Z_i + 79321}{39200Z_i^2} x^2 + \left( \frac{3(192Z_i^2 + 301Z_i\sqrt{2} + 178)}{700Z_i^2} \right)^2 ,
\end{aligned}
\qquad (82)
$$

where $x = \Omega_e/\nu_{ei}$ and the relations $\eta_1^e(x) = \eta_2^e(2x)$, $\eta_3^e(x) = \eta_4^e(2x)$. For the particular case of $Z_i = 1$, these electron viscosities become

$$ 
\begin{aligned}
\eta_0^e &= \frac{p_e}{\nu_{ei}} \frac{2040 + 1025\sqrt{2}}{2220 + 1806\sqrt{2}} ; \\
\eta_2^e &= \frac{p_e}{\nu_{ei}} \left[ \frac{3\sqrt{2} + 6}{5} x^2 + \frac{297987}{98000} \sqrt{2} + \frac{82311}{19600} \right] / \Delta ; \\
\eta_4^e &= \frac{p_e}{\nu_{ei}} x \left[ x^2 + \frac{12723}{4900} \sqrt{2} + \frac{166081}{39200} \right] / \Delta ; \\
\Delta &= x^4 + \left( \frac{22047}{4900} \sqrt{2} + \frac{291577}{39200} \right) x^2 + \left( \frac{1431459}{245000} + \frac{14319}{3500} \sqrt{2} \right) ,
\end{aligned}
\qquad (83)
$$

or with numerical values,

$$ 
\begin{aligned}
\eta_0^e &= 0.73094 \frac{p_e}{\nu_{ei}} ; \\
\eta_2^e &= \frac{p_e}{\nu_{ei}} (2.049x^2 + 8.500) / \Delta ; \\
\eta_4^e &= \frac{p_e}{\nu_{ei}} x(x^2 + 7.909) / \Delta ; \\
\Delta &= x^4 + 13.801x^2 + 11.628,
\end{aligned}
\qquad (84)
$$

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

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

$$ 
\begin{aligned}
\eta_0^{e(4)} &= \frac{p_e}{\nu_{ei}} \frac{35Z_i(552Z_i + 241\sqrt{2})}{6(192Z_i^2 + 301Z_i\sqrt{2} + 178)}; \\
\eta_2^{e(4)} &= \frac{p_e}{\nu_{ei}} \left[ \frac{33\sqrt{2} + 48Z_i}{10Z_i} x^2 + \frac{3(192Z_i^2 + 301Z_i\sqrt{2} + 178)(552Z_i + 241\sqrt{2})}{28000Z_i^3} \right] / \Delta; \\
\eta_4^{e(4)} &= \frac{p_e}{\nu_{ei}} x \left[ 7x^2 + \frac{173088Z_i^2 + 142032Z_i\sqrt{2} + 59989}{5600Z_i^2} \right] / \Delta,
\end{aligned} 
\quad (85) $$

where the denominator $\Delta$ is equivalent to (82). For the particular case of $Z_i = 1$, these electron viscosities become

$$ 
\begin{aligned}
\eta_0^{e(4)} &= \frac{p_e}{\nu_{ei}} \frac{35(241\sqrt{2} + 552)}{6(301\sqrt{2} + 370)}; \\
\eta_2^{e(4)} &= \frac{p_e}{\nu_{ei}} \left[ \frac{33\sqrt{2} + 48}{10} x^2 + \frac{382983}{14000}\sqrt{2} + \frac{523983}{14000} \right] / \Delta; \\
\eta_4^{e(4)} &= \frac{p_e}{\nu_{ei}} x \left[ 7x^2 + \frac{8877}{350}\sqrt{2} + \frac{233077}{5600} \right] / \Delta,
\end{aligned} 
\quad (86) $$

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

$$ 
\begin{aligned}
\eta_0^{e(4)} &= 6.546 \frac{p_e}{\nu_{ei}}; \\
\eta_2^{e(4)} &= \frac{p_e}{\nu_{ei}} (9.467x^2 + 76.114) / \Delta; \\
\eta_4^{e(4)} &= \frac{p_e}{\nu_{ei}} x (7x^2 + 77.489) / \Delta; \\
\Delta &= x^4 + 13.801x^2 + 11.628.
\end{aligned} 
\quad (87) $$

### 4. Generalized Electron Coefficients for Multispecies Plasmas

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

$$ Z_{\text{eff}} = \frac{\sum_i \nu_{ei}}{\sqrt{2} \nu_{ee}} = \frac{\sum_i n_i Z_i^2}{n_e}; \quad \langle \mathbf{u}_i \rangle_{\text{eff}} = \frac{\sum_i \nu_{ei} \mathbf{u}_i}{\sum_i \nu_{ei}}, \quad (88) $$

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

$$ \mathbf{R}_e = -\rho_e(\sum_i \nu_{ei})(\mathbf{u}_e - \langle \mathbf{u}_i \rangle_{\text{eff}}) + \frac{21}{10} \frac{\rho_e}{p_e}(\sum_i \nu_{ei})\mathbf{q}_e - \frac{3}{56} \frac{\rho_e^2}{p_e^2}(\sum_i \nu_{ei})\mathbf{X}_e^{(5)}; \quad (89) $$

$$ \mathbf{Q}_e^{(3)\prime} = +\frac{3}{2} p_e (\sum_i \nu_{ei}) (\mathbf{u}_e - \langle \mathbf{u}_i \rangle_{\text{eff}}) - \left[ \frac{\sqrt{2}}{Z_{\text{eff}}} + \frac{19}{4} \right] (\sum_i \nu_{ei}) \mathbf{q}_e + \left[ \frac{3}{70\sqrt{2}Z_{\text{eff}}} + \frac{69}{560} \right] (\sum_i \nu_{ei}) \frac{\rho_e}{p_e} \mathbf{X}_e^{(5)}; \quad (90) $$

$$ \mathbf{Q}_e^{(5)\prime} = +27 \frac{p_e^2}{\rho_e} (\sum_i \nu_{ei}) (\mathbf{u}_e - \langle \mathbf{u}_i \rangle_{\text{eff}}) - \left[ \frac{76}{5\sqrt{2}Z_{\text{eff}}} + \frac{759}{10} \right] (\sum_i \nu_{ei}) \frac{p_e}{\rho_e} \mathbf{q}_e - \left[ \frac{3}{35\sqrt{2}Z_{\text{eff}}} - \frac{533}{280} \right] (\sum_i \nu_{ei}) \mathbf{X}_e^{(5)}. \quad (91) $$

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

$$ Z_i \rightarrow Z_{\text{eff}}; \quad \nu_{ei} \rightarrow \sum_i \nu_{ei}; \quad \delta \mathbf{u} = \mathbf{u}_e - \mathbf{u}_i \rightarrow \mathbf{u}_e - \langle \mathbf{u}_i \rangle_{\text{eff}}. \quad (92) $$
If the evolution equations can be obtained with the transformation (92), then of course their solutions can be obtained with the same transformation as well. The same transformation applies for the viscous evolution Equations (80), (81) and their solutions. As an example, the generalized (thermal) electron heat of Braginskii (1965) for multispecies plasmas reads

$$ \boldsymbol{q}_e^T = - \kappa_{\parallel}^e \nabla_{\parallel} T_e - \kappa_{\perp}^e \nabla_{\perp} T_e + \kappa_{\times}^e \hat{\boldsymbol{b}} \times \nabla T_e; \qquad (93) $$

$$ \kappa_{\parallel}^{e} = \frac{P_e}{m_e(\sum_i \nu_{ei})} \gamma_0; \qquad \kappa_{\perp}^{e} = \frac{P_e}{m_e(\sum_i \nu_{ei})} \frac{\gamma_1' x^2 + \gamma_0'}{\Delta}; \qquad \kappa_{\times}^{e} = \frac{P_e}{m_e(\sum_i \nu_{ei})} \frac{\gamma_1'' x^3 + \gamma_0'' x}{\Delta}; \qquad (94) $$

$$ \gamma_0 = \frac{25Z_{\text{eff}}(433Z_{\text{eff}} + 180\sqrt{2})}{4(217Z_{\text{eff}}^2 + 604Z_{\text{eff}}\sqrt{2} + 288)}; \qquad \gamma_1' = \frac{13Z_{\text{eff}} + 4\sqrt{2}}{4Z_{\text{eff}}}; \qquad \gamma_1'' = \frac{5}{2}; $$

$$ \gamma_0' = \frac{(217Z_{\text{eff}}^2 + 604Z_{\text{eff}}\sqrt{2} + 288)(433Z_{\text{eff}} + 180\sqrt{2})}{78400Z_{\text{eff}}^3}; \qquad \gamma_0'' = \frac{320797Z_{\text{eff}}^2 + 202248Z_{\text{eff}}\sqrt{2} + 72864}{31360Z_{\text{eff}}^2}; $$

$$ \Delta = x^4 + \delta_1 x^2 + \delta_0; \qquad \delta_0 = \left( \frac{217Z_{\text{eff}}^2 + 604Z_{\text{eff}}\sqrt{2} + 288}{700Z_{\text{eff}}^2} \right)^2; $$

$$ \delta_1 = \frac{586601Z_{\text{eff}}^2 + 330152Z_{\text{eff}}\sqrt{2} + 106016}{78400Z_{\text{eff}}^2}, \qquad (95) $$

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

From the electron momentum equation, the electric field then becomes

$$ \boldsymbol{E} = -\frac{1}{c}\boldsymbol{u}_e \times \boldsymbol{B} - \frac{1}{en_e}\nabla \cdot \bar{\boldsymbol{\rho}}_e + \frac{m_e}{e}\boldsymbol{G} \\ + (\sum_i \nu_{ei}) \left[ +\frac{m_e}{e}(\langle \boldsymbol{u}_i \rangle_{\text{eff}} - \boldsymbol{u}_e) + \frac{21 m_e}{10 e p_e} \boldsymbol{q}_e - \frac{3}{56} \frac{\rho_e^2}{e n_e p_e^2} \boldsymbol{X}_e^{(5)} \right] - \frac{m_e}{e} \frac{d_e \boldsymbol{u}_e}{dt}, \qquad (96) $$

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

## 5. Generalization with Coupling of Stress Tensors and Heat Fluxes

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

$$ \frac{d_a}{dt} \bar{\boldsymbol{\Pi}}_a^{(2)} + \Omega_a (\hat{\boldsymbol{b}} \times \bar{\boldsymbol{\Pi}}_a^{(2)})^S + p_a \bar{\boldsymbol{W}}_a + \frac{2}{5} \left( (\nabla \boldsymbol{q}_a)^S - \frac{2}{3} \bar{\boldsymbol{I}} \nabla \cdot \boldsymbol{q}_a \right) \\ = \bar{\boldsymbol{Q}}_a^{(2)}' = \bar{\boldsymbol{Q}}_a^{(2)} - \frac{\bar{\boldsymbol{I}}}{3} \text{Tr} \bar{\boldsymbol{Q}}_a^{(2)}; \qquad (97) $$

$$ \frac{d_a}{dt} \boldsymbol{q}_a + \Omega_a \hat{\boldsymbol{b}} \times \boldsymbol{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) + \frac{1}{2} \nabla \cdot \bar{\boldsymbol{\Pi}}_a^{(4)} - \frac{5}{2} \frac{p_a}{\rho_a} \nabla \cdot \bar{\boldsymbol{\Pi}}_a^{(2)} \\ = \boldsymbol{Q}_a^{(3)}' = \frac{1}{2} \text{Tr} \bar{\boldsymbol{Q}}_a^{(3)} - \frac{5}{2} \frac{p_a}{\rho_a} \boldsymbol{R}_a; \qquad (98) $$

$$ \frac{d_a}{dt} \bar{\boldsymbol{\Pi}}_a^{(4)} + \Omega_a (\hat{\boldsymbol{b}} \times \bar{\boldsymbol{\Pi}}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{\boldsymbol{W}}_a + \frac{1}{5} \left[ (\nabla \boldsymbol{X}_a^{(5)})^S - \frac{2}{3} \bar{\boldsymbol{I}} (\nabla \cdot \boldsymbol{X}_a^{(5)}) \right] \\ = \bar{\boldsymbol{Q}}_a^{(4)}' = \text{Tr} \bar{\boldsymbol{Q}}_a^{(4)} - \frac{\bar{\boldsymbol{I}}}{3} \text{Tr} \text{Tr} \bar{\boldsymbol{Q}}_a^{(4)}; \qquad (99) $$

$$ \frac{d_a}{dt} \boldsymbol{X}_a^{(5)} + \Omega_a \hat{\boldsymbol{b}} \times \boldsymbol{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) + 18 \frac{p_a}{\rho_a} \nabla \cdot \bar{\boldsymbol{\Pi}}_a^{(4)} - 98 \frac{p_a^2}{\rho_a^2} \nabla \cdot \bar{\boldsymbol{\Pi}}_a^{(2)} \\ = \boldsymbol{Q}_a^{(5)}' = \text{Tr} \text{Tr} \bar{\boldsymbol{Q}}_a^{(5)} - 35 \frac{p_a^2}{\rho_a^2} \boldsymbol{R}_a. \qquad (100) $$

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

$$ \begin{gathered} \mathbf{q}_a = \mathbf{q}_a^{(1)} + \mathbf{q}_a^{(2)}; \quad \mathbf{X}_a^{(5)} = \mathbf{X}_a^{(5,1)} + \mathbf{X}_a^{(5,2)}; \\ \bar{\mathbf{\Pi}}_a^{(2)} = \bar{\mathbf{\Pi}}_a^{(2,1)} + \bar{\mathbf{\Pi}}_a^{(2,2)}; \quad \bar{\mathbf{\Pi}}_a^{(4)} = \bar{\mathbf{\Pi}}_a^{(4,1)} + \bar{\mathbf{\Pi}}_a^{(4,2)}, \end{gathered} \qquad (101) $$

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

$$ \frac{d_a \bar{\mathbf{\Pi}}_a^{(2,1)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(2)})^S + p_a \bar{\mathbf{W}}_a + \frac{2}{5} \left( (\nabla \mathbf{q}_a^{(1)})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \mathbf{q}_a^{(1)} \right) = \bar{\mathbf{Q}}_a^{(2) \prime}, \qquad (102) $$

$$ \frac{d_a \bar{\mathbf{\Pi}}_a^{(4,1)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{\mathbf{W}}_a + \frac{1}{5} \left[ (\nabla \mathbf{X}_a^{(5,1)})^S - \frac{2}{3} \bar{\mathbf{I}} (\nabla \cdot \mathbf{X}_a^{(5,1)}) \right] = \bar{\mathbf{Q}}_a^{(4) \prime}; \qquad (103) $$

$$ \frac{d_a \mathbf{q}_a^{(1)}}{dt} + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) + \frac{1}{2} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(4,1)} - \frac{5}{2} \frac{p_a}{\rho_a} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(2,1)} = \mathbf{Q}_a^{(3) \prime}; \qquad (104) $$

$$ \frac{d_a \mathbf{X}_a^{(5,1)}}{dt} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) + 18 \frac{p_a}{\rho_a} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(4,1)} - 98 \frac{p_a^2}{\rho_a^2} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(2,1)} = \mathbf{Q}_a^{(5) \prime}, \qquad (105) $$

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

## 6. Coupling for Unmagnetized One Ion–Electron Plasma

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

$$ \bar{\mathbf{Y}}_a^{(3,1)} = (\nabla \mathbf{q}_a^{(1)})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \mathbf{q}_a^{(1)}; \quad \bar{\mathbf{Y}}_a^{(5,1)} = (\nabla \mathbf{X}_a^{(5,1)})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \mathbf{X}_a^{(5,1)}, \qquad (106) $$

which are symmetric and traceless, analogous to the matrix $\bar{\mathbf{W}}_a$.

### 6.1. Ion Species (Self-collisions)

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

$$ \begin{align} \bar{\mathbf{\Pi}}_a^{(2,1)} &= -\frac{1025}{1068} \frac{p_a}{\nu_{aa}} \bar{\mathbf{W}}_a; \notag \\ \bar{\mathbf{\Pi}}_a^{(2,2)} &= -\frac{1}{\nu_{aa}} \left[ \frac{79}{534} \bar{\mathbf{Y}}_a^{(3,1)} + \frac{3}{178} \frac{\rho_a}{p_a} \bar{\mathbf{Y}}_a^{(5,1)} + \frac{395}{1068} \frac{\partial \bar{\mathbf{\Pi}}_a^{(2,1)}}{\partial t} + \frac{15}{178} \frac{\rho_a}{p_a} \frac{\partial \bar{\mathbf{\Pi}}_a^{(4,1)}}{\partial t} \right]; \notag \\ \bar{\mathbf{\Pi}}_a^{(4,1)} &= -\frac{8435}{1068} \frac{p_a^2}{\rho_a \nu_{aa}} \bar{\mathbf{W}}_a; \notag \\ \bar{\mathbf{\Pi}}_a^{(4,2)} &= +\frac{1}{\nu_{aa}} \left[ \frac{371}{534} \frac{p_a}{\rho_a} \bar{\mathbf{Y}}_a^{(3,1)} - \frac{49}{178} \bar{\mathbf{Y}}_a^{(5,1)} + \frac{1855}{1068} \frac{p_a}{\rho_a} \frac{\partial \bar{\mathbf{\Pi}}_a^{(2,1)}}{\partial t} - \frac{245}{178} \frac{\partial \bar{\mathbf{\Pi}}_a^{(4,1)}}{\partial t} \right], \tag{107} \end{align} $$
and the heat fluxes become

$$q_a^{(1)} = -\frac{125}{32} \frac{p_a}{m_a \nu_{aa}} \nabla T_a;$$

$$q_a^{(2)} = +\frac{1}{\nu_{aa}} \left[ +\frac{515}{96} \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(2,1)} - \frac{95}{96} \nabla \cdot \bar{\Pi}_a^{(4,1)} - \frac{5}{48} \frac{\partial q_a^{(1)}}{\partial t} - \frac{5}{96} \frac{\rho_a}{p_a} \frac{\partial X_a^{(5,1)}}{\partial t} \right];$$

$$X_a^{(5,1)} = -\frac{2975}{24} \frac{p_a^2}{\rho_a m_a \nu_{aa}} \nabla T_a;$$

$$X_a^{(5,2)} = +\frac{1}{\nu_{aa}} \left[ +\frac{p_a^2}{\rho_a^2} \frac{13825}{72} \nabla \cdot \bar{\Pi}_a^{(2,1)} - \frac{2485}{72} \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(4,1)} + \frac{665}{36} \frac{p_a}{\rho_a} \frac{\partial q_a^{(1)}}{\partial t} - \frac{175}{72} \frac{\partial X_a^{(5,1)}}{\partial t} \right]. \quad (108)$$

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

$$\bar{\Pi}_a^{(2,2)} = +\frac{45575}{17088} \frac{p_a}{m_a \nu_{aa}^2} \left[ 2\nabla\nabla T_a - \frac{2}{3}\bar{I}\nabla^2 T_a \right] + \frac{1164025}{1140624} \frac{p_a}{\nu_{aa}^2} \frac{\partial \bar{W}_a}{\partial t};$$

$$\bar{\Pi}_a^{(4,2)} = +\frac{536725}{17088} \frac{p_a^2}{\rho_a m_a \nu_{aa}^2} \left[ 2\nabla\nabla T_a - \frac{2}{3}\bar{I}\nabla^2 T_a \right] + \frac{10498075}{1140624} \frac{p_a^2}{\rho_a \nu_{aa}^2} \frac{\partial \bar{W}_a}{\partial t}, \quad (109)$$

and the heat flux corrections become

$$\mathbf{q}_a^{(2)} = +\frac{45575}{17088} \frac{p_a^2}{\rho_a \nu_{aa}^2} \nabla \cdot \bar{W}_a + \frac{31625}{4608} \frac{p_a}{m_a \nu_{aa}^2} \frac{\partial \nabla T_a}{\partial t};$$

$$X_a^{(5,2)} = +\frac{1131725}{12816} \frac{p_a^3}{\rho_a^2 \nu_{aa}^2} \nabla \cdot \bar{W}_a + \frac{791875}{3456} \frac{p_a^2}{\rho_a m_a \nu_{aa}^2} \frac{\partial \nabla T_a}{\partial t}. \quad (110)$$

## 6.2. Electron Species

For the electron species, it is useful to introduce the denominator

$$D_1 = 192Z_i^2 + 301\sqrt{2}Z_i + 178, \quad (111)$$

and the solutions for the stress tensors are:

$$\bar{\Pi}_e^{(2,1)} = -\frac{5Z_i(205\sqrt{2} + 408Z_i)}{6D_1} \frac{p_e}{\nu_{ei}} \bar{W}_e;$$

$$\bar{\Pi}_e^{(2,2)} = -\frac{1}{D_1 \nu_{ei}} \left[ \frac{Z_i}{3}(79\sqrt{2} - 96Z_i)\bar{Y}_e^{(3,1)} + 3Z_i(\sqrt{2} + 4Z_i)\frac{\rho_e}{p_e}\bar{Y}_e^{(5,1)} \right. \\ \left. + \frac{5}{6}Z_i(79\sqrt{2} - 96Z_i)\frac{\partial \bar{\Pi}_e^{(2,1)}}{\partial t} + 15Z_i(\sqrt{2} + 4Z_i)\frac{\rho_e}{p_e}\frac{\partial \bar{\Pi}_e^{(4,1)}}{\partial t} \right];$$

$$\bar{\Pi}_e^{(4,1)} = -\frac{35Z_i(241\sqrt{2} + 552Z_i)}{6D_1} \frac{p_e^2}{\rho_e \nu_{ei}} \bar{W}_e;$$

$$\bar{\Pi}_e^{(4,2)} = +\frac{1}{D_1 \nu_{ei}} \left[ \frac{7}{3} Z_i (53\sqrt{2} + 288Z_i) \frac{p_e}{\rho_e} \bar{Y}_e^{(3,1)} - 7Z_i (7\sqrt{2} + 20Z_i) \bar{Y}_e^{(5,1)} \right. \\ \left. + \frac{35}{6} Z_i (53\sqrt{2} + 288Z_i) \frac{p_e}{\rho_e} \frac{\partial \bar{\Pi}_e^{(2,1)}}{\partial t} - 35Z_i (7\sqrt{2} + 20Z_i) \frac{\partial \bar{\Pi}_e^{(4,1)}}{\partial t} \right], \quad (112)$$
with the matrices $\bar{\mathbf{Y}}_e$ defined by (106). For the heat fluxes, it is useful to define the denominator

$$D_2 = 217Z_i^2 + 604Z_i\sqrt{2} + 288, \qquad (113)$$

together with $\delta\mathbf{u} = \mathbf{u}_e - \mathbf{u}_i$, and the results read

$$\mathbf{q}_e^{(1)} = -\frac{25Z_i(180\sqrt{2} + 433Z_i)}{4D_2} \frac{p_e}{m_e \nu_{ei}} \nabla T_e + \frac{30Z_i(15\sqrt{2} + 11Z_i)}{D_2} p_e \delta \mathbf{u};$$

$$\mathbf{q}_e^{(2)} = +\frac{1}{D_2 \nu_{ei}} \left[ \frac{5}{4} Z_i (1236\sqrt{2} + 4097Z_i) \frac{p_e}{\rho_e} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(2,1)} - \frac{5}{4} Z_i (228\sqrt{2} + 709Z_i) \nabla \cdot \bar{\mathbf{\Pi}}_e^{(4,1)} \right. \\ \left. - \frac{5}{2} Z_i (12\sqrt{2} - 533Z_i) \frac{\partial \mathbf{q}_e^{(1)}}{\partial t} - \frac{15}{4} Z_i (4\sqrt{2} + 23Z_i) \frac{\rho_e}{p_e} \frac{\partial \mathbf{X}_e^{(5,1)}}{\partial t} \right];$$

$$\mathbf{X}_e^{(5,1)} = -\frac{175Z_i(204\sqrt{2} + 571Z_i)}{D_2} \frac{p_e^2}{\rho_e m_e \nu_{ei}} \nabla T_e + \frac{840Z_i(13\sqrt{2} + 12Z_i)}{D_2} \frac{p_e^2}{\rho_e} \delta \mathbf{u};$$

$$\mathbf{X}_e^{(5,2)} = \frac{1}{D_2 \nu_{ei}} \left[ +175Z_i(316\sqrt{2} + 1103Z_i) \frac{p_e^2}{\rho_e^2} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(2,1)} - 35Z_i(284\sqrt{2} + 951Z_i) \frac{p_e}{\rho_e} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(4,1)} \right. \\ \left. + 70Z_i(76\sqrt{2} + 759Z_i) \frac{p_e}{\rho_e} \frac{\partial \mathbf{q}_e^{(1)}}{\partial t} - 175Z_i(4\sqrt{2} + 19Z_i) \frac{\partial \mathbf{X}_e^{(5,1)}}{\partial t} \right]. \qquad (114)$$

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

$$\bar{\mathbf{\Pi}}_e^{(2,1)} = -0.7309 \frac{p_e}{\nu_{ei}} \bar{\mathbf{W}}_e;$$

$$\bar{\mathbf{\Pi}}_e^{(2,2)} = -\frac{1}{\nu_{ei}} \left[ 0.006587 \bar{\mathbf{Y}}_e^{(3,1)} + 0.02041 \frac{\rho_e}{p_e} \bar{\mathbf{Y}}_e^{(5,1)} + 0.01647 \frac{\partial \bar{\mathbf{\Pi}}_e^{(2,1)}}{\partial t} + 0.1021 \frac{\rho_e}{p_e} \frac{\partial \bar{\mathbf{\Pi}}_e^{(4,1)}}{\partial t} \right];$$

$$\bar{\mathbf{\Pi}}_e^{(4,1)} = -6.5455 \frac{p_e^2}{\rho_e \nu_{ei}} \bar{\mathbf{W}}_e;$$

$$\bar{\mathbf{\Pi}}_e^{(4,2)} = +\frac{1}{\nu_{ei}} \left[ 1.0644 \frac{p_e}{\rho_e} \bar{\mathbf{Y}}_e^{(3,1)} - 0.2630 \bar{\mathbf{Y}}_e^{(5,1)} + 2.6609 \frac{p_e}{\rho_e} \frac{\partial \bar{\mathbf{\Pi}}_e^{(2,1)}}{\partial t} - 1.3152 \frac{\partial \bar{\mathbf{\Pi}}_e^{(4,1)}}{\partial t} \right]; \qquad (115)$$

$$\mathbf{q}_e^{(1)} = -3.1616 \frac{p_e}{m_e \nu_{ei}} \nabla T_e + 0.7110 p_e \delta \mathbf{u};$$

$$\mathbf{q}_e^{(2)} = +\frac{1}{\nu_{ei}} \left[ 5.3754 \frac{p_e}{\rho_e} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(2,1)} - 0.9486 \nabla \cdot \bar{\mathbf{\Pi}}_e^{(4,1)} + 0.9492 \frac{\partial \mathbf{q}_e^{(1)}}{\partial t} - 0.07906 \frac{\rho_e}{p_e} \frac{\partial \mathbf{X}_e^{(5,1)}}{\partial t} \right];$$

$$\mathbf{X}_e^{(5,1)} = -110.664 \frac{p_e^2}{\rho_e m_e \nu_{ei}} \nabla T_e + 18.7783 \frac{p_e^2}{\rho_e} \delta \mathbf{u};$$

$$\mathbf{X}_e^{(5,2)} = \frac{1}{\nu_{ei}} \left[ 199.554 \frac{p_e^2}{\rho_e^2} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(2,1)} - 34.831 \frac{p_e}{\rho_e} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(4,1)} + 44.625 \frac{p_e}{\rho_e} \frac{\partial \mathbf{q}_e^{(1)}}{\partial t} - 3.1747 \frac{\partial \mathbf{X}_e^{(5,1)}}{\partial t} \right]. \qquad (116)$$
By further applying the quasilinear approximation, the corrections to the electron viscosities become

$$ \bar{\Pi}_e^{(2,2)} = + \frac{25Z_i^2(119520Z_i^2 + 101784Z_i\sqrt{2} + 46561) p_e}{18D_1^2} \frac{\partial \bar{W}_e}{\nu_{ei}^2 \partial t} \\ - \frac{10Z_i^2(11040Z_i^2 + 15557Z_i\sqrt{2} + 8922) p_e}{D_2 D_1} \frac{\partial \left[ (\nabla \delta \mathbf{u})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \delta \mathbf{u} \right]}{\nu_{ei}} \\ + \frac{25Z_i^2(53400Z_i^2 + 366451Z_i\sqrt{2} + 131256) p_e}{12D_2 D_1} \frac{\partial \left[ 2\nabla\nabla T_e - \frac{2}{3}\bar{\mathbf{I}}\nabla^2 T_e \right]}{\nu_{ei}^2 m_e} \quad (117) $$

$$ \bar{\Pi}_e^{(4,2)} = + \frac{175Z_i^2(173088Z_i^2 + 142032Z_i\sqrt{2} + 59989) p_e^2}{18D_1^2} \frac{\partial \bar{W}_e}{\nu_{ei}^2 \rho_e \partial t} \\ - \frac{70Z_i^2(16992Z_i^2 + 23993Z_i\sqrt{2} + 13698) p_e^2}{D_2 D_1} \frac{\partial \left[ (\nabla \delta \mathbf{u})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \delta \mathbf{u} \right]}{\nu_{ei} \rho_e} \\ + \frac{175Z_i^2(834576Z_i^2 + 603679Z_i\sqrt{2} + 220824) p_e^2}{12D_2 D_1} \frac{\partial \left[ 2\nabla\nabla T_e - \frac{2}{3}\bar{\mathbf{I}}\nabla^2 T_e \right]}{\nu_{ei}^2 m_e \rho_e} \quad (118) $$

together with the corrections for the heat fluxes:

$$ \begin{aligned} q_e^{(2)} = & \frac{25Z_i^2(534000Z_i^2 + 366451Z_i\sqrt{2} + 131256) p_e^2}{12D_1 D_2} \frac{\partial \bar{W}_e}{\rho_e \nu_{ei}^2} \nabla \cdot \bar{\mathbf{W}}_e \\ & - \frac{75Z_i^2(5729Z_i^2 + 6711Z_i\sqrt{2} + 4728) p_e}{D_2^2} \frac{\partial \delta \mathbf{u}}{\nu_{ei} \partial t} \\ & + \frac{125Z_i^2(320797Z_i^2 + 202248Z_i\sqrt{2} + 72864) p_e}{8D_2^2} \frac{\partial \nabla T_e}{\nu_{ei}^2 m_e \partial t}; \end{aligned} $$

$$ X_e^{(5,2)} = \frac{175Z_i^2(712272Z_i^2 + 463249Z_i\sqrt{2} + 155208) p_e^3}{3D_1 D_2} \frac{\partial \bar{W}_e}{\rho_e^2 \nu_{ei}^2} \nabla \cdot \bar{\mathbf{W}}_e \\ - \frac{2100Z_i^2(7611Z_i^2 + 8429Z_i\sqrt{2} + 5000) p_e^2}{D_2^2} \frac{\partial \delta \mathbf{u}}{\nu_{ei} \rho_e \partial t} \\ + \frac{875Z_i^2(430783Z_i^2 + 261672Z_i\sqrt{2} + 86880) p_e^2}{2D_2^2} \frac{\partial \nabla T_e}{\nu_{ei}^2 \rho_e m_e \partial t} \quad (119) $$

For the ion charge $Z_i = 1$, the numerical values read

$$ \bar{\Pi}_e^{(2,2)} = +0.6801 \frac{p_e}{\nu_{ei}^2} \frac{\partial \bar{W}_e}{\partial t} - 0.3880 \frac{p_e}{\nu_{ei}} \left[ (\nabla \delta \mathbf{u})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \delta \mathbf{u} \right] \\ + 2.2799 \frac{p_e}{\nu_{ei}^2 m_e} \left[ 2 \nabla \nabla T_e - \frac{2}{3} \bar{\mathbf{I}} \nabla^2 T_e \right]; \quad (120) $$

$$ \bar{\Pi}_e^{(4,2)} = +6.6638 \frac{p_e^2}{\nu_{ei}^2 \rho_e} \frac{\partial \bar{W}_e}{\partial t} - 4.1827 \frac{p_e^2}{\nu_{ei} \rho_e} \left[ (\nabla \delta \mathbf{u})^S - \frac{2}{3} \bar{\mathbf{I}} \nabla \cdot \delta \mathbf{u} \right] \\ + 25.7440 \frac{p_e^2}{\nu_{ei}^2 m_e \rho_e} \left[ 2 \nabla \nabla T_e - \frac{2}{3} \bar{\mathbf{I}} \nabla^2 T_e \right]; \quad (121) $$

together with

$$ q_e^{(2)} = 2.2799 \frac{p_e^2}{\rho_e \nu_{ei}^2} \nabla \cdot \bar{\mathbf{W}}_e - 0.8098 \frac{p_e}{\nu_{ei}} \frac{\partial \delta \mathbf{u}}{\partial t} + 5.7487 \frac{p_e}{\nu_{ei}^2 m_e} \frac{\partial \nabla T_e}{\partial t}; $$

$$ X_e^{(5,2)} = 82.1278 \frac{p_e^3}{\nu_{ei}^2 \nu_{ei}^2} \nabla \cdot \bar{\mathbf{W}}_e - 27.8859 \frac{p_e^2}{\nu_{ei} \rho_e} \frac{\partial \delta \mathbf{u}}{\partial t} + 210.2318 \frac{p_e^2}{\nu_{ei}^2 \rho_e m_e} \frac{\partial \nabla T_e}{\partial t}. \quad (122) $$

The rate-of-strain tensor $\bar{\mathbf{W}}_e$ obviously enters the electron heat fluxes, even in a quasistatic approximation.
## 6.3. Momentum Exchange Rates

The collisional momentum exchange rates $R_e = -R_i$, given by (48), can also be split into first- and second-order $R_e = R_e^{(1)} + R_e^{(2)}$, according to

$$ 
\begin{aligned}
\mathbf{R}_e^{(1)} &= -\rho_e \nu_{ei} \delta\mathbf{u} + \frac{21}{10} \frac{\rho_e}{p_e} \nu_{ei} \mathbf{q}_e^{(1)} - \frac{3}{56} \frac{\rho_e^2}{p_e^2} \nu_{ei} \mathbf{X}_e^{(5,1)}; \\
\mathbf{R}_e^{(2)} &= +\frac{21}{10} \frac{\rho_e}{p_e} \nu_{ei} \mathbf{q}_e^{(2)} - \frac{3}{56} \frac{\rho_e^2}{p_e^2} \nu_{ei} \mathbf{X}_e^{(5,2)}.
\end{aligned} 
\quad (123) $$

Then by using results given in the previous section,

$$ 
\begin{aligned}
\mathbf{R}_e^{(1)} &= -\nu_{ei} \rho_e \frac{(D_2 - 153Z_i^2 - 360Z_i\sqrt{2})}{D_2} \delta\mathbf{u} - \frac{30Z_i(15\sqrt{2} + 11Z_i)}{D_2} n_e \nabla T_e; \\
\mathbf{R}_e^{(2)} &= +\frac{6Z_i(47\sqrt{2} + 69Z_i)}{D_2} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(2,1)} - \frac{6Z_i(11\sqrt{2} + 13Z_i)}{D_2} \frac{\rho_e}{p_e} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(4,1)} \\
&\quad -\frac{12Z_i(29\sqrt{2} + 4Z_i)}{D_2} \frac{\rho_e}{p_e} \frac{\partial \mathbf{q}_e^{(1)}}{\partial t} + \frac{3Z_i(2\sqrt{2} - Z_i)}{D_2} \frac{\rho_e^2}{p_e^2} \frac{\partial \mathbf{X}_e^{(5,1)}}{\partial t},
\end{aligned} 
\quad (124) $$

or for a particular case of $Z_i = 1$,

$$ 
\begin{aligned}
\mathbf{R}_e^{(2)} &= +0.5980 \nabla \cdot \bar{\mathbf{\Pi}}_e^{(2,1)} - 0.1261 \frac{\rho_e}{p_e} \nabla \cdot \bar{\mathbf{\Pi}}_e^{(4,1)} \\
&\quad -0.3974 \frac{\rho_e}{p_e} \frac{\partial \mathbf{q}_e^{(1)}}{\partial t} + 0.004036 \frac{\rho_e^2}{p_e^2} \frac{\partial \mathbf{X}_e^{(5,1)}}{\partial t}.
\end{aligned} 
\quad (125) $$

Finally, at a semilinear level,

$$ 
\begin{aligned}
\mathbf{R}_e^{(2)} &= \frac{10Z_i^2(11040Z_i^2 + 15557Z_i\sqrt{2} + 8922)}{D_1 D_2} \frac{p_e}{\nu_{ei}} \nabla \cdot \bar{\mathbf{W}}_e \\
&\quad -\frac{720Z_i^2(64Z_i^2 + 151Z_i\sqrt{2} + 253)}{D_2^2} \rho_e \frac{\partial(\delta\mathbf{u})}{\partial t} \\
&\quad +\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},
\end{aligned} 
\quad (126) $$

and for $Z_i = 1$, the full momentum exchange rates become

$$ 
\mathbf{R}_e = -0.5129 \nu_{ei} \rho_e \delta\mathbf{u} - 0.7110 n_e \nabla T_e + 0.3880 \frac{p_e}{\nu_{ei}} \nabla \cdot \bar{\mathbf{W}}_e - 0.2068 \rho_e \frac{\partial(\delta\mathbf{u})}{\partial t} + 0.8098 \frac{n_e}{\nu_{ei}} \frac{\partial\nabla T_e}{\partial t}, 
\quad (127) $$

where $\delta\mathbf{u} = \mathbf{u}_e - \mathbf{u}_i$. Only the first two terms of (127) were considered by Braginskii (1965) and Spitzer & Härm (1953; with this latter having slightly different proportionality constants—see Appendix I). A further generalization, by keeping the full $\bar{\mathbf{\Pi}}_a^{(2)}$ and $\bar{\mathbf{\Pi}}_a^{(4)}$ viscosity tensors in the last terms of (104) and (105), brings another three terms to $R_e$ (not shown). Naturally, in a highly collisional regime ($\nu_{ei} \gg \omega$), all of the 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{\mathbf{W}}_e$ enters the momentum exchange rates (even at the linear level), with the contribution $\nabla \cdot \bar{\mathbf{W}}_e = \nabla^2 \mathbf{u}_e + (1/3)\nabla(\nabla \cdot \mathbf{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 quasistatic approximation). The evolution Equations (97)–(100) are of course well defined in the regime of low collisionality.

## 7. Multifluid 22-moment Model

Here we consider a natural generalization of the 21-moment model, by accounting for a fully contracted perturbation of the fourth-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^3v$. The fully contracted (scalar) moment is decomposed into its Maxwellian core and
a perturbation $\tilde{X}_a^{(4)}$ (denoted with tilde), according to

$$X_a^{(4)} = m_a \int |\mathbf{c}_a|^4 f_a d^3v = 15 \frac{p_a^2}{\rho_a} + \tilde{X}_a^{(4)}, \qquad (128)$$

meaning a definition of $\tilde{X}_a^{(4)} = m_a \int |\mathbf{c}_a|^4 (f_a - f_a^{(0)}) d^3v$, where $f_a^{(0)}$ is Maxwellian. The scalar perturbation $\tilde{X}_a^{(4)}$ enters the decomposition of the fourth-order moment

$$X_{ijkl}^{a(4)} = \frac{1}{15} \left( 15 \frac{p_a^2}{\rho_a} + \tilde{X}_a^{(4)} \right) (\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \\ + \frac{1}{7} \left[ \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} \right] + \sigma_{ijkl}^{a(4)\prime}, \qquad (129)$$

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

$$\begin{aligned} & \frac{d_a \mathbf{q}_a}{dt} + \frac{7}{5} \mathbf{q}_a \nabla \cdot \mathbf{u}_a + \frac{7}{5} \mathbf{q}_a \cdot \nabla \mathbf{u}_a + \frac{2}{5} (\nabla \mathbf{u}_a) \cdot \mathbf{q}_a + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) \\ & + \frac{1}{6} \nabla \tilde{X}_a^{(4)} + \frac{1}{2} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(4)} - \frac{5}{2} \frac{p_a}{\rho_a} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(2)} - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{\Pi}}_a) \cdot \bar{\mathbf{\Pi}}_a^{(2)} \\ & = \mathbf{Q}_a^{(3)\prime} \equiv \frac{1}{2} \mathrm{Tr} \bar{\mathbf{Q}}_a^{(3)} - \frac{5}{2} \frac{p_a}{\rho_a} \mathbf{R}_a - \frac{1}{\rho_a} \mathbf{R}_a \cdot \bar{\mathbf{\Pi}}_a^{(2)}; \end{aligned} \qquad (130)$$

$$\begin{aligned} & \frac{d_a \bar{\mathbf{\Pi}}_a^{(4)}}{dt} + \frac{1}{5} \left[ (\nabla \mathbf{X}_a^{(5)})^S - \frac{2}{3} \bar{\mathbf{\bar{I}}} (\nabla \cdot \mathbf{X}_a^{(5)}) \right] + \frac{9}{7} (\nabla \cdot \mathbf{u}_a) \bar{\mathbf{\Pi}}_a^{(4)} + \frac{9}{7} (\bar{\mathbf{\Pi}}_a^{(4)} \cdot \nabla \mathbf{u}_a)^S \\ & + \frac{2}{7} ((\nabla \mathbf{u}_a) \cdot \bar{\mathbf{\Pi}}_a^{(4)})^S - \frac{22}{21} \bar{\mathbf{\bar{I}}} (\bar{\mathbf{\Pi}}_a^{(4)} : \nabla \mathbf{u}_a) - \frac{14}{5 \rho_a} \left[ ((\nabla \cdot \bar{\mathbf{\Pi}}_a) \mathbf{q}_a)^S - \frac{2}{3} \bar{\mathbf{\bar{I}}} (\nabla \cdot \bar{\mathbf{\Pi}}_a) \cdot \mathbf{q}_a \right] \\ & + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(4)})^S + \frac{7}{15} \left( 15 \frac{p_a^2}{\rho_a} + \tilde{X}_a^{(4)} \right) \bar{\mathbf{W}}_a \\ & = \bar{\mathbf{Q}}_a^{(4)\prime} \equiv \mathrm{Tr} \bar{\mathbf{Q}}_a^{(4)} - \frac{\bar{\mathbf{\bar{I}}}}{3} \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{Q}}_a^{(4)} - \frac{14}{5 \rho_a} \left[ (\mathbf{R}_a \mathbf{q}_a)^S - \frac{2}{3} \bar{\mathbf{\bar{I}}} (\mathbf{R}_a \cdot \mathbf{q}_a) \right]; \end{aligned} \qquad (131)$$

$$\begin{aligned} & \frac{d_a \tilde{X}_a^{(4)}}{dt} + \nabla \cdot \mathbf{X}_a^{(5)} - 20 \frac{p_a}{\rho_a} \nabla \cdot \mathbf{q}_a + \frac{7}{3} \tilde{X}_a^{(4)} (\nabla \cdot \mathbf{u}_a) + 4 \left( \bar{\mathbf{\Pi}}_a^{(4)} - 5 \frac{p_a}{\rho_a} \bar{\mathbf{\Pi}}_a^{(2)} \right) : \nabla \mathbf{u}_a \\ & - \frac{8}{\rho_a} (\nabla \cdot \bar{\mathbf{\Pi}}_a) \cdot \mathbf{q}_a = \tilde{Q}_a^{(4)\prime} \equiv \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{Q}}_a^{(4)} - 20 \frac{p_a}{\rho_a} \mathbf{Q}_a - \frac{8}{\rho_a} \mathbf{R}_a \cdot \mathbf{q}_a; \end{aligned} \qquad (132)$$

$$\begin{aligned} & \frac{d_a \mathbf{X}_a^{(5)}}{dt} + \frac{1}{3} \nabla \tilde{X}_a^{(6)} + \nabla \cdot \bar{\mathbf{\Pi}}_a^{(6)} + \frac{9}{5} \mathbf{X}_a^{(5)} (\nabla \cdot \mathbf{u}_a) + \frac{9}{5} \mathbf{X}_a^{(5)} \cdot \nabla \mathbf{u}_a + \frac{4}{5} (\nabla \mathbf{u}_a) \cdot \mathbf{X}_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} \\ & + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) - 35 \frac{p_a^2}{\rho_a^2} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(2)} - \frac{7}{3 \rho_a} (\nabla \cdot \bar{\mathbf{\Pi}}_a^a) \tilde{X}_a^{(4)} - \frac{4}{\rho_a} (\nabla \cdot \bar{\mathbf{\Pi}}_a^a) \cdot \bar{\mathbf{\Pi}}_a^{(4)} \\ & = \mathbf{Q}_a^{(5)\prime} \equiv \mathbf{Q}_a^{(5)} - 35 \frac{p_a^2}{\rho_a^2} \mathbf{R}_a - \frac{7}{3 \rho_a} \mathbf{R}_a \tilde{X}_a^{(4)} - \frac{4}{\rho_a} \mathbf{R}_a \cdot \bar{\mathbf{\Pi}}_a^{(4)}. \end{aligned} \qquad (133)$$

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

$$\tilde{X}_a^{(6)} = m_a \int |\mathbf{c}_a|^6 (f_a - f_a^{(0)}) d^3v = 21 \frac{p_a}{\rho_a} \tilde{X}_a^{(4)}. \qquad (134)$$
In the semilinear approximations, the 22-moment model reads

$$ \frac{d_a \bar{\Pi}_a^{(2)}}{dt} + \Omega_a (\hat{\mathbf{b}} \times \bar{\Pi}_a^{(2)})^S + p_a \bar{W}_a + \frac{2}{5} \left( (\nabla \mathbf{q}_a)^S - \frac{2}{3} \bar{I} \nabla \cdot \mathbf{q}_a \right) = \bar{Q}_a^{(2)}; \quad (135) $$

$$ \frac{d_a \mathbf{q}_a + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) + \frac{1}{2} \nabla \cdot \bar{\Pi}_a^{(4)} - \frac{5}{2} \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(2)}}{dt} + \frac{1}{6} \nabla \tilde{X}_a^{(4)} = \mathcal{Q}_a^{(3)}; \quad (136) $$

$$ \frac{d_a \bar{\Pi}_a^{(4)} + \Omega_a (\hat{\mathbf{b}} \times \bar{\Pi}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{W}_a + \frac{1}{5} \left[ (\nabla \mathbf{X}_a^{(5)})^S - \frac{2}{3} \bar{I} (\nabla \cdot \mathbf{X}_a^{(5)}) \right]}{dt} = \bar{Q}_a^{(4)}; \quad (137) $$

$$ \frac{d_a \tilde{X}_a^{(4)} + \nabla \cdot \mathbf{X}_a^{(5)} - 20 \frac{p_a}{\rho_a} \nabla \cdot \mathbf{q}_a}{dt} = \tilde{Q}_a^{(4)} = \text{TrTr} \bar{Q}_a^{(4)} - 20 \frac{p_a}{\rho_a} Q_a; \quad (138) $$

$$ \frac{d_a \mathbf{X}_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) + 18 \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(4)} - 98 \frac{p_a^2}{\rho_a^2} \nabla \cdot \bar{\Pi}_a^{(2)}}{dt} + 7 \frac{p_a}{\rho_a} \nabla \tilde{X}_a^{(4)} = \mathcal{Q}_a^{(5)}; \quad (139) $$

As discussed in Section 2.4, in the semilinear approximation, we are neglecting terms such as $(\nabla p_a)\mathbf{q}_a$, which might become significant in the presence of large-scale gradients, together with the other terms that are neglected. In comparison to the 21-moment model given by (97)–(100), the evolution Equations (135) and (137) for stress tensors $\bar{\Pi}_a^{(2)}$ and $\bar{\Pi}_a^{(4)}$ remain unchanged. Importantly, the collisional contributions $R_a$, $\bar{Q}_a^{(2)}$, $\mathcal{Q}_a^{(3)}$, $\bar{Q}_a^{(4)}$, and $\mathcal{Q}_a^{(5)}$, given in Section 2.1, remain unchanged as well. The only differences are: (1) the scalar perturbations $\tilde{X}_a^{(4)}$ now enter the left-hand sides of the evolution Equations (136) and (139) for heat fluxes $\mathbf{q}_a$ and $\mathbf{X}_a^{(5)}$; (2) a new evolution Equation (138) for scalar $\tilde{X}_a^{(4)}$ is present, with collisional contributions $\tilde{Q}_a^{(4)}$ that need to be specified; and (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

$$ Q_a = \sum_{b \neq a} Q_{ab} = \sum_{b \neq a} \frac{\rho_a \nu_{ab}}{(m_a + m_b)} \left\{ 3(T_b - T_a) + \hat{P}_{ab(1)} \frac{\rho_a}{n_a p_a} \tilde{X}_a^{(4)} - \hat{P}_{ab(2)} \frac{\rho_b}{n_b p_b} \tilde{X}_b^{(4)} \right\}, \quad (140) $$

with the mass-ratio coefficients

$$ \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}; \quad \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}. \quad (141) $$

Interestingly, the scalar perturbations $\tilde{X}_a^{(4)}$ thus enter the energy exchange rates. For self-collisions, all of the contributions naturally disappear. As also discussed later, in Section 8.1, for multifluid models the conservation of energy $Q_{ab} + Q_{ba} = (\mathbf{u}_b - \mathbf{u}_a) \cdot \mathbf{R}_{ab}$ is only satisfied approximately, because in the semilinear approximation the differences in the drifts $\mathbf{u}_b - \mathbf{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 nonlinearly, with unrestricted 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} + (\mathbf{u}_b - \mathbf{u}_a) \cdot \mathbf{R}_{ab}$.

The collisional exchange rates entering evolution Equation (138) are given by

$$ \tilde{Q}_a^{(4)} = -\frac{4}{5} \nu_{aa} \tilde{X}_a^{(4)} + \sum_{b \neq a} \nu_{ab} \left\{ -\frac{p_a^2 (T_b - T_a)}{\rho_a T_a} \hat{S}_{ab(0)} - \tilde{X}_a^{(4)} \hat{S}_{ab(1)} + \frac{p_a^2 \rho_b}{P_b^2 \rho_a} \tilde{X}_b^{(4)} \hat{S}_{ab(2)} \right\}, \quad (142) $$
with the mass-ratio coefficients

$$ 
\begin{aligned}
\hat{S}_{ab(0)} &= \frac{36T_a m_a m_b}{(T_a m_b + T_b m_a)(m_b + m_a)}; \\
\hat{S}_{ab(1)} &= -\{m_a(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 \\
&\quad - 40T_b^3 m_a^3 - 84T_b^3 m_a^2 m_b)\}[10(T_a m_b + T_b m_a)^3(m_b + m_a)]^{-1}; \\
\hat{S}_{ab(2)} &= \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)},
\end{aligned} 
\quad (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

$$ 
\begin{aligned}
\hat{P}_{ab(1)} &= \frac{3m_b}{10(m_b + m_a)}; & \hat{P}_{ab(2)} &= \frac{3m_a}{10(m_b + m_a)}; \\
\hat{S}_{ab(0)} &= \frac{36m_a m_b}{(m_b + m_a)^2}; & \hat{S}_{ab(1)} &= \frac{2m_a(10m_a^2 + 8m_a m_b + 13m_b^2)}{5(m_b + m_a)^3}; & \hat{S}_{ab(2)} &= \frac{6m_a^2 m_b}{(m_b + m_a)^3},
\end{aligned} 
\quad (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 quasistatic approximation, the solution of Equation (138) becomes

$$ \tilde{X}_a^{(4)} = - \frac{5}{4\nu_{aa}} \left[ \nabla \cdot \mathbf{X}_a^{(5)} - 20 \frac{p_a}{\rho_a} \nabla \cdot \mathbf{q}_a \right]. \quad (145) $$

The quasistatic solution is thus completely determined by the heat fluxes $\mathbf{q}_a$ and $\mathbf{X}_a^{(5)}$, and for a magnetized plasma it has the following form:

$$ \tilde{X}_a^{(4)} = - \frac{5}{4\nu_{aa}} \left\{ \nabla \cdot \left[ \frac{p_a}{\rho_a} \left( -\kappa_{\parallel}^{a(5)} \nabla_{\parallel} T_a - \kappa_{\perp}^{a(5)} \nabla_{\perp} T_a + \kappa_{\times}^{a(5)} \hat{\mathbf{b}} \times \nabla T_a \right) \right] \right. \\ \left. - 20 \frac{p_a}{\rho_a} \nabla \cdot \left( -\kappa_{\parallel}^{a} \nabla_{\parallel} T_a - \kappa_{\perp}^{a} \nabla_{\perp} T_a + \kappa_{\times}^{a} \hat{\mathbf{b}} \times \nabla T_a \right) \right\}, \quad (146) $$

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

It feels natural to define the thermal conductivities (of the moment $\tilde{X}_a^{(4)}$):

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

and result (146) then transforms into

$$ \tilde{X}_a^{(4)} = - \frac{p_a}{\nu_{aa} \rho_a} \nabla \cdot [-\kappa_{\parallel}^{a(4)} \nabla_{\parallel} T_a - \kappa_{\perp}^{a(4)} \nabla_{\perp} T_a + \kappa_{\times}^{a(4)} \hat{\mathbf{b}} \times \nabla T_a] \\ - \frac{5}{4\nu_{aa}} (-\kappa_{\parallel}^{a(5)} \nabla_{\parallel} T_a - \kappa_{\perp}^{a(5)} \nabla_{\perp} T_a + \kappa_{\times}^{a(5)} \hat{\mathbf{b}} \times \nabla T_a) \cdot \nabla \left( \frac{p_a}{\rho_a} \right), \quad (148) $$

with the thermal conductivities

$$ 
\begin{aligned}
\kappa_{\parallel}^{a(4)} &= \frac{1375}{24} \frac{p_a}{\nu_{aa} m_a}; \\
\kappa_{\perp}^{a(4)} &= \frac{p_a}{\nu_{aa} m_a} \frac{5x^2 + (9504/245)}{x^4 + (3313/1225)x^2 + (20736/30625)}; \\
\kappa_{\times}^{a(4)} &= \frac{p_a}{\nu_{aa} m_a} \frac{25x^3 + (3810/49)x}{x^4 + (3313/1225)x^2 + (20736/30625)}.
\end{aligned} 
\quad (149) $$
The second term of (148) is strictly nonlinear, and may be neglected for simplicity. The solution for $\tilde{X}_a^{(4)}$ can thus 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 = \hat{b}\hat{b} \cdot \nabla$ and $\nabla_\perp = \hat{I}_\perp \cdot \nabla = -\hat{b} \times \hat{b} \times \nabla$.

The resulting (148) can be further simplified in the semilinear approximation, where one may use $\nabla \cdot (\hat{b} \times \nabla T_a) = 0$, and so

$$ \tilde{X}_a^{(4)} = + \frac{P_a}{\nu_{aa} \rho_a} [\kappa_\parallel^{a(4)} \nabla_\parallel^2 T_a + \kappa_\perp^{a(4)} \nabla_\perp^2 T_a], \quad (150) $$

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

$$ \tilde{X}_a^{(4)} = + \underbrace{\frac{1375}{24} \frac{p_a^2}{\nu_{aa}^2 \rho_a m_a}}_{57.292} \nabla^2 T_a. \quad (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

$$ Q_{ie} = \frac{\rho_i \nu_{ie}}{m_i} \left[ 3(T_e - T_i) + \frac{3}{10} m_e \left( \frac{\tilde{X}_i^{(4)}}{p_i} - \frac{\tilde{X}_e^{(4)}}{p_e} \right) \right]; \quad Q_{ei} = -Q_{ie} - (\mathbf{u}_e - \mathbf{u}_i) \cdot \mathbf{R}_{ei}. \quad (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 the collisional contributions (142) have a simple form:

$$ \tilde{Q}_e^{(4)} = - \frac{4}{5} \nu_{ee} \tilde{X}_e^{(4)}, \quad (153) $$

determined solely by the electron–electron collisions. A quasistatic solution of Equation (138) then becomes

$$ \tilde{X}_e^{(4)} = - \frac{5\sqrt{2}Z_i}{4\nu_{ei}} \left[ \nabla \cdot \mathbf{X}_e^{(5)} - 20 \frac{p_e}{\rho_e} \nabla \cdot \mathbf{q}_e \right], \quad (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 the thermal and frictional parts $\tilde{X}_e^{(4)} = \tilde{X}_e^{(4)T} + \tilde{X}_e^{(4)u}$, where

$$ \tilde{X}_e^{(4)T} = - \frac{5\sqrt{2}Z_i}{4\nu_{ei}} \left\{ \nabla \cdot \left[ \frac{p_e}{\rho_e} \left( -\kappa_\parallel^{e(5)} \nabla_\parallel T_e - \kappa_\perp^{e(5)} \nabla_\perp T_e + \kappa_\times^{e(5)} \hat{b} \times \nabla T_e \right) \right] \right. \\ \left. - 20 \frac{p_e}{\rho_e} \nabla \cdot \left( -\kappa_\parallel^{e} \nabla_\parallel T_e - \kappa_\perp^{e} \nabla_\perp T_e + \kappa_\times^{e} \hat{b} \times \nabla T_e \right) \right\}; \quad (155) $$

$$ \tilde{X}_e^{(4)u} = - \frac{5Z_i \sqrt{2}}{4\nu_{ei}} \left\{ \nabla \cdot \left[ \frac{p_e^2}{\rho_e} \left( \beta_0^{(5)} \delta\mathbf{u}_\parallel + \frac{\beta_1^{(5)\prime} x^2 + \beta_0^{(5)\prime}}{\Delta} \delta\mathbf{u}_\perp - \frac{\beta_1^{(5)\prime} x^3 + \beta_0^{(5)\prime}}{\Delta} \hat{b} \times \delta\mathbf{u} \right) \right] \right. \\ \left. - 20 \frac{p_e}{\rho_e} \nabla \cdot \left( \beta_0 p_e \delta\mathbf{u}_\parallel + p_e \delta\mathbf{u}_\perp \frac{\beta_1^{\prime} x^2 + \beta_0^{\prime}}{\Delta} - p_e \hat{b} \times \delta\mathbf{u} \frac{\beta_1^{\prime\prime} x^3 + \beta_0^{\prime\prime}}{\Delta} \right) \right\}, \quad (156) $$

with $\delta\mathbf{u} = \mathbf{u}_e - \mathbf{u}_i$. It is again natural to define the electron thermal conductivities (of the moment $\tilde{X}_e^{(4)}$):

$$ \kappa_\parallel^{e(4)} = \frac{5\sqrt{2}Z_i}{4}(\kappa_\parallel^{e(5)} - 20\kappa_\parallel^{e}); \quad \kappa_\perp^{e(4)} = \frac{5\sqrt{2}Z_i}{4}(\kappa_\perp^{e(5)} - 20\kappa_\perp^{e}); \quad \kappa_\times^{e(4)} = \frac{5\sqrt{2}Z_i}{4}(\kappa_\times^{e(5)} - 20\kappa_\times^{e}), \quad (157) $$
together with the transport coefficients

$$
\begin{align}
\beta_0^{(4)} &= \frac{5\sqrt{2}Z_i}{4}(\beta_0^{(5)} - 20\beta_0); & \beta_1^{(4)'} &= \frac{5\sqrt{2}Z_i}{4}(\beta_1^{(5)'} - 20\beta_1'); & \beta_0^{(4)'} &= \frac{5\sqrt{2}Z_i}{4}(\beta_0^{(5)'} - 20\beta_0'); \nonumber \\
\beta_1^{(4)''} &= \frac{5\sqrt{2}Z_i}{4}(\beta_1^{(5)''} - 20\beta_1''); & \beta_0^{(4)''} &= \frac{5\sqrt{2}Z_i}{4}(\beta_0^{(5)''} - 20\beta_0''); & & \tag{158}
\end{align}
$$

$$
\begin{align}
\gamma_0^{(4)} &= \frac{5\sqrt{2}Z_i}{4}(\gamma_0^{(5)} - 20\gamma_0); & \gamma_1^{(4)'} &= \frac{5\sqrt{2}Z_i}{4}(\gamma_1^{(5)'} - 20\gamma_1'); & \gamma_0^{(4)'} &= \frac{5\sqrt{2}Z_i}{4}(\gamma_0^{(5)'} - 20\gamma_0'); \nonumber \\
\gamma_1^{(4)''} &= \frac{5\sqrt{2}Z_i}{4}(\gamma_1^{(5)''} - 20\gamma_1''); & \gamma_0^{(4)''} &= \frac{5\sqrt{2}Z_i}{4}(\gamma_0^{(5)''} - 20\gamma_0''); & & \tag{159}
\end{align}
$$

The thermal and frictional parts then become

$$
\begin{equation}
\begin{split}
\tilde{\chi}_e^{(4)T} = & -\frac{P_e}{\nu_{ei} \rho_e} \nabla \cdot \left( -\kappa_{\parallel}^{e(4)} \nabla_{\parallel} T_e - \kappa_{\perp}^{e(4)} \nabla_{\perp} T_e + \kappa_{\times}^{e(4)} \hat{\mathbf{b}} \times \nabla T_e \right) \\
& - \frac{5\sqrt{2}Z_i}{4\nu_{ei}} \left( -\kappa_{\parallel}^{e(5)} \nabla_{\parallel} T_e - \kappa_{\perp}^{e(5)} \nabla_{\perp} T_e + \kappa_{\times}^{e(5)} \hat{\mathbf{b}} \times \nabla T_e \right) \cdot \nabla \left( \frac{p_e}{\rho_e} \right);
\end{split}
\tag{160}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
\tilde{\chi}_e^{(4)u} = & -\frac{P_e}{\nu_{ei} \rho_e} \nabla \cdot \left( \beta_0^{(4)} p_e \delta \mathbf{u}_{\parallel} + \frac{\beta_1^{(4)'} x^2 + \beta_0^{(4)'} \Delta}{\Delta} p_e \delta \mathbf{u}_{\perp} - \frac{\beta_1^{(4)''} x^3 + \beta_0^{(4)''} x \Delta}{\Delta} p_e \hat{\mathbf{b}} \times \delta \mathbf{u} \right) \\
& - \frac{5\sqrt{2}Z_i}{4\nu_{ei}} \nabla \cdot \left( \beta_0^{(5)} p_e \delta \mathbf{u}_{\parallel} + \frac{\beta_1^{(5)'} x^2 + \beta_0^{(5)'} \Delta}{\Delta} p_e \delta \mathbf{u}_{\perp} - \frac{\beta_1^{(5)''} x^3 + \beta_0^{(5)''} x \Delta}{\Delta} p_e \hat{\mathbf{b}} \times \delta \mathbf{u} \right) \cdot \nabla \left( \frac{p_e}{\rho_e} \right),
\end{split}
\tag{161}
\end{equation}
$$

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

$$
\kappa_{\parallel}^{e(4)} = \frac{P_e}{m_e \nu_{ei}} \gamma_0^{(4)}; \quad \kappa_{\perp}^{e(4)} = \frac{P_e}{m_e \nu_{ei}} \frac{\gamma_1^{(4)'} x^2 + \gamma_0^{(4)'} \Delta}{\Delta}; \quad \kappa_{\times}^{e(4)} = \frac{P_e}{m_e \nu_{ei}} \frac{\gamma_1^{(4)''} x^3 + \gamma_0^{(4)''} x \Delta}{\Delta}, \tag{162}
$$

and the transport coefficients become

$$
\begin{align}
\beta_0^{(4)} &= \frac{150Z_i^2 \sqrt{2} (16\sqrt{2} + 29Z_i)}{217Z_i^2 + 604Z_i \sqrt{2} + 288}; & \beta_1^{(4)'} &= -\frac{3\sqrt{2}(548\sqrt{2} + 1261Z_i)}{224}; \notag \\
\beta_0^{(4)'} &= \frac{3\sqrt{2}(217Z_i^2 + 604Z_i \sqrt{2} + 288)(16\sqrt{2} + 29Z_i)}{9800Z_i^2}; \notag \\
\beta_1^{(4)''} &= -\frac{15Z_i \sqrt{2}}{4}; & \beta_0^{(4)''} &= \frac{3\sqrt{2}(3079Z_i^2 + 3181Z_i \sqrt{2} + 1420)}{490Z_i}; \tag{163}
\end{align}
$$

$$
\begin{align}
\gamma_0^{(4)} &= \frac{250Z_i^2 \sqrt{2} (66\sqrt{2} + 229Z_i)}{217Z_i^2 + 604Z_i \sqrt{2} + 288}; & \gamma_1^{(4)'} &= \frac{5\sqrt{2}(4\sqrt{2} - 17Z_i)}{8}; \notag \\
\gamma_0^{(4)'} &= \frac{\sqrt{2}(217Z_i^2 + 604Z_i \sqrt{2} + 288)(66\sqrt{2} + 229Z_i)}{1960Z_i^2}; \notag \\
\gamma_1^{(4)''} &= 25Z_i \sqrt{2}; & \gamma_0^{(4)''} &= \frac{\sqrt{2}(176437Z_i^2 + 102558Z_i \sqrt{2} + 30480)}{784Z_i}; \tag{164}
\end{align}
$$

$$
\begin{align}
\Delta &= x^4 + \delta_1 x^2 + \delta_0; \notag \\
\delta_0 &= \left( \frac{217Z_i^2 + 604Z_i \sqrt{2} + 288}{700Z_i^2} \right)^2; & \delta_1 &= \frac{586601Z_i^2 + 330152Z_i \sqrt{2} + 106016}{78400Z_i^2}, \tag{165}
\end{align}
$$

and with numerical values for $Z_i = 1$,

$$
\begin{align}
\beta_0^{(4)} &= 8.0576; & \beta_1^{(4)'} &= -38.5624; & \beta_0^{(4)'} &= 30.3787; & \beta_1^{(4)''} &= -5.3033; & \beta_0^{(4)''} &= 77.9054; \notag \\
\gamma_0^{(4)} &= 83.8471; & \gamma_1^{(4)'} &= -10.0260; & \gamma_0^{(4)'} &= 316.1179; & \gamma_1^{(4)''} &= 35.3553; & \gamma_0^{(4)''} &= 634.8735; \notag \\
\delta_0 &= 3.7702; & \delta_1 &= 14.7898. & & & & \tag{166}
\end{align}
$$
At the semilinear level, the solution becomes

$$ \tilde{X}_e^{(4)T} = + \frac{p_e}{\nu_{ei} \rho_e} [\kappa_{\parallel}^{e(4)} \nabla_{\parallel}^2 T_e + \kappa_{\perp}^{e(4)} \nabla_{\perp}^2 T_e]; \quad (167) $$

$$ \tilde{X}_e^{(4)u} = - \frac{p_e^2}{\nu_{ei} \rho_e} \left[ \beta_0^{(4)} \nabla \cdot \delta \mathbf{u}_{\parallel} + \frac{\beta_1^{(4)'} x^2 + \beta_0^{(4)'} \Delta}{\Delta} \nabla \cdot \delta \mathbf{u}_{\perp} - \frac{\beta_1^{(4)''} x^3 + \beta_0^{(4)''} \Delta}{\Delta} \nabla \cdot (\hat{\mathbf{b}} \times \delta \mathbf{u}) \right], \quad (168) $$

and for zero magnetic field,

$$ \tilde{X}_e^{(4)} = \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 \mathbf{u}. \quad (169) $$

# 8. Discussion and Conclusions

Here we discuss various topics that we find to be of importance.

## 8.1. Energy Conservation

The collisional integrals were calculated in a semilinear approximation, where all quantities such as $q_a \cdot (u_b - u_a)$ or $|u_b - u_a|^2$ were neglected and considered small. This approach is typically used for calculations with Landau or Boltzmann collisional operators, and is used in the models of Burgers (1969) and Schunk (1977), for example. Importantly, an exact energy conservation $Q_{ab} + Q_{ba} = (u_b - u_a) \cdot R_{ab}$ cannot be achieved, because the collisional integrals would have to be calculated nonlinearly. 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}; \quad Q_{ei} = -Q_{ie} - (\mathbf{u}_e - \mathbf{u}_i) \cdot \mathbf{R}_{ei}, \quad (170) $$

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

The second particular case involves 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, multifluid calculations can be done for unrestricted drifts (see Burgers 1969; Schunk 1977; and our Appendix G.3), yielding

$$ \mathbf{R}_{ab} = \rho_a \nu_{ab} (\mathbf{u}_b - \mathbf{u}_a) \Phi_{ab}; \quad (171) $$

$$ Q_{ab} = \rho_a \nu_{ab} \left[ 3 \frac{T_b - T_a}{m_a + m_b} \Psi_{ab} + \frac{m_b}{m_a + m_b} |\mathbf{u}_b - \mathbf{u}_a|^2 \Phi_{ab} \right], \quad (172) $$

where one defines the functions

$$ \Psi_{ab} = e^{-\epsilon^2}; \quad \Phi_{ab} = \left( \frac{3}{4} \sqrt{\pi} \frac{\text{erf}(\epsilon)}{\epsilon^3} - \frac{3}{2} \frac{e^{-\epsilon^2}}{\epsilon^2} \right); \quad \epsilon = \frac{|\mathbf{u}_b - \mathbf{u}_a|}{\sqrt{\nu_{\text{tha}}^2 + \nu_{\text{thb}}^2}}, \quad (173) $$

the thermal velocities $\nu_{\text{tha}}^2 = 2T_a/m_a$, and the collisional frequencies (178). Because $\rho_a \nu_{ab} = \rho_b \nu_{ba}$ holds, both momentum and energy are conserved. The 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 the differences in the drift velocities $|u_b - u_a|$ become much smaller than the thermal velocities, so that $\epsilon \ll 1$, functions $\Phi_{ab} \to 1$ and $\Psi_{ab} \to 1$ and $R_{ab} = \rho_a \nu_{ab} (u_b - u_a)$. To correctly account for the small $|u_b - 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:

$$ Q_{ab} = \rho_a \nu_{ab} \left[ 3 \frac{T_b - T_a}{m_a + m_b} \left( 1 - \frac{|\mathbf{u}_b - \mathbf{u}_a|^2}{\frac{2T_a}{m_a} + \frac{2T_b}{m_b}} \right) + \frac{m_b}{m_b + m_a} |\mathbf{u}_b - \mathbf{u}_a|^2 \right]; \quad (174) $$

$$ Q_{ab} = \rho_a \nu_{ab} \left[ 3 \frac{T_b - T_a}{m_a + m_b} + \frac{3}{2} \left( \frac{T_a m_b}{T_a m_b + T_b m_a} - \frac{1}{3} \frac{m_b}{m_b + m_a} \right) |\mathbf{u}_b - \mathbf{u}_a|^2 \right]; \quad (175) $$
and 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} \left[ 3 \frac{T_b - T_a}{m_a + m_b} + \frac{m_b}{m_a + m_b} |\mathbf{u}_b - \mathbf{u}_a|^2 \right]. \quad (176)$$

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

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

$$Q_{ab} = \frac{\rho_a \nu_{ab}}{(m_a + m_b)} \left[ 3(T_b - T_a) + \hat{P}_{ab(1)} \frac{\rho_a}{n_a p_a} \tilde{X}_a^{(4)} - \hat{P}_{ab(2)} \frac{\rho_b}{n_b p_b} \tilde{X}_b^{(4)} \right]; \\ \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}; \quad \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}, \quad (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_{ab}^{-1} = \frac{16}{3} \sqrt{\pi} \frac{n_b e^4 Z_a^2 Z_b^2 \ln \Lambda}{m_a^2 (\nu_{tha}^2 + \nu_{thb}^2)^{3/2}} \left(1 + \frac{m_a}{m_b}\right), \quad (178)$$

where $\nu_{tha}^2 = 2T_a/m_a$, and $\rho_a \nu_{ab} = \rho_b \nu_{ba}$ holds. Equivalently, in the form of Burgers (1969) and Schunk (1977),

$$\nu_{ab} = \frac{16}{3} \sqrt{\pi} \left( \frac{\mu_{ab}}{2T_{ab}} \right)^{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}, \quad (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}}. \quad (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} m_a} \frac{\sqrt{\mu_{ab}}}{\sqrt{m_e}}, \quad (181)$$

which identifies with Equation (7.6) of Braginskii (1965; after one uses $\nu_{ab} = n_b \mu_{ab} \alpha_{ab}'/m_a$). For a particular case of a one ion–electron plasma, the collisional frequencies simplify into

$$\begin{aligned} \nu_{ii} &= \frac{4}{3} \sqrt{\pi} \frac{n_i e^4 Z_i^4 \ln \Lambda}{T_i^{3/2} \sqrt{m_i}}; & \nu_{ie} &= \frac{4}{3} \sqrt{2\pi} \frac{n_e e^4 Z_i^2 \ln \Lambda}{T_e^{3/2} \sqrt{m_i}} \sqrt{\frac{m_e}{m_i}}; \\ \nu_{ee} &= \frac{4}{3} \sqrt{\pi} \frac{n_e e^4 \ln \Lambda}{T_e^{3/2} \sqrt{m_e}}; & \nu_{ei} &= \frac{4}{3} \sqrt{2\pi} \frac{n_i e^4 Z_i^2 \ln \Lambda}{T_e^{3/2} \sqrt{m_e}}, \end{aligned} \quad (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$), but $\nu_{ee} \sim \nu_{ei}$, with the 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}$, $\nu_{ii}$ does not. Thus, comparing the Braginskii (1965) expressions (2.5i) and (2.5e) with the (182) definitions, Braginskii clearly uses

$$\tau_i = \tau_{ii}; \quad \tau_e = \tau_{ei}, \quad (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:

$$m_a \ll m_b: \quad \nu_{ab} = \tau_{ab}^{-1} = \frac{16}{3} \sqrt{\pi} \frac{n_b e^4 Z_a^2 Z_b^2 \ln \Lambda}{m_a^2 \nu_{tha}^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}}, \quad (184)$$
which, for example, agrees with the appendix of Helander & Sigmar (2002, p. 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 of obtaining (184) from the general (178). Importantly, $\rho_a\nu_{ab} \neq \rho_b\nu_{ba}$, and if one were to 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:

$$ 
\begin{aligned}
\nu_{ii} &= \frac{4}{3} \sqrt{2\pi} \frac{n_i e^4 Z_i^4 \ln \Lambda}{T_i^{3/2} \sqrt{m_i}}; & \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_e^{3/2} \sqrt{m_i}} \sqrt{\frac{m_e}{m_i}}; \\
\nu_{ee} &= \frac{4}{3} \sqrt{2\pi} \frac{n_e e^4 \ln \Lambda}{T_e^{3/2} \sqrt{m_e}}; & \nu_{ei} &= \frac{4}{3} \sqrt{2\pi} \frac{n_i e^4 Z_i^2 \ln \Lambda}{T_e^{3/2} \sqrt{m_e}}.
\end{aligned} 
\quad (185) $$

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

$$ \tau_i = \sqrt{2} \tau_{ii}; \quad \tau_e = \tau_{ei}. \quad (186) $$

Also, for $Z_i = 1$, the relation $\nu_{ee} = \nu_{ei}$ now holds. These definitions of the collisional frequencies are used in the 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 multispecies 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) for an ion-electron plasma also. 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 one needs to add $\nu_{ee} + \nu_{ei}$, for example. Calculating this addition according to (185) would be incorrect, and one has to use (182) instead. A comparison of the 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 nth 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

$$ 
\begin{aligned}
\mathbf{X}_a^{(2n+1)} &= m_a \int \mathbf{c}_a |\mathbf{c}_a|^{2n} f_a \, d^3v; \\
\bar{\mathbf{\Pi}}_a^{(2n)} &= m_a \int \left( \mathbf{c}_a \mathbf{c}_a - \frac{\bar{\mathbf{I}}}{3} |\mathbf{c}_a|^2 \right) |\mathbf{c}_a|^{2n-2} f_a \, d^3v; \\
\mathbf{X}_a^{(2n)} &= m_a \int |\mathbf{c}_a|^{2n} f_a \, d^3v = (2n+1)!! \frac{p_a^n}{\rho_a^{n-1}} + \tilde{\mathbf{X}}_a^{(2n)},
\end{aligned} 
\quad (187) $$

together with the collisional contributions

$$ 
\begin{aligned}
\bar{\mathcal{Q}}_a^{(2n+1)} &= m_a \int |\mathbf{c}_a|^{2n} \mathbf{c}_a C(f_a) d^3v; \\
\bar{\mathcal{Q}}_a^{(2n)} &= m_a \int |\mathbf{c}_a|^{2n-2} \mathbf{c}_a \mathbf{c}_a C(f_a) d^3v; \\
\mathcal{Q}_a^{(2n)} &= m_a \int |\mathbf{c}_a|^{2n} C(f_a) d^3v; & \mathcal{Q}_a &= \frac{m_a}{2} \int |\mathbf{c}_a|^2 C(f_a) d^3v, 
\end{aligned} 
\quad (188) $$

where, to prevent incompatibility with the previous notation, we use $\mathcal{Q}$ (mathcal{Q}) instead of $\mathcal{Q}$ for vectors and matrices. The new notation fixes the problem that, for example, $\mathcal{Q}_a^{(3)}$ was used for the right-hand side of the evolution equation of the heat flux $q_a$, and not for $\mathbf{X}_a^{(3)}$. It also clarifies that in the vector notation, the matrix $\bar{\mathcal{Q}}_a^{(2n)} = \text{TrTr} \dots \text{Tr} \bar{\mathcal{Q}}_a^{(2n)}$. Note that $\mathcal{Q}_a^{(2)} = 2\mathcal{Q}_a$, $\mathbf{X}_a^{(2)} = 3p_a$, $\tilde{\mathbf{X}}_a^{(2)} = 0$, $\mathbf{X}_a^{(3)} = 2q_a$, and $\mathbf{X}_a^{(1)} = 0$.

Fully nonlinear evolution equations are given in Appendix D; see (D13)–(D15). In the semilinear approximation, these simplify into evolution equations for vectors valid for $n \ge 1$:

$$ 
\begin{aligned}
& \frac{d}{dt} \mathbf{X}_a^{(2n+1)} + \frac{1}{3} \nabla \tilde{\mathbf{X}}_a^{(2n+2)} + \nabla \cdot \bar{\mathbf{\Pi}}_a^{(2n+2)} - \frac{(2n+3)!!}{3} \frac{p_a^n}{\rho_a^n} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(2)} + \Omega_a \hat{\mathbf{b}} \times \mathbf{X}_a^{(2n+1)} \\
& \qquad + (2n+3)!! \frac{(n)}{3} \frac{p_a^n}{\rho_a^{n-1}} \nabla \left( \frac{p_a}{\rho_a} \right) = \bar{\mathcal{Q}}_a^{(2n+1)}, \quad \bar{\mathcal{Q}}_a^{(2n+1)} = \frac{(2n+3)!!}{3} \frac{p_a^n}{\rho_a^n} \mathbf{R}_a,
\end{aligned} 
\quad (189) $$
stress tensors valid for $n \ge 1$:

$$ \frac{d_a \bar{\Pi}_a^{(2n)}}{dt} + \frac{1}{5} \left[ (\nabla X_a^{(2n+1)})^S - \frac{2}{3} \bar{I} \nabla \cdot X_a^{(2n+1)} \right] + \Omega_a (\hat{b} \times \bar{\Pi}_a^{(2n)})^S \\ + \frac{(2n+3)!!}{15} \frac{p_a^n}{\rho_a^{n-1}} \bar{W}_a = \bar{Q}_a^{(2n)}' = \bar{Q}_a^{(2n)} - \frac{\bar{I}}{3} Q_a^{(2n)}, \qquad (190) $$

and scalar perturbations valid for $n \ge 2$:

$$ \frac{d_a \tilde{\chi}_a^{(2n)}}{dt} + \nabla \cdot X_a^{(2n+1)} - (2n+1)!! \frac{(2n)}{3} \left( \frac{p_a}{\rho_a} \right)^{n-1} \nabla \cdot q_a \\ = \tilde{Q}_a^{(2n)}' = Q_a^{(2n)} - (2n+1)!! \frac{(2n)}{3} \left( \frac{p_a}{\rho_a} \right)^{n-1} Q_a, \qquad (191) $$

where $(n)$ without a 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 the previous notation, $\tilde{Q}_a^{(3)}' = 2Q_a^{(3)}'$, $\tilde{Q}_a^{(5)}' = Q_a^{(5)}'$, $\tilde{Q}_a^{(2)}' = \tilde{Q}_a^{(2)}'$, and $\tilde{Q}_a^{(4)}' = \tilde{Q}_a^{(4)}'$.

### 8.4. Reducible and Irreducible Hermite Polynomials

The irreducible Hermite polynomials $H(\tilde{c})$ (notation without tilde) are usually defined through Laguerre–Sonine polynomials $L(\tilde{c})$ (see, for example, Equation (G1.4.4) on p. 326 of Balescu 1988):

$$ H^{(2n)}(\tilde{c}) = L_n^{(1/2)}\left(\frac{\tilde{c}^2}{2}\right); \qquad H_i^{(2n+1)}(\tilde{c}) = \sqrt{\frac{3}{2}} \tilde{c}_i L_n^{(3/2)}\left(\frac{\tilde{c}^2}{2}\right); \\ H_{ij}^{(2n)}(\tilde{c}) = \sqrt{\frac{15}{8}} (\tilde{c}_i \tilde{c}_j - \frac{\tilde{c}^2}{3} \delta_{ij}) L_{n-1}^{(5/2)}\left(\frac{\tilde{c}^2}{2}\right), \qquad (192) $$

where we use tilde for the normalized fluctuating velocity $\tilde{c} = \sqrt{m_a/T_a c_a}$, with the species indices dropped. In our calculations, we find it more natural to use the reducible Hermite polynomials $\tilde{H}(\tilde{c})$ (notation with tilde) of Grad, defined according to

$$ \tilde{H}_{r_1 r_2 \dots r_n}^{(m)}(\tilde{c}) = (-1)^m e^{\frac{\tilde{c}^2}{2}} \frac{\partial}{\partial \tilde{c}_{r_1}} \frac{\partial}{\partial \tilde{c}_{r_2}} \dots \frac{\partial}{\partial \tilde{c}_{r_n}} e^{-\frac{\tilde{c}^2}{2}}. \qquad (193) $$

Applying a sufficient number of contractions then yields definitions of fully contracted scalars, vectors, and matrices:

$$ \tilde{H}^{(2n)} = \tilde{H}_{r_1 r_1 \dots r_n r_n}^{(2n)}; \qquad \tilde{H}_i^{(2n+1)} = \tilde{H}_{i r_1 r_1 \dots r_n r_n}^{(2n+1)}; \qquad \tilde{H}_{ij}^{(2n)} = \tilde{H}_{ij r_1 r_1 \dots r_{n-1} r_{n-1}}^{(2n)} \qquad (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)}. \qquad (195) $$

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

$$ H^{(2n)} = \left( \frac{1}{2^n n! (2n+1)!!} \right)^{1/2} \tilde{H}^{(2n)}; \qquad H_i^{(2n+1)} = \left( \frac{3}{2^n n! (2n+3)!!} \right)^{1/2} \tilde{H}_i^{(2n+1)}; \\ H_{ij}^{(2n)} = \left( \frac{15}{2^n (n-1)! (2n+3)!!} \right)^{1/2} \hat{H}_{ij}^{(2n)}, \qquad (196) $$

with both approaches using essentially the same polynomials, the only difference being the location of the normalization factors. The reducible Hermite polynomials are used to define the Hermite moments

$$ \tilde{h}^{(2n)} = \frac{1}{n_a} \int f_a \tilde{H}^{(2n)} d^3 c; \qquad \tilde{h}_i^{(2n+1)} = \frac{1}{n_a} \int f_a \tilde{H}_i^{(2n+1)} d^3 c; \qquad \hat{h}_{ij}^{(2n)} = \frac{1}{n_a} \int f_a \hat{H}_{ij}^{(2n)} d^3 c, \qquad (197) $$

and analogously for the irreducible ones. Note that the scalar $\tilde{h}^{(2)} = 0$, and we thus often use $\hat{h}_{ij}^{(2)} = \tilde{h}_{ij}^{(2)} = (1/n_a) \int f_a \tilde{H}_{ij}^{(2)} d^3 c$. Finally, by using orthogonality relations, one obtains the perturbation $\chi_a$ of the distribution function $f_a = f_a^{(0)}(1 + \chi_a)$ around the
<table><caption>**Table 1**
Summary of the Various Models with the Perturbation $\chi_a$ Given in Reducible Hermite Moments</caption><thead><tr><th>Model Name</th><th>Corresponding Perturbation of f<sub>a</sub> = f<sub>a</sub><sup>(0)</sup> (1 + χ<sub>a</sub>) in Hermite Moments</th></tr></thead><tbody><tr><td>5-moment</td><td>χ<sub>a</sub> = 0</td></tr><tr><td>8-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup></td></tr><tr><td>10-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup></td></tr><tr><td>13-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup> + <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup></td></tr><tr><td>20-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup> + <sup>1</sup>&frasl;<sub>6</sub>h&#x0302;<sub>ijk</sub><sup>(3)</sup>H&#x0302;<sub>ijk</sub><sup>(3)</sup></td></tr><tr><td>21-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup> + <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup> + <sup>1</sup>&frasl;<sub>28</sub>h&#x0302;<sub>ij</sub><sup>(4)</sup>H&#x0302;<sub>ij</sub><sup>(4)</sup> + <sup>1</sup>&frasl;<sub>280</sub>h&#x0302;<sub>i</sub><sup>(5)</sup>H&#x0302;<sub>i</sub><sup>(5)</sup></td></tr><tr><td>22-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup> + <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup> + <sup>1</sup>&frasl;<sub>28</sub>h&#x0302;<sub>ij</sub><sup>(4)</sup>H&#x0302;<sub>ij</sub><sup>(4)</sup> + <sup>1</sup>&frasl;<sub>120</sub>h&#x0302;<sub>i</sub><sup>(4)</sup>H&#x0302;<sub>i</sub><sup>(4)</sup> + <sup>1</sup>&frasl;<sub>280</sub>h&#x0302;<sub>i</sub><sup>(5)</sup>H&#x0302;<sub>i</sub><sup>(5)</sup></td></tr><tr><td>9-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup> + <sup>1</sup>&frasl;<sub>120</sub>h&#x0302;<sub>i</sub><sup>(4)</sup>H&#x0302;<sub>i</sub><sup>(4)</sup></td></tr><tr><td>11-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup> + <sup>1</sup>&frasl;<sub>280</sub>h&#x0302;<sub>i</sub><sup>(5)</sup>H&#x0302;<sub>i</sub><sup>(5)</sup></td></tr><tr><td>12-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>10</sub>h&#x0302;<sub>i</sub><sup>(3)</sup>H&#x0302;<sub>i</sub><sup>(3)</sup> + <sup>1</sup>&frasl;<sub>120</sub>h&#x0302;<sub>i</sub><sup>(4)</sup>H&#x0302;<sub>i</sub><sup>(4)</sup> + <sup>1</sup>&frasl;<sub>280</sub>h&#x0302;<sub>i</sub><sup>(5)</sup>H&#x0302;<sub>i</sub><sup>(5)</sup></td></tr><tr><td>15-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup> + <sup>1</sup>&frasl;<sub>28</sub>h&#x0302;<sub>ij</sub><sup>(4)</sup>H&#x0302;<sub>ij</sub><sup>(4)</sup></td></tr><tr><td>16-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2</sub>h&#x0302;<sub>ij</sub><sup>(2)</sup>H&#x0302;<sub>ij</sub><sup>(2)</sup> + <sup>1</sup>&frasl;<sub>28</sub>h&#x0302;<sub>ij</sub><sup>(4)</sup>H&#x0302;<sub>ij</sub><sup>(4)</sup> + <sup>1</sup>&frasl;<sub>120</sub>h&#x0302;<sub>i</sub><sup>(4)</sup>H&#x0302;<sub>i</sub><sup>(4)</sup></td></tr></tbody></table>

Note. 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.

Maxwellian $f_a^{(0)}$, in the following form:

$$ \chi_a = \sum_{n=1}^{N} \left[ \frac{15}{2^n(n-1)!(2n+3)!!} \hat{h}_{ij}^{(2n)} \hat{H}_{ij}^{(2n)} + \frac{1}{2^n n! (2n+1)!!} \tilde{h}^{(2n)} \tilde{H}^{(2n)} + \frac{3}{2^n n! (2n+3)!!} \tilde{h}_i^{(2n+1)} \tilde{H}_i^{(2n+1)} \right]; \quad (198) $$

$$ \chi_a = \sum_{n=1}^{N} [h_{ij}^{(2n)} H_{ij}^{(2n)} + h^{(2n)} H^{(2n)} + h_i^{(2n+1)} H_i^{(2n+1)}], \quad (199) $$

and the two approaches are equivalent. Alternatively, because $\hat{h}_{ij}^{(2n)}$ are traceless, it is possible to use $\hat{h}_{ij}^{(2n)} \hat{H}_{ij}^{(2n)} = \hat{h}_{ij}^{(2n)} \tilde{H}_{ij}^{(2n)}$. 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}_i^{(4)} \tilde{H}_i^{(4)} + \frac{1}{280} \tilde{h}_i^{(5)} \tilde{H}_i^{(5)}, \quad (200) $$

with Hermite polynomials

$$ \begin{aligned} \tilde{H}_i^{(3)} &= \tilde{c}_i(\tilde{c}^2 - 5); & \tilde{H}_i^{(5)} &= \tilde{c}_i(\tilde{c}^4 - 14\tilde{c}^2 + 35); \\ \hat{H}_{ij}^{(2)} &= (\tilde{c}_i\tilde{c}_j - \frac{\delta_{ij}}{3}\tilde{c}^2); & \hat{H}_{ij}^{(4)} &= (\tilde{c}_i\tilde{c}_j - \frac{\delta_{ij}}{3}\tilde{c}^2)(\tilde{c}^2 - 7); & \tilde{H}^{(4)} &= \tilde{c}^4 - 10\tilde{c}^2 + 15, \end{aligned} \quad (201) $$

and neglecting $\tilde{h}^{(4)} = 0$ (meaning $\tilde{X}^{(4)} = 0$) yields the 21-moment model. The transformation from Hermite to fluid moments is done according to

$$ \begin{aligned} \tilde{\mathbf{h}}_a^{(3)} &= \frac{2}{p_a} \sqrt{\frac{m_a}{T_a}} \mathbf{q}_a; & \tilde{\mathbf{h}}_a^{(5)} &= \frac{1}{p_a} \sqrt{\frac{m_a}{T_a}} \left( \frac{m_a}{T_a} \mathbf{X}_a^{(5)} - 28 \mathbf{q}_a \right); \\ \hat{\tilde{\mathbf{h}}}_a^{(2)} &= \hat{\tilde{\mathbf{h}}}_a^{(2)} = \frac{1}{p_a} \bar{\mathbf{\Pi}}_a^{(2)}; & \hat{\tilde{\mathbf{h}}}_a^{(4)} &= \frac{\rho_a}{p_a^2} \bar{\mathbf{\Pi}}_a^{(4)} - \frac{7}{p_a} \bar{\mathbf{\Pi}}_a^{(2)}; & \tilde{\mathbf{h}}_a^{(4)} &= \frac{\rho_a}{p_a^2} \tilde{X}_a^{(4)}. \end{aligned} \quad (202) $$

Various models are summarized in Tables 1 and 2. In Table 1, the perturbation $\chi_a$ is given in reducible Hermite moments, and in Table 2 the perturbation is given in fluid moments.
<table><caption>**Table 2**
Summary of the Various Models with the Perturbation $\chi_a$ Given in Fluid Moments</caption><thead><tr><th>Model Name</th><th>Corresponding Perturbation of f<sub>a</sub> = f<sub>a</sub><sup>(0)</sup>(1 + χ<sub>a</sub>) in Fluid Moments</th></tr></thead><tbody><tr><td>5-moment</td><td>χ<sub>a</sub> = 0</td></tr><tr><td>8-moment</td><td>χ<sub>a</sub> = - <sup>m<sub>a</sub></sup>&frasl;<sub>p<sub>a</sub>T<sub>a</sub></sub> (q<sub>a</sub> · c<sub>a</sub>) (1 - <sup>m<sub>a</sub></sup>&frasl;<sub>5T<sub>a</sub>c<sub>a</sub><sup>2</sup></sub>)</td></tr><tr><td>10-moment</td><td>χ<sub>a</sub> = <sup>m<sub>a</sub></sup>&frasl;<sub>2p<sub>a</sub>T<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c<sub>a</sub>c<sub>a</sub>)</td></tr><tr><td>13-moment</td><td>χ<sub>a</sub> = <sup>m<sub>a</sub></sup>&frasl;<sub>2p<sub>a</sub>T<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c<sub>a</sub>c<sub>a</sub>) - <sup>m<sub>a</sub></sup>&frasl;<sub>p<sub>a</sub>T<sub>a</sub></sub> (q<sub>a</sub> · c<sub>a</sub>) (1 - <sup>m<sub>a</sub></sup>&frasl;<sub>5T<sub>a</sub>c<sub>a</sub><sup>2</sup></sub>)</td></tr><tr><td>20-moment</td><td>χ<sub>a</sub> = <sup>m<sub>a</sub></sup>&frasl;<sub>2p<sub>a</sub>T<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c<sub>a</sub>c<sub>a</sub>) + <sup>m<sub>a</sub><sup>2</sup></sup>&frasl;<sub>6p<sub>a</sub>T<sub>a</sub><sup>2</sup></sub> (c<sub>a</sub> · q&#x0302;<sub>a</sub> : c<sub>a</sub>c<sub>a</sub>) - <sup>m<sub>a</sub></sup>&frasl;<sub>p<sub>a</sub>T<sub>a</sub></sub> (q<sub>a</sub> · c<sub>a</sub>)</td></tr><tr><td>21-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2p<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c&#x0302;<sub>a</sub>c&#x0302;<sub>a</sub>) + <sup>1</sup>&frasl;<sub>28</sub> [ <sup>p<sub>a</sub><sup>2</sup></sup>&frasl;<sub>p<sub>a</sub><sup>2</sup></sub> (Π&#x0302;<sub>a</sub><sup>(4)</sup> : c&#x0302;<sub>a</sub>c&#x0302;<sub>a</sub>) - <sup>7</sup>&frasl;<sub>p<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c&#x0302;<sub>a</sub>c&#x0302;<sub>a</sub>) ] (c&#x0302;<sub>a</sub><sup>2</sup> - 7) + <sup>1</sup>&frasl;<sub>5p<sub>a</sub></sub> √<sub>T<sub>a</sub></sub> (q<sub>a</sub> · c&#x0302;<sub>a</sub>) (c&#x0302;<sub>a</sub><sup>2</sup> - 5) + <sup>1</sup>&frasl;<sub>280p<sub>a</sub></sub> √<sub>T<sub>a</sub></sub> [ <sup>p<sub>a</sub></sup>&frasl;<sub>p<sub>a</sub></sub> (X<sub>a</sub><sup>(5)</sup> · c&#x0302;<sub>a</sub>) - 28(q<sub>a</sub> · c&#x0302;<sub>a</sub>) ] (c&#x0302;<sub>a</sub><sup>4</sup> - 14c&#x0302;<sub>a</sub><sup>2</sup> + 35)</td></tr><tr><td>22-moment</td><td>χ<sub>a</sub> = <sup>1</sup>&frasl;<sub>2p<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c&#x0302;<sub>a</sub>c&#x0302;<sub>a</sub>) + <sup>1</sup>&frasl;<sub>28</sub> [ <sup>p<sub>a</sub><sup>2</sup></sup>&frasl;<sub>p<sub>a</sub><sup>2</sup></sub> (Π&#x0302;<sub>a</sub><sup>(4)</sup> : c&#x0302;<sub>a</sub>c&#x0302;<sub>a</sub>) - <sup>7</sup>&frasl;<sub>p<sub>a</sub></sub> (Π&#x0302;<sub>a</sub><sup>(2)</sup> : c&#x0302;<sub>a</sub>c&#x0302;<sub>a</sub>) ] (c&#x0302;<sub>a</sub><sup>2</sup> - 7) + <sup>1</sup>&frasl;<sub>5p<sub>a</sub></sub> √<sub>T<sub>a</sub></sub> (q<sub>a</sub> · c&#x0302;<sub>a</sub>) (c&#x0302;<sub>a</sub><sup>2</sup> - 5) + <sup>1</sup>&frasl;<sub>280p<sub>a</sub></sub> √<sub>T<sub>a</sub></sub> [ <sup>p<sub>a</sub></sup>&frasl;<sub>p<sub>a</sub></sub> (X<sub>a</sub><sup>(5)</sup> · c&#x0302;<sub>a</sub>) - 28(q<sub>a</sub> · c&#x0302;<sub>a</sub>) ] (c&#x0302;<sub>a</sub><sup>4</sup> - 14c&#x0302;<sub>a</sub><sup>2</sup> + 35) + <sup>1</sup>&frasl;<sub>120</sub> <sup>p<sub>a</sub></sup>&frasl;<sub>p<sub>a</sub><sup>2</sup></sub> X&#x0302;<sub>a</sub><sup>(4)</sup> (c&#x0302;<sub>a</sub><sup>4</sup> - 10c&#x0302;<sub>a</sub><sup>2</sup> + 15)</td></tr></tbody></table>

Note. The results for the 21- and 22-moment models are written with normalized $\tilde{c}_a = \sqrt{m_a/T_a c_a}$.

## 8.5. Rosenbluth Potentials (22-moment Model)

Here we summarize the Rosenbluth potentials, defined according to

$$H_b(v) = \int \frac{f_b(v')}{|v' - v|} d^3v'; \quad G_b(v) = \int |v' - v| f_b(v') d^3v', \qquad (203)$$

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

$$H_b(v) = n_b \sqrt{\frac{m_b}{T_b}} \left\{ \frac{1}{\tilde{y}} \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) - \sqrt{\frac{2}{\pi}} \frac{e^{-\tilde{y}^2/2}}{10} \left( \tilde{y} \cdot \tilde{h}_b^{(3)} + (\tilde{y}^2 - 5) \frac{\tilde{y} \cdot \tilde{h}_b^{(5)}}{28} \right) \right. \\ \left. + \frac{1}{2} (\hat{\tilde{h}}_b^{(2)} : \tilde{y} \tilde{y}) \left[ \frac{3}{\tilde{y}^5} \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) - \sqrt{\frac{2}{\pi}} \left(\frac{1}{\tilde{y}^2} + \frac{3}{\tilde{y}^4}\right) e^{-\tilde{y}^2/2} \right] \right. \\ \left. - \frac{1}{28} (\hat{\tilde{h}}_b^{(4)} : \tilde{y} \tilde{y}) \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} + \frac{1}{120} \tilde{h}_b^{(4)} (3 - \tilde{y}^2) \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} \right\}; \qquad (204)$$

$$G_b(v) = n_b \sqrt{\frac{T_b}{m_b}} \left\{ \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} + \left(\tilde{y} + \frac{1}{\tilde{y}}\right) \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) \right. \\ \left. + \left(\frac{\operatorname{erf}(\tilde{y}/\sqrt{2})}{5\tilde{y}^3} - \sqrt{\frac{2}{\pi}} \frac{e^{-\tilde{y}^2/2}}{5\tilde{y}^2}\right) \tilde{y} \cdot \tilde{h}_b^{(3)} - \sqrt{\frac{2}{\pi}} \frac{e^{-\tilde{y}^2/2}}{140} \tilde{y} \cdot \tilde{h}_b^{(5)} \right. \\ \left. - \frac{1}{2} (\hat{\tilde{h}}_b^{(2)} : \tilde{y} \tilde{y}) \left[ \frac{3}{\tilde{y}^4} \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} + \left(\frac{1}{\tilde{y}^3} - \frac{3}{\tilde{y}^5}\right) \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) \right] \right. \\ \left. - \frac{1}{14} (\hat{\tilde{h}}_b^{(4)} : \tilde{y} \tilde{y}) \left[ \sqrt{\frac{2}{\pi}} \left(\frac{1}{\tilde{y}^2} + \frac{3}{\tilde{y}^4}\right) e^{-\tilde{y}^2/2} - \frac{3}{\tilde{y}^5} \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) \right] - \frac{1}{60} \tilde{h}_b^{(4)} \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} \right\}, \qquad (205)$$

where we use the variable

$$\tilde{y} = \sqrt{\frac{m_b}{T_b}} (v - u_b). \qquad (206)$$
<table><caption>**Table 3**
Summary of (MHD) Hermite Closures, Together with Corresponding Fluid Closures</caption><thead><tr><th>Hermite Closures</th><th>Fluid Closures</th></tr></thead><tbody><tr><td>h&#x0302;<sub>i</sub><sup>(3)</sup> = 0</td><td>X<sub>i</sub><sup>(3)</sup> = 0</td></tr><tr><td>h&#x0302;<sub>i</sub><sup>(4)</sup> = 0</td><td>X&#x0302;<sup>(4)</sup> = 0</td></tr><tr><td>h&#x0302;<sub>i</sub><sup>(5)</sup> = 0</td><td>X<sub>i</sub><sup>(5)</sup> = 14<sup>p</sup>&frasl;<sub>p</sub>X<sub>i</sub><sup>(3)</sup></td></tr><tr><td>h&#x0302;<sup>(6)</sup> = 0</td><td>X&#x0302;<sup>(6)</sup> = 21<sup>p</sup>&frasl;<sub>p</sub>X&#x0302;<sup>(4)</sup></td></tr><tr><td>h&#x0302;<sub>i</sub><sup>(7)</sup> = 0</td><td>X<sub>i</sub><sup>(7)</sup> = 27<sup>p</sup>&frasl;<sub>p</sub>X<sub>i</sub><sup>(5)</sup> - 189<sup>p<sup>2</sup></sup>&frasl;<sub>p<sup>2</sup></sub>X<sub>i</sub><sup>(3)</sup></td></tr><tr><td>h&#x0302;<sup>(8)</sup> = 0</td><td>X&#x0302;<sup>(8)</sup> = 36<sup>p</sup>&frasl;<sub>p</sub>X&#x0302;<sup>(6)</sup> - 378<sup>p<sup>2</sup></sup>&frasl;<sub>p<sup>2</sup></sub>X&#x0302;<sup>(4)</sup></td></tr><tr><td>h&#x0302;<sub>i</sub><sup>(9)</sup> = 0</td><td>X<sub>i</sub><sup>(9)</sup> = 44<sup>p</sup>&frasl;<sub>p</sub>X<sub>i</sub><sup>(7)</sup> - 594<sup>p<sup>2</sup></sup>&frasl;<sub>p<sup>2</sup></sub>X<sub>i</sub><sup>(5)</sup> + 2772<sup>p<sup>3</sup></sup>&frasl;<sub>p<sup>3</sup></sub>X<sub>i</sub><sup>(3)</sup></td></tr></tbody></table>

**Note.** Species indices “a” are dropped. The usual heat flux $q_i = X_i^{(3)}/2$. Note that beyond the fourth-order moment, both classes start to differ. It can be shown that erroneously prescribing closures at the last retained moment, such as $X_i^{(5)} = 0$ or $X&#x0302;<sup>(6)</sup> = 0$, leads to unphysical instabilities (unless one prescribes $X_i^{(3)}$ or $X&#x0302;<sup>(4)</sup> = 0$ as well), which is later demonstrated in Appendix B.8, Table 5. A general form for the closures corresponding to $\tilde{h}_i^{(2n+1)} = 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 6.

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

$$ \mathbf{A}_{ab}(\mathbf{v}) = 2 \frac{c_{ab}}{m_a^2} \left(1 + \frac{m_a}{m_b}\right) \frac{\partial H_b(\mathbf{v})}{\partial \mathbf{v}}; \quad \tilde{\mathbf{D}}_{ab}(\mathbf{v}) = 2 \frac{c_{ab}}{m_a^2} \frac{\partial^2 G_b(\mathbf{v})}{\partial \mathbf{v} \partial \mathbf{v}}; \quad c_{ab} = 2\pi e^4 Z_a^2 Z_b^2 \ln \Lambda; \\ C_{ab}(f_a, f_b) = -\frac{\partial}{\partial \mathbf{v}} \cdot \left[ \mathbf{A}_{ab} f_a - \frac{1}{2} \frac{\partial}{\partial \mathbf{v}} \cdot (\tilde{\mathbf{D}}_{ab} f_a) \right]. \tag{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 nonlinear and could potentially be useful for constructing more sophisticated models that could capture collisional effects beyond the semilinear approximation, or perhaps for exploring the runaway effect numerically. All of 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. The 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}_a^{(6)} = 0$ and $\tilde{h}_{ij}^{(6)} = 0$, which translate into fluid closures (134) and (13).

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

$$ \begin{aligned} X_a^{(2n+1)} &= \sum_{m=1}^{n-1} (-1)^{m+n+1} \left(\frac{p_a}{\rho_a}\right)^{n-m} \frac{n!}{m!(n-m)!} \frac{(2n+3)!!}{(2m+3)!!} X_a^{(2m+1)}; \\ \tilde{X}_a^{(2n)} &= \sum_{m=2}^{n-1} (-1)^{m+n+1} \left(\frac{p_a}{\rho_a}\right)^{n-m} \frac{n!}{m!(n-m)!} \frac{(2n+1)!!}{(2m+1)!!} \tilde{X}_a^{(2m)}, \end{aligned} \tag{208} $$

together with closures for stress-tensors

$$ \tilde{\mathbf{\Pi}}_a^{(2n)} = \sum_{m=0}^{n-2} (-1)^{m+n} \left(\frac{p_a}{\rho_a}\right)^{n-m-1} \frac{(n-1)!}{m!(n-m-1)!} \frac{(2n+3)!!}{(2m+5)!!} \tilde{\mathbf{\Pi}}_a^{(2m+2)}, \tag{209} $$

where the result is zero if the upper summation index is less than the lower summation index, yielding closures $X_a^{(3)} = 0$, $X&#x0302;<sup>(4)</sup> = 0$, and $\tilde{\mathbf{\Pi}}_a^{(2)} = 0$. The closures are summarized in Tables 3 and 4.

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, 2019b), it is claimed that Landau fluid closures are necessary to go beyond the fourth-order moment in the fluid hierarchy. This interpretation was obtained in the Chew, Goldberger, and Low (Chew et al. 1956; CGL)
<table><caption>**Table 4**
Similar to Table 3, but for Hermite Closures $\hat{h}_{ij}^{(2n)} = 0$</caption><thead><tr><th>Hermite Closures</th><th>Fluid Closures</th></tr></thead><tbody><tr><td>$\hat{h}_{ij}^{(2)} = 0$</td><td>$\Pi_{ij}^{(2)} = 0$</td></tr><tr><td>$\hat{h}_{ij}^{(4)} = 0$</td><td>$\Pi_{ij}^{(4)} = 7\frac{p}{\rho}\Pi_{ij}^{(2)}$</td></tr><tr><td>$\hat{h}_{ij}^{(6)} = 0$</td><td>$\Pi_{ij}^{(6)} = 18\frac{p}{\rho}\Pi_{ij}^{(4)} - 63\frac{p^2}{\rho^2}\Pi_{ij}^{(2)}$</td></tr><tr><td>$\hat{h}_{ij}^{(8)} = 0$</td><td>$\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)}$</td></tr></tbody></table>

Note. A general form for the closures corresponding to $\hat{h}_{ij}^{(2n)} = 0$ is given by (209).

framework for parallel moments by considering the closures at the last retained moment, $\tilde{X}_a^{(2n)} = 0$ and $X_{\parallel a}^{(2n+1)} = 0$. It was shown (see the detailed proof in Section 12.2 in Hunana et al. 2019b) that beyond the fourth-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 other words, beyond the fourth-order moment it is not possible to cut the fluid hierarchy by simply neglecting the next order moment with closures such as $X_a^{(5)} = 0$ or $\tilde{X}_a^{(6)} = 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 6). 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 $q_a$ and $\bar{\Pi}_a^{(2)}$ neglect all of the scalar perturbations $\tilde{X}_a^{(4)} = \dots = \tilde{X}_a^{(2n)} = 0$, together with neglecting the coupling between heat fluxes and stress tensors. In our formulation, this yields the system

$$ \frac{d_a X_a^{(2n+1)}}{dt} + \Omega_a \hat{b} \times X_a^{(2n+1)} + (2n+3)!! \frac{(n)}{3} \frac{p_a^n}{\rho_a^{n-1}} \nabla \left( \frac{p_a}{\rho_a} \right) = \bar{Q}_a^{(2n+1)} - \frac{(2n+3)!!}{3} \frac{p_a^n}{\rho_a^n} \mathbf{R}_a; \quad (210) $$

$$ \frac{d_a \bar{\Pi}_a^{(2n)}}{dt} + \Omega_a (\hat{b} \times \bar{\Pi}_a^{(2n)})^S + \frac{(2n+3)!!}{15} \frac{p_a^n}{\rho_a^{n-1}} \bar{W}_a = \bar{Q}_a^{(2n)} - \frac{\bar{I}}{3} Q_a^{(2n)}. \quad (211) $$

The closures (208), (209) are not required, because the equations are decoupled. We did not calculate the collisional contributions for higher-order moments, but in the semilinear approximation Equations (210)–(211) remain decoupled and represent two independent hierarchies. An essential feature of the Landau (or the Boltzmann) collisional operator is that the operator couples all of the heat fluxes together, and it also couples all of the stress tensors together. Thus, going higher and higher in the fluid hierarchy does not create new contributions in a quasistatic approximation, but yields increasingly precise transport coefficients for $q_a$ and $\bar{\Pi}_a^{(2)}$. Also, because the momentum exchange rates $R_a$ contain contributions from all of the heat fluxes $X_a^{(3)} \dots X_a^{(2n+1)}$, 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 decoupled system. A brief comparison of the various models is presented in Appendix I.

## 8.7. Inclusion of Gravity

We have not explicitly considered the force of gravity in our calculations in the appendix; nevertheless, its inclusion is trivial. With the gravitational acceleration **G** included, the Boltzmann equation reads

$$ \frac{\partial f_a}{\partial t} + \mathbf{v} \cdot \nabla f_a + \left[ \mathbf{G} + \frac{eZ_a}{m_a} (\mathbf{E} + \frac{1}{c} \mathbf{v} \times \mathbf{B}) \right] \cdot \nabla v f_a = C(f_a). \quad (212) $$

We use big **G** instead of small **g** to clearly distinguish it from the heat flux **q**. Gravity does not enter the collisional operator, and collisional integrals with the right-hand side are not affected. 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 the evolution equations for all other moments, analogously as the electric field does. An explicit
collisionless equation for a general nth-order moment with the electric field present is Equation (12.13) of Hunana et al. (2019b), for example. Because no Maxwell equations are used in deriving the fluid hierarchy, the presence of gravity can be accounted for by simply replacing

$$ \frac{eZ_a}{m_a} \mathbf{E} \rightarrow \mathbf{G} + \frac{eZ_a}{m_a} \mathbf{E}. \qquad (213) $$

Furthermore, such a hierarchy is not very useful, because the evolution equation for an nth-order moment is coupled with “n” momentum equations. Subtracting these momentum equations yields the 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 in the entire model is thus achieved trivially, by adding – **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 that contains gravity is the electric field Equation (96).

## 8.8. Precision of $m_e/m_i$ Expansions (Unmagnetized Proton–Electron Plasma)

The multifluid formulation is also an excellent tool for double-checking the precision of $m_e/m_i$ expansions. It is again possible to consider a one ion–electron plasma, but this time to 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 maintain $\nabla T_e \neq \nabla T_p$, however, because the gradients can be different. We first calculate the heat fluxes. For clarity, we are solving four coupled evolution equations, which are explicitly given in Appendix N; see Equations (N1)–(N4).

A 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 quasistatic approximation then yields the heat fluxes

$$ \begin{aligned} q_e &= [-3.159370\nabla T_e + 8.301 \times 10^{-6}\nabla T_p] \frac{p_e}{m_e \nu_{ep}} + 0.711046 p_e \delta u; \\ X_e^{(5)} &= [-110.5793\nabla T_e + 1.376 \times 10^{-3}\nabla T_p] \frac{p_e^2}{\rho_e m_e \nu_{ep}} + 18.78249 \frac{p_e^2}{\rho_e} \delta u; \\ q_p &= [-3.302411\nabla T_p + 0.2516 \times 10^{-3}\nabla T_e] \frac{p_p}{m_p \nu_{pp}} + 0.206535 \times 10^{-4} p_p \delta u; \\ X_p^{(5)} &= [-103.3984\nabla T_p + 0.7863 \times 10^{-2}\nabla T_e] \frac{p_p^2}{\rho_p m_p \nu_{pp}} + 0.646475 \times 10^{-3} \frac{p_p^2}{\rho_p} \delta u, \end{aligned} \qquad (214) $$

where $\delta u = (u_e - u_p)$. For the electron heat flux $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, as well as by using slightly different ratios of frequencies (a less-precise calculation, 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 $q_p$, the relatively large difference between the thermal conductivity 3.302 and the Braginskii self-collisional value 125/32 = 3.906 is caused by the proton–electron collisions. This is similarly the case for the $X_p^{(5)}$, where the self-collisional value is 2975/24 = 123.96. Calculating the coupled system exactly has a nice advantage, since one can calculate the momentum exchange rates in two different ways:

$$ \begin{aligned} \mathbf{R}_e &= \nu_{ep} \left\{ -\rho_e \delta \mathbf{u} + \frac{\mu_{ep}}{T_{ep}} \left[ V_{ep(1)} \mathbf{q}_e - V_{ep(2)} \frac{\rho_e}{\rho_p} \mathbf{q}_p \right] - \frac{3}{56} \left( \frac{\mu_{ep}}{T_{ep}} \right)^2 \left[ \mathbf{X}_e^{(5)} - \frac{\rho_e}{\rho_p} \mathbf{X}_p^{(5)} \right] \right\}; \\ \mathbf{R}_p &= \nu_{pe} \left\{ +\rho_p \delta \mathbf{u} + \frac{\mu_{ep}}{T_{ep}} \left[ V_{pe(1)} \mathbf{q}_p - V_{pe(2)} \frac{\rho_p}{\rho_e} \mathbf{q}_e \right] - \frac{3}{56} \left( \frac{\mu_{ep}}{T_{ep}} \right)^2 \left[ \mathbf{X}_p^{(5)} - \frac{\rho_p}{\rho_e} \mathbf{X}_e^{(5)} \right] \right\}, \end{aligned} \qquad (215) $$

and both options yield the same result:

$$ \mathbf{R}_e = - \mathbf{R}_p = - 0.711046 n_e \nabla T_e - 0.2065 \times 10^{-4} n_e \nabla T_p - 0.513306 \rho_e \nu_{ep} \delta u. \qquad (216) $$
The viscosities of the proton-electron plasma are (for clarity, we are solving four equations in four unknowns, explicitly given by (N5)–(N8))

$$
\begin{align}
\bar{\Pi}_e^{(2)} &= [-0.730622 \bar{\bar{W}}_e - 0.2800 \times 10^{-2} \bar{\bar{W}}_p] \frac{p_e}{\nu_{ei}}; \nonumber \\
\bar{\Pi}_e^{(4)} &= [-6.542519 \bar{\bar{W}}_e + 3.1509 \times 10^{-2} \bar{\bar{W}}_p] \frac{p_e^2}{\rho_e \nu_{ei}}; \nonumber \\
\bar{\Pi}_p^{(2)} &= [-0.892105 \bar{\bar{W}}_p - 0.4621 \times 10^{-4} \bar{\bar{W}}_e] \frac{p_p}{\nu_{pp}}; \nonumber \\
\bar{\Pi}_p^{(4)} &= [-7.250870 \bar{\bar{W}}_p - 0.3759 \times 10^{-3} \bar{\bar{W}}_e] \frac{p_p^2}{\rho_p \nu_{pp}}, \tag{217}
\end{align}
$$

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 other examples of coupling between the two species, and we calculate the heat fluxes and viscosities for protons and 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, or radiative transfer. Additionally, our approach shows that the coupling of stress tensors and heat fluxes should ideally be 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 have not used 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) loses its validity. A better approach then is 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 the expansions around an anisotropic bi-Maxwellian would look. The simplest anisotropic model is known as the CGL, after the pioneering work of Chew, Goldberger, and Low (Chew et al. 1956). The differences with our current approach start with the decomposition of the pressure tensor $p_{ij}^a$ defined in (A2), and the decomposition reads

$$
\text{isotropic:} \quad \bar{\bar{p}}_a = p_a \bar{\bar{I}} + \bar{\bar{\Pi}}_a^{(2)}; \tag{218}
$$

$$
\begin{align}
\text{anisotropic:} \quad \bar{\bar{p}}_a &= p_{\parallel a} \hat{b} \hat{b} + p_{\perp a} (\bar{\bar{I}} - \hat{b} \hat{b}) + \bar{\bar{\Pi}}_a^{(2)\text{CGL}} \nonumber \\
&= p_a \bar{\bar{I}} + (p_{\parallel a} - p_{\perp a}) \left( \hat{b} \hat{b} - \frac{\bar{\bar{I}}}{3} \right) + \bar{\bar{\Pi}}_a^{(2)\text{CGL}}, \tag{219}
\end{align}
$$

with scalar pressures

$$
p_{\parallel a} = \bar{\bar{p}}_a : \hat{b}\hat{b} = m_a \int c_{\parallel a}^2 f_a d^3 v; \quad p_{\perp a} = \bar{\bar{p}}_a : \bar{\bar{I}} / 2 = \frac{m_a}{2} \int |c_{\perp a}|^2 f_a d^3 v. \tag{220}
$$

Directly from the above definitions, the stress tensors satisfy

$$
\mathrm{Tr} \bar{\bar{\Pi}}_a^{(2)} = \mathrm{Tr} \bar{\bar{\Pi}}_a^{(2)\mathrm{CGL}} = 0; \quad \bar{\bar{\Pi}}_a^{(2)} : \hat{b}\hat{b} \neq 0; \quad \bar{\bar{\Pi}}_a^{(2)\mathrm{CGL}} : \hat{b}\hat{b} = 0, \tag{221}
$$

and while $\bar{\bar{\Pi}}_a^{(2)}$ has five independent components, $\bar{\bar{\Pi}}_a^{(2)\text{CGL}}$ has only four.

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:

$$
\mathbf{S}_a^{\parallel} = \bar{\bar{q}}_a : \hat{b}\hat{b} = m_a \int c_{\parallel a}^2 \mathbf{c}_a f_a d^3 v; \quad \mathbf{S}_a^{\perp} = \bar{\bar{q}}_a : \bar{\bar{I}} / 2 = \frac{m_a}{2} \int |c_{\perp a}|^2 \mathbf{c}_a f_a d^3 v. \tag{222}
$$
These heat flux vectors are further split by projecting them along the $\hat{b}$, which defines the gyrotropic (scalar) heat fluxes $q_{\parallel a}$ and $q_{\perp a}$, and the perpendicular projection defines the nongyrotropic heat flux vectors $S_{\perp a}^{\parallel}$ and $S_{\perp a}^{\perp}$, according to

$$ 
\begin{aligned}
q_{\parallel a} &= \hat{b} \cdot S_a^{\parallel} = m_a \int c_{\parallel a}^2 c_{\parallel a} f_a d^3 v; & q_{\perp a} &= \hat{b} \cdot S_a^{\perp} = \frac{m_a}{2} \int |c_{\perp a}|^2 c_{\parallel a} f_a d^3 v; \\
S_{\perp a}^{\parallel} &= \bar{\mathbf{I}}_{\perp} \cdot S_a^{\parallel} = m_a \int c_{\parallel a}^2 c_{\perp a} f_a d^3 v; & S_{\perp a}^{\perp} &= \bar{\mathbf{I}}_{\perp} \cdot S_a^{\perp} = \frac{m_a}{2} \int |c_{\perp a}|^2 c_{\perp a} f_a d^3 v.
\end{aligned} 
\quad (223) $$

The following relations then hold $S_a^{\parallel} = q_{\parallel a} \hat{b} + S_a^{\perp}$ and $S_a^{\perp} = q_{\perp a} \bar{\mathbf{I}} + S_a^{\perp}$, together with $\hat{b} \cdot S_{\perp a}^{\parallel} = 0$ and $\hat{b} \cdot S_{\perp a}^{\perp} = 0$. The two different decompositions of the entire heat flux tensor then read

$$ \text{isotropic:} \qquad \bar{\bar{q}}_a = \frac{2}{5} [q_a \bar{\mathbf{I}}]^S + \bar{\bar{\sigma}}_a', \quad (224) $$

$$ \begin{aligned}
\text{anisotropic:} \quad & \bar{\bar{q}}_a = q_{\parallel a} \hat{b} \hat{b} \hat{b} + q_{\perp a} [\hat{b} \bar{\mathbf{I}}]^{\parallel} + [S_{\perp a}^{\parallel} \hat{b} \hat{b}]^{\parallel} + \frac{1}{2} [S_{\perp a}^{\perp} \bar{\mathbf{I}}]^{\parallel} + \bar{\bar{\sigma}}_a; \\
& = q_{\perp a} [\hat{b} \bar{\mathbf{I}}]^{\parallel} + (q_{\parallel a} - 3q_{\perp a}) \hat{b} \hat{b} \hat{b} + \frac{1}{2} [S_{\perp a}^{\perp} \bar{\mathbf{I}}]^{\parallel} + \left[ \left( S_{\perp a}^{\parallel} - \frac{S_{\perp a}^{\perp}}{2} \right) \hat{b} \hat{b} \right]^{\parallel} + \bar{\bar{\sigma}}_a,
\end{aligned} 
\quad (225) $$

where both $\bar{\bar{\sigma}}_a'$ and $\bar{\bar{\sigma}}_a$ are traceless. Neglecting these traceless contributions, the isotropic approach accounts for three (out of 10) scalar components of $\bar{\bar{q}}_a$, and represents a 13-moment model (one density, three velocity, one scalar pressure, five stress tensor, and three heat flux $q_a$ components). The anisotropic approach accounts for six scalar components of $\bar{\bar{q}}_a$ and represents a 16-moment model, described by 16 scalar evolution Equations. (One density, three velocity, two scalar pressures, four stress tensor components, and three for each heat flux vector $S_a^{\parallel}$ and $S_a^{\perp}$.) Unfortunately, such a complicated decomposition of the heat flux is necessary in an arbitrary collisional regime, and we only used decomposition (224). For clarity, the direct relation with the usual heat flux vector $q_a$ reads

$$ q_a = \frac{1}{2} S_a^{\parallel} + S_a^{\perp} = \left( \frac{1}{2} q_{\parallel a} + q_{\perp a} \right) \hat{b} + \frac{1}{2} S_{\perp a}^{\parallel} + S_{\perp a}^{\perp}. \quad (226) $$

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

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), and 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 fourth- and fifth-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) and Cuperman & Dryer (1985) considered what seems like a mixture of anisotropic and isotropic moments, with anisotropic temperatures, an isotropic heat flux vector, and the parallel (anisotropic) perturbation of the fourth-order moment (which we call $\tilde{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, the collisionless linear kinetic theory expanded around a bi-Maxwellian with mean zero drifts in gyrotropic limit yields in Fourier space a 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)}} \left[ 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})} \right] + \Phi \frac{eZ_a}{T_{\parallel a}^{(0)}} \frac{k_{\parallel} v_{\parallel}}{(\omega - k_{\parallel} v_{\parallel})}, \quad (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_{\text{th}\parallel a} n_a^{(0)} T_{\parallel a}^{(0)} \text{sign}(k_{\parallel})(\zeta_a + 2\zeta_a^3 R(\zeta_a) - 3\zeta_a R(\zeta_a)) \left[ \frac{B_{\parallel}^{(1)} T_{\perp a}^{(0)}}{B_0 T_{\parallel a}^{(0)}} + \Phi \frac{eZ_a}{T_{\parallel a}^{(0)}} \right], \quad (228) $$

with variable $\zeta_a = \omega/(|k_{\parallel}|v_{\text{th}\parallel a})$, parallel thermal speed $v_{\text{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 + \text{erf}(i\zeta_a)]$. Such a kinetic answer can be expressed in fluid variables by searching for
<figure><img src="image_4.png" /><figcaption>**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: the imaginary part of $\hat{q}_{\parallel a}$. Right: the real part of $\hat{q}_{\parallel a}$. The colors are described in the text. Our Braginskii-type models do not contain these collisionless heat fluxes.</figcaption></figure>

Landau fluid closures, by, for example, replacing the $R(\zeta_a)$ with its three-pole Padé approximants:

$$R_{3,2}(\zeta_a) = \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} \qquad (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:

$$R_{3,2}(\zeta_a): \qquad q_{\parallel a}^{(1)} = -i \frac{2}{\sqrt{\pi}} n_a^{(0)} v_{\text{th}\parallel a} \text{sign}(k_{\parallel}) T_{\parallel a}^{(1)}; \qquad (230)$$

$$R_{3,1}(\zeta_a): \qquad q_{\parallel a}^{(1)} = \frac{3\pi - 8}{4 - \pi} p_{\parallel a}^{(0)} u_{\parallel a}^{(1)} - i \frac{\sqrt{\pi}}{4 - \pi} n_a^{(0)} v_{\text{th}\parallel a} \text{sign}(k_{\parallel}) T_{\parallel a}^{(1)}, \qquad (231)$$

where $T_{\parallel a}^{(1)}$ is the perturbed temperature and $u_{\parallel a}^{(1)}$ is the perturbed fluid velocity (a mean value of $u_{\parallel a}^{(0)} = 0$ is assumed). The heat flux closure (230) was obtained by Hammett & Perkins (1990) and 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

$$R_{3,2}(\zeta_a): \qquad q_{\parallel a}(z) = - \frac{2}{\pi^{3/2}} n_a^{(0)} v_{\text{th}\parallel a} \text{V.P.} \int_0^{\infty} \frac{T_{\parallel a}^{(1)}(z+z') - T_{\parallel a}^{(1)}(z-z')}{z'} dz'; \qquad (232)$$

$$R_{3,1}(\zeta_a): \qquad q_{\parallel a}(z) = \frac{3\pi - 8}{4 - \pi} p_{\parallel a}^{(0)} u_{\parallel a}^{(1)} - \frac{n_a^{(0)} v_{\text{th}\parallel a}}{\sqrt{\pi} (4 - \pi)} \text{V.P.} \int_0^{\infty} \frac{T_{\parallel a}^{(1)}(z+z') - T_{\parallel a}^{(1)}(z-z')}{z'} dz', \qquad (233)$$

where the nonlocality presents itself as an integral over the entire magnetic field line, where the temperatures everywhere along that field line matter in order to determine the heat flux at a specific spatial point. Note that the thermal part of (233) is almost two times larger than (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 quasistatic approximation. However, as is well known, these expressions are not very precise with respect to kinetic theory. For example, the precision can easily be compared by plotting the normalized heat fluxes $\hat{q}_{\parallel a} = \zeta_a + 2\zeta_a^3 R(\zeta_a) - 3\zeta_a R(\zeta_a)$, which are shown in Figure 1. A weakly damped regime with a real-valued $\zeta_a$ is considered. The left-hand panel shows the imaginary part of $\hat{q}_{\parallel a}$, and the right-hand panel shows the real part of $\hat{q}_{\parallel a}$. The exact kinetic heat flux is shown by the solid black line, the heat flux $R_{3,2}$ is shown by the dashed magenta line, and the heat flux $R_{3,1}$ is shown by the dashed cyan line. For comparison, higher-order fluid models with approximants $R_{5,3}$ (the dotted blue line) and $R_{7,5}$ (the dashed red line) are shown as well (see Equations (A11) and (A38) in Hunana et al. 2019a). The $R_{5,3}$ model represents a dynamic closure at the fourth-order moment, and the $R_{7,5}$ model represents a dynamic closure at the sixth-order moment, given by Equations (5) and (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 quasistatic 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$). Regimes $\zeta_a \ll 1$ can be viewed as isothermal, and regimes $\zeta_a \gg 1$ can be viewed as adiabatic. In the left-hand panel of Figure 1, the $R_{3,1}$ is more precise up to roughly $\zeta_a = 2.3$, and in the right-hand panel up to $\zeta_a = 1.6$. For larger $\zeta_a$ values than shown, the $R_{3,1}$ heat flux converges much more slowly to the correct zero values than the $R_{3,2}$, especially for the real part.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Normalized damping rate $\zeta_i = \omega_i/(|k_{\parallel}|v_{th_{\parallel}}|$) for a parallel propagating ion sound wave as a function of $k_{\parallel}\lambda_{mfp}$, where $\lambda_{mfp}$ is a mean-free path. Red line: Braginskii-type (isotropic) 13-moment model with quasistatic 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 the quasistatic heat fluxes of Snyder et al. (1997).</figcaption></figure>

The major obstacle to precision for the quasistatic heat fluxes of the Landau fluid models actually comes from the perpendicular heat flux $q_{\perp a}$ (which is along the $\hat{b}$), because only a closure of Snyder et al. (1997) with a crude Padé approximant $R_l(\zeta_a) = 1/(1 - i\sqrt{\pi}\zeta_a)$ is available. As a consequence, for large $\zeta_a$ values, the quasistatic 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 quasistatic $q_{\perp a}$ and consider its evolution equation, with a closure at the fourth-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 by Scudder 2021, for example). The closures are constructed universally for all of 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 nonlinear 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$ that 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 three-component plasma, which allows the net current to be zero by the ion-ion (or the electron-ion and electron-electron) acoustic instability, see Gary (1993, pp. 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 the 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 highly collisional to weakly collisional ones. In particular, let us consider the damping of a monochromatic ion sound wave of parallel wavenumber $k_{\parallel}$ 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., the discussion in Hunana et al. 2019a, p. 73), but it allows one to simplify the presentation with the goal of describing the general behavior, not providing 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. The x-axis shows $k_{\parallel}\lambda_{mfp}$, where $\lambda_{mfp} = v_{th_{\parallel}}/\nu$ is the ion mean-free path and $\nu$ is the collisional frequency, so that $k_{\parallel}\lambda_{mfp} \ll 1$ represents a highly collisional regime and $k_{\parallel}\lambda_{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_{th_{\parallel}})$. The usual isotropic 13-moment model (the green line) and the anisotropic 16-moment model (the blue line), with all of 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}$, and $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}$, and $q_{\perp}$ (we consider the 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_{\parallel}\lambda_{mfp} \sim 0.5 - 1$, and converge toward a zero damping rate in the collisionless regime (with only a small shift in $k_{\parallel}\lambda_{mfp}$ between them). In contrast, the red line, corresponding to the 13-moment model with the $\Pi_{zz}$ and $q_z$ taken in the quasistatic approximation, shows that the damping rate does not reach a maximum and instead continues to increase in a weakly collisional regime, while around $k_{\parallel}\lambda_{mfp} \sim 6.3$ the sound mode stops existing (it becomes nonpropagating, with zero real frequency). This is consequence of the quasistatic 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