1 Introduction
Modelling of gas dynamics has been attracting people’s attention for centuries. Even for the simplest single-species, monatomic gas, while the classical continuum models such as the Euler equations and the Navier–Stokes–Fourier equations work well in most circumstances, people do find them inadequate when we care about some ‘extreme cases’, such as the low-density regime and a large velocity slip or temperature jump at a solid wall due to gas–surface interaction. In these cases, the interaction between gas molecules are either insufficient or severely ruined by gas–surface interactions, resulting in the failure of continuum models. Although some microscopic models, such as the Boltzmann equation, Enskog equation or even molecular dynamics, have been validated to be accurate for most applications, they are usually too expensive to solve due to the high dimensionality. Despite the fast computers developed nowadays, a full three-dimensional simulation of the Boltzmann equation still requires a huge amount of computational resources (Dimarco et al. Reference Dimarco, Loubére, Narski and Rey2018), and therefore lower-dimensional models are preferable for slip or early transitional flows. Since there is a large gap between the continuum models and the microscopic models, researchers have been trying to find models sitting in between, which are cheaper to simulate than the kinetic models.
Since the Euler equations and the Navier–Stokes equations can be considered as zeroth-order and first-order approximations of the Boltzmann equation in the continuum limit (Struchtrup Reference Struchtrup2005a), various attempts have been made to derive higher-order approximations. For example, by Chapman–Enskog expansion (Chapman Reference Chapman1916; Enskog Reference Enskog1921), one obtains the Burnett equations and the super-Burnett equations as third- and fourth-order approximations (Burnett Reference Burnett1936; Reinecke & Kremer Reference Reinecke and Kremer1990; Shavaliyev Reference Shavaliyev1993); by Grad’s expansion, equations for the stress tensor and heat fluxes can be derived to provide better closure than the Navier–Stokes and Fourier laws so that the models are suitable for a wider range of Knudsen numbers (Grad Reference Grad1949, Reference Grad1958); by the assumption of maximum entropy, the Euler equations can be extended to include 14 (or more) moments (Dreyer Reference Dreyer1987; Levermore Reference Levermore1996; Müller & Ruggeri Reference Müller and Ruggeri1998). However, these attempts show that going beyond Navier–Stokes is quite non-trivial: the Burnett and super-Burnett equations are linearly unstable (Bobylev Reference Bobylev1982); Grad’s method has the hyperbolicity problem and the convergence problem (Müller & Ruggeri Reference Müller and Ruggeri1998; Cai, Fan & Li Reference Cai, Fan and Li2014); equations by maximum entropy are still difficult to solve numerically due to an ill-posed optimization problem hidden in the equations (Tallec & Perlat Reference Tallec and Perlat1997). These deficiencies have been severely restricting the applications of these models.
In spite of this, a number of new thoughts have been introduced for this classical modelling problem. In the current century, all these classical models have been re-studied. The Burnett equations have been fixed to regain linear stability (Bobylev Reference Bobylev2006, Reference Bobylev2008); the hyperbolicity problem of Grad’s equations is fixed in Cai, Fan & Li (Reference Cai, Fan and Li2015); approximations of the maximum-entropy equations have been proposed which have explicit analytical expressions (McDonald & Torrilhon Reference McDonald and Torrilhon2013). At the same time, Grad’s old idea hidden in his notes (Grad Reference Grad1958) has been picked up to build new models, called regularized moment equations (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003). Although all these models are quite new, based on current studies, we find the regularized moment equations to be interesting due to the relatively complete theory (boundary conditions (Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2008), H-theorem (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2007; Torrilhon Reference Torrilhon2012)) and a number of numerical studies (Torrilhon Reference Torrilhon2006). However, the complete regularized 13 (R13)-moment equations have been derived only for Maxwell molecules. In Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2004), it has been demonstrated by the example of plane shock structure that the equations derived for Maxwell molecules are not directly applicable to the hard-sphere model, which indicates that collision models need to be taken into account during the model derivation. This inspires us to go beyond Maxwell molecules, and study more realistic interaction models directly. In this work, we focus on inverse power potentials, which cover both Maxwell molecules and hard-sphere molecules (as the limit), and have been verified by experiments to be realistic for a number of gases (Torrens Reference Torrens1972, Table 8.1).
As far as we know, the only regularized moment model derived for non-Maxwell monatomic molecules is by Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013), which consists of the fully linearized equations for hard-sphere molecules. Such a model cannot be applied in nonlinear regimes such as shock waves. In this work, we extend this work (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2013) and write down equations for all inverse power law (IPL) models linearized about the local Maxwellian. The long derivation of the R13 equations is done by our automated Mathematica code. In Gupta & Torrilhon (Reference Gupta and Torrilhon2012), the authors already used computer algebra systems to derive complicated moment equations, which turns out to be much more efficient than using pen and paper. Plane shock wave structures will be computed based on these 13-moment models to show better results than a simple alteration of the Maxwell model.
The rest of this paper is organized as follows. In § 2, we first introduce the explicit expressions for the R13-moment equations, and then show the linear stability analysis. In § 3, the derivation of the R13-moment equations is presented. Some numerical experiments verifying the capability of the R13 system are carried out in § 4 and some concluding remarks are made in § 5. A brief introduction to the Boltzmann equation, the expressions for the infinite moment equations and the concrete form of the right-hand side of the R13 equations are given in the appendices.
2 R13-moment equations for linearized IPL model
In this section, we are going to present the R13-moment equations for the IPL model, followed by their linear stability and dispersion relations. Before that, we start with a quick review of some properties of the IPL model.
2.1 A brief review of the IPL model
The IPL model contains a class of potentials that are frequently studied in the gas kinetic theory (cf. Harris Reference Harris1971; Chapman & Cowling Reference Chapman and Cowling1990; Bird Reference Bird1994; Struchtrup Reference Struchtrup2005b). It assumes that the potential between two molecules is proportional to an inverse power of the distance between them
$$\begin{eqnarray}\unicode[STIX]{x1D711}(r)=\frac{\unicode[STIX]{x1D705}}{1-\unicode[STIX]{x1D702}}r^{1-\unicode[STIX]{x1D702}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D705}$ specifies the intensity of the force between the particles. Based on this assumption, the viscosity coefficient of the gas in equilibrium is proportional to a certain power of the temperature of the gas, which is usually written as 
$\unicode[STIX]{x1D707}_{ref}(\unicode[STIX]{x1D703}/\unicode[STIX]{x1D703}_{ref})^{\unicode[STIX]{x1D714}}$, where 
$\unicode[STIX]{x1D703}$ is the temperature represented in units of specific energy
$$\begin{eqnarray}\unicode[STIX]{x1D703}=\frac{k_{B}}{\mathfrak{m}}T,\end{eqnarray}$$ with 
$T$ being the temperature in Kelvin, 
$k_{B}$ being the Boltzmann constant and 
$\unicode[STIX]{x1D707}_{ref}$ being the reference viscosity coefficient at temperature 
$\unicode[STIX]{x1D703}_{ref}$. The notation 
$\mathfrak{m}$ is the mass of a single molecule, and the viscosity index 
$\unicode[STIX]{x1D714}$ is related to 
$\unicode[STIX]{x1D702}$ by 
$\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D702}+3)/(2\unicode[STIX]{x1D702}-2)$. When 
$\unicode[STIX]{x1D702}=5$, the model is Maxwell molecules, whose viscosity index is 
$1$; when 
$\unicode[STIX]{x1D702}\rightarrow \infty$, the IPL model reduces to the hard-sphere model, whose viscosity index is 
$1/2$. A detailed introduction to the IPL model based on kinetic models is presented in appendix A.
In the following we use the symbol 
$\unicode[STIX]{x1D707}$ to denote a more familiar ‘first approximation’ of the viscosity coefficient, which is obtained by the truncated series expansion using Sonine polynomials (Vincenti, Kruger & Teichmann Reference Vincenti, Kruger and Teichmann1966). For IPL models, the value of 
$\unicode[STIX]{x1D707}$ can be obtained by the following formula (Bird Reference Bird1994):
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{5\mathfrak{m}(k_{B}T/(\mathfrak{m}\unicode[STIX]{x03C0}))^{1/2}(2k_{B}T/\unicode[STIX]{x1D705})^{2/(\unicode[STIX]{x1D702}-1)}}{8A_{2}(\unicode[STIX]{x1D702})\unicode[STIX]{x1D6E4}[4-2/(\unicode[STIX]{x1D702}-1)]},\end{eqnarray}$$ with 
$A_{2}(\unicode[STIX]{x1D702})$ being a constant depending only on 
$\unicode[STIX]{x1D702}$. Some of the values of this constant are given in table 1.
Table 1. Coefficients 
$A_{2}(\unicode[STIX]{x1D702})$ for different 
$\unicode[STIX]{x1D702}$.
2.2 R13-moment equations
As the main result of this paper, the R13-moment equations for general IPL models will be presented in this section. For convenience, the equations are to be written down using ‘primitive variables’, which are density 
$\unicode[STIX]{x1D70C}$, velocity 
$v_{i}$, temperature 
$\unicode[STIX]{x1D703}$, trace-free stress tensor 
$\unicode[STIX]{x1D70E}_{ij}$ and heat flux 
$q_{i}$. All the indices run from 
$1$ to 
$3$. Owing to the constraint 
$\unicode[STIX]{x1D70E}_{11}+\unicode[STIX]{x1D70E}_{22}+\unicode[STIX]{x1D70E}_{33}=0$, these variables amount to 13 quantities, as in Grad (Reference Grad1949). In the following we use the Einstein summation convection without using superscripts. For instance, the above constraint will be written as 
$\unicode[STIX]{x1D70E}_{ii}=0$.
With these 13 variables, the equations for 
$\unicode[STIX]{x1D70C}$, 
$v_{i}$ and 
$\unicode[STIX]{x1D703}$ can be written as
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\displaystyle \frac{\text{d}\unicode[STIX]{x1D70C}}{\text{d}t}}+\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}=0,\\ \displaystyle \unicode[STIX]{x1D70C}{\displaystyle \frac{\text{d}v_{i}}{\text{d}t}}+\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}}+\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}}+{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ik}}{\unicode[STIX]{x2202}x_{k}}}=0,\\ \displaystyle \frac{3}{2}\unicode[STIX]{x1D70C}{\displaystyle \frac{\text{d}\unicode[STIX]{x1D703}}{\text{d}t}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}+{\displaystyle \frac{\unicode[STIX]{x2202}q_{k}}{\unicode[STIX]{x2202}x_{k}}}+\unicode[STIX]{x1D70E}_{kl}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{l}}}=0,\end{array}\right\}\end{eqnarray}$$ which are, in fact, the conservation laws of mass, momentum and total energy represented by primitive variables. More precisely, from the above equations, one can derive the equations for the momentum density 
$\unicode[STIX]{x1D70C}v_{i}$ with momentum flux 
$\unicode[STIX]{x1D70C}(v_{i}v_{k}+\unicode[STIX]{x1D703}\unicode[STIX]{x1D6FF}_{ik})+\unicode[STIX]{x1D70E}_{ik}$, and for the energy density 
$\frac{1}{2}\unicode[STIX]{x1D70C}(v_{i}v_{i}+3\unicode[STIX]{x1D703})$ with energy flux 
$\frac{1}{2}\unicode[STIX]{x1D70C}v_{k}(v_{i}v_{i}+5\unicode[STIX]{x1D703})+\unicode[STIX]{x1D70E}_{ik}v_{i}+q_{k}$. To close the above system, the evolution of the stress tensor 
$\unicode[STIX]{x1D70E}_{kl}$ and the heat flux 
$q_{k}$ needs to be specified. The system (2.2) turns out to be Euler equations if 
$\unicode[STIX]{x1D70E}_{ik}$ and 
$q_{k}$ are set to zero. Finer models based on Chapman–Enskog expansion, such as the Navier–Stokes–Fourier equations, the Burnett equations and the super-Burnett equations, represent 
$\unicode[STIX]{x1D70E}_{ik}$ and 
$q_{k}$ using derivatives of 
$\unicode[STIX]{x1D70C}$, 
$v_{i}$ and 
$\unicode[STIX]{x1D703}$. Following Grad (Reference Grad1949), the 13-moment equations describe the evolution of 
$\unicode[STIX]{x1D70E}_{ik}$ and 
$q_{k}$ by supplementing (2.2) with additional equations. Here, we adopt the form used in Struchtrup (Reference Struchtrup2005b, equations (6.5) and (6.6)) and write down these equations as
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\displaystyle \frac{\text{d}\unicode[STIX]{x1D70E}_{ij}}{\text{d}t}}+\unicode[STIX]{x1D70E}_{ij}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}+\frac{4}{5}{\displaystyle \frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}+2\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}+2\unicode[STIX]{x1D70E}_{k\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}}+{\displaystyle \frac{\unicode[STIX]{x2202}m_{ijk}^{(\unicode[STIX]{x1D702})}}{\unicode[STIX]{x2202}x_{k}}}=\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},1)}+\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},2)},\\ \displaystyle \begin{array}{@{}rcl@{}}\; & \; & \displaystyle {\displaystyle \frac{\text{d}q_{i}}{\text{d}t}}+\frac{5}{2}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}}+\frac{5}{2}\unicode[STIX]{x1D70E}_{ik}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}}+\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ik}}{\unicode[STIX]{x2202}x_{k}}}-\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ik}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k}}}+\frac{7}{5}q_{k}{\displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{k}}}\\ \; & \; & \displaystyle \quad +\,\frac{2}{5}q_{k}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{i}}}+\frac{7}{5}q_{i}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}-\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x2202}x_{k}}}\\ \; & \; & \displaystyle \quad +\,C^{(\unicode[STIX]{x1D702})}\left(\unicode[STIX]{x1D70E}_{ik}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}}+\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ik}}{\unicode[STIX]{x2202}x_{k}}}\right)+\frac{1}{2}{\displaystyle \frac{\unicode[STIX]{x2202}R_{ik}^{(\unicode[STIX]{x1D702})}}{\unicode[STIX]{x2202}x_{k}}}+\frac{1}{6}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}}{\unicode[STIX]{x2202}x_{i}}}+m_{ijk}^{(\unicode[STIX]{x1D702})}{\displaystyle \frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}}=Q_{i}^{(\unicode[STIX]{x1D702},1)}+Q_{i}^{(\unicode[STIX]{x1D702},2)},\end{array}\end{array}\right\}\end{eqnarray}$$ where the angular brackets represent the symmetric and trace-free part of a tensor (
$\unicode[STIX]{x1D61B}_{\langle ij\rangle }=\frac{1}{2}(\unicode[STIX]{x1D61B}_{ij}+\unicode[STIX]{x1D61B}_{ji})-\frac{1}{3}\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x1D61B}_{kk}$ for any two-tensor 
$\unicode[STIX]{x1D61B}$), and some values of the constant 
$C^{(\unicode[STIX]{x1D702})}$ are listed in table 2. Note that Struchtrup (Reference Struchtrup2005b, equations (6.5) and (6.6)) uses the notation 
$u_{ijk}^{0}$, 
$u_{ik}^{1}$ and 
$w^{2}$, while we use the notation 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$, 
$R_{ij}^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$. The relations are
$$\begin{eqnarray}m_{ijk}^{(\unicode[STIX]{x1D702})}=u_{ijk}^{0},\quad R_{ij}^{(\unicode[STIX]{x1D702})}=u_{ij}^{1}-(7+2C^{(\unicode[STIX]{x1D702})})\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij},\quad \unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}=w^{2}.\end{eqnarray}$$ In (2.3), the newly introduced variables 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$, 
$R_{ik}^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ are the moments contributing to second- and higher-order terms in the Chapman–Enskog expansion, and the right-hand sides 
$\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},1)},\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},2)}$ and 
$Q_{i}^{(\unicode[STIX]{x1D702},1)},Q_{i}^{(\unicode[STIX]{x1D702},2)}$ come from the collisions between gas molecules. (This holds for any molecular potential. We refer the reader to Struchtrup (Reference Struchtrup2005b, equation (8,14)), where our 
$R_{ij}^{(\unicode[STIX]{x1D702})}$ is denoted by 
$w_{ij}^{1}$.) Here we assume that the collision is linearized about the local Maxwellian. Up to now, the system is exact but still not closed. The main contribution of this work is to close the system by providing expressions for 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$, 
$R_{ik}^{(\unicode[STIX]{x1D702})}$, 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ and the right-hand sides using the 13 moments. Note that the moments 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$, 
$R_{ik}^{(\unicode[STIX]{x1D702})}$, 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ are quantities appearing in Grad’s 26-moment theory. More specifically, 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$ is the three-tensor formed by all trace-free third-order moments, and 
$R_{ik}^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ are fourth-order moments. In the following we first provide the expressions for 
$\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},1)}$ and 
$Q_{i}^{(\unicode[STIX]{x1D702},1)}$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \begin{array}{@{}rcl@{}}\displaystyle \unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},1)}\; & =\; & \displaystyle D_{0}^{(\unicode[STIX]{x1D702})}\frac{\unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D70E}_{ij}+D_{1}^{(\unicode[STIX]{x1D702})}\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{k}}}\unicode[STIX]{x1D70E}_{j\!\rangle k}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{\langle \!i}}}\unicode[STIX]{x1D70E}_{j\!\rangle k}\right)+D_{2}^{(\unicode[STIX]{x1D702})}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}\unicode[STIX]{x1D70E}_{ij}+D_{3}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}},\\ \; & \; & \displaystyle +\,D_{4}^{(\unicode[STIX]{x1D702})}q_{\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}+D_{5}^{(\unicode[STIX]{x1D702})}q_{\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}+D_{6}^{(\unicode[STIX]{x1D702})}\frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }},\end{array}\\ \displaystyle \begin{array}{@{}rcl@{}}\displaystyle Q_{i}^{(\unicode[STIX]{x1D702},1)}\; & =\; & \displaystyle E_{0}^{(\unicode[STIX]{x1D702})}\frac{\unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}q_{i}+E_{1}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D70E}_{ik}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}}+E_{2}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ik}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k}}}+E_{3}^{(\unicode[STIX]{x1D702})}q_{k}\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{i}}}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{k}}}\right)\\ \; & \; & \displaystyle +\,E_{4}^{(\unicode[STIX]{x1D702})}q_{i}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}+E_{5}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ki}}{\unicode[STIX]{x2202}x_{k}}}+E_{6}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}},\end{array}\end{array}\right\}\end{eqnarray}$$ where the coefficients 
$D_{i}^{(\unicode[STIX]{x1D702})}$ and 
$E_{i}^{(\unicode[STIX]{x1D702})}$ are partially tabulated in table 3. Note that we have intentionally split the right-hand sides in (2.3) into two parts, so that the ‘generalized Grad 13-moment (GG13) equations’ can be extracted from (2.3) by setting
$$\begin{eqnarray}m_{ijk}^{(\unicode[STIX]{x1D702})}=R_{ik}^{(\unicode[STIX]{x1D702})}=\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}=\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},2)}=Q_{i}^{(\unicode[STIX]{x1D702},2)}=0.\end{eqnarray}$$The GG13 equations are introduced (Struchtrup Reference Struchtrup2005a) by the order of magnitude method, and its fully linearized version for hard spheres has been derived in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013).
Table 2. Coefficients 
$C^{(\unicode[STIX]{x1D702})}$ for different 
$\unicode[STIX]{x1D702}$.
Table 3. Coefficient of 
$A_{i}^{(\unicode[STIX]{x1D702})}$, 
$B_{i}^{(\unicode[STIX]{x1D702})}$, 
$C_{i}^{(\unicode[STIX]{x1D702})}$, 
$D_{i}^{(\unicode[STIX]{x1D702})}$ and 
$E_{i}^{(\unicode[STIX]{x1D702})}$ for different 
$\unicode[STIX]{x1D702}$.
To give the R13 equations, we need to close (2.3) by specifying 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$, 
$R_{ik}^{(\unicode[STIX]{x1D702})}$, 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$, 
$\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},2)}$ and 
$Q_{i}^{(\unicode[STIX]{x1D702},2)}$. The closure depends on the specific form of the collision model. Here, we again assume that the collision is linearized about the local Maxwellian. Thus, the R13 theory gives the following closure:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle m_{ijk}^{(\unicode[STIX]{x1D702})}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}}\left(A_{1}^{(\unicode[STIX]{x1D702})}q_{\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k\!\rangle }}}+A_{2}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{\langle \!ij}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k\!\rangle }}}+A_{3}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D70E}_{\langle \!ij}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k\!\rangle }}}+A_{4}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\langle \!ij}}{\unicode[STIX]{x2202}x_{k}\!\rangle }}\right),\\ \displaystyle \unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}}\left(B_{1}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}q_{k}}{\unicode[STIX]{x2202}x_{k}}}+B_{2}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}{\displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}}+B_{3}^{(\unicode[STIX]{x1D702})}q_{k}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}}+B_{4}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}q_{k}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k}}}\right),\\ \displaystyle \begin{array}{@{}rcl@{}}\displaystyle R_{ij}^{(\unicode[STIX]{x1D702})}\; & =\; & \displaystyle C_{0}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}+\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D703}\unicode[STIX]{x1D70C}}\left(C_{1}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}+C_{2}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{k\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}}+C_{3}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}^{2}\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}\right.\\ \; & \; & \displaystyle +\left.C_{4}^{(\unicode[STIX]{x1D702})}q_{\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}+C_{5}^{(\unicode[STIX]{x1D702})}\unicode[STIX]{x1D703}q_{\langle \!i}{\displaystyle \frac{\unicode[STIX]{x2202}\ln \unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j\!\rangle }}}\right).\end{array}\end{array}\right\}\end{eqnarray}$$ Some values of the coefficients 
$A_{i}^{(\unicode[STIX]{x1D702})}$, 
$B_{i}^{(\unicode[STIX]{x1D702})}$, and 
$C_{i}^{(\unicode[STIX]{x1D702})}$ are given in table 3. The full expressions of 
$\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},2)}$ and 
$Q_{i}^{(\unicode[STIX]{x1D702},2)}$ are quite lengthy and we provide them in appendix C. A simple case is 
$\unicode[STIX]{x1D702}=5$, for which we have
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{ij}^{(5,2)}=0,\quad Q_{i}^{(5,2)}=0,\end{eqnarray}$$ and the corresponding model matches that derived for Maxwell molecules by Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) (with terms nonlinear in 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$ removed since we use the linearized collision model).
2.3 Discussion on the order of accuracy
One possible way to describe the accuracy of the moment models in the near-continuum regime is to use the notion of ‘order of accuracy’ (Struchtrup Reference Struchtrup2005a,Reference Struchtrupb). In such a regime, the Knudsen number 
$Kn$, i.e. the ratio of the mean free path to the characteristic length of the problem, is regarded as a small number. Thus, Chapman–Enskog expansion can be applied, and all non-equilibrium moments are expanded into power series of 
$Kn$, e.g.
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70E}_{ij}=Kn\unicode[STIX]{x1D70E}_{ij}^{(1)}+Kn^{2}\unicode[STIX]{x1D70E}_{ij}^{(2)}+Kn^{3}\unicode[STIX]{x1D70E}_{ij}^{(3)}+\cdots \,,\\ q_{i}=Knq_{i}^{(1)}+Kn^{2}q_{i}^{(2)}+Kn^{3}q_{i}^{(3)}+\cdots \,.\end{array}\right\}\end{eqnarray}$$ By asymptotic analysis, all these terms can be represented by the conservative variables and their derivatives. Truncating the above series up to the term 
$Kn^{k}$ and inserting the result into (2.2), one obtains moment equations with 
$k$th order of accuracy. By this approach, the models derived from Chapman–Enskog expansion from zeroth to third order are, respectively, the Euler equations, the Navier–Stokes–Fourier equations, the Burnett equations and the super-Burnett equations. These equations contain only equilibrium variables: density, velocity and temperature.
In the 13-moment model, one can also assume 
$Kn$ is small and apply the expansion (2.8) to obtain models including only equilibrium variables. Suppose the second-order Chapman–Enskog expansion of a moment model agrees with the Burnett equations, while its third-order Chapman–Enskog expansion differs from super-Burnett equations, then we say that the moment model has the second-order accuracy (or Burnett order). For example, Grad’s 13-moment equations have the first-order accuracy for general IPL potentials, but have second-order accuracy for Maxwell molecules; GG13 equations are extensions to Grad’s 13-moment theory to achieve second-order accuracy for all molecule potentials. In general, the expansion (2.8) is usually obtained by multiplying the equations of 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$ in the moment system by 
$Kn$. Therefore, in most cases, a 13-moment model has 
$k$th-order accuracy if the equations for 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$ (2.3) are accurate up to the 
$(k-1)$th order. R13 equations have third-order accuracy, as the closure (2.6) provides the equations (2.3) exact second-order contributions.
In the literature, there exist similar 13-moment models obtained by other approaches to including second-order derivatives in the equations for 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$. For instance, the relaxed Burnett equations (Jin & Slemrod Reference Jin and Slemrod2001) are also derived for arbitrary interaction potentials, and as mentioned in Jin & Slemrod (Reference Jin and Slemrod2001) and Struchtrup (Reference Struchtrup2005b), these equations have second-order accuracy. Another similar model is the nonlinear coupled constitutive relation equation model (Myong Reference Myong1999). These equations do not include information from the Burnett order, and therefore they have first-order accuracy and distinguish different interaction models by viscosity and heat conductivity coefficients. The full R13 models for Maxwell molecules and the Bhatnagar–Gross–Krook (BGK) model, which have third-order accuracy, have been derived in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) and Struchtrup (Reference Struchtrup2005b), and the linear R13 equations for the hard-sphere model have been derived in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013).
2.4 Linear stability and dispersion
For the newly proposed R13 equations for IPL models, we are going to check some of its basic properties in this work. In this section, we focus on the linear properties, including stability in time and space, and the dispersion and damping of sound waves.
Following Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003, Reference Struchtrup and Torrilhon2013) and Struchtrup (Reference Struchtrup2005b), we apply the analysis to one-dimensional linear dimensionless equations. The linearization is performed about a global equilibrium state with density 
$\unicode[STIX]{x1D70C}_{0}$, zero velocity and temperature 
$\unicode[STIX]{x1D703}_{0}$. The derivation of the one-dimensional linear dimensionless equations consists of the following steps.
(i) Introduce the small dimensionless variables
$\hat{\unicode[STIX]{x1D70C}}$, $\hat{\unicode[STIX]{x1D703}}$, $\hat{v}_{i}$, $\hat{\unicode[STIX]{x1D70E}}_{ij}$ and $\hat{q}_{i}$ by (2.9a-e)$$\begin{eqnarray}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(1+\hat{\unicode[STIX]{x1D70C}}),\quad \unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{0}(1+\hat{\unicode[STIX]{x1D703}}),\quad v_{i}=\sqrt{\unicode[STIX]{x1D703}_{0}}\hat{v}_{i},\quad \unicode[STIX]{x1D70E}_{ij}=\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D703}_{0}\hat{\unicode[STIX]{x1D70E}}_{ij},\quad q_{i}=\unicode[STIX]{x1D70C}_{0}\sqrt{\unicode[STIX]{x1D703}_{0}}^{3}\hat{q}_{i}.\end{eqnarray}$$(ii) Let
$L$ be the characteristic length, and define the dimensionless space and time variables by (2.10a,b)$$\begin{eqnarray}x_{i}=L\hat{x}_{i},\quad t=\frac{L}{\sqrt{\unicode[STIX]{x1D703}_{0}}}\hat{t}.\end{eqnarray}$$(iii) Substitute (2.9) and (2.10) into the R13 equations (2.2), (2.3) and (2.6), and drop all the terms nonlinear in the variables with hats introduced in (2.9).
(iv) Reduce the resulting equations to a one-dimensional system by dropping all the terms with derivatives with respect to
$x_{2}$ and $x_{3}$, and setting (2.11)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\hat{v}_{1}=\hat{v},\quad \hat{q}_{1}=\hat{q},\quad \hat{\unicode[STIX]{x1D70E}}_{11}=\hat{\unicode[STIX]{x1D70E}},\quad \hat{\unicode[STIX]{x1D70E}}_{22}=\hat{\unicode[STIX]{x1D70E}}_{33}=-{\textstyle \frac{1}{2}}\hat{\unicode[STIX]{x1D70E}},\\ \hat{v}_{2}=\hat{v}_{3}=\hat{q}_{2}=\hat{q}_{3}=\hat{\unicode[STIX]{x1D70E}}_{12}=\hat{\unicode[STIX]{x1D70E}}_{23}=\hat{\unicode[STIX]{x1D70E}}_{13}=0.\end{array}\right\}\end{eqnarray}$$
The resulting equations can be written down more neatly if we introduce the Knudsen number
$$\begin{eqnarray}Kn=\frac{\unicode[STIX]{x1D707}_{0}\sqrt{\unicode[STIX]{x1D703}_{0}}}{\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D703}_{0}L},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D707}_{0}$ is the viscosity coefficient at temperature 
$\unicode[STIX]{x1D703}_{0}$. For all IPL models, the one-dimensional linear dimensionless equations have the form
where 
$\unicode[STIX]{x1D6FC}_{i}^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6FD}_{i}^{(\unicode[STIX]{x1D702})}$ depend only on 
$\unicode[STIX]{x1D702}$, and their values for some choices of 
$\unicode[STIX]{x1D702}$ are listed in table 4. In (
), if we replace all the terms with double underlines by zero, we obtain the linearized GG13 equations. Furthermore, if we set all the terms with both single and double underlines to zero, then the result is the linearized Navier–Stokes–Fourier equations. The left-hand side of (
2.13) comes from the advection, and the right-hand side comes from the collision. By table 4, it can be clearly seen that when 
$\unicode[STIX]{x1D702}=5$ (Maxwell molecules), due to the simplicity of the collision operator, all the underlined terms on the right-hand side disappear.
Table 4. Coefficients of the linearized R13 system for different 
$\unicode[STIX]{x1D702}$.
For simplicity, we omit the hats on the variables hereafter. In general, the linear GG13 or R13 system has the form
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}u_{A}}{\unicode[STIX]{x2202}t}}+{\mathcal{A}}_{1}^{(\unicode[STIX]{x1D702})}{\displaystyle \frac{\unicode[STIX]{x2202}u_{A}}{\unicode[STIX]{x2202}x}}+{\mathcal{A}}_{2}^{(\unicode[STIX]{x1D702})}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}u_{A}}{\unicode[STIX]{x2202}x^{2}}}+{\mathcal{A}}_{3}^{(\unicode[STIX]{x1D702})}u_{A}=0,\end{eqnarray}$$ where 
$u_{A}=(\unicode[STIX]{x1D70C},v,\unicode[STIX]{x1D703},\unicode[STIX]{x1D70E},q)^{\text{T}}$ and the matrices 
${\mathcal{A}}_{i}^{(\unicode[STIX]{x1D702})}$ are constant matrices which can be observed from (2.13). To study the linear waves, we consider the plane wave solution
$$\begin{eqnarray}u_{A}(x,t)=\tilde{u} _{A}\exp [\text{i}(\unicode[STIX]{x1D6FA}t-kx)],\end{eqnarray}$$ where 
$\tilde{u} _{A}$ is the initial amplitude of the wave, 
$\unicode[STIX]{x1D6FA}$ is frequency and 
$k$ is the wavenumber. Inserting the above solution into (2.14) yields
$$\begin{eqnarray}{\mathcal{G}}^{(\unicode[STIX]{x1D702})}\tilde{u} _{A}=0,\quad \text{where }{\mathcal{G}}^{(\unicode[STIX]{x1D702})}=(\text{i}\unicode[STIX]{x1D6FA}-\text{i}k{\mathcal{A}}_{1}^{(\unicode[STIX]{x1D702})}-k^{2}{\mathcal{A}}_{2}^{(\unicode[STIX]{x1D702})}+{\mathcal{A}}_{3}^{(\unicode[STIX]{x1D702})}),\end{eqnarray}$$and the existence of non-trivial solutions requires
$$\begin{eqnarray}\det [{\mathcal{G}}^{(\unicode[STIX]{x1D702})}]=0.\end{eqnarray}$$ From (2.17), we can obtain the relation between 
$\unicode[STIX]{x1D6FA}$ and 
$k$, and thus all the desired properties such as the amplification and dispersion of the linear waves can naturally be obtained. In the following analysis, we choose 
$Kn=1$ to obtain quantitative results.
Remark 1. Equations (2.13) in the case 
$\unicode[STIX]{x1D702}=\infty$ can be used to compare with the results in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013) for cross-checking. Small deviation between our coefficients and the coefficients in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013) can be observed. For example, in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013), the value of 
$\unicode[STIX]{x1D6FC}_{2}^{(\infty )}$ is 
$-0.98632$, while our analysis gives 
$\unicode[STIX]{x1D6FC}_{2}^{(\infty )}=-0.9842$. We believe that such discrepancies are due to a different truncation when inverting the collision operator during the derivation. According to the method reported in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013, equation (23)), our result is probably more accurate since we preserve more terms in the truncation. Details are given in § 3.
2.4.1 Linear stability in time and space
We first discuss the linear stability of the GG13 and R13 systems in time and space. For the time stability, we require that the norm of the amplitude decreases with time for any given wavenumber 
$k\in \mathbb{R}$. Precisely, if we assume 
$\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6FA}_{r}(k)+\text{i}\unicode[STIX]{x1D6FA}_{i}(k)$, the time stability requires 
$\unicode[STIX]{x1D6FA}_{i}(k)\geqslant 0$. Figure 1 shows possible values of 
$\unicode[STIX]{x1D6FA}_{i}(k)$ on the complex plane. Note that 
$\unicode[STIX]{x1D6FA}_{i}(k)$ is a multi-valued function since (2.17) may have multiple solutions for a given 
$k$. It is observed that, for all choices 
$\unicode[STIX]{x1D702}$, the values of 
$\unicode[STIX]{x1D6FA}(k)$ are always located on the upper half of the complex plane, indicating linear stability for both R13 and GG13 equations.
Figure 1. Damping coefficients 
$\unicode[STIX]{x1D6FA}_{i}(k)$ of the GG13 and R13 systems for different 
$\unicode[STIX]{x1D702}$.
For the stability in space, we require that, for a given wave frequency, the amplitude should not increase along the direction of wave propagation. Now, we assume that 
$\unicode[STIX]{x1D6FA}\in \mathbb{R}$ is given, and let 
$k=k_{r}(\unicode[STIX]{x1D6FA})+\text{i}k_{i}(\unicode[STIX]{x1D6FA})$. Then the wave is stable in space if 
$k_{r}(\unicode[STIX]{x1D6FA})k_{i}(\unicode[STIX]{x1D6FA})\leqslant 0$. Figure 2 shows the values of 
$k$ in the complex plane with 
$\unicode[STIX]{x1D6FA}$ as the parameter. The results show that all the curves do not enter the upper right or the lower left quadrant for both the R13 and GG13 equations, showing the spatial stability for both models. Again, such a stability result holds for all 
$\unicode[STIX]{x1D702}$ considered in our experiments.
Figure 2. The solutions 
$k(\unicode[STIX]{x1D6FA})$ of the dispersion relation in the complex plane with 
$\unicode[STIX]{x1D6FA}$ as the parameter of the GG13 and R13 systems for different 
$\unicode[STIX]{x1D702}$.
2.4.2 Dispersion and damping
We proceed by discussing the phase speeds as functions of frequency for the R13 and GG13 systems. For a given wave frequency 
$\unicode[STIX]{x1D6FA}$, we define the damping rate 
$\unicode[STIX]{x1D6FC}$ and the wave speed 
$v_{ph}$ by
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}=-k_{i}(\unicode[STIX]{x1D6FA}),\quad v_{ph}=\frac{\unicode[STIX]{x1D6FA}}{k_{r}(\unicode[STIX]{x1D6FA})}.\end{eqnarray}$$ For the Euler equations, where 
$\unicode[STIX]{x1D70E}=q=0$ in (2.13), the absolute phase velocity is 
$|v_{ph}|=\sqrt{5/3}$ for all wave frequencies 
$\unicode[STIX]{x1D6FA}$. For R13 and GG13 equations, 
$v_{ph}$ depends on 
$\unicode[STIX]{x1D6FA}$, causing the dispersion of sound waves. Here, we define the dimensionless phase speed 
$c_{ph}=v_{ph}/\sqrt{5/3}$ and plot 
$c_{ph}$ as a function of the frequency 
$\unicode[STIX]{x1D6FA}$ in figure 3 for the R13 and GG13 equations. Note that 
$c_{ph}(\unicode[STIX]{x1D6FA})$ is also a multi-valued function, and in figure 3, we only plot the positive phase velocities. For the GG13 equations, the phase velocity has an upper limit, indicating the hyperbolic nature of the system, while the R13 equations can achieve infinitely large phase velocities. In general, the phase speeds do not change much as 
$\unicode[STIX]{x1D702}$ varies, which predicts similar behaviour of sound waves in different monatomic gases.
Figure 3. Phase speed 
$c_{ph}$ over frequency 
$\unicode[STIX]{x1D6FA}$ of the GG13 and R13 systems for different 
$\unicode[STIX]{x1D702}$.
To study the phase speeds for large frequency waves, we plot the inverse wave speed 
$1/c_{ph}$ as a function of the inverse frequency 
$1/\unicode[STIX]{x1D6FA}$ in figure 4(a), and the reduced damping rate 
$\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6FA}$ is plotted in figure 4(b) also as a function of 
$1/\unicode[STIX]{x1D6FA}$. In these figures, only the mode with the weakest damping is given. Figure 4(a) shows that, for the GG13 equations, the phase velocity increases monotonically as the frequency increases, while for the R13 equations, the wave slows down as 
$\unicode[STIX]{x1D6FA}$ reaches a value close to 
$1$. Such an observation agrees with the results in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003, Reference Struchtrup and Torrilhon2013), while the experimental results for argon (Meyer & Sessler Reference Meyer and Sessler1957) suggest the monotonicity of the phase velocity, which is closer to the prediction of the GG13 equations. Note that, here, we set the Knudsen number 
$Kn$ to be 
$1$. For another Knudsen number, the actual frequency of the wave should be 
$\unicode[STIX]{x1D6FA}/Kn$. This means that, if we consider a wave with a fixed actual frequency travelling in the gas with a low Knudsen number, we need to focus on large values of 
$\unicode[STIX]{x1D6FA}^{-1}$. Indeed, when 
$\unicode[STIX]{x1D6FA}^{-1}>1$, R13 models give a better approximation of the phase velocity, while for small 
$\unicode[STIX]{x1D6FA}^{-1}$, which corresponds to large Knudsen numbers, there is no guarantee as to whether the R13 or GG13 is superior, and the reason why GG13 provides better prediction requires further investigation. Nevertheless, for the damping rate shown in figure 4(b), R13 equations give significantly better agreement with the experimental data for the whole range of frequencies.
Figure 4. Inverse phase speed and damping over frequency 
$\unicode[STIX]{x1D6FA}$ of the GG13 and R13 systems for different 
$\unicode[STIX]{x1D702}$. The bullets are the experimental results for argon Meyer & Sessler (Reference Meyer and Sessler1957).
3 Derivation of moment systems
In this section, we provide the detailed procedure to derive the GG13 and R13 equations. In general, both models can be derived from infinite moment equations by the method of order of magnitude. Following Grad (Reference Grad1949), the moments can be considered as the coefficients in the series expansion of the distribution function in the gas kinetic theory. The distribution function is a function of position 
$\boldsymbol{x}$, particle velocity 
$\unicode[STIX]{x1D743}$ and time 
$t$, and is a mesoscopic description of fluid states in statistical physics. The moment method proposed by Grad (Reference Grad1949) is one of the methods to derive macroscopic models from the kinetic theory. Our starting point is the same as Grad (Reference Grad1949), but we adopt the form used in Kumar (Reference Kumar1966a), which expands the distribution function 
$f(\boldsymbol{x},\unicode[STIX]{x1D743},t)$ as
$$\begin{eqnarray}f(\boldsymbol{x},\unicode[STIX]{x1D743},t)=\mathop{\sum }_{l=0}^{+\infty }\mathop{\sum }_{m=-l}^{l}\mathop{\sum }_{n=0}^{+\infty }f_{lmn}(\boldsymbol{x},t)\unicode[STIX]{x1D713}_{lmn}(\boldsymbol{x},\unicode[STIX]{x1D743},t),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D713}_{lmn}(\cdot )$ is the basis function based on Sonine polynomials and spherical harmonics, the detailed form of which is listed in appendix B. Here, 
$\boldsymbol{v}=(v_{1},v_{2},v_{3})^{\text{T}}$ is the velocity vector. The coefficients 
$f_{lmn}$ satisfy 
$\overline{f_{lmn}}=(-1)^{m}f_{l,-m,n}$, and they are related to Grad’s 13 moments by
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}f_{000}=\unicode[STIX]{x1D70C},\quad f_{1m0}=0,\quad m=-1,0,1,\quad f_{001}=0,\\ \unicode[STIX]{x1D70E}_{11}=\sqrt{2}\text{Re}(f_{220})-f_{200}/\sqrt{3},\quad \unicode[STIX]{x1D70E}_{12}=-\sqrt{2}\text{Im}(f_{220}),\quad \unicode[STIX]{x1D70E}_{13}=-\sqrt{2}\text{Re}(f_{210}),\\ \unicode[STIX]{x1D70E}_{22}=-\sqrt{2}\text{Re}(f_{220})-f_{200}/\sqrt{3},\quad \unicode[STIX]{x1D70E}_{23}=\sqrt{2}\text{Im}(f_{210}),\quad \unicode[STIX]{x1D70E}_{33}=2f_{200}/\sqrt{3},\\ q_{1}=\sqrt{5}\text{Re}(f_{111}),\quad q_{2}=-\sqrt{5}\text{Im}(f_{111}),\quad q_{3}=-\sqrt{5/2}f_{101}.\end{array}\right\}\end{eqnarray}$$These relations indicate the equivalence between Grad’s 13 moments and the following 13 variables:
$$\begin{eqnarray}f_{000},v_{1},v_{2},v_{3},\unicode[STIX]{x1D703},f_{220},f_{210},f_{200},f_{2,-1,0},f_{2,-2,0},f_{111},f_{101},f_{1,-1,1}.\end{eqnarray}$$In the following, we focus only on the derivation of equations for these quantities.
The exact evolution equations for 
$f_{lmn}$ have been derived from the Boltzmann equation with a linearized collision operator in Cai & Torrilhon (Reference Cai and Torrilhon2018). In general, the equations for other 
$f_{lmn}$ have the form
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}f_{lmn}}{\unicode[STIX]{x2202}t}}+S_{lmn}+T_{lmn}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}f_{lmn^{\prime }},\end{eqnarray}$$ where 
$S_{lmn}$ contains time derivatives and 
$T_{lmn}$ contains spatial derivatives. More precisely, 
$S_{lmn}$ is the linear combination of the terms
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}f_{l,m^{\prime },n-1},\quad {\displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}t}}f_{l-1,m^{\prime },n},\quad \text{and}\quad {\displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}t}}f_{l+1,m^{\prime },n-1},\end{eqnarray}$$ with 
$i=1,2,3$ and 
$m^{\prime }=m-1,m,m+1$, and 
$T_{lmn}$ has the form
$$\begin{eqnarray}T_{lmn}=\mathop{\sum }_{l^{\prime },m^{\prime },n^{\prime }}\left(\unicode[STIX]{x1D6FC}_{lmn}^{l^{\prime }m^{\prime }n^{\prime }}(\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{v},\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D703},\boldsymbol{v},\unicode[STIX]{x1D703})f_{l^{\prime }m^{\prime }n^{\prime }}+\mathop{\sum }_{i=1}^{3}\unicode[STIX]{x1D6FD}_{lmn,i}^{l^{\prime }m^{\prime }n^{\prime }}(\boldsymbol{v},\unicode[STIX]{x1D703}){\displaystyle \frac{\unicode[STIX]{x2202}f_{l^{\prime }m^{\prime }n^{\prime }}}{\unicode[STIX]{x2202}x_{i}}}\right),\end{eqnarray}$$ which shows that 
$T_{lmn}$ is linear in all the coefficients 
$f_{l^{\prime }m^{\prime }n^{\prime }}$, while the linear coefficients are nonlinear functions of 
$\boldsymbol{v}$, 
$\unicode[STIX]{x1D703}$ and their spatial derivatives. The convection term 
$T_{lmn}$ also has the following properties.
- (P1)
The differential operator appears only once in each coefficient
$\unicode[STIX]{x1D6FC}_{lmn}^{l^{\prime }m^{\prime }n^{\prime }}(\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{v},\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D703},\boldsymbol{v},\unicode[STIX]{x1D703})$.- (P2)
The coefficients
$\unicode[STIX]{x1D6FC}_{lmn}^{l^{\prime }m^{\prime }n^{\prime }}(\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{v},\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D703},\boldsymbol{v},\unicode[STIX]{x1D703})$ and $\unicode[STIX]{x1D6FD}_{lmn,i}^{l^{\prime }m^{\prime }n^{\prime }}(\boldsymbol{v},\unicode[STIX]{x1D703})$ are non-zero only if $$\begin{eqnarray}l+2n-3\leqslant l^{\prime }+2n^{\prime }\leqslant l+2n+1,\quad l-2\leqslant l^{\prime }\quad \text{and}\quad n-2\leqslant n^{\prime }\leqslant n+1.\end{eqnarray}$$
The second property (P2) will play an important role in the derivation of the moment equations.
The precise expressions of 
$S_{lmn}$ and 
$T_{lmn}$ are given in appendix B, and on the right-hand side of (3.4), 
$a_{lnn^{\prime }}$ are pure numbers for all IPL models. Note that (3.4) has already included the conservation laws (2.2), which can be obtained by setting 
$(l,m,n)$ to be
$$\begin{eqnarray}(0,0,0),\quad (1,-1,0),\quad (1,0,0),\quad (1,1,0),\quad (0,0,1).\end{eqnarray}$$The subsequent derivation may include tedious formulae, and in our implementation, all the calculations are done by the computer algebra system Wolfram Mathematica. In the following, we only describe the algorithm we use in the Mathematica code, and do not write out the lengthy intermediate results in the calculation process.
Remark 2. Note that complex basis functions 
$\unicode[STIX]{x1D713}_{lmn}$ are introduced in the expansion (3.1), resulting in complex coefficients 
$f_{lmn}$, which seems to complicate the derivation. In Grad’s original formulation (Grad Reference Grad1949), basis functions based on Cartesian coordinates are considered, so that all the coefficients are real. However, in Grad’s original expansion using Hermite polynomials, one cannot find 13 coefficients that match exactly all the moments in his 13-moment equations. Our basis functions, which are similar to those in Kumar (Reference Kumar1966b), correspond to the 13 moments very well (see (3.3)). Meanwhile, our basis functions are based on spherical coordinates, so that the rotational invariance of the collision operator can be easily utilized to simplify the calculation. Based on spherical coordinates, complex basis functions can also be avoided by using real spherical harmonics instead of complex spherical harmonics (Bayin Reference Bayin2018, § 1.5.2). However, real spherical harmonics do not have simple recurrence formulae such as Arfken, Weber & Harris (Reference Arfken, Weber and Harris2013, equations (15.150) (15.151)), which will result in a more complicated form of 
$S_{lmn}$ and 
$T_{lmn}$.
3.1 General idea for the moment closure
Before providing the details of the derivation, we would first like to explain the general methodology of the moment closure. To close the advection term, one can observe from (2.3) that we need to provide expressions for the moments 
$m_{ijk}$, 
$R_{ik}$ and 
$\unicode[STIX]{x1D6E5}$, which correspond to the coefficients 
$f_{3m0}$, 
$f_{1m1}$ and 
$f_{002}$, respectively. To close the collision term, it can be seen from (3.4) that we have to provide expressions for infinite terms 
$f_{1mn}$ and 
$f_{2mn}$ for any positive integer 
$n$. In our implementation, this is done by a truncation of the infinite series in (3.4), whose fast convergence has already been demonstrated in Cai & Torrilhon (Reference Cai and Torrilhon2015), so that only a finite number of coefficients need to be considered.
The derivation of the R13 equations is similar to the Chapman–Enskog expansion. However, there are two key differences.
- (K1)
All the moments are to be represented by the 13 moments and their derivatives, while only 5 equilibrium moments are involved in the Chapman–Enskog expansion.
- (K2)
We expect that the equations have the Burnett order and involve only second-order derivatives, while third-order derivatives are involved in Burnett equations.
Due to the first difference, there are in principle infinite versions of the R13 equations in the Burnett order, since the leading-order terms in 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$ can be represented by equilibrium variables
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{ij}\approx -\frac{2\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D6FC}_{2}^{(\unicode[STIX]{x1D702})}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}},\quad q_{i}\approx -\frac{5\unicode[STIX]{x1D707}}{2\unicode[STIX]{x1D6FD}_{5}^{(\unicode[STIX]{x1D702})}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}},\end{eqnarray}$$ as already derived in many different versions of R13 equations for Maxwell molecules (Timokhin et al. Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017). In our derivation, since 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$ are already included in the system, we will avoid using (3.6) to make any replacement except in the derivation of the first-order expressions. More precisely, according to the Chapman–Enskog expansion, the first-order term of the distribution function can be completely represented by the following quantities:
$$\begin{eqnarray}\unicode[STIX]{x1D70C},v_{i},\unicode[STIX]{x1D703},{\displaystyle \frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}},\quad {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}}.\end{eqnarray}$$ At this step, we apply (3.6) to replace the derivatives of 
$v_{i}$ and 
$\unicode[STIX]{x1D703}$ by 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$. By such replacement, one order of derivative can be eliminated, so that (K2) can be automatically achieved.
3.2 Chapman–Enskog expansion of the coefficients
This section is devoted to the details of the asymptotic analysis. As mentioned in the previous section, the idea of a Chapman–Enskog expansion is utilized here to derive models with different orders of accuracy. To begin with, we introduce the scaling 
$t=t^{\prime }/\unicode[STIX]{x1D716}$ and 
$x_{i}=x_{i}^{\prime }/\unicode[STIX]{x1D716}$, and rewrite equations (3.4) with time and spatial variables 
$t^{\prime }$ and 
$x_{i}^{\prime }$. In the resulting equations, a factor 
$\unicode[STIX]{x1D716}^{-1}$ is introduced to the right-hand side of (3.4). In this section, we work on the scaled equations, and the prime symbol on 
$t^{\prime }$ and 
$x_{i}^{\prime }$ is omitted. Based on such a transform, we write down the asymptotic expansion for all the moments as
$$\begin{eqnarray}f_{lmn}=f_{lmn}^{(0)}+\unicode[STIX]{x1D716}f_{lmn}^{(1)}+\unicode[STIX]{x1D716}^{2}f_{lmn}^{(2)}+\unicode[STIX]{x1D716}^{3}f_{lmn}^{(3)}+\cdots \,.\end{eqnarray}$$ In the original Chapman–Enskog expansion of the distribution function 
$f=f^{(0)}+\unicode[STIX]{x1D716}f^{(1)}$, we require that 
$f^{(0)}$ is the local Maxwellian 
${\mathcal{M}}$. The corresponding assumption for the coefficients is
$$\begin{eqnarray}f_{lmn}^{(0)}=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70C}, & \text{if }(l,m,n)=(0,0,0),\\ 0, & \text{otherwise}.\end{array}\right.\end{eqnarray}$$ By (3.7), the terms 
$S_{lmn}$ and 
$T_{lmn}$ can be expanded correspondingly
$$\begin{eqnarray}S_{lmn}=S_{lmn}^{(0)}+\unicode[STIX]{x1D716}S_{lmn}^{(1)}+\unicode[STIX]{x1D716}^{2}S_{lmn}^{(2)}+\unicode[STIX]{x1D716}^{3}S_{lmn}^{(3)}+\cdots \,,\quad T_{lmn}=T_{lmn}^{(0)}+\unicode[STIX]{x1D716}T_{lmn}^{(1)}+\unicode[STIX]{x1D716}^{2}T_{lmn}^{(2)}+\unicode[STIX]{x1D716}^{3}T_{lmn}^{(3)}+\cdots \,.\end{eqnarray}$$ Note that the above expansion is straightforward since both 
$S_{lmn}$ and 
$T_{lmn}$ are linear in all the coefficients 
$f_{lmn}$. Thus, the moment equations turn out to be
$$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}(f_{lmn}^{(0)}+\unicode[STIX]{x1D716}f_{lmn}^{(1)}+\unicode[STIX]{x1D716}^{2}f_{lmn}^{(2)}+\unicode[STIX]{x1D716}^{3}f_{lmn}^{(3)}+\cdots \,)}{\unicode[STIX]{x2202}t}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,(S_{lmn}^{(0)}+\unicode[STIX]{x1D716}S_{lmn}^{(1)}+\unicode[STIX]{x1D716}^{2}S_{lmn}^{(2)}+\unicode[STIX]{x1D716}^{3}S_{lmn}^{(3)}+\cdots \,)+(T_{lmn}^{(0)}+\unicode[STIX]{x1D716}T_{lmn}^{(1)}+\unicode[STIX]{x1D716}^{2}T_{lmn}^{(2)}+\unicode[STIX]{x1D716}^{3}T_{lmn}^{(3)}+\cdots \,)\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{\unicode[STIX]{x1D716}}\left(\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}(f_{lmn^{\prime }}^{(0)}+\unicode[STIX]{x1D716}f_{lmn^{\prime }}^{(1)}+\unicode[STIX]{x1D716}^{2}f_{lmn^{\prime }}^{(2)}+\unicode[STIX]{x1D716}^{3}f_{lmn^{\prime }}^{(3)}+\cdots \,)\right).\end{eqnarray}$$ Matching the terms with the same orders with respect to 
$\unicode[STIX]{x1D716}$, one obtains
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}f_{lmn}^{(k)}}{\unicode[STIX]{x2202}t}}+S_{lmn}^{(k)}+T_{lmn}^{(k)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}f_{lmn^{\prime }}^{(k+1)},\quad k\geqslant 0.\end{eqnarray}$$ In Chapman–Enskog expansion, due to the assumption (3.8), equation (3.11) is only applied to the case 
$(l,n)\neq (0,0),(1,0),(0,1)$, in which the right-hand side of (3.11) is non-zero and (3.11) can provide us expressions for 
$f_{lmn}^{(k+1)}$. For more details about the Chapman–Enskog expansion, we refer the reader to textbooks such as Struchtrup (Reference Struchtrup2005a). In what follows, we introduce a generalized version of the Chapman–Enskog expansion involving 13 moments in the assumption, which is carried out by studying each 
$k$ incrementally. Note that the idea of the following method is also applicable for more general collision models.
3.2.1 First order ($k=0$)
When 
$(l,n)\neq (0,0),(1,0),(0,1)$, using (3.8) and the fact that 
$S_{lmn}$ is a linear combination of (3.5), one can see that 
$f_{lmn}^{(0)}=S_{lmn}^{(0)}=0$. Thus, when 
$k=0$, equation (3.11) becomes
$$\begin{eqnarray}T_{lmn}^{(0)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}f_{lmn^{\prime }}^{(1)}.\end{eqnarray}$$ If 
$l=0$, 
$n>1$ or 
$l\geqslant 3$, the property (P2) and (3.8) show that 
$T_{lmn}^{(0)}=0$, meaning that
$$\begin{eqnarray}f_{lmn}^{(1)}=0,\quad \text{if }l\geqslant 3\text{ or }l=0.\end{eqnarray}$$ In the following, we focus on the cases 
$l=1$ and 
$l=2$. Note that 
$f_{1m0}=0$, by which we can rewrite (3.12) as
$$\begin{eqnarray}T_{1mn}^{(0)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=1}^{+\infty }a_{1nn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}f_{1mn^{\prime }}^{(1)},\quad T_{2mn}^{(0)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{2nn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}f_{2mn^{\prime }}^{(1)}.\end{eqnarray}$$ Again by the property (P2) and (3.8), we see that 
$T_{1mn}^{(0)}=0$ for all 
$n>1$ and 
$T_{2mn}^{(0)}=0$ for all 
$n>0$. For any given 
$m$, the values of 
$f_{1mn^{\prime }}^{(1)}$ and 
$f_{2mn^{\prime }}^{(1)}$ can be solved from (3.14), and the result has the form
$$\begin{eqnarray}f_{1mn}^{(1)}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}A_{1n}\unicode[STIX]{x1D703}^{n-1}T_{1m1}^{(0)},\quad f_{2mn}^{(1)}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}A_{2n}\unicode[STIX]{x1D703}^{n}T_{2m0}^{(0)},\end{eqnarray}$$ where 
$A_{1n}$ and 
$A_{2n}$ are pure numbers. In principle, obtaining (3.15) requires solving an infinite matrix. In our implementation, this is approximated by a cutoff of the right-hand sides of (3.14) up to 
$n^{\prime }\leqslant 9$.
Until now, our calculation is completely the same as the classical Chapman–Enskog expansion. In (3.15), all first-order quantities can be represented by the conservative quantities and their derivatives (hidden in the expressions for 
$T_{1m1}^{(0)}$ and 
$T_{2m0}^{(0)}$). However, to derive 13-moment equations, we are required to represent the distribution functions using more moments and less derivatives. For example, in Grad’s 13-moment theory, no derivatives are included in the ansatz of the distribution function. In our derivation, this can be achieved by writing (3.15) as
$$\begin{eqnarray}f_{1mn}^{(1)}=\frac{A_{1n}}{A_{11}}\unicode[STIX]{x1D703}^{n-1}f_{1m1}^{(1)},\quad f_{2mn}^{(1)}=\frac{A_{2n}}{A_{20}}\unicode[STIX]{x1D703}^{n}f_{2m0}^{(1)}.\end{eqnarray}$$ Since 
$f_{2m0}$ and 
$f_{1m1}$ are included in the 13-moment theory, we mimic the assumption of Chapman–Enskog expansion (3.8) and set
$$\begin{eqnarray}f_{1m1}=\unicode[STIX]{x1D716}f_{1m1}^{(1)},\quad f_{2m0}=\unicode[STIX]{x1D716}f_{2m0}^{(1)},\quad f_{1m1}^{(k)}=f_{2m0}^{(k)}=0,\quad k\geqslant 2.\end{eqnarray}$$ By (3.16) and (3.17), we can write down the approximation of the distribution function up to first order in 
$\unicode[STIX]{x1D716}$ using the 13 moments
$$\begin{eqnarray}\displaystyle f(\boldsymbol{x},\unicode[STIX]{x1D743},t) & {\approx} & \displaystyle {\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743},t)+\mathop{\sum }_{m=-1}^{1}f_{1m1}(\boldsymbol{x},t)\mathop{\sum }_{n=1}^{+\infty }\frac{A_{1n}}{A_{11}}[\unicode[STIX]{x1D703}(\boldsymbol{x},t)]^{n-1}\unicode[STIX]{x1D713}_{1mn}(\boldsymbol{x},\unicode[STIX]{x1D743},t)\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{m=-2}^{2}f_{2m0}(\boldsymbol{x},t)\mathop{\sum }_{n=0}^{+\infty }\frac{A_{2n}}{A_{20}}[\unicode[STIX]{x1D703}(\boldsymbol{x},t)]^{n}\unicode[STIX]{x1D713}_{2mn}(\boldsymbol{x},\unicode[STIX]{x1D743},t),\end{eqnarray}$$ where we have written all the parameters 
$t,\boldsymbol{x}$ and 
$\unicode[STIX]{x1D743}$ for clarification. Note that, in contrast to the Chapman–Enskog expansion of the distribution function, no derivatives are involved in the above expression. Equation (3.18) is also different from Grad’s ansatz for the 13-moment theory, since (3.18) can be written equivalently as
$$\begin{eqnarray}\displaystyle & & \displaystyle f(\unicode[STIX]{x1D743})\approx \left[1-\mathop{\sum }_{i=1}^{3}\frac{(\unicode[STIX]{x1D709}_{i}-v_{i})q_{i}}{\unicode[STIX]{x1D703}^{2}}\mathop{\sum }_{n=1}^{+\infty }\frac{A_{1n}}{A_{11}}\sqrt{\frac{3\unicode[STIX]{x03C0}^{1/2}n!}{10\unicode[STIX]{x1D6E4}(n+5/2)}}L_{n}^{(3/2)}\left(\frac{\unicode[STIX]{x1D743}-\boldsymbol{v}}{2\unicode[STIX]{x1D703}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\mathop{\sum }_{i=1}^{3}\mathop{\sum }_{j=1}^{3}\frac{\unicode[STIX]{x1D70E}_{ij}(\unicode[STIX]{x1D709}_{i}-v_{i})(\unicode[STIX]{x1D709}_{j}-v_{j})}{\unicode[STIX]{x1D703}^{2}}\mathop{\sum }_{n=0}^{+\infty }\frac{A_{2n}}{A_{20}}\sqrt{\frac{15\unicode[STIX]{x03C0}^{1/2}n!}{32\unicode[STIX]{x1D6E4}(n+7/2)}}L_{n}^{(5/2)}\left(\frac{\unicode[STIX]{x1D743}-\boldsymbol{v}}{2\unicode[STIX]{x1D703}}\right)\right]{\mathcal{M}}(\unicode[STIX]{x1D743}).\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ Here, 
$L_{n}^{(\unicode[STIX]{x1D6FC})}$ is the Laguerre polynomial defined in (B 2), and we have omitted the variables 
$t$ and 
$\boldsymbol{x}$ for conciseness. Grad’s ansatz for the 13-moment theory can be obtained by truncating both infinite sums in (3.19) by preserving only their first terms, and therefore Grad’s 13-moment theory does not fully represent the first-order term of the distribution function for general collision models. Due to the assumption (3.17), when we apply (3.11), equations for all the 13 moments should be excluded, i.e. 
$(l,n)\neq (0,0),(0,1),(1,0),(1,1),(2,0)$.
Remark 3. The solvability of (3.14) relies on the existence of the spectral gap for the linearized Boltzmann collision operator. For IPL models, this has been proven in Mouhot & Strain (Reference Mouhot and Strain2007). In particular, when 
$\unicode[STIX]{x1D702}=5$, we have 
$A_{1n}=\unicode[STIX]{x1D6FF}_{1n}/a_{111}$ and 
$A_{2n}=\unicode[STIX]{x1D6FF}_{0n}/a_{200}$. In this case, the first two orders (3.18) are exactly the ansatz in Grad’s 13-moment theory.
3.2.2 Second order ($k=1$)
Now we set 
$k=1$ in (3.11). Since 
$f_{000}=f_{001}=f_{1m0}=0$ and 
$f_{2m0}^{(k)}=f_{1m1}^{(k)}=0$ for 
$k\geqslant 2$, the result can be written as
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}f_{lmn}^{(1)}}{\unicode[STIX]{x2202}t}}+S_{lmn}^{(1)}+T_{lmn}^{(1)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n^{\prime }=n_{0}(l)}^{+\infty }a_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}f_{lmn^{\prime }}^{(2)},\end{eqnarray}$$where
$$\begin{eqnarray}n_{0}(l)=\left\{\begin{array}{@{}ll@{}}2, & \text{if }l=0,1,\\ 1, & \text{if }l=2,\\ 0, & \text{if }l\geqslant 3.\end{array}\right.\end{eqnarray}$$ Since 
$f_{lmn}^{(1)}$ has been fully obtained in the previous section (see (3.13) and (3.17)), the expressions of 
$S_{lmn}^{(1)}$ and 
$T_{lmn}^{(1)}$ can be naturally obtained. Thus, the left-hand side of (3.20) can already be represented by the 13 moments. To obtain the second-order contributions 
$f_{lmn}^{(2)}$, for any 
$l$ and 
$m$, we just need to solve the infinite linear system (3.20), and the general result is
$$\begin{eqnarray}f_{lmn}^{(2)}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}\mathop{\sum }_{n^{\prime }=n_{0}(l)}^{+\infty }B_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}\left({\displaystyle \frac{\unicode[STIX]{x2202}f_{lmn^{\prime }}^{(1)}}{\unicode[STIX]{x2202}t}}+S_{lmn^{\prime }}^{(1)}+T_{lmn^{\prime }}^{(1)}\right),\end{eqnarray}$$ where 
$B_{lnn^{\prime }}$ are all pure numbers. In our implementation, we again truncate the system (3.20) at 
$n^{\prime }=\lfloor 10-l/2\rfloor$. For the purpose of deriving R13 equations, we only need 
$f_{lmn}^{(2)}$ up to 
$l=4$.
The right-hand side of (3.22) still contains time derivatives, which are not desired. In general, they can all be replaced by spatial derivatives. Note that 
$\unicode[STIX]{x2202}_{t}\,f_{lmn}^{(1)}$ is non-zero only if 
$l=1$ or 
$l=2$. In these cases, we have
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}f_{1mn^{\prime }}^{(1)}}{\unicode[STIX]{x2202}t}=\frac{A_{1n^{\prime }}}{A_{11}}\left(\unicode[STIX]{x1D703}^{n^{\prime }-1}\frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}+(n^{\prime }-1)\unicode[STIX]{x1D703}^{n^{\prime }-2}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}f_{1m1}^{(1)}\right),\\ \displaystyle \frac{\unicode[STIX]{x2202}f_{2mn^{\prime }}^{(1)}}{\unicode[STIX]{x2202}t}=\frac{A_{2n^{\prime }}}{A_{20}}\left(\unicode[STIX]{x1D703}^{n^{\prime }}\frac{\unicode[STIX]{x2202}f_{2m0}^{(1)}}{\unicode[STIX]{x2202}t}+n^{\prime }\unicode[STIX]{x1D703}^{n^{\prime }-1}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}f_{2m0}^{(1)}\right).\end{array}\right\}\end{eqnarray}$$ We first focus on 
$\unicode[STIX]{x2202}_{t}\,f_{1m1}^{(1)}$ and 
$\unicode[STIX]{x2202}_{t}\,f_{2m0}^{(1)}$. Taking 
$\unicode[STIX]{x2202}_{t}\,f_{1m1}^{(1)}$ as an example, we perform the following calculation:
from which it can be solved that
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t} & = & \displaystyle \left(1-\mathop{\sum }_{n^{\prime }=2}^{+\infty }\unicode[STIX]{x1D706}_{1n^{\prime }}\right)^{-1}\left[\frac{1}{\unicode[STIX]{x1D716}}\left(\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\frac{f_{1m1}^{(1)}}{A_{11}}-T_{1m1}^{(0)}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\mathop{\sum }_{n^{\prime }=2}^{+\infty }\frac{(n^{\prime }-1)\unicode[STIX]{x1D706}_{1n^{\prime }}}{\unicode[STIX]{x1D703}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}f_{1m1}^{(1)}+\mathop{\sum }_{n^{\prime }=1}^{+\infty }C_{1n^{\prime }}\unicode[STIX]{x1D703}^{1-n^{\prime }}(S_{1mn^{\prime }}^{(1)}+T_{1mn^{\prime }}^{(1)})\right]+O(\unicode[STIX]{x1D716}),\qquad\end{eqnarray}$$ where we have used 
$S_{1m1}^{(0)}=0$, and the newly introduced constants are
$$\begin{eqnarray}C_{ln^{\prime }}=\left\{\begin{array}{@{}ll@{}}-1, & \text{if }n^{\prime }=n_{0}(l)-1,\\ \displaystyle \mathop{\sum }_{n=n_{0}(l)}^{+\infty }a_{l,n_{0}(l)-1,n}B_{lnn^{\prime }}, & \text{if }n^{\prime }\geqslant n_{0}(l),\end{array}\right.\quad \unicode[STIX]{x1D706}_{ln^{\prime }}=\frac{C_{ln^{\prime }}A_{ln^{\prime }}}{A_{l,n_{0}(l)-1}}.\end{eqnarray}$$Similarly, we can obtain
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f_{2m0}^{(1)}}{\unicode[STIX]{x2202}t} & = & \displaystyle \left(1-\mathop{\sum }_{n^{\prime }=1}^{+\infty }\unicode[STIX]{x1D706}_{2n^{\prime }}\right)^{-1}\left[\frac{1}{\unicode[STIX]{x1D716}}\left(\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\frac{f_{2m0}^{(1)}}{A_{20}}-T_{2m0}^{(0)}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\mathop{\sum }_{n^{\prime }=1}^{+\infty }\frac{n^{\prime }\unicode[STIX]{x1D706}_{2n^{\prime }}}{\unicode[STIX]{x1D703}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}f_{2m0}^{(1)}+\mathop{\sum }_{n^{\prime }=0}^{+\infty }C_{2n^{\prime }}\unicode[STIX]{x1D703}^{-n^{\prime }}(S_{2mn^{\prime }}^{(1)}+T_{2mn^{\prime }}^{(1)})\right]+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$As a summary, the time derivatives in (3.22) can be replaced with spatial derivatives by the following operations.
(i) If
$l\neq 1$ and $l\neq 2$, set the derivative $\unicode[STIX]{x2202}_{t}\,f_{lmn^{\prime }}^{(1)}$ to zero. If $l=1$ or $l=2$, replace $\unicode[STIX]{x2202}_{t}\,f_{1mn^{\prime }}^{(1)}$ by (3.23), (3.25) and (3.27).(ii) After replacement, the results still include the time derivatives of
$\boldsymbol{v}$ and $\unicode[STIX]{x1D703}$. These terms can be replaced by conservation laws (2.2). In fact, we can use the fact that $\unicode[STIX]{x1D70E}_{kl}$ and $q_{k}$ are $O(\unicode[STIX]{x1D716})$ terms to rewrite the conservation laws of momentum and energy as (3.28a,b)and use these equations for substitution.$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}t}=-v_{k}{\displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{k}}}-\frac{\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D70C}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}}+O(\unicode[STIX]{x1D716}),\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}=-\frac{2}{3}\unicode[STIX]{x1D703}{\displaystyle \frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}}-v_{k}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}}+O(\unicode[STIX]{x1D716}),\end{eqnarray}$$(iii) After the above replacements, the
$O(\unicode[STIX]{x1D716})$ terms can be dropped.
By now, we have presented all the coefficients 
$f_{lmn}^{(2)}$ by the 13 moments and their spatial derivatives. When 
$l=1,2$, the results include a coefficient 
$1/\unicode[STIX]{x1D716}$ coming from (3.25) and (3.27). More precisely, 
$f_{1mn}^{(2)}$ and 
$f_{2mn}^{(2)}$ have the form
$$\begin{eqnarray}\displaystyle & \displaystyle f_{1mn}^{(2)}=\frac{1}{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D703}^{n-1}\left(\mathop{\sum }_{n^{\prime }=2}^{+\infty }B_{1nn^{\prime }}\frac{A_{1n^{\prime }}}{A_{11}}\right)\left(1-\mathop{\sum }_{n^{\prime }=2}^{+\infty }\unicode[STIX]{x1D706}_{1n^{\prime }}\right)^{-1}\left(\frac{f_{1m1}^{(1)}}{A_{11}}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}T_{1m1}^{(0)}\right)+W_{1mn}^{(2)},\quad n\geqslant 2, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle f_{2mn}^{(2)}=\frac{1}{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D703}^{n}\left(\mathop{\sum }_{n^{\prime }=1}^{+\infty }B_{2nn^{\prime }}\frac{A_{2n^{\prime }}}{A_{20}}\right)\left(1-\mathop{\sum }_{n^{\prime }=1}^{+\infty }\unicode[STIX]{x1D706}_{2n^{\prime }}\right)^{-1}\left(\frac{f_{2m0}^{(1)}}{A_{20}}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}T_{2m0}^{(0)}\right)+W_{2mn}^{(2)},\quad n\geqslant 1, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ where 
$W_{1mn}^{(2)}$ and 
$W_{2mn}^{(2)}$ are terms independent of 
$\unicode[STIX]{x1D716}$. Note that only first-order derivatives have been introduced into 
$f_{lmn}^{(2)}$, while in the original Chapman–Enskog expansion, the second-order (Burnett-order) term 
$f^{(2)}$ includes second-order derivatives.
3.2.3 Third order ($k=2$)
Similar to the case 
$k=1$, when 
$k=2$, we can solve the linear system (3.11) to obtain
$$\begin{eqnarray}f_{lmn}^{(3)}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}\mathop{\sum }_{n^{\prime }=n_{0}(l)}^{+\infty }B_{lnn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}\left({\displaystyle \frac{\unicode[STIX]{x2202}f_{lmn^{\prime }}^{(2)}}{\unicode[STIX]{x2202}t}}+S_{lmn^{\prime }}^{(2)}+T_{lmn^{\prime }}^{(2)}\right),\quad n\geqslant n_{0}(l).\end{eqnarray}$$ After inserting the expression for 
$f_{lmn}^{(2)}$ into the above equation, we again need to deal with the time derivatives. The time derivatives appearing in 
$S_{lmn^{\prime }}^{(2)}$ can again be replaced by (3.28) without 
$O(\unicode[STIX]{x1D716})$ terms. In the following we focus only the time derivative of 
$f_{lmn^{\prime }}^{(2)}$.
When 
$l\neq 1$ and 
$l\neq 2$, the second-order term 
$f_{lmn^{\prime }}^{(2)}$ does not include the coefficient 
$1/\unicode[STIX]{x1D716}$, and therefore 
$\unicode[STIX]{x2202}_{t}\,f_{lmn^{\prime }}^{(2)}$ can be computed by inserting the expression for 
$f_{lmn^{\prime }}^{(2)}$, expanding the time derivative and then replacing the time derivative of each moment in (3.3) by (3.28), (3.25) and (3.27) and the continuity equation 
$\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}=-\text{div}(\unicode[STIX]{x1D70C}\boldsymbol{v})$. After replacement, we can also safely drop the 
$O(\unicode[STIX]{x1D716})$ terms. The treatment for 
$l=1,2$ is more complicated, and the process is described in detail.
When 
$l=1$ and 
$n\geqslant 2$ (the case 
$l=2,n\geqslant 1$ is similar), by (3.29),
$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{1mn^{\prime }}^{(2)}}{\unicode[STIX]{x2202}t}} & = & \displaystyle \frac{1}{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D703}^{n^{\prime }-1}\frac{D_{1n^{\prime }}}{A_{11}}{\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}}-\frac{1}{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D703}^{n^{\prime }-1}D_{1n^{\prime }}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\left(\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}T_{1m1}^{(0)}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{\unicode[STIX]{x1D716}}(n^{\prime }-1)\unicode[STIX]{x1D703}^{n^{\prime }-2}D_{1n^{\prime }}\left(\frac{f_{1m1}^{(1)}}{A_{11}}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}T_{1m1}^{(0)}\right){\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{\unicode[STIX]{x2202}W_{1mn^{\prime }}^{(2)}}{\unicode[STIX]{x2202}t}},\end{eqnarray}$$where
$$\begin{eqnarray}D_{1n}=\left(\mathop{\sum }_{n^{\prime }=2}^{+\infty }B_{1nn^{\prime }}\frac{A_{1n^{\prime }}}{A_{11}}\right)\left(1-\mathop{\sum }_{n^{\prime }=2}^{+\infty }\unicode[STIX]{x1D706}_{1n^{\prime }}\right)^{-1}.\end{eqnarray}$$ The last term 
$\unicode[STIX]{x2202}_{t}W_{1mn^{\prime }}^{(2)}$ is independent of 
$\unicode[STIX]{x1D716}$, and therefore the time derivative can also be replaced by conservation laws and (3.25) and (3.27) without 
$O(\unicode[STIX]{x1D716})$ terms. However, for the first three terms on the right-hand side of (3.32), due to the existence of 
$1/\unicode[STIX]{x1D716}$, when replaced by spatial derivatives, the 
$O(\unicode[STIX]{x1D716})$ terms have to be taken into account to capture the 
$O(1)$ contribution. Note that 
$T_{1m1}^{(0)}$ also appears in (3.15), where the first equation for 
$n=1$ represents Fourier’s law. We know that 
$(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703})^{-1}T_{1m1}^{(0)}$ is essentially the spatial derivative of 
$\unicode[STIX]{x1D703}$ multiplied by a constant. Therefore, the second and third terms involve only the time derivative of density 
$\unicode[STIX]{x1D70C}$ and temperature 
$\unicode[STIX]{x1D703}$, and they can be replaced exactly by the conservation laws (2.2). Thus, the only troublesome term is again 
$\unicode[STIX]{x2202}_{t}\,f_{1m1}^{(1)}$. Before discussing this term, we first implement all the aforementioned replacements and write the result as
$$\begin{eqnarray}f_{1mn}^{(3)}=\frac{1}{\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x1D703}^{n-1}}{A_{11}}\left(\mathop{\sum }_{n^{\prime }=2}^{+\infty }B_{1nn^{\prime }}D_{1n^{\prime }}\right){\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}}+R_{1mn}^{(3)},\end{eqnarray}$$ where 
$R_{1mn}^{(3)}$ is the collection of terms which does not include any time derivatives.
By now, we can carry out a calculation similar to (3.24)
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D716}}(S_{1m1}^{(0)}+T_{1m1}^{(0)})+(S_{1m1}^{(1)}+T_{1m1}^{(1)})+\unicode[STIX]{x1D716}(S_{1m1}^{(2)}+T_{1m1}^{(2)})\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=1}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}f_{1mn}^{(1)}+\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=1}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}f_{1mn}^{(2)}+\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=1}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}f_{1mn}^{(3)}+O(\unicode[STIX]{x1D716}^{2})\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{\unicode[STIX]{x1D716}}\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\frac{f_{1m1}^{(1)}}{A_{11}}+\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=2}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}f_{1mn}^{(2)}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\mathop{\sum }_{n=2}^{+\infty }\frac{a_{11n}}{A_{11}}\mathop{\sum }_{n^{\prime }=2}^{+\infty }B_{1nn^{\prime }}D_{1n^{\prime }}{\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D716}\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=2}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}R_{1mn}^{(3)}+O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$ Solving 
$\unicode[STIX]{x2202}_{t}\,f_{1m1}^{(1)}$ from the above equation, we obtain
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t} & = & \displaystyle \left(1-\frac{1}{A_{11}}\mathop{\sum }_{n^{\prime }=2}^{+\infty }C_{1n^{\prime }}D_{1n^{\prime }}\right)^{-1}\left[\frac{1}{\unicode[STIX]{x1D716}}\left(\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\frac{f_{1m1}^{(1)}}{A_{11}}-T_{1m1}^{(0)}\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left(\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=2}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}f_{1mn}^{(2)}-(S_{1m1}^{(1)}+T_{1m1}^{(1)})\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\left.\unicode[STIX]{x1D716}\left(\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\mathop{\sum }_{n=2}^{+\infty }a_{11n}\unicode[STIX]{x1D703}^{1-n}R_{1mn}^{(3)}-(S_{1m1}^{(2)}+T_{1m1}^{(2)})\right)\right]+O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$ In (3.36), the time derivatives in 
$S_{1m1}^{(2)}$ can be replaced by (3.28) with 
$O(\unicode[STIX]{x1D716})$ terms dropped, while the time derivatives in 
$S_{1m1}^{(1)}$ have to be replaced by complete conservation laws (2.2). The last step is to substitute the above equation into (3.34), and the 
$O(\unicode[STIX]{x1D716}^{2})$ term in (3.36) can now be discarded. This completes the calculation of 
$f_{1mn}^{(3)}$.
The calculation of 
$f_{2mn}^{(3)}$ follows exactly the same procedure, and the details are omitted. The above procedure shows that the expression for 
$f_{lmn}^{(3)}$ includes second-order derivatives of the 13 moments.
3.3 The 13-moment equations
With all the moments up to third order calculated, we are ready to write down the 13-moment equations. Note that all the 13-moment equations include the conservation laws (2.2). Therefore, we focus only on the equations for 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{j}$ or, equivalently, 
$f_{2m0}$ and 
$f_{1m1}$. Since 
$f_{2m0}=\unicode[STIX]{x1D716}f_{2m0}^{(1)}$ and 
$f_{1m1}=\unicode[STIX]{x1D716}f_{1m1}^{(1)}$, these equations are to be obtained by truncation of (3.10). Two different truncations are considered in the following, which correspond to GG13 equations and R13 equations, respectively.
3.3.1 GG13-moment equations
The derivation of GG13 equations is basically a truncation of (3.10) up to the first order. The result reads
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D716}{\displaystyle \frac{\unicode[STIX]{x2202}f_{2m0}^{(1)}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D716}S_{2m0}^{(1)}+T_{2m0}^{(0)}+\unicode[STIX]{x1D716}T_{2m0}^{(1)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\left(\frac{1}{A_{20}}f_{2m0}^{(1)}+\unicode[STIX]{x1D716}\mathop{\sum }_{n^{\prime }=1}^{+\infty }a_{20n^{\prime }}\unicode[STIX]{x1D703}^{-n^{\prime }}f_{2mn^{\prime }}^{(2)}\right),\\ \displaystyle \unicode[STIX]{x1D716}{\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D716}S_{1m1}^{(1)}+T_{1m1}^{(0)}+\unicode[STIX]{x1D716}T_{1m1}^{(1)}=\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\left(\frac{1}{A_{11}}f_{1m1}^{(1)}+\unicode[STIX]{x1D716}\mathop{\sum }_{n^{\prime }=2}^{+\infty }a_{11n^{\prime }}\unicode[STIX]{x1D703}^{1-n^{\prime }}f_{1mn^{\prime }}^{(2)}\right),\end{array}\right\}\end{eqnarray}$$where we have used
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{20n^{\prime }}\unicode[STIX]{x1D703}^{-n^{\prime }}f_{2mn^{\prime }}^{(1)}=\mathop{\sum }_{n^{\prime }=0}^{+\infty }a_{20n^{\prime }}\frac{A_{2n^{\prime }}}{A_{20}}f_{2m0}^{(1)}=\frac{1}{A_{20}}f_{2m0}^{(1)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{n^{\prime }=1}^{+\infty }a_{11n^{\prime }}\unicode[STIX]{x1D703}^{1-n^{\prime }}f_{1mn^{\prime }}^{(1)}=\mathop{\sum }_{n^{\prime }=1}^{+\infty }a_{11n^{\prime }}\frac{A_{1n^{\prime }}}{A_{11}}f_{1m1}^{(1)}=\frac{1}{A_{11}}f_{1m1}^{(1)}. & \displaystyle\end{eqnarray}$$ All the terms in (3.37) have been represented by the 13 moments (3.3) and their derivatives in § 3.2, and the time derivatives in 
$S_{2m0}^{(1)}$ can be replaced by (3.28) with 
$O(\unicode[STIX]{x1D716})$ terms discarded. The final step is to revert to the scaling of space and time introduced in the beginning of § 3.2, which can be simply achieved by setting 
$\unicode[STIX]{x1D716}$ to 
$1$.
This set of equations is called the GG13-moment equations, as proposed in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013). Similar to Grad’s 13-moment equations, the GG13 equations are also first-order quasi-linear equations. However, they are different from Grad’s 13-moment equations in two ways: (i) in Grad’s 13-moment equations, the terms 
$T_{1m1}^{(1)}$ and 
$T_{2m0}^{(1)}$ come directly from the truncation of the distribution function, while in the GG13-moment equations, they include information from inversion of the linearized collision operator; (ii) in Grad’s 13-moment equations, the collision term does not include information from 
$f_{1mn^{\prime }}$ with 
$n^{\prime }>1$ or 
$f_{2mn^{\prime }}$ with 
$n^{\prime }>0$. Due to such differences, as mentioned in § 2.3, Grad’s 13-moment equations are accurate only up to first order for general IPL models.
Now we compare equations (3.37) with equations (2.3). On the left-hand side, we can see from (3.37) that no second-order terms are involved. Since the moments 
$m_{ijk}$ (
$f_{3m0}$), 
$R_{ik}$ (
$f_{2m1}$) and 
$\unicode[STIX]{x1D6E5}$ (
$f_{002}$) are all 
$O(\unicode[STIX]{x1D716}^{2})$ terms, they are simply set to zero in the GG13 equations. The right-hand side is much more complicated owing to the involved formulae for second-order terms. We would just like to point out that the coefficient 
$D_{0}^{(\unicode[STIX]{x1D702})}$ in (2.4) does not equal the coefficient 
$1/A_{20}$ in (3.37), since the same term also appears in 
$f_{2mn}^{(2)}$, as is shown in (3.30). For a similar reason, the coefficient 
$E_{0}^{(\unicode[STIX]{x1D702})}$ in (2.4) does not equal 
$1/A_{11}$ in (3.37).
3.3.2 R13-moment equations
To gain one more order of accuracy, we need to keep the second-order terms in (3.10), and the result is
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \begin{array}{@{}rcl@{}}\; & \; & \displaystyle \unicode[STIX]{x1D716}{\displaystyle \frac{\unicode[STIX]{x2202}f_{2m0}^{(1)}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D716}S_{2m0}^{(1)}+\unicode[STIX]{x1D716}^{2}S_{2m0}^{(2)}+T_{2m0}^{(0)}+\unicode[STIX]{x1D716}T_{2m0}^{(1)}+\unicode[STIX]{x1D716}^{2}T_{2m0}^{(2)}\\ \; & \; & \displaystyle \quad =\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\left(\frac{1}{A_{20}}f_{2m0}^{(1)}+\mathop{\sum }_{n^{\prime }=1}a_{2nn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}(\unicode[STIX]{x1D716}f_{2mn^{\prime }}^{(2)}+\unicode[STIX]{x1D716}^{2}f_{2mn^{\prime }}^{(3)})\right),\end{array}\\ \displaystyle \begin{array}{@{}rcl@{}}\; & \; & \displaystyle \unicode[STIX]{x1D716}{\displaystyle \frac{\unicode[STIX]{x2202}f_{1m1}^{(1)}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D716}S_{1m1}^{(1)}+\unicode[STIX]{x1D716}^{2}S_{1m1}^{(2)}+T_{1m1}^{(0)}+\unicode[STIX]{x1D716}T_{1m1}^{(1)}+\unicode[STIX]{x1D716}^{2}T_{1m1}^{(2)}\\ \; & \; & \displaystyle \quad =\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D707}}\left(\frac{1}{A_{11}}f_{1m1}^{(1)}+\mathop{\sum }_{n^{\prime }=2}a_{1nn^{\prime }}\unicode[STIX]{x1D703}^{n-n^{\prime }}(\unicode[STIX]{x1D716}f_{1mn^{\prime }}^{(2)}+\unicode[STIX]{x1D716}^{2}f_{1mn^{\prime }}^{(3)})\right).\end{array}\end{array}\right\}\end{eqnarray}$$ Again, the time derivatives for velocity and temperature in 
$S_{2m0}^{(i)}$ and 
$S_{1m1}^{(i)},i=1,2$ need to be replaced by spatial derivatives. In order to preserve the second-order terms, the time derivatives in 
$S_{2m0}^{(1)}$ and 
$S_{1m1}^{(1)}$ need to be substituted by the complete conservation laws (2.2), where, as in 
$S_{2m0}^{(2)}$ and 
$S_{1m1}^{(2)}$, the replacement of time derivatives is done by using (3.28) and discarding 
$O(\unicode[STIX]{x1D716})$ terms. Afterwards, we set 
$\unicode[STIX]{x1D716}$ to 
$1$, and the result is the R13-moment equations. The Mathematica code for the deduction of the R13-moment equations is available in the supplementary material, which is available at https://doi.org/10.1017/jfm.2020.251.
By replacing all the coefficients with the primitive variables 
$\unicode[STIX]{x1D70E}_{ij}$ and 
$q_{i}$, equations (2.3) with (2.4) and (2.6) and (C 1) and (C 2) can be obtained. Compared with the linear R13 equations obtained in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2013), much more information is included in this nonlinear version. For instance, in (2.6), one can see that all of the first three terms in 
$m_{ijk}^{(\unicode[STIX]{x1D702})}$ are nonlinear and all of the last three terms in 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ are nonlinear. Clearly, these terms cannot be ignored in a problem such as the structure of plane shock waves, which will be studied in the following section.
4 Numerical examples
In this section, we test the behaviour of the nonlinear R13 equations by computing the structure of one-dimensional plane shock waves, which is a benchmark problem in gas kinetic theory. It involves strong non-equilibrium, but does not have any boundary condition, which makes it suitable for testing the ability to describe non-equilibrium processes for our models.
For one-dimensional flow, the moments satisfy
$$\begin{eqnarray}v_{2}=v_{3}=0,\quad \unicode[STIX]{x1D70E}_{12}=\unicode[STIX]{x1D70E}_{13}=\unicode[STIX]{x1D70E}_{23}=0,\quad \unicode[STIX]{x1D70E}_{22}=\unicode[STIX]{x1D70E}_{33}=-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70E}_{11},\quad q_{2}=q_{3}=0.\end{eqnarray}$$ For simplicity, we use the notation 
$v=v_{1}$, 
$\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{11}$ and 
$q=q_{1}$. Then 13 moments are reduced to five moments 
$\unicode[STIX]{x1D70C}$, 
$v$, 
$\unicode[STIX]{x1D703}$, 
$\unicode[STIX]{x1D70E}$ and 
$q$. In the following we write the R13 model for these quantities in the form of balance laws
where
$$\begin{eqnarray}\displaystyle \text{FQ} & = & \displaystyle \frac{16}{5}vq+\frac{1}{2}\unicode[STIX]{x1D70C}v^{4}+4v^{2}\unicode[STIX]{x1D703}+\frac{5}{2}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}+\frac{5}{2}v^{2}\unicode[STIX]{x1D70E}+\left(\frac{7}{2}-\sqrt{\frac{14}{3}}C_{1D}^{(\unicode[STIX]{x1D702})}\right)\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}\nonumber\\ \displaystyle & & \displaystyle +\,m^{(\unicode[STIX]{x1D702})}v+\frac{1}{2}R^{(\unicode[STIX]{x1D702})}+\frac{1}{6}\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})},\end{eqnarray}$$ with 
$m^{(\unicode[STIX]{x1D702})},R^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ being 
$m_{111}^{(\unicode[STIX]{x1D702})},R_{11}^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6E5}^{(\unicode[STIX]{x1D702})}$ in (2.6) substituted by (4.1). On the right-hand side, 
$\unicode[STIX]{x1D6F4}_{1D}^{(\unicode[STIX]{x1D702})}$ and 
$Q_{1D}^{(\unicode[STIX]{x1D702})}$ are, respectively, 
$\unicode[STIX]{x1D6F4}_{11}^{(\unicode[STIX]{x1D702},1)}+\unicode[STIX]{x1D6F4}_{11}^{(\unicode[STIX]{x1D702},2)}$ and 
$Q_{1}^{(\unicode[STIX]{x1D702},1)}+Q_{1}^{(\unicode[STIX]{x1D702},2)}$ in (2.3) subject to (4.1). The constant 
$C_{1D}^{(\unicode[STIX]{x1D702})}$ depends only on 
$\unicode[STIX]{x1D702}$, and some of its values are listed in table 5.
Table 5. Coefficient 
$C_{1D}^{(\unicode[STIX]{x1D702})}$ for different 
$\unicode[STIX]{x1D702}$.
For simplicity, we adopt the dimensionless setting in our numerical tests. Let 
$\unicode[STIX]{x1D70C}_{0}$ and 
$\unicode[STIX]{x1D703}_{0}$ be, respectively, the density and temperature in front of the shock. We scale the moments as follows:
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{0},\quad \hat{v}=v/\sqrt{\unicode[STIX]{x1D703}},\quad \hat{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D703}/\unicode[STIX]{x1D703}_{0},\quad \hat{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70E}/(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D703}_{0}),\quad \hat{q}=q/(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D703}_{0}^{3/2}).\end{eqnarray}$$ The non-dimensionalization of the viscosity 
$\unicode[STIX]{x1D707}$ is 
$\hat{\unicode[STIX]{x1D707}}/\unicode[STIX]{x1D707}_{0}$ with 
$\unicode[STIX]{x1D707}_{0}$ being the viscosity at density 
$\unicode[STIX]{x1D70C}_{0}$ and temperature 
$\unicode[STIX]{x1D703}_{0}$ computed using (2.1). The dimensionless time and spatial variables are given by 
$\hat{x}=x/\unicode[STIX]{x1D706}_{0}$ and 
$\hat{t}=t/(\unicode[STIX]{x1D706}_{0}/\sqrt{\unicode[STIX]{x1D703}_{0}})$, where 
$\unicode[STIX]{x1D706}_{0}$ is the mean free path defined by
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{0}=\frac{8}{15}\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{(\unicode[STIX]{x1D702}-2)(3\unicode[STIX]{x1D702}-5)}{(\unicode[STIX]{x1D702}-1)^{2}}\frac{\unicode[STIX]{x1D707}_{0}}{\unicode[STIX]{x1D70C}_{0}\sqrt{\unicode[STIX]{x1D703}_{0}}}.\end{eqnarray}$$ After non-dimensionalization, we remove the ‘hats’ on the variables, and the resulting equations still have the form (4.2), with 
$\unicode[STIX]{x1D707}$ defined by
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{15}{8}\sqrt{\frac{\unicode[STIX]{x03C0}}{2}}\frac{(\unicode[STIX]{x1D702}-1)^{2}}{(\unicode[STIX]{x1D702}-2)(3\unicode[STIX]{x1D702}-5)}\unicode[STIX]{x1D703}^{(\unicode[STIX]{x1D702}-5)/(2(\unicode[STIX]{x1D702}-1))}.\end{eqnarray}$$ Based on this dimensionless setting, the structure of plane shock waves with Mach number 
$Ma$ can be obtained by choosing the initial data to be
$$\begin{eqnarray}(\unicode[STIX]{x1D70C},v,\unicode[STIX]{x1D703},\unicode[STIX]{x1D70E},q)=\left\{\begin{array}{@{}ll@{}}(\unicode[STIX]{x1D70C}_{l},v_{l},\unicode[STIX]{x1D703}_{l},0,0), & x<0,\\ (\unicode[STIX]{x1D70C}_{r},v_{r},\unicode[STIX]{x1D703}_{r},0,0), & x>0,\end{array}\right.\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{l}=1,\quad v_{l}=\sqrt{5/3}Ma,\quad \unicode[STIX]{x1D703}_{l}=1, & \displaystyle \nonumber\\ \displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{r}=\frac{4Ma^{2}}{Ma^{2}+3},\quad v_{r}=\sqrt{\frac{5}{3}}\frac{Ma^{2}+3}{4Ma},\quad \unicode[STIX]{x1D703}_{r}=\frac{5Ma^{2}-1}{4\unicode[STIX]{x1D70C}_{r}}. & \displaystyle \nonumber\end{eqnarray}$$To solve (4.2) and (4.4) numerically, the finite volume method is adopted. Since the left-hand side of (4.2) has the form of a conservation law, we apply the Harten–Lax–van Leer (HLL) scheme in the discretization. The right-hand side provides the non-conservative part, for which a central difference method is used to approximate both the first and second derivatives. For the time discretization, we use the classical forward Euler method in all the examples. The direct simulation Monte Carlo (DSMC) results for variable hard-sphere models with the same viscosity index are used as reference solutions (Bird Reference Bird1994).
Figure 5. Normalized density, velocity and temperature of shock structures for the Maxwell molecules model and Mach numbers 
$Ma=1.55,9$. DSMC solutions for the variable hard-sphere model are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path. The solid lines are the numerical results for the linearized collision model and the dashed lines are those for the quadratic collision model.
4.1 Shock structure for Maxwell molecules
This section is devoted to shock structure computation of Maxwell molecules. The same computation has been carried out in Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017), and the main purpose of this section is the verification of our numerical method and the comparison between linearized and quadratic collision terms. Note that the R13 equations for Maxwell molecules with quadratic collision terms are already available (Struchtrup Reference Struchtrup2005b). Two Mach numbers 
$1.55$ and 
$9.0$ are tested, and we show the numerical results for all the moments in figures 5 and 6, where the density, velocity and temperature are given in their normalized form
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D70C}}=\frac{\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{l}}{|\unicode[STIX]{x1D70C}_{r}-\unicode[STIX]{x1D70C}_{l}|},\quad \bar{v}=\frac{v-v_{r}}{|v_{l}-v_{r}|},\quad \bar{\unicode[STIX]{x1D703}}=\frac{\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{l}}{|\unicode[STIX]{x1D703}_{r}-\unicode[STIX]{x1D703}_{l}|}.\end{eqnarray}$$ For the small Mach number 
$1.55$, both linearized and quadratic collision terms provide good agreement with DSMC results, except for a slight underestimation of the peak heat flux. Surprisingly, when the Mach number reaches 
$9.0$, the R13 equations with linearized and quadratic collision terms still provide an almost identical shock structure. Similar to the results in Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017), our profiles also show some typical structures for R13 solutions with high Mach numbers, such as kinks in the profiles and a too fast decay in the low-density region, which indicates the correctness of our simulation. As is stated in Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017), the underpredicted heat flux is due to the loss of some fourth-order terms in the regularized moment equations. Nevertheless, these results indicate that quadratic collision terms do not contribute too much to the shock structure for a Mach number lower than 
$9.0$.
Figure 6. The stress 
$\unicode[STIX]{x1D70E}$ and heat flux 
$q$ of shock structures for the Maxwell molecules model and Mach numbers 
$Ma=1.55,9$. DSMC results for the variable hard-sphere model are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path. The left 
$y$-axis corresponds to the stress and the right 
$y$-axis corresponds to the heat flux. The solid lines are the numerical results for the linearized collision model and the dashed lines are those for the quadratic collision model.
These numerical tests again show how the R13 equations lose their validity as the Mach number grows. For large Mach numbers, an unphysical kink is predicted in the low-density region, and, as shown in Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017), it becomes more obvious as the Mach number grows. Timokhin et al. (Reference Timokhin, Struchtrup, Kokhanchik and Bondar2017) predict 
$Ma\approx 5.5$ as a formal upper boundary of applicability of the R13 system of Maxwell molecules. In the following experiments, we show that the kink-like structure will also appear in other collision models, while the range of applicability may change when the parameter 
$\unicode[STIX]{x1D702}$ varies.
4.2 Shock structures for different Mach numbers
In this experiment, we test the approximability of the R13 model by varying the Mach number for a non-Maxwell collision model. Four Mach numbers 
$Ma=1.55,3.0,6.5,9.0$ are taken into account, and we consider the IPL model with 
$\unicode[STIX]{x1D702}=10$ and the hard-sphere model (
$\unicode[STIX]{x1D702}=\infty$) in our tests.
Figures 7 and 8 show the comparison between the R13 results and the DSMC results for 
$\unicode[STIX]{x1D702}=10$. The profiles of all the five quantities have been plotted. Although DSMC uses a variable hard-sphere model as an approximation of the IPL model, the simulation results in Hu & Cai (Reference Hu and Cai2019) show that, for the shock structure problem, the variable hard-sphere model and the IPL model show almost identical results for Mach numbers 
$6.5$ and 
$9.0$, which means that it is reliable to use DSMC results to check the quality of R13 results. For increasing Mach number, it is generally harder for macroscopic models to accurately capture the non-equilibrium effects. This can be clearly observed from figure 8, which shows that the heat flux is underestimated in the low-density region. Note that the shock structure in the high-density region is well captured for all Knudsen numbers, since the high density and temperature in this region result in distribution functions close to the local Maxwellians, which can be relatively easier to represent using the Chapman–Enskog expansion.
Figures 9 and 10 show the shock structures for the same Mach numbers for the hard-sphere model. Similarly, the normalized density 
$\bar{\unicode[STIX]{x1D70C}}$, velocity 
$\bar{v}$ and temperature 
$\bar{\unicode[STIX]{x1D703}}$ are plotted in figure 9. It is interesting that, when the Mach number increases from 
$3.0$ to 
$9.0$, there is no significant decrement of the general quality of the R13 approximation. Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2004) calculated the shock structure for the hard-sphere model using the R13 equations for Maxwell molecules where in these R13 equations the expression for the viscosity is changed to match the hard-sphere model. At Mach number 
$3.0$, such a method already shows significant deviation in the profile of the heat flux. After switching to the ‘true’ R13 equations for hard spheres, much better agreement can be obtained. Note that the peak of the heat flux is again underestimated in all results (including the model with 
$\unicode[STIX]{x1D702}=10$). In general, up to a Mach number of 
$9.0$, the R13 results show quite satisfactory agreement with the reference solutions for both models.
Figure 7. Normalized density, velocity and temperature of shock structures for the IPL model with 
$\unicode[STIX]{x1D702}=10$ and Mach numbers 
$Ma=1.55,3,6.5,9$. DSMC solutions for the variable hard-sphere model are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path.
Figure 8. The stress 
$\unicode[STIX]{x1D70E}$ and heat flux 
$q$ of shock structures for the IPL model with 
$\unicode[STIX]{x1D702}=10$ and Mach numbers 
$Ma=1.55,3,6.5,9$. DSMC results for the variable hard-sphere model are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path. The left 
$y$-axis corresponds to the stress and the right 
$y$-axis corresponds to the heat flux.
Figure 9. Normalized density, velocity and temperature of shock structures for the hard-sphere model and Mach numbers 
$Ma=1.55,3,6.5,9$. DSMC solutions for the variable hard-sphere model are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path.
Figure 10. The stress 
$\unicode[STIX]{x1D70E}$ and heat flux 
$q$ of shock structures for the hard-sphere model and Mach numbers 
$Ma=1.55,3,6.5,9$. DSMC results for the variable hard-sphere model are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path. The left 
$y$-axis corresponds to the stress and the right 
$y$-axis corresponds to the heat flux.
4.3 Shock structures for different indices $\unicode[STIX]{x1D702}$
Now we perform the tests by fixing the Mach number as 
$Ma=6.5$ and changing the parameter 
$\unicode[STIX]{x1D702}$. Here, we focus only on hard potentials with 
$\unicode[STIX]{x1D702}=7,10,17$ and 
$\infty$ (hard-sphere model). The results for all the five quantities are plotted in figures 11 and 12. Similarly, the normalized density 
$\bar{\unicode[STIX]{x1D70C}}$, velocity 
$\bar{v}$ and temperature 
$\bar{\unicode[STIX]{x1D703}}$ are plotted in figure 11. Both R13 results and DSMC results show that the shock structure differs for different collision models, and it can be observed that better agreement between two results can be achieved for larger 
$\unicode[STIX]{x1D702}$. The possible reason is that larger 
$\unicode[STIX]{x1D702}$ gives a smaller viscosity index, which brings the distribution function closer to the local Maxwellian.
Again, the most obvious deviation between R13 and DSMC results appears in the plots of heat fluxes in the low-density region. In general, the distribution function inside a shock wave is similar to the superposition of two Maxwellians: a narrow one coming from the front of the shock wave and a wide one from the back of the shock wave (Valentini & Schwartzentruber Reference Valentini and Schwartzentruber2009). In the low-density region, the portion of the wide Maxwellian is quite small. However, when evaluating high-order moments, the contribution of this small portion of wide Maxwellian becomes obvious due to its slow decay at infinity. For the 13-moment approximation, it can be expected that the contribution of the tail may be underestimated, since the decay rate of the distribution function in the Chapman–Enskog expansion is mainly set by the local temperature, which is significantly faster than the wide Maxwellian in the low-density region.
As a summary, we observe that R13 models predicts reasonable shock structures both qualitatively and quantitatively, although the derivation of the models does not involve any special consideration for this specific problem. For all the IPL models tested, the kink structure exists when 
$Ma$ is large. For a fixed Mach number, such a problem looks milder when 
$\unicode[STIX]{x1D702}$ is larger. This indicates that the range of Mach numbers for which the R13 equations are valid expands as 
$\unicode[STIX]{x1D702}$ increases. This can also be observed by comparing figure 5(b) with figure 9(d), from which one can observe that for Maxwell molecules, the temperature profile is highly inaccurate for Mach number 
$9.0$, while for hard-sphere molecules, the model seems to remain valid. Current numerical studies are insufficient to determine a specific range for all IPL models, and we leave this for future work. Nevertheless, these numerical results indicate the potential use of the R13 model not only for the low-Knudsen-number case, but also for high-speed rarefied gas flows.
Figure 11. Normalized density, velocity and temperature of shock structures for IPL models with 
$\unicode[STIX]{x1D702}=7,10,17,\infty$ (hard sphere) and Mach numbers 
$Ma=6.5$. DSMC solutions for the variable hard-sphere models are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path.
Figure 12. The stress 
$\unicode[STIX]{x1D70E}$ and heat flux 
$q$ of shock structures for IPL models with
$\unicode[STIX]{x1D702}=7,10,17,\infty$ (hard sphere) and Mach numbers 
$Ma=6.5$. DSMC results for the variable hard-sphere models are provided as reference results. The horizontal axis is 
$x/\unicode[STIX]{x1D706}$ with 
$\unicode[STIX]{x1D706}$ being the mean free path. The left 
$y$-axis corresponds to the stress and the right 
$y$-axis corresponds to the heat flux.
5 Conclusion
In this work, we have derived the R13-moment equations for all IPL models and the hard-sphere model. This work can be considered as a generalization of Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) to a much more general class of gas molecules. It also generalizes the methodology of Struchtrup (Reference Struchtrup2005a), which proposed the derivation of GG13 equations for general collision models, to one more order of accuracy. The derivation follows a systematic routine that can, in principle, be applied to all collision models. In the numerical experiment for shock structures, these new models show good agreement with the kinetic model in strong non-equilibrium regimes. To better understand how the R13 equations describe the distribution functions, one may plot the distribution function predicted by the R13 models inside the shock structure, which is to be considered in future work.
A significant drawback of these models is the high complexity of the collision terms, which are given in appendix C. This may cause difficulties in both understanding the models and designing the numerical methods. One possible way to simplify the equations is to linearize the regularization terms, as in Torrilhon (Reference Torrilhon2006) and Cai, Li & Wang (Reference Cai, Li and Wang2012), which still needs further justification. We are currently also working on the derivation of R13-moment equations for the Boltzmann equation with binary collision terms.
Acknowledgements
Z.C. is supported by the National University of Singapore Startup Fund under grant no. R-146-000-241-133. Y.W. is supported by the National Natural Scientific Foundation of China (grant no. 91630310 and U1930402).
Declaration of interests
The authors report no conflict of interest.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2020.251.
Appendix A. Introduction to the Boltzmann equation and the linearized IPL model
As introduced in the beginning of § 3, both the GG13 equations and R13 equations are derived from the kinetic equation, which governs the distribution function 
$f(\boldsymbol{x},\unicode[STIX]{x1D743},t)$. The relation between the distribution function and the moments has been demonstrated in (3.1) and (3.2). An equivalent but more straightforward way to write down the relationship is the following:
where 
$\mathfrak{m}$ is the mass of a single molecule. For monatomic gases, the governing equation for the distribution function is the Boltzmann equation, which reads
$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D743}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}f=C[f],\end{eqnarray}$$ where 
$C(f)$ is the collision term. Here, we only focus on the linearized collision term, which can be expressed as (Harris Reference Harris1971)
$$\begin{eqnarray}\displaystyle & & \displaystyle C[f](\boldsymbol{x},\unicode[STIX]{x1D743},t)=\int _{\mathbb{R}^{3}}\int _{\boldsymbol{n}\bot \boldsymbol{g}}B(\boldsymbol{g},\unicode[STIX]{x1D712}){\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743},t){\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743}_{1},t)\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\left[\frac{f(\boldsymbol{x},\unicode[STIX]{x1D743}_{1}^{\prime },t)}{{\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743}_{1}^{\prime },t)}+\frac{f(\boldsymbol{x},\unicode[STIX]{x1D743}^{\prime },t)}{{\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743}^{\prime },t)}-\frac{f(\boldsymbol{x},\unicode[STIX]{x1D743}_{1},t)}{{\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743}_{1},t)}-\frac{f(\boldsymbol{x},\unicode[STIX]{x1D743},t)}{{\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743},t)}\right]\,\text{d}\unicode[STIX]{x1D712}\,\text{d}\boldsymbol{n}\,\text{d}\unicode[STIX]{x1D743}_{1},\qquad\end{eqnarray}$$ where 
${\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743},t)$ is the local Maxwellian
$$\begin{eqnarray}{\mathcal{M}}(\boldsymbol{x},\unicode[STIX]{x1D743},t)=\frac{\unicode[STIX]{x1D70C}(\boldsymbol{x},t)}{\mathfrak{m}(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D703}(\boldsymbol{x},t))^{3/2}}\exp \left(-\frac{|\unicode[STIX]{x1D743}-\boldsymbol{v}(\boldsymbol{x},t)|^{2}}{2\unicode[STIX]{x1D703}(\boldsymbol{x},t)}\right),\end{eqnarray}$$ which satisfies 
$C[{\mathcal{M}}]=0$. In (A 3), 
$\boldsymbol{g}=\unicode[STIX]{x1D743}-\unicode[STIX]{x1D743}_{1}$ and 
$\boldsymbol{n}$ is a unit vector. The post-collisional velocities 
$\unicode[STIX]{x1D743}^{\prime }$ and 
$\unicode[STIX]{x1D743}_{1}^{\prime }$ are
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D743}^{\prime }=\cos ^{2}(\unicode[STIX]{x1D712}/2)\unicode[STIX]{x1D743}+\sin ^{2}(\unicode[STIX]{x1D712}/2)\unicode[STIX]{x1D743}_{1}-|\boldsymbol{g}|\cos (\unicode[STIX]{x1D712}/2)\sin (\unicode[STIX]{x1D712}/2)\boldsymbol{n},\\ \unicode[STIX]{x1D743}_{1}^{\prime }=\cos ^{2}(\unicode[STIX]{x1D712}/2)\unicode[STIX]{x1D743}_{1}+\sin ^{2}(\unicode[STIX]{x1D712}/2)\unicode[STIX]{x1D743}+|\boldsymbol{g}|\cos (\unicode[STIX]{x1D712}/2)\sin (\unicode[STIX]{x1D712}/2)\boldsymbol{n}.\end{array}\right\}\end{eqnarray}$$ Now, the only unexplained term in the collision term (A 3) is the collision kernel 
$B(|\boldsymbol{g}|,\unicode[STIX]{x1D712})$, which is a non-negative function determined by the force between gas molecules. For the IPL model, it has the form (Bird Reference Bird1994)
$$\begin{eqnarray}B(|\boldsymbol{g}|,\unicode[STIX]{x1D712})=\left(\frac{2\unicode[STIX]{x1D705}}{\mathfrak{m}}\right)^{2/(\unicode[STIX]{x1D702}-1)}|\boldsymbol{g}|^{(\unicode[STIX]{x1D702}-5)/(\unicode[STIX]{x1D702}-1)}W_{0}\left|\frac{\text{d}W_{0}}{\text{d}\unicode[STIX]{x1D712}}\right|,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D702}$ is the same as the parameter used throughout this paper. The dimensionless impact parameter 
$W_{0}$ is related to the angle 
$\unicode[STIX]{x1D712}$ by the following two equations:
$$\begin{eqnarray}\unicode[STIX]{x1D712}=\unicode[STIX]{x03C0}-2\int _{0}^{W_{1}}\left[1-W^{2}-\frac{2}{\unicode[STIX]{x1D702}-1}\left(\frac{W}{W_{0}}\right)^{\unicode[STIX]{x1D702}-1}\right]^{-1/2}\,\text{d}W,\quad 1-W_{1}^{2}-\frac{2}{\unicode[STIX]{x1D702}-1}\left(\frac{W_{1}}{W_{0}}\right)^{\unicode[STIX]{x1D702}-1}=0.\end{eqnarray}$$We also refer the readers to Struchtrup (Reference Struchtrup2005b) for more information on the Boltzmann equation and the collision models.
Appendix B. Basis functions and moment equations
In this appendix, we explain the basis functions used in the expansion (3.1) and the terms in the moment equations (3.4). Here, we follow Cai & Torrilhon (Reference Cai and Torrilhon2015) to define the basis function 
$\unicode[STIX]{x1D713}_{lmn}(\boldsymbol{x},\unicode[STIX]{x1D743},t)$ as
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D713}_{lmn}(\boldsymbol{x},\unicode[STIX]{x1D743},t)\nonumber\\ \displaystyle & & \displaystyle \quad =[\unicode[STIX]{x1D703}(\boldsymbol{x},t)]^{-(2n+l)/2}p_{lmn}\left(\frac{\unicode[STIX]{x1D743}-\boldsymbol{v}(\boldsymbol{x},t)}{\sqrt{\unicode[STIX]{x1D703}(\boldsymbol{x},t)}}\right)\cdot [2\unicode[STIX]{x03C0}\unicode[STIX]{x1D703}(\boldsymbol{x},t)]^{-3/2}\exp \left(-\frac{|\unicode[STIX]{x1D743}-\boldsymbol{v}(\boldsymbol{x},t)|^{2}}{2\unicode[STIX]{x1D703}(\boldsymbol{x},t)}\right),\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ where 
$p_{lmn}(\cdot )$ is an orthogonal polynomial in 
$\mathbb{R}^{3}$
$$\begin{eqnarray}p_{lmn}(\unicode[STIX]{x1D743})=\sqrt{\frac{2^{1-l}\unicode[STIX]{x03C0}^{3/2}n!}{\unicode[STIX]{x1D6E4}(n+l+3/2)}}L_{n}^{(l+1/2)}\left(\frac{|\unicode[STIX]{x1D743}|^{2}}{2}\right)|\unicode[STIX]{x1D743}|^{l}Y_{l}^{m}\left(\frac{\unicode[STIX]{x1D743}}{|\unicode[STIX]{x1D743}|}\right),\quad l,n\in \mathbb{N},\quad m=-l,\ldots ,l,\end{eqnarray}$$where we have used Laguerre polynomials
$$\begin{eqnarray}L_{n}^{(\unicode[STIX]{x1D6FC})}(x)=\frac{x^{-\unicode[STIX]{x1D6FC}}\exp (x)}{n!}\frac{\text{d}^{n}}{\text{d}x^{n}}[x^{n+\unicode[STIX]{x1D6FC}}\exp (-x)],\end{eqnarray}$$and spherical harmonics
$$\begin{eqnarray}Y_{l}^{m}(\boldsymbol{n})=\sqrt{\frac{2l+1}{4\unicode[STIX]{x03C0}}\frac{(l-m)!}{(l+m)!}}P_{l}^{m}(\cos \unicode[STIX]{x1D703})\exp (\text{i}m\unicode[STIX]{x1D719}),\quad \boldsymbol{n}=(\sin \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719},\sin \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719},\cos \unicode[STIX]{x1D703})^{\text{T}},\end{eqnarray}$$ with 
$P_{l}^{m}$ being the associate Legendre polynomial
$$\begin{eqnarray}P_{l}^{m}(x)=\frac{(-1)^{m}}{2^{l}l!}(1-x^{2})^{m/2}\frac{\text{d}^{l+m}}{\text{d}x^{l+m}}[(x^{2}-1)^{l}].\end{eqnarray}$$ The orthogonality of the polynomials 
$p_{lmn}$ is
$$\begin{eqnarray}\frac{1}{(2\unicode[STIX]{x03C0})^{3/2}}\int _{\mathbb{R}^{3}}\overline{p_{l_{1}m_{1}n_{1}}(\unicode[STIX]{x1D743})}p_{l_{2}m_{2}n_{2}}(\unicode[STIX]{x1D743})\exp \left(-\frac{|\unicode[STIX]{x1D743}|^{2}}{2}\right)\text{d}\unicode[STIX]{x1D743}=\unicode[STIX]{x1D6FF}_{l_{1}l_{2}}\unicode[STIX]{x1D6FF}_{m_{1}m_{2}}\unicode[STIX]{x1D6FF}_{n_{1}n_{2}}.\end{eqnarray}$$ Next, we show the expressions for 
$S_{lmn}$ and 
$T_{lmn}$ in (3.4), which have been obtained by Cai & Torrilhon (Reference Cai and Torrilhon2018). For simplicity, we introduce the following velocities:
$$\begin{eqnarray}V_{-1}={\textstyle \frac{1}{2}}(v_{1}-\text{i}v_{2}),\quad V_{0}=v_{3},\quad V_{1}=-{\textstyle \frac{1}{2}}(v_{1}+\text{i}v_{2}).\end{eqnarray}$$Then, we have
$$\begin{eqnarray}\displaystyle & & \displaystyle S_{lmn}=-\sqrt{n(n+l+1/2)}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}}f_{l,m,n-1}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\sqrt{2}\mathop{\sum }_{\unicode[STIX]{x1D707}=-1}^{1}{\displaystyle \frac{\unicode[STIX]{x2202}V_{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}t}}[(-1)^{\unicode[STIX]{x1D707}}\sqrt{n+l+1/2}\unicode[STIX]{x1D6FE}_{l,m+\unicode[STIX]{x1D707}}^{-\unicode[STIX]{x1D707}}\,f_{l-1,m+\unicode[STIX]{x1D707},n}-\sqrt{n}\unicode[STIX]{x1D6FE}_{-l-1,m+\unicode[STIX]{x1D707}}^{-\unicode[STIX]{x1D707}}\,f_{l+1,m+\unicode[STIX]{x1D707},n-1}],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FE}_{lm}^{\unicode[STIX]{x1D707}}$ are constants defined by
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{lm}^{\unicode[STIX]{x1D707}}=\sqrt{\frac{[l+(2\unicode[STIX]{x1D6FF}_{1,\unicode[STIX]{x1D707}}-1)m+\unicode[STIX]{x1D6FF}_{1,\unicode[STIX]{x1D707}}][l-(2\unicode[STIX]{x1D6FF}_{-1,\unicode[STIX]{x1D707}}-1)m+\unicode[STIX]{x1D6FF}_{-1,\unicode[STIX]{x1D707}}]}{(2l-1)(2l+1)}}.\end{eqnarray}$$ To introduce 
$T_{lmn}$, we first define the following operators:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}X_{-1}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{1}}+\text{i}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{2}},\quad \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}X_{0}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{3}},\quad \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}X_{1}}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{1}}+\text{i}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{2}},\end{eqnarray}$$ using which we can write 
$T_{lmn}$ as
$$\begin{eqnarray}\displaystyle T_{lmn} & = & \displaystyle \mathop{\sum }_{\unicode[STIX]{x1D707}=-1}^{1}\left(V_{\unicode[STIX]{x1D707}}F_{lmn\unicode[STIX]{x1D707}}+\frac{1}{2^{|\unicode[STIX]{x1D707}|}}[\!\sqrt{2(n+l)+1}\unicode[STIX]{x1D6FE}_{l,m-\unicode[STIX]{x1D707}}^{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D703}F_{l-1,m-\unicode[STIX]{x1D707},n,\unicode[STIX]{x1D707}}\right.\nonumber\\ \displaystyle & & \displaystyle -\,\sqrt{2(n+1)}\unicode[STIX]{x1D6FE}_{l,m-\unicode[STIX]{x1D707}}^{\unicode[STIX]{x1D707}}F_{l-1,m-\unicode[STIX]{x1D707},n+1,\unicode[STIX]{x1D707}}\nonumber\\ \displaystyle & & \displaystyle +\left.(-1)^{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D6FE}_{-l-1,m-\unicode[STIX]{x1D707}}^{\unicode[STIX]{x1D707}}(\sqrt{2(n+l)+3}F_{l+1,m-\unicode[STIX]{x1D707},n,\unicode[STIX]{x1D707}}-\sqrt{2n}\unicode[STIX]{x1D703}F_{l+1,m-\unicode[STIX]{x1D707},n-1,\unicode[STIX]{x1D707}})\!]\right),\qquad\end{eqnarray}$$where
$$\begin{eqnarray}\displaystyle F_{lmn\unicode[STIX]{x1D707}} & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}f_{lmn}}{\unicode[STIX]{x2202}X_{\unicode[STIX]{x1D707}}}}-\sqrt{n(n+l+1/2)}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}X_{\unicode[STIX]{x1D707}}}}\,f_{l,m,n-1}\nonumber\\ \displaystyle & & \displaystyle +\,\sqrt{2}\mathop{\sum }_{\unicode[STIX]{x1D708}=-1}^{1}{\displaystyle \frac{\unicode[STIX]{x2202}V_{\unicode[STIX]{x1D708}}}{\unicode[STIX]{x2202}X_{\unicode[STIX]{x1D707}}}}[(-1)^{\unicode[STIX]{x1D708}}\sqrt{n+l+1/2}\unicode[STIX]{x1D6FE}_{l,m+\unicode[STIX]{x1D708}}^{-\unicode[STIX]{x1D708}}\,f_{l-1,m+\unicode[STIX]{x1D708},n}-\sqrt{n}\unicode[STIX]{x1D6FE}_{-l-1,m+\unicode[STIX]{x1D708}}^{-\unicode[STIX]{x1D708}}\,f_{l+1,m+\unicode[STIX]{x1D708},n-1}].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$Table 6. Coefficients 
$\unicode[STIX]{x1D6EF}_{i}^{(\unicode[STIX]{x1D702})}$ in (C 1) for different 
$\unicode[STIX]{x1D702}$.
Table 7. Coefficients 
$\unicode[STIX]{x1D6EC}_{i}^{(\unicode[STIX]{x1D702})}$ in (C 2) for different 
$\unicode[STIX]{x1D702}$.
Appendix C. R13 collision terms
In this appendix, we provide the explicit forms of 
$\unicode[STIX]{x1D6F4}_{ij}^{(\unicode[STIX]{x1D702},2)}$ and 
$Q_{i}^{(\unicode[STIX]{x1D702},2)}$ and tabulate some values of the coefficients for some choices of 
$\unicode[STIX]{x1D702}$. In expressions (
) and (
C 2), 
$\unicode[STIX]{x1D6EF}_{k}^{(\unicode[STIX]{x1D702})}$ and 
$\unicode[STIX]{x1D6EC}_{k}^{(\unicode[STIX]{x1D702})}$ are constants depending only on 
$\unicode[STIX]{x1D702}$, whose values are given in tables 6 and 7. These tables are to be read horizontally. For example, in table 6, the row below ‘
$\unicode[STIX]{x1D702}=10$’ gives the values of 
$\unicode[STIX]{x1D6EF}_{0}^{(10)},\unicode[STIX]{x1D6EF}_{1}^{(10)},\ldots ,\unicode[STIX]{x1D6EF}_{9}^{(10)}$, and the next row gives the values of 
$\unicode[STIX]{x1D6EF}_{10}^{(10)},\unicode[STIX]{x1D6EF}_{11}^{(10)},\ldots ,\unicode[STIX]{x1D6EF}_{19}^{(10)}$.
