damping is preferred against a quantity that blows up in a weakly collisional regime, all three models are technically incorrect, because the Landau damping provides a significant contribution to the damping rate as the plasma becomes weakly collisional. To illustrate the importance of the Landau damping, the magenta line displays the damping rate obtained with a Landau fluid model that contains the evolution equations for $\rho$, $u_z$, $p_\parallel$, and $p_\perp$, but where the quasistatic $q_\parallel$ and $q_\perp$ are given by the collisionally modified 3+1 closures of Snyder et al. (1997), i.e., their Equations (48)–(49), which for the isotropic mean temperatures $T_{\parallel a}^{(0)} = T_{\perp a}^{(0)}$ considered here are equivalent to (our thermal speed contains a factor of 2, which is not the case in that paper)

$$q_{\parallel a}^{(1)} = - \frac{\frac{4}{3\pi - 8} n_a^{(0)} v_{\text{th}\parallel a}^2}{\bar{\nu}_a + \frac{2\sqrt{\pi}}{3\pi - 8} v_{\text{th}\parallel a} |k_\parallel|} ik_\parallel T_{\parallel a}^{(1)}; \qquad q_{\perp a}^{(1)} = - \frac{\frac{1}{2} n_a^{(0)} v_{\text{th}\parallel a}^2}{\bar{\nu}_a + \frac{\sqrt{\pi}}{2} v_{\text{th}\parallel a} |k_\parallel|} ik_\parallel T_{\perp a}^{(1)} \quad (234)$$

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

### 8.9.4. Large Gradients and Large Drifts

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

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

$$q_e^{\text{Hollweg}} = \frac{3}{2} p_e \mathbf{u}_{\text{sw}} \alpha, \quad (235)$$

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

### 8.9.5. Comments on the Positivity of the Perturbed Distribution Function

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

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

### 8.10. Conclusions

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

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

$$ \mathbf{q}_e = \frac{\mathbf{X}_e^{(3)}}{2} = -3.2 \frac{p_e}{m_e \nu_{ei}} \nabla T_e; \quad \mathbf{X}_e^{(5)} = -110.7 \frac{p_e^2}{\rho_e m_e \nu_{ei}} \nabla T_e, \qquad (236) $$

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

$$ \tilde{X}_e^{(4)} = +83.8 \frac{p_e^2}{\nu_{ei}^2 \rho_e m_e} \nabla^2 T_e. \qquad (237) $$

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

$$ \text{Burgers-Schunk:} \quad \mathbf{R}_e^T = + \frac{3}{5} \frac{\rho_e}{p_e} \nu_{ei} \mathbf{q}_e; \qquad (238) $$

$$ \text{Killie et al.:} \quad \mathbf{R}_e^T = + \frac{6}{35} \frac{\rho_e}{p_e} \nu_{ei} \mathbf{q}_e; \qquad (239) $$

$$ \text{present paper:} \quad \mathbf{R}_e^T = + \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)}. \qquad (240) $$

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

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

## Appendix A
General Evolution Equations

We consider the Boltzmann equation (in cgs units)

$$ \frac{\partial f_a}{\partial t} + \mathbf{v} \cdot \nabla f_a + \frac{eZ_a}{m_a} \left( \mathbf{E} + \frac{1}{c} \mathbf{v} \times \mathbf{B} \right) \cdot \nabla_v f_a = C(f_a), \quad (A1) $$

where “a” is a species index and $C(f_a) = \sum_b C_{ab}(f_a, f_b)$ is the Landau collisional operator, so Equation (A1) is called the Landau equation. One defines the usual number density $n_a = \int f_a d^3v$, density $\rho_a = m_a n_a$, fluid velocity $\mathbf{u}_a = (1/n_a) \int v f_a d^3v$, and fluctuating velocity $\mathbf{c}_a = \mathbf{v} - \mathbf{u}_a$, and further defines the pressure tensor $\bar{\mathbf{p}}_a$, heat flux tensor $\bar{\mathbf{q}}_a$, fourth-order moment $\bar{\mathbf{r}}_a$, and fifth-order and sixth-order moments $\bar{\mathbf{X}}_a^{(5)}$, $\bar{\mathbf{X}}_a^{(6)}$, respectively, according to

$$ \bar{\mathbf{p}}_a = m_a \int \mathbf{c}_a \mathbf{c}_a f_a d^3v; \quad \bar{\mathbf{q}}_a = m_a \int \mathbf{c}_a \mathbf{c}_a f_a d^3v; \quad \bar{\mathbf{r}}_a = m_a \int \mathbf{c}_a \mathbf{c}_a f_a d^3v; \quad (A2) $$

$$ \bar{\mathbf{X}}_a^{(5)} = m_a \int \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a f_a d^3v; \quad \bar{\mathbf{X}}_a^{(6)} = m_a \int \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a \mathbf{c}_a f_a d^3v. \quad (A3) $$

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

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

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

$$ \frac{\partial n_a}{\partial t} + \nabla \cdot (n_a \mathbf{u}_a) = 0; \qquad (A4) $$

$$ \frac{\partial \mathbf{u}_a}{\partial t} + \mathbf{u}_a \cdot \nabla \mathbf{u}_a + \frac{1}{\rho_a} \nabla \cdot \bar{\mathbf{p}}_a - \frac{eZ_a}{m_a} \left( \mathbf{E} + \frac{1}{c} \mathbf{u}_a \times \mathbf{B} \right) = \frac{\mathbf{R}_a}{\rho_a}; \qquad (A5) $$

$$ \frac{\partial \bar{\mathbf{p}}_a}{\partial t} + \nabla \cdot (\bar{\mathbf{q}}_a + \mathbf{u}_a \bar{\mathbf{p}}_a) + [\bar{\mathbf{p}}_a \cdot \nabla \mathbf{u}_a + \Omega_a \hat{\mathbf{b}} \times \bar{\mathbf{p}}_a]^S = \bar{\mathbf{Q}}_a^{(2)}; \qquad (A6) $$

$$ \frac{\partial \bar{\mathbf{q}}_a}{\partial t} + \nabla \cdot (\bar{\mathbf{r}}_a + \mathbf{u}_a \bar{\mathbf{q}}_a) + \left[ \bar{\mathbf{q}}_a \cdot \nabla \mathbf{u}_a + \Omega_a \hat{\mathbf{b}} \times \bar{\mathbf{q}}_a - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \bar{\mathbf{p}}_a \right]^S = \bar{\mathbf{Q}}_a^{(3)} - \frac{1}{\rho_a} [\mathbf{R}_a \bar{\mathbf{p}}_a]^S; \qquad (A7) $$

$$ \frac{\partial \bar{\mathbf{r}}_a}{\partial t} + \nabla \cdot (\bar{\mathbf{X}}_a^{(5)} + \mathbf{u}_a \bar{\mathbf{r}}_a) + \left[ \bar{\mathbf{r}}_a \cdot \nabla \mathbf{u}_a + \Omega_a \hat{\mathbf{b}} \times \bar{\mathbf{r}}_a - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \bar{\mathbf{q}}_a \right]^S = \bar{\mathbf{Q}}_a^{(4)} - \frac{1}{\rho_a} [\mathbf{R}_a \bar{\mathbf{q}}_a]^S; \qquad (A8) $$

$$ \frac{\partial \bar{\mathbf{X}}_a^{(5)}}{\partial t} + \nabla \cdot (\bar{\mathbf{X}}_a^{(6)} + \mathbf{u}_a \bar{\mathbf{X}}_a^{(5)}) + \left[ \bar{\mathbf{X}}_a^{(5)} \cdot \nabla \mathbf{u}_a + \Omega_a \hat{\mathbf{b}} \times \bar{\mathbf{X}}_a^{(5)} - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \bar{\mathbf{r}}_a \right]^S = \bar{\mathbf{Q}}_a^{(5)} - \frac{1}{\rho_a} [\mathbf{R}_a \bar{\mathbf{r}}_a]^S, \qquad (A9) $$

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

$$ X_{r_1 r_2 \dots r_n}^{(n)} = m \int c_{r_1} c_{r_2} \dots c_{r_n} f d^3 v; \qquad Q_{r_1 r_2 \dots r_n}^{(n)} = m \int c_{r_1} c_{r_2} \dots c_{r_n} C(f) d^3 v, \qquad (A10) $$

together with a symmetric operator “S” that cycles around all of its indices:

$$ [X^{(n)}]_{r_1 r_2 \dots r_n}^S = X_{r_1 r_2 \dots r_n}^{(n)} + X_{r_2 r_3 \dots r_n}^{(n)} + X_{r_3 \dots r_n r_1 r_2}^{(n)} + \dots + X_{r_n r_1 r_2 \dots r_{n-1}}^{(n)}, \qquad (A11) $$

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

$$ \frac{\partial}{\partial t} \bar{\mathbf{X}}_a^{(n)} + \nabla \cdot (\bar{\mathbf{X}}_a^{(n+1)} + \mathbf{u}_a \bar{\mathbf{X}}_a^{(n)}) + \left[ \bar{\mathbf{X}}_a^{(n)} \cdot \nabla \mathbf{u}_a + \Omega_a \hat{\mathbf{b}} \times \bar{\mathbf{X}}_a^{(n)} - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \bar{\mathbf{X}}_a^{(n-1)} \right]^S = \bar{\mathbf{Q}}_a^{(n)} - \frac{1}{\rho_a} [\mathbf{R}_a \bar{\mathbf{X}}_a^{(n-1)}]^S, \qquad (A12) $$

valid for $n \ge 2$. The left-hand side of (A12) is equal to the collisionless Equation (12.16) of Hunana et al. (2019b). Evolution Equations (A6)–(A9) can then easily be obtained by the evaluation of (A12). Note that definition (A10) yields $\bar{\mathbf{X}}^{(2)} = \bar{\mathbf{p}}$, $\bar{\mathbf{X}}^{(3)} = \bar{\mathbf{q}}$, and $\bar{\mathbf{X}}^{(4)} = \bar{\mathbf{r}}$, but $\bar{\mathbf{X}}^{(1)} = 0$.

As has already been pointed out by Grad (1949a, 1949b), who developed the moment approach considering rarefied gases, because fluid moments are symmetric in all of their indices, a general nth-order moment $\bar{\mathbf{X}}_a^{(n)}$ contains $(\frac{n+2}{n}) = (n+1)(n+2)/2$ distinct (scalar) components. So the density has 1, the velocity has 3, the pressure tensor has 6, the heat flux tensor has 10, $\bar{\mathbf{X}}_a^{(4)}$ has 15, and $\bar{\mathbf{X}}_a^{(5)}$ has 21 scalar components. The system (A4)–(A9) thus represents a 56-moment model.

## Appendix B
Tensorial Hermite Decomposition

In the famous work of Grad (1949a, 1949b, 1958), the so-called tensorial Hermite decomposition is used, which is a generalization of the one-dimensional version. The one-dimensional Hermite polynomials of order “m” are defined as

$$ H^{(m)}(x) = (-1)^m e^{\frac{x^2}{2}} \frac{d^m}{dx^m} e^{-\frac{x^2}{2}}, \qquad (B1) $$

and evaluated step by step as $H^{(0)} = 1$, $H^{(1)} = x$, $H^{(2)} = x^2 - 1$, $H^{(3)} = x^3 - 3x$, $H^{(4)} = x^4 - 6x^2 + 3$, and $H^{(5)} = x^5 - 10x^3 + 15x$. So polynomials of even order contain only even powers of x and polynomials of odd order contain only odd powers of x. These polynomials are orthogonal to each other by

$$ \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} H^{(n)}(x) H^{(m)}(x) e^{-\frac{x^2}{2}} dx = n! \delta_{nm}. \qquad (B2) $$

Note that the “weight” $\exp(-x^2/2)$ was used by Grad instead of the quantum-mechanical $\exp(-x^2)$. Of course, it is important to use the correct weights with both classes of Hermite polynomials. Curiously, if the weight is accidentally mismatched (i.e., by using $\exp(-x^2)$ in our (B2) or $\exp(-x^2/2)$ in the quantum version), in addition to naturally wrong numerical constants, the even–even and odd–odd couples of polynomials are not orthogonal any more! The generalization to tensors for the isotropic Maxwellian distribution
reads

$$ \tilde{H}_{r_1, r_2, \dots, r_m}^{(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_m}} e^{-\frac{\tilde{c}^2}{2}}. \quad (B3) $$

We use the same notation as Balescu (1988), where *reducible* Hermite polynomials are denoted with tilde, and *irreducible* polynomials have no tilde. We have added tilde on normalized $\tilde{c}$ to make transitioning to usual fluid moments straightforward. Then explicit evaluation step by step gives

$$ \begin{aligned} \tilde{H}^{(0)}(\tilde{c}) &= 1; \\ \tilde{H}_i^{(1)}(\tilde{c}) &= \tilde{c}_i; \\ \tilde{H}_{ij}^{(2)}(\tilde{c}) &= \tilde{c}_i \tilde{c}_j - \delta_{ij}; \\ \tilde{H}_{ijk}^{(3)}(\tilde{c}) &= \tilde{c}_i \tilde{c}_j \tilde{c}_k - (\delta_{ij} \tilde{c}_k + \delta_{jk} \tilde{c}_i + \delta_{ik} \tilde{c}_j); \\ \tilde{H}_{ijkl}^{(4)}(\tilde{c}) &= \tilde{c}_i \tilde{c}_j \tilde{c}_k \tilde{c}_l - (\delta_{ij} \tilde{c}_k \tilde{c}_l + \delta_{jk} \tilde{c}_i \tilde{c}_l + \delta_{kl} \tilde{c}_i \tilde{c}_j + \delta_{li} \tilde{c}_j \tilde{c}_k + \delta_{ik} \tilde{c}_j \tilde{c}_l + \delta_{jl} \tilde{c}_i \tilde{c}_k) \\ &\quad + \delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}, \end{aligned} \quad (B4) $$

and quickly starts to grow:

$$ \begin{aligned} \tilde{H}_{ijklm}^{(5)}(\tilde{c}) ={}& \tilde{c}_i \tilde{c}_j \tilde{c}_k \tilde{c}_l \tilde{c}_m - (\delta_{ij} \tilde{c}_k \tilde{c}_l \tilde{c}_m + \delta_{jk} \tilde{c}_l \tilde{c}_i \tilde{c}_m + \delta_{kl} \tilde{c}_i \tilde{c}_j \tilde{c}_m + \delta_{li} \tilde{c}_j \tilde{c}_k \tilde{c}_m + \delta_{ik} \tilde{c}_j \tilde{c}_l \tilde{c}_m \\ & + \delta_{jl} \tilde{c}_i \tilde{c}_k \tilde{c}_m + \delta_{im} \tilde{c}_j \tilde{c}_k \tilde{c}_l + \delta_{jm} \tilde{c}_i \tilde{c}_j \tilde{c}_l + \delta_{km} \tilde{c}_i \tilde{c}_j \tilde{c}_l + \delta_{lm} \tilde{c}_i \tilde{c}_j \tilde{c}_k) \\ & + \delta_{ij} \delta_{kl} \tilde{c}_m + \delta_{ik} \delta_{jl} \tilde{c}_m + \delta_{il} \delta_{jk} \tilde{c}_m + \delta_{ij} \delta_{km} \tilde{c}_l + \delta_{ij} \delta_{lm} \tilde{c}_k \\ & + \delta_{jk} \delta_{lm} \tilde{c}_i + \delta_{jk} \delta_{im} \tilde{c}_l + \delta_{kl} \delta_{im} \tilde{c}_j + \delta_{kl} \delta_{jm} \tilde{c}_i + \delta_{li} \delta_{jm} \tilde{c}_k \\ & + \delta_{li} \delta_{km} \tilde{c}_j + \delta_{ik} \delta_{jm} \tilde{c}_l + \delta_{ik} \delta_{lm} \tilde{c}_j + \delta_{jl} \delta_{im} \tilde{c}_k + \delta_{jl} \delta_{km} \tilde{c}_i. \end{aligned} \quad (B5) $$

The choice of Grad with exp(-x²/2) has a great benefit, because no numerical constants are present in the entire hierarchy of Hermite polynomials, which is not the case for the weight exp(-x²). Here, numerical factors appear only after one applies contractions (traces) at the above expressions. Similar to the one-dimensional case, polynomials of even order contain only terms with even numbers of velocities $\tilde{c}$, and polynomials of odd order contain only terms with odd numbers of $\tilde{c}$.

For a Maxwellian distribution, the normalized velocity is

$$ \tilde{c} = \sqrt{\frac{m_a}{T_a}} (v - u_a) = \sqrt{\frac{m_a}{T_a}} c, \quad (B6) $$

where, for simplicity, we suppress the writing of species index “a” for velocity c in the expressions that follow, and for many other variables as well (the Hermite decomposition is done independently for each species, and the species variable “a” just makes the expressions more blurry). It is possible to work in both normalized and physical units. The entire distribution function can be written as

$$ f_a = f_a^{(0)} (1 + \chi_a) = n_a \left(\frac{m_a}{T_a}\right)^{3/2} \phi^{(0)}(1 + \chi_a); \quad \phi^{(0)} = \frac{e^{-\frac{\tilde{c}^2}{2}}}{(2\pi)^{3/2}}, \quad (B7) $$

where $\chi_a$ represents the wanted perturbation. One can go quickly between the physical and normalized units by

$$ \int f_a(c) d^3 c = n_a \int \phi^{(0)}(1 + \chi_a(\tilde{c})) d^3 \tilde{c}. \quad (B8) $$

The tensorial polynomials are again orthogonal to each other, where, by using “weight” $\phi^{(0)}$:

$$ \begin{aligned} \int \phi^{(0)} \tilde{H}^{(0)} \tilde{H}^{(0)} d^3 \tilde{c} &= 1; \\ \int \phi^{(0)} \tilde{H}_i^{(1)} \tilde{H}_j^{(1)} d^3 \tilde{c} &= \delta_{ij}; \\ \int \phi^{(0)} \tilde{H}_{ij}^{(2)} \tilde{H}_{kl}^{(2)} d^3 \tilde{c} &= \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}; \\ \int \phi^{(0)} \tilde{H}_{r_1 r_2 r_3}^{(3)} \tilde{H}_{s_1 s_2 s_3}^{(3)} d^3 \tilde{c} &= \delta_{r_1 s_1} \delta_{r_2 s_2} \delta_{r_3 s_3} + \delta_{r_1 s_1} \delta_{r_2 s_3} \delta_{r_3 s_2} + \delta_{r_1 s_2} \delta_{r_2 s_1} \delta_{r_3 s_3} \\ &\quad + \delta_{r_1 s_2} \delta_{r_2 s_3} \delta_{r_3 s_1} + \delta_{r_1 s_3} \delta_{r_2 s_1} \delta_{r_3 s_2} + \delta_{r_1 s_3} \delta_{r_2 s_2} \delta_{r_3 s_1}, \end{aligned} \quad (B9) $$
and the expressions quickly become long:

$$ \int \phi^{(0)} \tilde{H}_{r_1 r_2 r_3 r_4}^{(4)} \tilde{H}_{s_1 s_2 s_3 s_4}^{(4)} d^3 \tilde{c} = + \delta_{r_1 s_1} \delta_{r_2 s_2} \delta_{r_3 s_3} \delta_{r_4 s_4} + \delta_{r_1 s_1} \delta_{r_2 s_2} \delta_{r_3 s_4} \delta_{r_4 s_3} + \delta_{r_1 s_1} \delta_{r_2 s_3} \delta_{r_3 s_2} \delta_{r_4 s_4} \\ + \delta_{r_1 s_1} \delta_{r_2 s_3} \delta_{r_3 s_4} \delta_{r_4 s_2} + \delta_{r_1 s_1} \delta_{r_2 s_4} \delta_{r_3 s_2} \delta_{r_4 s_3} + \delta_{r_1 s_1} \delta_{r_2 s_4} \delta_{r_3 s_3} \delta_{r_4 s_2} \\ + \delta_{r_1 s_2} \delta_{r_2 s_1} \delta_{r_3 s_3} \delta_{r_4 s_4} + \delta_{r_1 s_2} \delta_{r_2 s_1} \delta_{r_3 s_4} \delta_{r_4 s_3} + \delta_{r_1 s_2} \delta_{r_2 s_3} \delta_{r_3 s_1} \delta_{r_4 s_4} \\ + \delta_{r_1 s_2} \delta_{r_2 s_3} \delta_{r_3 s_4} \delta_{r_4 s_1} + \delta_{r_1 s_2} \delta_{r_2 s_4} \delta_{r_3 s_1} \delta_{r_4 s_3} + \delta_{r_1 s_2} \delta_{r_2 s_4} \delta_{r_3 s_3} \delta_{r_4 s_1} \\ + \delta_{r_1 s_3} \delta_{r_2 s_1} \delta_{r_3 s_2} \delta_{r_4 s_4} + \delta_{r_1 s_3} \delta_{r_2 s_1} \delta_{r_3 s_4} \delta_{r_4 s_2} + \delta_{r_1 s_3} \delta_{r_2 s_2} \delta_{r_3 s_1} \delta_{r_4 s_4} \\ + \delta_{r_1 s_3} \delta_{r_2 s_2} \delta_{r_3 s_4} \delta_{r_4 s_1} + \delta_{r_1 s_3} \delta_{r_2 s_4} \delta_{r_3 s_1} \delta_{r_4 s_2} + \delta_{r_1 s_3} \delta_{r_2 s_4} \delta_{r_3 s_3} \delta_{r_4 s_1} \\ + \delta_{r_1 s_4} \delta_{r_2 s_1} \delta_{r_3 s_2} \delta_{r_4 s_3} + \delta_{r_1 s_4} \delta_{r_2 s_1} \delta_{r_3 s_3} \delta_{r_4 s_2} + \delta_{r_1 s_4} \delta_{r_2 s_2} \delta_{r_3 s_1} \delta_{r_4 s_3} \\ + \delta_{r_1 s_4} \delta_{r_2 s_2} \delta_{r_3 s_3} \delta_{r_4 s_1} + \delta_{r_1 s_4} \delta_{r_2 s_3} \delta_{r_3 s_1} \delta_{r_4 s_2} + \delta_{r_1 s_4} \delta_{r_2 s_3} \delta_{r_3 s_2} \delta_{r_4 s_1}. \qquad (B10) $$

The general orthogonality can be written by introducing the multi-indices $r = r_1...r_n$ and $s = s_1...s_m$:

$$ \int \phi^{(0)} \tilde{H}_{r}^{(n)} \tilde{H}_{s}^{(m)} d^3 \tilde{c} = \delta_{mn} \delta_{rs}^{(n)}, \qquad (B11) $$

where the new symbol $\delta_{rs}^{(n)}$ is equal to one, if the indices $r_1...r_n$ are a permutation of $s_1...s_n$, but otherwise is zero. In other words, for $n = m$, the right-hand side contains $n!$ terms, where each of these terms has the form $\delta_{r_1s_1} \delta_{r_2s_2} ... \delta_{r_n s_n}$, and to calculate the other terms it is necessary to keep the $r$-indices fixed, and do all the possible permutations with $s$-indices (or vice versa). A particular case of (B11) reads

$$ m \neq 0: \qquad \int \phi^{(0)} \tilde{H}_{s}^{(m)} d^3 \tilde{c} = 0, \qquad (B12) $$

i.e., the integral over a single Hermite polynomial with weight $\phi^{(0)}$ is zero.

The goal of the Hermite expansion is to find the perturbation of the distribution function $\chi_a$ in (B7). For the most general decomposition, one can chose to express the perturbation $\chi_a$ as a sum of Hermite polynomials:

$$ \chi_a = \sum_{m=1}^{\infty} A_{r_1 r_2 \dots r_m}^{(m)} \tilde{H}_{r_1 r_2 \dots r_m}^{(m)} \\ = A_{r_1}^{(1)} \tilde{H}_{r_1}^{(1)} + A_{r_1 r_2}^{(2)} \tilde{H}_{r_1 r_2}^{(2)} + A_{r_1 r_2 r_3}^{(3)} \tilde{H}_{r_1 r_2 r_3}^{(3)} + A_{r_1 r_2 r_3 r_4}^{(4)} \tilde{H}_{r_1 r_2 r_3 r_4}^{(4)} + \dots, \qquad (B13) $$

where the coefficients $A_{r_1 r_2 \dots r_m}^{(m)}$ need to be found. Note that full contractions over all indices are present and the result is a scalar. Multiplying (B13) by weight $\phi^{(0)}$ and polynomial $\tilde{H}_{s_1 s_2 \dots s_n}^{(n)}$, and integrating over $d^3 \tilde{c}$ by using orthogonality (B11), then yields

$$ \int \chi_a \phi^{(0)} \tilde{H}_s^{(n)} d^3 \tilde{c} = A_r^{(n)} \delta_{rs}^{(n)} = n! A_s^{(n)}, \qquad (B14) $$

where the last equality holds because coefficient $A_s^{(n)}$ is a fluid variable and symmetric in all of its indices. The coefficients $A_s^{(n)}$ are thus found according to

$$ A_s^{(n)} = \frac{1}{n!} \int \chi_a \phi^{(0)} \tilde{H}_s^{(n)}(\tilde{c}) d^3 \tilde{c} = \frac{1}{n!} \int (1 + \chi_a) \phi^{(0)} \tilde{H}_s^{(n)}(\tilde{c}) d^3 \tilde{c} = \frac{1}{n!} \left[ \frac{\int_a \tilde{H}_s^{(n)}(\tilde{c}) d^3 \tilde{c}}{\tilde{h}_s^{(n)}} \right], \qquad (B15) $$

where we have used the orthogonality relation (B12) and changed the integration variable to $d^3 c$ with (B8). The quantities in the brackets of (B15) are called Hermite moments $\tilde{h}_s^{(n)}$. The entire Hermite expansion then can be summarized into two easy steps.

(1) Calculate the Hermite moments:

$$ \tilde{h}_{r_1 r_2 \dots r_m}^{a(m)} = \frac{1}{n_a} \int f_a \tilde{H}_{r_1 r_2 \dots r_m}^{a(m)}(\tilde{c}) d^3 c, \qquad (B16) $$

(2) the final perturbation is:

$$ \chi_a = \sum_{m=1}^{\infty} \frac{1}{m!} \tilde{h}_{r_1 r_2 \dots r_m}^{a(m)} \tilde{H}_{r_1 r_2 \dots r_m}^{a(m)}(\tilde{c}). \qquad (B17) $$

It is useful to omit writing the species indices “a” on both $\tilde{h}$ and $\tilde{H}$, as well as on the fluid moments, so we will keep the species index only for $n_a$, $m_a$, $T_a$, and $p_a$. The final perturbations will be written in a full form.
By using definitions of general fluid moments, one straightforwardly calculates the Hermite moments:

$$
\begin{align}
\tilde{h}_i^{(1)} &= \frac{1}{n_a} \int f_a \tilde{H}_i^{(1)} d^3 c = 0; \\
\tilde{h}_{ij}^{(2)} &= \frac{1}{n_a} \int f_a \tilde{H}_{ij}^{(2)} d^3 c = \frac{1}{p_a} \Pi_{ij}^{(2)}; \\
\tilde{h}_{ijk}^{(3)} &= \frac{1}{n_a} \int f_a \tilde{H}_{ijk}^{(3)} d^3 c = \frac{1}{p_a} \sqrt{\frac{m_a}{T_a}} q_{ijk}, \tag{B18}
\end{align}
$$

together with

$$
\begin{align}
\tilde{h}_{ijkl}^{(4)} &= \frac{1}{n_a} \int f_a \tilde{H}_{ijkl}^{(4)}(\mathbf{c}) d^3 c = \frac{\rho_a}{p_a^2} r_{ijkl} + \delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} \nonumber \\
&\quad -\frac{1}{p_a}(\delta_{ij}p_{kl} + \delta_{jk}p_{li} + \delta_{kl}p_{ij} + \delta_{li}p_{jk} + \delta_{ik}p_{jl} + \delta_{jl}p_{ik}) \nonumber \\
&= \frac{\rho_a}{p_a^2} r_{ijkl} - (\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \nonumber \\
&\quad -\frac{1}{p_a}(\delta_{ij}\Pi_{kl}^{(2)} + \delta_{jk}\Pi_{li}^{(2)} + \delta_{kl}\Pi_{ij}^{(2)} + \delta_{li}\Pi_{jk}^{(2)} + \delta_{ik}\Pi_{jl}^{(2)} + \delta_{jl}\Pi_{ik}^{(2)}), \tag{B19}
\end{align}
$$

and

$$
\begin{equation}
\begin{split}
\tilde{h}_{ijklm}^{(5)} = \frac{1}{n_a} \int f_a \tilde{H}_{ijklm}^{(5)}(\mathbf{c}) d^3 c ={}& \frac{\rho_a^{3/2}}{p_a^{5/2}} X_{ijklm}^{(5)} - \frac{\rho_a^{1/2}}{p_a^{3/2}} (\delta_{ij} q_{klm} + \delta_{jk} q_{lim} + \delta_{kl} q_{ijm} + \delta_{li} q_{jkm} \\
& \qquad + \delta_{ik} q_{jlm} + \delta_{jl} q_{ikm} + \delta_{im} q_{jkl} + \delta_{jm} q_{ikl} + \delta_{km} q_{ijl} + \delta_{lm} q_{ijk}).
\end{split}
\tag{B20}
\end{equation}
$$

B.1. Usual Perturbations of Grad

B.1.1. 20-moment Model

Using the definition of the perturbation (B17) and cutting the hierarchy at

$$
\chi_a = \tilde{h}_i^{(1)} \tilde{H}_i^{(1)} + \frac{1}{2} \tilde{h}_{ij}^{(2)} \tilde{H}_{ij}^{(2)} + \frac{1}{6} \tilde{h}_{ijk}^{(3)} \tilde{H}_{ijk}^{(3)}, \quad (B21)
$$

yields the 20-moment perturbation of Grad:

$$
\text{20-moment:} \quad \chi_a = \frac{m_a}{2p_a T_a} (\bar{\Pi}_a^{(2)} : c_a c_a) + \frac{m_a^2}{6p_a T_a^2} (c_a \cdot \bar{q}_a : c_a c_a) - \frac{m_a}{p_a T_a} (q_a \cdot c_a), \quad (B22)
$$

where one defines vector $\mathbf{q}_a = (1/2)\text{Tr}\bar{q}_a$.

B.1.2. 13-moment Model

To quickly obtain the simplified 13-moment model of Grad, one can use $\bar{q} = (2/5)(q\bar{I})^S + \sigma'$ with $\sigma'$ neglected (the validity of this equation is shown below), and calculating $\mathbf{c} \cdot \bar{\mathbf{q}}: \mathbf{cc} = (6/5)(\mathbf{q} \cdot \mathbf{c})c^2$ yields the 13-moment model:

$$
\text{13-moment:} \quad \chi_a = \frac{m_a}{2p_a T_a} (\bar{\Pi}_a^{(2)} : c_a c_a) - \frac{m_a}{p_a T_a} (q_a \cdot c_a) \left(1 - \frac{m_a}{5T_a} c_a^2\right). \quad (B23)
$$

Rederiving the heat flux contribution of the 13-moment model from scratch can be done by using a contracted Hermite polynomial:

$$
\tilde{H}_i^{(3)} \equiv \delta_{jk} \tilde{H}_{ijk}^{(3)} = \tilde{c}_i (\tilde{c}^2 - 5). \quad (B24)
$$

However, one needs to be careful about the normalization constant, because applying the contractions $\delta_{r_1 r_2}$ and $\delta_{s_1 s_2}$ on (B9) yields

$$
\int \phi^{(0)} \tilde{H}_i^{(3)} \tilde{H}_j^{(3)} d^3 \tilde{c} = 10 \delta_{ij}, \quad (B25)
$$

which can be also verified by direct calculation. (Note that it is important to apply the contractions on (B9) as stated above, and not accidentally as $\delta_{r_1 s_1} \delta_{r_2 s_2}$, which would yield an erroneous coefficient 20, as the contraction must satisfy definition (B24)). Then, one
and calculating the heat flux moment again, this time with the Hermite variables:

$$
\begin{align*}
q_{r_1 r_2 r_3} &= m_a \int c_{r_1} c_{r_2} c_{r_3} f_a^{(0)} (1 + \chi_a) d^3 c \\
&= \frac{p_a}{6} \sqrt{\frac{T_a}{m_a}} \tilde{h}_{s_1 s_2 s_3}^{(3)} \int \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{H}_{s_1 s_2 s_3}^{(3)} \phi^{(0)} d^3 \tilde{c} \\
&= \frac{p_a}{6} \sqrt{\frac{T_a}{m_a}} \tilde{h}_{s_1 s_2 s_3}^{(3)} \int \tilde{H}_{r_1 r_2 r_3}^{(3)} \tilde{H}_{s_1 s_2 s_3}^{(3)} \phi^{(0)} d^3 \tilde{c} \\
&= \frac{p_a}{6} \sqrt{\frac{T_a}{m_a}} \tilde{h}_{s_1 s_2 s_3}^{(3)} \delta_{(r_1 r_2 r_3)}^{(3)} (s_1 s_2 s_3) \\
&= p_a \sqrt{\frac{T_a}{m_a}} \tilde{h}_{r_1 r_2 r_3}^{(3)}. \tag{B34}
\end{align*}
$$

In the derivation, we did not use the complicated right-hand side of (B9), we only used $\tilde{h}_r^{(n)} \delta_{rs}^{(n)} = n! \tilde{h}_s^{(n)}$, and the factor of 3! canceled out as well.

Similarly, using the same perturbation (B33), one can derive the fifth-order fluid moment $\bar{\mathbf{X}}^{(5)}$, by using the Hermite polynomial $\tilde{H}_{ijklm}^{(5)}$, Equation (B5), according to

$$
\begin{align*}
X_{r_1 r_2 r_3 r_4 r_5}^{(5)} &= m_a \int c_{r_1} c_{r_2} c_{r_3} c_{r_4} c_{r_5} f_a^{(0)} (1 + \chi_a) d^3 c \\
&= \frac{p_a}{6} \left(\frac{T_a}{m_a}\right)^{3/2} \tilde{h}_{s_1 s_2 s_3}^{(3)} \int \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{c}_{r_4} \tilde{c}_{r_5} \tilde{H}_{s_1 s_2 s_3}^{(3)} \phi^{(0)} d^3 \tilde{c} \\
&= \frac{p_a}{6} \left(\frac{T_a}{m_a}\right)^{3/2} \tilde{h}_{s_1 s_2 s_3}^{(3)} \int \left[ \delta_{r_1 r_2} \tilde{c}_{r_3} \tilde{c}_{r_4} \tilde{c}_{r_5} + \delta_{r_2 r_3} \tilde{c}_{r_4} \tilde{c}_{r_1} \tilde{c}_{r_5} + \delta_{r_3 r_4} \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_5} \right. \\
&\qquad \left. + \delta_{r_4 r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{c}_{r_5} + \delta_{r_1 r_3} \tilde{c}_{r_2} \tilde{c}_{r_4} \tilde{c}_{r_5} + \delta_{r_2 r_4} \tilde{c}_{r_1} \tilde{c}_{r_3} \tilde{c}_{r_5} + \delta_{r_1 r_5} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{c}_{r_4} \right. \\
&\qquad \left. + \delta_{r_2 r_5} \tilde{c}_{r_1} \tilde{c}_{r_3} \tilde{c}_{r_4} + \delta_{r_3 r_5} \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_4} + \delta_{r_4 r_5} \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \right] \tilde{H}_{s_1 s_2 s_3}^{(3)} \phi^{(0)} d^3 \tilde{c} \\
&= p_a \left(\frac{T_a}{m_a}\right)^{3/2} \left[ \delta_{r_1 r_2} \tilde{h}_{r_3 r_4 r_5}^{(3)} + \delta_{r_2 r_3} \tilde{h}_{r_4 r_1 r_5}^{(3)} + \delta_{r_3 r_4} \tilde{h}_{r_1 r_2 r_5}^{(3)} \right. \\
&\qquad \left. + \delta_{r_4 r_1} \tilde{h}_{r_2 r_3 r_5}^{(3)} + \delta_{r_1 r_3} \tilde{h}_{r_2 r_4 r_5}^{(3)} + \delta_{r_2 r_4} \tilde{h}_{r_1 r_3 r_5}^{(3)} + \delta_{r_1 r_5} \tilde{h}_{r_2 r_3 r_4}^{(3)} \right. \\
&\qquad \left. + \delta_{r_2 r_5} \tilde{h}_{r_1 r_3 r_4}^{(3)} + \delta_{r_3 r_5} \tilde{h}_{r_1 r_2 r_4}^{(3)} + \delta_{r_4 r_5} \tilde{h}_{r_1 r_2 r_3}^{(3)} \right]. \tag{B35}
\end{align*}
$$

Or this can be rewritten with the heat fluxes according to (B18) and using the usual indices:

$$
X_{ijklm}^{(5)} = \frac{p_a}{\rho_a} \left[ \begin{aligned}
& \delta_{ij} q_{klm} + \delta_{jk} q_{lim} + \delta_{kl} q_{ijm} + \delta_{li} q_{jkm} + \delta_{ik} q_{jlm} \\
& + \delta_{jl} q_{ikm} + \delta_{im} q_{jkl} + \delta_{jm} q_{ikl} + \delta_{km} q_{ijl} + \delta_{lm} q_{ijk}
\end{aligned} \right], \tag{B36}
$$

and by using the heat flux decomposition (B29) with $\sigma'$ neglected:

$$
X_{ijklm}^{(5)} = \frac{4}{5} \frac{p_a}{\rho_a} \left[ \begin{aligned}
& q_i (\delta_{jk} \delta_{lm} + \delta_{kl} \delta_{jm} + \delta_{jl} \delta_{km}) + q_j (\delta_{kl} \delta_{im} + \delta_{il} \delta_{km} + \delta_{ik} \delta_{lm}) \\
& + q_k (\delta_{ij} \delta_{lm} + \delta_{il} \delta_{jm} + \delta_{jl} \delta_{im}) + q_l (\delta_{ij} \delta_{km} + \delta_{jk} \delta_{im} + \delta_{ik} \delta_{jm}) \\
& + q_m (\delta_{ij} \delta_{kl} + \delta_{jk} \delta_{il} + \delta_{ik} \delta_{jl})
\end{aligned} \right]. \tag{B37}
$$

Applying the contractions at (B36) yields

$$
\begin{align}
[\mathrm{Tr}\,\bar{\mathbf{X}}^{(5)}]_{ijk} &= \frac{P_a}{\rho_a} [2(\bar{\mathbf{I}}\mathbf{q})^S + 9\bar{\mathbf{q}}]_{ijk} = \frac{28}{5} \frac{P_a}{\rho_a} (\bar{\mathbf{I}}\mathbf{q})_{ijk}^S + 9 \frac{P_a}{\rho_a} \sigma_{ijk}'; \nonumber \\
\mathbf{X}^{(5)} &= \mathrm{Tr}\mathrm{Tr}\,\bar{\mathbf{X}}^{(5)} = 28 \frac{P_a}{\rho_a} \mathbf{q}. \tag{B38}
\end{align}
$$
## B.2. Higher-order Perturbations (Full $\bar{\mathbf{X}}^{(4)}$ and $\bar{\mathbf{X}}^{(5)}$ Moments)

By using the technique described above, it is possible to use the following higher-order perturbation:

$$ \chi_a = \frac{1}{2!} \tilde{h}_{s_1 s_2}^{(2)} \tilde{H}_{s_1 s_2}^{(2)} + \frac{1}{3!} \tilde{h}_{s_1 s_2 s_3}^{(3)} \tilde{H}_{s_1 s_2 s_3}^{(3)} + \frac{1}{4!} \tilde{h}_{s_1 s_2 s_3 s_4}^{(4)} \tilde{H}_{s_1 s_2 s_3 s_4}^{(4)} + \frac{1}{5!} \tilde{h}_{s_1 s_2 s_3 s_4 s_5}^{(5)} \tilde{H}_{s_1 s_2 s_3 s_4 s_5}, \quad (B39) $$

and directly calculate the fluid moments (we use $\bar{\mathbf{X}}^{(4)}$ instead of $\bar{\mathbf{r}}$ from now on):

$$ [\bar{\mathbf{X}}^{(4)}]_{r_1 r_2 r_3 r_4} = \frac{p_a^2}{\rho_a} [\tilde{h}_{r_1 r_2 r_3 r_4}^{(4)} + \delta_{r_1 r_2} \tilde{h}_{r_3 r_4}^{(2)} + \delta_{r_2 r_3} \tilde{h}_{r_1 r_4}^{(2)} + \delta_{r_3 r_4} \tilde{h}_{r_1 r_2}^{(2)} \\ + \delta_{r_1 r_4} \tilde{h}_{r_2 r_3}^{(2)} + \delta_{r_1 r_3} \tilde{h}_{r_2 r_4}^{(2)} + \delta_{r_2 r_4} \tilde{h}_{r_1 r_3}^{(2)} \\ + \delta_{r_1 r_2} \delta_{r_3 r_4} + \delta_{r_1 r_3} \delta_{r_2 r_4} + \delta_{r_2 r_3} \delta_{r_1 r_4}], \quad (B40) $$

and

$$ [\bar{\mathbf{X}}^{(5)}]_{r_1 r_2 r_3 r_4 r_5} = p_a \left(\frac{T_a}{m_a}\right)^{3/2} [\tilde{h}_{r_1 r_2 r_3 r_4 r_5}^{(5)} + \delta_{r_1 r_2} \tilde{h}_{r_3 r_4 r_5}^{(3)} + \delta_{r_2 r_3} \tilde{h}_{r_4 r_1 r_5}^{(3)} + \delta_{r_3 r_4} \tilde{h}_{r_1 r_2 r_5}^{(3)} \\ + \delta_{r_3 r_1} \tilde{h}_{r_2 r_3 r_5}^{(3)} + \delta_{r_1 r_3} \tilde{h}_{r_2 r_4 r_5}^{(3)} + \delta_{r_2 r_4} \tilde{h}_{r_1 r_3 r_5}^{(3)} + \delta_{r_1 r_5} \tilde{h}_{r_2 r_3 r_4}^{(3)} \\ + \delta_{r_2 r_5} \tilde{h}_{r_1 r_3 r_4}^{(3)} + \delta_{r_3 r_5} \tilde{h}_{r_1 r_2 r_4}^{(3)} + \delta_{r_4 r_5} \tilde{h}_{r_1 r_2 r_3}^{(3)}]. \quad (B41) $$

Both results contain new contributions, represented by the $\tilde{h}_{r_1 r_2 r_3 r_4}^{(4)}$ and $\tilde{h}_{r_1 r_2 r_3 r_4 r_5}^{(5)}$.

It is useful to introduce notation where, by applying a contraction at a tensor, the contracted indices will be suppressed, so, for example, $\tilde{h}_i^{(3)} \equiv \tilde{h}_{ijk}^{(3)}$, or $X_{ij}^{(4)} \equiv X_{ijk}^{(4)}$ and $X^{(4)} \equiv X_{ijk}^{(4)}$. We define all of the contractions without any additional factors, with the sole exception of the heat flux vector $\mathbf{q}$, where the additional factor of 1/2 is present, to match its usual definition. To emphasize this difference, in the index notation we thus keep an arrow on the components of the heat flux vector $\mathbf{q}_i$, to clearly distinguish it from the contracted tensor $\mathbf{q}_{ijk}$.

Applying the contractions at (B40), (B41) then yields

$$ \begin{aligned} X_{ij}^{(4)} &= \frac{p_a^2}{\rho_a} [\tilde{h}_{ij}^{(4)} + 7\tilde{h}_{ij}^{(2)} + 5\delta_{ij}]; \\ X_{ij}^{(4)} &= \frac{p_a^2}{\rho_a} [\tilde{h}^{(4)} + 15]; \\ X_{ijk}^{(5)} &= \frac{p_a^2}{\rho_a} \sqrt{\frac{T_a}{m_a}} [\tilde{h}_{ijk}^{(5)} + \delta_{ij} \tilde{h}_k^{(3)} + \delta_{jk} \tilde{h}_i^{(3)} + \delta_{ik} \tilde{h}_j^{(3)} + 9\tilde{h}_{ijk}^{(3)}]; \\ X_i^{(5)} &= \frac{p_a^2}{\rho_a} \sqrt{\frac{T_a}{m_a}} [\tilde{h}_i^{(5)} + 14\tilde{h}_i^{(3)}], \end{aligned} \quad (B42) $$

and applying the contractions at the Hermite moments (B18)–(B20) yields

$$ \begin{aligned} \tilde{h}_{ij}^{(2)} &= \frac{1}{p_a} \Pi_{ij}^{(2)}; & \tilde{h}_i^{(3)} &= \frac{2}{p_a} \sqrt{\frac{m_a}{T_a}} \vec{q}_i; \\ \tilde{h}_{ij}^{(4)} &= \frac{\rho_a}{p_a^2} X_{ij}^{(4)} - 5\delta_{ij} - \frac{7}{p_a} \Pi_{ij}^{(2)}; \\ \tilde{h}^{(4)} &= \frac{\rho_a}{p_a^2} X^{(4)} - 15; \\ \tilde{h}_{ijk}^{(5)} &= \frac{1}{p_a} \sqrt{\frac{m_a}{T_a}} \left[ \frac{\rho_a}{p_a} X_{ijk}^{(5)} - (2\delta_{ij}\vec{q}_k + 2\delta_{jk}\vec{q}_i + 2\delta_{ik}\vec{q}_j + 9q_{ijk}) \right]; \\ \tilde{h}_i^{(5)} &= \frac{1}{p_a} \sqrt{\frac{m_a}{T_a}} \left( \frac{\rho_a}{p_a} X_i^{(5)} - 28\vec{q}_i \right). \end{aligned} \quad (B43) $$
### B.2.1. Viscosity $\Pi_{ij}^{(4)}$ of the Fourth-order Moment $X_{ij}^{(4)}$

The usual viscosity tensor is defined as a traceless matrix:

$$ \Pi_{ij}^{(2)} = m_a \int \left( c_i c_j - \frac{1}{3} \delta_{ij} c^2 \right) f_a d^3 c. \qquad (B44) $$

Similarly, it is beneficial to introduce a traceless viscosity tensor of the fourth-order fluid moment:

$$ \Pi_{ij}^{(4)} = m_a \int \left( c_i c_j - \frac{1}{3} \delta_{ij} c^2 \right) c^2 f_a d^3 c. \qquad (B45) $$

In other words, the moment $X_{ij}^{(4)}$ is decomposed as

$$ X_{ij}^{(4)} = \frac{\delta_{ij}}{3} X^{(4)} + \Pi_{ij}^{(4)}, \qquad (B46) $$

where the fully contracted $X^{(4)} = m_a \int c^4 f_a d^3 c$. Scalar $X^{(4)}$ is further decomposed to its “core” Maxwellian part, and the additional perturbation $\tilde{X}^{(4)}$ (with wide tilde), according to

$$ X^{(4)} = 15 \frac{P_a^2}{\rho_a} + \tilde{X}^{(4)}, \qquad (B47) $$

and the corresponding Hermite moments thus become

$$ \begin{aligned} \tilde{h}_{ij}^{(4)} &= \frac{\rho_a}{P_a^2} \frac{\delta_{ij}}{3} \tilde{X}^{(4)} + \frac{\rho_a}{P_a^2} \Pi_{ij}^{(4)} - \frac{7}{P_a} \Pi_{ij}^{(2)}; \\ \tilde{h}^{(4)} &= \frac{\rho_a}{P_a^2} \tilde{X}^{(4)}. \end{aligned} \qquad (B48) $$

It is important to emphasize that depending on the choice of perturbation $\chi_a$, in general $\tilde{X}^{(4)}$ is nonzero. However, this perturbation is not required to derive the model of Braginskii (1965), and Balescu (1988), for example, prescribes irreducible $h^{(4)} = 0$. In the next section, we will consider simplified perturbations and derive the above results in a more direct manner; nevertheless, the more general case (B39) is a very useful guide in demonstrating that it is possible to consider perturbations with nonzero $h^{(4)}$.

Finally, because the reducible matrix $\tilde{h}_{ij}^{(4)}$ is not traceless in general (unless one prescribes the Hermite closure $\tilde{h}^{(4)} = 0$, which makes it traceless by definition), it is useful to introduce traceless

$$ \tilde{h}_{ij}^{(4)} = \tilde{h}_{ij}^{(4)} - \frac{\delta_{ij}}{3} \tilde{h}^{(4)} = \frac{\rho_a}{P_a^2} \Pi_{ij}^{(4)} - \frac{7}{P_a} \Pi_{ij}^{(2)}, \qquad (B49) $$

where we use hat instead of tilde.

### B.2.2. Simplified Perturbations (21-moment Model)

Instead of working with very complicated perturbations (B39), it was shown by Balescu (1988) that to obtain the model of Braginskii (1965), it is sufficient to work with simplified

$$ \chi_a = h_{ij}^{(2)} H_{ij}^{(2)} + h_i^{(3)} H_i^{(3)} + h_{ij}^{(4)} H_{ij}^{(4)} + h_i^{(5)} H_i^{(5)}. \qquad (B50) $$

The perturbation (B50) is written with irreducible Hermite polynomials (notation without tilde), as discussed in the next section. This perturbation represents the 21-moment model, and recovers both the stress tensor and the heat flux of Braginskii. However, the connection between irreducible and reducible Hermite polynomials can be very blurry at first, and we continue with reducible Hermite polynomials.

Applying the contractions at the hierarchy of the reducible polynomials (B5) yields

$$ \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); \\ \tilde{H}_{ij}^{(2)} &= \tilde{c}_i \tilde{c}_j - \delta_{ij}; & \tilde{H}_{ij}^{(4)} &= \tilde{c}_i \tilde{c}_j (\tilde{c}^2 - 7) - \delta_{ij} (\tilde{c}^2 - 5). \end{aligned} \qquad (B51) $$

By using these polynomials, the Hermite moments then calculate

$$ \begin{aligned} \tilde{h}_i^{(3)} &= \frac{2}{P_a} \sqrt{\frac{m_a}{T_a}} \vec{q}_i; & \tilde{h}_i^{(5)} &= \frac{1}{P_a} \sqrt{\frac{m_a}{T_a}} \left( \frac{\rho_a}{P_a} X_i^{(5)} - 28 \vec{q}_i \right); \\ \tilde{h}_{ij}^{(2)} &= \frac{1}{P_a} \Pi_{ij}^{(2)}; & \hat{h}_{ij}^{(4)} &= \frac{\rho_a}{P_a^2} \Pi_{ij}^{(4)} - \frac{7}{P_a} \Pi_{ij}^{(2)}, \end{aligned} $$
of course recovering previous results. The reducible Hermite polynomials satisfy the following orthogonality relations:

$$
\begin{align}
\int \tilde{H}_i^{(3)} \tilde{H}_j^{(3)} \phi^{(0)} d^3 \tilde{c} &= 10 \delta_{ij}; & \int \tilde{H}_i^{(5)} \tilde{H}_j^{(5)} \phi^{(0)} d^3 \tilde{c} &= 280 \delta_{ij}; \nonumber \\
\int \tilde{H}_{ij}^{(2)} \tilde{H}_{kl}^{(2)} \phi^{(0)} d^3 \tilde{c} &= \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}; \nonumber \\
\int \tilde{H}_{ij}^{(4)} \tilde{H}_{kl}^{(4)} \phi^{(0)} d^3 \tilde{c} &= 14(\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}) + 4 \delta_{ij} \delta_{kl}, \tag{B52}
\end{align}
$$

and because the Hermite moments $\tilde{h}_{kl}^{(2)}$, $\hat{h}_{kl}^{(4)}$ are symmetric and traceless,

$$
\tilde{h}_{kl}^{(2)} \int \tilde{H}_{ij}^{(2)} \tilde{H}_{kl}^{(2)} \phi^{(0)} d^3 \tilde{c} = 2 \tilde{h}_{ij}^{(2)}; \quad \hat{h}_{kl}^{(4)} \int \tilde{H}_{ij}^{(4)} \tilde{H}_{kl}^{(4)} \phi^{(0)} d^3 \tilde{c} = 28 \hat{h}_{ij}^{(4)}. \qquad (B53)
$$

Thus, a perturbation that can be directly derived from the hierarchy of the reducible Hermite polynomials (with no reference to irreducible Hermite polynomials or Laguerre–Sonine polynomials) reads

$$
\chi_a = \frac{1}{2} \tilde{h}_{ij}^{(2)} \tilde{H}_{ij}^{(2)} + \frac{1}{10} \tilde{h}_i^{(3)} \tilde{H}_i^{(3)} + \frac{1}{28} \hat{h}_{ij}^{(4)} \tilde{H}_{ij}^{(4)} + \frac{1}{280} \tilde{h}_i^{(5)} \tilde{H}_i^{(5)}, \quad (B54)
$$

where each term is calculated as

$$
\begin{align*}
\frac{1}{10} \tilde{h}_i^{(3)} \tilde{H}_i^{(3)} &= \frac{1}{5p_a} \sqrt{\frac{m_a}{T_a}} (\mathbf{q}_a \cdot \tilde{\mathbf{c}}_a)(\tilde{c}_a^2 - 5); \\
\frac{1}{280} \tilde{h}_i^{(5)} \tilde{H}_i^{(5)} &= \frac{1}{280p_a} \sqrt{\frac{m_a}{T_a}} \left[ \frac{\rho_a}{p_a} (\mathbf{X}^{a(5)} \cdot \tilde{\mathbf{c}}_a) - 28(\mathbf{q}_a \cdot \tilde{\mathbf{c}}_a) \right] (\tilde{c}_a^4 - 14\tilde{c}_a^2 + 35); \\
\frac{1}{2} \tilde{h}_{ij}^{(2)} \tilde{H}_{ij}^{(2)} &= \frac{1}{2p_a} (\tilde{\mathbf{\Pi}}_a^{(2)} : \tilde{\mathbf{c}}_a \tilde{\mathbf{c}}_a); \\
\frac{1}{28} \hat{h}_{ij}^{(4)} \tilde{H}_{ij}^{(4)} &= \frac{1}{28} \left[ \frac{\rho_a}{p_a^2} (\tilde{\mathbf{\Pi}}_a^{(4)} : \tilde{\mathbf{c}}_a \tilde{\mathbf{c}}_a) - \frac{7}{p_a} (\tilde{\mathbf{\Pi}}_a^{(2)} : \tilde{\mathbf{c}}_a \tilde{\mathbf{c}}_a) \right] (\tilde{c}_a^2 - 7), \tag{B55}
\end{align*}
$$

with normalized velocity $\tilde{\mathbf{c}}_a = \sqrt{m_a/T_a} \mathbf{c}_a$. Below, we show that the perturbation (B54), (B55) is equivalent to the perturbation of Balescu (B50), obtained with irreducible polynomials. The final perturbation of the 21-moment model that recovers Braginskii (1965) thus reads

$$
\chi_a = \frac{1}{2p_a} (\bar{\Pi}_a^{(2)} : \tilde{\mathbf{c}}_a \tilde{\mathbf{c}}_a) + \frac{1}{28} \left[ \frac{\rho_a}{p_a^2} (\bar{\Pi}_a^{(4)} : \tilde{\mathbf{c}}_a \tilde{\mathbf{c}}_a) - \frac{7}{p_a} (\bar{\Pi}_a^{(2)} : \tilde{\mathbf{c}}_a \tilde{\mathbf{c}}_a) \right] (\tilde{c}_a^2 - 7) \\
+ \frac{1}{5p_a} \sqrt{\frac{m_a}{T_a}} (\mathbf{q}_a \cdot \tilde{\mathbf{c}}_a)(\tilde{c}_a^2 - 5) + \frac{1}{280p_a} \sqrt{\frac{m_a}{T_a}} \left[ \frac{\rho_a}{p_a} (\mathbf{X}^{a(5)} \cdot \tilde{\mathbf{c}}_a) - 28(\mathbf{q}_a \cdot \tilde{\mathbf{c}}_a) \right] (\tilde{c}_a^4 - 14\tilde{c}_a^2 + 35). \quad (B56)
$$

Finally, because $\hat{h}_{ij}^{(4)}$ is traceless, its double contraction with $\tilde{H}_{ij}^{(4)}$ makes the part of this polynomial proportional to $\delta_{ij}$ redundant in the final perturbation. It is possible to define another traceless polynomial (with hat instead of tilde):

$$
\hat{H}_{ij}^{(4)} = \tilde{H}_{ij}^{(4)} - \frac{\delta_{ij}}{3} \tilde{H}^{(4)} = \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) (\tilde{c}^2 - 7), \quad (B57)
$$

and replace the following term in the perturbation (B54):

$$
\hat{h}_{ij}^{(4)} \tilde{H}_{ij}^{(4)} = \hat{h}_{ij}^{(4)} \hat{H}_{ij}^{(4)}, \qquad (B58)
$$

where the part of (B57) proportional to $\delta_{ij}$ is still suppressed in the final perturbation. However, the traceless definition (B57) now makes it possible to directly define the traceless Hermite moment $\hat{h}_{ij}^{(4)}$ as an integral over $\hat{H}_{ij}^{(4)}$:

$$
\hat{h}_{ij}^{(4)} = \frac{1}{n_a} \int f_a \hat{H}_{ij}^{(4)} d^3 c. \qquad (B59)
$$

This is the main motivation behind irreducible Hermite polynomials, as is further clarified below.
### B.3. Irreducible Hermite Polynomials

In the work of Balescu (1988), the irreducible Hermite polynomials are defined through Laguerre–Sonine polynomials, according to (see Equation (G1.4.4) on p. 326 of Balescu)

$$ 
\begin{aligned}
H^{(2n)}(\tilde{c}) &= L_n^{1/2}\left(\frac{\tilde{c}^2}{2}\right); \\
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).
\end{aligned}
\quad (B60) $$

To recover the Braginskii (1965) model, one only needs (see Table 4.1 on p. 327 of Balescu)

$$ 
\begin{aligned}
H_i^{(3)} &= \frac{1}{\sqrt{10}}\tilde{c}_i(\tilde{c}^2 - 5); & H_i^{(5)} &= \frac{1}{2\sqrt{70}}\tilde{c}_i(\tilde{c}^4 - 14\tilde{c}^2 + 35); \\
H_{ij}^{(2)} &= \frac{1}{\sqrt{2}}(\tilde{c}_i\tilde{c}_j - \frac{1}{3}\tilde{c}^2\delta_{ij}); & H_{ij}^{(4)} &= \frac{1}{2\sqrt{7}}(\tilde{c}_i\tilde{c}_j - \frac{1}{3}\tilde{c}^2\delta_{ij})(\tilde{c}^2 - 7),
\end{aligned}
\quad (B61) $$

yielding Hermite moments

$$ 
\begin{aligned}
h_i^{(3)} &= \sqrt{\frac{2}{5}}\frac{1}{p_a}\sqrt{\frac{m_a}{T_a}}\vec{q}_i; & h_i^{(5)} &= \frac{1}{2\sqrt{70}}\frac{p_a}{T_a}\sqrt{\frac{m_a}{T_a}}\left[\frac{\rho_a}{p_a}X_i^{(5)} - 28\vec{q}_i\right]; \\
h_{ij}^{(2)} &= \frac{1}{\sqrt{2}}\Pi_{ij}^{(2)}; & h_{ij}^{(4)} &= \frac{1}{2\sqrt{7}}\Pi_{ij}^{(4)} - 7\Pi_{ij}^{(2)}.
\end{aligned}
\quad (B62) $$

Furthermore, the orthogonal relations are

$$ \int \phi^{(0)} H_i^{(2n+1)} H_j^{(2n+1)} d^3c = \delta_{ij}; \quad h_{kl}^{(2n)} \int \phi^{(0)} H_{ij}^{(2n)} H_{kl}^{(2n)} d^3c = h_{ij}^{(2n)}, \quad (B63) $$

yielding perturbation (B50), which then recovers perturbation (B54), (B55) obtained with reducible polynomials. Both approaches are therefore equivalent, which is further addressed in Appendix B.7.

#### B.3.1. Higher-order Tensorial “Anisotropies”

It is useful to clarify what contributions are obtained by using the irreducible Hermite polynomials:

$$ 
\begin{aligned}
H_{ijk}^{(3)}(\tilde{c}) &= \tilde{c}_i \tilde{c}_j \tilde{c}_k - \frac{1}{5} \tilde{c}^2 (\delta_{ij} \tilde{c}_k + \delta_{jk} \tilde{c}_i + \delta_{ik} \tilde{c}_j); \\
H_{ijkl}^{(4)}(\tilde{c}) &= \tilde{c}_i \tilde{c}_j \tilde{c}_k \tilde{c}_l - \frac{1}{7} \tilde{c}^2 (\delta_{ij} \tilde{c}_k \tilde{c}_l + \delta_{jk} \tilde{c}_i \tilde{c}_k + \delta_{kl} \tilde{c}_i \tilde{c}_j + \delta_{li} \tilde{c}_j \tilde{c}_k + \delta_{ik} \tilde{c}_j \tilde{c}_l + \delta_{jl} \tilde{c}_i \tilde{c}_k) \\
&\quad + \frac{1}{35} \tilde{c}^4 (\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}),
\end{aligned}
\quad (B64) $$

which Balescu (1988) calls “anisotropies” (even though they are valid as a perturbation for isotropic Maxwellians). Importantly, applying a trace on (B64) yields zero. The corresponding Hermite moments calculate

$$ h_{ijk}^{(3)} = \frac{1}{p_a} \sqrt{\frac{m_a}{T_a}} \left[ \bar{q} - \frac{2}{5} (\bar{I} \bar{q})^S \right]_{ijk} = \frac{1}{p_a} \sqrt{\frac{m_a}{T_a}} \sigma'_{ijk}; \quad h_{ijkl}^{(4)} = \frac{\rho_a}{p_a^2} \sigma_{ijkl}^{(4)\prime}, \quad (B65) $$

and directly yield the highest-order irreducible parts.

### B.4. Decomposition of $X_{ijkl}^{(4)}$

We continue with the reducible Hermite polynomials. To decompose the full fourth-order fluid moment $X_{ijkl}^{(4)}$, it is necessary to consider the following perturbation (i.e., the 16-moment model):

$$ \chi_a = \frac{1}{2} \tilde{h}_{ij}^{(2)} \tilde{H}_{ij}^{(2)} + \frac{1}{28} \hat{h}_{ij}^{(4)} \tilde{H}_{ij}^{(4)} + \frac{1}{120} \tilde{h}^{(4)} \tilde{H}^{(4)}, \quad (B66) $$

and, by using this perturbation, to calculate $X_{ijkl}^{(4)}$. In comparison to the previous perturbation of the 21-moment model, the last term with Hermite polynomial $\tilde{H}^{(4)}$ is new. It is derived with orthogonality relation $\int \phi^{(0)} H^{(4)} H^{(4)} d^3c = 120$. We will need the following
integrals. Applying the contraction $\delta_{r_3 r_4}$ at the orthogonality relation (B10) yields

$$ \begin{aligned} \int \tilde{H}_{r_1 r_2}^{(4)} \tilde{c}_{s_1} \tilde{c}_{s_2} \tilde{c}_{s_3} \tilde{c}_{s_4} \phi^{(0)} d^3 \tilde{c} &= \int \tilde{H}_{r_1 r_2}^{(4)} \tilde{H}_{s_1 s_2 s_3 s_4}^{(4)} \phi^{(0)} d^3 \tilde{c} \\ &= + 2\delta_{r_1 s_1} \delta_{r_2 s_2} \delta_{s_3 s_4} + 2\delta_{r_1 s_1} \delta_{r_2 s_3} \delta_{s_2 s_4} + 2\delta_{r_1 s_1} \delta_{r_2 s_4} \delta_{s_2 s_3} \\ &\quad + 2\delta_{r_1 s_2} \delta_{r_2 s_1} \delta_{s_3 s_4} + 2\delta_{r_1 s_2} \delta_{r_2 s_3} \delta_{s_1 s_4} + 2\delta_{r_1 s_2} \delta_{r_2 s_4} \delta_{s_1 s_3} \\ &\quad + 2\delta_{r_1 s_3} \delta_{r_2 s_1} \delta_{s_2 s_4} + 2\delta_{r_1 s_3} \delta_{r_2 s_2} \delta_{s_1 s_4} + 2\delta_{r_1 s_3} \delta_{r_2 s_4} \delta_{s_1 s_2} \\ &\quad + 2\delta_{r_1 s_4} \delta_{r_2 s_1} \delta_{s_2 s_3} + 2\delta_{r_1 s_4} \delta_{r_2 s_2} \delta_{s_1 s_3} + 2\delta_{r_1 s_4} \delta_{r_2 s_3} \delta_{s_1 s_2}, \end{aligned} \quad (B67) $$

and, further, applying traceless $\hat{h}^{(4)}$ at (B67) leads to

$$ \begin{aligned} \hat{h}_{r_1 r_2}^{(4)} \int \tilde{H}_{r_1 r_2}^{(4)} \tilde{c}_{s_1} \tilde{c}_{s_2} \tilde{c}_{s_3} \tilde{c}_{s_4} \phi^{(0)} d^3 \tilde{c} &= 4[\hat{h}_{s_1 s_2}^{(4)} \delta_{s_3 s_4} + \hat{h}_{s_1 s_3}^{(4)} \delta_{s_2 s_4} + \hat{h}_{s_1 s_4}^{(4)} \delta_{s_2 s_3} + \hat{h}_{s_2 s_3}^{(4)} \delta_{s_1 s_4} + \hat{h}_{s_2 s_4}^{(4)} \delta_{s_1 s_3} + \hat{h}_{s_3 s_4}^{(4)} \delta_{s_1 s_2}]. \end{aligned} \quad (B68) $$

Applying the contraction $\delta_{r_1 r_2}$ at (B67) and multiplying by $\tilde{h}^{(4)}$ yields

$$ \tilde{h}^{(4)} \int \tilde{H}^{(4)} \tilde{c}_{s_1} \tilde{c}_{s_2} \tilde{c}_{s_3} \tilde{c}_{s_4} \phi^{(0)} d^3 \tilde{c} = 8 \tilde{h}^{(4)} [\delta_{s_1 s_2} \delta_{s_3 s_4} + \delta_{s_1 s_3} \delta_{s_2 s_4} + \delta_{s_1 s_4} \delta_{s_2 s_3}]. \quad (B69) $$

Similarly,

$$ \begin{aligned} \tilde{h}_{r_1 r_2}^{(2)} \int \tilde{H}_{r_1 r_2}^{(2)} \tilde{c}_{s_1} \tilde{c}_{s_2} \tilde{c}_{s_3} \tilde{c}_{s_4} \phi^{(0)} d^3 \tilde{c} &= 2[\tilde{h}_{s_1 s_2}^{(2)} \delta_{s_3 s_4} + \tilde{h}_{s_1 s_3}^{(2)} \delta_{s_2 s_4} + \tilde{h}_{s_1 s_4}^{(2)} \delta_{s_2 s_3} + \tilde{h}_{s_2 s_3}^{(2)} \delta_{s_1 s_4} + \tilde{h}_{s_2 s_4}^{(2)} \delta_{s_1 s_3} + \tilde{h}_{s_3 s_4}^{(2)} \delta_{s_1 s_2}]. \end{aligned} \quad (B70) $$

The results (B68), (B69), (B70) allow us to calculate the $X_{ijkl}^{(4)}$ moment, which becomes

$$ \begin{aligned} X_{ijkl}^{(4)} &= m_a \int f_a^{(0)} \left[ 1 + \frac{1}{2} \tilde{h}_{r_1 r_2}^{(2)} \tilde{H}_{r_1 r_2}^{(2)} + \frac{1}{28} \hat{h}_{r_1 r_2}^{(4)} \tilde{H}_{r_1 r_2}^{(4)} + \frac{1}{120} \tilde{h}^{(4)} \tilde{H}^{(4)} \right] c_i c_j c_k c_l d^3 c \\ &= + \frac{1}{15} \frac{p_a^2}{\rho_a} (15 + \tilde{h}^{(4)})(\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{jk} \delta_{il}) \\ &\quad + \frac{1}{7} \frac{p_a^2}{\rho_a} [(\hat{h}_{ij}^{(4)} + 7 \tilde{h}_{ij}^{(2)}) \delta_{kl} + (\hat{h}_{ik}^{(4)} + 7 \tilde{h}_{ik}^{(2)}) \delta_{jl} + (\hat{h}_{il}^{(4)} + 7 \tilde{h}_{il}^{(2)}) \delta_{jk} \\ &\quad + (\hat{h}_{jk}^{(4)} + 7 \tilde{h}_{jk}^{(2)}) \delta_{il} + (\hat{h}_{jl}^{(4)} + 7 \tilde{h}_{jl}^{(2)}) \delta_{ik} + (\hat{h}_{kl}^{(4)} + 7 \tilde{h}_{kl}^{(2)}) \delta_{ij}]. \end{aligned} \quad (B71) $$

Form (B71) nicely shows how various parts of perturbation (B66) contribute to the decomposition, including the new $\tilde{h}^{(4)}$. Prescribing the Hermite closures $\hat{h}_{ij}^{(4)} = 0$, $\tilde{h}^{(4)} = 0$ recovers decomposition (B30) used in the Burgers–Schunk model. Finally, rewritten with the fluid moments

$$ \Pi_{ij}^{(4)} = \frac{p_a^2}{\rho_a} (\hat{h}_{ij}^{(4)} + 7 \tilde{h}_{ij}^{(2)}); \quad \tilde{X}^{(4)} = \frac{p_a^2}{\rho_a} \tilde{h}^{(4)}; \quad X^{(4)} = 15 \frac{p_a^2}{\rho_a} + \tilde{X}^{(4)}, \quad (B72) $$

and representing all of the other terms that were not obtained from (B66) by traceless $\sigma_{ijkl}^{(4)'}$ (which represents the highest-order irreducible part of $X_{ijkl}^{(4)}$), the decomposition becomes

$$ X_{ijkl}^{(4)} = \frac{1}{15} X^{(4)} (\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}) \\ + \frac{1}{7} [\Pi_{ij}^{(4)} \delta_{kl} + \Pi_{ik}^{(4)} \delta_{jl} + \Pi_{il}^{(4)} \delta_{jk} + \Pi_{jk}^{(4)} \delta_{il} + \Pi_{jl}^{(4)} \delta_{ik} + \Pi_{kl}^{(4)} \delta_{ij}] + \sigma_{ijkl}^{(4)'} \quad (B73) $$

or equivalently,

$$ X_{ijkl}^{(4)} = -\frac{1}{35} X^{(4)} (\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}) \\ + \frac{1}{7} [X_{ij}^{(4)} \delta_{kl} + X_{ik}^{(4)} \delta_{jl} + X_{il}^{(4)} \delta_{jk} + X_{jk}^{(4)} \delta_{il} + X_{jl}^{(4)} \delta_{ik} + X_{kl}^{(4)} \delta_{ij}] + \sigma_{ijkl}^{(4)'} \quad (B74) $$

Decomposition (B73) is equivalent to Equation (30.22) of Grad (1958). Essentially, any tensorial moment can be decomposed by subtracting all of the possible contractions of that moment. Note that simply prescribing the closure $\Pi_{ij}^{(4)} = 0$ in (B73) would be erroneous, unless one also prescribes $\Pi_{ij}^{(2)} = 0$ as well. The correct simplification of (B73) is obtained by prescribing the Hermite closure $\hat{h}_{ij}^{(4)} = 0$, meaning by prescribing the fluid closure $\Pi_{ij}^{(4)} = 7(p_a/\rho_a)\Pi_{ij}^{(2)}$. Additionally, one can also prescribe the Hermite closure $\hat{h}_{ij}^{(4)} = 0$, meaning by prescribing the fluid closure $\Pi_{ij}^{(4)} = 7(p_a/\rho_a)\Pi_{ij}^{(2)}$.
closure $\tilde{h}^{(4)} = 0$, which is equivalent to the fluid closure $\tilde{X}^{(4)} = 0$.

## B.5. Decomposition of $X_{ijklm}^{(5)}$

We only use the simplified perturbation

$$ \chi_a = \frac{1}{10} \tilde{h}_i^{(3)} \tilde{H}_i^{(3)} + \frac{1}{280} \tilde{h}_i^{(5)} \tilde{H}_i^{(5)}. \qquad (B75) $$

By using this perturbation, it is possible to calculate the fifth-order fluid moment:

$$ \begin{aligned} X_{s_1 s_2 s_3 s_4 s_5}^{(5)} = & \frac{1}{35} \frac{p_a^{5/2}}{\rho_a^{3/2}} [ (\tilde{h}_{s_1}^{(5)} + 14\tilde{h}_{s_1}^{(3)})(\delta_{s_2 s_3} \delta_{s_4 s_5} + \delta_{s_2 s_4} \delta_{s_3 s_5} + \delta_{s_2 s_5} \delta_{s_3 s_4}) \\ & + (\tilde{h}_{s_2}^{(5)} + 14\tilde{h}_{s_2}^{(3)})(\delta_{s_1 s_3} \delta_{s_4 s_5} + \delta_{s_1 s_4} \delta_{s_3 s_5} + \delta_{s_1 s_5} \delta_{s_3 s_4}) \\ & + (\tilde{h}_{s_3}^{(5)} + 14\tilde{h}_{s_3}^{(3)})(\delta_{s_1 s_2} \delta_{s_4 s_5} + \delta_{s_1 s_4} \delta_{s_2 s_5} + \delta_{s_1 s_5} \delta_{s_2 s_4}) \\ & + (\tilde{h}_{s_4}^{(5)} + 14\tilde{h}_{s_4}^{(3)})(\delta_{s_1 s_2} \delta_{s_3 s_5} + \delta_{s_1 s_3} \delta_{s_2 s_5} + \delta_{s_1 s_5} \delta_{s_2 s_3}) \\ & + (\tilde{h}_{s_5}^{(5)} + 14\tilde{h}_{s_5}^{(3)})(\delta_{s_1 s_2} \delta_{s_3 s_4} + \delta_{s_1 s_3} \delta_{s_2 s_4} + \delta_{s_1 s_4} \delta_{s_2 s_3})]. \end{aligned} \qquad (B76) $$

Because we consider the simplified perturbation (B75), we do not consider full decomposition with $\sigma^{(5)}$. Prescribing the Hermite closure $\tilde{h}_i^{(5)} = 0$ yields the previously obtained decomposition (B37). Finally, by switching from Hermite to fluid moments,

$$ \tilde{h}_i^{(3)} = 2 \frac{\rho_a^{1/2}}{p_a^{3/2}} \tilde{q}_i; \qquad \tilde{h}_i^{(5)} = \frac{\rho_a^{1/2}}{p_a^{3/2}} \left( \frac{\rho_a}{p_a} X_i^{(5)} - 28 \tilde{q}_i \right); \qquad \tilde{h}_i^{(5)} + 14 \tilde{h}_i^{(3)} = \frac{\rho_a^{3/2}}{p_a^{5/2}} X_i^{(5)}, \qquad (B77) $$

the decomposition becomes

$$ \begin{aligned} X_{s_1 s_2 s_3 s_4 s_5}^{(5)} = & \frac{1}{35} [X_{s_1}^{(5)} (\delta_{s_2 s_3} \delta_{s_4 s_5} + \delta_{s_2 s_4} \delta_{s_3 s_5} + \delta_{s_2 s_5} \delta_{s_3 s_4}) \\ & + X_{s_2}^{(5)} (\delta_{s_1 s_3} \delta_{s_4 s_5} + \delta_{s_1 s_4} \delta_{s_3 s_5} + \delta_{s_1 s_5} \delta_{s_3 s_4}) \\ & + X_{s_3}^{(5)} (\delta_{s_1 s_2} \delta_{s_4 s_5} + \delta_{s_1 s_4} \delta_{s_2 s_5} + \delta_{s_1 s_5} \delta_{s_2 s_4}) \\ & + X_{s_4}^{(5)} (\delta_{s_1 s_2} \delta_{s_3 s_5} + \delta_{s_1 s_3} \delta_{s_2 s_5} + \delta_{s_1 s_5} \delta_{s_2 s_3}) \\ & + X_{s_5}^{(5)} (\delta_{s_1 s_2} \delta_{s_3 s_4} + \delta_{s_1 s_3} \delta_{s_2 s_4} + \delta_{s_1 s_4} \delta_{s_2 s_3})]. \end{aligned} \qquad (B78) $$

As a double-check, applying contraction $\delta_{s_4 s_5}$ at the last expression yields

$$ X_{s_1 s_2 s_3}^{(5)} = \frac{1}{5} [X_{s_1}^{(5)} \delta_{s_2 s_3} + X_{s_2}^{(5)} \delta_{s_1 s_3} + X_{s_3}^{(5)} \delta_{s_1 s_2}], \qquad (B79) $$

and applying another contraction yields an identity. Note that it is not possible to perform the closure $X^{(5)} = 0$, as such a closure would be erroneous (unless $q = 0$ is prescribed as well). Instead, one needs to perform the closure at the Hermite moment $\tilde{h}_i^{(5)} = 0$, or, in other words, the correct closure is $X^{(5)} = 28(p_a/\rho_a)q$.

## B.6. Table of Useful Integrals

The Hermite polynomials allow one to build the hierarchy of the following integrals. One introduces the weight

$$ \phi^{(0)} = \frac{e^{-\tilde{c}^2/2}}{(2\pi)^{3/2}}, \qquad (B80) $$

and for any odd "m," the following integral holds:

$$ m = \text{odd}: \qquad \int \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \dots \tilde{c}_{r_m} \phi^{(0)} d^3 \tilde{c} = 0. \qquad (B81) $$

The validity of (B81) can be shown by using “common-sense” symmetries and Gaussian integration, or by rewriting the integral with pairs of Hermite polynomials, one of even order and one of odd order, $\tilde{H}^{(r_m+1)/2} \tilde{H}^{(r_m-1)/2}$ (where the result of the integration is zero), and a hierarchy of lower-order integrals that will also be odd-even pairs, yielding zero.
A particular case of the orthogonality theorem is that for any $m \neq 0$, an integral over any single Hermite polynomial with weight $\phi^{(0)}$ is zero:

$$m \neq 0: \quad \int \tilde{H}_{r_1 r_2 r_3 \dots r_m}^{(m)} \phi^{(0)} d^3 \tilde{c} = 0. \qquad (B82)$$

The two rules (B81), (B82) allow one to calculate the integrals for any even “m” number of velocities $\tilde{c}$, such as $\tilde{c}\tilde{c}\tilde{c}\tilde{c}\tilde{c}$, which would otherwise be very difficult to do just by using “common-sense” symmetries and Gaussian integration. Actually, for “m” being even, quicker than using (B82) is rewriting the integrals into $\tilde{H}^{(m/2)}\tilde{H}^{(m/2)}$ and using the orthogonality relations (B9). A very useful integral also reads

$$m < n: \quad \int \tilde{c}_{r_1} \dots \tilde{c}_{r_m} H_{s_1 \dots s_n}^{(n)} \phi^{(0)} d^3 \tilde{c} = 0, \qquad (B83)$$

the validity of which is easily shown by rewriting $\tilde{c}_{r_1} \dots \tilde{c}_{r_m}$ with $H_{r_1 \dots r_m}^{(m)}$ (where the result of the integration is zero) and a hierarchy of lower-order Hermite polynomials where the result of the integration is also zero.

It is possible to build the following table when “m” is even:

$$\begin{align*}
& \int \phi^{(0)} d^3 \tilde{c} = 1; \\
& \int \tilde{c}_i \tilde{c}_j \phi^{(0)} d^3 \tilde{c} = \delta_{ij}; \\
& \int \tilde{c}_i \tilde{c}_j \tilde{c}_k \phi^{(0)} d^3 \tilde{c} = \delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{jk} \delta_{il}; \tag{B84}
\end{align*}$$

$$\begin{align}
\int \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{c}_{s_1} \tilde{c}_{s_2} \tilde{c}_{s_3} \phi^{(0)} d^3 \tilde{c} ={}& \delta_{r_1 s_1} (\delta_{r_2 s_2} \delta_{r_3 s_3} + \delta_{r_2 s_3} \delta_{r_3 s_2}) \nonumber \\
& + \delta_{r_1 s_2} (\delta_{r_2 s_1} \delta_{r_3 s_3} + \delta_{r_2 s_3} \delta_{r_3 s_1}) + \delta_{r_1 s_3} (\delta_{r_2 s_1} \delta_{r_3 s_2} + \delta_{r_2 s_2} \delta_{r_3 s_1}) \nonumber \\
& + \delta_{r_1 r_2} (\delta_{s_1 s_2} \delta_{r_3 s_3} + \delta_{s_2 s_3} \delta_{r_3 s_1} + \delta_{s_3 s_1} \delta_{r_3 s_2}) \nonumber \\
& + \delta_{r_1 r_3} (\delta_{s_1 s_2} \delta_{r_2 s_3} + \delta_{s_2 s_3} \delta_{r_2 s_1} + \delta_{s_3 s_1} \delta_{r_2 s_2}) \nonumber \\
& + \delta_{r_2 r_3} (\delta_{s_1 s_2} \delta_{r_1 s_3} + \delta_{s_2 s_3} \delta_{r_1 s_1} + \delta_{s_3 s_1} \delta_{r_1 s_2}). \tag{B85}
\end{align}$$

These integrals can be used to obtain other useful integrals, for example, those that are valid for any (three-dimensional) vector $q$:

$$\int \tilde{c}_i \tilde{c}_j \tilde{c}_k (\tilde{c} \cdot q) \phi^{(0)} d^3 \tilde{c} = \delta_{ij} \tilde{q}_k + \delta_{jk} \tilde{q}_i + \delta_{ki} \tilde{q}_j = [\bar{\mathbf{I}} \mathbf{q}]_{ijk}^S; \qquad (B86)$$

$$\begin{align}
\int \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{c}_{s_1} \tilde{c}_{s_2} (\tilde{c} \cdot q) \phi^{(0)} d^3 \tilde{c} ={}& \delta_{r_1 s_1} (\delta_{r_2 s_2} \tilde{q}_{r_3} + \tilde{q}_{r_2} \delta_{r_3 s_2}) \nonumber \\
& + \delta_{r_1 s_2} (\delta_{r_2 s_1} \tilde{q}_{r_3} + \tilde{q}_{r_2} \delta_{r_3 s_1}) + \tilde{q}_{r_1} (\delta_{r_2 s_1} \delta_{r_3 s_2} + \delta_{r_2 s_2} \delta_{r_3 s_1}) \nonumber \\
& + \delta_{r_1 r_2} (\delta_{s_1 s_2} \tilde{q}_{r_3} + \tilde{q}_{s_2} \delta_{r_3 s_1} + \tilde{q}_{s_1} \delta_{r_3 s_2}) \nonumber \\
& + \delta_{r_1 r_3} (\delta_{s_1 s_2} \tilde{q}_{r_2} + \tilde{q}_{s_2} \delta_{r_2 s_1} + \tilde{q}_{s_1} \delta_{r_2 s_2}) \nonumber \\
& + \delta_{r_2 r_3} (\delta_{s_1 s_2} \tilde{q}_{r_1} + \tilde{q}_{s_2} \delta_{r_1 s_1} + \tilde{q}_{s_1} \delta_{r_1 s_2}). \tag{B87}
\end{align}$$

and by further contractions,

$$\int \tilde{c}_i \tilde{c}_j \tilde{c}_k \tilde{c}^2 (\tilde{c} \cdot q) \phi^{(0)} d^3 \tilde{c} = 7(\delta_{ij} \tilde{q}_k + \delta_{jk} \tilde{q}_i + \delta_{ki} \tilde{q}_j) = 7[\bar{\mathbf{I}} \mathbf{q}]_{ijk}^S; \qquad (B88)$$

$$\int \tilde{c}_i \tilde{c}_j \tilde{c}_k \left(1 - \frac{\tilde{c}^2}{5}\right) (\tilde{c} \cdot q) \phi^{(0)} d^3 \tilde{c} = -\frac{2}{5} [\bar{\mathbf{I}} \mathbf{q}]_{ijk}^S. \qquad (B89)$$

As a quick double-check of the above results, by performing the further contractions

$$\int \tilde{c}_i \tilde{c}^2 (\tilde{c} \cdot q) \phi^{(0)} d^3 \tilde{c} = 5 \tilde{q}_i; \quad \int \tilde{c}_i \tilde{c}^4 (\tilde{c} \cdot q) \phi^{(0)} d^3 \tilde{c} = 35 \tilde{q}_i, \qquad (B90)$$

which is easy to verify directly.

Similarly, for a triple contraction with any fully symmetric third-rank tensor $\bar{q}$:

$$\int \tilde{c}_i \tilde{c}_j \tilde{c}_k (\tilde{c} \tilde{c}: \bar{\mathbf{q}} \cdot \tilde{c}) \phi^{(0)} d^3 \tilde{c} = 6(q_{ijk} + \tilde{q}_i \delta_{jk} + \tilde{q}_j \delta_{ik} + \tilde{q}_k \delta_{ij}); \qquad (B91)$$

$$\int \tilde{c}_i \tilde{c}^2 (\tilde{c} \tilde{c}: \bar{\mathbf{q}} \cdot \tilde{c}) \phi^{(0)} d^3 \tilde{c} = 42 \tilde{q}_i, \qquad (B92)$$
where one defines vector $q = (1/2)\text{Tr}\bar{q}$. Finally, for any $(3 \times 3)$ matrix $\bar{\mathbf{A}}$:

$$ \int \tilde{c}_i \tilde{c}_j (\tilde{c} \tilde{c}: \bar{\mathbf{A}}) \phi^{(0)} d^3 \tilde{c} = A_{ij} + A_{ji} + (\text{Tr} \bar{\mathbf{A}}) \delta_{ij}; \qquad (B93) $$

$$ \begin{aligned} \int \tilde{c}_{r_1} \tilde{c}_{r_2} \tilde{c}_{r_3} \tilde{c}_{s_1} (\tilde{c} \tilde{c}: \bar{\mathbf{A}}) \phi^{(0)} d^3 \tilde{c} = & \delta_{r_1 s_1} (A_{r_2 r_3} + A_{r_3 r_2}) \\ & + \delta_{r_2 s_1} (A_{r_1 r_3} + A_{r_3 r_1}) + \delta_{r_3 s_1} (A_{r_1 r_2} + A_{r_2 r_1}) \\ & + \delta_{r_1 r_2} (A_{r_3 s_1} + A_{s_1 r_3} + (\text{Tr} \bar{\mathbf{A}}) \delta_{r_3 s_1}) \\ & + \delta_{r_1 r_3} (A_{r_2 s_1} + A_{s_1 r_2} + (\text{Tr} \bar{\mathbf{A}}) \delta_{r_2 s_1}) \\ & + \delta_{r_2 r_3} (A_{r_1 s_1} + A_{s_1 r_1} + (\text{Tr} \bar{\mathbf{A}}) \delta_{r_1 s_1}), \end{aligned} \qquad (B94) $$

and by further contractions,

$$ \int \tilde{c}^2 (\tilde{c} \tilde{c}: \bar{\mathbf{A}}) \phi^{(0)} d^3 \tilde{c} = 5 \text{Tr} \bar{\mathbf{A}}; \qquad (B95) $$

$$ \int \tilde{c}_i \tilde{c}_j \tilde{c}^2 (\tilde{c} \tilde{c}: \bar{\mathbf{A}}) \phi^{(0)} d^3 \tilde{c} = 7(A_{ij} + A_{ji} + (\text{Tr} \bar{\mathbf{A}}) \delta_{ij}); \qquad (B96) $$

$$ \int \tilde{c}^4 (\tilde{c} \tilde{c}: \bar{\mathbf{A}}) \phi^{(0)} d^3 \tilde{c} = 35 \text{Tr} \bar{\mathbf{A}}, \qquad (B97) $$

and so for symmetric traceless matrix $\bar{\mathbf{\Pi}}$:

$$ \int \tilde{c}_i \tilde{c}_j (\tilde{c} \tilde{c}: \bar{\mathbf{\Pi}}) \phi^{(0)} d^3 \tilde{c} = 2 \Pi_{ij}; \qquad (B98) $$

$$ \int \tilde{c}_i \tilde{c}_j \tilde{c}^2 (\tilde{c} \tilde{c}: \bar{\mathbf{\Pi}}) \phi^{(0)} d^3 \tilde{c} = 14 \Pi_{ij}; \qquad (B99) $$

$$ \begin{aligned} \int \tilde{c}_i \tilde{c}_j \tilde{c}_k (\tilde{c} \tilde{c}: \bar{\mathbf{\Pi}}) \phi^{(0)} d^3 \tilde{c} \\ = 2(\delta_{ij} \Pi_{kl} + \delta_{ik} \Pi_{jl} + \delta_{il} \Pi_{jk} + \delta_{jk} \Pi_{il} + \delta_{jl} \Pi_{ik} + \delta_{kl} \Pi_{ij}). \end{aligned} \qquad (B100) $$

A curious reader might find the following integrals useful:

$$ \begin{aligned} \int \tilde{c}_i \tilde{c}_j \tilde{c}_k \tilde{c}_l \tilde{c}^2 \phi^{(0)} d^3 \tilde{c} &= 7[\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}]; \\ \int \tilde{c}_i \tilde{c}_j \tilde{c}_k \tilde{c}_l (\tilde{c}^2 - 7) \phi^{(0)} d^3 \tilde{c} &= 0; \\ \int \tilde{c}_i \tilde{c}_j (\tilde{c}^2 - 5) \phi^{(0)} d^3 \tilde{c} &= 0; \\ \int (\tilde{c}^2 - 3) \phi^{(0)} d^3 \tilde{c} &= 0. \end{aligned} \qquad (B101) $$

## B.7. General nth-order Perturbation

The hierarchy of the simplified reducible Hermite polynomials (with tilde) can be calculated directly from (B3) as

$$ \begin{aligned} \tilde{H}_i^{(1)} &= \tilde{c}_i; \\ \tilde{H}_{ij}^{(2)} &= \tilde{c}_i \tilde{c}_j - \delta_{ij}; \\ \tilde{H}_i^{(3)} &= \tilde{c}_i (\tilde{c}^2 - 5); \\ \tilde{H}_{ij}^{(4)} &= \tilde{c}_i \tilde{c}_j (\tilde{c}^2 - 7) - \delta_{ij} (\tilde{c}^2 - 5); \\ \tilde{H}_i^{(5)} &= \tilde{c}_i (\tilde{c}^4 - 14\tilde{c}^2 + 35); \\ \tilde{H}_{ij}^{(6)} &= \tilde{c}_i \tilde{c}_j (\tilde{c}^4 - 18\tilde{c}^2 + 63) - \delta_{ij} (\tilde{c}^4 - 14\tilde{c}^2 + 35); \\ \tilde{H}_i^{(7)} &= \tilde{c}_i (\tilde{c}^6 - 27\tilde{c}^4 + 189\tilde{c}^2 - 315); \\ \tilde{H}_{ij}^{(8)} &= \tilde{c}_i \tilde{c}_j (\tilde{c}^6 - 33\tilde{c}^4 + 297\tilde{c}^2 - 693) - \delta_{ij} (\tilde{c}^6 - 27\tilde{c}^4 + 189\tilde{c}^2 - 315); \\ \tilde{H}_i^{(9)} &= \tilde{c}_i (\tilde{c}^8 - 44\tilde{c}^6 + 594\tilde{c}^4 - 2772\tilde{c}^2 + 3465), \end{aligned} \qquad (B102) $$
and the fully contracted ones for the even orders are

$$
\begin{align*}
\tilde{H}^{(2)} &= \tilde{c}^2 - 3; \\
\tilde{H}^{(4)} &= \tilde{c}^4 - 10\tilde{c}^2 + 15; \\
\tilde{H}^{(6)} &= \tilde{c}^6 - 21\tilde{c}^4 + 105\tilde{c}^2 - 105; \\
\tilde{H}^{(8)} &= \tilde{c}^8 - 36\tilde{c}^6 + 378\tilde{c}^4 - 1260\tilde{c}^2 + 945. \tag{B103}
\end{align*}
$$

The even order polynomials $\tilde{H}_{ij}^{(2n)}$ can be rewritten into

$$
\begin{align*}
\tilde{H}_{ij}^{(2)} &= \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) + \frac{\delta_{ij}}{3} \tilde{H}^{(2)}; \\
\tilde{H}_{ij}^{(4)} &= \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) (\tilde{c}^2 - 7) + \frac{\delta_{ij}}{3} \tilde{H}^{(4)}; \\
\tilde{H}_{ij}^{(6)} &= \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) (\tilde{c}^4 - 18\tilde{c}^2 + 63) + \frac{\delta_{ij}}{3} \tilde{H}^{(6)}; \\
\tilde{H}_{ij}^{(8)} &= \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) (\tilde{c}^6 - 33\tilde{c}^4 + 297\tilde{c}^2 - 693) + \frac{\delta_{ij}}{3} \tilde{H}^{(8)}. \tag{B104}
\end{align*}
$$

The orthogonality relations can be calculated as

$$
\begin{align*}
\int \phi^{(0)} \tilde{H}_i^{(1)} \tilde{H}_j^{(1)} d^3 \tilde{c} &= \delta_{ij}; & \int \phi^{(0)} \tilde{H}^{(2)} \tilde{H}^{(2)} d^3 \tilde{c} &= 6; \\
\int \phi^{(0)} \tilde{H}_i^{(3)} \tilde{H}_j^{(3)} d^3 \tilde{c} &= 10 \delta_{ij}; & \int \phi^{(0)} \tilde{H}^{(4)} \tilde{H}^{(4)} d^3 \tilde{c} &= 120; \\
\int \phi^{(0)} \tilde{H}_i^{(5)} \tilde{H}_j^{(5)} d^3 \tilde{c} &= 280 \delta_{ij}; & \int \phi^{(0)} \tilde{H}^{(6)} \tilde{H}^{(6)} d^3 \tilde{c} &= 5040; \tag{B105} \\
\int \phi^{(0)} \tilde{H}_i^{(7)} \tilde{H}_j^{(7)} d^3 \tilde{c} &= 15120 \delta_{ij}; & \int \phi^{(0)} \tilde{H}^{(8)} \tilde{H}^{(8)} d^3 \tilde{c} &= 362880; \\
\int \phi^{(0)} \tilde{H}_i^{(9)} \tilde{H}_j^{(9)} d^3 \tilde{c} &= 1330560 \delta_{ij},
\end{align*}
$$

together with

$$
\begin{align}
\hat{h}_{kl}^{(2)} \int \phi^{(0)} \tilde{H}_{ij}^{(2)} \tilde{H}_{kl}^{(2)} d^3 \tilde{c} &= 2 \tilde{h}_{ij}^{(2)}; & \hat{h}_{kl}^{(4)} \int \phi^{(0)} \tilde{H}_{ij}^{(4)} \tilde{H}_{kl}^{(4)} d^3 \tilde{c} &= 28 \hat{h}_{ij}^{(4)}; \nonumber \\
\hat{h}_{kl}^{(6)} \int \phi^{(0)} \tilde{H}_{ij}^{(6)} \tilde{H}_{kl}^{(6)} d^3 \tilde{c} &= 1008 \hat{h}_{ij}^{(6)}; & \hat{h}_{kl}^{(8)} \int \phi^{(0)} \tilde{H}_{ij}^{(8)} \tilde{H}_{kl}^{(8)} d^3 \tilde{c} &= 66528 \hat{h}_{ij}^{(8)}, \tag{B106}
\end{align}
$$

where we used traceless Hermite moments (with hat):

$$
\hat{h}_{ij}^{(2n)} = \tilde{h}_{ij}^{(2n)} - \frac{1}{3} \delta_{ij} \tilde{h}^{(2n)}, \qquad (B107)
$$

with $\tilde{h}^{(2)} = 0$ (so that $\hat{h}_{ij}^{(2)} = \tilde{h}_{ij}^{(2)}$). The perturbation of the distribution function then becomes

$$
\begin{align}
\chi_a ={}& \frac{1}{2} \tilde{h}_{ij}^{(2)} \tilde{H}_{ij}^{(2)} + \frac{1}{10} \tilde{h}_i^{(3)} \tilde{H}_i^{(3)} + \frac{1}{28} \hat{h}_{ij}^{(4)} \tilde{H}_{ij}^{(4)} + \frac{1}{120} \tilde{h}_i^{(4)} \tilde{H}_i^{(4)} + \frac{1}{280} \tilde{h}_i^{(5)} \tilde{H}_i^{(5)} \nonumber \\
& + \frac{1}{1008} \hat{h}_{ij}^{(6)} \tilde{H}_{ij}^{(6)} + \frac{1}{5040} \tilde{h}_i^{(6)} \tilde{H}_i^{(6)} + \frac{1}{15120} \tilde{h}_i^{(7)} \tilde{H}_i^{(7)} \nonumber \\
& + \frac{1}{66528} \hat{h}_{ij}^{(8)} \tilde{H}_{ij}^{(8)} + \frac{1}{362880} \tilde{h}_i^{(8)} \tilde{H}_i^{(8)} + \frac{1}{1330560} \tilde{h}_i^{(9)} \tilde{H}_i^{(9)} + \dots \tag{B108}
\end{align}
$$

The corresponding perturbation with the irreducible polynomials reads

$$
\begin{align}
\chi_a ={}& h_{ij}^{(2)} H_{ij}^{(2)} + h_i^{(3)} H_i^{(3)} + h_{ij}^{(4)} H_{ij}^{(4)} + h^{(4)} H^{(4)} + h_i^{(5)} H_i^{(5)} \nonumber \\
& + h_{ij}^{(6)} H_{ij}^{(6)} + h^{(6)} H^{(6)} + h_i^{(7)} H_i^{(7)} + h_{ij}^{(8)} H_{ij}^{(8)} + h^{(8)} H^{(8)} + h_i^{(9)} H_i^{(9)} + \dots, \tag{B109}
\end{align}
$$

i.e., no normalization constants are explicitly present. Now one can clearly see the motivation behind the definition of irreducible
polynomials of Balescu (1988), where the direct relation between irreducible (no tilde) and reducible (tilde) Hermite polynomials can
be shown to be

$$
\begin{align*}
H^{(2n)} &= \left( \frac{1}{2^n n! (2n+1)!!} \right)^{1/2} \tilde{H}^{(2n)}; \\
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} \left( \tilde{H}_{ij}^{(2n)} - \frac{1}{3} \delta_{ij} \tilde{H}^{(2n)} \right). \tag{B110}
\end{align*}
$$

Up to the normalization constants (which can be viewed as coming from the orthogonality relations), the scalar and vector
polynomials are equivalent to each other. The only difference is for the matrices $H_{ij}^{(2n)}$, where the irreducible polynomials are defined
as traceless. Multiplying (B110) by $f_a/n_\alpha$ and integrating over $d^3c$ yields analogous relations for the Hermite moments:

$$
\begin{align*}
h^{(2n)} &= \left( \frac{1}{2^n n! (2n + 1)!!} \right)^{1/2} \tilde{h}^{(2n)}; \\
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} \frac{\left( \tilde{h}_{ij}^{(2n)} - \frac{1}{3} \delta_{ij} \tilde{h}^{(2n)} \right)}{\tilde{h}_{ij}^{(2n)}}. \tag{B111}
\end{align*}
$$

Importantly, because $\hat{h}_{ij}^{(2n)}$ is traceless, multiplying (B110) and (B111) yields

$$
h_{ij}^{(2n)} H_{ij}^{(2n)} = \frac{15}{2^n (n-1)!(2n+3)!!} \hat{h}_{ij}^{(2n)} \tilde{H}_{ij}^{(2n)}. \quad (B112)
$$

The two approaches, with reducible and irreducible polynomials, thus yield the same result, with the only difference being the
location of the normalization constants. Furthermore, it feels natural to define traceless polynomials (with hat instead of tilde):

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

and on the right-hand side of (B112) replace

$$
\hat{h}_{ij}^{(2n)} \tilde{H}_{ij}^{(2n)} = \hat{h}_{ij}^{(2n)} \hat{H}_{ij}^{(2n)}, \quad (B114)
$$

which holds because $\hat{h}_{ij}^{(2n)}$ is traceless. The main advantage of introducing the polynomials $\hat{H}_{ij}^{(2n)}$ is that, instead of calculating $\hat{h}_{ij}^{(2n)}$ from its definition (B107), one can directly define

$$
\hat{h}_{ij}^{(2n)} = \frac{1}{n_a} \int f_a \hat{H}_{ij}^{(2n)} d^3 c. \qquad (B115)
$$

Then the two approaches are indeed equivalent, because the same polynomials are used, with the location of the normalization constants being an ad hoc choice.

From the Appendix of Balescu (1988), one can guess and then verify the following generalizations for the reducible polynomials:

$$
\tilde{H}^{(2n)} = \sum_{m=0}^{n} (-1)^{m+n} \frac{n!}{m!(n-m)!} \frac{(2n+1)!!}{(2m+1)!!} \tilde{c}^{2m}; \quad (B116)
$$

$$
\tilde{H}_i^{(2n+1)} = \tilde{c}_i \sum_{m=0}^{n} (-1)^{m+n} \frac{n!}{m!(n-m)! (2m+3)!!} \tilde{c}^{2m}; \quad (B117)
$$

$$
\begin{align}
\tilde{H}_{ij}^{(2n)} &= \tilde{c}_i \tilde{c}_j \left( \sum_{m=0}^{n-1} (-1)^{m+n-1} \frac{(n-1)!}{m!(n-m-1)!} \frac{(2n+3)!!}{(2m+5)!!} \tilde{c}^{2m} \right) - \delta_{ij} \frac{\tilde{H}_k^{(2n-1)}}{\tilde{c}_k} \nonumber \\
&= \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) \left( \sum_{m=0}^{n-1} (-1)^{m+n-1} \frac{(n-1)!}{m!(n-m-1)!} \frac{(2n+3)!!}{(2m+5)!!} \tilde{c}^{2m} \right) + \frac{\delta_{ij}}{3} \tilde{H}^{(2n)}; \tag{B118}
\end{align}
$$

$$
\hat{H}_{ij}^{(2n)} = \left( \tilde{c}_i \tilde{c}_j - \frac{\delta_{ij}}{3} \tilde{c}^2 \right) \left( \sum_{m=0}^{n-1} (-1)^{m+n-1} \frac{(n-1)!}{m!(n-m-1)!} \frac{(2n+3)!!}{(2m+5)!!} \tilde{c}^{2m} \right). \quad (B119)
$$
Applying a trace at (B118) yields (B116). Similarly, the orthogonal relations are

$$ \int \phi^{(0)} \tilde{H}^{(2n)} \tilde{H}^{(2m)} d^3 \tilde{c} = 2^n n! (2n + 1)!! \delta_{nm}; \qquad (B120) $$

$$ \int \phi^{(0)} \tilde{H}_{i}^{(2n+1)} \tilde{H}_{j}^{(2m+1)} d^3 \tilde{c} = \frac{2^n n! (2n + 3)!!}{3} \delta_{ij} \delta_{nm}; \qquad (B121) $$

$$ \int \phi^{(0)} \tilde{H}_{ij}^{(2n)} \tilde{H}_{kl}^{(2m)} d^3 \tilde{c} = \frac{2^{n-1}(n-1)!(2n+1)!!}{15} [(2n+3)(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) + 2(n-1)\delta_{ij}\delta_{kl}] \delta_{nm}; \qquad (B122) $$

$$ \int \phi^{(0)} \hat{H}_{ij}^{(2n)} \hat{H}_{kl}^{(2m)} d^3 \tilde{c} = \frac{2^{n-1}(n-1)!(2n+3)!!}{15} \left[ \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} - \frac{2}{3}\delta_{ij}\delta_{kl} \right] \delta_{nm}, \qquad (B123) $$

and applying $\delta_{ij}\delta_{kl}$ on (B122) recovers (B120). Note that the orders of the Hermite moments “m” and “n” are one-dimensional and $\delta_{nn} = 1$. In contrast, for the indices, $\delta_{ii} = 3$ applies. Also note that $n! = n!(n-1)!!$ and $2^n n! = (2n)!!$, implying $2^n n! (2n+1)!! = (2n+1)!$. Applying traceless $\hat{H}_{kl}^{(2n)}$ on (B122) or (B123) yields the orthogonal relation

$$ \begin{aligned} \hat{h}_{kl}^{(2n)} \int \phi^{(0)} \hat{H}_{ij}^{(2n)} \hat{H}_{kl}^{(2n)} d^3 \tilde{c} &= \hat{h}_{kl}^{(2n)} \int \phi^{(0)} \tilde{H}_{ij}^{(2n)} \tilde{H}_{kl}^{(2n)} d^3 \tilde{c} \\ &= \frac{2^n (n-1)!(2n+3)!!}{15} \hat{h}_{ij}^{(2n)}. \end{aligned} \qquad (B124) $$

Finally, the general perturbation then can be written as

$$ \chi_a = \sum_{n=1}^{\infty} \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], \qquad (B125) $$

where for the first term $\tilde{h}^{(2)} = 0$ (and so $\hat{h}_{ij}^{(2)} = \tilde{h}_{ij}^{(2)}$). Alternatively, $\hat{h}_{ij}^{(2n)} \hat{H}_{ij}^{(2n)} = \hat{h}_{ij}^{(2n)} \tilde{H}_{ij}^{(2n)}$. The perturbation (B125) is equivalent to a perturbation with irreducible polynomials:

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

where again $h^{(2)} = 0$.

## B.8. Hierarchy of MHD Hermite Closures

Let us use the third-order moment $X_i^{(3)} = 2q_i$ instead of the heat flux, so that no additional factors are present (also note that $X^{(2)} = 3p$). The even order moments are decomposed according to

$$ X_{ij}^{(2n)} = \frac{\delta_{ij}}{3} X^{(2n)} + \Pi_{ij}^{(2n)}, \qquad (B127) $$

where the scalar part $X^{(2n)}$ is further decomposed into its Maxwellian “core” and perturbation $\tilde{X}^{(2n)}$ (with wide tilde), as

$$ X^{(2n)} = (2n + 1)!! \frac{p^n}{\rho^{n-1}} + \tilde{X}^{(2n)}, \qquad (B128) $$

so, for example,

$$ X^{(4)} = 15 \frac{p^2}{\rho} + \tilde{X}^{(4)}; \qquad X^{(6)} = 105 \frac{p^3}{\rho^2} + \tilde{X}^{(6)}; \qquad X^{(8)} = 945 \frac{p^4}{\rho^3} + \tilde{X}^{(8)}. \qquad (B129) $$
Then, by using the Hermite polynomials (B102)–(B104), one calculates the hierarchy of the Hermite moments:

$$
\begin{align*}
\tilde{h}_i^{(3)} &= \frac{\rho^{1/2}}{p^{3/2}} X_i^{(3)}; & \tilde{h}^{(4)} &= \frac{\rho}{p} \tilde{X}^{(4)}; \\
\tilde{h}_i^{(5)} &= \frac{\rho^{1/2}}{p^{3/2}} \left[ \frac{\rho}{p} X_i^{(5)} - 14X_i^{(3)} \right]; & \tilde{h}^{(6)} &= \frac{\rho}{p^2} \left[ \frac{\rho}{p} \tilde{X}^{(6)} - 21\tilde{X}^{(4)} \right]; \\
\tilde{h}_i^{(7)} &= \frac{\rho^{1/2}}{p^{3/2}} \left[ \frac{\rho^2}{p^2} X_i^{(7)} - 27\frac{\rho}{p}X_i^{(5)} + 189X_i^{(3)} \right]; \\
\tilde{h}^{(8)} &= \frac{\rho}{p^2} \left[ \frac{\rho^2}{p^2} \tilde{X}^{(8)} - 36\frac{\rho}{p}\tilde{X}^{(6)} + 378\tilde{X}^{(4)} \right]; \\
\tilde{h}_i^{(9)} &= \frac{\rho^{1/2}}{p^{3/2}} \left[ \frac{\rho^3}{p^3} X_i^{(9)} - 44\frac{\rho^2}{p^2} X_i^{(7)} + 594\frac{\rho}{p} X_i^{(5)} - 2772X_i^{(3)} \right], \tag{B130}
\end{align*}
$$

together with

$$
\begin{align*}
\hat{h}_{ij}^{(4)} &= \frac{1}{p} \left[ \frac{\rho}{p} \Pi_{ij}^{(4)} - 7\Pi_{ij}^{(2)} \right]; \\
\hat{h}_{ij}^{(6)} &= \frac{1}{p} \left[ \frac{\rho^2}{p^2} \Pi_{ij}^{(6)} - 18\frac{\rho}{p} \Pi_{ij}^{(4)} + 63\Pi_{ij}^{(2)} \right]; \\
\hat{h}_{ij}^{(8)} &= \frac{1}{p} \left[ \frac{\rho^3}{p^3} \Pi_{ij}^{(8)} - 33\frac{\rho^2}{p^2} \Pi_{ij}^{(6)} + 297\frac{\rho}{p} \Pi_{ij}^{(4)} - 693\Pi_{ij}^{(2)} \right]. \tag{B131}
\end{align*}
$$

Prescribing the last retained Hermite moment to be zero then yields the corresponding fluid closures that are summarized in Section 8.6, Tables 3 and 4.

B.8.1. Propagation along the B-field (Ion-acoustic Mode)

For a propagation parallel to the mean magnetic field that is applied in the z-direction, linearized equations without collisions read

$$
\begin{align*}
& \frac{\partial \rho}{\partial t} + \rho_0 \partial_z u_z = 0; && \frac{\partial u_z}{\partial t} + \frac{1}{\rho_0} \partial_z p = 0; \\
& \frac{\partial p}{\partial t} + \frac{5}{3} p_0 \partial_z u_z + \frac{2}{3} \partial_z q_z = 0; \\
& \frac{\partial q_z}{\partial t} + \frac{1}{6} \partial_z X^{(4)} - \frac{5}{2} \frac{p_0}{\rho_0} \partial_z p = 0; \\
& \frac{\partial X^{(4)}}{\partial t} + \partial_z X_z^{(5)} + \frac{7}{3} X_0^{(4)} \partial_z u_z = 0; \\
& \frac{\partial X_z^{(5)}}{\partial t} + \frac{1}{3} \partial_z X^{(6)} - \frac{7}{3} \frac{X_0^{(4)}}{\rho_0} \partial_z p = 0, \tag{B132}
\end{align*}
$$

where all of the variables are scalars. We are neglecting collisions and viscosities, to make direct comparison with the CGL model in the next section. The even order moments are decomposed into a Maxwellian “core” and tilde perturbations with (B129), so that their mean values are $X_0^{(4)} = 15p_0^2/\rho_0$ and $X_0^{(6)} = 105p_0^3/\rho_0^2$. These moments are thus linearized according to

$$
X^{(4)} \stackrel{\text{lin.}}{=} X_0^{(4)} \left( 2 \frac{p}{p_0} - \frac{\rho}{\rho_0} \right) + \tilde{X}^{(4)}; \quad X^{(6)} \stackrel{\text{lin.}}{=} X_0^{(6)} \left( 3 \frac{p}{p_0} - 2 \frac{\rho}{\rho_0} \right) + \tilde{X}^{(6)}, \qquad \text{(B133)}
$$
<table><caption>**Table 5**
Summary of the Hermite Closures and Corresponding Dispersion Relations for the Parallel Propagating Ion-acoustic Mode (Electrons Are Cold), Where $\zeta = \omega/(|k_{\parallel}|v_{th})$</caption><thead><tr><th>Closure</th><th>Dispersion Relation</th><th>Solution ±ζ =</th></tr></thead><tbody><tr><td>$\tilde{h}_z^{(3)} = 0$</td><td>$\zeta^2 - 5/6 = 0$</td><td>0.913</td></tr><tr><td>$\tilde{h}_z^{(4)} = 0$</td><td>$\zeta^4 - (5/3)\zeta^2 + (5/12) = 0$</td><td>0.553; 1.166</td></tr><tr><td>$\tilde{h}_z^{(5)} = 0$</td><td>$\zeta^4 - (7/3)\zeta^2 + (35/36) = 0$</td><td>0.737; 1.338</td></tr><tr><td>$\tilde{h}_z^{(6)} = 0$</td><td>$\zeta^6 - (7/2)\zeta^4 + (35/12)\zeta^2 - (35/72) = 0$</td><td>0.471; 0.966; 1.531</td></tr><tr><td>$X_z^{(5)} = 0$</td><td>$\zeta^4 - (35/36) = 0$</td><td>0.99; 0.99 i</td></tr><tr><td>$\tilde{X}_z^{(6)} = 0$</td><td>$\zeta^6 - (35/12)\zeta^2 + (35/36) = 0$</td><td>0.59; 1.23; 1.36 i</td></tr></tbody></table>

**Note.** With the Hermite closures (the upper half of the table), no spurious instabilities are present. Unphysical instabilities appear if one prescribes erroneous fluid closures at the last retained moment $X_z^{(5)} = 0$ or $\tilde{X}_z^{(6)} = 0$ (the lower half of the table). However, if one prescribes at the same time, $\tilde{X}_z^{(6)} = 0$ and $\tilde{X}_z^{(4)} = 0$, the system is again well defined, with a dispersion relation equivalent to the closure $\tilde{h}_z^{(4)} = 0$.

and the last three equations of (B132) then become

$$ 
\begin{aligned}
& \frac{\partial q_z}{\partial t} + \frac{1}{6} \partial_z \tilde{X}^{(4)} + \frac{5}{2} \frac{p_0}{\rho_0} \left( \partial_z p - \frac{p_0}{\rho_0} \partial_z \rho \right) = 0; \\
& \frac{\partial \tilde{X}^{(4)}}{\partial t} + \partial_z X_z^{(5)} - 20 \frac{p_0}{\rho_0} \partial_z q_z = 0; \\
& \frac{\partial X_z^{(5)}}{\partial t} + \frac{1}{3} \partial_z \tilde{X}^{(6)} + 70 \frac{p_0^2}{\rho_0^2} \left( \partial_z p - \frac{p_0}{\rho_0} \partial_z \rho \right) = 0.
\end{aligned}
\quad (B134) $$

Prescribing a closure at the last retained moment yields the dispersion relations in the variable $\zeta = \omega/(|k_{\parallel}|v_{th})$ that are summarized in Table 5.

The example clearly demonstrates that Landau fluid closures are actually not required to go beyond the fourth-order moment, which contradicts a claim in the last paragraph of Hunana et al. (2018) and also claims in various parts of Hunana et al. (2019a, 2019b; see, e.g., Section 12.2 in Part 1). Obviously, the closures $X_z^{(5)} = 0$ or $\tilde{X}_z^{(6)} = 0$ are not allowed by the fluid hierarchy (unless $q_z = 0$ or $\tilde{X}_z^{(4)} = 0$ as well). Instead, for moments of order $n \ge 5$, one needs to construct “classical” closures at the Hermite moments. Nevertheless, all of the Landau fluid closures reported in the above papers are constructed correctly.

Out of curiosity, prescribing closures with a free parameter “a” as $X_z^{(5)} = 28a(p_0/\rho_0)q_z$ or $\tilde{X}_z^{(6)} = 21a(p_0/\rho_0)\tilde{X}_z^{(4)}$ yields the following dispersion relations:

$$ \zeta^4 - \frac{7a}{3}\zeta^2 + \frac{35a}{18} - \frac{35}{36} = 0; \quad (B135) $$

$$ \zeta^6 - \frac{7a}{2}\zeta^4 + \left(\frac{35a}{6} - \frac{35}{12}\right)\zeta^2 - \frac{35a}{24} + \frac{35}{36} = 0. \quad (B136) $$

The $X_z^{(5)}$ closure with dispersion relation (B135) yields an instability for $a < 1/2$, and the $\tilde{X}_z^{(6)}$ closure with (B136) yields an instability for $a < 2/3$. There are therefore a lot of closures that do not create these unphysical instabilities.

Finally, the situation is saved by completely decoupling the odd and even moments, prescribing $\tilde{X}_z^{(6)} = \tilde{X}_z^{(4)} = 0$, for example, so that Equations (B134) are replaced by

$$ 
\begin{aligned}
& \frac{\partial q_z}{\partial t} + \frac{5}{2} \frac{p_0}{\rho_0} \left( \partial_z p - \frac{p_0}{\rho_0} \partial_z \rho \right) = 0; \\
& \frac{\partial X_z^{(5)}}{\partial t} + 70 \frac{p_0^2}{\rho_0^2} \left( \partial_z p - \frac{p_0}{\rho_0} \partial_z \rho \right) = 0.
\end{aligned}
\quad (B137) $$

The dispersion relation of this model is equivalent to the closure $\tilde{h}_z^{(4)} = 0$.
## B.9. Hierarchy of CGL (Parallel) Hermite Closures

The hierarchy of one-dimensional Hermite polynomials calculates (with weight $\exp(-\tilde{c}^2/2)$)

$$
\begin{aligned}
H^{(1)} &= \tilde{c}; \\
H^{(2)} &= \tilde{c}^2 - 1; \\
H^{(3)} &= \tilde{c}(\tilde{c}^2 - 3); \\
H^{(4)} &= \tilde{c}^4 - 6\tilde{c}^2 + 3; \\
H^{(5)} &= \tilde{c}(\tilde{c}^4 - 10\tilde{c}^2 + 15); \\
H^{(6)} &= \tilde{c}^6 - 15\tilde{c}^4 + 45\tilde{c}^2 - 15; \\
H^{(7)} &= \tilde{c}(\tilde{c}^6 - 21\tilde{c}^4 + 105\tilde{c}^2 - 105); \\
H^{(8)} &= \tilde{c}^8 - 28\tilde{c}^6 + 210\tilde{c}^4 - 420\tilde{c}^2 + 105; \\
H^{(9)} &= \tilde{c}(\tilde{c}^8 - 36\tilde{c}^6 + 378\tilde{c}^4 - 1260\tilde{c}^2 + 945),
\end{aligned}
\qquad (B138)
$$

further yielding the following hierarchy of Hermite moments:

$$
\begin{aligned}
h^{(1)} &= 0; & h^{(2)} &= 0; \\
h^{(3)} &= \frac{\rho^{1/2}}{p^{3/2}} X^{(3)}; & h^{(4)} &= \frac{\rho}{p^2} \tilde{X}^{(4)}; \\
h^{(5)} &= \frac{\rho^{1/2}}{p^{3/2}} \left( \frac{\rho}{p} X^{(5)} - 10X^{(3)} \right); & h^{(6)} &= \frac{\rho}{p^2} \left( \frac{\rho}{p} \tilde{X}^{(6)} - 15\tilde{X}^{(4)} \right); \\
h^{(7)} &= \frac{\rho^{1/2}}{p^{3/2}} \left( \frac{\rho^2}{p^2} X^{(7)} - 21\frac{\rho}{p} X^{(5)} + 105X^{(3)} \right); \\
h^{(8)} &= \frac{\rho}{p^2} \left( \frac{\rho^2}{p^2} \tilde{X}^{(8)} - 28\frac{\rho}{p} \tilde{X}^{(6)} + 210\tilde{X}^{(4)} \right); \\
h^{(9)} &= \frac{\rho^{1/2}}{p^{3/2}} \left( \frac{\rho^3}{p^3} X^{(9)} - 36\frac{\rho^2}{p^2} X^{(7)} + 378\frac{\rho}{p} X^{(5)} - 1260X^{(3)} \right),
\end{aligned}
\qquad (B139)
$$

where the even moments were separated into

$$
\begin{aligned}
X^{(4)} &= 3 \frac{p^2}{\rho} + \tilde{X}^{(4)}; & X^{(6)} &= 15 \frac{p^3}{\rho^2} + \tilde{X}^{(6)}; & X^{(8)} &= 105 \frac{p^4}{\rho^3} + \tilde{X}^{(8)}; \\
X^{(2n)} &= (2n - 1)!! \frac{p^n}{\rho^{n-1}} + \tilde{X}^{(2n)}.
\end{aligned}
\qquad (B140)
$$

This yields the hierarchy of Hermite closures summarized in Table 6. Note the difference between (B140) and the isotropic (MHD) decomposition (B128) (in the three-dimensional CGL geometry, one typically uses the notation $X^{(4)} = r_{\parallel}$).

Hermite polynomials (B138) can be written in a general form:

$$
\begin{aligned}
H^{(2n+1)} &= \sum_{m=0}^{n} (-1)^{n-m} \frac{(2n+1)!}{2^{n-m}(2m+1)!(n-m)!} \tilde{c}^{2m+1}; \\
H^{(2n)} &= \sum_{m=0}^{n} (-1)^{n-m} \frac{(2n)!}{2^{n-m}(2m)!(n-m)!} \tilde{c}^{2m}.
\end{aligned}
\qquad (B141)
$$

Then it can be shown that prescribing the Hermite closure $h^{(2n+1)} = 0$ or $h^{(2n)} = 0$ is equivalent to prescribing the fluid closure

$$
\begin{aligned}
X^{(2n+1)} &= \sum_{m=1}^{n-1} (-1)^{n-m+1} \frac{(2n+1)!}{2^{n-m}(2m+1)!(n-m)!} \left(\frac{p}{\rho}\right)^{n-m} X^{(2m+1)}; \\
\tilde{X}^{(2n)} &= \sum_{m=2}^{n-1} (-1)^{n-m+1} \frac{(2n)!}{2^{n-m}(2m)!(n-m)!} \left(\frac{p}{\rho}\right)^{n-m} \tilde{X}^{(2m)}.
\end{aligned}
\qquad (B142)
$$

By using Equations (12.49)–(12.54) from Hunana et al. (2019b), we calculate the corresponding dispersion relations, which are summarized in Table 7.
<table><caption>**Table 6**
Summary of Hermite Closures for Parallel CGL Moments, Together with Corresponding Fluid Closures</caption><thead><tr><th>Hermite Closures</th><th>Fluid Closures</th></tr></thead><tbody><tr><td>h<sup>(3)</sup> = 0</td><td>X<sup>(3)</sup> = 0</td></tr><tr><td>h<sup>(4)</sup> = 0</td><td>X&#771;<sup>(4)</sup> = 0</td></tr><tr><td>h<sup>(5)</sup> = 0</td><td>X<sup>(5)</sup> = 10<sup>p</sup>/<sub>p</sub>X<sup>(3)</sup></td></tr><tr><td>h<sup>(6)</sup> = 0</td><td>X&#771;<sup>(6)</sup> = 15<sup>p</sup>/<sub>p</sub>X&#771;<sup>(4)</sup></td></tr><tr><td>h<sup>(7)</sup> = 0</td><td>X<sup>(7)</sup> = 21<sup>p</sup>/<sub>p</sub>X<sup>(5)</sup> - 105<sup>p<sup>2</sup></sup>/<sub>p<sup>2</sup></sub>X<sup>(3)</sup></td></tr><tr><td>h<sup>(8)</sup> = 0</td><td>X&#771;<sup>(8)</sup> = 28<sup>p</sup>/<sub>p</sub>X&#771;<sup>(6)</sup> - 210<sup>p<sup>2</sup></sup>/<sub>p<sup>2</sup></sub>X&#771;<sup>(4)</sup></td></tr><tr><td>h<sup>(9)</sup> = 0</td><td>X<sup>(9)</sup> = 36<sup>p</sup>/<sub>p</sub>X<sup>(7)</sup> - 378<sup>p<sup>2</sup></sup>/<sub>p<sup>2</sup></sub>X<sup>(5)</sup> + 1260<sup>p<sup>3</sup></sup>/<sub>p<sup>3</sup></sub>X<sup>(3)</sup></td></tr></tbody></table>

**Note.** The usual parallel heat flux $q_\parallel = X^{(3)}$. Note that beyond the fourth-order moment, both classes start to differ. A general form corresponding to $h^{(2n+1)} = 0$ and $h^{(2n)} = 0$ is given by (B142).

<table><caption>**Table 7**
Summary of Hermite Closures and Corresponding Dispersion Relations for the Parallel Propagating Ion-acoustic Mode (Electrons Are Cold), Where $\zeta = \omega/(k_\parallel v_\text{th})$</caption><thead><tr><th>Closure</th><th>Dispersion Relation</th><th>Solution ±ζ =</th></tr></thead><tbody><tr><td>h<sup>(3)</sup> = 0</td><td>ζ<sup>2</sup> - 3/2 = 0</td><td>1.225</td></tr><tr><td>h<sup>(4)</sup> = 0</td><td>ζ<sup>4</sup> - 3ζ<sup>2</sup> + 3/4 = 0</td><td>0.525; 1.651</td></tr><tr><td>h<sup>(5)</sup> = 0</td><td>ζ<sup>4</sup> - 5ζ<sup>2</sup> + 15/4 = 0</td><td>0.959; 2.020</td></tr><tr><td>h<sup>(6)</sup> = 0</td><td>ζ<sup>6</sup> - (15/2)ζ<sup>4</sup> + (45/4)ζ<sup>2</sup> - 15/8 = 0</td><td>0.436; 1.336; 2.351</td></tr><tr><td>h<sup>(7)</sup> = 0</td><td>ζ<sup>6</sup> - (21/2)ζ<sup>4</sup> + (105/4)ζ<sup>2</sup> - 105/8 = 0</td><td>0.816; 1.674; 2.652</td></tr><tr><td>h<sup>(8)</sup> = 0</td><td>ζ<sup>8</sup> - 14ζ<sup>6</sup> + (105/2)ζ<sup>4</sup> - (105/2)ζ<sup>2</sup> + 105/16 = 0</td><td>0.381; 1.157; 1.982; 2.931</td></tr><tr><td>h<sup>(9)</sup> = 0</td><td>ζ<sup>8</sup> - 18ζ<sup>6</sup> + (189/2)ζ<sup>4</sup> - (315/2)ζ<sup>2</sup> + 945/16 = 0</td><td>0.724; 1.469; 2.267; 3.191</td></tr></tbody></table>

**Note.** No spurious instabilities are present. Spurious instabilities occur if one prescribes the closures $\tilde{X}^{(2n)} = 0$ or $X^{(2n+1)} = 0$ at the last retained moment.

Curiously, from Hunana et al. (2019a), the not “well-behaved” Padé approximants of the plasma dispersion function $R(\zeta)$ that contain no Landau damping read

$$ 
\begin{aligned}
R_{4,5}(\zeta) &= \frac{1 - (2/3)\zeta^2}{1 - 4\zeta^2 + (4/3)\zeta^4}; \\
R_{6,9}(\zeta) &= \frac{1 - (8/5)\zeta^2 + (4/15)\zeta^4}{1 - 6\zeta^2 + 4\zeta^4 - (8/15)\zeta^6}; \\
R_{8,13}(\zeta) &= \frac{1 - (94/35)\zeta^2 + (20/21)\zeta^4 - (8/105)\zeta^6}{1 - 8\zeta^2 + 8\zeta^4 - (32/15)\zeta^6 + (16/105)\zeta^8}.
\end{aligned} 
\quad (B143) $$

Comparing (B143) with Table 7, one comes to the nonobvious observation that the denominators of the above approximants are equal to the dispersion relations obtained with the Hermite closures $h^{(4)} = 0$, $h^{(6)} = 0$, and $h^{(8)} = 0$. This observation is analogous to the Landau fluid closures when electrons are cold; see Equation (3.358) of Hunana et al. (2019a). Thus, it is expected that for proton-electron plasma with finite temperatures (and with electron inertia retained), these three dispersion relations will be equivalent to

$$ \frac{T_{\parallel e}^{(0)}}{T_{\parallel p}^{(0)}} R_{n,n'}(\zeta_p) + R_{n,n'}(\zeta_e) = 0, \quad (B144) $$

which we did not verify.

## Appendix C
Evolution Equations for the 22-moment Model

Here we use evolution Equations (A4)–(A9), and by applying the contractions at these equations we obtain the 22-moment model in detail. The pressure tensor is decomposed as $p_{ij}^a = p_a \delta_{ij} + \Pi_{ij}^{(2)}$, where the scalar pressure $p_a = p_{ii}^a/3$. Instead of considering the full moments $X_{ijk}^{(3)}$, $X_{ijkl}^{(4)}$, $X_{ijklm}^{(5)}$, and $X_{ijklmn}^{(6)}$, one only considers the contracted vectors and matrices

$$ X_i^{a(3)} = X_{ij}^{a(3)}; \quad X_{ij}^{a(4)} = X_{ijk}^{a(4)}; \quad X_i^{a(5)} = X_{ijk}^{a(5)}; \quad X_{ij}^{a(6)} = X_{ijk}^{a(6)}. \quad (C1) $$
The even order moments are decomposed by separating the traceless viscosity tensors $\Pi_{ij}^{(2n)}$:

$$X_{ij}^{a(4)} = \frac{\delta_{ij}}{3} X^{a(4)} + \Pi_{ij}^{a(4)}; \quad X_{ij}^{a(6)} = \frac{\delta_{ij}}{3} X^{a(6)} + \Pi_{ij}^{a(6)}, \qquad (C2)$$

where the fully contracted (scalsars) $X^{a(4)} = X_{ijij}^{a(4)}$, $X^{a(6)} = X_{ijijkk}^{a(6)}$. The scalars are further decomposed into their “Maxwellian core” and a perturbation around this core (which is denoted by wide tilde):

$$X_a^{(4)} = \mathrm{Tr}\mathrm{Tr}\bar{\bar{\mathbf{r}}}_a = 15 \frac{p_a^2}{\rho_a} + \tilde{\mathbf{X}}_a^{(4)}; \quad X_a^{(6)} = \mathrm{Tr}\mathrm{Tr}\bar{\bar{\mathbf{X}}}_a = 105 \frac{p_a^3}{\rho_a^2} + \tilde{\mathbf{X}}_a^{(6)}. \qquad (C3)$$

As in Braginskii (1965), we use notation with the Boltzmann constant $k_B = 1$, and the temperature is defined as $T_a = p_a/n_a$. Note that $m_a/T_a = \rho_a/p_a$.

### C.1. Decomposition of Moments

The heat flux tensor $q_{ijk}$ and moments $X_{ijkl}^{(4)}$, $X_{ijklm}^{(5)}$ are decomposed according to (see Appendix B)

$$q_{ijk}^{a} = \frac{2}{5} [\bar{\bar{\mathbf{I}}} \mathbf{q}^{a}]_{ijk}^{S}; \qquad (C4)$$

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

$$X_{s_1 s_2 s_3 s_4 s_5}^{a(5)} = \frac{1}{35} [X_{s_1}^{a(5)} (\delta_{s_2 s_3} \delta_{s_4 s_5} + \delta_{s_2 s_4} \delta_{s_3 s_5} + \delta_{s_2 s_5} \delta_{s_3 s_4}) \\ + X_{s_2}^{a(5)} (\delta_{s_1 s_3} \delta_{s_4 s_5} + \delta_{s_1 s_4} \delta_{s_3 s_5} + \delta_{s_1 s_5} \delta_{s_3 s_4}) \\ + X_{s_3}^{a(5)} (\delta_{s_1 s_2} \delta_{s_4 s_5} + \delta_{s_1 s_4} \delta_{s_2 s_5} + \delta_{s_1 s_5} \delta_{s_2 s_4}) \\ + X_{s_4}^{a(5)} (\delta_{s_1 s_2} \delta_{s_3 s_5} + \delta_{s_1 s_3} \delta_{s_2 s_5} + \delta_{s_1 s_5} \delta_{s_2 s_3}) \\ + X_{s_5}^{a(5)} (\delta_{s_1 s_2} \delta_{s_3 s_4} + \delta_{s_1 s_3} \delta_{s_2 s_4} + \delta_{s_1 s_4} \delta_{s_2 s_3})], \qquad (C6)$$

where the highest-order irreducible parts of moments (C4)–(C6), denoted as $\sigma_{ijk}^{(3)\prime}$, $\sigma_{ijkl}^{(4)\prime}$, and $\sigma_{ijklm}^{(5)\prime}$, are neglected (which provides the reduction from the 56-moment model to the 22-moment model).

### C.2. Evolution Equation for the Scalar Pressure $p_a$

By using the decomposition $\bar{\bar{\mathbf{p}}}_a = p_a \bar{\bar{\mathbf{I}}} + \bar{\bar{\mathbf{\Pi}}}_a^{(2)}$, the evolution equation for the scalar pressure $p_a$ is obtained by applying (1/3)Tr on the pressure tensor Equation (A6), yielding

$$\frac{\partial p_a}{\partial t} + \mathbf{u}_a \cdot \nabla 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{\mathbf{\Pi}}}_a^{(2)} : (\nabla \mathbf{u}_a) = \frac{1}{3} \mathrm{Tr} \bar{\bar{\mathbf{Q}}}_a^{(2)} = \frac{2}{3} Q_a. \qquad (C7)$$

Alternatively, using temperature $T_a = p_a/n_a$ yields the following equation:

$$\frac{3}{2} n_a \frac{d_a T_a}{dt} + p_a \nabla \cdot \mathbf{u}_a + \nabla \cdot \mathbf{q}_a + \bar{\bar{\mathbf{\Pi}}}_a^{(2)} : (\nabla \mathbf{u}_a) = \frac{1}{2} \mathrm{Tr} \bar{\bar{\mathbf{Q}}}_a^{(2)} = Q_a, \qquad (C8)$$

which identifies with Equation (2.3) of Braginskii (1965). The collisional energy exchange rates

$$Q_a = \frac{m_a}{2} \int |c_a|^2 C(f_a) d^3 v. \qquad (C9)$$
### C.3. Evolution Equation for the Viscosity Tensor $\bar{\Pi}_a^{(2)}$

The evolution equation for the usual viscosity tensor is obtained by subtracting $\bar{\mathcal{I}}$ times (C7) from (A6), yielding

$$ \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{\mathcal{I}} (\bar{\Pi}_a^{(2)} : \nabla \mathbf{u}_a) + \nabla \cdot \bar{\mathbf{q}}_a - \frac{2}{3} \bar{\mathcal{I}} \nabla \cdot \mathbf{q}_a \\ + p_a \left[ (\nabla \mathbf{u}_a)^S - \frac{2}{3} \bar{\mathcal{I}} \nabla \cdot \mathbf{u}_a \right] = \bar{\mathcal{Q}}_a^{(2)} - \frac{\bar{\mathcal{I}}}{3} \mathrm{Tr} \bar{\mathcal{Q}}_a^{(2)}. \qquad (C10) $$

It is possible to define the well-known rate-of-strain tensor:

$$ \bar{W}_a = (\nabla \mathbf{u}_a)^S - \frac{2}{3} \bar{\mathcal{I}} \nabla \cdot \mathbf{u}_a. \qquad (C11) $$

Equation (C10) is exact. Using the heat flux decomposition (C4) yields $\nabla \cdot \bar{\mathbf{q}}_a = (2/5)((\nabla \mathbf{q}_a)^S + \bar{\mathcal{I}} \nabla \cdot \mathbf{q}_a)$, and so Equation (C10) becomes

$$ \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{\mathcal{I}} (\bar{\Pi}_a^{(2)} : \nabla \mathbf{u}_a) + \frac{2}{5} \left[ (\nabla \mathbf{q}_a)^S - \frac{2}{3} \bar{\mathcal{I}} \nabla \cdot \mathbf{q}_a \right] \\ + p_a \bar{W}_a = \bar{\mathcal{Q}}_a^{(2)} = \bar{\mathcal{Q}}_a^{(2)}' \equiv \bar{\mathcal{Q}}_a^{(2)} - \frac{\bar{\mathcal{I}}}{3} \mathrm{Tr} \bar{\mathcal{Q}}_a^{(2)}, \qquad (C12) $$

which, for example, identifies with Equations (39)–(40) of Schunk (1977). It is possible to define

$$ \bar{W}_a^q = \frac{2}{5} \left[ (\nabla \mathbf{q}_a)^S - \frac{2}{3} \bar{\mathcal{I}} \nabla \cdot \mathbf{q}_a \right], \qquad (C13) $$

where we use a heat flux superscript “q” to differentiate it from (C11). As a double check, applying a trace on (C12) shows that both sides are zero.

### C.4. Evolution Equation for the Heat Flux Vector $\mathbf{q}_a$

The evolution equation for $\mathbf{q}_a$ is obtained by applying (1/2)Tr on (A7), yielding

$$ \frac{d_a \mathbf{q}_a}{dt} + \mathbf{q}_a \nabla \cdot \mathbf{u}_a + \mathbf{q}_a \cdot \nabla \mathbf{u}_a + \bar{\mathbf{q}}_a : \nabla \mathbf{u}_a + \Omega_a \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{1}{2} \mathrm{Tr} \nabla \cdot \bar{\mathbf{r}}_a - \frac{1}{\rho_a} \left[ \frac{3}{2} p_a \nabla \cdot \bar{\mathbf{p}}_a + (\nabla \cdot \bar{\mathbf{p}}_a) \cdot \bar{\mathbf{p}}_a \right] \\ = \frac{1}{2} \mathrm{Tr} \left[ \bar{\mathcal{Q}}_a^{(3)} - \frac{p_a}{\rho_a} (\mathbf{R}_a \bar{\mathcal{I}})^S \right] - \frac{1}{\rho_a} \mathbf{R}_a \cdot \bar{\mathbf{\Pi}}_a^{(2)}, \qquad (C14) $$

where $\mathrm{Tr}(\mathbf{R}_a \bar{\mathcal{I}})^S = 5 \mathbf{R}_a$. This equation is exact. Using the heat flux decomposition (C4) yields

$$ \bar{\mathbf{q}}_a : \nabla \mathbf{u}_a = (2/5) [\mathbf{q}_a \cdot \nabla \mathbf{u}_a + (\nabla \mathbf{u}_a) \cdot \mathbf{q}_a + \mathbf{q}_a \nabla \cdot \mathbf{u}_a], \qquad (C15) $$

and applying a trace at decomposition (C5) yields

$$ \mathrm{Tr} \bar{\mathbf{r}}_a = 5 \frac{p_a^2}{\rho_a} \bar{\mathcal{I}} + \frac{\bar{\mathcal{I}}}{3} \tilde{\mathbf{X}}_a^{(4)} + \bar{\mathbf{\Pi}}_a^{(4)}, \qquad (C16) $$

which is of course equivalent to decomposition (C2), (C3). Note that the closure $\mathrm{Tr} \bar{\mathbf{r}}_a = 5 \frac{p_a^2}{\rho_a} \bar{\mathcal{I}}$ can be viewed as an isotropic analogy of the anisotropic bi-Maxwellian “normal” closure $r_{\parallel a} = \frac{3p_a^2}{\rho}$, $r_{\perp a} = \frac{p_{\parallel a} p_{\perp a}}{\rho_a}$, $r_{\perp a} = \frac{2p_a^2}{\rho_a}$, with $p_{\parallel a} = p_{\perp a} = p_a$, because the following general identity holds for any gyrotropic distribution function: $\mathrm{Tr} \bar{\mathbf{r}}_a^g = r_{\parallel a} \hat{\mathbf{b}} \hat{\mathbf{b}} + r_{\perp a} (\bar{\mathcal{I}} + \hat{\mathbf{b}} \hat{\mathbf{b}}) + 2r_{\perp a} (\bar{\mathcal{I}} - \hat{\mathbf{b}} \hat{\mathbf{b}})$. Then one calculates

$$ \frac{1}{2} \mathrm{Tr} \nabla \cdot \bar{\mathbf{r}}_a = \frac{5}{2} \nabla \left( \frac{p_a^2}{\rho_a} \right) + \frac{1}{6} \nabla \tilde{\mathbf{X}}_a^{(4)} + \frac{1}{2} \nabla \cdot \bar{\mathbf{\Pi}}_a^{(4)}, \qquad (C17) $$

together with

$$ \frac{1}{2} \mathrm{Tr} \nabla \cdot \bar{\mathbf{r}}_a - \frac{1}{\rho_a} \left[ \frac{3}{2} p_a \nabla \cdot \bar{\mathbf{p}}_a + (\nabla \cdot \bar{\mathbf{p}}_a) \cdot \bar{\mathbf{p}}_a \right] = \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) + \frac{1}{6} \nabla \tilde{\mathbf{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{p}}_a) \cdot \bar{\mathbf{\Pi}}_a^{(2)}, \qquad (C18) $$
and evolution Equation (C14) becomes

$$ 
\begin{aligned}
& \frac{d_a \boldsymbol{q}_a}{dt} + \frac{7}{5} \boldsymbol{q}_a \nabla \cdot \boldsymbol{u}_a + \frac{7}{5} \boldsymbol{q}_a \cdot \nabla \boldsymbol{u}_a + \frac{2}{5} (\nabla \boldsymbol{u}_a) \cdot \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) \\
& \qquad + \frac{1}{6} \nabla \tilde{\boldsymbol{X}}_a^{(4)} + \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)} - \frac{1}{\rho_a} (\nabla \cdot \bar{\boldsymbol{\bar{p}}}_a) \cdot \bar{\boldsymbol{\Pi}}_a^{(2)} \\
& = \boldsymbol{Q}_a^{(3)'} \equiv \frac{1}{2} \mathrm{Tr} \bar{\boldsymbol{Q}}_a^{(3)} - \frac{5}{2} \frac{p_a}{\rho_a} \boldsymbol{R}_a - \frac{1}{\rho_a} \boldsymbol{R}_a \cdot \bar{\boldsymbol{\Pi}}_a^{(2)}.
\end{aligned}
\quad (C19) $$

As a double check, reducing the 22-moment model into the 13-moment model with the closures $\tilde{\boldsymbol{X}}_a^{(4)} = 0$ and $\bar{\boldsymbol{\Pi}}_a^{(4)} = 7(p_a/\rho_a)\bar{\boldsymbol{\Pi}}_a^{(2)}$ yields

$$ \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)} \xrightarrow{13-m} \frac{p_a}{\rho_a} \nabla \cdot \bar{\boldsymbol{\Pi}}_a^{(2)} + \frac{7}{2} \bar{\boldsymbol{\Pi}}_a^{(2)} \cdot \nabla \left( \frac{p_a}{\rho_a} \right), \quad (C20) $$

then evolution Equation (C19) recovers Equations (39)–(40) of Schunk (1977).

### C.5. Evolution Equation for the Viscosity Tensor $\bar{\boldsymbol{\Pi}}_a^{(4)}$

The nonlinear evolution equation for the fourth-order moment $r_{ijkl}^a = X_{ijkl}^{a(4)}$ is given by (A8). First, we need to obtain the evolution equation for the matrix $(\mathrm{Tr}\,\bar{\boldsymbol{r}}_a^a)_{ij} = X_{ij}^{a(4)}$, which is further decomposed into (C2) and (C3). Applying a trace at (A8) yields

$$ 
\begin{aligned}
& \frac{d_a}{dt} \mathrm{Tr} \bar{\boldsymbol{r}}_a + \nabla \cdot (\mathrm{Tr} \bar{\boldsymbol{X}}_a^{(5)}) + (\nabla \cdot \boldsymbol{u}_a) \mathrm{Tr} \bar{\boldsymbol{r}}_a + 2 \bar{\boldsymbol{r}}_a: \nabla \boldsymbol{u}_a \\
& \qquad + \left[ (\mathrm{Tr} \bar{\boldsymbol{r}}_a) \cdot \nabla \boldsymbol{u}_a + \Omega_a \hat{\boldsymbol{b}} \times (\mathrm{Tr} \bar{\boldsymbol{r}}_a) - \frac{2}{\rho_a} (\nabla \cdot \bar{\boldsymbol{\bar{p}}}_a) \boldsymbol{q}_a \right]^S - \frac{2}{\rho_a} (\nabla \cdot \bar{\boldsymbol{\bar{p}}}_a) \cdot \bar{\boldsymbol{\bar{q}}}_a \\
& = \mathrm{Tr} \bar{\boldsymbol{Q}}_a^{(4)} - \frac{2}{\rho_a} [(\boldsymbol{R}_a \boldsymbol{q}_a)^S + \boldsymbol{R}_a \cdot \bar{\boldsymbol{\bar{q}}}_a].
\end{aligned}
\quad (C21) $$

As a quick double check, Equation (C21) appears equivalent to Equation (3.4.35) of Balescu (1988, p. 154; after accounting for the different normalization constants of 1/2 and adding a missing “s” index to his fourth-order moment $S_{rsm}$). Applying another trace at (C21) yields

$$ 
\begin{aligned}
& \frac{d_a}{dt} X_a^{(4)} + \nabla \cdot X_a^{(5)} + (\nabla \cdot \boldsymbol{u}_a) X_a^{(4)} + 4 (\mathrm{Tr} \, \bar{\boldsymbol{r}}_a): \nabla \boldsymbol{u}_a - \frac{8}{\rho_a} (\nabla \cdot \bar{\boldsymbol{\bar{p}}}_a) \cdot \boldsymbol{q}_a \\
& = \mathrm{Tr} \, \mathrm{Tr} \, \bar{\boldsymbol{Q}}_a^{(4)} - \frac{8}{\rho_a} \boldsymbol{R}_a \cdot \boldsymbol{q}_a.
\end{aligned}
\quad (C22) $$

To obtain the evolution equation for the matrix $\Pi_{ij}^{a(4)}$, we need to subtract $(\bar{\boldsymbol{\bar{I}}}/3)$ times (C22) from (C21). For example, we need to calculate

$$ 
\begin{aligned}
X_{ijk}^{a(5)} &= \frac{1}{5} (X_i^{a(5)} \delta_{jk} + X_j^{a(5)} \delta_{ik} + X_k^{a(5)} \delta_{ij}); \\
\partial_k X_{kij}^{a(5)} &= \frac{1}{5} (\partial_j X_i^{a(5)} + \partial_i X_j^{a(5)} + \delta_{ij} \partial_k X_k^{a(5)}); \\
(\partial_k X_{kij}^{a(5)}) - \frac{\delta_{ij}}{3} \partial_k X_k^{a(5)} &= \frac{1}{5} \left( \partial_j X_i^{a(5)} + \partial_i X_j^{a(5)} - \frac{2}{3} \delta_{ij} \partial_k X_k^{a(5)} \right),
\end{aligned}
\quad (C23) $$
together with

$$
\begin{aligned}
\bar{\bar{r}}_a: \nabla u_a &= \frac{1}{15} X_a^{(4)} ((\nabla u_a)^S + \bar{\bar{I}}(\nabla \cdot u_a)) \\
&\quad + \frac{1}{7} [\bar{\bar{\Pi}}_a^{(4)}(\nabla \cdot u_a) + \bar{\bar{I}}(\bar{\bar{\Pi}}_a^{(4)}: \nabla u_a) + (\bar{\bar{\Pi}}_a^{(4)} \cdot \nabla u_a)^S + ((\nabla u_a) \cdot \bar{\bar{\Pi}}_a^{(4)})^S];
\end{aligned}
$$

$$
\begin{aligned}
\operatorname{Tr}\bar{\bar{r}}_a: \nabla u_a &= \frac{1}{3} X_a^{(4)} (\nabla \cdot u_a) + \bar{\bar{\Pi}}_a^{(4)}: \nabla u_a; \\
2\bar{\bar{r}}_a: \nabla u_a - \frac{\bar{\bar{I}}}{3} 4 \operatorname{Tr}\bar{\bar{r}}_a: \nabla u_a &= \frac{2}{15} X_a^{(4)} \left( (\nabla u_a)^S - \frac{7}{3} \bar{\bar{I}} (\nabla \cdot u_a) \right) - \frac{22}{21} \bar{\bar{I}} (\bar{\bar{\Pi}}_a^{(4)}: \nabla u_a) \\
&\quad + \frac{2}{7} [\bar{\bar{\Pi}}_a^{(4)} (\nabla \cdot u_a) + (\bar{\bar{\Pi}}_a^{(4)} \cdot \nabla u_a)^S + ((\nabla u_a) \cdot \bar{\bar{\Pi}}_a^{(4)})^S],
\end{aligned}
\qquad \text{(C24)}
$$

and the useful identities are

$$
\begin{aligned}
&[(\operatorname{Tr}\bar{\bar{r}}_a) \cdot \nabla u_a]^S = \frac{1}{3} X_a^{(4)} (\nabla u_a)^S + [\bar{\bar{\Pi}}_a^{(4)} \cdot \nabla u]^S; \\
&[\hat{\mathbf{b}} \times (\operatorname{Tr}\bar{\bar{r}}_a)]^S = [\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(4)}]^S.
\end{aligned}
\qquad \text{(C25)}
$$

The heat flux contributions calculate

$$
\begin{aligned}
& (\nabla \cdot \bar{\bar{p}}_a) \cdot \bar{\bar{q}}_a = \frac{2}{5} [((\nabla \cdot \bar{\bar{p}}_a) \mathbf{q}_a)^S + \bar{\bar{I}}(\nabla \cdot \bar{\bar{p}}_a) \cdot \mathbf{q}_a], \\
& \text{so the heat fluxes are added as} \\
& -2((\nabla \cdot \bar{\bar{p}}_a) \mathbf{q}_a)^S - 2(\nabla \cdot \bar{\bar{p}}_a) \cdot \bar{\bar{q}}_a + \bar{\bar{I}} \frac{8}{3} (\nabla \cdot \bar{\bar{p}}_a) \cdot \mathbf{q}_a \\
& = -\frac{14}{5} \left[ ((\nabla \cdot \bar{\bar{p}}_a) \mathbf{q}_a)^S - \frac{2}{3} \bar{\bar{I}} (\nabla \cdot \bar{\bar{p}}_a) \cdot \mathbf{q}_a \right].
\end{aligned}
\qquad \text{(C26)}
$$

The fully nonlinear evolution equation for the matrix $\bar{\bar{\Pi}}_a^{(4)}$ thus reads

$$
\begin{aligned}
& \frac{d_a}{dt} \bar{\bar{\Pi}}_a^{(4)} + \frac{1}{5} \left[ (\nabla X_a^{(5)})^S - \frac{2}{3} \bar{\bar{I}} (\nabla \cdot X_a^{(5)}) \right] + \frac{9}{7} (\nabla \cdot u_a) \bar{\bar{\Pi}}_a^{(4)} + \frac{9}{7} (\bar{\bar{\Pi}}_a^{(4)} \cdot \nabla u_a)^S \\
& + \frac{2}{7} ((\nabla u_a) \cdot \bar{\bar{\Pi}}_a^{(4)})^S - \frac{22}{21} \bar{\bar{I}} (\bar{\bar{\Pi}}_a^{(4)}: \nabla u_a) - \frac{14}{5 \rho_a} \left[ ((\nabla \cdot \bar{\bar{p}}_a) \mathbf{q}_a)^S - \frac{2}{3} \bar{\bar{I}} (\nabla \cdot \bar{\bar{p}}_a) \cdot \mathbf{q}_a \right] \\
& + \Omega_a (\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(4)})^S + \frac{7}{15} \left( 15 \frac{p_a^2}{\rho_a} + \tilde{X}_a^{(4)} \right) \left[ (\nabla u_a)^S - \frac{2}{3} \bar{\bar{I}} (\nabla \cdot u_a) \right] \\
& = \bar{\bar{Q}}_a^{(4)'} \equiv \operatorname{Tr} \bar{\bar{Q}}_a^{(4)} - \frac{\bar{\bar{I}}}{3} \operatorname{Tr} \operatorname{Tr} \bar{\bar{Q}}_a^{(4)} - \frac{14}{5 \rho_a} \left[ (\mathbf{R}_a \mathbf{q}_a)^S - \frac{2}{3} \bar{\bar{I}} (\mathbf{R}_a \cdot \mathbf{q}_a) \right].
\end{aligned}
\qquad \text{(C27)}
$$

At the semilinear level (while keeping the $d/dt$), evolution Equation (C27) simplifies into

$$
\begin{aligned}
& \frac{d_a}{dt} \bar{\bar{\Pi}}_a^{(4)} + \frac{1}{5} \left[ (\nabla X_a^{(5)})^S - \frac{2}{3} \bar{\bar{I}} (\nabla \cdot X_a^{(5)}) \right] + \Omega_a (\hat{\mathbf{b}} \times \bar{\bar{\Pi}}_a^{(4)})^S \\
& + 7 \frac{p_a^2}{\rho_a} \left[ (\nabla u_a)^S - \frac{2}{3} \bar{\bar{I}} (\nabla \cdot u_a) \right] = \bar{\bar{Q}}_a^{(4)'} = \operatorname{Tr} \bar{\bar{Q}}_a^{(4)} - \frac{\bar{\bar{I}}}{3} \operatorname{Tr} \operatorname{Tr} \bar{\bar{Q}}_a^{(4)}.
\end{aligned}
\qquad \text{(C28)}
$$

Finally, neglecting the coupling between heat fluxes and viscosities (which is the choice of Braginskii), the simplest evolution equation reads

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

where $\bar{\bar{W}}_a = (\nabla u_a)^S - (2/3)\bar{\bar{I}}(\nabla \cdot u_a)$ is the usual rate-of-strain tensor.
### C.6. Evolution Equation for the Perturbation $\tilde{X}_a^{(4)}$

The fully nonlinear evolution Equation (C22) for $X_a^{(4)}$ reads

$$ \begin{aligned} & \frac{d_a X_a^{(4)}}{dt} + \nabla \cdot X_a^{(5)} + \frac{7}{3} X_a^{(4)} (\nabla \cdot \mathbf{u}_a) + 4 \bar{\mathbf{\Pi}}_a^{(4)} : \nabla \mathbf{u}_a - \frac{8}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot \mathbf{q}_a \\ & = \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{Q}}_a^{(4)} - \frac{8}{\rho_a} \mathbf{R}_a \cdot \mathbf{q}_a. \end{aligned} \quad (C30) $$

Then, by using $X_a^{(4)} = 15(p_a^2/\rho_a) + \tilde{X}_a^{(4)}$ with

$$ \frac{d_a}{dt} \left( \frac{p_a^2}{\rho_a} \right) = \frac{p_a}{\rho_a} \left[ -\frac{7}{3} p_a \nabla \cdot \mathbf{u}_a - \frac{4}{3} \nabla \cdot \mathbf{q}_a - \frac{4}{3} \bar{\mathbf{\Pi}}_a^{(2)} : \nabla \mathbf{u}_a + \frac{4}{3} Q_a \right], \quad (C31) $$

one obtains the fully nonlinear evolution equation for $\tilde{X}_a^{(4)}$:

$$ \begin{aligned} & \frac{d_a}{dt} \tilde{X}_a^{(4)} + \nabla \cdot 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 \\ & \quad - \frac{8}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_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} Q_a - \frac{8}{\rho_a} \mathbf{R}_a \cdot \mathbf{q}_a, \end{aligned} \quad (C32) $$

and, at the semilinear level,

$$ \frac{d_a}{dt} \tilde{X}_a^{(4)} + \nabla \cdot X_a^{(5)} - 20 \frac{p_a}{\rho_a} \nabla \cdot \mathbf{q}_a = \tilde{Q}_a^{(4) \prime} = \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{Q}}_a^{(4)} - 20 \frac{p_a}{\rho_a} Q_a. \quad (C33) $$

The collisional contributions can be found in Section 7.1; see Equation (142).

### C.7. Evolution Equation for the Heat Flux Vector $X_a^{(5)}$

Applying a trace twice at (A9) yields

$$ \begin{aligned} & \frac{\partial}{\partial t} \mathrm{Tr} \bar{\mathbf{X}}^{a(5)} + \nabla \cdot (\mathrm{Tr} \bar{\mathbf{X}}^{a(6)}) + \nabla \cdot (\mathbf{u}^a \mathrm{Tr} \bar{\mathbf{X}}^{a(5)}) + (\mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{a(5)} \cdot \nabla) \mathbf{u}^a \\ & \quad + 4(\mathrm{Tr} \bar{\mathbf{X}}^{a(5)}): \nabla \mathbf{u}^a + \Omega_a \hat{\mathbf{b}} \times (\mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{a(5)}) - \frac{1}{\rho_a} [(\nabla \cdot \bar{\mathbf{p}}^a) \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{a(4)} + 4(\nabla \cdot \bar{\mathbf{p}}^a) \cdot \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{a(4)}] \\ & = \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{Q}}^{a(5)} - \frac{1}{\rho_a} [\mathbf{R}^a \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{a(4)} + 4 \mathbf{R}^a \cdot \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{a(4)}]. \end{aligned} \quad (C34) $$

By using the definitions of the vectors $X^{(5)} = \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{X}}^{(5)}$, $Q^{(5)} = \mathrm{Tr} \mathrm{Tr} \bar{\mathbf{Q}}^{(5)}$, and

$$ X_{ijk}^{(5)} = \frac{1}{5} [X_i^{(5)} \delta_{jk} + X_j^{(5)} \delta_{ik} + X_k^{(5)} \delta_{ij}]; $$

$$ X_{ijk}^{(5)} \partial_j u_k = \frac{1}{5} [X_i^{(5)} \nabla \cdot \mathbf{u}_a + X_j^{(5)} \partial_j u_i^a + X_k^{(5)} \partial_i u_k^a], \quad (C35) $$

together with the decompositions (C2) and (C3), the fully nonlinear evolution equation becomes

$$ \begin{aligned} & \frac{d_a}{dt} X_a^{(5)} + \frac{1}{3} \nabla \tilde{X}_a^{(6)} + \nabla \cdot \bar{\mathbf{\Pi}}_a^{(6)} + \frac{9}{5} X_a^{(5)} (\nabla \cdot \mathbf{u}_a) + \frac{9}{5} X_a^{(5)} \nabla \mathbf{u}_a + \frac{4}{5} (\nabla \mathbf{u}_a) \cdot X_a^{(5)} + \Omega_a \hat{\mathbf{b}} \times X_a^{(5)} \\ & \quad + 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{p}}^a) \tilde{X}_a^{(4)} - \frac{4}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}^a) \cdot \bar{\mathbf{\Pi}}_a^{(4)} \\ & = Q_a^{(5) \prime} \equiv 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} \quad (C36) $$
Because we do not go higher in the hierarchy, the model is closed with the closures (see Equations (B130) and (B131) or Section 8.6 with Tables 3 and 4)

$$ \tilde{X}_a^{(6)} = 21 \frac{P_a}{\rho_a} \tilde{X}_a^{(4)}; \quad \bar{\Pi}_a^{(6)} = 18 \frac{P_a}{\rho_a} \bar{\Pi}_a^{(4)} - 63 \frac{P_a^2}{\rho_a^2} \bar{\Pi}_a^{(2)}. \qquad (C37) $$

At a semilinear level, Equation (C36) becomes

$$ \frac{d_a}{dt} \mathbf{X}_a^{(5)} + 7 \frac{P_a}{\rho_a} \nabla \tilde{X}_a^{(4)} + 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)} \\ + \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} = \mathbf{Q}_a^{(5)} - 35 \frac{P_a^2}{\rho_a^2} \mathbf{R}_a. \qquad (C38) $$

## Appendix D
### Simplified General Fluid Hierarchy

Previously, we introduced a full fluid hierarchy in Section A, which contains n-dimensional moments $X_{ijk...n}^{(n)}$. By applying contractions at these moments in Appendix C, we derived evolution equations for the 22-moment model. Instead of doing that, it is of course possible to obtain evolution equations for contracted moments directly from the Boltzmann equation. This simplified hierarchy is formulated with heat fluxes (vectors) and stress tensors (matrices)

$$ \mathbf{X}_a^{(2n+1)} = m_a \int \mathbf{c}_a |\mathbf{c}_a|^{2n} f_a d^3v; \quad \bar{\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, \qquad (D1) $$

together with fully contracted scalars, which are decomposed into a Maxwellian core and perturbation (notation with tilde)

$$ 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{X}_a^{(2n)}, \qquad (D2) $$

meaning a definition of $\tilde{X}_a^{(2n)} = m_a \int |\mathbf{c}_a|^{2n} (f_a - f_a^{(0)}) d^3v$, where $f_a^{(0)}$ is Maxwellian. In another words, one considers the matrices

$$ X_{ij}^{a(2n)} = m_a \int |\mathbf{c}_a|^{2n-2} c_i^a c_j^a f_a d^3v = \frac{\delta_{ij}}{3} X^{a(2n)} + \Pi_{ij}^{a(2n)}, \qquad (D3) $$

which are decomposed into fully contracted scalars and stress tensors. Note that $X_a^{(1)} = 0$ and $\tilde{X}_a^{(2)} = 0$.

Unfortunately, the traditional definition of the heat flux vector $\mathbf{q}_a = (1/2)\text{Tr}\bar{\mathbf{q}}_a$, which contains a factor of 1/2, goes against the general ideology that no additional factors are introduced by contractions. Also, we have previously reserved vector $\mathbf{Q}_a^{(3)\prime}$ for the right-hand side of the heat flux $\mathbf{q}_a$ evolution equation, and not for $\tilde{X}_a^{(3)}$. Obviously, our previous notation is not ideal for generalization to nth-order moments. To circumvent all of the problems with the previous definitions, we define new collisional contributions for the heat fluxes and stress tensors with $\mathcal{Q}$ (mathcal{Q}), as the vectors and matrices

$$ \begin{aligned} Q_i^{a(2n+1)} &= m_a \int |\mathbf{c}_a|^{2n} c_i^a C(f_a) d^3v; \\ Q_{ij}^{a(2n)} &= m_a \int |\mathbf{c}_a|^{2n-2} c_i^a c_j^a C(f_a) d^3v; \end{aligned} \qquad (D4) $$

together with the fully contracted

$$ Q_a^{(2n)} = m_a \int |\mathbf{c}_a|^{2n} C(f_a) d^3v; \quad Q_a = \frac{m_a}{2} \int |\mathbf{c}_a|^2 C(f_a) d^3v. \qquad (D5) $$

The energy exchange rates $Q_a$ contain the traditional factor of 1/2, and $Q_a^{(2)} = 2Q_a$. The momentum exchange rates $\mathbf{R}_a = m_a \int \mathbf{v} C(f_a) d^3v$. In the vector notation, matrix $\bar{\mathbf{Q}}_a^{(2n)} = \text{TrTr} \dots \text{Tr} \bar{\mathbf{Q}}_a^{(2n)}$.
Then, the direct integration of the Boltzmann equation and the subtraction of the momentum equations yields the evolution equations for the scalars

$$ \frac{\partial}{\partial t} X^{a(2n)} + \partial_k (u_k^a X^{a(2n)}) + \partial_k X_k^{a(2n+1)} + (2n) X_{ik}^{a(2n)} \partial_k u_i^a \\ - \frac{(2n)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}^a)_k X_k^{a(2n-1)} = Q^{a(2n)} - \frac{(2n)}{\rho_a} R_k^a X_k^{a(2n-1)}, \qquad (D6) $$

where (n) without a species index should not be confused with the number density, the evolution equations for the vectors

$$ \frac{\partial}{\partial t} X_i^{a(2n+1)} + \partial_k (u_k^a X_i^{a(2n+1)}) + \partial_k X_{ki}^{a(2n+2)} + X_k^{a(2n+1)} \partial_k u_i^a + (2n) X_{ijk}^{a(2n+1)} \partial_k u_j^a \\ - \frac{(2n)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}^a)_k X_{ki}^{a(2n)} - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}^a)_i X^{a(2n)} + \Omega_a (\hat{\mathbf{b}} \times \mathbf{X}^{a(2n+1)})_i \\ = Q_i^{a(2n+1)} - \frac{1}{\rho_a} R_i^a X^{a(2n)} - \frac{(2n)}{\rho_a} R_k^a X_{ki}^{a(2n)}, \qquad (D7) $$

and the matrices

$$ \frac{\partial}{\partial t} X_{ij}^{a(2n)} + \partial_k (u_k^a X_{ij}^{a(2n)}) + \partial_k X_{kij}^{a(2n+1)} + (2n-2) X_{ijkl}^{a(2n)} (\partial_k u_l^a) \\ + \left[ X_{ik}^{a(2n)} \partial_k u_j^a + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{X}}^{a(2n)})_{ij} - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}^a)_i X_j^{a(2n-1)} \right]^S - \frac{(2n-2)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}^a)_k X_{kij}^{a(2n-1)} \\ = Q_{ij}^{a(2n)} - \frac{1}{\rho_a} [R_i^a X_j^{a(2n-1)}]^S - \frac{(2n-2)}{\rho_a} R_k^a X_{kij}^{a(2n-1)}, \qquad (D8) $$

which are valid for $n \ge 1$. For example, evaluating (D6) for $n = 1$ yields the evolution equation for the scalar pressure $p_a$. Applying a trace at (D8) recovers (D6).

The matrices $X_{ij}^{a(2n)}$ are then decomposed according to (D3), where the stress tensors $\Pi_{ij}^{a(2n)}$ are traceless, and higher-order tensors are decomposed according to (where the tensors $\sigma$ are neglected, which is the core of the hierarchy simplification)

$$ X_{ijk}^{a(2n+1)} = \frac{1}{5} [X_{i}^{a(2n+1)} \delta_{jk} + X_{j}^{a(2n+1)} \delta_{ik} + X_{k}^{a(2n+1)} \delta_{ij}], \qquad (D9) $$

$$ X_{ijkl}^{a(2n)} = \frac{1}{15} X^{a(2n)} (\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}) \\ + \frac{1}{7} [\Pi_{ij}^{a(2n)} \delta_{kl} + \Pi_{ik}^{a(2n)} \delta_{jl} + \Pi_{il}^{a(2n)} \delta_{jk} + \Pi_{jk}^{a(2n)} \delta_{il} + \Pi_{jl}^{a(2n)} \delta_{ik} + \Pi_{kl}^{a(2n)} \delta_{ij}]. \qquad (D10) $$

Applying a trace at (D9) yields the identity, and applying a trace at (D10) yields the decomposition (D3). The evolution equations for the fully contracted moments (scalsors) then become

$$ \frac{d_a}{dt} X_a^{(2n)} + \nabla \cdot X_a^{(2n+1)} + \frac{(2n+3)}{3} X_a^{(2n)} \nabla \cdot u_a + (2n) \bar{\Pi}_a^{(2n)} : \nabla u_a \\ - \frac{(2n)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot X_a^{(2n-1)} = Q_a^{(2n)} - \frac{(2n)}{\rho_a} \mathbf{R}_a \cdot X_a^{(2n-1)}, \qquad (D11) $$

for the heat fluxes (vectors)

$$ \frac{d_a}{dt} X_a^{(2n+1)} + \frac{(2n+5)}{5} [X_a^{(2n+1)} \nabla \cdot \mathbf{u}_a + X_a^{(2n+1)} \cdot \nabla \mathbf{u}_a] + \frac{(2n)}{5} (\nabla \mathbf{u}_a) \cdot X_a^{(2n+1)} \\ + \frac{1}{3} \nabla X_a^{(2n+2)} + \nabla \cdot \bar{\Pi}_a^{(2n+2)} - \frac{(2n+3)}{3\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) X_a^{(2n)} - \frac{(2n)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot \bar{\Pi}_a^{(2n)} \\ + \Omega_a \hat{\mathbf{b}} \times X_a^{(2n+1)} = \bar{Q}_a^{(2n+1)} - \frac{(2n+3)}{3\rho_a} \mathbf{R}_a X_a^{(2n)} - \frac{(2n)}{\rho_a} \mathbf{R}_a \cdot \bar{\Pi}_a^{(2n)}, \qquad (D12) $$
and for the stress tensors (matrices)

$$
\begin{aligned}
& \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] + \frac{(2n+5)}{7} \bar{\Pi}_a^{(2n)} (\nabla \cdot \mathbf{u}_a) \\
& \qquad + \left[ \frac{(2n+5)}{7} (\bar{\Pi}_a^{(2n)}) \cdot \nabla \mathbf{u}_a \right]^S + \frac{(2n-2)}{7} ((\nabla \mathbf{u}_a) \cdot \bar{\Pi}_a^{(2n)})^S - \frac{2(4n+3)}{21} \bar{I} (\bar{\Pi}_a^{(2n)} : \nabla \mathbf{u}_a) \\
& \qquad - \frac{(2n+3)}{5\rho_a} \left[ ((\nabla \cdot \bar{\mathbf{p}}_a) X_a^{(2n-1)})^S - \frac{2}{3} \bar{I} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot X_a^{(2n-1)} \right] \\
& \qquad + \Omega_a (\hat{\mathbf{b}} \times \bar{\Pi}_a^{(2n)})^S + \frac{(2n+3)}{15} X_a^{(2n)} \bar{W}_a \\
& = \bar{Q}_a^{(2n)} \quad \equiv \bar{Q}_a^{(2n)} - \frac{\bar{I}}{3} Q_a^{(2n)} - \frac{(2n+3)}{5\rho_a} \left[ (\mathbf{R}_a X_a^{(2n-1)})^S - \frac{2}{3} \bar{I} \mathbf{R}_a \cdot X_a^{(2n-1)} \right].
\end{aligned}
\qquad (D13)
$$

By applying a trace at Equation (D13), it can be verified that it is traceless.

The fully contracted scalar variables are then decomposed into a Maxwellian core and perturbation (with tilde) according to (D2), yielding the evolution equation for scalars

$$
\begin{aligned}
& \frac{d_a \tilde{X}_a^{(2n)}}{dt} + \nabla \cdot X_a^{(2n+1)} + \frac{(2n+3)}{3} \tilde{X}_a^{(2n)} \nabla \cdot \mathbf{u}_a + (2n) \bar{\Pi}_a^{(2n)} : \nabla \mathbf{u}_a \\
& \qquad - (2n+1)!! \frac{(2n)}{3} \left(\frac{p_a}{\rho_a}\right)^{n-1} \left[\nabla \cdot \mathbf{q}_a + \bar{\Pi}_a^{(2)} : \nabla \mathbf{u}_a\right] - \frac{(2n)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot X_a^{(2n-1)} \\
& = \tilde{Q}_a^{(2n)} \quad \equiv Q_a^{(2n)} - (2n+1)!! \frac{(2n)}{3} \left(\frac{p_a}{\rho_a}\right)^{n-1} Q_a - \frac{(2n)}{\rho_a} \mathbf{R}_a \cdot X_a^{(2n-1)},
\end{aligned}
\qquad (D14)
$$

and for heat fluxes

$$
\begin{aligned}
& \frac{d_a X_a^{(2n+1)}}{dt} + \frac{(2n+5)}{5} [X_a^{(2n+1)} \nabla \cdot \mathbf{u}_a + X_a^{(2n+1)} \cdot \nabla \mathbf{u}_a] + \frac{(2n)}{5} (\nabla \mathbf{u}_a) \cdot X_a^{(2n+1)} \\
& \qquad + \frac{1}{3} \nabla \tilde{X}_a^{(2n+2)} + \nabla \cdot \bar{\Pi}_a^{(2n+2)} - \frac{(2n+3)}{3\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \tilde{X}_a^{(2n)} - \frac{(2n)}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot \bar{\Pi}_a^{(2n)} \\
& \qquad + (2n+3)!! \frac{(n)}{3} \frac{p_a^n}{\rho_a^{n-1}} \nabla \left( \frac{p_a}{\rho_a} \right) - \frac{(2n+3)!!}{3} \frac{p_a^n}{\rho_a^n} \nabla \cdot \bar{\Pi}_a^{(2)} + \Omega_a \hat{\mathbf{b}} \times X_a^{(2n+1)} \\
& = \tilde{Q}_a^{(2n+1)} \quad \equiv \tilde{Q}_a^{(2n+1)} - \frac{(2n+3)}{3\rho_a} \mathbf{R}_a \tilde{X}_a^{(2n)} - \frac{(2n+3)!!}{3} \frac{p_a^n}{\rho_a^n} \mathbf{R}_a - \frac{(2n)}{\rho_a} \mathbf{R}_a \cdot \bar{\Pi}_a^{(2n)}.
\end{aligned}
\qquad (D15)
$$

The evolution equation for the stress tensors (D13) contains only one trivial term with $X_a^{(2n)}$, where

$$
\frac{(2n + 3)}{15} X_a^{(2n)} \bar{W}_a = \frac{(2n + 3)!!}{15} \frac{p_a^n}{\rho_a^{n-1}} \bar{W}_a + \frac{(2n + 3)}{15} \tilde{X}_a^{(2n)} \bar{W}_a,
$$

and we do not rewrite the full equation. Equations (D13)–(D15) are valid for $n \ge 1$, where for $n = 1$ (D14) reduces to zero, so this equation is meaningful only for $n \ge 2$. In the semilinear approximation, the hierarchy simplifies into (189)–(191).

## Appendix E
### BGK Collisional Operator

Before the calculations with the Landau collisional operator, it is beneficial to first become familiar with the heuristic relaxation-type operator known as BGK, after Bhatnagar-Gross-Krook (Bhatnagar et al. 1954; Gross & Krook 1956), written in the following form:

$$C(f_a) = \sum_b C_{ab}(f_a) = - \sum_b \nu_{ab}(f_a - f_{ab}^{(0)}). \qquad (E1)$$

The Maxwellian $f_{ab}^{(0)}$ has two indices and is defined as

$$f_{ab}^{(0)} = n_a \left( \frac{m_a}{2\pi T_a} \right)^{3/2} \exp \left( - \frac{m_a |\mathbf{v} - \mathbf{u}_b|^2}{2T_a} \right). \qquad (E2)$$
Note that only the velocity $\mathbf{u}_b$ has the index “b” and that the temperature, mass, and density have the index “a.” Accounting for different temperatures is possible by considering the generalized BGK operators of Haack et al. (2017). The simple BGK operator yields momentum and energy exchange rates

$$ \mathbf{R}_{ab} = \rho_a \nu_{ab} (\mathbf{u}_b - \mathbf{u}_a); \quad Q_{ab} = \frac{1}{2} \rho_a \nu_{ab} |\mathbf{u}_b - \mathbf{u}_a|^2, \qquad (E3) $$

where both the momentum and energy are conserved (note that for heuristic operators it is advisable to directly calculate both $\mathbf{R}_{ab}$ and $\mathbf{R}_{ba}$ together with $Q_{ab}$ and $Q_{ba}$ to verify that they are well-defined). This BGK operator also satisfies the Boltzmann H-theorem, which for multispecies plasmas has the general form

$$ \int C_{ab}(f_a) \ln f_a d^3v + \int C_{ba}(f_b) \ln f_b d^3v \le 0, \qquad (E4) $$

where the equality is true only if $f_a$ and $f_b$ are Maxwellians. For the BGK operator, each part of the H-theorem (E4) is satisfied independently. It can be shown that $\int (f_a - f_{ab}^{(0)}) \ln f_{ab}^{(0)} d^3v = 0$, and subtracting this integral from the first term of (E4) yields

$$ \begin{aligned} \int C_{ab}(f_a) \ln f_a d^3v &= \nu_{ab} \int (f_{ab}^{(0)} - f_a) \ln f_a d^3v - \underbrace{\nu_{ab} \int (f_{ab}^{(0)} - f_a) \ln f_{ab}^{(0)} d^3v}_{0} \\ &= \nu_{ab} \int (f_{ab}^{(0)} - f_a) \ln \left( \frac{f_a}{f_{ab}^{(0)}} \right) d^3v \le 0, \end{aligned} \qquad (E5) $$

where in the last step one uses that for any real numbers $a > 0$ and $b > 0$, the following identity holds: $(a - b)\ln(b/a) \le 0$ (the identity is easily verified, because for $a > b$ the first term is positive and the logarithm is negative, and for $a < b$ the first term is negative and the logarithm is positive; the identity is equal to zero only if $a = b$).

The BGK collisional contributions calculate

$$ \bar{\mathbf{Q}}_{ab}^{(2)} = m_a \int \mathbf{c}_a \mathbf{c}_a C_{ab}(f_a) d^3v = - \nu_{ab} \bar{\mathbf{\Pi}}_a^{(2)} + \nu_{ab} \rho_a \delta \mathbf{u} \delta \mathbf{u}; \qquad (E6) $$

$$ \bar{\mathbf{Q}}_{ab}^{(3)} = m_a \int \mathbf{c}_a \mathbf{c}_a C_{ab}(f_a) d^3v = - \nu_{ab} \bar{\mathbf{q}}_a + \nu_{ab} p_a [\delta \mathbf{u} \bar{\mathbf{I}}]^S + \nu_{ab} \rho_a \delta \mathbf{u} \delta \mathbf{u} \delta \mathbf{u}, \qquad (E7) $$

where $\delta \mathbf{u} = \mathbf{u}_b - \mathbf{u}_a$.

### E.1. Viscosity Tensor $\bar{\mathbf{\Pi}}_a^{(2)}$

The collisional contributions that enter the right-hand side of evolution Equation (C12) are

$$ \bar{\mathbf{Q}}_a^{(2)} \equiv \bar{\mathbf{Q}}_a^{(2)} - \frac{\bar{\mathbf{I}}}{3} \mathrm{Tr} \bar{\mathbf{Q}}_a^{(2)} = - \bar{\nu}_a \bar{\mathbf{\Pi}}_a^{(2)} - \bar{\mathbf{W}}_a^{\mathrm{frict}}, \qquad (E8) $$

where we define

$$ \bar{\nu}_a = \sum_b \nu_{ab}; \qquad (E9) $$

$$ \bar{\mathbf{W}}_a^{\mathrm{frict}} = - \rho_a \sum_b \nu_{ab} \left( \delta \mathbf{u} \delta \mathbf{u} - \frac{\bar{\mathbf{I}}}{3} |\delta \mathbf{u}|^2 \right), \qquad (E10) $$

and where the superscript “frict” means the frictional contributions due to $\delta \mathbf{u}$. The frictional contributions are only nonlinear, but we keep them to show that it is possible to take them into account. Using the quasistatic approximation, evolution Equation (C12) can be simplified into

$$ (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(2)})^S + \frac{\bar{\nu}_a}{\Omega_a} \bar{\mathbf{\Pi}}_a^{(2)} = - \frac{1}{\Omega_a} (p_a \bar{\mathbf{W}}_a + \bar{\mathbf{W}}_a^q + \bar{\mathbf{W}}_a^{\mathrm{frict}}), \qquad (E11) $$

where the matrices $\bar{\mathbf{W}}_a$ and $\bar{\mathbf{W}}_a^q$ are given by (C11), (C13). Equation (E11) can be directly solved. Nevertheless, the stress tensor of Braginskii does not contain heat flux contributions or frictional contributions. To understand the solution of Braginskii more clearly, let us first solve the above equation only with the matrix $\bar{\mathbf{W}}_a$.

The simplest quasistatic $\bar{\mathbf{\Pi}}_a^{(2)}$ is thus obtained by solving

$$ (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(2)})^S + \frac{\bar{\nu}_a}{\Omega_a} \bar{\mathbf{\Pi}}_a^{(2)} = - \frac{p_a}{\Omega_a} \bar{\mathbf{W}}_a. \qquad (E12) $$
For any traceless and symmetric matrix $\bar{\bar{W}}_a$, the solution of (E12) reads (see the details in Appendix E.4)

$$
\begin{align}
\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; \nonumber \\
\bar{\bar{W}}_0 &= \frac{3}{2} (\bar{\bar{W}}_a \cdot \hat{b}\hat{b}) \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right); \nonumber \\
\bar{\bar{W}}_1 &= \bar{\bar{I}}_{\perp} \cdot \bar{\bar{W}}_a \cdot \bar{\bar{I}}_{\perp} + \frac{1}{2} (\bar{\bar{W}}_a \cdot \hat{b}\hat{b}) \bar{\bar{I}}_{\perp}; \nonumber \\
\bar{\bar{W}}_2 &= (\bar{\bar{I}}_{\perp} \cdot \bar{\bar{W}}_a \cdot \hat{b}\hat{b})^S; \nonumber \\
\bar{\bar{W}}_3 &= \frac{1}{2} (\hat{b} \times \bar{\bar{W}}_a \cdot \bar{\bar{I}}_{\perp})^S; \nonumber \\
\bar{\bar{W}}_4 &= (\hat{b} \times \bar{\bar{W}}_a \cdot \hat{b}\hat{b})^S, \tag{E13}
\end{align}
$$

with BGK viscosity coefficients

$$
\eta_0^a = \frac{p_a}{\bar{\nu}_a}; \quad \eta_1^a = \frac{p_a \bar{\nu}_a}{4\Omega_a^2 + \bar{\nu}_a^2}; \quad \eta_2^a = \frac{p_a \bar{\nu}_a}{\Omega_a^2 + \bar{\nu}_a^2}; \quad \eta_3^a = \frac{2p_a \Omega_a}{4\Omega_a^2 + \bar{\nu}_a^2}; \quad \eta_4^a = \frac{p_a \Omega_a}{\Omega_a^2 + \bar{\nu}_a^2}. \tag{E14}
$$

The coefficient $\eta_0$ is called the parallel viscosity, $\eta_1$, $\eta_2$ are perpendicular viscosities, and $\eta_3$, $\eta_4$ are gyroviscosities. Importantly, the BGK solution (E13) is identical to the form of the Braginskii (1965) viscosity tensor, his Equations (4.41)–(4.42), only his viscosities are different. A comparison is presented in the next section. All four matrices $\bar{\bar{W}}_0, \dots, \bar{\bar{W}}_4$ are traceless and $\bar{\bar{W}}_0 + \bar{\bar{W}}_1 + \bar{\bar{W}}_2 = \bar{\bar{W}}_a$.

When the magnetic field is zero, so that $\Omega_a = 0$ and $\eta_0^a = \eta_1^a = \eta_2^a$, the stress tensor (E13) simplifies into $\bar{\bar{\Pi}}_a^{(2)} = - \eta_0^a \bar{\bar{W}}_a$ and contributes to the momentum equations in a familiar form:

$$
\mathbf{B} = 0: \quad \nabla \cdot \bar{\bar{\Pi}}_a^{(2)} = - \nabla \cdot (\eta_0^a \bar{\bar{W}}_a) = - \eta_0^a \left( \nabla^2 \mathbf{u}_a + \frac{1}{3} \nabla (\nabla \cdot \mathbf{u}_a) \right) - (\nabla \eta_0^a) \cdot \bar{\bar{W}}_a, \quad (\text{E15})
$$

analogously to the viscosity of the Navier–Stokes equations (the last term can be neglected if $\eta_0^a$ is spatially independent).

If the mean magnetic field is sufficiently strong that its curvature can be neglected, (E13) can be evaluated with respect to $\hat{b}_0 = (0, 0, 1)$, yielding

$$
\begin{align}
\Pi_{xx}^{(2)a} &= -\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; \nonumber \\
\Pi_{xy}^{(2)a} &= \frac{\eta_3^a}{2}(W_{xx}^a - W_{yy}^a) - \eta_1^a W_{xy}^a; \nonumber \\
\Pi_{xz}^{(2)a} &= -\eta_4^a W_{yz}^a - \eta_2^a W_{xz}^a; \nonumber \\
\Pi_{yy}^{(2)a} &= -\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; \nonumber \\
\Pi_{yz}^{(2)a} &= \eta_4^a W_{xz}^a - \eta_2^a W_{yz}^a; \nonumber \\
\Pi_{zz}^{(2)a} &= -\eta_0^a W_{zz}^a, \tag{E16}
\end{align}
$$

which is Equation (2.21) of Braginskii (1965). As a double check, adding $\Pi_{xx}^{(2)a} + \Pi_{yy}^{(2)a} + \Pi_{zz}^{(2)a} = - \eta_0^a(W_{xx}^a + W_{yy}^a + W_{zz}^a) = 0$, so the stress tensor is indeed traceless (even though all of the diagonal components are nonzero). For strong magnetic field $\Omega_a \gg \bar{\nu}_a$, the viscosities (E14) simplify into

$$
\eta_0^a = \frac{p_a}{\bar{\nu}_a}; \quad \eta_1^a = \frac{1}{4} \frac{p_a \bar{\nu}_a}{\Omega_a^2}; \quad \eta_2^a = 4\eta_1^a; \quad \eta_3^a = \frac{p_a}{2\Omega_a}; \quad \eta_4^a = 2\eta_3^a. \tag{E17}
$$

Considering only self-collisions, the BGK viscosity coefficients (E14) were first recovered by Kaufman (1960), even though he does not write them explicitly, and one needs to get them by rearranging his Equations (12)–(15) into form (E16). The same results for $\eta_0 - \eta_3$ can also be found in Helander & Sigmar (2002, p. 86), for example; see also Zank (2014, p. 164), although the $\eta_4$ coefficient is erroneously related to $\eta_3 = 2\eta_4$, which is a valid relation only in the limit when $x = \Omega_a/\bar{\nu}_a$ is small (i.e., a weak magnetic field). The correct relations are $\eta_3^a(x) = \eta_4^a(2x)$ and $\eta_1^a(x) = \eta_2^a(2x)$, which are valid for both the BGK and Braginskii solutions.

Now one can consider more general (E11), with the heat flux contributions $\bar{\bar{W}}_a^q$ and the frictional contributions $\bar{\bar{W}}_a^{\text{frict}}$. The solution of (E11) is analogous to (E13), because all of the matrices on the right-hand side are traceless and symmetric. However, it is useful to rewrite the solution in a different form by defining the new matrix

$$
\tilde{\tilde{W}}_a = (\nabla \mathbf{u}_a)^S + \frac{2}{5p_a} (\nabla \mathbf{q}_a)^S, \quad (\text{E18})
$$
and the stress tensor then reads

$$
\begin{align}
\bar{\Pi}_a^{(2)} &= -\eta_0^a \bar{W}_0 - \eta_1^a \bar{W}_1 - \eta_2^a \bar{W}_2 + \eta_3^a \bar{W}_3 + \eta_4^a \bar{W}_4; \\
\bar{W}_0 &= \left[ \frac{3}{2} (\tilde{W}_a : \hat{b}\hat{b}) - \nabla \cdot \mathbf{u}_a - \frac{2}{5p_a} \nabla \cdot \mathbf{q}_a \right] \left( \hat{b}\hat{b} - \frac{\bar{I}}{3} \right) \nonumber \\
&\quad - \frac{3}{2} \frac{\rho_a}{p_a} \left[ \sum_b \nu_{ab} \left( \delta u_\parallel^2 - \frac{1}{3} |\delta u|^2 \right) \right] \left( \hat{b}\hat{b} - \frac{\bar{I}}{3} \right); \\
\bar{W}_1 &= \bar{I}_\perp \cdot \tilde{W}_a \cdot \bar{I}_\perp + \left[ \frac{1}{2} (\tilde{W}_a : \hat{b}\hat{b}) - \nabla \cdot \mathbf{u}_a - \frac{2}{5p_a} \nabla \cdot \mathbf{q}_a \right] \bar{I}_\perp \nonumber \\
&\quad - \frac{\rho_a}{p_a} \sum_b \nu_{ab} \left( \delta \mathbf{u}_\perp \delta \mathbf{u}_\perp - \frac{\bar{I}_\perp}{2} |\delta \mathbf{u}_\perp|^2 \right); \\
\bar{W}_2 &= (\bar{I}_\perp \cdot \tilde{W}_a \cdot \hat{b}\hat{b})^S - \frac{\rho_a}{p_a} \sum_b \nu_{ab} [\delta u_\parallel \hat{b} \delta u_\perp]^S; \\
\bar{W}_3 &= \frac{1}{2} (\hat{b} \times \tilde{W}_a \cdot \bar{I}_\perp)^S - \frac{\rho_a}{2p_a} \sum_b \nu_{ab} [(\hat{b} \times \delta \mathbf{u}) \delta \mathbf{u}_\perp]^S; \\
\bar{W}_4 &= (\hat{b} \times \tilde{W}_a \cdot \hat{b}\hat{b})^S - \frac{\rho_a}{p_a} \sum_b \nu_{ab} [(\hat{b} \times \delta \mathbf{u}) \delta u_\parallel \hat{b}]^S, \tag{E19}
\end{align}
$$

with viscosities (E14). Prescribing $\mathbf{q}_a = 0$ and $\delta \mathbf{u} = 0$ of course recovers (E13).

### E.2. Heat Flux Vector $\mathbf{q}_a$

We consider the 13-moment model, where evolution Equation (C19) becomes

$$
\begin{align}
& \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{b} \times \mathbf{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) \nonumber \\
& \qquad + \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(2)} + \frac{7}{2} \bar{\Pi}_a^{(2)} \cdot \nabla \left( \frac{p_a}{\rho_a} \right) - \frac{1}{\rho_a} (\nabla \cdot \bar{\mathbf{p}}_a) \cdot \bar{\Pi}_a^{(2)}, \nonumber \\
& = \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)}, \tag{E20}
\end{align}
$$

and the BGK collisional contributions calculate

$$ \frac{1}{2} \mathrm{Tr} \bar{Q}_{ab}^{(3)} - \frac{5}{2} \frac{p_a}{\rho_a} \mathbf{R}_{ab} = - \nu_{ab} \mathbf{q}_a + \frac{\nu_{ab}}{2} \rho_a \delta \mathbf{u} |\delta \mathbf{u}|^2. \quad (E21) $$

In a quasistatic approximation, (E20) can be simplified into

$$
\begin{equation}
\begin{split}
\hat{\mathbf{b}} \times \mathbf{q}_a + \frac{\bar{\nu}_a}{\Omega_a} \mathbf{q}_a = -\frac{1}{\Omega_a} \Biggl[ & \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) + \frac{p_a}{\rho_a} \nabla \cdot \bar{\Pi}_a^{(2)} + \frac{7}{2} \bar{\Pi}_a^{(2)} \cdot \nabla \left( \frac{p_a}{\rho_a} \right) - \frac{1}{\rho_a} \left( \nabla \cdot \bar{\mathbf{p}}_a \right) \cdot \bar{\Pi}_a^{(2)} \\
& + \frac{1}{\rho_a} \mathbf{R}_a \cdot \bar{\Pi}_a^{(2)} - \sum_b \frac{\nu_{ab}}{2} \rho_a \delta \mathbf{u} |\delta \mathbf{u}|^2 \Biggr].
\end{split}
\tag{E22}
\end{equation}
$$

A general vector equation (where $\mathbf{a}$ is an unspecified vector, unrelated to the species index)

$$ \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{\bar{\nu}_a}{\Omega_a} \mathbf{q}_a = - \frac{\mathbf{a}}{\Omega_a}, \quad (E23) $$

has the following exact solution (splitting the equation into parallel and perpendicular parts $\mathbf{q}_a = \mathbf{q}_{\parallel a} + \mathbf{q}_{\perp a}$ and $\mathbf{a} = \mathbf{a}_{\parallel} + \mathbf{a}_{\perp}$ with $\hat{\mathbf{b}} \times \mathbf{q}_{\parallel a} = 0$; applying $\hat{\mathbf{b}} \times$ on the perpendicular part; using $\hat{\mathbf{b}} \times (\hat{\mathbf{b}} \times \mathbf{q}_{\perp a}) = - \mathbf{q}_{\perp a}$; and solving the two coupled perpendicular equations by eliminating $\hat{\mathbf{b}} \times \mathbf{q}_{\perp a}$):

$$ \mathbf{q}_a = - \frac{1}{\bar{\nu}_a} (\mathbf{a} \cdot \hat{\mathbf{b}}) \hat{\mathbf{b}} + \frac{\Omega_a}{\Omega_a^2 + \bar{\nu}_a^2} \hat{\mathbf{b}} \times \mathbf{a} - \frac{\bar{\nu}_a}{\Omega_a^2 + \bar{\nu}_a^2} \mathbf{a}_{\perp}. \quad (E24) $$

Note that $\hat{\mathbf{b}} \times \mathbf{a} = \hat{\mathbf{b}} \times \mathbf{a}_{\perp}$. Result (E24) represents the solution of Equation (E22). For zero magnetic field, $\mathbf{q}_a = - \mathbf{a}/\bar{\nu}_a$. The BGK frictional contributions due to $\delta \mathbf{u}$ are only nonlinear; in contrast, the electron heat flux of Braginskii contains frictional $\delta \mathbf{u}$
contributions that are linear. At the semilinear level, (E22) simplifies into

$$ \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{\bar{\nu}_a}{\Omega_a} \mathbf{q}_a = - \frac{1}{\Omega_a} \left[ \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) + \frac{p_a}{\rho_a} \nabla \cdot \bar{\bar{\Pi}}_a^{(2)} \right], \qquad (E25) $$

with a solution again given by (E24). The BGK operator can therefore account for linear (!) contributions of the stress tensor $\bar{\bar{\Pi}}_a^{(2)}$ that enters the heat flux $\mathbf{q}_a$, similar to the previous result for (E19), where the heat $\mathbf{q}_a$ flux entered the stress tensor $\bar{\bar{\Pi}}_a^{(2)}$. Such a coupling is typically neglected with the Landau collisional operator.

The simplest BGK heat flux is a solution of the equation

$$ \hat{\mathbf{b}} \times \mathbf{q}_a + \frac{\bar{\nu}_a}{\Omega_a} \mathbf{q}_a = - \frac{5}{2} \frac{p_a}{\Omega_a m_a} \nabla T_a, \qquad (E26) $$

and the solution reads:

$$ \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, \qquad (E27) $$

with thermal conductivities

$$ \kappa_{\parallel}^{a} = \frac{5}{2} \frac{p_a}{\bar{\nu}_a m_a}; \qquad \kappa_{\perp}^{a} = \frac{5}{2} \frac{p_a}{m_a} \frac{\bar{\nu}_a}{(\Omega_a^2 + \bar{\nu}_a^2)}; \qquad \kappa_{\times}^{a} = \frac{5}{2} \frac{p_a}{m_a} \frac{\Omega_a}{(\Omega_a^2 + \bar{\nu}_a^2)}. \qquad (E28) $$

We use the Braginskii notation with vector $\nabla_{\parallel} = \hat{\mathbf{b}}\hat{\mathbf{b}} \cdot \nabla$. If the magnetic field is zero, so that $\Omega_a = 0$ and $\kappa_{\parallel}^a = \kappa_{\perp}^a$, the solution simplifies into $\mathbf{q}_a = - \kappa_{\parallel}^a \nabla T_a$.

### E.3. BGK versus Braginskii Comparison

Here we compare the BGK viscosities and heat conductivities with those of Braginskii (1965) for a one ion-electron plasma with ion charge $Z_i = 1$. The BGK viscosities (E14) contain $\bar{\nu}_a = \sum_b \nu_{ab}$, and in general should be added according to

$$ \begin{aligned} \bar{\nu}_i &= \nu_{ii} + \nu_{ie} = \nu_{ii} \left( 1 + \frac{\sqrt{2}}{Z_i} \sqrt{\frac{m_e}{m_i}} \left( \frac{T_i}{T_e} \right)^{3/2} \right); \\ \bar{\nu}_e &= \nu_{ee} + \nu_{ei} = \nu_{ei} \left( 1 + \frac{1}{Z_i \sqrt{2}} \right). \end{aligned} \qquad (E29) $$

However, for the ion species, Braginskii neglects ion-electron collisions, and thus $\bar{\nu}_i = \nu_{ii}$ and $\bar{\nu}_e = 1.707\nu_{ei}$; see Section 8.2. Using Braginskii's notation with one-index $\nu_i = \nu_{ii}$ and $\nu_e = \nu_{ei}$ then implies

$$ \bar{\nu}_a = \alpha_a \nu_a; \quad \text{where } \alpha_i = 1; \quad \alpha_e = 1.707, \qquad (E30) $$

and introducing the quantity $x = \Omega_a/\nu_a$, the BGK viscosities (E14) then become

$$ \eta_0^a = \frac{p_a}{\alpha_a \nu_a}; \qquad \eta_1^a = \frac{p_a}{\nu_a 4x^2 + \alpha_a^2}; \qquad \eta_2^a = \frac{p_a}{\nu_a x^2 + \alpha_a^2}; \qquad \eta_3^a = \frac{p_a}{\nu_a 4x^2 + \alpha_a^2}; \qquad \eta_4^a = \frac{p_a}{\nu_a x^2 + \alpha_a^2}. \qquad (E31) $$

Note that $\eta_1^a(x) = \eta_2^a(2x)$ and $\eta_3^a(x) = \eta_4^a(2x)$. Similarly, the BGK heat conductivities (E28) become

$$ \kappa_{\parallel}^{a} = \frac{5}{2\alpha_a} \frac{p_a}{\nu_a m_a}; \qquad \kappa_{\perp}^{a} = \frac{5}{2} \frac{p_a}{m_a} \frac{\alpha_a}{(x^2 + \alpha_a^2)}; \qquad \kappa_{\times}^{a} = \frac{5}{2} \frac{p_a}{m_a} \frac{x}{(x^2 + \alpha_a^2)}. \qquad (E32) $$

The viscosities and heat conductivities for Braginskii are given in the main text. The ion viscosities are compared in Figure 3, the electron viscosities in Figure 4, and the heat conductivities in Figure 5. A small value of $x$ represents a weak magnetic field and a large value of $x$ represents a strong magnetic field.

### E.4. Nonlinear Stress Tensor Decomposition

Here we want to consider the BGK equation for the stress tensor (E12),

$$ (\hat{\mathbf{b}} \times \bar{\bar{\Pi}})^S + \frac{\nu}{\Omega} \bar{\bar{\Pi}} = - \frac{p}{\Omega} \bar{\bar{W}}, \qquad (E33) $$

and clarify the solution (E13). Species indices are dropped, and both $\bar{\bar{\Pi}}$ and $\bar{\bar{W}}$ are symmetric and traceless. First, we need to learn how to decompose any general matrix. It is useful for a moment to consider the undefined matrices $\bar{\bar{W}}$ and $\bar{\bar{\Pi}}$, which are not necessarily symmetric or traceless.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Ion viscosities of the BGK model (red) and of the Braginskii model (black), normalized as $\hat{\eta}^i = \eta^i \nu_{ii} / p_i$ vs. the ratio $x = \Omega_{ii} / \nu_{ii}$. Left: perpendicular viscosities $\eta_1^i, \eta_2^i$. Right: gyroviscosities $\eta_3^i, \eta_4^i$. For large values of $x$, the BGK asymptotic profiles for $\hat{\eta}_3^i = 1/(2x)$ and $\hat{\eta}_4^i = 1/x$ become independent of collisional frequencies and match the asymptotic profiles of Braginskii exactly. The BGK asymptotic profiles for $\hat{\eta}_1^i = 1/(4x^2)$ and $\hat{\eta}_2^i = 1/x^2$ have the correct functional dependence, but differ from the Braginskii asymptotes by a proportionality constant. The BGK operator reproduces the ion viscosity of Braginskii with surprisingly good accuracy.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Electron viscosities, normalized as $\hat{\eta}^e = \eta^e \nu_{ei} / p_e$ vs. the ratio $|x| = |\Omega_e| / \nu_{ei}$. The results are less precise than for the ions in Figure 3, especially for small values of $x$; nevertheless, the same conclusions are obtained.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Heat conductivities $\kappa_{\perp}^a$ and $\kappa_{\parallel}^a$. Left: ion species, normalized as $\kappa^i m_i \nu_{ii} / p_i$. Right: electron species, normalized as $\kappa^e m_e \nu_{ei} / p_e$. For large values of $x$, the BGK asymptotic profiles $\kappa_{\perp}^a$ (dashed lines) match the Braginskii results exactly, whereas for $\kappa_{\parallel}^a$ (solid lines) the results differ by a proportionality constant.</figcaption></figure>
### E.4.1. Decomposition of a Matrix

We will work both in the reference frame of the magnetic field lines ($\hat{b}_0 = (0, 0, 1)$), which nicely guide and clarify the calculations, and also in a laboratory reference frame with general $\hat{b}$. In the reference frame of the magnetic field lines, one uses the matrices

$$ \hat{b}\hat{b} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}; \quad \bar{\mathcal{I}}_{\perp} = \bar{\mathcal{I}} - \hat{b}\hat{b} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}; \quad \bar{\mathcal{I}}^{\times} = \begin{pmatrix} 0 & -1 & 0 \\ +1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (E34) $$

where the last matrix is defined as $\hat{b} \times \bar{\mathcal{W}} = (\bar{\mathcal{I}}^{\times}) \cdot \bar{\mathcal{W}}$. Then, one takes a general matrix $\bar{\mathcal{W}}$, and starts multiplying it with the matrices $\hat{b}\hat{b}$ and $\bar{\mathcal{I}}_{\perp}$ from the left and right, yielding a general decomposition:

$$ \begin{aligned} \bar{\mathcal{W}} &= \bar{\mathcal{W}}_0' + \bar{\mathcal{W}}_1' + \bar{\mathcal{W}}_2; \\ \bar{\mathcal{W}}_0' &= \hat{b}\hat{b} \cdot \bar{\mathcal{W}} \cdot \hat{b}\hat{b} = (\bar{\mathcal{W}}: \hat{b}\hat{b})\hat{b}\hat{b}; \\ \bar{\mathcal{W}}_1' &= \bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \bar{\mathcal{I}}_{\perp}; \\ \bar{\mathcal{W}}_2' &= \bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \hat{b}\hat{b} + \hat{b}\hat{b} \cdot \bar{\mathcal{W}} \cdot \bar{\mathcal{I}}_{\perp} = (\bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \hat{b}\hat{b})^s. \end{aligned} \qquad (E35) $$

In the reference frame of the magnetic field lines,

$$ \bar{\mathcal{W}}_0' = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & W_{zz} \end{pmatrix}; \quad \bar{\mathcal{W}}_1' = \begin{pmatrix} W_{xx} & W_{xy} & 0 \\ W_{yx} & W_{yy} & 0 \\ 0 & 0 & 0 \end{pmatrix}; \quad \bar{\mathcal{W}}_2 = \begin{pmatrix} 0 & 0 & W_{xz} \\ 0 & 0 & W_{yz} \\ W_{zx} & W_{zy} & 0 \end{pmatrix}, \qquad (E36) $$

and adding these matrices together obviously yields the full matrix $\bar{\mathcal{W}}$. However, the decomposition (E35) also works in the laboratory reference frame with general $\hat{b}$, as can be verified by adding the general matrices together. It is possible to consider an alternative decomposition, according to

$$ \begin{aligned} \bar{\mathcal{W}} &= \bar{\mathcal{W}}_0 + \bar{\mathcal{W}}_1 + \bar{\mathcal{W}}_2; \\ \bar{\mathcal{W}}_0 &= (\bar{\mathcal{W}}: \hat{b}\hat{b})\hat{b}\hat{b} + \frac{1}{2}(\bar{\mathcal{W}}: \bar{\mathcal{I}}_{\perp})\bar{\mathcal{I}}_{\perp}; \\ \bar{\mathcal{W}}_1 &= \bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \bar{\mathcal{I}}_{\perp} - \frac{1}{2}(\bar{\mathcal{W}}: \bar{\mathcal{I}}_{\perp})\bar{\mathcal{I}}_{\perp}; \\ \bar{\mathcal{W}}_2 &= (\bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \hat{b}\hat{b})^s, \end{aligned} \qquad (E37) $$

where in the reference frame of the magnetic field lines,

$$ \bar{\mathcal{W}}_0 = \begin{pmatrix} \frac{1}{2}(W_{xx} + W_{yy}) & 0 & 0 \\ 0 & \frac{1}{2}(W_{xx} + W_{yy}) & 0 \\ 0 & 0 & W_{zz} \end{pmatrix}; \quad \bar{\mathcal{W}}_1 = \begin{pmatrix} \frac{1}{2}(W_{xx} - W_{yy}) & W_{xy} & 0 \\ W_{yx} & -\frac{1}{2}(W_{xx} - W_{yy}) & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (E38) $$

The decomposition (E37) again works for general $\hat{b}$, and in comparison to the previous decomposition $\bar{\mathcal{W}}_0' + \bar{\mathcal{W}}_1' = \bar{\mathcal{W}}_0 + \bar{\mathcal{W}}_1$. The advantage is that if $\bar{\mathcal{W}}$ is traceless, then all three matrices are traceless. It is useful to rearrange (E37), by separating the trace of $\bar{\mathcal{W}}$ with $(\bar{\mathcal{W}}: \bar{\mathcal{I}}_{\perp})\bar{\mathcal{I}}_{\perp} = (\bar{\mathcal{W}}: \bar{\mathcal{I}})\bar{\mathcal{I}}_{\perp} - (\bar{\mathcal{W}}: \hat{b}\hat{b})\bar{\mathcal{I}}_{\perp}$, yielding the decomposition

$$ \begin{aligned} \bar{\mathcal{W}} &= \bar{\mathcal{W}}_0 + \bar{\mathcal{W}}_1 + \bar{\mathcal{W}}_2; \\ \bar{\mathcal{W}}_0 &= \frac{3}{2}(\bar{\mathcal{W}}: \hat{b}\hat{b})\left(\hat{b}\hat{b} - \frac{\bar{\mathcal{I}}}{3}\right) + \frac{1}{2}(\bar{\mathcal{W}}: \bar{\mathcal{I}})\bar{\mathcal{I}}_{\perp}; \\ \bar{\mathcal{W}}_1 &= \bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \bar{\mathcal{I}}_{\perp} + \frac{1}{2}(\bar{\mathcal{W}}: \hat{b}\hat{b})\bar{\mathcal{I}}_{\perp} - \frac{1}{2}(\bar{\mathcal{W}}: \bar{\mathcal{I}})\bar{\mathcal{I}}_{\perp}; \\ \bar{\mathcal{W}}_2 &= (\bar{\mathcal{I}}_{\perp} \cdot \bar{\mathcal{W}} \cdot \hat{b}\hat{b})^s. \end{aligned} \qquad (E39) $$
The same decomposition is used for the stress tensor $\bar{\bar{\Pi}}$:

$$ 
\begin{aligned}
\bar{\bar{\Pi}} &= \bar{\bar{\Pi}}_0 + \bar{\bar{\Pi}}_1 + \bar{\bar{\Pi}}_2; \\
\bar{\bar{\Pi}}_0 &= \frac{3}{2}(\bar{\bar{\Pi}}: \hat{b}\hat{b}) \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right) + \frac{1}{2}(\bar{\bar{\Pi}}: \bar{\bar{I}})\bar{\bar{I}}_{\perp}; \\
\bar{\bar{\Pi}}_1 &= \bar{\bar{I}}_{\perp} \cdot \bar{\bar{\Pi}} \cdot \bar{\bar{I}}_{\perp} + \frac{1}{2}(\bar{\bar{\Pi}}: \hat{b}\hat{b})\bar{\bar{I}}_{\perp} - \frac{1}{2}(\bar{\bar{\Pi}}: \bar{\bar{I}})\bar{\bar{I}}_{\perp}; \\
\bar{\bar{\Pi}}_2 &= (\bar{\bar{I}}_{\perp} \cdot \bar{\bar{\Pi}} \cdot \hat{b}\hat{b})^S.
\end{aligned}
\quad (E40) $$

Let us solve for $\bar{\bar{\Pi}}_0$. Applying $: \hat{b}\hat{b}$ and $: \bar{\bar{I}}$ at Equation (E33) and using the identities

$$ (\hat{b} \times \bar{\bar{\Pi}})^S: \hat{b}\hat{b} = 0; \quad (\hat{b} \times \bar{\bar{\Pi}})^S: \bar{\bar{I}} = 0, \quad (E41) $$

yields

$$ \bar{\bar{\Pi}}: \hat{b}\hat{b} = -\frac{p}{\nu}\bar{\bar{W}}: \hat{b}\hat{b}; \quad \bar{\bar{\Pi}}: \bar{\bar{I}} = -\frac{p}{\nu}\bar{\bar{W}}: \bar{\bar{I}}, \quad (E42) $$

and plugging these results into (E40) yields the final solution for the parallel stress tensor:

$$ \bar{\bar{\Pi}}_0 = -\frac{p}{\nu} \left[ \frac{3}{2} (\bar{\bar{W}}: \hat{b}\hat{b}) \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right) + \frac{1}{2} (\bar{\bar{W}}: \bar{\bar{I}}) \bar{\bar{I}}_{\perp} \right] = -\frac{p}{\nu} \bar{\bar{W}}_0. \quad (E43) $$

The solution is valid for any general matrix $\bar{\bar{W}}$ (not necessarily symmetric or traceless). If this result is compared with the expression (4.42) of Braginskii (1965), given below by (E46), one notices

$$ \bar{\bar{W}}_{0}^{\text{BR}} = (E46) = \frac{3}{2}(\bar{\bar{W}}: \hat{b}\hat{b})\left(\hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3}\right) + \frac{1}{2}(\bar{\bar{W}}: \bar{\bar{I}})\left(\frac{\bar{\bar{I}}}{3} - \hat{b}\hat{b}\right) \neq \bar{\bar{W}}_{0}, \quad (E44) $$

and his result is valid only if $\bar{\bar{W}}$ is traceless (which it is). The reason why Braginskii left his result in form (E46), and did not simplify it with $\bar{\bar{W}}: \bar{\bar{I}} = 0$, is likely an alternative form (E47).

### E.4.2. Symmetric and Traceless Matrices

We further consider only the symmetric and traceless matrices $\bar{\bar{W}}$ and $\bar{\bar{\Pi}}$, so all of the previous expressions are simplified with $\bar{\bar{W}}: \bar{\bar{I}} = 0$, $\bar{\bar{\Pi}}: \bar{\bar{I}} = 0$, and the BGK parallel stress tensor $\bar{\bar{\Pi}}_0 = -(p/\nu)\bar{\bar{W}}_0$. For clarity, it is useful to write several possible forms for

$$ \bar{\bar{W}}_0 = \frac{3}{2} (\bar{\bar{W}}: \hat{b}\hat{b}) \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right); \quad (E45) $$

$$ = \frac{3}{2} \left[ \bar{\bar{W}}: \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right) \right] \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right); \quad (E46) $$

$$ = \frac{3}{2} \left[ (\nabla u)^S: \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right) \right] \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right); \quad (E47) $$

$$ = 3 \left[ (\nabla u): \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right) \right] \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right); \quad (E48) $$

$$ = 3 \left[ \hat{b} \cdot (\nabla u) \cdot \hat{b} - \frac{1}{3} \nabla \cdot u \right] \left( \hat{b}\hat{b} - \frac{\bar{\bar{I}}}{3} \right). \quad (E49) $$

Braginskii uses (E46), Fitzpatrick (2015), for example, uses (E48), and we use (E45). In the reference frame of the magnetic field lines,

$$ \bar{\bar{W}}_0 = \frac{3}{2} W_{zz} \begin{pmatrix} -1/3 & 0 & 0 \\ 0 & -1/3 & 0 \\ 0 & 0 & +2/3 \end{pmatrix}; \quad \bar{\bar{\Pi}}_0 = \frac{p}{\nu} W_{zz} \begin{pmatrix} 1/2 & 0 & 0 \\ 0 & 1/2 & 0 \\ 0 & 0 & -1 \end{pmatrix}. \quad (E50) $$
To solve Equation (E33), it is beneficial to introduce two other matrices $\bar{W}_3$, $\bar{W}_4$, by decomposing

$$ 
\begin{aligned}
(\hat{b} \times \bar{W})^S &= 2\bar{W}_3 + \bar{W}_4; \\
2\bar{W}_3 &= (\hat{b} \times \bar{W} \cdot \bar{I}_{\perp})^S; \\
\bar{W}_4 &= (\hat{b} \times \bar{W} \cdot \hat{b}\hat{b})^S,
\end{aligned}
\qquad (E51) $$

where in the reference frame of the magnetic field lines,

$$ 2\bar{W}_3 = \begin{pmatrix} -2W_{xy} & W_{xx} - W_{yy} & 0 \\ W_{xx} - W_{yy} & 2W_{xy} & 0 \\ 0 & 0 & 0 \end{pmatrix}; \qquad \bar{W}_4 = \begin{pmatrix} 0 & 0 & -W_{yz} \\ 0 & 0 & W_{xz} \\ -W_{yz} & W_{xz} & 0 \end{pmatrix}. \qquad (E52) $$

The decomposition (E51) is again valid for general $\hat{b}$, which is easily verified using $\bar{I}_{\perp} + \hat{b}\hat{b} = \bar{I}$, and the stress tensor is decomposed in the same way:

$$ 
\begin{aligned}
(\hat{b} \times \bar{\Pi})^S &= 2\bar{\Pi}_3 + \bar{\Pi}_4; \\
2\bar{\Pi}_3 &= (\hat{b} \times \bar{\Pi} \cdot \bar{I}_{\perp})^S; \\
\bar{\Pi}_4 &= (\hat{b} \times \bar{\Pi} \cdot \hat{b}\hat{b})^S.
\end{aligned}
\qquad (E53) $$

Finally, applying $\hat{b} \times$ at the matrices $\bar{W}_0$ ... $\bar{W}_4$ yields the following identities:

$$ 
\begin{aligned}
(\hat{b} \times \bar{W}_0)^S &= 0; & \hat{b} \times \bar{W}_1 &= \bar{W}_3; & (\hat{b} \times \bar{W}_2)^S &= \bar{W}_4; \\
\hat{b} \times \bar{W}_3 &= -\bar{W}_1; & (\hat{b} \times \bar{W}_4)^S &= -\bar{W}_2,
\end{aligned}
\qquad (E54) $$

which are easy to verify in a general reference frame with $\hat{b}$. The same identities hold for the stress tensor

$$ 
\begin{aligned}
(\hat{b} \times \bar{\Pi}_0)^S &= 0; & \hat{b} \times \bar{\Pi}_1 &= \bar{\Pi}_3; & (\hat{b} \times \bar{\Pi}_2)^S &= \bar{\Pi}_4; \\
\hat{b} \times \bar{\Pi}_3 &= -\bar{\Pi}_1; & (\hat{b} \times \bar{\Pi}_4)^S &= -\bar{\Pi}_2.
\end{aligned}
\qquad (E55) $$

### E.4.3. Final Solution

Now we are ready to solve Equation (E33), which is rewritten as

$$ 2\bar{\Pi}_3 + \bar{\Pi}_4 + \frac{\nu}{\Omega}(\bar{\Pi}_0 + \bar{\Pi}_1 + \bar{\Pi}_2) = -\frac{p}{\Omega}(\bar{W}_0 + \bar{W}_1 + \bar{W}_2). \qquad (E56) $$

One solution, $\bar{\Pi}_0 = -(p/\nu)\bar{W}_0$, has already been obtained, and can be eliminated from (E56). For the rest of the equation, the most straightforward approach is to be guided by the reference frame of the magnetic field lines, which shows that the system (E56) can be directly split into two independent equations:

$$ 2\bar{\Pi}_3 + \frac{\nu}{\Omega}\bar{\Pi}_1 = -\frac{p}{\Omega}\bar{W}_1; \qquad (E57) $$

$$ \bar{\Pi}_4 + \frac{\nu}{\Omega}\bar{\Pi}_2 = -\frac{p}{\Omega}\bar{W}_2. \qquad (E58) $$

In the general reference frame, the split can be achieved by applying $\bar{I}_{\perp}$ from left and right at (E56), for example, which, using identities $\bar{I}_{\perp} \cdot \bar{\Pi}_4 \cdot \bar{I}_{\perp} = 0$, $\bar{I}_{\perp} \cdot \bar{\Pi}_2 \cdot \bar{I}_{\perp} = 0$, and $\bar{I}_{\perp} \cdot \bar{W}_2 \cdot \bar{I}_{\perp} = 0$, yields (E57) and subsequently (E58). The split significantly simplifies the “inversion procedure.”

Furthermore, applying $\hat{b} \times$ at (E57), applying $\hat{b} \times$ together with the symmetric operator at (E58), and using the identities (E54)–(E55) then gives

$$ -2\bar{\Pi}_1 + \frac{\nu}{\Omega}\bar{\Pi}_3 = -\frac{p}{\Omega}\bar{W}_3; \qquad (E59) $$

$$ -\bar{\Pi}_2 + \frac{\nu}{\Omega}\bar{\Pi}_4 = -\frac{p}{\Omega}\bar{W}_4. \qquad (E60) $$

Equations (E57), (E59) are coupled and can be treated as two equations in two unknowns, and similarly Equations (E58), (E60), finally yielding the solutions

$$ \bar{\Pi}_1 = -\frac{p\nu}{4\Omega^2 + \nu^2}\bar{W}_1 + \frac{2p\Omega}{4\Omega^2 + \nu^2}\bar{W}_3; \qquad (E61) $$

$$ \bar{\Pi}_2 = -\frac{p\nu}{\Omega^2 + \nu^2}\bar{W}_2 + \frac{p\Omega}{\Omega^2 + \nu^2}\bar{W}_4. \qquad (E62) $$
The entire solution for the stress tensor $\bar{\bar{\Pi}} = \bar{\bar{\Pi}}_0 + \bar{\bar{\Pi}}_1 + \bar{\bar{\Pi}}_2$ thus reads

$$ \bar{\bar{\Pi}} = -\frac{p}{\nu}\bar{\bar{W}}_0 - \frac{p\nu}{4\Omega^2 + \nu^2}\bar{\bar{W}}_1 - \frac{p\nu}{\Omega^2 + \nu^2}\bar{\bar{W}}_2 + \frac{2p\Omega}{4\Omega^2 + \nu^2}\bar{\bar{W}}_3 + \frac{p\Omega}{\Omega^2 + \nu^2}\bar{\bar{W}}_4. \quad (E63) $$

### E.5. BGK Operator and Electric Field

The BGK operator is also an excellent tool for clarifying various processes in fully ionized or partially ionized plasmas. Here we want to clarify the ohmic (magnetic) diffusion together with the ambipolar diffusion, both caused by the momentum exchange rates

$$ \mathbf{R}_a = \sum_{b \neq a} \rho_a \nu_{ab} (\mathbf{u}_b - \mathbf{u}_a). \quad (E64) $$

From the BGK perspective, one does not need to worry about complicated Landau and Boltzmann operators, and simply “adopts” the correct collisional frequencies; see, for example, Appendix C of Schunk (1977). The momentum exchange rates (E64) are actually the correct answer if the relative drift velocities are small and one considers the 5-moment model (i.e., if the heat flux is neglected).

We restrict our focus to spatial scales much longer than the Debye length. The displacement current is neglected, the Gauss law $\nabla \cdot \mathbf{E} = 4\pi e \sum_a Z_a n_a$ is replaced by the charge neutrality, and no condition is placed on $\nabla \cdot \mathbf{E}$. The Maxwell equations then read

$$ \sum_a Z_a n_a = 0; \quad \mathbf{j} = \sum_a e Z_a n_a \mathbf{u}_a = \frac{c}{4\pi} \nabla \times \mathbf{B}; \quad (E65) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \nabla \times \mathbf{E}; \quad \nabla \cdot \mathbf{B} = 0. \quad (E66) $$

By focusing on the spatial and temporal scales of the ion and neutral species, we do not need to resolve the electron motion. In the electron momentum equation, the electron inertia represented by $d_e \mathbf{u}_e/dt$ is neglected (which does not mean that $m_e = 0$; relations $\rho_a \nu_{ab} = \rho_b \nu_{ba}$ still hold), and the electric field is expressed as

$$ \mathbf{E} = -\frac{1}{c}\mathbf{u}_e \times \mathbf{B} - \frac{1}{en_e}\nabla \cdot \bar{\bar{\mathbf{p}}}_e + \frac{\mathbf{R}_e}{en_e}. \quad (E67) $$

The momentum equations for the ions then become

$$ \rho_i \frac{d_i \mathbf{u}_i}{dt} + \nabla \cdot \bar{\bar{\mathbf{p}}}_i + \frac{Z_i n_i}{n_e} \nabla \cdot \bar{\bar{\mathbf{p}}}_e - \frac{e Z_i n_i}{c} (\mathbf{u}_i - \mathbf{u}_e) \times \mathbf{B} = \mathbf{R}_i + \frac{Z_i n_i}{n_e} \mathbf{R}_e. \quad (E68) $$

Also, by using (E65), the electron density $n_e$ and electron velocity $\mathbf{u}_e$ are expressed as

$$ n_e = \sum_i Z_i n_i; \quad \mathbf{u}_e = \frac{1}{n_e} \sum_i Z_i n_i \mathbf{u}_i - \frac{\mathbf{j}}{en_e}; \quad \mathbf{j} = \frac{c}{4\pi} \nabla \times \mathbf{B}, \quad (E69) $$

where the summations are over ion species. The electron density equation $\partial n_e/\partial t + \nabla \cdot (n_e \mathbf{u}_e) = 0$ becomes redundant, because multiplying all of the density equations for the charges (including electrons) by $Z_a$ and summing them together yields a requirement $\nabla \cdot (\sum_a Z_a n_a \mathbf{u}_a) = 0$, which is automatically satisfied by $\nabla \cdot \mathbf{j} = 0$ in (E69). Expressions (E69) and (E67) can then be substituted to all of the other equations (which is easy to do numerically), and the occurrences of $\mathbf{E}$, $\mathbf{u}_e$, and $n_e$ in the entire model are thus eliminated.

For a particular case of $\mathbf{R}_e$ given by (E64), the electric field (E67) then becomes

$$ \begin{aligned} \mathbf{E} = & -\frac{1}{cn_e} \left( \sum_i Z_i n_i \mathbf{u}_i \right) \times \mathbf{B} + \frac{\mathbf{j} \times \mathbf{B}}{cen_e} - \frac{1}{en_e} \nabla \cdot \bar{\bar{\mathbf{p}}}_e + \frac{m_e}{e^2 n_e} \mathbf{j} \sum_{a \neq e} \nu_{ea} \\ & + \frac{m_e}{e} \left[ \left( \sum_{a \neq e} \nu_{ea} \mathbf{u}_a \right) - \frac{1}{n_e} \left( \sum_{a \neq e} \nu_{ea} \right) \left( \sum_i Z_i n_i \mathbf{u}_i \right) \right]. \end{aligned} \quad (E70) $$

The summations over “a” include both ions and neutrals. The terms on the right-hand-side can be called the convective term, the Hall term, the electron pressure term, the ohmic term, and a mixed collisional term due to ion and neutral velocities, respectively. When (E70) is used in the induction equation, the ohmic term ($\sim \mathbf{j}$) becomes directly diffusive through the identity $\nabla \times (\eta_B \nabla \times \mathbf{B}) = -\eta_B \nabla^2 \mathbf{B} + \nabla(\eta_B) \times (\nabla \times \mathbf{B})$, where one defines a coefficient of magnetic diffusion $\eta_B = (\sum_{a \neq e} \nu_{ea}) m_e c^2 / (4\pi e^2 n_e)$. In contrast, no other term in (E70) is directly diffusive in this sense. Nevertheless, the so-called ambipolar diffusion due to the differences in the velocities $\mathbf{u}_a$ between different species is still present implicitly, which can be shown by solving the dispersion relations. The explicit presence of ambipolar diffusion caused by $\sim (\mathbf{j} \times \mathbf{B}) \times \mathbf{B} = \mathbf{j} \perp |\mathbf{B}|^2$ is revealed by the construction of a single-fluid model, formulated with respect to the center-of-mass velocity of all of the species. In general, ambipolar diffusion between two species with indices (a, b) exists if

$$ \frac{Z_a}{m_a} \neq \frac{Z_b}{m_b}, \quad (E71) $$

which is demonstrated in Appendix E.6.
In partially ionized solar plasmas, one often focuses on a two-fluid model, formulated with center-of-mass velocities for the ion species $\langle \mathbf{u}_i \rangle = (\sum_i \rho_i \mathbf{u}_i)/\sum_i \rho_i$ and for the neutral species $\langle \mathbf{u}_n \rangle = (\sum_n \rho_n \mathbf{u}_n)/\sum_n \rho_n$. The velocities for each species are thus decomposed into $\mathbf{u}_i = \langle \mathbf{u}_i \rangle + \mathbf{w}_i$, $\mathbf{u}_n = \langle \mathbf{u}_n \rangle + \mathbf{w}_n$, where $\mathbf{w}_i$, $\mathbf{w}_n$ represent drifts; and because $\langle \mathbf{u}_i \rangle$, $\langle \mathbf{u}_n \rangle$ can be pulled out in front of the summations, the electric field (E70) transforms into

$$ 
\begin{aligned}
\mathbf{E} = & -\frac{1}{c} \langle \mathbf{u}_i \rangle \times \mathbf{B} - \frac{1}{cn_e} \left( \sum_i Z_i n_i \mathbf{w}_i \right) \times \mathbf{B} + \frac{\mathbf{j} \times \mathbf{B}}{cen_e} - \frac{1}{en_e} \nabla \cdot \bar{\mathbf{p}}_e + \frac{m_e}{e^2 n_e} \mathbf{j} \sum_{a \neq e} \nu_{ea} \\
& + (\langle \mathbf{u}_n \rangle - \langle \mathbf{u}_i \rangle) \frac{m_e}{e} \sum_n \nu_{en} \\
& + \frac{m_e}{e} \left[ \left( \sum_{a \neq e} \nu_{ea} \mathbf{w}_a \right) - \frac{1}{n_e} \left( \sum_{a \neq e} \nu_{ea} \right) \left( \sum_i Z_i n_i \mathbf{w}_i \right) \right].
\end{aligned}
\quad (E72) $$

The electric field (E72) still represents a multifluid electric field, where one considers separate evolution equations for all of the drifts $\mathbf{w}_a$. To obtain a two-fluid electric field, these drifts have to be somehow eliminated, which is of course not straightforward to justify. In partially ionized solar plasmas, the usual justification is that (1) one takes into account only the first ionization degree, with all of the ions having $Z_i = 1$; (2) one prescribes that on average $\sum_i n_i \mathbf{w}_i = 0$ (which, for example, eliminates ambipolar diffusion between different ions), together with $\sum_n n_n \mathbf{w}_n = 0$; (3) all of the species have roughly the same temperature, which, using the collisional frequencies $\nu_{ei} = n_e f(T)/\sqrt{m_e}$, yields $\sum_i \nu_{ei} \mathbf{w}_i = 0$; and (4) all of the neutrals have roughly the same cross sections (radii $r_n$), which, using $\nu_{en} = n_n f(T) r_n^2 / \sqrt{m_e}$, yields $\sum_n \nu_{en} \mathbf{w}_n = 0$. The two-fluid electric field thus reads

$$ 
\begin{aligned}
\mathbf{E} = & -\frac{1}{c} \langle \mathbf{u}_i \rangle \times \mathbf{B} + \frac{\mathbf{j} \times \mathbf{B}}{cen_e} - \frac{1}{en_e} \nabla \cdot \bar{\mathbf{p}}_e + \frac{m_e}{e^2 n_e} \mathbf{j} \sum_{a \neq e} \nu_{ea} \\
& + (\langle \mathbf{u}_n \rangle - \langle \mathbf{u}_i \rangle) \frac{m_e}{e} \sum_n \nu_{en}.
\end{aligned}
\quad (E73) $$

The center-of-mass velocity for the ions $\langle \mathbf{u}_i \rangle$ can be freely replaced by the center-of-mass velocity for all the charges $\langle \mathbf{u}_e \rangle$ (which include electrons). Then, the electric field (E73) is almost identical to Equation (115) of Khomenko et al. (2014), except that the $\sum_n \nu_{en}$ in the last term of (E73) is replaced by $(\sum_n \nu_{en}) - (\sum_i \nu_{in})$ in that paper. The difference arises from the alternative approach in that paper, where the electron inertia is not neglected from the beginning, but instead the electric field is derived by first summing momentum equations for all of the species together, and prescribing quasistatic current $j$. Then, the subsequent expansion in the mass ratios retains the contributions from $\mathbf{R}_i$. Nevertheless, the missing contributions are small $\nu_{in} \ll \nu_{en}$, explaining the small difference between these two approaches.

For a particular case of only one ion species and one neutral species, so that $n_e = Z_i n_i$ and $\mathbf{u}_e = \mathbf{u}_i - \mathbf{j}/(en_e)$, the electric field (E70) simplifies into

$$ 
\begin{aligned}
\mathbf{E} = & -\frac{1}{c} \mathbf{u}_i \times \mathbf{B} + \frac{\mathbf{j} \times \mathbf{B}}{cen_e} - \frac{1}{en_e} \nabla \cdot \bar{\mathbf{p}}_e + \frac{m_e}{e^2 n_e} \mathbf{j} (\nu_{ei} + \nu_{en}) \\
& + \frac{m_e}{e} \nu_{en} (\mathbf{u}_n - \mathbf{u}_i);
\end{aligned}
\quad (E74) $$

$$ 
\begin{aligned}
\frac{\partial \mathbf{B}}{\partial t} = & \nabla \times (\mathbf{u}_i \times \mathbf{B}) - \nabla \times \left( \frac{\mathbf{j}}{en_e} \times \mathbf{B} \right) + \frac{c}{e} \nabla \times \left( \frac{1}{n_e} \nabla \cdot \bar{\mathbf{p}}_e \right) \\
& - \nabla \times (\eta_B \nabla \times \mathbf{B}) - \nabla \times \left[ \frac{cm_e}{e} \nu_{en} (\mathbf{u}_n - \mathbf{u}_i) \right],
\end{aligned}
\quad (E75) $$

with the coefficient of magnetic diffusion $\eta_B = (\nu_{ei} + \nu_{en})m_e c^2 / (4\pi e^2 n_e)$.

## E.6. Ambipolar Diffusion of Two Ion Species

Here we consider a two-fluid model, consisting of two different ion species with species indices $(i, j)$, so that the charge neutrality reads $n_e = Z_i n_i + Z_j n_j$. A particular case consisting of one ion and one neutral species can be obtained by prescribing $Z_j = 0$ and index $j = n$ (or $Z_i = 0$ and $i = n$). The momentum equations are

$$ \rho_i \frac{d_i \mathbf{u}_i}{dt} + \nabla \cdot \bar{\mathbf{p}}_i + \frac{Z_i n_i}{n_e} \nabla \cdot \bar{\mathbf{p}}_e - \frac{e Z_i n_i Z_j n_j}{c n_e} (\mathbf{u}_i - \mathbf{u}_j) \times \mathbf{B} - \frac{Z_i n_i}{cn_e} \mathbf{j} \times \mathbf{B} = \mathbf{R}_i + \frac{Z_i n_i}{n_e} \mathbf{R}_e; \quad (E76) $$

$$ \rho_j \frac{d_j \mathbf{u}_j}{dt} + \nabla \cdot \bar{\mathbf{p}}_j + \frac{Z_j n_j}{n_e} \nabla \cdot \bar{\mathbf{p}}_e + \frac{e Z_j n_j Z_j n_j}{c n_e} (\mathbf{u}_i - \mathbf{u}_j) \times \mathbf{B} - \frac{Z_j n_j}{cn_e} \mathbf{j} \times \mathbf{B} = \mathbf{R}_j + \frac{Z_j n_j}{n_e} \mathbf{R}_e, \quad (E77) $$
with the collisional right-hand sides being

$$
\begin{align}
\mathbf{R}_i + \frac{Z_i n_i}{n_e} \mathbf{R}_e &= -(\mathbf{u}_i - \mathbf{u}_j) \left[ \rho_i \nu_{ij} + \rho_e \nu_{ei} \left( \frac{Z_j n_j}{n_e} \right)^2 + \rho_e \nu_{ej} \left( \frac{Z_i n_i}{n_e} \right)^2 \right] \nonumber \\
&\quad -j \frac{m_e}{en_e} (\nu_{ei} Z_j n_j - \nu_{ej} Z_i n_i); \nonumber \\
\mathbf{R}_j + \frac{Z_j n_j}{n_e} \mathbf{R}_e &= (\mathbf{u}_i - \mathbf{u}_j) \left[ \rho_i \nu_{ij} + \rho_e \nu_{ei} \left( \frac{Z_j n_j}{n_e} \right)^2 + \rho_e \nu_{ej} \left( \frac{Z_i n_i}{n_e} \right)^2 \right] \nonumber \\
&\quad +j \frac{m_e}{en_e} (\nu_{ei} Z_j n_j - \nu_{ej} Z_i n_i), \tag{E78}
\end{align}
$$

and the electric field (which determines the induction equation) reading

$$
\begin{equation}
\begin{aligned}
\mathbf{E} = & -\frac{1}{cn_e}(Z_i n_i \mathbf{u}_i + Z_j n_j \mathbf{u}_j) \times \mathbf{B} + \frac{j \times \mathbf{B}}{cen_e} - \frac{1}{en_e} \nabla \cdot \bar{\mathbf{p}}_e + \frac{m_e}{e^2 n_e} \mathbf{j}(\nu_{ei} + \nu_{ej}) \\
& + \frac{m_e}{en_e}(\mathbf{u}_i - \mathbf{u}_j)(Z_j n_j \nu_{ei} - Z_i n_i \nu_{ej}).
\end{aligned}
\tag{E79}
\end{equation}
$$

The ambipolar diffusion term $-j \times B \times B = j_{\perp} |B|^2$ is not directly present in the electric field, and the only term that directly causes magnetic diffusion in the induction equation is the ohmic term ($\sim j$). Nevertheless, the ambipolar diffusion is still present implicitly, which can be shown by solving the dispersion relations or by constructing a single-fluid model.

Using the same notation as Zaqarashvili et al. (2011), and introducing the center-of-mass velocity $V = (\rho_i \mathbf{u}_i + \rho_j \mathbf{u}_j)/\rho$, where the total density $\rho = \rho_i + \rho_j$ and the difference in velocities $\mathbf{w} = \mathbf{u}_i - \mathbf{u}_j$, so that $\mathbf{u}_i = V + (\rho_j/\rho)\mathbf{w}$, $\mathbf{u}_j = V - (\rho_i/\rho)\mathbf{w}$, yields the momentum equations

$$
\rho \frac{\partial \mathbf{V}}{\partial t} + \rho \mathbf{V} \cdot \nabla \mathbf{V} + \nabla \cdot (\bar{\mathbf{p}}_i + \bar{\mathbf{p}}_j + \bar{\mathbf{p}}_e) - \frac{1}{c} \mathbf{j} \times \mathbf{B} + \nabla \cdot \left( \frac{\rho_i \rho_j}{\rho} \mathbf{w} \mathbf{w} \right) = 0; \quad (\text{E80})
$$

$$
\begin{align}
& \frac{\partial \mathbf{w}}{\partial t} + \mathbf{w} \cdot \nabla \mathbf{V} + \mathbf{V} \cdot \nabla \mathbf{w} + \frac{\rho_j}{\rho} \mathbf{w} \cdot \nabla \mathbf{w} - \mathbf{w} \cdot \nabla \left( \frac{\rho_i}{\rho} \mathbf{w} \right) - \frac{e Z_i n_i Z_j n_j}{c n_e} \frac{\rho}{\rho_i \rho_j} \mathbf{w} \times \mathbf{B} \nonumber \\
& \qquad + \frac{1}{\rho_i} \nabla \cdot \bar{\mathbf{p}}_i - \frac{1}{\rho_j} \nabla \cdot \bar{\mathbf{p}}_j + \frac{1}{n_e} \left( \frac{Z_i n_i}{\rho_i} - \frac{Z_j n_j}{\rho_j} \right) \left( \nabla \cdot \bar{\mathbf{p}}_e - \frac{1}{c} \mathbf{j} \times \mathbf{B} \right) \nonumber \\
& = - \mathbf{w} \cdot \frac{\rho}{\rho_i \rho_j} \left[ \rho_i \nu_{ij} + \rho_e \nu_{ei} \left( \frac{Z_j n_j}{n_e} \right)^2 + \rho_e \nu_{ej} \left( \frac{Z_i n_i}{n_e} \right)^2 \right] - \mathbf{j} \cdot \frac{\rho}{\rho_i \rho_j} \frac{m_e}{en_e} (\nu_{ei} Z_j n_j - \nu_{ej} Z_i n_i), \tag{E81}
\end{align}
$$

with the electric field

$$
\begin{equation}
\begin{aligned}
\mathbf{E} = & -\frac{1}{c} \mathbf{V} \times \mathbf{B} - \frac{1}{cn_e \rho} (Z_i n_i \rho_j - Z_j n_j \rho_i) \mathbf{w} \times \mathbf{B} + \frac{j \times B}{cen_e} - \frac{1}{en_e} \nabla \cdot \bar{\mathbf{p}}_e + \frac{m_e}{e^2 n_e} \mathbf{j}(\nu_{ei} + \nu_{ej}) \\
& + \frac{m_e}{en_e} \mathbf{w} (Z_j n_j \nu_{ei} - Z_i n_i \nu_{ej}).
\end{aligned}
\tag{E82}
\end{equation}
$$

The system (E80)–(E82) is of course equivalent to (E76)–(E79). However, in a particular case when the collisions are very frequent, the right-hand side of (E81) becomes very large, and neglecting all of the “inertial” terms in the first line of (E81) with $\mathbf{w}$ allows one to obtain an explicit expression for the velocity difference

$$
\mathbf{w} = \frac{1}{D} \left[ -j \frac{m_e}{en_e} (\nu_{ei} Z_j n_j - \nu_{ej} Z_i n_i) - \frac{\rho_j}{\rho} \nabla \cdot \bar{\mathbf{p}}_i + \frac{\rho_i}{\rho} \nabla \cdot \bar{\mathbf{p}}_j - \frac{1}{\rho n_e} (Z_i n_i \rho_j - Z_j n_j \rho_i) \left( \nabla \cdot \bar{\mathbf{p}}_e - \frac{1}{c} j \times \mathbf{B} \right) \right], \quad (\text{E83})
$$

where we define the denominator

$$
D = \left[ \rho_i \nu_{ij} + \rho_e \nu_{ei} \left( \frac{Z_j n_j}{n_e} \right)^2 + \rho_e \nu_{ej} \left( \frac{Z_i n_i}{n_e} \right)^2 \right]. \quad (\text{E84})
$$

For frequent collisions, only the first term in (E83) ~j is finite, and all of the other terms are small. Nevertheless, the sought-after term is the last term in (E83) ~j × B, because when (E83) is used in (E82) it creates the ambipolar term ~-j × B × B. The single-
fluid electric field reads

$$ 
\begin{aligned}
E = & -\frac{1}{c} \mathbf{V} \times \mathbf{B} + j \frac{m_e}{e^2 n_e} \left[ \nu_{ei} + \nu_{ej} - \frac{m_e}{n_e D} (\nu_{ei} Z_j n_j - \nu_{ej} Z_i n_i)^2 \right] \\
& + \frac{j \times \mathbf{B}}{c n_e} \left[ 1 + \frac{2m_e}{n_e \rho D} (Z_i n_i \rho_j - Z_j n_j \rho_i)(\nu_{ei} Z_j n_j - \nu_{ej} Z_i n_i) \right] \\
& - \frac{j \times \mathbf{B} \times \mathbf{B}}{c^2 n_e^2 \rho^2 D} (Z_i n_i \rho_j - Z_j n_j \rho_i)^2 \\
& - \frac{\nabla \cdot \bar{\mathbf{p}}_e}{e n_e} \left[ 1 + \frac{m_e}{\rho n_e D} (Z_j n_j \nu_{ei} - Z_i n_i \nu_{ej})(Z_i n_i \rho_j - Z_j n_j \rho_i) \right] \\
& + \frac{m_e}{e n_e \rho D} (Z_j n_j \nu_{ei} - Z_i n_i \nu_{ej}) \left[ -\rho_j \nabla \cdot \bar{\mathbf{p}}_i + \rho_i \nabla \cdot \bar{\mathbf{p}}_j \right] \\
& - \frac{1}{c n_e \rho^2 D} (Z_i n_i \rho_j - Z_j n_j \rho_i) \left[ -\rho_j \nabla \cdot \bar{\mathbf{p}}_i + \rho_i \nabla \cdot \bar{\mathbf{p}}_j - \frac{1}{n_e} (Z_i n_i \rho_j - Z_j n_j \rho_i) \nabla \cdot \bar{\mathbf{p}}_e \right] \times \mathbf{B}.
\end{aligned}
\quad (E85) $$

Importantly, the sign in front of the ambipolar term is negative, and because $-j \times B \times B = +j_{\perp}|B|^2$, the term indeed creates diffusion in the induction equation. It is possible to define a coefficient of ambipolar diffusion:

$$ \eta_A = \frac{|\mathbf{B}|^2}{4\pi\rho} A = V_A^2 A; \quad \text{where} \quad A = \frac{(Z_i n_i \rho_j - Z_j n_j \rho_i)^2}{n_e^2 \rho D}, \quad (E86) $$

and $V_A$ is the Alfvén speed. As a double check, prescribing zero charge for one of the species, the electric field (E85) identifies with Equation (A.10) of Zaqarashvili et al. (2011; for example, our denominator simplifies to $D = \alpha_{in} + \alpha_{en} = \alpha_n$). Also, $\eta_A = |\mathbf{B}|^2 \rho_n^2 / (4\pi\rho^2(\rho_i\nu_{in} + \rho_e\nu_{en}))$ identifies with the usual coefficient of ambipolar diffusion; see, for example, Equation (20) in Khomenko & Collados (2012; after switching to cgs units with $\mu_0 \to 4\pi$). The ambipolar diffusion exists if

$$ \frac{Z_i}{m_i} \neq \frac{Z_j}{m_j}. \quad (E87) $$

It is important to emphasize that the reduction to a single-fluid model was obtained by assuming that collisions are sufficiently frequent, and that the ambipolar diffusion (as well as other terms) now contains a denominator $D$, which can be simplified into $D = \rho_i\nu_{ij}$. So, when the collisional frequencies $\nu_{ij}$ become small, this leads to an artificial “explosion” of the ambipolar diffusion. This is nicely demonstrated in the figures of Zaqarashvili et al. (2011) plotted with respect to a wavenumber $\bar{k} \sim k/\nu$, where it is shown that for a single-fluid description, the ambipolar diffusion in a collisionless regime (when $\bar{k}$ becomes large) yields cutoff frequencies for waves. The mechanism is completely analogous to the “explosion” of the Braginskii stress tensor or the heat flux vector in a collisionless regime. In contrast, as they show in their two-fluid figures, no “explosion” of the ambipolar diffusion is present. The effect is further discussed in Zaqarashvili et al. (2012).

### E.6.1. Damping of Alfvén Waves

For example, considering Alfvén waves at long wavelengths, and focusing only on the ambipolar diffusion (with the Hall-term, ohmic terms, and pressure terms neglected), the induction equation reads

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{V} \times \mathbf{B}) + \nabla \times [\eta_A (\nabla \times \mathbf{B})_{\perp}], \quad (E88) $$

with the coefficient of ambipolar diffusion (E86). This yields the following dispersion relation for Alfvén waves:

$$ \omega^2 + i\omega V_A^2 k_{\parallel}^2 A - V_A^2 k_{\parallel}^2 = 0, \quad (E89) $$

with solutions

$$ \omega = \pm V_A |k_{\parallel}| \sqrt{1 - \left(\frac{V_A k_{\parallel} A}{2}\right)^2} - i \frac{V_A^2 k_{\parallel}^2 A}{2}. \quad (E90) $$

Obviously, the Alfvén waves are damped, and for wavenumbers $k_{\parallel} \ge 2/(V_A A)$ the real part of the frequency even becomes zero, so the wave stops existing (i.e., a cutoff wavenumber). For the particular case of one species being neutral, the quantity $A = \rho_n^2 / (\rho\alpha_n)$, which can be approximated as $A = \rho_n^2 / (\rho\alpha_{in}) = \zeta_n^2 / (\zeta_i\nu_{in})$. Then, expressions (E89), (E90) identify with Equations (44)–(47) of Zaqarashvili et al. (2011); however, one needs to use their definition $\nu_{in} = \alpha_{in}/\rho$ instead of the more logical (and correct) $\nu_{in} = \alpha_{in}/\rho_i$.