3.1 Counting moments in $d$ dimensions

Before deriving the various moment systems, it is worthwhile knowing how many moments a Grad moment system contains in a general dimension $d$. As it is more convenient to work with symmetric–traceless moments, the number of moments in a Grad moment system in a general dimension can be determined by knowing the number of independent components in a symmetric $r$-rank tensor and the number of traces in this tensor. Indeed, the number of independent components of a fully symmetric $r$-rank tensor in $d$ dimensions is $\binom{d+r-1}{r}=[(d+r-1)!]/[r!(d-1)!]$, and the number of traces in this tensor is $0$ for $r\in \{0,1\}$ while $\binom{d+r-3}{r-2}=[(d+r-3)!]/[(r-2)!(d-1)!]$ for $r\in \mathbb{N}\,\setminus \,\{1\}$. Consequently, the number of independent components of a fully symmetric–traceless $r$-rank tensor ($r\in \mathbb{N}\,\setminus \,\{1\}$) in $d$ dimensions is

$$\begin{eqnarray}\binom{d+r-1}{r}-\binom{d+r-3}{r-2}=\frac{d+2r-2}{d+r-2}\binom{d+r-2}{r}.\end{eqnarray}$$

Notably, any symmetric–traceless $r$-rank tensor ($r\in \mathbb{N}$) in two dimensions has only two independent components while any symmetric–traceless $r$-rank tensor ($r\in \mathbb{N}$) in three dimensions has $2r+1$ independent components. The counting of number of moments in some of the Grad moment systems considered in this paper is illustrated in table 1. Notwithstanding, any Grad moment system considered in this paper henceforth will be referred by its number of moments in three dimensions since Grad moment systems with the number of moments in three dimensions are more familiar to us (see e.g. Jenkins & Richman Reference Jenkins and Richman1985a; Levermore Reference Levermore1996; Struchtrup Reference Struchtrup2005; Kremer & Marques Jr. Reference Kremer and Marques2011). For instance, the Grad $(d^{2}+5d+2)/2!$-, $(d+1)(d^{2}+8d+6)/3!$- or $(d+3)(d^{2}+6d+2)/3!$-moment systems will simply be referred to as the Grad $13$-, $26$- or $29$-moment systems, respectively, which we are more acquainted with.

3.2 The system of the 29-moment equations

The system of the 29-moment equations includes the governing equations for the third rank tensor, for one and full traces of the fourth rank tensor and for full trace of the fifth rank tensor along with the governing equations for the well-known 13 moments. In other words, the system of the 29-moment equations consists of the governing equations for the moments $n$, $v_{i}$, $T$, $\unicode[STIX]{x1D70E}_{ij}$, $q_{i}$, $u_{ijk}^{0}$, $u^{2}$, $u_{ij}^{1}$, $u_{i}^{2}$, and is obtained by multiplying the Boltzmann equation (2.1) with $1$, $c_{i}$, $(1/d)mC^{2}$, $mC_{\langle \!i}C_{j\!\rangle }$, $\frac{1}{2}mC^{2}C_{i}$, $mC_{\langle \!i}C_{j}C_{k\!\rangle }$, $mC^{4}$, $mC^{2}C_{\langle \!i}C_{j\!\rangle }$ and $mC^{4}C_{i}$, and integrating each of the resulting equations over the velocity space successively. The detailed derivation of the 29-moment equations is provided as supplementary material available at https://doi.org/10.1017/jfm.2020.20. Here, they are presented directly. The system of the 29-moment equations consists of the mass, momentum and energy balance equations (2.8)–(2.10) and other higher-order moment equations, which on using the abbreviations

(3.3)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle m_{ijk}:=u_{ijk}^{0},\quad \unicode[STIX]{x1D6E5}:=\frac{u^{2}}{d(d+2)\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}-1,\\ \displaystyle R_{ij}:=u_{ij}^{1}-(d+4)\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij},\quad \unicode[STIX]{x1D711}_{i}:=u_{i}^{2}-4(d+4)\unicode[STIX]{x1D703}q_{i}\end{array}\right\}\end{eqnarray}$$

read

(3.4)$$\begin{eqnarray}\frac{\text{D}\unicode[STIX]{x1D70E}_{ij}}{\text{D}t}+\frac{\unicode[STIX]{x2202}m_{ijk}}{\unicode[STIX]{x2202}x_{k}}+\frac{4}{d+2}\frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}+2\unicode[STIX]{x1D70E}_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}+2\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}={\mathcal{P}}_{ij}^{0},\end{eqnarray}$$

(3.5)$$\begin{eqnarray}\displaystyle \frac{\text{D}q_{i}}{\text{D}t} & + & \displaystyle \frac{1}{2}\frac{\unicode[STIX]{x2202}R_{ij}}{\unicode[STIX]{x2202}x_{j}}+\frac{d+2}{2}\left[\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}+(1+2\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j}}\right]\nonumber\\ \displaystyle & + & \displaystyle \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}-\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j}}\right)+m_{ijk}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}\nonumber\\ \displaystyle & + & \displaystyle \frac{d+4}{d+2}q_{i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}+\frac{d+4}{d+2}q_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{2}{d+2}q_{j}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{i}}=\frac{1}{2}{\mathcal{P}}_{i}^{1},\end{eqnarray}$$

(3.6)$$\begin{eqnarray}\displaystyle \frac{\text{D}m_{ijk}}{\text{D}t} & + & \displaystyle \frac{\unicode[STIX]{x2202}u_{ijkl}^{0}}{\unicode[STIX]{x2202}x_{l}}+\frac{3}{d+4}\frac{\unicode[STIX]{x2202}R_{\langle \!ij}}{\unicode[STIX]{x2202}x_{k\!\rangle }}+3\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\langle \!ij}}{\unicode[STIX]{x2202}x_{k\!\rangle }}-3\frac{\unicode[STIX]{x1D70E}_{\langle \!ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{k\!\rangle l}}{\unicode[STIX]{x2202}x_{l}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k\!\rangle }}\right)\nonumber\\ \displaystyle & + & \displaystyle m_{ijk}\frac{\unicode[STIX]{x2202}v_{l}}{\unicode[STIX]{x2202}x_{l}}+3m_{l\langle \!ij}\frac{\unicode[STIX]{x2202}v_{k\!\rangle }}{\unicode[STIX]{x2202}x_{l}}+\frac{12}{d+2}q_{\langle \!i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k\!\rangle }}={\mathcal{P}}_{ijk}^{0},\end{eqnarray}$$

(3.7)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}\unicode[STIX]{x1D6E5}}{\text{D}t}+\frac{8}{d(d+2)}\frac{1}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}\left(1-\frac{d+2}{2}\unicode[STIX]{x1D6E5}\right)\left(\frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{d(d+2)}\frac{1}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{i}}{\unicode[STIX]{x2202}x_{i}}+4(d+2)q_{i}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}-8\frac{q_{i}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}\right)+4R_{ij}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{d(d+2)}\frac{1}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}[{\mathcal{P}}^{2}-2(d+2)(1+\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D703}{\mathcal{P}}^{1}],\end{eqnarray}$$

(3.8)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}R_{ij}}{\text{D}t}+\frac{2}{d+2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+\frac{4(d+4)}{d+2}\left(\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+q_{\langle \!i}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j\!\rangle }}-\frac{q_{\langle \!i}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{j\rangle k}}{\unicode[STIX]{x2202}x_{k}}-\frac{\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D70C}}q_{\langle \!i}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j\!\rangle }}\right)\nonumber\\ \displaystyle & & \displaystyle \phantom{\frac{\text{D}R_{ij}}{\text{D}t}}+\,4\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+4\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}-\frac{8}{d}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}-\frac{2(d+4)}{d}\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}q_{k}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D70E}_{kl}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{l}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x2202}u_{ijk}^{1}}{\unicode[STIX]{x2202}x_{k}}-(d+4)\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}m_{ijk}}{\unicode[STIX]{x2202}x_{k}}-2\frac{m_{ijk}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{kl}}{\unicode[STIX]{x2202}x_{l}}+\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k}}\right)+2u_{ijkl}^{0}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{l}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{d+6}{d+4}\left(R_{ij}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}+2R_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}\right)+\frac{4}{d+4}R_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+2(d+4)\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}\nonumber\\ \displaystyle & & \displaystyle \quad ={\mathcal{P}}_{ij}^{1}-(d+4)\unicode[STIX]{x1D703}{\mathcal{P}}_{ij}^{0}-\frac{d+4}{d}\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}{\mathcal{P}}^{1},\end{eqnarray}$$

(3.9)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}\unicode[STIX]{x1D711}_{i}}{\text{D}t}-\frac{8(d+4)}{d}\frac{q_{i}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}q_{j}}{\unicode[STIX]{x2202}x_{j}}+\unicode[STIX]{x1D70E}_{jk}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}\right)+\frac{\unicode[STIX]{x2202}u_{ij}^{2}}{\unicode[STIX]{x2202}x_{j}}+\frac{1}{d}\frac{\unicode[STIX]{x2202}u^{3}}{\unicode[STIX]{x2202}x_{i}}-2(d+4)\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}R_{ij}}{\unicode[STIX]{x2202}x_{j}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,4R_{ij}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j}}-(d+4)[(d+6)+(d+2)\unicode[STIX]{x1D6E5}]\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}-2(d+4)^{2}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,(d+2)(d+4)\left[2\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x2202}x_{i}}+(1+3\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D703}^{3}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}+(3+5\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad -\,4\frac{R_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j}}\right)+4u_{ijk}^{1}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}-4(d+4)\unicode[STIX]{x1D703}m_{ijk}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{8(d+4)}{d+2}\unicode[STIX]{x1D703}\left(q_{i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}+q_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}+q_{j}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{i}}\right)+\frac{d+6}{d+2}\unicode[STIX]{x1D711}_{i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}+\frac{d+6}{d+2}\unicode[STIX]{x1D711}_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{4}{d+2}\unicode[STIX]{x1D711}_{j}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{i}}={\mathcal{P}}_{i}^{2}-2(d+4)\unicode[STIX]{x1D703}{\mathcal{P}}_{i}^{1}-\frac{4(d+4)}{d}\frac{q_{i}}{\unicode[STIX]{x1D70C}}{\mathcal{P}}^{1}.\end{eqnarray}$$

The abbreviations (3.3) are introduced in such a way that $m_{ijk}$, $\unicode[STIX]{x1D6E5}$, $R_{ij}$ and $\unicode[STIX]{x1D711}_{i}$ vanish if computed with the well-known G13 distribution function

(3.10)$$\begin{eqnarray}f_{|\text{G13}}=f_{M}\left[1+\frac{1}{2}\frac{\unicode[STIX]{x1D70E}_{ij}C_{i}C_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}+\frac{q_{i}C_{i}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}\left(\frac{1}{d+2}\frac{C^{2}}{\unicode[STIX]{x1D703}}-1\right)\right],\end{eqnarray}$$

where

(3.11)$$\begin{eqnarray}f_{M}\equiv f_{M}(t,\boldsymbol{x},\boldsymbol{c})=n\left(\frac{1}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D703}}\right)^{d/2}\exp \left(-\frac{C^{2}}{2\unicode[STIX]{x1D703}}\right)\end{eqnarray}$$

is the Maxwellian distribution function (Garzó Reference Garzó2013). In general, the computation of the collisional production terms ${\mathcal{P}}_{i_{1}i_{2}\ldots i_{r}}^{a}$ requires the knowledge of the distribution function and is not easy for particles interacting with a general interaction potential. Nevertheless, for IMM (considered in this work), the collisional production terms can be evaluated easily – indeed, without the knowledge of the explicit form of the distribution function. A strategy for computing them for IMM in an automated way using the computer algebra software Mathematica^® is demonstrated in appendix B. Using this strategy, the production terms associated with the G29 equations for $d$-dimensional IMM have been computed. They turn out to be

(3.12)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}^{1}=-\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D708}d\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$

(3.13)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}_{ij}^{0}=-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D708}\unicode[STIX]{x1D70E}_{ij}, & \displaystyle\end{eqnarray}$$

(3.14)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}_{i}^{1}=-2\unicode[STIX]{x1D708}_{q}^{\ast }\unicode[STIX]{x1D708}q_{i}, & \displaystyle\end{eqnarray}$$

(3.15)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}_{ijk}^{0}=-\unicode[STIX]{x1D708}_{m}^{\ast }\unicode[STIX]{x1D708}m_{ijk}, & \displaystyle\end{eqnarray}$$

(3.16)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}^{2}=-\unicode[STIX]{x1D708}\left[(\unicode[STIX]{x1D6FC}_{0}+\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}+d(d+2)\unicode[STIX]{x1D70D}_{0}\frac{\unicode[STIX]{x1D70E}_{ij}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}\right], & \displaystyle\end{eqnarray}$$

(3.17)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}_{ij}^{1}=-\unicode[STIX]{x1D708}\left[\unicode[STIX]{x1D708}_{R}^{\ast }R_{ij}+\unicode[STIX]{x1D6FC}_{2}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}+\unicode[STIX]{x1D70D}_{1}\frac{\unicode[STIX]{x1D70E}_{k\langle \!i}\unicode[STIX]{x1D70E}_{j\!\rangle k}}{\unicode[STIX]{x1D70C}}\right], & \displaystyle\end{eqnarray}$$

(3.18)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{P}}_{i}^{2}=-\unicode[STIX]{x1D708}\left[\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }\unicode[STIX]{x1D711}_{i}+\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D703}q_{i}+\unicode[STIX]{x1D70D}_{2}\frac{\unicode[STIX]{x1D70E}_{ij}q_{j}}{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D70D}_{3}\frac{m_{ijk}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x1D70C}}\right], & \displaystyle\end{eqnarray}$$

where the coefficients $\unicode[STIX]{x1D701}_{0}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }$, $\unicode[STIX]{x1D708}_{q}^{\ast }$, $\unicode[STIX]{x1D708}_{m}^{\ast }$, $\unicode[STIX]{x1D708}_{R}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }$, $\unicode[STIX]{x1D6FC}_{0}$, $\unicode[STIX]{x1D6FC}_{1}$, $\unicode[STIX]{x1D6FC}_{2}$, $\unicode[STIX]{x1D6FC}_{3}$, $\unicode[STIX]{x1D70D}_{0}$, $\unicode[STIX]{x1D70D}_{1}$, $\unicode[STIX]{x1D70D}_{2}$ and $\unicode[STIX]{x1D70D}_{3}$ depend only on the dimension $d$ and coefficient of restitution $e$, and are relegated to appendix C for better readability. Collisional production terms (3.12)–(3.17) for IMM agree with those obtained in Garzó & Santos (Reference Garzó and Santos2007), wherein they have been computed till fourth order. Moreover, the coefficients $\unicode[STIX]{x1D701}_{0}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }$, $\unicode[STIX]{x1D708}_{q}^{\ast }$, $\unicode[STIX]{x1D708}_{R}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }$ and $\unicode[STIX]{x1D6FC}_{1}$ relate to the collisional rate $\unicode[STIX]{x1D708}_{2r|s}$ – associated with the Ikenberry polynomial $Y_{2r|i_{1}i_{2}\ldots i_{s}}(\boldsymbol{C})$ – given in Santos & Garzó (Reference Santos and Garzó2012) for $s\in \{0,1,2\}$ via $\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{2|0}$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{0|2}$, $\unicode[STIX]{x1D708}_{q}^{\ast }\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{2|1}$, $\unicode[STIX]{x1D708}_{R}^{\ast }\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{2|2}$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{4|1}$ and $\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D708}=d(d+2)\unicode[STIX]{x1D708}_{4|0}$. I could not find the full expression for the collisional production term (3.18) for granular gases in the existing literature. Nonetheless, for monatomic gases (i.e. for $d=3$ and $e=1$), it can be found, for instance, in Gu & Emerson (Reference Gu and Emerson2009) – although not explicitly. The source code for computing the above collisional production terms is provided as supplementary material with the present paper. The collisional production terms associated with the G26 equations for three-dimensional IHS can be found in Gupta & Torrilhon (Reference Gupta and Torrilhon2012) and Gupta et al. (Reference Gupta, Shukla and Torrilhon2018).

The relation ${\mathcal{P}}^{1}=-dnT\unicode[STIX]{x1D701}$, on exploiting (3.12), gives the cooling rate for IMM,

(3.19)$$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D708},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}_{0}^{\ast }=(d+2)(1-e^{2})/(4d)$ (see (C 1)). The cooling rate (3.19) is the same as that obtained in Santos (Reference Santos2003), Garzó & Santos (Reference Garzó and Santos2011) and Khalil et al. (Reference Khalil, Garzó and Santos2014), and vanishes identically for monatomic gases (i.e. for $e=1$), guaranteeing the conservation of energy for them. It is important to note from (3.19) that the cooling rate for IMM neither depends on the gradients of any field nor on any higher-order moment (in contrast to the cooling rate for IHS that also depends on the scalar fourth moment $\unicode[STIX]{x1D6E5}$; see Gupta et al. (Reference Gupta, Shukla and Torrilhon2018)).

3.3 Grad 29-moment closure

The system of the 29-moment equations for IMM ((2.8)–(2.10) and (3.4)–(3.9) along with collisional production terms (3.12)–(3.18)) is still not closed as it possesses the additional unknown moments $u_{ijkl}^{0}$, $u_{ijk}^{1}$, $u_{ij}^{2}$, $u^{3}$. The system is closed with the Grad distribution function based on the considered 29 moments, which is referred to as the G29 distribution function. The ($d$-dimensional) G29 distribution function $f_{|\text{G29}}$ reads

(3.20)$$\begin{eqnarray}\displaystyle f_{|\text{G29}} & = & \displaystyle f_{M}\left[1+\frac{1}{2}\frac{\unicode[STIX]{x1D70E}_{ij}C_{i}C_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}+\frac{q_{i}C_{i}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}\left(\frac{1}{d+2}\frac{C^{2}}{\unicode[STIX]{x1D703}}-1\right)+\frac{1}{6}\frac{m_{ijk}C_{i}C_{j}C_{k}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\frac{d(d+2)\unicode[STIX]{x1D6E5}}{8}\left(1-\frac{2}{d}\frac{C^{2}}{\unicode[STIX]{x1D703}}+\frac{1}{d(d+2)}\frac{C^{4}}{\unicode[STIX]{x1D703}^{2}}\right)+\frac{1}{4}\frac{R_{ij}C_{i}C_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}}\left(\frac{1}{d+4}\frac{C^{2}}{\unicode[STIX]{x1D703}}-1\right)\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{1}{8}\frac{\unicode[STIX]{x1D711}_{i}C_{i}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}}\left(1-\frac{2}{d+2}\frac{C^{2}}{\unicode[STIX]{x1D703}}+\frac{1}{(d+2)(d+4)}\frac{C^{4}}{\unicode[STIX]{x1D703}^{2}}\right)\right].\end{eqnarray}$$

The details of computing the G29 distribution function (3.20) can be found in appendix D. Insertion of the G29 distribution function (3.20) into the definitions of unknown moments $u_{ijkl}^{0}$, $u_{ijk}^{1}$, $u_{ij}^{2}$ and $u^{3}$ expresses them in terms of the considered 29 moments:

(3.21a)$$\begin{eqnarray}\displaystyle & \displaystyle u_{ijkl|\text{G29}}^{0}=0, & \displaystyle\end{eqnarray}$$

(3.21b)$$\begin{eqnarray}\displaystyle & \displaystyle u_{ijk|\text{G29}}^{1}=(d+6)\unicode[STIX]{x1D703}m_{ijk}, & \displaystyle\end{eqnarray}$$

(3.21c)$$\begin{eqnarray}\displaystyle & \displaystyle u_{ij|\text{G29}}^{2}=(d+6)\unicode[STIX]{x1D703}[2R_{ij}+(d+4)\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}], & \displaystyle\end{eqnarray}$$

(3.21d)$$\begin{eqnarray}\displaystyle & \displaystyle u_{|\text{G29}}^{3}=d(d+2)(d+4)(1+3\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}, & \displaystyle\end{eqnarray}$$

$|\text{G29}$

3.4 The G29 system for IMM

Equations (2.8)–(2.10) and (3.4)–(3.9) closed with (3.21) and (3.12)–(3.18) form the system of the G29 equations for $d$-dimensional IMM. Combining all of them, the system of the G29 equations for $d$-dimensional IMM reads

(3.22)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}n}{\text{D}t}+n\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{i}}=0, & \displaystyle\end{eqnarray}$$

(3.23)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}v_{i}}{\text{D}t}+\frac{1}{mn}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}(nT)}{\unicode[STIX]{x2202}x_{i}}\right]-F_{i}=0, & \displaystyle\end{eqnarray}$$

(3.24)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}T}{\text{D}t}+\frac{2}{d}\frac{1}{n}\left[\frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}+nT\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{i}}\right]=-\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D708}T, & \displaystyle\end{eqnarray}$$

(3.25)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\unicode[STIX]{x1D70E}_{ij}}{\text{D}t}+\frac{\unicode[STIX]{x2202}m_{ijk}}{\unicode[STIX]{x2202}x_{k}}+\frac{4}{d+2}\frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}+2\unicode[STIX]{x1D70E}_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}+2\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}=-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D708}\unicode[STIX]{x1D70E}_{ij}, & \displaystyle\end{eqnarray}$$

(3.26)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}q_{i}}{\text{D}t}+\frac{1}{2}\frac{\unicode[STIX]{x2202}R_{ij}}{\unicode[STIX]{x2202}x_{j}}+\frac{d+2}{2}\left[\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x0394}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}+(1+2\unicode[STIX]{x1D6E5})\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j}}\right]\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}-\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x2202}x_{k}}-\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j}}\right)+m_{ijk}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{d+4}{d+2}q_{i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}+\frac{d+4}{d+2}q_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{2}{d+2}q_{j}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{i}}=-\unicode[STIX]{x1D708}_{q}^{\ast }\unicode[STIX]{x1D708}q_{i},\end{eqnarray}$$

(3.27)$$\begin{eqnarray}\displaystyle \frac{\text{D}m_{ijk}}{\text{D}t} & + & \displaystyle \frac{3}{d+4}\frac{\unicode[STIX]{x2202}R_{\langle \!ij}}{\unicode[STIX]{x2202}x_{k\!\rangle }}+3\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\langle \!ij}}{\unicode[STIX]{x2202}x_{k\!\rangle }}-3\frac{\unicode[STIX]{x1D70E}_{\langle \!ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{k\!\rangle l}}{\unicode[STIX]{x2202}x_{l}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k\!\rangle }}\right)\nonumber\\ \displaystyle & + & \displaystyle m_{ijk}\frac{\unicode[STIX]{x2202}v_{l}}{\unicode[STIX]{x2202}x_{l}}+3m_{l\langle \!ij}\frac{\unicode[STIX]{x2202}v_{k\!\rangle }}{\unicode[STIX]{x2202}x_{l}}+\frac{12}{d+2}q_{\langle \!i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k\!\rangle }}=-\unicode[STIX]{x1D708}_{m}^{\ast }\unicode[STIX]{x1D708}m_{ijk},\end{eqnarray}$$

(3.28)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}\unicode[STIX]{x1D6E5}}{\text{D}t}+\frac{8}{d(d+2)}\frac{1}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}}\left(1-\frac{d+2}{2}\unicode[STIX]{x1D6E5}\right)\left(\frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}x_{i}}+\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{d(d+2)}\frac{1}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{i}}{\unicode[STIX]{x2202}x_{i}}+4(d+2)q_{i}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}-8\frac{q_{i}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}\right)+4R_{ij}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D708}\left[\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6E5}}^{\ast }\left\{\unicode[STIX]{x1D6E5}-\frac{6(1-e)^{2}}{4d-7+6e-3e^{2}}\right\}+\text{}\underline{\unicode[STIX]{x1D70D}_{0}\frac{\unicode[STIX]{x1D70E}_{ij}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x1D703}^{2}}}\right],\end{eqnarray}$$

(3.29)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}R_{ij}}{\text{D}t}+\frac{2}{d+2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D711}_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+\frac{4(d+4)}{d+2}\left(\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}q_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+q_{\langle \!i}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j\!\rangle }}-\frac{q_{\langle \!i}}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{j\!\rangle k}}{\unicode[STIX]{x2202}x_{k}}-\frac{\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D70C}}q_{\langle \!i}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j\!\rangle }}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,4\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+4\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}-\frac{8}{d}\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}-\frac{2(d+4)}{d}\frac{\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}q_{k}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D70E}_{kl}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{l}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,2\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}m_{ijk}}{\unicode[STIX]{x2202}x_{k}}+(d+4)m_{ijk}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{k}}-2\frac{m_{ijk}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{kl}}{\unicode[STIX]{x2202}x_{l}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{k}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{d+6}{d+4}\left(R_{ij}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{k}}+2R_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{j\!\rangle }}{\unicode[STIX]{x2202}x_{k}}\right)+\frac{4}{d+4}R_{k\langle \!i}\frac{\unicode[STIX]{x2202}v_{k}}{\unicode[STIX]{x2202}x_{j\!\rangle }}+2(d+4)\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D708}\left(\unicode[STIX]{x1D708}_{R}^{\ast }R_{ij}-\unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}+\text{}\underline{\unicode[STIX]{x1D70D}_{1}\frac{\unicode[STIX]{x1D70E}_{k\langle \!i}\unicode[STIX]{x1D70E}_{j\!\rangle k}}{\unicode[STIX]{x1D70C}}}\right),\end{eqnarray}$$

(3.30)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{D}\unicode[STIX]{x1D711}_{i}}{\text{D}t}-\frac{8(d+4)}{d}\frac{q_{i}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}q_{j}}{\unicode[STIX]{x2202}x_{j}}+\unicode[STIX]{x1D70E}_{jk}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}\right)+4\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}R_{ij}}{\unicode[STIX]{x2202}x_{j}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,(d+2)(d+4)\unicode[STIX]{x1D703}^{2}\left[\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E5}}{\unicode[STIX]{x2202}x_{i}}+4\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}-\unicode[STIX]{x1D6E5}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}x_{j}}\right]\nonumber\\ \displaystyle & & \displaystyle \quad -\,4\frac{R_{ij}}{\unicode[STIX]{x1D70C}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x2202}x_{k}}+\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j}}\right)+2(d+4)R_{ij}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j}}+4(d+4)\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{j}}+8\unicode[STIX]{x1D703}m_{ijk}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{k}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{8(d+4)}{d+2}\unicode[STIX]{x1D703}\left(q_{i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}+q_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}+q_{j}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{i}}\right)+\frac{d+6}{d+2}\unicode[STIX]{x1D711}_{i}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{j}}+\frac{d+6}{d+2}\unicode[STIX]{x1D711}_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{4}{d+2}\unicode[STIX]{x1D711}_{j}\frac{\unicode[STIX]{x2202}v_{j}}{\unicode[STIX]{x2202}x_{i}}\nonumber\\ \displaystyle & & \displaystyle =-\unicode[STIX]{x1D708}\left(\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }\unicode[STIX]{x1D711}_{i}-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }\unicode[STIX]{x1D703}q_{i}+\text{}\underline{\unicode[STIX]{x1D70D}_{2}\frac{\unicode[STIX]{x1D70E}_{ij}q_{j}}{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D70D}_{3}\frac{m_{ijk}\unicode[STIX]{x1D70E}_{jk}}{\unicode[STIX]{x1D70C}}}\right),\end{eqnarray}$$

where

(3.31)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6E5}}^{\ast }=\frac{(1+e)^{2}(4d-7+6e-3e^{2})}{16d},\\ \displaystyle \unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }=\frac{3(1+e)^{2}(1-e)(d+2-2e)}{4d},\\ \displaystyle \unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }=\frac{3(1+e)^{2}(1-e)[5(d+2)-(d+14)e]}{4d},\end{array}\right\}\end{eqnarray}$$

and the other coefficients $\unicode[STIX]{x1D701}_{0}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }$, $\unicode[STIX]{x1D708}_{q}^{\ast }$, $\unicode[STIX]{x1D708}_{m}^{\ast }$, $\unicode[STIX]{x1D708}_{R}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }$, $\unicode[STIX]{x1D70D}_{0}$, $\unicode[STIX]{x1D70D}_{1}$, $\unicode[STIX]{x1D70D}_{2}$ and $\unicode[STIX]{x1D70D}_{3}$ appearing on the right-hand sides of the G29 equations (3.22)–(3.30) depend only on the dimension $d$ and coefficient of restitution $e$ (see appendix C for their expressions). For $d=3$ and $e=1$, these coefficients become $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }=1$, $\unicode[STIX]{x1D708}_{q}^{\ast }=2/3$, $\unicode[STIX]{x1D708}_{m}^{\ast }=3/2$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6E5}}^{\ast }=2/3$, $\unicode[STIX]{x1D708}_{R}^{\ast }=7/6$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }=1$, $\unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }=0$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }=0$, $\unicode[STIX]{x1D70D}_{0}=2/45$, $\unicode[STIX]{x1D70D}_{1}=2/3$, $\unicode[STIX]{x1D70D}_{2}=28/15$ and $\unicode[STIX]{x1D70D}_{3}=2/3$, which are the same as the respective coefficients for monatomic gases of Maxwell molecules; see e.g. Struchtrup (Reference Struchtrup2005) and Gu & Emerson (Reference Gu and Emerson2009). In particular, the vanishing coefficients $\unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }=0$ and $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }=0$ make the right-hand sides of the linearised G29 equations for monatomic gases of Maxwell molecules completely decoupled, which is not the case for granular gases. Furthermore, the underlined nonlinear terms in (3.28)–(3.30) will be discarded for simplicity while investigating the HCS of a granular gas in § 5.

3.5 Various Grad moment systems

The abbreviations (3.3) have been introduced in such a way that the smaller systems of the Grad moment equations can be obtained directly from the G29 system (3.22)–(3.30). The other Grad moment systems considered in this paper are as follows:

1. (i) The G13 system: the system of the 13-moment equations contains the governing equations for variables $n$, $v_{i}$, $T$, $\unicode[STIX]{x1D70E}_{ij}$ and $q_{i}$, i.e. it consists of (3.22)–(3.26). However, equations (3.22)–(3.26) contain additional unknowns $m_{ijk}$, $\unicode[STIX]{x1D6E5}$ and $R_{ij}$ that vanish on being computed with the G13 distribution function (3.10). Thus, the G13 system for $d$-dimensional IMM consists of (3.22)–(3.26) with $m_{ijk}=\unicode[STIX]{x1D6E5}=R_{ij}=0$.

2. (ii) The G14 system: the system of the 14-moment equations contains the governing equations for variables $n$, $v_{i}$, $T$, $\unicode[STIX]{x1D70E}_{ij}$, $q_{i}$ and $\unicode[STIX]{x1D6E5}$, i.e. it consists of (3.22)–(3.26) and (3.28). However, equations (3.22)–(3.26) and (3.28) contain additional unknowns $m_{ijk}$, $R_{ij}$ and $\unicode[STIX]{x1D711}_{i}$ that also vanish on being computed with the G14 distribution function

   (3.32)$$\begin{eqnarray}\displaystyle f_{|\text{G14}} & = & \displaystyle f_{M}\left[1+\frac{1}{2}\frac{\unicode[STIX]{x1D70E}_{ij}C_{i}C_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}+\frac{q_{i}C_{i}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}\left(\frac{1}{d+2}\frac{C^{2}}{\unicode[STIX]{x1D703}}-1\right)\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\frac{d(d+2)\unicode[STIX]{x1D6E5}}{8}\left(1-\frac{2}{d}\frac{C^{2}}{\unicode[STIX]{x1D703}}+\frac{1}{d(d+2)}\frac{C^{4}}{\unicode[STIX]{x1D703}^{2}}\right)\right],\end{eqnarray}$$
   which can be obtained easily by following a similar procedure presented in appendix D. Thus, the G14 system for $d$-dimensional IMM consists of (3.22)–(3.26) and (3.28) with $m_{ijk}=R_{ij}=\unicode[STIX]{x1D711}_{i}=0$.

3. (iii) The G26 system: the system of the 26-moment equations contains the governing equations for variables $n$, $v_{i}$, $T$, $\unicode[STIX]{x1D70E}_{ij}$, $q_{i}$, $m_{ijk}$, $\unicode[STIX]{x1D6E5}$ and $R_{ij}$, i.e. it consists of (3.22)–(3.26) and (3.6)–(3.8) with the right-hand sides computed using the collisional production terms (3.12)–(3.18). However, equations (3.6)–(3.8) contain additional unknowns $u_{ijkl}^{0}$, $\unicode[STIX]{x1D711}_{i}$ and $u_{ijk}^{1}$ that are computed with the G26 distribution function

   (3.33)$$\begin{eqnarray}\displaystyle & & \displaystyle f_{|\text{G26}}=f_{M}\left[1+\frac{1}{2}\frac{\unicode[STIX]{x1D70E}_{ij}C_{i}C_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}+\frac{q_{i}C_{i}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{2}}\left(\frac{1}{d+2}\frac{C^{2}}{\unicode[STIX]{x1D703}}-1\right)+\frac{1}{6}\frac{m_{ijk}C_{i}C_{j}C_{k}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}}\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{d(d+2)\unicode[STIX]{x1D6E5}}{8}\left(1-\frac{2}{d}\frac{C^{2}}{\unicode[STIX]{x1D703}}+\frac{1}{d(d+2)}\frac{C^{4}}{\unicode[STIX]{x1D703}^{2}}\right)+\frac{1}{4}\frac{R_{ij}C_{i}C_{j}}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}^{3}}\left(\frac{1}{d+4}\frac{C^{2}}{\unicode[STIX]{x1D703}}-1\right)\right],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
   which can also be obtained easily by following a similar procedure presented in appendix D. With the G26 distribution function (3.33), $u_{ijkl}^{0}$ and $\unicode[STIX]{x1D711}_{i}$ vanish and $u_{ijk}^{1}$ turns out to be $u_{ijk|\text{G26}}^{1}=(d+6)\unicode[STIX]{x1D703}m_{ijk}$, which is exactly the same as the value of $u_{ijk}^{1}$ obtained with the G29 distribution function (3.20). Therefore, inserting the G26 closure (i.e. $u_{ijkl}^{0}=\unicode[STIX]{x1D711}_{i}=0$ and $u_{ijk}^{1}=(d+6)\unicode[STIX]{x1D703}m_{ijk}$), equations (3.6)–(3.8) turn to (3.27)–(3.29) in which $\unicode[STIX]{x1D711}_{i}=0$. Thus, the G26 system for $d$-dimensional IMM consists of (3.22)–(3.29) with $\unicode[STIX]{x1D711}_{i}=0$.

It is worthwhile noting that the G13 and G26 theories belong to the category of ordered moment theories, which always include the neglected fluxes of a moment theory at the previous level (Torrilhon Reference Torrilhon2015). Also, there are other moment theories, which consider complete (traces and traceless) moments of a given order; such moment theories are referred to as full moment theories (Torrilhon Reference Torrilhon2015). The first few examples of full moment theories are the Grad 10-, 20- and 35-moment theories (in three dimensions). In this sense, the G14 and G29 theories considered in the present work neither belong to the category of ordered moment theories nor to that of full moment theories.

4.1 Zeroth-order contributions in spatial gradients

To zeroth order in the spatial gradients, the mass, momentum and energy balance equations (3.22)–(3.24) reduce to

(4.2a,b)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}=0,\quad \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}t}=0\quad \text{and}\quad \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D708}T,\end{eqnarray}$$

while the balance equations for the other higher moments (3.25)–(3.30) reduce to

(4.3)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70E}_{ij}=0,\quad q_{i}=0,\quad m_{ijk}=0,\quad \unicode[STIX]{x1D6E5}-\frac{6(1-e)^{2}}{4d-7+6e-3e^{2}}=0,\\ \displaystyle \unicode[STIX]{x1D708}_{R}^{\ast }R_{ij}-\unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D703}\unicode[STIX]{x1D70E}_{ij}=0,\quad \unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }\unicode[STIX]{x1D711}_{i}-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }\unicode[STIX]{x1D703}q_{i}=0.\end{array}\right\}\end{eqnarray}$$

Equations in (4.3) readily imply that

(4.4a,b)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{ij}=q_{i}=m_{ijk}=R_{ij}=\unicode[STIX]{x1D711}_{i}=0\quad \text{and}\quad \unicode[STIX]{x1D6E5}=a_{2},\end{eqnarray}$$

where

(4.5)$$\begin{eqnarray}a_{2}=\frac{6(1-e)^{2}}{4d-7+6e-3e^{2}}\end{eqnarray}$$

is the same as the value of the fourth cumulant for IMM reported in previous studies (Santos Reference Santos2003; Khalil et al. Reference Khalil, Garzó and Santos2014). Thus, to zeroth order in spatial gradients, $\unicode[STIX]{x1D70E}_{ij}$, $q_{i}$, $m_{ijk}$, $R_{ij}$ and $\unicode[STIX]{x1D711}_{i}$ are zero while $\unicode[STIX]{x1D6E5}=a_{2}$.

4.2 First-order contributions in spatial gradients

Now, the terms having first-order spatial derivatives are also retained in the moment equations. To first order in spatial gradients, moment equations (3.22)–(3.26), read

(4.6)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}t}=-v_{i}\frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}x_{i}}-n\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{i}}, & \displaystyle\end{eqnarray}$$

(4.7)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}t}=-v_{j}\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{j}}-\frac{1}{mn}\frac{\unicode[STIX]{x2202}(nT)}{\unicode[STIX]{x2202}x_{i}}, & \displaystyle\end{eqnarray}$$

(4.8)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}=-v_{i}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{i}}-\frac{2}{d}T\frac{\unicode[STIX]{x2202}v_{i}}{\unicode[STIX]{x2202}x_{i}}-\unicode[STIX]{x1D701}T, & \displaystyle\end{eqnarray}$$

(4.9)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}t}=-2\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D708}\unicode[STIX]{x1D70E}_{ij}, & \displaystyle\end{eqnarray}$$

(4.10)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}t}=-\frac{d+2}{2}\left[a_{2}\unicode[STIX]{x1D703}^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{i}}+(1+2a_{2})\unicode[STIX]{x1D70C}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x_{i}}\right]-\unicode[STIX]{x1D708}_{q}^{\ast }\unicode[STIX]{x1D708}q_{i}. & \displaystyle\end{eqnarray}$$

Notice that, unlike IHS (see Gupta et al. Reference Gupta, Shukla and Torrilhon2018), here, none of the balance equations (3.27)–(3.30) are required to determine the transport coefficients for IMM up to first-order accuracy in spatial gradients, since the stress and heat flux balance equations ((4.9) and (4.10)) have no coupling with the higher moments. The balance equations (3.27)–(3.30) will only be needed for computing the transport coefficients beyond the first-order accuracy in spatial gradients, which is not the focus of the present work.

The time derivatives of the stress and heat flux in (4.9) and (4.10) are computed using dimensional analysis and using the zeroth-order accurate mass, momentum and energy balance equations (4.2). They turn out to be (Garzó Reference Garzó2013; Gupta et al. Reference Gupta, Shukla and Torrilhon2018)

(4.11a,b)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ij}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D702}\unicode[STIX]{x1D701}\frac{\unicode[STIX]{x2202}v_{\langle \!i}}{\unicode[STIX]{x2202}x_{j\!\rangle }}\quad \text{and}\quad \frac{\unicode[STIX]{x2202}q_{i}}{\unicode[STIX]{x2202}t}=2\unicode[STIX]{x1D705}\unicode[STIX]{x1D701}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{i}}+\left(\unicode[STIX]{x1D705}\frac{T}{n}+\frac{3}{2}\unicode[STIX]{x1D706}\right)\unicode[STIX]{x1D701}\frac{\unicode[STIX]{x2202}n}{\unicode[STIX]{x2202}x_{i}}.\end{eqnarray}$$

Now, in the first-order accurate stress and heat flux balance equations ((4.9) and (4.10)), one replaces $\unicode[STIX]{x1D70E}_{ij}$ and $q_{i}$ using (4.1) and their time derivatives using (4.11). Subsequent comparison of the coefficients of each hydrodynamic gradient in both the resulting equations leads to the transport coefficients in the NSF laws (4.1):

(4.12a-c)$$\begin{eqnarray}\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{0}\unicode[STIX]{x1D702}^{\ast },\quad \unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}\unicode[STIX]{x1D705}^{\ast }\quad \text{and}\quad \unicode[STIX]{x1D706}=\frac{\unicode[STIX]{x1D705}_{0}T}{n}\unicode[STIX]{x1D706}^{\ast },\end{eqnarray}$$

where

(4.13a,b)$$\begin{eqnarray}\unicode[STIX]{x1D702}_{0}=\frac{nT}{\unicode[STIX]{x1D708}}\quad \text{and}\quad \unicode[STIX]{x1D705}_{0}=\frac{d(d+2)}{2(d-1)m}\unicode[STIX]{x1D702}_{0}\end{eqnarray}$$

are the shear viscosity and thermal conductivity, respectively, in the elastic limit; and $\unicode[STIX]{x1D702}^{\ast }$, $\unicode[STIX]{x1D705}^{\ast }$ and $\unicode[STIX]{x1D706}^{\ast }$ are the reduced shear viscosity, reduced thermal conductivity and reduced Dufour-like coefficient, respectively. These reduced transport coefficients are given by

(4.14a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}^{\ast }=\frac{1}{\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }-\frac{1}{2}\unicode[STIX]{x1D701}_{0}^{\ast }}, & \displaystyle\end{eqnarray}$$

(4.14b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D705}^{\ast }=\frac{d-1}{d}\frac{1+2a_{2}}{\unicode[STIX]{x1D708}_{q}^{\ast }-2\unicode[STIX]{x1D701}_{0}^{\ast }}, & \displaystyle\end{eqnarray}$$

(4.14c)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}^{\ast }=\frac{\displaystyle \unicode[STIX]{x1D705}^{\ast }\unicode[STIX]{x1D701}_{0}^{\ast }+\frac{d-1}{d}a_{2}}{\displaystyle \unicode[STIX]{x1D708}_{q}^{\ast }-\frac{3}{2}\unicode[STIX]{x1D701}_{0}^{\ast }}=\frac{\unicode[STIX]{x1D705}^{\ast }}{1+2a_{2}}\frac{\unicode[STIX]{x1D701}_{0}^{\ast }+\unicode[STIX]{x1D708}_{q}^{\ast }a_{2}}{\displaystyle \unicode[STIX]{x1D708}_{q}^{\ast }-\frac{3}{2}\unicode[STIX]{x1D701}_{0}^{\ast }}. & \displaystyle\end{eqnarray}$$

$a_{2}$

$\unicode[STIX]{x1D701}_{0}^{\ast }$

$\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }$

$\unicode[STIX]{x1D708}_{q}^{\ast }$

$\unicode[STIX]{x1D705}^{\ast }$

$\unicode[STIX]{x1D706}^{\ast }$

$a_{2}$

$\unicode[STIX]{x1D701}_{0}^{\ast }$

$\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }$

$\unicode[STIX]{x1D708}_{q}^{\ast }$

(4.15a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}^{\ast }=\frac{8d}{(1+e)[3d+2+(d-2)e]}, & \displaystyle\end{eqnarray}$$

(4.15b)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D705}^{\ast }=\frac{8(d-1)(4d+5-18e+9e^{2})}{(1+e)(d-4+3de)(4d-7+6e-3e^{2})}, & \displaystyle\end{eqnarray}$$

(4.15c)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}^{\ast }=\frac{16(1-e)[2d^{2}+8d-1-6(d+2)e+9e^{2}]}{(1+e)^{2}(d-4+3de)(4d-7+6e-3e^{2})}. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D705}^{\ast }$

$\unicode[STIX]{x1D706}^{\ast }$

$e=(4-d)/(3d)$

$\unicode[STIX]{x1D705}^{\ast }$

$\unicode[STIX]{x1D706}^{\ast }$

$e=1/3$

$d=2$

$e=1/9$

$d=3$

$e<(4-d)/(3d)$

$\unicode[STIX]{x1D705}^{\ast }$

$\unicode[STIX]{x1D706}^{\ast }$

$e\leqslant (4-d)/(3d)$

$T$

$n\sqrt{T}$

(4.16)$$\begin{eqnarray}q_{i}=-\unicode[STIX]{x1D705}^{\prime }\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{i}}-\frac{\unicode[STIX]{x1D706}}{\sqrt{T}}\frac{\unicode[STIX]{x2202}(n\sqrt{T})}{\unicode[STIX]{x2202}x_{i}},\quad \text{where}\quad \unicode[STIX]{x1D705}^{\prime }=\unicode[STIX]{x1D705}-\unicode[STIX]{x1D706}\frac{n}{2T}\end{eqnarray}$$

is referred to as the modified thermal conductivity (Garzó et al. Reference Garzó, Santos and Montanero2007). The reduced modified thermal conductivity $\unicode[STIX]{x1D705}^{\prime \ast }=\unicode[STIX]{x1D705}^{\prime }/\unicode[STIX]{x1D705}_{0}$, using (4.14)–(4.16), is given by

(4.17)$$\begin{eqnarray}\unicode[STIX]{x1D705}^{\prime \ast }=\frac{d-1}{d}\frac{1+\frac{3}{2}a_{2}}{\unicode[STIX]{x1D708}_{q}^{\ast }-\frac{3}{2}\unicode[STIX]{x1D701}_{0}^{\ast }}=\frac{8(2d+1-6e+3e^{2})}{(1+e)^{2}(4d-7+6e-3e^{2})}.\end{eqnarray}$$

The reduced modified thermal conductivity $\unicode[STIX]{x1D705}^{\prime \ast }$ does not possess the above singularity and hence is finite for all $0\leqslant e\leqslant 1$ and for $d=2,3$.

4.3 Comparison with existing theories and computer simulations

I have not found any simulation data on the transport coefficients for IMM. Therefore, in this subsection, I compare the reduced transport coefficients for IMM obtained above with those for IHD ($d=2$) and IHS ($d=3$) obtained through various theoretical and simulation methods.

Figure 1. Reduced shear viscosity $\unicode[STIX]{x1D702}^{\ast }$ plotted over the coefficient of restitution $e$ for (a) two and (b) three dimensions. The thick solid (red) line represents the result for IMM. All other lines or symbols are the results for IHD ($d=2$) or IHS ($d=3$). The thin solid (green) and dash-dotted (magenta) lines delineate the results from the first Sonine and modified first Sonine approximations, respectively (Garzó et al. Reference Garzó, Santos and Montanero2007). The dashed (black) line depicts the result obtained with the G14 distribution function (Garzó Reference Garzó2013) while the solid cyan line depicts that obtained with the G26 equations (Gupta et al. Reference Gupta, Shukla and Torrilhon2018). The squares are the results from the theoretical expressions obtained via the computer-aided method devised by Noskowicz et al. (Reference Noskowicz, Bar-Lev, Serero and Goldhirsch2007). The circles are the DSMC results from Montanero, Santos & Garzó (Reference Montanero, Santos and Garzó2005) and Garzó et al. (Reference Garzó, Santos and Montanero2007).

Figure 2. Reduced thermal conductivity $\unicode[STIX]{x1D705}^{\ast }$ plotted over the coefficient of restitution $e$ for (a) two and (b) three dimensions. The circles are the DSMC results from Brey & Ruiz-Montero (Reference Brey and Ruiz-Montero2004) for $d=2$ and in Brey et al. (Reference Brey, Ruiz-Montero, Maynar and García de Soria2005) for $d=3$. The lines and squares are the same as described in figure 1.

The reduced transport coefficients $\unicode[STIX]{x1D702}^{\ast }$, $\unicode[STIX]{x1D705}^{\ast }$, $\unicode[STIX]{x1D706}^{\ast }$ and $\unicode[STIX]{x1D705}^{\prime \ast }$ for a dilute granular gas are plotted over the coefficient of restitution $e$ in figures 1–4, respectively. Panels (a,b) in each figure exhibit the results for $d=2$ and $d=3$, respectively. The thick solid (red) line in each figure denotes the result for IMM obtained from (4.14) or (4.15) and (4.17), which have been obtained in this paper through the Grad moment equations. Recall that the reduced transport coefficients for IMM obtained through the moment method above are exactly the same as those obtained at first order of the CE expansion (see Santos Reference Santos2003; Garzó & Santos Reference Garzó and Santos2011). Therefore, the thick solid (red) line in each figure also represents the results for IMM from the first-order CE expansion. The remaining lines and symbols in figures 1–4 depict the results for IHD (in case of $d=2$) or for IHS (in case of $d=3$). The thin solid (green) and dash-dotted (magenta) lines are the plots for the reduced transport coefficients from Garzó et al. (Reference Garzó, Santos and Montanero2007) obtained at the first Sonine and modified first Sonine approximations, respectively, in the CE expansion. The dashed (black) lines depict the reduced transport coefficients obtained through the G14 distribution function for $d$-dimensional IHS in Garzó (Reference Garzó2013). Figure 1(b) also displays a solid cyan line, which is not present in the other figures. This solid cyan line is the result for the reduced shear viscosity obtained with the G26 equations for IHS very recently by Gupta et al. (Reference Gupta, Shukla and Torrilhon2018). Indeed, the other transport coefficients from the G14 or G26 equations remain the same; consequently, the solid cyan line in panel (b) of each of figures 2–4 coincides with the dashed black line, and hence has not been shown separately. The squares are the results from the theoretical expressions derived via the so-called computer-aided method devised by Noskowicz et al. (Reference Noskowicz, Bar-Lev, Serero and Goldhirsch2007) and the circles are the numerical solutions of the Boltzmann equation obtained via the direct simulation Monte Carlo (DSMC) method in Brey & Ruiz-Montero (Reference Brey and Ruiz-Montero2004), Montanero et al. (Reference Montanero, Santos and Garzó2005), Brey et al. (Reference Brey, Ruiz-Montero, Maynar and García de Soria2005) and Montanero, Santos & Garzó (Reference Montanero, Santos and Garzó2007). The paper by Garzó et al. (Reference Garzó, Santos and Montanero2007) summarises and presents the DSMC results in the aforementioned references that are computed with two approaches: (i) through Green–Kubo relations in Brey & Ruiz-Montero (Reference Brey and Ruiz-Montero2004) (for $d=2$) and Brey et al. (Reference Brey, Ruiz-Montero, Maynar and García de Soria2005) (for $d=3$), and (ii) by implementing an external force in Montanero et al. (Reference Montanero, Santos and Garzó2005) and Montanero et al. (Reference Montanero, Santos and Garzó2007). The external force in the latter compensates for collisional cooling and yields somewhat better results. Therefore, the DSMC results in figures 1–4 are shown with the latter for whichever coefficient they are available else they are shown with the former – figures 1 and 4(b) display the DSMC results with the latter while figures 2 and 3 show the DSMC results with the former. The DSMC results for the reduced shear viscosity from the latter approach were obtained by Montanero et al. (Reference Montanero, Santos and Garzó2005) for $e=0.6,0.7,0.8,0.9,1$ in the case of $d=3$ while that for $e=0.2,0.3,0.4,0.5$ in the case of $d=3$ and that for all $e$ in the case of $d=2$ were obtained by Garzó et al. (Reference Garzó, Santos and Montanero2007). The DSMC data for $\unicode[STIX]{x1D705}^{\prime \ast }$ in two dimensions (figure 4a) are apparently unavailable.

Figure 3. Reduced coefficient $\unicode[STIX]{x1D706}^{\ast }$ plotted over the coefficient of restitution $e$ for (a) two and (b) three dimensions. The circles are the DSMC results from Brey & Ruiz-Montero (Reference Brey and Ruiz-Montero2004) for $d=2$ and in Brey et al. (Reference Brey, Ruiz-Montero, Maynar and García de Soria2005) for $d=3$. The lines and squares are the same as described in figure 1.

Figure 4. Reduced modified thermal conductivity $\unicode[STIX]{x1D705}^{\prime \ast }$ plotted over the coefficient of restitution $e$ for (a) two and (b) three dimensions. The circles in panel (b) are the DSMC results from Montanero et al. (Reference Montanero, Santos and Garzó2007). The lines and symbols are the same as described in figure 1.

It is clear from figures 1–4 that the IMM model overpredicts all the transport coefficients significantly in comparison to the IHS model, despite the transport coefficients for IMM computed from the Grad moment method in the present work being exactly the same as those obtained through CE expansion for IMM, for example, in Santos (Reference Santos2003) and Garzó & Santos (Reference Garzó and Santos2011). The discrepancies between the results for IMM and IHS are apparently linked to the choice of the effective collision frequency ${\unicode[STIX]{x1D708}\unicode[STIX]{x0030A}}$ (see (2.3a)) in the IMM model (Santos Reference Santos2003), which is chosen in such a way that the cooling rates from the Boltzmann equation for IHS and IMM remain exactly the same. Furthermore, the reduced transport coefficients $\unicode[STIX]{x1D705}^{\ast }$ and $\unicode[STIX]{x1D706}^{\ast }$ for IMM (shown by the thick solid red lines in figures 2 and 3) diverge at $e=1/3$ for $d=2$ (figures 2a and 3a) and at $e=1/9$ for $d=3$ (figures 2b and 3b), and remain unphysical below these values of the coefficient of restitution. On the other hand, the reduced modified thermal conductivity $\unicode[STIX]{x1D705}^{\prime \ast }$ for IMM remains positive for all values of the coefficient of restitution in both dimensions (see figure 4). Nevertheless, $\unicode[STIX]{x1D705}^{\prime \ast }$ for IMM is also much higher than that for IHS. In the case of $d=2$ (figure 4a), $\unicode[STIX]{x1D705}^{\prime \ast }$ for IHS from any theory first decreases then increases on increasing the coefficient of restitution (although the profiles of $\unicode[STIX]{x1D705}^{\prime \ast }$ from the modified first Sonine approximation and first Sonine approximation/G14 theory differ significantly) whereas that for IMM decreases monotonically on increasing the coefficient of restitution. However, as the DSMC data are not available in this case, it is difficult to discern which theory for IHS yields better results in this case.

Among fully theoretical methods, the modified version of the first Sonine approximation (dash-dotted magenta lines) proposed by Garzó et al. (Reference Garzó, Santos and Montanero2007) for IHS seems to be the best model, which captures all the transport coefficient very well, although the G26 model of Gupta et al. (Reference Gupta, Shukla and Torrilhon2018) was able to capture the coefficient of the reduced shear viscosity (but not the other transport coefficients) better than the modified first Sonine approximation.

5.1 Haff’s law

Equations (5.2) and (5.3) with the initial conditions of the HCS imply $n_{\ast }(t_{\ast })=1$ and $v_{i}^{\ast }(t_{\ast })=0$. Therefore, equation (5.4) using the initial conditions of the HCS yields Haff’s law (Haff Reference Haff1983) for the evolution of the granular temperature:

(5.11a,b)$$\begin{eqnarray}T_{\ast }(t_{\ast })=\frac{1}{(1+t_{\ast }/\unicode[STIX]{x1D70F}_{\ast })^{2}}\quad \text{or}\quad T(t)=\frac{T_{0}}{(1+t/\unicode[STIX]{x1D70F}_{0})^{2}},\end{eqnarray}$$

where

(5.12a,b)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\ast }=\frac{2}{\unicode[STIX]{x1D701}_{0}^{\ast }}\quad \text{and}\quad \unicode[STIX]{x1D70F}_{0}=\frac{\unicode[STIX]{x1D70F}_{\ast }}{\unicode[STIX]{x1D708}_{0}}=\frac{2}{\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D708}_{0}}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D70F}_{0}$ is the time scale in Haff’s law for IMM and $\unicode[STIX]{x1D70F}_{\ast }$ is the corresponding dimensionless time scale. Haff’s law (5.11) with time scale (5.12) is exactly the same as that obtained in Garzó & Santos (Reference Garzó and Santos2011) for IMM. It is worthwhile noting that, unlike the energy balance equation in the case of IHS that also contains the scalar fourth moment $\unicode[STIX]{x1D6E5}$ (see e.g. Kremer & Marques Jr. Reference Kremer and Marques2011; Gupta et al. Reference Gupta, Shukla and Torrilhon2018), equation (5.4) does not contain any other moment except $n_{\ast }$ and $T_{\ast }$. Consequently, Haff’s law for IMM does not depend on higher moments; or in other words, Haff’s law remains unchanged for IMM, no matter how large a moment system it is determined from.

Note that the dimensionless time scale $\unicode[STIX]{x1D70F}_{\ast }$ in Haff’s law (5.11) for $d$-dimensional IHS obtained at first approximation of the Sonine expansion is given by (van Noije & Ernst Reference van Noije and Ernst1998)

(5.13)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{\ast }^{(\text{IHS})}=\frac{8d}{(d+2)(1-e^{2})}\left(1+\frac{3}{16}a_{2}^{(\text{IHS})}\right)^{-1},\end{eqnarray}$$

with an excellent estimate for the fourth cumulant (van Noije & Ernst Reference van Noije and Ernst1998)

(5.14)$$\begin{eqnarray}a_{2}^{(\text{IHS})}=\frac{16(1-e)(1-2e^{2})}{24d+9+e(8d-41)+30e^{2}(1-e)}.\end{eqnarray}$$

Figure 5 illustrates the decay of the dimensionless granular temperature in the HCS via Haff’s law (5.11a) for the coefficients of restitution $e=0.75$ (depicted by solid lines and squares) and $e=0.95$ (depicted by dotted lines and circles). The lines denote Haff’s law for IMM, where the associated (dimensionless) time scale $\unicode[STIX]{x1D70F}_{\ast }$ for the decay is given by (5.12a), while the symbols denote that for IHS at first approximation of the Sonine expansion, where the associated (dimensionless) time scale $\unicode[STIX]{x1D70F}_{\ast }$ is given by (5.13). Although the time scale $\unicode[STIX]{x1D70F}_{\ast }$ for IMM does not contain the fourth cumulant $a_{2}$ while that for IHS does contain it, Haff’s law from IMM (lines) is still in very good agreement with that from IHS (symbols) in both two and three dimensions. The granular temperature relaxes faster with increasing inelasticity due to the fact that more inelastic particles dissipate more energy during a collision in comparison to less inelastic ones.

Figure 5. Relaxation of the granular temperature $T_{\ast }$ in the HCS via Haff’s law (5.11a) for (a) $d=2$ and (b) $d=3$ with the initial conditions $n_{\ast }(0)=T_{\ast }(0)=1$. The lines depict granular temperature profiles for IMM, i.e. with $\unicode[STIX]{x1D70F}_{\ast }$ as in (5.12a), while the symbols denote those for IHS, i.e. with $\unicode[STIX]{x1D70F}_{\ast }$ as in (5.13).

5.2 Relaxation of other moments in the HCS

For monatomic gases, it can be shown through the order of magnitude method devised by Struchtrup (Reference Struchtrup2004) that all the non-equilibrium moments ($\unicode[STIX]{x1D70E}_{ij}$, $q_{i}$, $m_{ijk}$, $\unicode[STIX]{x1D6E5}$, $R_{ij}$ and beyond) are at least of first order in spatial gradients (see Struchtrup Reference Struchtrup2004, Reference Struchtrup2005). In contrast, the order of magnitude method in the case of granular gases is not straightforward due to non-conservation of energy and, to the best of my knowledge, has never been attempted so far. Notwithstanding, I would expect that all the higher vectorial and tensorial moments ($\unicode[STIX]{x1D70E}_{ij}$, $q_{i}$, $m_{ijk}$, $R_{ij}$, $\unicode[STIX]{x1D711}_{i}$ and beyond) for granular gases are also at least of first order in spatial gradients; this conjecture is well known for $\unicode[STIX]{x1D70E}_{ij}$ and $q_{i}$ in the case of granular gases as well. This means that all the vectorial and tensorial moments remain zero in the HCS because spatial gradients are zero in this state. Nonetheless, it is interesting to analyse how these higher-order moments relax in the HCS if started with non-vanishing initial conditions.

Equations (5.5)–(5.10) are coupled with (5.2) and (5.4), but can be solved analytically. In order to compare the decay rates of the moments, the initial conditions for all the variables in (5.5)–(5.10) are taken as unity. With these initial conditions, equations (5.5)–(5.10) yield the following solution for the other variables:

(5.15a)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}^{\ast }(t_{\ast })=\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D70F}_{\ast }}, & \displaystyle\end{eqnarray}$$

(5.15b)$$\begin{eqnarray}\displaystyle & \displaystyle q_{i}^{\ast }(t_{\ast })=\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{q}^{\ast }\unicode[STIX]{x1D70F}_{\ast }}, & \displaystyle\end{eqnarray}$$

(5.15c)$$\begin{eqnarray}\displaystyle & \displaystyle m_{ijk}^{\ast }(t_{\ast })=\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{m}^{\ast }\unicode[STIX]{x1D70F}_{\ast }}, & \displaystyle\end{eqnarray}$$

(5.15d)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E5}(t_{\ast })=a_{2}+(1-a_{2})\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6E5}}^{\ast }\unicode[STIX]{x1D70F}_{\ast }}, & \displaystyle\end{eqnarray}$$

(5.15e)$$\begin{eqnarray}\displaystyle & \displaystyle R_{ij}^{\ast }(t_{\ast })=(1-\unicode[STIX]{x1D718}_{R})\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{R}^{\ast }\unicode[STIX]{x1D70F}_{\ast }}+\unicode[STIX]{x1D718}_{R}\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D70F}_{\ast }-2}, & \displaystyle\end{eqnarray}$$

(5.15f)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{i}^{\ast }(t_{\ast })=(1-\unicode[STIX]{x1D718}_{\unicode[STIX]{x1D711}})\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }\unicode[STIX]{x1D70F}_{\ast }}+\unicode[STIX]{x1D718}_{\unicode[STIX]{x1D711}}\left(1+\frac{t_{\ast }}{\unicode[STIX]{x1D70F}_{\ast }}\right)^{-\unicode[STIX]{x1D708}_{q}^{\ast }\unicode[STIX]{x1D70F}_{\ast }-2}, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D718}_{R}=\unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }\unicode[STIX]{x1D70F}_{\ast }/[(\unicode[STIX]{x1D708}_{R}^{\ast }-\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast })\unicode[STIX]{x1D70F}_{\ast }-2]$

$\unicode[STIX]{x1D718}_{\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }\unicode[STIX]{x1D70F}_{\ast }/[(\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }-\unicode[STIX]{x1D708}_{q}^{\ast })\unicode[STIX]{x1D70F}_{\ast }-2]$

$d=3$

$e=1$

(5.16)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70E}_{ij}^{\ast }(t_{\ast })=\text{e}^{-t_{\ast }},\quad q_{i}^{\ast }(t_{\ast })=\text{e}^{-(2/3)t_{\ast }},\quad m_{ijk}^{\ast }(t_{\ast })=\text{e}^{-(3/2)t_{\ast }},\\ \unicode[STIX]{x1D6E5}(t_{\ast })=\text{e}^{-(2/3)t_{\ast }},\quad R_{ij}^{\ast }(t_{\ast })=\text{e}^{-(7/6)t_{\ast }},\quad \unicode[STIX]{x1D711}_{i}^{\ast }(t_{\ast })=\text{e}^{-t_{\ast }}.\end{array}\right\}\end{eqnarray}$$

From solution (5.16), it is clear that, for dilute monatomic gases of Maxwell molecules, the third-order moment $m_{ijk}^{\ast }$ relaxes faster than all other higher-order moments considered in the present work; $R_{ij}^{\ast }$ relaxes slower than $m_{ijk}^{\ast }$ but faster than the remaining moments ($\unicode[STIX]{x1D70E}_{ij}^{\ast }$, $q_{i}^{\ast }$, $\unicode[STIX]{x1D6E5}$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$); $\unicode[STIX]{x1D70E}_{ij}^{\ast }$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$ relax with equal relaxation rates but faster than $q_{i}^{\ast }$ and $\unicode[STIX]{x1D6E5}$, which also relax with equal relaxation rates.

Figure 6. Relaxation of the other moments – $\unicode[STIX]{x1D70E}_{ij}^{\ast }$, $q_{i}^{\ast }$, $m_{ijk}^{\ast }$, $\unicode[STIX]{x1D6E5}$, $R_{ij}^{\ast }$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$ – for IMM in the HCS for (a) $d=2$ and (b) $d=3$ evaluated using analytical solutions (5.15) and $\unicode[STIX]{x1D70F}_{\ast }$ from (5.12a). The coefficient of restitution is taken as $e=0.95$ (main panels) and $e=0.25$ (insets). Initial conditions are taken as $n_{\ast }(0)=T_{\ast }(0)=\unicode[STIX]{x1D70E}_{ij}^{\ast }(0)=q_{i}^{\ast }(0)=m_{ijk}^{\ast }(0)=\unicode[STIX]{x1D6E5}(0)=R_{ij}^{\ast }(0)=\unicode[STIX]{x1D711}_{i}^{\ast }(0)=1$.

Figure 6 illustrates the relaxation of the other (dimensionless) moments – $\unicode[STIX]{x1D70E}_{ij}^{\ast }$, $q_{i}^{\ast }$, $m_{ijk}^{\ast }$, $\unicode[STIX]{x1D6E5}$, $R_{ij}^{\ast }$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$ over the (dimensionless) time $t_{\ast }$ (via analytical solutions (5.15)) in two and three dimensions for granular gases (i.e. for $e\neq 1$): $e=0.95$ (main panels) and $e=0.25$ (insets). It turns out that all these moments relax with time much faster than the granular temperature. It is interesting to note that all the vectorial and tensorial moments relax to zero whereas the scalar moment $\unicode[STIX]{x1D6E5}$ relaxes to a non-zero value $a_{2}$ for all $e\neq 1$. These can also be verified from analytical solutions (5.15) in the limit $t_{\ast }\rightarrow \infty$, since all the exponents in (5.15) are negative for all $e\neq 1$ (note that $\unicode[STIX]{x1D70F}_{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D70E}}^{\ast }$, $\unicode[STIX]{x1D708}_{q}^{\ast }$, $\unicode[STIX]{x1D708}_{m}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6E5}}^{\ast }$, $\unicode[STIX]{x1D708}_{R}^{\ast }$, $\unicode[STIX]{x1D708}_{R\unicode[STIX]{x1D70E}}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}}^{\ast }$, $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D711}q}^{\ast }$, $\unicode[STIX]{x1D6FC}_{0}$, $\unicode[STIX]{x1D6FC}_{1}$, $\unicode[STIX]{x1D6FC}_{2}$ and $\unicode[STIX]{x1D6FC}_{3}$ are positive for $e\neq 1$). For large values of the coefficient of restitution (main panels), all these moments decay monotonically and their decay rates are, apparently, proportional to their tensorial orders, i.e. the third-order moment ($m_{ijk}^{\ast }$) decays faster than the second-order moments ($\unicode[STIX]{x1D70E}_{ij}^{\ast }$ and $R_{ij}^{\ast }$), which decay faster than the vectorial moment ($q_{i}^{\ast }$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$), which decay faster than the scalar moment ($\unicode[STIX]{x1D6E5}$). However, for sufficiently small values of the coefficient of restitution (insets), the higher-order moments ($R_{ij}^{\ast }$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$) do not decay monotonically. This is attributed to the fact that higher-order (sixth order and beyond) moments for IMM are prone to diverge for sufficiently small values of the coefficient of restitution in the HCS (Santos & Garzó Reference Santos and Garzó2012) and it is manifested already through the non-monotonic relaxation of $R_{ij}^{\ast }$ and $\unicode[STIX]{x1D711}_{i}^{\ast }$ for small coefficients of restitution (insets), although they themselves do not diverge.

6.1 Eigenmodes from the longitudinal and transverse systems

The nine roots of dispersion relation (6.16a) lead to nine eigenmodes for the longitudinal system (6.15a) associated with the G29 equations while the six roots of dispersion relation (6.16b) yield six eigenmodes for the transverse system (6.15b) associated with the G29 equations. The real part of the frequency ($\text{Re}(\unicode[STIX]{x1D714})$) associated with each eigenmode and its growth rate ($\text{Im}(\unicode[STIX]{x1D714})$) from both the longitudinal and transverse systems are illustrated in figures 7–10 for $d=2$ (figures 7 and 9), $d=3$ (figures 8 and 10) – with figures 7 and 8 being for the longitudinal system (6.15a) and figures 9 and 10 for the transverse system (6.15b). Panels (a,b) and (c,d) in each figure depict the eigenmodes for the inelastic ($e=0.75$) and elastic ($e=1$) cases, respectively, while panels (a,c) and (b,d) in each figure delineate the real part of the frequency and growth rate, respectively. It is apparent from figures 7(a,c) and 8(a,c) that four pairs out of the nine eigenmodes from the longitudinal system have non-zero $\text{Re}(\unicode[STIX]{x1D714})$, i.e. the four pairs of associated eigenmodes are travelling whereas one eigenmode is purely imaginary, i.e. it has $\text{Re}(\unicode[STIX]{x1D714})=0$ for all wavenumbers and hence always remains stationary. Similarly, it is clear from figures 9(a,c) and 10(a,c) that two pairs out of the six eigenmodes from the transverse system have non-zero $\text{Re}(\unicode[STIX]{x1D714})$, i.e. the two pairs of associated eigenmodes are travelling, whereas two eigenmodes are purely imaginary and hence remain stationary for all wavenumbers. A travelling eigenmode is commonly referred to as a sound mode and a stationary eigenmode as a heat mode (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004; Garzó Reference Garzó2005).

6.2 Comparison among various Grad moment theories

As discussed in § 3.5, a lower-level Grad moment system can be obtained from the G29 equations by discarding the appropriate field variables. Accordingly, the longitudinal and transverse problems associated with the G13, G14 and G26 systems can be obtained by eliminating the appropriate variables and corresponding rows and columns of the matrices $\mathscr{L}$ and $\mathscr{T}$ in (6.15). For comparison purposes, I shall also include the results obtained from the linear stability analysis of the system of the NSF equations for IMM along with these Grad moment systems. The linear–dimensionless NSF equations in the perturbed field variables are (6.4)–(6.6) with

(6.17a,b)$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70E}}_{ij}=-2\unicode[STIX]{x1D702}^{\ast }\frac{\unicode[STIX]{x2202}\tilde{v}_{\langle \!i}}{\unicode[STIX]{x2202}\tilde{x}_{j\!\rangle }}\quad \text{and}\quad \tilde{q}_{i}=-\frac{d(d+2)}{2(d-1)}\left(\unicode[STIX]{x1D705}^{\ast }\frac{\unicode[STIX]{x2202}\tilde{T}}{\unicode[STIX]{x2202}\tilde{x}_{i}}+\unicode[STIX]{x1D706}^{\ast }\frac{\unicode[STIX]{x2202}{\tilde{n}}}{\unicode[STIX]{x2202}\tilde{x}_{i}}\right),\end{eqnarray}$$

where the reduced transport coefficients $\unicode[STIX]{x1D702}^{\ast }$, $\unicode[STIX]{x1D705}^{\ast }$ and $\unicode[STIX]{x1D706}^{\ast }$ for IMM are given by (4.14). From the linear–dimensionless NSF equations, it is straightforward to obtain the longitudinal and transverse problems associated with the system of the NSF equations by following a similar procedure as above. The longitudinal and transverse problems associated with the system of the NSF equations read

(6.18a,b)$$\begin{eqnarray}\mathscr{L}_{NSF}\left[\begin{array}{@{}c@{}}{\check{n}}\\ \check{v}_{x}\\ {\check{T}}\end{array}\right]=\left[\begin{array}{@{}c@{}}0\\ 0\\ 0\end{array}\right]\quad \text{and}\quad \left(-\unicode[STIX]{x1D714}-i\unicode[STIX]{x1D702}^{\ast }k^{2}+i{\displaystyle \frac{\unicode[STIX]{x1D701}_{0}^{\ast }}{2}}\right)\check{v}_{y}=0,\end{eqnarray}$$

respectively, where the matrix $\mathscr{L}_{NSF}\equiv \mathscr{L}_{NSF}(k,\unicode[STIX]{x1D714},d,e)$ is also presented in appendix F. As far as the linear stability of the HCS is concerned, the NSF theory and all Grad moment theories – although not shown here explicitly for the NSF, G13, G14 and G26 equations – predict a similar behaviour for moderate to large values of the coefficient of restitution in the sense that one heat mode from the longitudinal system associated with each theory and one shear mode from the transverse system associated with each theory are unstable below some critical wavenumbers for granular gases while all the modes remain stable in the elastic case ($e=1$). The stability of a (longitudinal or transverse) system is regulated by its least stable eigenmode. Therefore, to analyse the stability of the longitudinal and transverse systems associated with different moment theories, the critical wavenumbers for the least stable mode (unstable shear mode for the longitudinal system and unstable heat mode for the transverse system) from each moment theory are plotted in the $(e,k_{c})$-plane in figure 11. The figure also includes the critical wavenumber profiles (shown by thin dashed black lines) for the least stable modes of the longitudinal and transverse systems associated with the NSF theory for IMM. Figures 11(a,b) and 11(c,d) display the critical wavenumbers for the longitudinal and transverse systems, respectively, while figures 11(a,c) and 11(b,d) exhibit the results for $d=2$ and $d=3$, respectively. Since the critical wavenumber is that wavenumber where the growth rate ($\text{Im}(\unicode[STIX]{x1D714})$) changes its sign, the curves in figure 11 are essentially the zero contours – in the $(e,k_{c})$-plane – of the growth rate of the least stable mode in each system. Each curve for the critical wavenumber corresponds to $\text{Im}(\unicode[STIX]{x1D714})=0$ and hence divides the $(e,k_{c})$-plane into two parts demarcating the stable and unstable regions: the region below a curve is unstable as $\text{Im}(\unicode[STIX]{x1D714})>0$ in this region whereas that above this curve is stable as $\text{Im}(\unicode[STIX]{x1D714})<0$ in this region. In general, the NSF and all the moment theories considered here predict qualitatively similar critical wavenumber profiles. In particular, the critical wavenumber is zero in the elastic ($e=1$) case since both the longitudinal and transverse systems are stable in this case, and it increases with increasing inelasticity, in general. For nearly elastic granular gases ($0.9\leqslant e\leqslant 1$), the respective critical wavenumber profiles from the NSF and all moment theories coincide for both the longitudinal and transverse systems.

Figure 11. Critical wavenumbers in the $(e,k_{c})$-plane from the NSF and different Grad moment theories. Panels (a,b) and (c,d) exhibit the critical wavenumbers for the longitudinal and transverse systems, respectively, while (a,c) and (b,d) display the results for $d=2$ and $d=3$, respectively. Each system is unstable (stable) below (above) its corresponding curve.

Table 2. Threshold values of the coefficient of restitution $e_{th}$ below which the longitudinal and transverse systems for IMM remain always unstable.

For the longitudinal system (figure 11a,b), the NSF and all the moment theories yield smooth critical wavenumber profiles for moderate to large values of the coefficient of restitution discerning the stability regions. However, each theory gives a threshold value of the coefficient of restitution below which the longitudinal system remains unstable (because for $e<e_{th}$, the growth rates of some of the eigenmodes of the longitudinal system remain positive even for large wavenumbers), which is not true (see e.g. Gupta et al. Reference Gupta, Shukla and Torrilhon2018). This simply means that the IMM model is not adequate for granular gases with moderate to large inelasticity. The threshold values of the coefficient of restitution $e_{th}$ for the longitudinal systems associated with the NSF and different Grad moment theories are given in table 2. Furthermore, the G13, G14 and G26 theories lead to kinks at $e\approx 0.7,0.648,0.608$ in the case of $d=2$ (at $e\approx 0.562,0.522,0.446$ in the case of $d=3$), respectively, indicating that the region of applicability increases on increasing the number of moments.

For the transverse system (figure 11c,d), the critical wavenumbers for the G13 and G14 theories are exactly the same due to the fact that the scalar fourth moment $\check{\unicode[STIX]{x1D6E5}}$ does not enter the transverse system associated with any moment theory (see (6.15b)). Hence, the critical wavenumbers from both the theories are depicted by a single dot-dashed magenta line. Among the theories considered, the transverse systems associated with the NSF and G13 (or G14) theories give the critical wavenumbers for all coefficients of restitution whereas those associated with the G26 and G29 theories again lead to threshold values of the coefficient of restitution below which the transverse systems associated with them remain unstable for all wavenumbers (due to the same reason as above). This again restricts the employability of the IMM model to moderately to nearly elastic granular gases. The threshold values of the coefficient of restitution $e_{th}$ for the transverse systems associated with the G26 and G29 theories are also given in table 2. The critical wavenumbers from the G26 (shown by dashed blue line) and G29 (shown by solid red line) theories closely follow each other for $0.7\lesssim e\leqslant 1$ in the case of $d=2$ and for $0.6\lesssim e\leqslant 1$ in the case of $d=3$. Similarly, the critical wavenumbers from the NSF (shown by thin dashed black line) and G13 or G14 (shown by dot-dashed magenta line) theories closely follow each other for $0.7\lesssim e\leqslant 1$ in the case of $d=2$ and for all values of $e$ in the case of $d=3$.

From figures 7–11, it is concluded that the instabilities of the longitudinal and transverse systems above – for moderately to nearly elastic granular gases – are confined to small wavenumbers (or long wavelengths), i.e. these instabilities are long-wave instabilities. Thus, there is a minimum system size, referred to as the critical system size, such that the instabilities will not appear in a system having size smaller than the critical system size.

6.3 Critical system size

It is well established – through the linear stability analysis of hydrodynamic models, through the DSMC method as well as through molecular dynamics (MD) simulations – that the HCS of a freely cooling granular gas is unstable but a minimum critical system size is necessary for the onset of instabilities (see e.g. Brey et al. Reference Brey, Ruiz-Montero and Moreno1998b; Brilliantov & Pöschel Reference Brilliantov and Pöschel2004; Garzó Reference Garzó2005; Gupta et al. Reference Gupta, Shukla and Torrilhon2018). Moreover, it is also known that during the instability phenomenon of the HCS, the instability of the unstable shear mode first engenders the formation of vortices in the system through linear effects (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004; Garzó Reference Garzó2005) and subsequently clustering set in due to the instability of the shear mode through nonlinear effects (Brey, Ruiz-Montero & Cubero Reference Brey, Ruiz-Montero and Cubero1999; Goldhirsch Reference Goldhirsch2003). In order to verify through the moment theories that it is the unstable shear mode which is responsible for the onset of instabilities in a freely cooling granular gas, the critical wavenumbers for the longitudinal and transverse systems are plotted together in the $(e,k_{c})$-plane (but only for that range of $e$, which apparently leads to meaningful critical wavenumber profiles) in figure 12. Denoting the critical wavenumbers associated with the unstable heat and shear modes for any moment system by $k_{h}$ and $k_{s}$, respectively, it can be easily deduced that a heat (shear) mode of the longitudinal (transverse) system for wavenumbers $k>k_{h}$ ($k>k_{s}$) will always decay while that for wavenumbers $k<k_{h}$ ($k<k_{s}$) will grow exponentially. For the ranges of the coefficient of restitution shown in figure 12, $k_{s}\geqslant k_{h}$. Since the wavelength, and hence the critical system size, is inversely proportional to the wavenumber, it is the unstable shear mode from the transverse system which becomes unstable first. From the above discussion, the critical system size can be determined with the knowledge of the critical wavenumber for the transverse system itself. Indeed, it is possible to obtain the analytical expressions for the critical wavenumbers from the moment theories as follows.

Figure 12. Critical wavenumbers for the longitudinal and transverse systems – associated with the NSF and different Grad moment theories (G13, G14, G26, G29) – plotted together in the $(e,k_{c})$-plane. Panels (a,c,e,g) and (b,d,f,h) display the results for $d=2$ and $d=3$, respectively.

It has been verified (although shown here only for the G29 system in panels (a,b) of figures 7–10 in the case of $e=0.75$ for the sake of succinctness) that the real part of the frequency for the least stable eigenmode is either always zero (for longitudinal systems associated with the NSF, G14 and G26 theories and for transverse systems associated with all the theories) or is non-zero only for wavenumbers above the critical wavenumber (for longitudinal systems associated with the G13 and G29 theories). Therefore, for the least stable mode, $\unicode[STIX]{x1D714}=0$ at critical wavenumber. Consequently, the critical wavenumbers for the longitudinal and transverse systems associated with a moment theory can also be determined by substituting $\unicode[STIX]{x1D714}=0$ in their dispersion relations and solving the resulting equations for $k$. For instance, the critical wavenumbers for the longitudinal and transverse systems associated with the G29 theory can also be determined by inserting $\unicode[STIX]{x1D714}=0$ in (6.16) and solving the resulting equations, namely $a_{9}=0$ and $b_{6}=0$. Hence, the roots of $a_{9}=0$ and $b_{6}=0$, respectively, yield the critical wavenumbers for the longitudinal and transverse systems associated with the G29 equations. It is worthwhile noting that the coefficients $a_{9}$ and $b_{6}$ are (even-degree) polynomials of degree eight and four in $k$, respectively. Similarly, the coefficients of $\unicode[STIX]{x1D714}^{0}$ in the dispersion relations for the longitudinal (transverse) systems associated with the NSF, G13, G14 and G26 are (even-degree) polynomials of degree four, four, four and six (two, two, two and four) in $k$, respectively. Consequently, each of them leads to more than one value of the critical wavenumber. Nevertheless, only one of them in each case is meaningful (positive) and that is the analytical expression for the corresponding critical wavenumber. After some algebra, the explicit expressions for $k_{h}$ and $k_{s}$ for each of the NSF and Grad moment systems, in a compact form, can be written as

(6.19)$$\begin{eqnarray}\displaystyle & \displaystyle k_{h|\text{NSF}}=\sqrt{\frac{d-1}{2(d+2)}}\sqrt{\frac{\unicode[STIX]{x1D701}_{0}^{\ast }}{\unicode[STIX]{x1D705}^{\ast }-\unicode[STIX]{x1D706}^{\ast }}}, & \displaystyle\end{eqnarray}$$

(6.20)$$\begin{eqnarray}\displaystyle & \displaystyle k_{s|\text{NSF}}=\sqrt{\frac{\unicode[STIX]{x1D701}_{0}^{\ast }}{2\unicode[STIX]{x1D702}^{\ast }}}, & \displaystyle\end{eqnarray}$$

(6.21)$$\begin{eqnarray}\displaystyle & \displaystyle k_{h|\text{G13}}=\sqrt{\frac{d(d+2)}{2}}\sqrt{\frac{\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D709}_{q}}{\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D709}_{3}+(d+2)^{2}(1+a_{2})\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D70E}}}}, & \displaystyle\end{eqnarray}$$

(6.22)$$\begin{eqnarray}\displaystyle & \displaystyle k_{h|\text{G14}}=\sqrt{\frac{d(d+2)}{2}}\sqrt{\frac{\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D709}_{q}\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6E5}}^{\ast }}{\unicode[STIX]{x1D709}_{4}}}, & \displaystyle\end{eqnarray}$$

(6.23)$$\begin{eqnarray}\displaystyle & \displaystyle k_{s|\text{G13}}=k_{s|\text{G14}}=\sqrt{\frac{d+2}{2}}\sqrt{\frac{\unicode[STIX]{x1D701}_{0}^{\ast }\unicode[STIX]{x1D709}_{\unicode[STIX]{x1D70E}}\unicode[STIX]{x1D709}_{q}}{(d+2)\unicode[STIX]{x1D709}_{q}-\unicode[STIX]{x1D701}_{0}^{\ast }}}, & \displaystyle\end{eqnarray}$$

(6.24)$$\begin{eqnarray}\displaystyle & \displaystyle k_{h|\text{G26}}=\sqrt{\frac{\sqrt{\unicode[STIX]{x1D717}_{12}^{2}+\unicode[STIX]{x1D717}_{11}\unicode[STIX]{x1D717}_{13}}-\unicode[STIX]{x1D717}_{12}}{\unicode[STIX]{x1D717}_{11}}}, & \displaystyle\end{eqnarray}$$

(6.25)$$\begin{eqnarray}\displaystyle & \displaystyle k_{s|\text{G26}}=\sqrt{\frac{\sqrt{\unicode[STIX]{x1D717}_{22}^{2}+\unicode[STIX]{x1D717}_{21}\unicode[STIX]{x1D717}_{23}}-\unicode[STIX]{x1D717}_{22}}{\unicode[STIX]{x1D717}_{21}}}, & \displaystyle\end{eqnarray}$$

(6.26)$$\begin{eqnarray}\displaystyle & \displaystyle k_{h|\text{G29}}=\sqrt{\frac{1}{\unicode[STIX]{x1D717}_{31}}\left(\frac{\unicode[STIX]{x1D717}_{33}}{\unicode[STIX]{x1D717}_{34}}+\unicode[STIX]{x1D717}_{34}-\unicode[STIX]{x1D717}_{32}\right)}, & \displaystyle\end{eqnarray}$$

(6.27)$$\begin{eqnarray}\displaystyle & \displaystyle k_{s|\text{G29}}=\sqrt{\frac{\sqrt{\unicode[STIX]{x1D717}_{42}^{2}+\unicode[STIX]{x1D717}_{41}\unicode[STIX]{x1D717}_{43}}-\unicode[STIX]{x1D717}_{42}}{\unicode[STIX]{x1D717}_{41}}}, & \displaystyle\end{eqnarray}$$

where the coefficients appearing in these expressions are relegated to appendix G for better readability; and ‘$|\text{NSF}$’ and ‘$|\text{G}\cdots \,$’ in subscripts denote the NSF and Grad moment systems which the critical wavenumbers belong to. The plots of $k_{h}$ and $k_{s}$ from the analytical expressions agree completely with those in figure 12 (at least for the depicted ranges of the coefficient of restitution; for instance, $k_{h|\text{G29}}$ becomes complex for smaller values of $e$).

The critical system size $L_{c}$ is estimated by $\tilde{L}_{c}=2\unicode[STIX]{x03C0}/\max \{k_{h},k_{s}\}$ (Garzó Reference Garzó2005; Gupta et al. Reference Gupta, Shukla and Torrilhon2018), where $\tilde{L}_{c}:=L_{c}/\ell$ is the dimensionless critical system size and $\ell$ is the length scale defined in (6.2). From figure 12, $k_{s}>k_{h}$ for the NSF and all the moment theories considered in this work (for the shown ranges of $e$ in the figure); therefore, the critical system size $L_{c}$ is given by

(6.28)$$\begin{eqnarray}L_{c}=\frac{2\unicode[STIX]{x03C0}}{k_{s}}\times \ell =\frac{2\unicode[STIX]{x03C0}}{k_{s}}\times \frac{d+2}{4\sqrt{2}}\frac{\displaystyle \unicode[STIX]{x1D6E4}\left(\frac{d}{2}\right)}{\displaystyle \unicode[STIX]{x1D6E4}\left(\frac{d+1}{2}\right)}\ell _{0},\end{eqnarray}$$

where $\ell _{0}=[\unicode[STIX]{x1D6E4}\{(d+1)/2\}]/[\sqrt{2}\,\unicode[STIX]{x1D70B}^{(d-1)/2}n_{0}d^{d-1}]$ is the mean free path of a dilute hard-sphere gas.

Figure 13. Critical system size in units of the mean free path $\ell _{0}$ plotted over the coefficient of restitution $e$ for (a) $d=2$ and (b) $d=3$. The three-dimensional MD simulation results of Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011) (depicted by triangles) in panel (b) are at solid fraction $\unicode[STIX]{x1D719}=0.1$ while the dilute limit refers to $\unicode[STIX]{x1D719}\rightarrow 0$. The cyan line in panel (b) represents the critical system size at solid fraction $\unicode[STIX]{x1D719}=0.1$ computed from the theoretical expression of Garzó (Reference Garzó2005), and is included only to show the good agreement between the theoretical results of Garzó (Reference Garzó2005) and MD simulations results of Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011).

The critical system size in units of the mean free path ($L_{c}/\ell _{0}$) is illustrated in figure 13 as a function of the coefficient of restitution $e$ for $d=2$ (panel a) and $d=3$ (panel b). The dot-dashed (magenta), dashed (blue) and solid (red) lines in the figure depict the critical system size from the G13 (or G14), G26 and G29 theories, respectively, computed with formula (6.28) using the analytical expressions for $k_{s}$ given in (6.23), (6.25) and (6.27), respectively. For comparison purposes, the figure also includes thin solid black lines depicting the critical system size from the NSF theory for IMM computed with formula (6.28) using the analytical expressions for $k_{s}$ given in (6.20). The dashed (green) lines with symbols delineate the critical system size computed from the theoretical expression given in Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b), which was obtained via the linear stability analysis of a kinetic model for granular gases of IHS due to Brey, Dufty & Santos (Reference Brey, Dufty and Santos1997). It should be noted that the $l_{0}$ used in Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b) relates to the mean free path $\ell _{0}$ used in the present work via

$$\begin{eqnarray}l_{0}=\frac{2\sqrt{2}}{C}\frac{\unicode[STIX]{x03C0}^{(d/2)-1}}{\displaystyle \unicode[STIX]{x1D6E4}\left(\frac{d+1}{2}\right)}\ell _{0}\quad \text{with }C\simeq \left\{\begin{array}{@{}ll@{}}2\quad & \text{for}\quad d=2,\\ {\textstyle \frac{16}{5}}\quad & \text{for}\quad d=3.\end{array}\right.\end{eqnarray}$$

The dotted (black) line in figure 13(b) denotes the critical system size determined from the theoretical expression obtained by the linear stability analysis of the granular NSF equations for IHS in Garzó (Reference Garzó2005). The (red) circles in figure 13(a) are the results from two-dimensional DSMC simulations carried out by Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b). The triangles in figure 13(b) delineate the results from MD simulations of IHS carried out by Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011) at solid fraction $\unicode[STIX]{x1D719}=0.1$ and are included only for qualitative comparison, since the present work deals with dilute granular gases for which $\unicode[STIX]{x1D719}\rightarrow 0$. Note that the critical system size in Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011) is scaled with the diameter of a particle while that in the present work is scaled with the mean free path. Therefore, the MD simulations results of Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011) have been multiplied by a factor $6\sqrt{2}\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}(\unicode[STIX]{x1D719})$ while displaying them in figure 13. Here, $\unicode[STIX]{x1D712}(\unicode[STIX]{x1D719})=(2-\unicode[STIX]{x1D719})/[2(1-\unicode[STIX]{x1D719})^{3}]$ is the pair correlation function. Figure 13(b) also illustrates a solid cyan line, which depicts the critical system size for $\unicode[STIX]{x1D719}=0.1$ computed with the theoretical expressions derived in Garzó (Reference Garzó2005), and is included just to show the good agreement between the theoretical results of Garzó (Reference Garzó2005) and MD simulations results of Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011).

Figure 13 reveals that the critical system size from all the theories and simulations decreases with increasing inelasticity. For disk flows ($d=2$, figure 13a), the critical system size from the NSF theory for IMM (thin solid black line) agrees well with that from the theoretical expression in Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b). However, the critical system size from the G13 or G14 theory (dashed magenta line) seems to be slightly better than that from the NSF theory and agrees perfectly with that from the theoretical expression in Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b) (dashed green line with symbols), and also agrees reasonably well with the DSMC results of Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b) (red circles) – for $0.65\leqslant e\leqslant 1$. Nevertheless, the G26 and G29 theories somewhat underpredict the critical system size for all coefficients of restitution $e\gtrsim 0.65$, although their predictions are also close to the other theories for $e\gtrsim 0.9$.

For sphere flows ($d=3$, figure 13b), I could not find any simulation data for the dilute limit, i.e. for solid fraction $\unicode[STIX]{x1D719}\rightarrow 0$. Therefore, the data from MD simulations carried out by Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011) for solid fraction $\unicode[STIX]{x1D719}=0.1$ are included for comparison. It is important to note that the results from MD simulations by Mitrano et al. (Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011) are in good agreement with those from theoretical expression of Garzó (Reference Garzó2005) not only for $\unicode[STIX]{x1D719}=0.1$ but also for $\unicode[STIX]{x1D719}=0.4$ (see Mitrano et al. Reference Mitrano, Dahl, Cromer, Pacella and Hrenya2011, figure 9). Therefore, in the dilute limit ($\unicode[STIX]{x1D719}\rightarrow 0$), the results from the theoretical expression of Garzó (Reference Garzó2005), shown by the dotted black line in figure 13, can be treated as a benchmark. Additionally, in this limit, the results from the theoretical expressions of Garzó (Reference Garzó2005) and Brey et al. (Reference Brey, Ruiz-Montero and Moreno1998b) are also in good agreement with each other. Clearly, the critical system size from the NSF theory for IMM is again in reasonably good agreement with that from Garzó (Reference Garzó2005). Moreover, the critical system sizes from all the moment theories are also in good agreement with that from Garzó (Reference Garzó2005) for $e\gtrsim 0.85$, but deviate slightly from the results of Garzó (Reference Garzó2005) for $0.55\lesssim e\lesssim 0.85$, where the G26 and G29 theories again underpredict the critical system size while the G13 or G14 theory overpredicts it. Between the G26 and G29 theories, the latter seems to perform slightly better at moderate values of the coefficient of restitution for both $d=2$ and $d=3$.

From figure 13, it is apparent that the critical system size obtained from the NSF and Grad moment theories for IMM is in qualitatively good agreement with that obtained from the NSF-level theories and simulations for IHD/IHS, although it is also noticeable from the figure that some lower-order Grad moment theories (e.g. the G13 or G14 theory) perform better than some higher-order Grad moment theories (e.g. the G26 theory). A possible reason for this could be the choice of the effective collision frequency ${\unicode[STIX]{x1D708}\unicode[STIX]{x0030A}}$ in the IMM model (see (2.3a)) that was chosen in such a way that the cooling rates for IHS and IMM remain exactly the same while the collisional production terms in the other moment equations for IMM follow accordingly based on this choice of the effective collision frequency. Consequently, with this choice of the effective collision frequency, even the NSF equations for IMM seem to perform better than some higher-order moment models. Therefore, it would be interesting to explore other possible choices for the effective collision frequency in the future in such a way that the results from the Boltzmann equations for IHS and IMM agree in an optimal way so that a higher-order moment model would perform better than a lower-order moment model in the case of IMM, similarly to moment models for IHS.