Now, by using (L36), (L42), we are ready to calculate the collisional integrals (L32). The first integral in (L32) calculates

$$
\begin{equation}
\begin{aligned}
& m_a \int f_a [(A^{ab} \mathbf{c}_a)^S |\mathbf{c}_a|^2 + 2(A^{ab} \cdot \mathbf{c}_a) \mathbf{c}_a \mathbf{c}_a] d^3v \\
&= \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \left\{ -\bar{\mathbf{I}} \frac{4(2T_a m_b + 5T_b m_a)}{(T_a m_b + T_b m_a)} - \bar{\mathbf{\tilde{h}}}^{(2)} \frac{6(3T_a m_b - 7T_b m_a)T_b m_a}{5(T_a m_b + T_b m_a)^2} \right. \\
& \qquad \left. + \hat{\mathbf{\tilde{h}}}^{(4)} \frac{3(T_a m_b - T_b m_a)T_b^2 m_a^2}{(T_a m_b + T_b m_a)^3} \right. \\
& \qquad \left. - \bar{\mathbf{\tilde{h}}}^{(2)} \frac{4(8T_a^2 m_b^2 + 28T_a T_b m_a m_b + 35T_b^2 m_a^2)}{5(T_a m_b + T_b m_a)^2} \right. \\
& \qquad \left. + \hat{\mathbf{\tilde{h}}}^{(4)} \frac{2(8T_a^3 m_b^3 + 36T_a^2 T_b m_a m_b^2 + 63T_a T_b^2 m_a^2 m_b - 70T_b^3 m_a^3)}{35(T_a m_b + T_b m_a)^3} \right\}.
\end{aligned}
\tag{L43}
\end{equation}
$$

The second integral in (L32) calculates

$$
\begin{equation}
\begin{aligned}
& m_a \int f_a [(\operatorname{Tr} \bar{\mathbf{D}}^{ab}) \mathbf{c}_a \mathbf{c}_a + \bar{\mathbf{D}}^{ab} |\mathbf{c}_a|^2 + 2((\bar{\mathbf{D}}^{ab} \cdot \mathbf{c}_a) \mathbf{c}_a)^S] d^3 v \\
&= \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \left\{ \bar{\mathbf{I}} \frac{4(2T_a m_b + 5T_b m_a)}{T_a(m_b + m_a)} + \bar{\mathbf{\tilde{h}}}^{(2)} \frac{2(11T_a m_b - 7T_b m_a)T_b m_a}{5T_a(T_a m_b + T_b m_a)(m_b + m_a)} \right. \\
& \qquad \left. - \hat{\mathbf{\tilde{h}}}^{(4)} \frac{3(23T_a m_b - 7T_b m_a)T_b^2 m_a^2}{35T_a(T_a m_b + T_b m_a)^2(m_b + m_a)} \right. \\
& \qquad \left. + \bar{\mathbf{\tilde{h}}}^{(2)} \frac{2(4T_a^2 m_b^2 + 21T_a T_b m_a m_b + 35T_b^2 m_a^2)}{5T_a(T_a m_b + T_b m_a)(m_b + m_a)} \right. \\
& \qquad \left. - \hat{\mathbf{\tilde{h}}}^{(4)} \frac{m_b(T_a m_b + 7T_b m_a)(4T_a m_b + 19T_b m_a)}{35(T_a m_b + T_b m_a)^2(m_b + m_a)} \right\}.
\end{aligned}
\tag{L44}
\end{equation}
$$

Adding (L43) and (L44) together then yields the collisional contributions

$$
\mathrm{Tr} \, \bar{\mathcal{Q}}_{ab}^{(4)} = \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \left\{ + \bar{\mathbf{I}} \frac{4(2T_a m_b + 5T_b m_a)m_a}{(T_a m_b + T_b m_a)(m_b + m_a)} \frac{T_b - T_a}{T_a} \right. \\
\left. - M_{ab(1)} \bar{\mathbf{\tilde{h}}}^{a(2)} + M_{ab(2)} \bar{\mathbf{\tilde{h}}}^{b(2)} - N_{ab(1)} \hat{\mathbf{\tilde{h}}}^{a(4)} - N_{ab(2)} \hat{\mathbf{\tilde{h}}}^{b(4)} \right\}, \quad \text{(L45)}
$$

with mass-ratio coefficients

$$
\begin{align}
M_{ab(1)} &= \{2(16T_a^3 m_a m_b^2 + 12T_a^3 m_b^3 + 56T_a^2 T_b m_a^2 m_b + 31T_a^2 T_b m_a m_b^2 + 70T_a T_b^2 m_a^3 \nonumber \\
&\quad +14T_a T_b^2 m_a^2 m_b - 35T_b^3 m_a^3)\}[5T_a(T_a m_b + T_b m_a)^2(m_b + m_a)]^{-1}; \nonumber \\
M_{ab(2)} &= -\frac{2T_b m_a (9T_a^2 m_a m_b - 2T_a^2 m_b^2 - 21T_a T_b m_a^2 - 25T_a T_b m_a m_b + 7T_b^2 m_a^2)}{5(T_a m_b + T_b m_a)^2 T_a (m_b + m_a)}; \nonumber \\
N_{ab(1)} &= \{16T_a^3 m_a m_b^3 + 12T_a^3 m_b^4 + 72T_a^2 T_b m_a^2 m_b^2 + 21T_a^2 T_b m_a m_b^3 + 126T_a T_b^2 m_a^3 m_b \nonumber \\
&\quad -54T_a T_b^2 m_a^2 m_b^2 - 140T_b^3 m_a^4 - 273T_b^3 m_a^3 m_b\}[35(T_a m_b + T_b m_a)^3(m_b + m_a)]^{-1}; \nonumber \\
N_{ab(2)} &= -\frac{3T_b^2 m_a^2 (35T_a^2 m_a m_b + 12T_a^2 m_b^2 - 35T_a T_b m_a^2 - 51T_a T_b m_a m_b + 7T_b^2 m_a^2)}{35(T_a m_b + T_b m_a)^3 T_a (m_b + m_a)}. \tag{L46}
\end{align}
$$

For a particular case of small temperature differences between species, the mass-ratio coefficients simplify into

$$
\begin{align}
M_{ab(1)} &= \frac{2(35m_a^2 + 35m_a m_b + 12m_b^2)}{5(m_b + m_a)^2}; \qquad M_{ab(2)} = \frac{4m_a(7m_a + m_b)}{5(m_b + m_a)^2}; \nonumber \\
N_{ab(1)} &= \frac{140m_a^3 + 7m_a^2 m_b - 25m_a m_b^2 - 12m_b^3}{35(m_b + m_a)^3}; \qquad N_{ab(2)} = \frac{12m_a^2(7m_a - 3m_b)}{35(m_b + m_a)^3}, \tag{L47}
\end{align}
$$

and for self-collisions, $M_{aa(1)} = 41/5$, $M_{aa(2)} = 8/5$, $N_{aa(1)} = 11/28$, and $N_{aa(2)} = 6/35$.
L.7. Exchange Rates $\tilde{Q}_a^{(4)}$,

Applying trace at (L45) yields the scalar

$$ \mathrm{TrTr} \bar{\mathcal{Q}}_{ab}^{(4)} = \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \left\{ +3 \frac{4(2T_a m_b + 5T_b m_a)m_a}{(T_a m_b + T_b m_a)(m_b + m_a)} \frac{T_b - T_a}{T_a} \right\}, \quad (L48) $$

and thus

$$ \begin{aligned} \bar{\mathcal{Q}}_{ab}^{(4)}' &= \mathrm{Tr} \bar{\mathcal{Q}}_{ab}^{(4)} - \frac{\bar{\mathcal{I}}}{3} \mathrm{Tr} \mathrm{Tr} \bar{\mathcal{Q}}_{ab}^{(4)} \\ &= \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} [-M_{ab(1)} \bar{\mathcal{h}}^{a(2)} + M_{ab(2)} \bar{\mathcal{h}}^{b(2)} - N_{ab(1)} \hat{\mathcal{h}}^{a(4)} - N_{ab(2)} \hat{\mathcal{h}}^{b(4)}]. \end{aligned} \quad (L49) $$

Finally, introducing summation over all of the “b” species, and rewritten with fluid moments,

$$ \begin{aligned} \bar{\mathcal{Q}}_{a}^{(4)}' = &-\frac{53}{20}\nu_{aa}\frac{p_a}{\rho_a}\bar{\mathcal{H}}_{a}^{(2)} - \frac{79}{140}\nu_{aa}\bar{\mathcal{H}}_{a}^{(4)} + \sum_{b \neq a}\nu_{ab}\left[-(M_{ab(1)} - 7N_{ab(1)})\frac{p_a}{\rho_a}\bar{\mathcal{H}}_{a}^{(2)}\right. \\ &\left.\quad+(M_{ab(2)} + 7N_{ab(2)})\frac{p_a^2}{\rho_a p_b}\bar{\mathcal{H}}_{b}^{(2)} - N_{ab(1)}\bar{\mathcal{H}}_{a}^{(4)} - N_{ab(2)}\frac{p_a^2}{p_b^2}\frac{\rho_b}{\rho_a}\bar{\mathcal{H}}_{b}^{(4)}\right]. \end{aligned} \quad (L50) $$

It is useful to define (introducing tilde)

$$ \hat{M}_{ab(1)} = M_{ab(1)} - 7N_{ab(1)}; \quad \hat{M}_{ab(2)} = M_{ab(2)} + 7N_{ab(2)}, \quad (L51) $$

and the final mass-ratio coefficients are given by (25).

## Appendix M
### Collisional Contributions for Scalar $\tilde{\mathcal{X}}^{(4)}$

Here we consider the perturbation

$$ \chi_b = \frac{1}{120} \tilde{h}^{b(4)} \tilde{H}^{b(4)}, \quad (M1) $$

with the Hermite polynomial $\tilde{H}^{(4)} = \tilde{c}^4 - 10\tilde{c}^2 + 15$ and the Hermite moment $\tilde{h}^{b(4)} = \frac{\rho_b}{p_b^2} \tilde{\mathcal{X}}^{b(4)}$. The Rosenbluth potentials become

$$ H_b(\mathbf{v}) = n_b \sqrt{\frac{m_b}{T_b}} \left\{ \frac{1}{\tilde{y}} \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) + \frac{1}{120} \tilde{h}^{b(4)} (3 - \tilde{y}^2) \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} \right\}; \quad (M2) $$

$$ G_b(\mathbf{v}) = n_b \sqrt{\frac{T_b}{m_b}} \left\{ \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} + \left(\tilde{y} + \frac{1}{\tilde{y}}\right) \operatorname{erf}\left(\frac{\tilde{y}}{\sqrt{2}}\right) - \frac{1}{60} \tilde{h}^{b(4)} \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} \right\}, \quad (M3) $$

and the dynamical friction vector and the diffusion tensor become

$$ \begin{aligned} A^{ab} = 2 \frac{c_{ab}}{m_a^2} \left(1 + \frac{m_a}{m_b}\right) \frac{n_b m_b}{T_b} \left\{ \tilde{y} \left( \sqrt{\frac{2}{\pi}} \frac{e^{-\tilde{y}^2/2}}{\tilde{y}^2} - \frac{\operatorname{erf}(\tilde{y}/\sqrt{2})}{\tilde{y}^3} \right) \right. \\ \left. - \tilde{y} \frac{\tilde{h}^{b(4)}}{120} \sqrt{\frac{2}{\pi}} (5 - y^2) e^{-\tilde{y}^2/2} \right\}, \end{aligned} \quad (M4) $$

$$ \begin{aligned} \bar{\mathcal{D}}^{ab} = 2 \frac{c_{ab}}{m_a^2} n_b \sqrt{\frac{m_b}{T_b}} \left\{ \bar{\mathcal{I}} \frac{\bar{\mathcal{A}}_1'}{\bar{\mathcal{I}}} \frac{\bar{\mathcal{A}}_1}{\bar{\mathcal{I}}} + \frac{\bar{\mathcal{y}} \bar{\mathcal{y}}}{\bar{\mathcal{y}}^2} \left( \bar{\mathcal{A}}_1'' - \frac{\bar{\mathcal{A}}_1'}{\bar{\mathcal{I}}} \right) \right. \\ \left. + (\bar{\mathcal{I}} - \bar{\mathcal{y}} \bar{\mathcal{y}}) \frac{\tilde{h}^{b(4)}}{60} \sqrt{\frac{2}{\pi}} e^{-\tilde{y}^2/2} \right\}. \end{aligned} \quad (M5) $$

The perturbation $\chi_a = (\tilde{h}^{a(4)}/120)(\alpha^4 \tilde{y}^4 - 10\alpha^2 \tilde{y}^2 + 15)$.
### M.1. Pressure Tensor Exchange Rates

It is sufficient to consider $c_a = \sqrt{T_a/m_a} \alpha \tilde{\mathbf{y}}$, and so

$$[A^{ab}c_a]^S = 2 \frac{c_{ab}}{m_a^2} \left(1 + \frac{m_a}{m_b}\right) \frac{n_b m_b}{T_b} \sqrt{\frac{T_a}{m_a}} \left\{ 2 \alpha \tilde{\mathbf{y}} \tilde{\mathbf{y}} \left( \sqrt{\frac{2}{\pi}} \frac{e^{-\tilde{\mathbf{y}}^2/2}}{\tilde{\mathbf{y}}^2} - \frac{\text{erf}(\tilde{\mathbf{y}}/\sqrt{2})}{\tilde{\mathbf{y}}^3} \right) \right. \\ \left. - 2 \alpha \tilde{\mathbf{y}} \tilde{\mathbf{y}} \frac{\tilde{\mathbf{h}}^{b(4)}}{120} \sqrt{\frac{2}{\pi}} (5 - y^2) e^{-\tilde{\mathbf{y}}^2/2} \right\}, \quad (M6)$$

which further integrates

$$m_a \int f_a [A^{ab} c_a]^S d^3v = \rho_a \nu_{ab} \bar{\mathbf{I}} \left[ -2 \frac{T_a}{m_a} - \tilde{\mathbf{h}}^{b(4)} \frac{T_a T_b^2 m_a}{4(T_a m_b + T_b m_a)^2} \right. \\ \left. - \tilde{\mathbf{h}}^{a(4)} \frac{m_b T_a^2 (T_a m_b - 4 T_b m_a)}{20 m_a (T_a m_b + T_b m_a)^2} \right], \quad (M7)$$

together with

$$m_a \int f_a \bar{\mathbf{D}}_{ab} d^3v = \frac{\rho_a \nu_{ab}}{m_a + m_b} \bar{\mathbf{I}} \left[ \frac{2}{m_a} (T_a m_b + T_b m_a) + \tilde{\mathbf{h}}^{b(4)} \frac{T_b^2 m_a}{20(T_a m_b + T_b m_a)} \right. \\ \left. + \tilde{\mathbf{h}}^{a(4)} \frac{m_b^2 T_a^2}{20 m_a (T_a m_b + T_b m_a)} \right]. \quad (M8)$$

Adding the last two results together yields the collisional contributions

$$\bar{\mathbf{Q}}_{ab}^{(2)} = \frac{\rho_a \nu_{ab}}{m_a + m_b} \bar{\mathbf{I}} \left[ +2(T_b - T_a) - T_b \tilde{\mathbf{h}}^{b(4)} \frac{T_b m_a (5T_a m_a + 4T_a m_b - T_b m_a)}{20(T_a m_b + T_b m_a)^2} \right. \\ \left. + T_a \tilde{\mathbf{h}}^{a(4)} \frac{T_a m_b (5T_b m_b + 4T_b m_a - T_a m_b)}{20(T_a m_b + T_b m_a)^2} \right], \quad (M9)$$

which can be written as

$$\bar{\mathbf{Q}}_{ab}^{(2)} = \frac{\rho_a \nu_{ab}}{m_a + m_b} \bar{\mathbf{I}} \left[ +2(T_b - T_a) + P_{ab(1)} T_a \tilde{\mathbf{h}}^{a(4)} - P_{ab(2)} T_b \tilde{\mathbf{h}}^{b(4)} \right], \quad (M10)$$

with mass-ratio coefficients

$$P_{ab(1)} = \frac{T_a m_b (5 T_b m_b + 4 T_b m_a - T_a m_b)}{20 (T_a m_b + T_b m_a)^2}; \quad P_{ab(2)} = \frac{T_b m_a (5 T_a m_a + 4 T_a m_b - T_b m_a)}{20 (T_a m_b + T_b m_a)^2}, \quad (M11)$$

or for the particular case of small temperature differences,

$$P_{ab(1)} = \frac{m_b}{5(m_b + m_a)}; \quad P_{ab(2)} = \frac{m_a}{5(m_b + m_a)}. \quad (M12)$$

The pressure tensor exchange rates (M10) are rewritten to fluid variables, according to

$$\bar{\mathbf{Q}}_{ab}^{(2)} = \frac{\rho_a \nu_{ab}}{m_a + m_b} \bar{\mathbf{I}} \left[ +2(T_b - T_a) + P_{ab(1)} \frac{\rho_a}{n_a p_a} \tilde{\mathbf{X}}_a^{(4)} - P_{ab(2)} \frac{\rho_b}{n_b p_b} \tilde{\mathbf{X}}_b^{(4)} \right]. \quad (M13)$$

The energy exchange rates then become

$$Q_{ab} = \frac{1}{2} \text{Tr} \bar{\mathbf{Q}}_{ab}^{(2)} = \frac{\rho_a \nu_{ab}}{(m_a + m_b)} \left[ +3(T_b - T_a) + \frac{3}{2} P_{ab(1)} \frac{\rho_a}{n_a p_a} \tilde{\mathbf{X}}_a^{(4)} - \frac{3}{2} P_{ab(2)} \frac{\rho_b}{n_b p_b} \tilde{\mathbf{X}}_b^{(4)} \right], \quad (M14)$$

and the collisional contributions for the stress tensor are

$$\bar{\mathbf{Q}}_{ab}^{(2)'} = \bar{\mathbf{Q}}_{ab}^{(2)} - \frac{\bar{\mathbf{I}}}{3} \text{Tr} \bar{\mathbf{Q}}_{ab}^{(2)} = 0. \quad (M15)$$

The scalar perturbations $\tilde{\mathbf{X}}_a^{(4)}$ and $\tilde{\mathbf{X}}_b^{(4)}$ thus do not modify the $\bar{\mathbf{Q}}_{ab}^{(2)'}$, but they enter the conservation of energy. The final model uses $\hat{P}_{ab(1)} = (3/2)P_{ab(1)}$ and $\hat{P}_{ab(2)} = (3/2)P_{ab(2)}$, and the result is written in Section 7.1, Equation (140). The result is also shown in the discussion, Equation (177).
## M.2. Fourth-order Moment Exchange Rates

It is straightforward to calculate

$$ (A^{ab}c_a)^S|c_a|^2 + 2(A^{ab} \cdot c_a)c_a c_a \\ = 2\frac{c_{ab}}{m_a^2}\left(1 + \frac{m_a}{m_b}\right)\frac{n_b m_b}{T_b}\left(\frac{T_a}{m_a}\right)^{3/2}\alpha^3\left\{4\tilde{y}\tilde{\tilde{y}}\left(\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^2/2} - \frac{\operatorname{erf}(\tilde{y}/\sqrt{2})}{\tilde{y}}\right) - 4\tilde{y}\tilde{\tilde{y}}\tilde{y}^2(5 - \tilde{y}^2)\frac{\tilde{h}^{b(4)}}{120}\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^2/2}\right\}, \quad (M16) $$

together with

$$ (\operatorname{Tr}\tilde{\bar{D}}^{ab})c_a c_a + \tilde{\bar{D}}^{ab}|c_a|^2 + 2[(\tilde{\bar{D}}^{ab} \cdot c_a)c_a]^S \\ = 2\frac{c_{ab}}{m_a^2}n_b\sqrt{\frac{m_b}{T_b}}\frac{T_a}{m_a}\alpha^2\left\{\tilde{\bar{I}}\tilde{\tilde{y}}\tilde{A}_1' + \tilde{\tilde{y}}\tilde{\tilde{y}}\left(6\tilde{A}_1'' + \frac{\tilde{A}_1'}{\tilde{\tilde{y}}}\right) + [\tilde{\bar{I}}\tilde{\tilde{y}}^2 + \tilde{\tilde{y}}\tilde{\tilde{y}}(7 - 6\tilde{y}^2)]\frac{\tilde{h}^{b(4)}}{60}\sqrt{\frac{2}{\pi}}e^{-\tilde{y}^2/2}\right\}, \quad (M17) $$

and to integrate

$$ m_a \int f_a [(A^{ab}c_a)^S|c_a|^2 + 2(A^{ab} \cdot c_a)c_a c_a] d^3v \\ = \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \tilde{\bar{I}} \left\{ -\frac{4(2T_a m_b + 5T_b m_a)}{(T_a m_b + T_b m_a)} + \tilde{h}^{b(4)} \frac{m_a^2 T_b^2 (2T_a m_b - 5T_b m_a)}{2(T_a m_b + T_b m_a)^3} + \tilde{h}^{a(4)} \frac{2T_a^3 m_b^3 + 9T_a^2 T_b m_a m_b^2 + 72T_a T_b^2 m_a^2 m_b - 40T_b^3 m_a^3}{30(T_a m_b + T_b m_a)^3} \right\}, \quad (M18) $$

together with

$$ m_a \int f_a [(\operatorname{Tr}\tilde{\bar{D}}^{ab})c_a c_a + \tilde{\bar{D}}^{ab}|c_a|^2 + 2((\tilde{\bar{D}}^{ab} \cdot c_a)c_a)^S] d^3v \\ = \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \tilde{\bar{I}} \left\{ \frac{4(2T_a m_b + 5T_b m_a)}{T_a(m_b + m_a)} - \tilde{h}^{b(4)} \frac{T_b^2 m_a^2 (2T_a m_b - T_b m_a)}{2T_a(T_a m_b + T_b m_a)^2 (m_b + m_a)} - \tilde{h}^{a(4)} \frac{m_b(2T_a^2 m_b^2 + T_a T_b m_a m_b + 44T_b^2 m_a^2)}{30(T_a m_b + T_b m_a)^2 (m_b + m_a)} \right\}. \quad (M19) $$

Adding the last two results together then yields the collisional contributions

$$ \operatorname{Tr} \tilde{\bar{Q}}_{ab}^{(4)} = \rho_a \nu_{ab} \frac{p_a^2}{\rho_a^2} \tilde{\bar{I}} \left\{ +S_{ab(0)} \frac{(T_b - T_a)}{T_a} - S_{ab(1)} \tilde{h}^{a(4)} - S_{ab(2)} \tilde{h}^{b(4)} \right\}, \quad (M20) $$

with mass-ratio coefficients

$$ S_{ab(0)} = \frac{4m_a(2T_a m_b + 5T_b m_a)}{(T_a m_b + T_b m_a)(m_b + m_a)}; \\ S_{ab(1)} = -\{m_a(2T_a^3 m_b^3 + 9T_a^2 T_b m_a m_b^2 + 6T_a^2 T_b m_b^3 + 72T_a T_b^2 m_a^2 m_b + 27T_a T_b^2 m_a m_b^2 - 40T_b^3 m_a^3 - 84T_b^3 m_a^2 m_b)\}[30(T_a m_b + T_b m_a)^3(m_b + m_a)]^{-1}; \\ S_{ab(2)} = -\frac{T_b^2 m_a^3 (2T_a^2 m_b - 5T_a T_b m_a - 6T_a T_b m_b + T_b^2 m_a)}{2T_a(T_a m_b + T_b m_a)^3(m_b + m_a)}. \quad (M21) $$

For the particular case of small temperature differences,

$$ S_{ab(0)} = \frac{4m_a(5m_a + 2m_b)}{(m_b + m_a)^2}; \\ S_{ab(1)} = \frac{2m_a(10m_a^2 - 7m_a m_b - 2m_b^2)}{15(m_b + m_a)^3}; \quad S_{ab(2)} = \frac{2m_a^3}{(m_b + m_a)^3}, \quad (M22) $$
and for self-collisions, $S_{aa(1)} = 1/60$ and $S_{aa(2)} = 1/4$. Applying a trace at (M20) and changing to fluid moments yields

$$ \mathrm{TrTr} \bar{\mathcal{Q}}_{ab}^{(4)} = 3\nu_{ab} \frac{p_a^2}{\rho_a} \left\{ +S_{ab(0)} \frac{(T_b - T_a)}{T_a} - S_{ab(1)} \frac{\rho_a}{p_a^2} \tilde{X}_a^{(4)} - S_{ab(2)} \frac{\rho_b}{p_b^2} \tilde{X}_b^{(4)} \right\}, \qquad (M23) $$

and the collisional contributions for the stress tensor $\bar{\Pi}_a^{(4)}$ are

$$ \bar{\mathcal{Q}}_{ab}^{(4)}' \equiv \mathrm{Tr} \bar{\mathcal{Q}}_{ab}^{(4)} - \frac{\bar{\mathcal{I}}}{3} \mathrm{Tr} \mathrm{Tr} \bar{\mathcal{Q}}_{ab}^{(4)} = 0. \qquad (M24) $$

### M.3. Collisional Contributions $\tilde{\mathcal{Q}}_a^{(4)}$,

The collisional contributions for the evolution equation $\tilde{X}_a^{(4)}$, Equation (C33), then become

$$ \begin{aligned} \tilde{\mathcal{Q}}_{ab}^{(4)}' &= \mathrm{TrTr} \bar{\mathcal{Q}}_{ab}^{(4)} - 20 \frac{p_a}{\rho_a} Q_{ab} \\ &= \nu_{ab} \left\{ + \frac{p_a^2 (T_b - T_a)}{T_a} \left( 3S_{ab(0)} - \frac{60m_a}{m_a + m_b} \right) - \tilde{X}_a^{(4)} \left( 3S_{ab(1)} + \frac{30m_a}{m_a + m_b} P_{ab(1)} \right) \right. \\ & \qquad \left. - \frac{p_a^2 \rho_b}{p_b^2 \rho_a} \tilde{X}_b^{(4)} \left( 3S_{ab(2)} - \frac{30m_a}{(m_a + m_b)} \frac{T_b}{T_a} P_{ab(2)} \right) \right\}. \end{aligned} \qquad (M25) $$

It is useful to define

$$ \begin{aligned} \hat{S}_{ab(0)} &= - \left( 3S_{ab(0)} - \frac{60m_a}{m_a + m_b} \right); \\ \hat{S}_{ab(1)} &= 3S_{ab(1)} + \frac{30m_a}{m_a + m_b} P_{ab(1)}; \\ \hat{S}_{ab(2)} &= - \left( 3S_{ab(2)} - \frac{30m_a}{(m_a + m_b)} \frac{T_b}{T_a} P_{ab(2)} \right), \end{aligned} \qquad (M26) $$

and the final model then reads

$$ \begin{aligned} \tilde{\mathcal{Q}}_{ab}^{(4)}' &= \mathrm{TrTr} \bar{\mathcal{Q}}_{ab}^{(4)} - 20 \frac{p_a}{\rho_a} Q_{ab} \\ &= \nu_{ab} \left\{ -\frac{p_a^2 (T_b - T_a)}{T_a} \hat{S}_{ab(0)} - \tilde{X}_a^{(4)} \hat{S}_{ab(1)} + \frac{p_a^2 \rho_b}{p_b^2 \rho_a} \tilde{X}_b^{(4)} \hat{S}_{ab(2)} \right\}, \end{aligned} \qquad (M27) $$

with mass-ratio coefficients given by (143).

## Appendix N
### Coupling of Two Species

Here we would like to emphasize the usefulness of the multifluid formulation, which makes the calculation of transport coefficients straightforward. We consider two species with indices “a” and “b.” The evolution equations for the heat fluxes “a” become

$$ \begin{aligned} \frac{d}{dt} \boldsymbol{q}_a + \Omega_a \hat{\boldsymbol{b}} \times \boldsymbol{q}_a + \frac{5}{2} p_a \nabla \left( \frac{p_a}{\rho_a} \right) &= -[2\nu_{aa} + \nu_{ab} \hat{D}_{ab(1)}] \boldsymbol{q}_a + \nu_{ab} \hat{D}_{ab(2)} \frac{\rho_a}{\rho_b} \boldsymbol{q}_b \\ &\quad + \left[ \frac{3}{70} \nu_{aa} + \nu_{ab} \hat{E}_{ab(1)} \right] \frac{\rho_a}{p_a} \boldsymbol{X}_a^{(5)} - \nu_{ab} \hat{E}_{ab(2)} \frac{\rho_a}{p_b} \boldsymbol{X}_b^{(5)} - p_a \nu_{ab} (\boldsymbol{u}_b - \boldsymbol{u}_a) \hat{U}_{ab(1)}; \end{aligned} \qquad (N1) $$

$$ \begin{aligned} \frac{d}{dt} \boldsymbol{X}_a^{(5)} + \Omega_a \hat{\boldsymbol{b}} \times \boldsymbol{X}_a^{(5)} + 70 \frac{p_a^2}{\rho_a} \nabla \left( \frac{p_a}{\rho_a} \right) &= - \left[ \frac{76}{5} \nu_{aa} + \nu_{ab} \hat{F}_{ab(1)} \right] \frac{p_a}{\rho_a} \boldsymbol{q}_a + \nu_{ab} \hat{F}_{ab(2)} \frac{p_a}{\rho_b} \boldsymbol{q}_b \\ &\quad - \left[ \frac{3}{35} \nu_{aa} + \nu_{ab} \hat{G}_{ab(1)} \right] \boldsymbol{X}_a^{(5)} - \nu_{ab} \hat{G}_{ab(2)} \frac{p_a}{p_b} \boldsymbol{X}_b^{(5)} - \frac{p_a^2}{\rho_a} \nu_{ab} (\boldsymbol{u}_b - \boldsymbol{u}_a) \hat{U}_{ab(2)}, \end{aligned} \qquad (N2) $$
together with the evolution equations for the heat fluxes “b”:

$$ \frac{d_b}{dt} \mathbf{q}_b + \Omega_b \hat{\mathbf{b}} \times \mathbf{q}_b + \frac{5}{2} p_b \nabla \left( \frac{p_b}{\rho_b} \right) = - \left[ 2\nu_{bb} + \nu_{ba} \hat{D}_{ba(1)} \right] \mathbf{q}_b + \nu_{ba} \hat{D}_{ba(2)} \frac{\rho_b}{p_a} \mathbf{q}_a \\ + \left[ \frac{3}{70} \nu_{bb} + \nu_{ba} \hat{E}_{ba(1)} \right] \frac{\rho_b}{p_b} \mathbf{X}_b^{(5)} - \nu_{ba} \hat{E}_{ba(2)} \frac{\rho_b}{p_a} \mathbf{X}_a^{(5)} + p_b \nu_{ba} (\mathbf{u}_b - \mathbf{u}_a) \hat{U}_{ba(1)}; \qquad (N3) $$

$$ \frac{d_b}{dt} \mathbf{X}_b^{(5)} + \Omega_b \hat{\mathbf{b}} \times \mathbf{X}_b^{(5)} + 70 \frac{p_b^2}{\rho_b} \nabla \left( \frac{p_b}{\rho_b} \right) = - \left[ \frac{76}{5} \nu_{bb} + \nu_{ba} \hat{F}_{ba(1)} \right] \frac{p_b}{\rho_b} \mathbf{q}_b + \nu_{ba} \hat{F}_{ba(2)} \frac{p_b}{\rho_a} \mathbf{q}_a \\ - \left[ \frac{3}{35} \nu_{bb} + \nu_{ba} \hat{G}_{ba(1)} \right] \mathbf{X}_b^{(5)} - \nu_{ba} \hat{G}_{ba(2)} \frac{p_b}{p_a} \mathbf{X}_a^{(5)} + \frac{p_b^2}{\rho_b} \nu_{ba} (\mathbf{u}_b - \mathbf{u}_a) \hat{U}_{ba(2)}, \qquad (N4) $$

where, for similar temperatures, the mass-ratio coefficients are given by (27), (28), and for arbitrary temperatures by (19), (21). The system is fully specified and after prescribing a quasistatic approximation it can be solved. Unfortunately, the general analytic solution is too long to write, even for the unmagnetized case. It is beneficial to consider a specific example. Nevertheless, the above system is a very powerful tool, which allows one to obtain the transport coefficients between two different species, it being a two ion plasma, or precise solutions for a specific one ion-electron plasma, without neglecting $m_e/m_i$ (see Section 8.8).

Similarly, the viscosity between two species is described by the evolution equations for the viscosity tensors of species “a”:

$$ \frac{d_a}{dt} \bar{\mathbf{\Pi}}_a^{(2)} + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(2)})^S + p_a \bar{\mathbf{W}}_a = - \frac{21}{10} \nu_{aa} \bar{\mathbf{\Pi}}_a^{(2)} + \frac{9}{70} \nu_{aa} \frac{\rho_a}{p_a} \bar{\mathbf{\Pi}}_a^{(4)} \\ + \frac{\rho_a}{m_a} \frac{\nu_{ab}}{m_b} \left[ - \hat{K}_{ab(1)} \frac{1}{n_a} \bar{\mathbf{\Pi}}_a^{(2)} + \hat{K}_{ab(2)} \frac{1}{n_b} \bar{\mathbf{\Pi}}_b^{(2)} + L_{ab(1)} \frac{\rho_a}{n_a p_a} \bar{\mathbf{\Pi}}_a^{(4)} - L_{ab(2)} \frac{\rho_b}{n_b p_b} \bar{\mathbf{\Pi}}_b^{(4)} \right]; \qquad (N5) $$

$$ \frac{d_a}{dt} \bar{\mathbf{\Pi}}_a^{(4)} + \Omega_a (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_a^{(4)})^S + 7 \frac{p_a^2}{\rho_a} \bar{\mathbf{W}}_a = - \frac{53}{20} \nu_{aa} \frac{p_a}{\rho_a} \bar{\mathbf{\Pi}}_a^{(2)} - \frac{79}{140} \nu_{aa} \bar{\mathbf{\Pi}}_a^{(4)} \\ + \nu_{ab} \left[ - \hat{M}_{ab(1)} \frac{p_a}{\rho_a} \bar{\mathbf{\Pi}}_a^{(2)} + \hat{M}_{ab(2)} \frac{p_a^2}{\rho_a p_b} \bar{\mathbf{\Pi}}_b^{(2)} - N_{ab(1)} \bar{\mathbf{\Pi}}_a^{(4)} - N_{ab(2)} \frac{p_b^2}{p_b^2} \bar{\mathbf{\Pi}}_b^{(4)} \right], \qquad (N6) $$

together with the evolution equations for the viscosity tensors of species “b”:

$$ \frac{d_b}{dt} \bar{\mathbf{\Pi}}_b^{(2)} + \Omega_b (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_b^{(2)})^S + p_b \bar{\mathbf{W}}_b = - \frac{21}{10} \nu_{bb} \bar{\mathbf{\Pi}}_b^{(2)} + \frac{9}{70} \nu_{bb} \frac{\rho_b}{p_b} \bar{\mathbf{\Pi}}_b^{(4)} \\ + \frac{\rho_b}{m_a + m_b} \frac{\nu_{ba}}{m_b} \left[ - \hat{K}_{ba(1)} \frac{1}{n_b} \bar{\mathbf{\Pi}}_b^{(2)} + \hat{K}_{ba(2)} \frac{1}{n_a} \bar{\mathbf{\Pi}}_a^{(2)} + L_{ba(1)} \frac{\rho_b}{n_b p_b} \bar{\mathbf{\Pi}}_b^{(4)} - L_{ba(2)} \frac{\rho_a}{n_a p_a} \bar{\mathbf{\Pi}}_a^{(4)} \right]; \qquad (N7) $$

$$ \frac{d_b}{dt} \bar{\mathbf{\Pi}}_b^{(4)} + \Omega_b (\hat{\mathbf{b}} \times \bar{\mathbf{\Pi}}_b^{(4)})^S + 7 \frac{p_b^2}{\rho_b} \bar{\mathbf{W}}_b = - \frac{53}{20} \nu_{bb} \frac{p_b}{\rho_b} \bar{\mathbf{\Pi}}_b^{(2)} - \frac{79}{140} \nu_{bb} \bar{\mathbf{\Pi}}_b^{(4)} \\ + \nu_{ba} \left[ - \hat{M}_{ba(1)} \frac{p_b}{\rho_b} \bar{\mathbf{\Pi}}_b^{(2)} + \hat{M}_{ba(2)} \frac{p_b^2}{\rho_b p_a} \bar{\mathbf{\Pi}}_a^{(2)} - N_{ba(1)} \bar{\mathbf{\Pi}}_b^{(4)} - N_{ba(2)} \frac{p_b^2}{p_b^2} \bar{\mathbf{\Pi}}_a^{(4)} \right]. \qquad (N8) $$

Here the heat fluxes (N1)–(N4) and viscosities (N5)–(N8) are decoupled, but one can consider more precise solutions with coupling between heat fluxes and viscosities, similar to Section 6.

## N.1. Protons and Alpha Particles (Unmagnetized)

As an example, we consider collisions between protons and alpha particles (fully ionized Helium with proton mass 4). The protons will be the “a” species and the alpha particles will be “b” species. For the ion coefficients, the collisions with electrons are neglected in an analogous fashion to Braginskii (1965). By prescribing mass $m_b = 4m_a$, the mass-ratio coefficients with equal temperatures
$T_a = T_b$ become

$$
\begin{align*}
\hat{D}_{ab(1)} &= \frac{499}{125}; & \hat{D}_{ab(2)} &= \frac{396}{125}; & \hat{E}_{ab(1)} &= \frac{87}{875}; & \hat{E}_{ab(2)} &= \frac{9}{175}; & \hat{U}_{ab(1)} &= \frac{6}{5}; \\
\hat{F}_{ab(1)} &= \frac{7624}{125}; & \hat{F}_{ab(2)} &= \frac{4848}{125}; & \hat{G}_{ab(1)} &= -\frac{171}{125}; & \hat{G}_{ab(2)} &= \frac{12}{25}; & \hat{U}_{ab(2)} &= 24, \\
\hat{D}_{ba(1)} &= \frac{2011}{500}; & \hat{D}_{ba(2)} &= \frac{117}{250}; & \hat{E}_{ba(1)} &= \frac{897}{14000}; & \hat{E}_{ba(2)} &= \frac{9}{700}; & \hat{U}_{ba(1)} &= \frac{3}{10}; \\
\hat{F}_{ba(1)} &= \frac{979}{50}; & \hat{F}_{ba(2)} &= \frac{1383}{125}; & \hat{G}_{ba(1)} &= \frac{8907}{7000}; & \hat{G}_{ba(2)} &= \frac{3}{10}; & \hat{U}_{ba(2)} &= \frac{39}{5}. \tag{N9}
\end{align*}
$$

By specifying the charges $Z_a = 1$, $Z_b = 2$, the four different collisional frequencies are related by

$$
\nu_{ba} = \frac{\rho_a}{\rho_b} \nu_{ab}; \quad \nu_{ab} = 8 \frac{n_b}{n_a} \sqrt{\frac{2}{5}} \nu_{aa}; \quad \nu_{bb} = 8 \frac{n_b}{n_a} \nu_{aa}, \qquad \text{(N10)}
$$

and we choose $\nu_{aa}$ as the reference frequency. Furthermore, applying the charge neutrality $n_a + 2n_b = n_e$, we choose as a reference the normalized density $N_a \equiv n_a/n_e$ and express $n_b/n_e = (1 - N_a)/2$. We also prescribe $\nabla T_a = \nabla T_b$.

Solving the system then yields the (parallel) thermal heat fluxes $\mathbf{q}_a^T = -\kappa_a \nabla T_a$, $\mathbf{q}_b^T = -\kappa_b \nabla T_b$, with thermal conductivities

$$
\kappa_a = \frac{T_a n_a}{m_a \nu_{aa}} \hat{\kappa}_a; \quad \kappa_b = \frac{T_b n_b}{m_b \nu_{bb}} \hat{\kappa}_b, \qquad \text{(N11)}
$$

and with normalized fully analytic values

$$
\hat{\kappa}_a = N_a \left\{ \left( -\frac{17989001}{10557600} \sqrt{10} + \frac{292708195}{54054912} \right) N_a^3 + \left( \frac{2129490299}{675686400} \sqrt{10} - \frac{1032644005}{108109824} \right) N_a^2 \right. \\
\left. + \left( -\frac{98252949}{45045760} \sqrt{10} + \frac{8035835}{1689216} \right) N_a + \frac{51625}{70384} \sqrt{10} + \frac{3425}{140768} \right\} / \Delta_1; \qquad \text{(N12)}
$$

$$
\hat{\kappa}_b = 32(1 - N_a) \left\{ \frac{125}{1024} + \left( \frac{128513167}{2162196480} \sqrt{10} - \frac{166007075}{864878592} \right) N_a^3 \right. \\
\left. + \left( -\frac{67953383}{540549120} \sqrt{10} + \frac{386788475}{864878592} \right) N_a^2 + \left( \frac{15671599}{216219648} \sqrt{10} - \frac{1540025}{4504576} \right) N_a \right\} / \Delta_1; \quad \text{(N13)}
$$

$$
\Delta_1 = \left[ 1 + \left( -\frac{722521001}{563072000} \sqrt{10} + \frac{14274588957}{3519200000} \right) N_a^4 + \left( \frac{1043512703}{337843200} \sqrt{10} - \frac{8606493541}{879800000} \right) N_a^3 \right. \\
\left. + \left( -\frac{23828129}{8798000} \sqrt{10} + \frac{15644893541}{1759600000} \right) N_a^2 + \left( \frac{23828129}{26394000} \sqrt{10} - 4 \right) N_a \right], \qquad \text{(N14)}
$$

or with numerical values

$$
\hat{\kappa}_a = N_a [2.3438 + 0.02684 N_a^3 + 0.4144 N_a^2 - 2.1404 N_a] / \Delta_1; \qquad \text{(N15)}
$$

$$
\hat{\kappa}_b = 32(1 - N_a)[0.1221 - 0.003988 N_a^3 + 0.04968 N_a^2 - 0.1127 N_a] / \Delta_1; \qquad \text{(N16)}
$$

$$
\Delta_1 = 1 - 0.001559 N_a^4 - 0.01485 N_a^3 + 0.3266 N_a^2 - 1.1451 N_a. \qquad \text{(N17)}
$$

Note that $n_a/\nu_{aa}$ is independent of $n_a$, and that is why the definitions (N11) were chosen. For the “b” species (the alpha particles), the results are written in a form so that it is easy to use $32n_b/(\nu_{bb}m_b) = n_a/(\nu_{aa}m_a)$. As a double check, prescribing

$$
\begin{align}
N_a = 1; \quad &\Rightarrow \quad \kappa_a = \frac{T_a n_a}{\nu_{aa} m_a} \frac{125}{32}; \quad \kappa_b = 0; \notag \\
N_a = 0; \quad &\Rightarrow \quad \kappa_a = 0; \quad \kappa_b = \frac{T_a n_b}{\nu_{bb} m_b} \frac{125}{32}, \tag{N18}
\end{align}
$$

as it should be. In general, the thermal conductivities of single-ion plasmas compare as $\kappa_a/\kappa_b = \sqrt{m_b/m_a}(Z_b/Z_a)^4$. In our case, the thermal conductivity of pure alpha particles is 32 times smaller than that of pure protons. The thermal conductivities $\hat{\kappa}_a, \hat{\kappa}_b$ are plotted in the left panel of Figure 7.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left: proton thermal conductivity $\hat{\kappa}_a$ (red), given by (N15), and alpha particle thermal conductivity $\hat{\kappa}_b$ (black), given by (N16). Middle: proton frictional heat flux given by $\beta_{0a}$ (N22). Right: alpha particle frictional heat flux given by $\beta_{0b}$ (N23). Note the surprising change of sign of $\beta_{0b}$ for $N_a > 0.57$. We have verified that the same effect is present in the simplified 13-moment model of Burgers (1969)–Schunk (1977).</figcaption></figure>

The frictional heat fluxes read

$$ q_a^u = -T_a n_e (\mathbf{u}_b - \mathbf{u}_a) \beta_{0a}; \quad q_b^u = -T_a n_e (\mathbf{u}_b - \mathbf{u}_a) \beta_{0b}; \qquad (N19) $$

$$ \beta_{0a} = N_a(1 - N_a) \left\{ \left( -\frac{150058601}{43990000} + \frac{1522393}{1407680} \sqrt{10} \right) N_a^3 + \left( +\frac{258658601}{43990000} - \frac{199422}{109975} \sqrt{10} \right) N_a^2 + \left( -\frac{16290}{4399} + \frac{99711}{109975} \sqrt{10} \right) N_a + \frac{5430}{4399} \right\} / \Delta_1; \qquad (N20) $$

$$ \beta_{0b} = N_a(1 - N_a) \left\{ \frac{7351}{1407680} \sqrt{10} + \left( -\frac{54551}{22522880} \sqrt{10} + \frac{264247}{35192000} \right) N_a^3 \right. \\ \left. + \left( \frac{289783}{22522880} \sqrt{10} - \frac{2663863}{70384000} \right) N_a^2 + \left( -\frac{22053}{1407680} \sqrt{10} + \frac{2663863}{140768000} \right) N_a \right\} / \Delta_1, \qquad (N21) $$

where the denominator $\Delta_1$ is identical to (N14), and with numerical values

$$ \beta_{0a} = N_a(1 - N_a)[1.2344 + 0.008776N_a^3 + 0.1457N_a^2 - 0.8360N_a] / \Delta_1; \qquad (N22) $$

$$ \beta_{0b} = N_a(1 - N_a)[0.01651 - 0.0001504N_a^3 + 0.002839N_a^2 - 0.03062N_a] / \Delta_1. \qquad (N23) $$

In both limits $N_a = 0, 1$, the frictional heat fluxes disappear. The frictional heat fluxes are plotted in the middle and right panels of Figure 7.

## N.2. Viscosities

One first calculates the required viscosity mass-ratio coefficients, which for the protons ("a") and alpha particles ("b") become

$$ \begin{aligned} \hat{K}_{ab(1)} &= \frac{398}{25}; & \hat{K}_{ab(2)} &= \frac{32}{25}; & L_{ab(1)} &= \frac{228}{175}; & L_{ab(2)} &= \frac{12}{175}; \\ \hat{M}_{ab(1)} &= \frac{934}{125}; & \hat{M}_{ab(2)} &= \frac{32}{125}; & N_{ab(1)} &= -\frac{8}{35}; & N_{ab(2)} &= -\frac{12}{875}; \\ \hat{K}_{ba(1)} &= \frac{323}{100}; & \hat{K}_{ba(2)} &= \frac{68}{25}; & L_{ba(1)} &= \frac{93}{700}; & L_{ba(2)} &= \frac{48}{175}; \\ \hat{M}_{ba(1)} &= -\frac{368}{125}; & \hat{M}_{ba(2)} &= \frac{1424}{125}; & N_{ba(1)} &= \frac{256}{125}; & N_{ba(2)} &= \frac{192}{175}, \end{aligned} \qquad (N24) $$

and which enter evolution Equations (N5)–(N8). For an unmagnetized plasma, a quasistatic solution of these equations then yields the viscosity tensors

$$ \begin{aligned} \bar{\Pi}_a^{(2)} &= -\frac{P_a}{\nu_{aa}} [\hat{\eta}_{aa} \bar{\bar{W}}_a + \hat{\eta}_{ab} \bar{\bar{W}}_b]; \\ \bar{\Pi}_b^{(2)} &= -\frac{P_b}{\nu_{bb}} [8 \hat{\eta}_{ab} \bar{\bar{W}}_a + \hat{\eta}_{bb} \bar{\bar{W}}_b], \end{aligned} \qquad (N25) $$
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Normalized viscosities of the proton and alpha particle plasma, according to (N26). Collisions with electrons are neglected, in an analogous fashion to Braginskii. Left: proton viscosity $\hat{\eta}_{aa}$ (red) and alpha particle viscosity $\hat{\eta}_{bb}$ (black). Right: “cross-viscosity” $\hat{\eta}_{ab}$.</figcaption></figure>

with numerical values

$$ 
\begin{aligned}
\hat{\eta}_{aa} &= N_a(-0.05464N_a^3 + 0.3704N_a^2 - 0.7717N_a + 0.5173)/\Delta; \\
\hat{\eta}_{ab} &= N_a(1 - N_a)(0.001874N_a^2 - 0.008142N_a + 0.01248)/\Delta; \\
\hat{\eta}_{bb} &= 8(1 - N_a)(-0.01150N_a^3 + 0.07862N_a^2 - 0.1729N_a + 0.11997)/\Delta; \\
\Delta &= 1 + 0.03923N_a^4 - 0.3759N_a^3 + 1.2959N_a^2 - 1.8953N_a.
\end{aligned}
\qquad (N26) $$

Note that $p_a/\nu_{aa} = 8p_b/\nu_{bb}$, and the chosen form (N25) emphasizes that the “cross-viscosities” $\hat{\eta}_{ab}$ are directly related. In general, the viscosities of a pure single ion species compare as $\eta_a/\eta_b = \sqrt{m_a/m_b}(Z_b/Z_a)^4$, so in our case the viscosity of the pure alpha particles is eight times smaller than that of the pure protons. We provide only numerical values for the solutions (N26), nevertheless it can be shown that for $N_a = 1$, the proton viscosity $\hat{\eta}_{aa} = 1025/1068$, and the same result is obtained for the alpha particle viscosity $\hat{\eta}_{bb}$ if $N_a = 0$. The “cross-viscosity” $\hat{\eta}_{ab}$ becomes zero for both $N_a = 1$ and $N_a = 0$. The results are plotted in Figure 8.

### N.3. Deuterium and Tritium Plasma (Unmagnetized)

Here we calculate another example of deuterium-tritium plasma, also considered by Simakov & Molvig (2016b). Plasma consisting of deuterium-tritium is probably the most efficient way of achieving plasma fusion. It is, for example, being used in the JET machine (see e.g., Joffrin et al. 2019), and it will be used in ITER.<fnref n="9" /> Of course, we do not consider the peculiar complications associated with the neoclassical toroidal geometry, as our calculation is classical. The deuterium core consists of one proton and one neutron. The tritium core consists of one proton and two neutrons. Deuterium will be the “a” species, and tritium will be the “b” species. The collisions with electrons are neglected. By prescribing $m_b = (3/2)m_a$, the mass-ratio coefficients with equal temperatures $T_b = T_a$ become

$$ 
\begin{aligned}
\hat{D}_{ab(1)} &= \frac{1989}{500}; & \hat{D}_{ab(2)} &= \frac{324}{125}; & \hat{E}_{ab(1)} &= \frac{189}{2000}; & \hat{E}_{ab(2)} &= \frac{81}{1400}; & \hat{U}_{ab(1)} &= \frac{9}{10}; \\
\hat{F}_{ab(1)} &= \frac{13543}{250}; & \hat{F}_{ab(2)} &= \frac{5022}{125}; & \hat{G}_{ab(1)} &= -\frac{1373}{1400}; & \hat{G}_{ab(2)} &= \frac{81}{100}; & \hat{U}_{ab(2)} &= \frac{99}{5}, \\
\hat{D}_{ba(1)} &= \frac{521}{125}; & \hat{D}_{ba(2)} &= \frac{189}{125}; & \hat{E}_{ba(1)} &= \frac{78}{875}; & \hat{E}_{ba(2)} &= \frac{27}{700}; & \hat{U}_{ba(1)} &= \frac{3}{5}; \\
\hat{F}_{ba(1)} &= \frac{5832}{125}; & \hat{F}_{ba(2)} &= \frac{3672}{125}; & \hat{G}_{ba(1)} &= -\frac{307}{875}; & \hat{G}_{ba(2)} &= \frac{18}{25}; & \hat{U}_{ba(2)} &= \frac{72}{5}.
\end{aligned}
\qquad (N27) $$

Further specifying $Z_a = Z_b = 1$, the collisional frequencies are related by

$$ \nu_{ab} = \frac{n_b}{n_a} \sqrt{\frac{6}{5}} \nu_{aa}; \qquad \nu_{bb} = \frac{n_b}{n_a} \sqrt{\frac{2}{3}} \nu_{aa}, \qquad (N28) $$

and the charge neutrality $n_a + n_b = n_e$ implies $n_b/n_e = 1 - N_a$, where $N_a = n_a/n_e$. These mass-ratio coefficients and collisional frequencies are used in the system (N1)–(N4). We present quasistatic solutions only for the unmagnetized case, and we assume $\nabla T_a = \nabla T_b$. The thermal heat fluxes $\mathbf{q}_a^T = -\kappa_a \nabla T_a; \mathbf{q}_b^T = -\kappa_b \nabla T_b$ are given by

$$ \kappa_a = \frac{T_a n_a}{m_a \nu_{aa}} \hat{\kappa}_a; \qquad \kappa_b = \frac{T_a n_b}{m_b \nu_{bb}} \hat{\kappa}_b, \qquad (N29) $$

<footnote n="9">⁹ www.iter.org/sci/FusionFuels</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Left: deuterium thermal conductivity $\hat{\kappa}_a$ (red) and tritium thermal conductivity $\hat{\kappa}_b$ (black), given by (N30). Middle: deuterium frictional heat flux, given by $\beta_{0a}$ (N32). Right: tritium frictional heat flux given by $\beta_{0b}$ (N33). Note that the frictional heat fluxes $q_b^u$ are defined with opposite signs in (N31) and (N19).</figcaption></figure>

and with numerical values

$$ 
\begin{aligned}
\hat{\kappa}_a &= N_a(4.2135 - 0.009780N_a^3 + 0.06292N_a^2 + 1.4992N_a) / \Delta; \\
\hat{\kappa}_b &= \sqrt{3/2}(1 - N_a)(3.1894 - 0.001385N_a^3 + 0.04936N_a^2 + 0.9845N_a) / \Delta; \\
\Delta &= 1 - 0.0021475N_a^4 - 0.01543N_a^3 + 0.01753N_a^2 + 0.4761N_a,
\end{aligned} 
\quad (N30) $$

where one can also use $\sqrt{3/2}n_b/(m_b\nu_{bb}) = n_a/(\nu_{aa}m_a)$. The frictional heat fluxes are given by

$$ q_a^u = -T_a n_e (\mathbf{u}_b - \mathbf{u}_a) \beta_{0a}; \quad q_b^u = +T_a n_e (\mathbf{u}_b - \mathbf{u}_a) \beta_{0b}; \quad (N31) $$

$$ \beta_{0a} = N_a(1 - N_a)[0.81156 + 0.010099N_a^3 + 0.098815N_a^2 + 0.50235N_a] / \Delta; \quad (N32) $$

$$ \beta_{0b} = N_a(1 - N_a)[0.26178 + 0.0088461N_a^3 + 0.069351N_a^2 + 0.24742N_a] / \Delta. \quad (N33) $$

## N.4. Viscosities

The required viscosity mass-ratio coefficients for deuterium (“a”) and tritium (“b”) become

$$ 
\begin{aligned}
\hat{K}_{ab(1)} &= \frac{397}{50}; & \hat{K}_{ab(2)} &= \frac{44}{25}; & L_{ab(1)} &= \frac{207}{350}; & L_{ab(2)} &= \frac{24}{175}; \\
\hat{M}_{ab(1)} &= \frac{166}{25}; & \hat{M}_{ab(2)} &= \frac{184}{125}; & N_{ab(1)} &= \frac{86}{875}; & N_{ab(2)} &= \frac{48}{875}; \\
\hat{K}_{ba(1)} &= \frac{124}{25}; & \hat{K}_{ba(2)} &= \frac{56}{25}; & L_{ba(1)} &= \frac{54}{175}; & L_{ba(2)} &= \frac{36}{175}; \\
\hat{M}_{ba(1)} &= \frac{444}{125}; & \hat{M}_{ba(2)} &= \frac{24}{5}; & N_{ba(1)} &= \frac{702}{875}; & N_{ba(2)} &= \frac{324}{875},
\end{aligned} 
\quad (N34) $$

and enter evolution Equations (N5)–(N8). For an unmagnetized plasma, the solutions read

$$ 
\begin{aligned}
\bar{\bar{\Pi}}_a^{(2)} &= -\frac{P_a}{\nu_{aa}}[\hat{\eta}_{aa}\bar{\bar{W}}_a + \hat{\eta}_{ab}\bar{\bar{W}}_b]; \\
\bar{\bar{\Pi}}_b^{(2)} &= -\frac{P_b}{\nu_{bb}}\left[\sqrt{\frac{2}{3}}\hat{\eta}_{ab}\bar{\bar{W}}_a + \hat{\eta}_{bb}\bar{\bar{W}}_b\right],
\end{aligned} 
\quad (N35) $$

with numerical values

$$ 
\begin{aligned}
\hat{\eta}_{aa} &= N_a(0.0046589N_a^3 + 0.0064481N_a^2 + 0.17316N_a + 0.85048) / \Delta; \\
\hat{\eta}_{ab} &= N_a(1 - N_a)(0.0049729N_a^2 + 0.028578N_a + 0.16621) / \Delta; \\
\hat{\eta}_{bb} &= \sqrt{2/3}(1 - N_a)(-0.0057061N_a^3 - 0.047294N_a^2 - 0.10519N_a + 1.17543) / \Delta; \\
\Delta &= 1 + 0.00017711N_a^4 - 0.00044516N_a^3 - 0.020987N_a^2 + 0.099409N_a.
\end{aligned} 
\quad (N36) $$
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Viscosities of deuterium and tritium plasma, according to (N36). Left: deuterium viscosity $\hat{\eta}_{aa}$ (red) and tritium viscosity $\hat{\eta}_{bb}$ (black). Right: “cross-viscosity” $\hat{\eta}_{ab}$.</figcaption></figure>

The solutions are written in a form so that one can directly use $\sqrt{2/3}p_b/\nu_{bb} = p_a/\nu_{aa}$, and are plotted in Figure 10. To obtain more precise solutions, one should include the collisions with electrons (i.e., consider coupling between three species). Nevertheless, the self-collisional values 1025/1068 = 0.96 will only change to roughly 0.89 (see, for example, Equation (217)), and the plotted viscosity profiles will not change much.

## ORCID iDs

P. Hunana https://orcid.org/0000-0002-9860-9759

T. Passot https://orcid.org/0000-0002-6845-9729

E. Khomenko https://orcid.org/0000-0003-3812-620X

D. Martínez-Gómez https://orcid.org/0000-0002-2653-7020

M. Collados https://orcid.org/0000-0002-6210-9648

A. Tenerani https://orcid.org/0000-0003-2880-6084

G. P. Zank https://orcid.org/0000-0002-4642-6192

Y. Maneva https://orcid.org/0000-0002-2422-6050

M. L. Goldstein https://orcid.org/0000-0002-5317-988X

G. M. Webb https://orcid.org/0000-0002-0617-9502

## References

Bale, S. D., Pulupa, M., Salem, C., Chen, C. H. K., & Quataert, E. 2013, ApJ, 769, L22

Balescu, R. 1988, Transport Processes in Plasmas, Vol. 1 and 2 (Amsterdam: North-Holland)

Barakat, A. R., & Schunk, R. W. 1981, JPhD, 14, 421

Barakat, A. R., & Schunk, R. W. 1982, PPh, 24, 389

Bhatnagar, P. L., Gross, E. P., & Krook, M. 1954, PhRv, 94, 511

Braginskii, S. I. 1958, JETP, 6, 358

Braginskii, S. I. 1965, RvPP, 1, 205

Burgers, J. M. 1969, Flow Equations for Composite Gases (New York: Academic)

Canullo, M. V., Costa, A., & Ferro-Fontan, C. 1996, ApJ, 462, 1005

Chapman, S., & Cowling, T. G. 1939, The Mathematical Theory of Non-Uniform Gases (Cambridge: Cambridge Univ. Press)

Chen, J. G., Xu, X. Q., & Lei, Y. A. 2019, CoPhC, 236, 128

Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, RSPSA, 236, 112

Chodura, R., & Pohl, F. 1971, PPh, 13, 645

Cohen, R. S., Spitzer, L., & Routly, P. M. 1950, PhRv, 80, 230

Cranmer, S. R., & Schiff, A. J. 2021, JGRA, 126, e29666

Cuperman, S., & Dryer, M. 1985, ApJ, 298, 414

Cuperman, S., Weiss, I., & Dryer, M. 1980, ApJ, 239, 345

Cuperman, S., Weiss, I., & Dryer, M. 1981, ApJ, 251, 297

Demars, H. G., & Schunk, R. W. 1979, JPhD, 12, 1051

Dreicer, H. 1959, PhRv, 115, 238

Epperlein, E. M., & Haines, M. G. 1986, PPhI, 29, 1029

Fitzpatrick, R. 2015, Plasma Physics. An Introduction (London, Boca Raton, FL: Taylor and Francis, CRC Press)

Gary, S. P. 1993, Theory of Space Plasma Microinstabilities (New York: Cambridge Univ. Press)

Goswami, P., Passot, T., & Sulem, P. L. 2005, PhPl, 12, 102109

Grad, H. 1949a, CPAM, 2, 325

Grad, H. 1949b, CPAM, 2, 331

Grad, H. 1958, HDP, 3, 205

Gross, E. P., & Krook, M. 1956, PhRv, 102, 593

Gurnett, D. A., & Bhattacharjee, A. 2005, Introduction to Plasma Physics: With Space and Laboratory Applications (Cambridge: Cambridge Univ. Press)

Haack, J. R., Hauck, C. D., & Murillo, M. S. 2017, JSP, 168, 826

Halekas, J. S., Whittlesey, P. L., Larson, D. E., et al. 2021, A&A, 650, A15

Hammett, G., Dorland, W., & Perkins, F. 1992, PhFIB, 4, 2052

Hammett, G. W., & Perkins, W. F. 1990, PhRvL, 64, 3019

Helander, P., & Sigmar, D. J. 2002, Collisional Transport in Magnetized Plasmas (Cambridge: Cambridge Univ. Press)

Hellinger, P., & Trávníček, P. M. 2009, PhPl, 16, 054501

Hinton, F. L. 1983, in Handbook of Plasma Physics, Vol. 1: Basic Plasma Physics 1, ed. M. N. Rosenbluth & R. Z. Sagdeev (Amsterdam: North-Holland), 147

Hirshman, S., & Sigmar, D. 1981, NucFu, 21, 1079

Hirshman, S. P., & Sigmar, D. J. 1977, PhFl, 20, 418

Hollweg, J. V. 1974, JGR, 79, 3845

Hollweg, J. V. 1976, JGR, 81, 1649

Hunana, P., Tenerani, A., Zank, G. P., et al. 2019a, JPIPh, 85, 205850603

Hunana, P., Tenerani, A., Zank, G. P., et al. 2019b, JPIPh, 85, 205850602

Hunana, P., & Zank, G. P. 2017, ApJ, 839, 13

Hunana, P., Zank, G. P., Laurenza, M., et al. 2018, PhRvL, 121, 135101

Ji, J.-Y., & Held, E. D. 2006, PhPl, 13, 102103

Ji, J.-Y., & Held, E. D. 2013, PhPl, 20, 042114

Ji, J.-Y., & Held, E. D. 2015, PhPl, 22, 062114

Ji, J.-Y., & Joseph, I. 2018, PhPl, 25, 032117

Joffrin, E., Abduallev, S., Abhangi, M., et al. 2019, NucFu, 59, 112021

Joseph, I., & Dimits, A. M. 2016, CoPP, 56, 504

Kaneko, B. S. 1960, JPSJ, 15, 1685

Kaneko, S., & Taguchi, M. 1978, JPSJ, 45, 1380

Kaneko, S., & Yamao, A. 1980, JPSJ, 48, 2098

Kaufman, A. N. 1960, PhFl, 3, 610

Khomenko, E., & Collados, M. 2012, ApJ, 747, 87

Khomenko, E., Collados, M., Díaz, A., & Vitas, N. 2014, PhPl, 21, 092901

Killie, M. A., Janse, A. M., Lie-Svendsen, O., & Leer, E. 2004, ApJ, 604, 842

Krommes, J. A. 2018a, JPIPh, 84, 925840401

Krommes, J. A. 2018b, JPIPh, 84, 905840601

Kunz, M. W. 2021, Lecture Notes on Irreversible Processes in Plasmas, www.astro.princeton.edu/~kunz/Site/AST554/AST554_lecture_notes_Kunz.pdf

Lancia, L., Albertazzi, B., Boniface, C., et al. 2014, PhRvL, 113, 235001

Landau, L. 1936, Phys. Z. Sowjetunion, 10, 154

Landau, L. 1937, J. Exptl. Theoret. Phys. U.S.S.R., 7, 203

Landi, S., Matteini, L., & Pantellini, F. 2014, ApJL, 790, L12

Landshoff, R. 1949, PhRv, 76, 904

Landshoff, R. 1951, PhRv, 82, 442
Luciani, J. F., Mora, P., & Virmont, J. 1983, PhRvL, 51, 1664

Macmahon, A. 1965, PhFl, 8, 1840

Meyer, E., & Sessler, G. 1957, ZPhy, 149, 15

Meyrand, R., Kanekar, A., Dorland, W., & Schekochihin, A. A. 2019, PNAS, 116, 1185

Mikhailovskii, A. B., & Smolyakov, A. I. 1985, Sov. Phys. JETP, 61, 109

Oraevskii, V., Chodura, R., & Feneberg, W. 1968, PPh, 10, 819

Passot, T., & Sulem, P. L. 2007, PhPl, 14, 082502

Passot, T., Sulem, P. L., & Hunana, P. 2012, PhPl, 19, 082113

Perrone, D., Passot, T., Laveder, D., et al. 2018, PhPl, 25, 052302

Ramos, J. J. 2005, PhPl, 12, 052102

Salem, C., Hubert, D., Lacombe, C., et al. 2003, ApJ, 585, 1147

Schunk, R. W. 1975, P&SS, 23, 437

Schunk, R. W. 1977, RvGSP, 15, 429

Schurtz, G. P., Nicolaï, P. D., & Busquet, M. 2000, PhPl, 7, 4238

Scudder, J. D. 2021, ApJ, 907, 90

Simakov, A. N., & Molvig, K. 2014, PhPl, 21, 024503

Simakov, A. N., & Molvig, K. 2016a, PhPl, 23, 032115

Simakov, A. N., & Molvig, K. 2016b, PhPl, 23, 032116

Snyder, P. B., & Hammett, G. W. 2001, PhPl, 8, 3199

Snyder, P. B., Hammett, G. W., & Dorland, W. 1997, PhPl, 4, 3974

Spitzer, L. 1962, Physics of Fully Ionized Gases (2nd ed.; New York: Wiley)

Spitzer, L., & Härm, R. 1953, PhRv, 89, 977

Stubbe, P. 1994, PhFl, 6, 3164

Stubbe, P., & Sukhorukov, A. I. 1999, PhPl, 6, 2976

Sulem, P. L., & Passot, T. 2015, JPIPh, 81, 325810103

Tanenbaum, B. S. 1967, Plasma Physics (New York: McGraw-Hill)

Verscharen, D., Klein, K. G., & Maruca, B. A. 2019, LRSP, 16, 5

Wang, L., Zhu, B., Xu, X.-q., & Li, B. 2019, AIPA, 9, 015217

Zank, G. P. 2014, Transport Processes in Space Physics and Astrophysics
(Berlin: Springer)

Zaqarashvili, T. V., Carbonell, M., Ballester, J. L., & Khodachenko, M. L. 2012,
A&A, 544, A143

Zaqarashvili, T. V., Khodachenko, M. L., & Rucker, H. O. 2011, A&A,
529, A82

Zhdanov, V. M. 2002, PPCF, 44, 2283