Exact Transport Coefficients from the Inelastic Rough Maxwell Model of a Granular Gas

Abstract

Granular gases demand models capable of capturing their distinct characteristics. The widely employed inelastic hard-sphere model (IHSM) introduces complexities that are compounded when incorporating realistic features like surface roughness and rotational degrees of freedom, resulting in the more intricate inelastic rough hard-sphere model (IRHSM). This paper focuses on the inelastic rough Maxwell model (IRMM), presenting a more tractable alternative to the IRHSM and enabling exact solutions. Building on the foundation of the inelastic Maxwell model (IMM) applied to granular gases, the IRMM extends the mathematical representation to encompass surface roughness and rotational degrees of freedom. The primary objective is to provide exact expressions for the Navier–Stokes–Fourier transport coefficients within the IRMM, including the shear and bulk viscosities, the thermal and diffusive heat conductivities, and the cooling-rate transport coefficient. In contrast to earlier approximations in the IRHSM, our study unveils inherent couplings, such as shear viscosity to spin viscosity and heat conductivities to counterparts associated with a torque-vorticity vector. These exact findings provide valuable insights into refining the Sonine approximation applied to the IRHSM, contributing to a deeper understanding of the transport properties in granular gases with realistic features.

1 Introduction

While it is well known that a hydrodynamic description can be applied to granular gases [1,2,3], the derivation of the associated Navier–Stokes–Fourier (NSF) transport coefficients faces important challenges. To begin with, the spatially uniform and isotropic base state is not that of equilibrium (as in the case of molecular gases) but either the so-called homogeneous cooling state (HCS) in the absence of any external driving or a nonequilibrium steady state if energy is injected to the granular gas by a certain thermostat. In top of that, one usually needs to leave out some realistic properties of the grains and focus on their most basic features by choosing a physically relevant model. In the most widely used model, the grains are assumed to be represented by perfectly smooth spheres characterized by a constant coefficient of normal restitution, $\alpha$ [4]. Despite that, the derivation of the NSF transport coefficients in this inelastic hard-sphere model (IHSM) requires the use of nontrivial approximations [3, 5,6,7]. Those approximations become even more problematic if the hard-sphere model is made more realistic by assuming surface roughness, so that the rotational degrees of freedom need to be included in the description. The resulting inelastic rough hard-sphere model (IRHSM) is characterized by a coefficient of tangential restitution, $\beta$ [8,9,10,11,12,13,14,15,16], in addition to $\alpha$. Needless to say, the NSF transport coefficients obtained from the IRHSM by an extension of the approximations used for the IHSM are much more involved than in the latter case [17,18,19,20].

The importance of exactly solvable models in statistical physics cannot be overemphasized, as they play a pivotal role in our understanding of complex physical systems and the development of fundamental principles. These models function as essential benchmarks, offering precise solutions that can serve as a reference for more realistic but less tractable systems Moreover, they are invaluable for testing theoretical concepts, verifying numerical simulations, and confirming experimental observations. In this context, Maxwell models applied to the Boltzmann equation serve as essential tools in the development of the kinetic theory of gases.

As is well known, the original Maxwell molecules interact via a continuous repulsive interaction potential decaying as $r^{-4}$ with the intermolecular distance r [21,22,23,24]. A notable consequence of this behavior is the independence of the collision rate on the relative velocity of the colliding pair. Consequently, the velocity moments of the Boltzmann collision term can be expressed as bilinear combinations of moments of the distribution function of equal or lower degree, facilitating the derivation of insightful exact results. The extension of this velocity-independence of the collision rate allows for the formulation of Maxwellian mathematical models that do not rely on a specific interaction potential [25].

While Maxwell models were originally designed for molecular gases, the exploration of granular gases has led to the development of the inelastic Maxwell model (IMM), accompanied by the derivation of various exact results, particularly focusing on the presence of fat high-energy tails [26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52]. Within the IMM framework, precise evaluations of the exact NSF transport coefficients [53,54,55,56,57,58] and the Burnett coefficients [59] have been conducted. Additionally, several rheological properties have been determined [60,61,62,63,64,65,66,67,68,69,70,71]. The interest in the IMM extends beyond its mathematical tractability, encompassing its ability to account for experimental results involving magnetic grains [72].

Since the IMM represents a mathematical model for granular gases much more tractable than the conventional IHSM, thus allowing for exact results, it is then natural to construct a Maxwell version of the more realistic IRHSM. In this sense, we recently proposed the inelastic rough Maxwell model (IRMM) [73] and derived the exact collisional production rates of first and second degree, as well as the most relevant ones of third and fourth degree. In addition, the results were applied to the evaluation of the rotational-to-translational temperature ratio and the velocity cumulants in the HCS.

The aim of this paper is to derive the exact expressions for the NSF transport coefficients predicted by the IRMM, namely the shear and bulk viscosities, the thermal and diffusive heat conductivities, and the cooling-rate transport coefficient, that is, the same coefficients as previously obtained from the IRHSM by the application of a standard Sonine approximation [17, 19]. However, in contrast to the approximate treatment in the IRHSM, the study carried out here unveils an inherent coupling of the shear viscosity to a spin viscosity (associated with the spin-spin tensor defined below). Likewise, the heat conductivities are coupled to their counterparts associated with a torque-vorticity vector defined below. Thus, the exact results derived here can shed light on improvements of the Sonine approximation as applied to the IRHSM.

Our study is motivated by the challenge of deriving exact NSF transport coefficients for granular gases, focusing on the IRMM. Unlike molecular gases, granular gases present unknown nonequilibrium base states, complicating the derivation process. Existing models rely on nontrivial approximations, but the IRMM avoids such complexities and uncovers previously overlooked couplings, thus enhancing our understanding of granular gas dynamics. This contributes to refining approximation methods in similar models, serving as a valuable benchmark in both theory and practical applications.

This paper is organized as follows. The IRMM is presented in Sect. 2, where also the necessary collisional production rates are provided, the expressions of the associated coefficients being found in the Appendix. Next, in Sect. 3, the Chapman–Enskog method is applied to the explicit evaluation of the NSF transport coefficients. The results are extensively discussed in Sect. 4, including a stability analysis of the HCS and an assessment of the impact of the couplings mentioned above. Finally, the paper is closed in Sect. 5 with some concluding remarks.

2 The Inelastic Rough Maxwell Model. Exact Collisional Production Rates

Let us consider a granular gas made of spherical particles of diameter $\sigma$, mass m, and moment of inertia I. The translational and angular velocities of a particle will be denoted by $\textbf{v}$ and $\varvec{\omega }$, respectively. According to the IRMM, the Boltzmann equation for the gas reads [73]

$$\begin{aligned} \frac{\partial f}{\partial t}+\textbf{v}\cdot \nabla f=J[\varvec{\Gamma }\vert f,f], \end{aligned}$$

(1)

where $\varvec{\Gamma }\equiv \{\textbf{v},\varvec{\omega }\}$, $d\varvec{\Gamma }\equiv d\textbf{v}d\varvec{\omega }$, $f(\textbf{r},\varvec{\Gamma }, t)$ is the one-particle velocity distribution function (VDF), and $J[\varvec{\Gamma }\vert f,f]$ is the bilinear collision operator

$$\begin{aligned} J[\varvec{\Gamma }_1\vert f,f]=\frac{\nu }{4\pi n}\int d\varvec{\Gamma }_2\int d\widehat{\varvec{\sigma }}\,\left( \frac{{\hat{b}}_{12}^{-1}}{\alpha \beta ^2} -1\right) f(\varvec{\Gamma }_1)f(\varvec{\Gamma }_2). \end{aligned}$$

(2)

Here, $\nu$ is an effective mean-field collision frequency, $n=\int d\varvec{\Gamma }\, f(\varvec{\Gamma })$ is the number density, and $\widehat{\varvec{\sigma }}$ is the intercenter unit vector at contact. The operator ${\hat{b}}_{12}$ expresses the postcollision velocities as functions of the precollisional velocities, the unit vector $\widehat{\varvec{\sigma }}$, the coefficient of normal restitution $\alpha$, the coefficient of tangential restitution $\beta$, and the reduced moment of inertia $\kappa \equiv 4I/m\sigma ^2$ [3, 17, 19]. For a brief description of the collision rules, the reader is referred to Sect. 2.1 of Ref. [73].

The range of the coefficient of normal restitution is $0\le \alpha \le 1$, where $\alpha =1$ corresponds to a perfectly elastic collision (in the sense that the normal component of the relative velocity just changes its sign upon collision). In the case of the coefficient of tangential restitution, the range is $-1\le \beta \le 1$, where $\beta =-1$ defines smooth particles (the tangential component of the relative velocity is preserved after a collision) and $\beta =1$ corresponds to perfectly rough particles (the tangential component of the relative velocity just changes its sign upon collision). Only if $\alpha =\mid \!\!\beta \!\!\mid =1$ is the kinetic energy conserved by collisions. As for the reduced moment of inertia, $\kappa =0$ refers to a concentration of the mass at the center of the sphere, $\kappa =\frac{2}{5}$ means a uniform mass distribution, and $\kappa =\frac{2}{3}$ if the mass is concentrated on the surface.

Given an arbitrary function $\Psi (\varvec{\Gamma })$, we define its average value as

$$\begin{aligned} \langle \Psi \rangle =\frac{1}{n}\int d\varvec{\Gamma }\, \Psi (\varvec{\Gamma })f(\varvec{\Gamma }). \end{aligned}$$

(3)

In particular,

$$\begin{aligned} \textbf{u}= & {} \left\langle \textbf{v}\right\rangle ,\quad \varvec{\Omega }=\left\langle \varvec{\omega }\right\rangle , \end{aligned}$$

(4a)

$$\begin{aligned} T_t= & {} \frac{m}{3}\left\langle V^2\right\rangle ,\quad T_r=\frac{I}{3}\left\langle \omega ^2\right\rangle ,\quad T=\frac{1}{2}\left( T_t+T_r\right) , \end{aligned}$$

(4b)

$$\begin{aligned} \Upsilon _{ij}= & {} n\left\langle \omega _iV_j\right\rangle ,\quad \Pi _{ij}=n m\left\langle V_iV_j-\frac{1}{3}V^2\delta _{ij}\right\rangle , \end{aligned}$$

(4c)

$$\begin{aligned} \Omega _{ij}= & {} nI\left\langle \omega _i\omega _j-\frac{1}{3}\omega ^2\delta _{ij}\right\rangle , \quad \textbf{Q}=\frac{nI}{2}\left\langle (\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\right\rangle , \end{aligned}$$

(4d)

$$\begin{aligned} \textbf{q}_t= & {} \frac{nm}{2}\left\langle V^2\textbf{V}\right\rangle ,\quad \textbf{q}_r=\frac{nI}{2}\left\langle \omega ^2\textbf{V}\right\rangle ,\quad \textbf{q}=\textbf{q}_t+\textbf{q}_r. \end{aligned}$$

(4e)

Here, $\textbf{u}$ is the flow velocity, $\textbf{V}\equiv \textbf{v}-\textbf{u}$ is the translational peculiar velocity, $T_t$ is the translational granular temperature, $\Pi _{ij}$ is the (traceless) stress tensor, and $\textbf{q}_t$ is the translational heat flux vector. Note that the pressure tensor is

$$\begin{aligned} P_{ij}=nm\left\langle V_iV_j\right\rangle =\Pi _{ij}+nT_t\delta _{ij}. \end{aligned}$$

(5)

The quantities $T_t$, $\Pi _{ij}$, $P_{ij}$, and $\textbf{q}_t$ are defined in the same way as in the case of smooth particles. However, the rotational degrees of freedom make it necessary to define additional quantities like the mean spin vector $\varvec{\Omega }$, the rotational granular temperature $T_r$, the (traceless) spin-spin tensor $\Omega _{ij}$, and the rotational heat flux vector $\textbf{q}_r$. The mean granular temperature and the total heat flux are T and $\textbf{q}$, respectively. Finally, the asymmetric quantity $\Upsilon _{ij}$ represents the couple stress tensor [74, 75], while $\textbf{Q}$ can be seen as a torque-vorticity vector.

The collisional production rate of an arbitrary function $\Psi (\varvec{\Gamma })$ is

$$\begin{aligned} {\mathcal {J}}[\Psi ]\equiv&\frac{1}{n}\int d\varvec{\Gamma }\,\Psi (\varvec{\Gamma })J[\varvec{\Gamma }\vert f,f]\nonumber \\ =&\frac{\nu }{8\pi n^2}\int d\varvec{\Gamma }_1\int d\varvec{\Gamma }_2\int d\widehat{\varvec{\sigma }}\,f(\varvec{\Gamma }_1)f(\varvec{\Gamma }_2) \left( {\hat{b}}_{12}-1\right) \left[ \Psi (\varvec{\Gamma }_1)+\Psi (\varvec{\Gamma }_2)\right] . \end{aligned}$$

(6)

In Ref. [73], it was proved that the collisional production rates associated with the quantities defined in Eqs. (4) are ${\mathcal {J}}\left[ \textbf{v}\right] =0$ and

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ \varvec{\omega }\right] }= \varphi _{01\mid 01} \varvec{\Omega },\quad -\nu ^{-1}{{\mathcal {J}}\left[ n\omega _i V_j\right] }= \psi _{11\mid 11} \Upsilon _{ij}, \end{aligned}$$

(7a)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ \frac{m}{3}V^2\right] }=\chi _{20\mid 20}T_t +\frac{4}{\kappa }\chi _{20\mid 02}T_r, \end{aligned}$$

(7b)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ \frac{I}{3}\omega ^2\right] }= \frac{\kappa }{4}\chi _{02\mid 20}T_t+\chi _{02\mid 02}T_r, \end{aligned}$$

(7c)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ nm\left( V_iV_j-\frac{1}{3}V^2\delta _{ij}\right) \right] }= \psi _{20\mid 20}\Pi _{ij} +\frac{4}{\kappa }\psi _{20\mid 02}\Omega _{ij}, \end{aligned}$$

(7d)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ nI\left( \omega _i\omega _j-\frac{1}{3}\omega ^2\delta _{ij}\right) \right] }= \frac{\kappa }{4}\psi _{02\mid 20}\Pi _{ij}+\psi _{02\mid 02}\Omega _{ij}, \end{aligned}$$

(7e)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ n\frac{m}{2}V^2\textbf{V}\right] }=\varphi _{30\mid 30}\textbf{q}_t+\frac{4}{\kappa }\varphi _{30\mid 12}\left( 2\textbf{q}_r-\textbf{Q}\right) , \end{aligned}$$

(7f)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ n\frac{I}{2}\omega ^2\textbf{V}\right] }=\varphi _{12\mid 12}^{(1)}\textbf{q}_r+\varphi _{12\mid 12}^{(2)}\textbf{Q}+\frac{\kappa }{4}\varphi _{12\mid 30}\textbf{q}_t, \end{aligned}$$

(7g)

$$\begin{aligned}{} & {} -\nu ^{-1}{{\mathcal {J}}\left[ n\frac{I}{2}(\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\right] }=\overline{\varphi }_{12\mid 12}^{(1)}\textbf{Q}+ \overline{\varphi }_{12\mid 12}^{(2)}\textbf{q}_r. \end{aligned}$$

(7h)

It should be noted that in Eqs. (7b)–(7h) it has been assumed that $\varvec{\Omega }=\Upsilon _{ij}=0$ (see below). The explicit expressions for the 17 coefficients in Eqs. (7) can be found in the Appendix.

The cooling rate $\zeta$ is defined as

$$\begin{aligned} \int d\varvec{\Gamma }\left( mV^2+I\omega ^2\right) J[\varvec{\Gamma }\vert f,f]=-6\zeta nT. \end{aligned}$$

(8)

Thus,

$$\begin{aligned} \zeta =\frac{\nu }{2T}\left[ \left( \chi _{20\mid 20}+\frac{\kappa }{4}\chi _{02\mid 20}\right) T_t+\left( \frac{4}{\kappa }\chi _{20\mid 02}+\chi _{02\mid 02}\right) T_r\right] . \end{aligned}$$

(9)

3 Chapman–Enskog Expansion

3.1 General Scheme

The hydrodynamic balance equations of a granular gas are [3, 17, 19]

$$\begin{aligned}{} & {} {{{\mathcal {D}}}}_tn+n\nabla \cdot \textbf{u}=0, \end{aligned}$$

(10a)

$$\begin{aligned}{} & {} m n{{{\mathcal {D}}}}_t \textbf{u}+\nabla \cdot \textsf{P} =0, \end{aligned}$$

(10b)

$$\begin{aligned}{} & {} {{{\mathcal {D}}}}_t T +{1\over 3n}\left( \nabla \cdot \textbf{q}+\textsf{P}:\nabla \textbf{u}\right) + T \zeta =0. \end{aligned}$$

(10c)

In the above balance equations, ${\mathcal {D}}_t=\partial _t+\textbf{u}\cdot \nabla$ denotes the material time derivative. Equations (10) are formally exact but they do not make a closed set unless constitutive equations expressing the pressure tensor $\textsf{P}$, the heat flux $\textbf{q}$, and the cooling rate $\zeta$ in terms of gradients of the hydrodynamic fields (n, $\textbf{u}$, and T) are proposed. The form of the NSF constitutive equations (to first order in gradients) is

$$\begin{aligned} P_{ij}= & {} n T_t^{(0)}\delta _{ij}-\eta \Delta _{ijk\ell }\nabla _k u_\ell -\eta _b\nabla \cdot \textbf{u}, \end{aligned}$$

(11a)

$$\begin{aligned} \textbf{q}= & {} -\lambda \nabla T-\mu \nabla n, \end{aligned}$$

(11b)

$$\begin{aligned} \zeta= & {} \zeta ^{(0)}-\xi \nabla \cdot \textbf{u}, \end{aligned}$$

(11c)

where the superscript (0) denotes quantities in the absence of gradients, $\eta$ is the shear viscosity, $\eta _b$ is the bulk viscosity, $\lambda$ is the thermal heat conductivity, $\mu$ is a Dufour-like coefficient that will be referred to as the diffusive heat conductivity [3, 76, 77], and $\xi$ is a cooling-rate transport coefficient. In Eq. (11a),

$$\begin{aligned} \Delta _{ijk\ell }\equiv \delta _{ik}\delta _{j\ell }+\delta _{i\ell }\delta _{jk}-\frac{2}{3}\delta _{ij}\delta _{k\ell }. \end{aligned}$$

(12)

Although the spin-spin tensor $\varvec{\Omega }$, the partial heat fluxes $\textbf{q}_t$ and $\textbf{q}_r$, and the torque-vorticity vector $\textbf{Q}$ do not enter into the hydrodynamic balance equations, Eqs. (10), they are intrinsically coupled to $\textbf{P}$ and $\textbf{q}$, as will be seen below.

Moreover, multiplying both sides of Eq. (1) by $\omega _i$ and integrating over $\varvec{\Gamma }$, one gets the following balance equation:

$$\begin{aligned} n{{{\mathcal {D}}}}_t\Omega _i+\nabla _j\Upsilon _{ij}=-n\nu \varphi _{01\mid 01}\Omega _i, \end{aligned}$$

(13)

where use has been made of Eq. (7a). This shows that, in a homogeneous state, $\partial _t\varvec{\Omega }=-\nu \varphi _{01\mid 01}\varvec{\Omega }$, so that $\lim _{t\rightarrow \infty }\varvec{\Omega }(t)=0$. As we will see later, $\varvec{\Omega }\rightarrow 0$ even in the presence of hydrodynamic gradients to NSF order.

Our main goal is to derive the exact expressions for the NSF transport coefficients within the IRMM. To that end, we assume that the VDF depends on space and time only through the slow hydrodynamic fields (n, $\textbf{u}$, and T), and apply the Chapman–Enskog expansion method [3, 17, 19].

In the method, one introduces a bookkeeping parameter $\epsilon$, which is used as a small parameter, so that

$$\begin{aligned}{} & {} \nabla \rightarrow \epsilon \nabla ,\quad f = f^{(0)} + \epsilon f^{(1)} + \epsilon ^2 f^{(2)}+\cdots , \end{aligned}$$

(14a)

$$\begin{aligned}{} & {} {\mathcal {D}}_t={\mathcal {D}}_t^{(0)}+\epsilon {\mathcal {D}}_t^{(1)}+\epsilon ^2{\mathcal {D}}_t^{(2)}+\cdots . \end{aligned}$$

(14b)

Thus, the Boltzmann equation, Eq. (1), decouples into a hierarchy of equations of orders $k=0,1,2,\ldots$. To make a qualitative contact with the case of hard spheres, we will assume $\nu \propto n\sqrt{T_t}$ and, therefore, the collision frequency must also be expanded as

$$\begin{aligned} \nu =\nu ^{(0)}+\epsilon \nu ^{(1)}+\epsilon ^2 \nu ^{(2)}+\cdots , \end{aligned}$$

(15)

where

$$\begin{aligned} \nu ^{(0)}\propto n\sqrt{T_t^{(0)}},\quad \nu ^{(1)}=\nu ^{(0)}\frac{T_t^{(1)}}{2T_t^{(0)}}. \end{aligned}$$

(16)

3.2 Zeroth-Order Distribution

The zeroth-order Boltzmann equation is

$$\begin{aligned} {\mathcal {D}}_t^{(0)} f^{(0)}(\varvec{\Gamma }_1) =\frac{\nu ^{(0)}}{4\pi n}\int d\varvec{\Gamma }_2\int d\widehat{\varvec{\sigma }}\,\left( \frac{{\hat{b}}_{12}^{-1}}{\alpha \beta ^2} -1\right) f^{(0)}(\varvec{\Gamma }_1) f^{(0)}(\varvec{\Gamma }_2). \end{aligned}$$

(17)

This shows that the zeroth-order VDF $f^{(0)}$ is the local version of the HCS VDF. Thus, according to the results derived in Ref. [73],

$$\begin{aligned} T_t^{(0)}=\tau _t T,\quad T_r^{(0)}=\tau _r T, \end{aligned}$$

(18)

where

$$\begin{aligned} \tau _t=2-\tau _r=\frac{2}{1+\theta },\quad \theta =h+\sqrt{1+h^2}, \end{aligned}$$

(19a)

$$\begin{aligned} h\equiv&\frac{\kappa }{8}\frac{\chi _{02\mid 02}-\chi _{20\mid 20}}{\chi _{20\mid 02}}\nonumber \\ =&\frac{1+\kappa }{2\kappa (1+\beta )}\left[ \frac{1+\kappa }{2}\frac{1-\alpha ^2}{1+\beta }-(1-\kappa )(1-\beta )\right] . \end{aligned}$$

(19b)

Moreover, the zeroth-order cooling rate is

$$\begin{aligned} \zeta ^{(0)}=&\frac{\nu ^{(0)}}{2}\left[ \left( \chi _{20\mid 20}+\frac{\kappa }{4}\chi _{02\mid 20}\right) \tau _t +\left( \frac{4}{\kappa }\chi _{20\mid 02}+\chi _{02\mid 02}\right) \tau _r\right] \nonumber \\ =&\frac{1}{6}\frac{\nu ^{(0)}}{1+\theta }\left[ 1-\alpha ^2+2\frac{1-\beta ^2}{1+\kappa }\theta \left( \frac{\kappa }{\theta }+1\right) \right] . \end{aligned}$$

(20)

Note also that, in the HCS,

$$\begin{aligned} \varvec{\Omega }^{(0)}=\Upsilon _{ij}^{(0)}=\Pi _{ij}^{(0)}=\Omega _{ij}^{(0)}=\textbf{q}_t^{(0)}=\textbf{q}_r^{(0)}=\textbf{Q}^{(0)}=0. \end{aligned}$$

(21)

3.3 First-Order Distribution

The integral equation for the first-order distribution function $f^{(1)}$ is

$$\begin{aligned} \left( {{{\mathcal {D}}}}_t^{(0)} +{\mathcal {L}}\right) f^{(1)}=&-\left( {{{\mathcal {D}}}}_t^{(1)}+\textbf{V}\cdot \nabla -\frac{T_t^{(1)}}{2\tau _tT}{{{\mathcal {D}}}}_t^{(0)}\right) f^{(0)}\nonumber \\ =&\textbf{A}\cdot \nabla \ln T\!+\!\textbf{B}\cdot \nabla \ln n{+C_{ij}\nabla _j u_i} \!+\!E\nabla \cdot \textbf{u}\!+\!\left[ \!\zeta ^{(1)}\!-\!\frac{\nu ^{(1)}\zeta ^{(0)}}{\nu ^{(0)}}\!\right] \! T \partial _T f^{(0)}, \end{aligned}$$

(22)

where

$$\begin{aligned} {\mathcal {L}}f^{(1)}(\varvec{\Gamma }_1)=&-\frac{\nu ^{(0)}}{4\pi n}\int d\varvec{\Gamma }_2\int d\widehat{\varvec{\sigma }}\,\left( \frac{{\hat{b}}_{12}^{-1}}{\alpha \beta ^2} -1\right) \left[ f^{(0)}(\varvec{\Gamma }_1)f^{(1)}(\varvec{\Gamma }_2) +f^{(0)}(\varvec{\Gamma }_2)f^{(1)}(\varvec{\Gamma }_1)\right] \end{aligned}$$

(23)

is the linearized collision operator and

$$\begin{aligned}{} & {} \textbf{A}=-T\left( \textbf{V}\partial _T+\frac{\tau _t}{m}\partial _{\textbf{V}}\right) f^{(0)}, \end{aligned}$$

(24a)

$$\begin{aligned}{} & {} \textbf{B}=-\left( \textbf{V}+\frac{\tau _t T}{m}\partial _{\textbf{V}}\right) f^{(0)}, \end{aligned}$$

(24b)

$$\begin{aligned}{} & {} {C_{ij}=\left( \partial _{V_i}V_j-\frac{1}{3}\delta _{ij}\partial _{\textbf{V}}\cdot \textbf{V}\right) f^{(0)}}, \end{aligned}$$

(24c)

$$\begin{aligned}{} & {} E=\frac{1}{3}\left( \partial _\textbf{V}\cdot \textbf{V}+\tau _tT\partial _T\right) f^{(0)}, \end{aligned}$$

(24d)

$$\begin{aligned}{} & {} T\partial _T f^{(0)}=-\frac{1}{2}\left( \partial _{\textbf{V}}\cdot \textbf{V}+\partial _{\varvec{\omega }}\cdot \varvec{\omega }\right) f^{(0)}. \end{aligned}$$

(24e)

The solution to Eq. (22) has the form

$$\begin{aligned} f^{(1)}= \varvec{{\mathcal {A}}}\cdot \nabla \ln T+\varvec{{\mathcal {B}}}\cdot \nabla \ln n+{\mathcal {C}}_{ij}{\nabla _j u_i} +{\mathcal {E}}\nabla \cdot \textbf{u}, \end{aligned}$$

(25)

where the vectors $\varvec{{\mathcal {A}}}$ and $\varvec{{\mathcal {B}}}$, the traceless tensor ${\mathcal {C}}_{ij}$, and the scalar ${\mathcal {E}}$ are the solutions of a set of linear integral equations [17, 19]. Whereas in the case of the IRHSM one needs to solve those integral equations by a Sonine approximation, the main advantage of the IRMM is that it is not necessary to get any explicit expression for $f^{(1)}$ (or, equivalently, $\varvec{{\mathcal {A}}}$, $\varvec{{\mathcal {B}}}$, ${\mathcal {C}}_{ij}$ and ${\mathcal {E}}$) since, as shown below, one can exactly obtain the transport coefficients by just taking velocity moments in both sides of Eq. (22).

From Eq. (23), it follows that

$$\begin{aligned} \int d\varvec{\Gamma }\,\Psi (\varvec{\Gamma }){\mathcal {L}}f^{(1)}(\varvec{\Gamma })=-n\nu ^{(0)}\nu ^{-1}{\mathcal {J}}[\Psi ] \end{aligned}$$

(26)

for the functions $\Psi (\varvec{\Gamma })$ in Eqs. (7a) and (7d)–(7h). On the other hand, from Eqs. (7b) and (7c), we have

$$\begin{aligned} \int d\varvec{\Gamma }\begin{Bmatrix} \frac{m}{3}V^2\\ \frac{I}{3}\omega ^2 \end{Bmatrix} {\mathcal {L}}f^{(1)}(\varvec{\Gamma })=n\nu ^{(0)} \begin{Bmatrix} \chi _{20\mid 20} -\frac{4}{\kappa }\chi _{20\mid 02}\\ \frac{\kappa }{4}\chi _{02\mid 20}-\chi _{02\mid 02} \end{Bmatrix}{T_t^{(1)}}, \end{aligned}$$

(27)

where we have taken into account that $T_t^{(1)}+T_r^{(1)}=0$ since T is a hydrodynamic variable. From Eq. (9), we obtain the following expression for the first-order cooling rate,

$$\begin{aligned} \zeta ^{(1)}=\nu ^{(0)}\left[ \frac{\zeta ^{(0)}}{\tau _t\nu ^{(0)}}+\chi _{20\mid 20} -\frac{4\chi _{20\mid 02}}{\kappa }-\chi _{02\mid 02} +\frac{\kappa \chi _{02\mid 20}}{4}\right] \frac{T_t^{(1)}}{2T}, \end{aligned}$$

(28)

where we have taken into account that $\nu ^{(1)}/\nu ^{(0)}={T_t^{(1)}}/{2\tau _t T}$ [see Eq. (16)].

As we will show below, Eqs. (26) and (27), together with the collision integrals in Eqs. (7), allow us to derive exactly the NSF transport coefficients within the IRMM.

3.3.1 Mean Spin Vector and Couple Stress Tensor

Here we prove that $\varvec{\Omega }^{(1)}=\Upsilon _{ij}^{(1)}=0$. Suppose that $\varvec{\Omega }^{(1)}$ and $\Upsilon _{ij}^{(1)}$ existed; then, by symmetry, one should have

$$\begin{aligned} n\varvec{\Omega }^{(1)}= & {} \varsigma _1\nabla T+\varsigma _2\nabla n, \end{aligned}$$

(29a)

$$\begin{aligned} \Upsilon _{ij}^{(1)}= & {} \varsigma _{ijk\ell }\nabla _ku_\ell , \end{aligned}$$

(29b)

where $\varsigma _1$, $\varsigma _2$, and $\varsigma _{ijk\ell }$ would be the associated transport coefficients. By dimensional analysis, $\varsigma _1\propto n/\sqrt{mT}$, $\varsigma _2\propto \sqrt{T/m}$, $\varsigma _{ijk\ell }\propto n\sqrt{T/m}$. Therefore,

$$\begin{aligned} {\mathcal {D}}_t^{(0)}n\varvec{\Omega }^{(1)}= & {} \zeta ^{(0)}\left[ -\varsigma _1\nabla T+\left( \frac{\varsigma _1 T}{n}-\varsigma _2\right) \nabla n\right] , \end{aligned}$$

(30a)

$$\begin{aligned} {\mathcal {D}}_t^{(0)}\Upsilon _{ij}^{(1)}= & {} -\frac{\zeta ^{(0)}}{2}\varsigma _{ijk\ell }\nabla _ku_\ell . \end{aligned}$$

(30b)

In Eq. (30a), use has been made of the property

$$\begin{aligned} {\mathcal {D}}_t^{(0)}\nabla T=-\nabla \left[ \zeta ^{(0)}T\right] =-\zeta ^{(0)}\left( \frac{T}{n}\nabla n+\frac{3}{2}\nabla T\right) . \end{aligned}$$

(31)

The key point is that, as can be easily checked from Eqs. (24),

$$\begin{aligned} \int d\varvec{\Gamma }\begin{Bmatrix} \varvec{\omega }\\ \omega _iV_j \end{Bmatrix} A_k(\varvec{\Gamma })=&\int d\varvec{\Gamma }\begin{Bmatrix} \varvec{\omega }\\ \omega _iV_j \end{Bmatrix} B_k(\varvec{\Gamma }) =\int d\varvec{\Gamma }\begin{Bmatrix} \varvec{\omega }\\ \omega _iV_j \end{Bmatrix} C_{k\ell }(\varvec{\Gamma })\nonumber \\ =&\int d\varvec{\Gamma }\begin{Bmatrix} \varvec{\omega }\\ \omega _iV_j \end{Bmatrix} E(\varvec{\Gamma }) =\int d\varvec{\Gamma }\begin{Bmatrix} \varvec{\omega }\\ \omega _iV_j \end{Bmatrix} \partial _T f^{(0)}(\varvec{\Gamma })=0. \end{aligned}$$

(32)

As a consequence, multiplying both sides of Eq. (22) by $\varvec{\omega }$ or $\omega _iV_j$, integrating over $\varvec{\Gamma }$, and using Eqs. (7a), (26), and (30), one gets

$$\begin{aligned} \varsigma _1\zeta ^{(0)}= & {} \varsigma _1\nu ^{(0)} \varphi _{01\mid 01}, \end{aligned}$$

(33a)

$$\begin{aligned} \left( \varsigma _2-\frac{\varsigma _1 T}{n}\right) \zeta ^{(0)}= & {} \varsigma _2\nu ^{(0)} \varphi _{01\mid 01}, \end{aligned}$$

(33b)

$$\begin{aligned} \varsigma _{ijk\ell }\frac{\zeta ^{(0)}}{2}= & {} \varsigma _{ijk\ell }\nu ^{(0)} \psi _{11\mid 11}. \end{aligned}$$

(33c)

The solution of those homogeneous linear equations is the trivial one, i.e.,

$$\begin{aligned} \varsigma _1=\varsigma _2=\varsigma _{ijk\ell }=0. \end{aligned}$$

(34)

Therefore, $\varvec{\Omega }^{(1)}=\Upsilon _{ij}^{(1)}=0$, thus justifying the assumption made in Eqs. (7b)–(7h).

The observation that both the mean spin vector and the couple stress tensor vanish to NSF order in the IRMM contrasts with expectations based on the micropolar fluids framework [74, 75, 78, 79]. This contrast is likely attributed to the low-density regime where the Boltzmann equation is applicable and the absence of boundary effects in the NSF description [75]. However, in the case of dense gases, a dependency does emerge, with the stress tensor correlating with the anti-symmetric part of flow velocity and mean spin gradients, and the heat flux linking with the curl of the mean spin vector [80,81,82,83].

3.3.2 Stress and Spin–Spin Tensors

Now we turn to the first-order translational temperature $T_t^{(1)}$, stress tensor $\Pi _{ij}^{(1)}$, and spin-spin tensor $\Omega _{ij}^{(1)}$. First, we need to make use of the integrals

$$\begin{aligned}{} & {} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} m V_iV_j\\ I \omega _i\omega _j \end{array} \right\} \textbf{A}(\varvec{\Gamma })=\int d\varvec{\Gamma }\, \left\{ \begin{array}{c} mV_iV_j\\ I\omega _i\omega _j \end{array} \right\} \textbf{B}(\varvec{\Gamma })=0, \end{aligned}$$

(35a)

$$\begin{aligned}{} & {} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} mV_iV_j\\ I\omega _i\omega _j \end{array} \right\} {E}(\varvec{\Gamma })=\frac{n\tau _t\tau _rT}{3}\left\{ \begin{array}{c} -1\\ 1 \end{array} \right\} \delta _{ij}, \end{aligned}$$

(35b)

$$\begin{aligned}{} & {} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} mV_iV_j\\ I\omega _i\omega _j \end{array} \right\} T\partial _T f^{(0)}(\varvec{\Gamma })={nT}\left( \begin{array}{c} \tau _t\\ \tau _r \end{array} \right) \delta _{ij}, \end{aligned}$$

(35c)

$$\begin{aligned}{} & {} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} mV_iV_j\\ I\omega _i\omega _j \end{array} \right\} C_{k\ell }(\varvec{\Gamma })=-{n\tau _tT}\left\{ \begin{array}{c} \Delta _{ijk\ell }\\ 0 \end{array} \right\} . \end{aligned}$$

(35d)

Multiplying both sides of Eq. (22) by $mV^2/3$ and integrating over $\varvec{\Gamma }$, one gets

$$\begin{aligned} {\mathcal {D}}_t^{(0)} T_t^{(1)} \!+\!\nu ^{(0)} \left( \!\chi _{20\mid 20} \!-\!\frac{4}{\kappa }\chi _{20\mid 02}\!\right) \!{T_t^{(1)}} \!=&\!-\frac{\tau _t \tau _r}{3} T\nabla \cdot \textbf{u}\!+\!\left[ \zeta ^{(1)}\!-\!\frac{\nu ^{(1)}\zeta ^{(0)}}{\nu ^{(0)}}\right] \tau _t T. \end{aligned}$$

(36)

The same result is obtained by multiplication of both sides of Eq. (22) by $I\omega ^2/3$.

Since, on physical grounds, $T_t^{(1)}$ must be proportional to $T(\nabla \cdot \textbf{u})/\nu ^{(0)}$, it turns out that ${\mathcal {D}}_t^{(0)}T_t^{(1)}=\frac{1}{2}[T_t^{(1)}/T]{\mathcal {D}}_t^{(0)}T=-\frac{1}{2}\zeta ^{(0)}T_t^{(1)}$. Thus, $T_t^{(1)}$ is given by

$$\begin{aligned} {T_t^{(1)}}=-\frac{\eta _b}{n} \nabla \cdot \textbf{u}, \end{aligned}$$

(37)

where

$$\begin{aligned} \eta _b=\frac{2nT\tau _t\tau _r}{3\nu ^{(0)}}\left[ \left( \chi _{20\mid 20} -\frac{4\chi _{20\mid 02}}{\kappa }\right) \tau _r +\left( \chi _{02\mid 02} -\frac{\kappa \chi _{02\mid 20}}{4}\right) \tau _t-\frac{\zeta ^{(0)}}{\nu ^{(0)}}\right] ^{-1} \end{aligned}$$

(38)

is the bulk viscosity introduced in Eq. (11a). Combination of Eqs. (28), (37), and (38) yields the cooling-rate transport coefficient introduced in Eq. (11c), namely

$$\begin{aligned} \xi =\left[ \frac{\zeta ^{(0)}}{\tau _t\nu ^{(0)}}+\chi _{20\mid 20} -\frac{4\chi _{20\mid 02}}{\kappa }-\chi _{02\mid 02} +\frac{\kappa \chi _{02\mid 20}}{4}\right] \frac{\eta _b}{2nT/\nu ^{(0)}}. \end{aligned}$$

(39)

Next, from Eq. (22) one obtains

$$\begin{aligned} {\mathcal {D}}_t^{(0)} \begin{Bmatrix} \Pi _{ij}^{(1)}\\ \Omega _{ij}^{(1)} \end{Bmatrix} +\nu ^{(0)} \begin{Bmatrix} \psi _{20\mid 20}\Pi _{ij}^{(1)} +\frac{4}{\kappa }\psi _{20\mid 02}\Omega _{ij}^{(1)}\\ \psi _{02\mid 02}\Omega _{ij}^{(1)} +\frac{\kappa }{4}\psi _{02\mid 20}\Pi _{ij}^{(1)} \end{Bmatrix} =\nabla _k u_\ell \begin{Bmatrix} -n\tau _t T\Delta _{ijk\ell }\\ 0 \end{Bmatrix}. \end{aligned}$$

(40)

Both $\Pi _{ij}^{(1)}$ and $\Omega _{ij}^{(1)}$ must be proportional to $n T\Delta _{ijk\ell }\nabla _k u_\ell$, so that ${\mathcal {D}}_t^{(0)}\Pi _{ij}^{(1)}=-\frac{1}{2}\zeta ^{(0)}\Pi _{ij}^{(1)}$ and ${\mathcal {D}}_t^{(0)}\Omega _{ij}^{(1)}=-\frac{1}{2}\zeta ^{(0)}\Omega _{ij}^{(1)}$. We then write

$$\begin{aligned} \Pi _{ij}^{(1)}=-\eta \Delta _{ijk\ell }\nabla _k u_\ell ,\quad \Omega _{ij}^{(1)}=-\eta _\Omega \Delta _{ijk\ell }\nabla _k u_\ell , \end{aligned}$$

(41)

where $\eta$ is the shear viscosity introduced in Eq. (11a) and $\eta _\Omega$ is a new transport coefficient here called spin viscosity. Inserting this into Eq. (40), we obtain the set of linear equations

$$\begin{aligned} \begin{bmatrix} \nu ^{(0)}\psi _{20\mid 20}-\frac{1}{2}\zeta ^{(0)}&{}\nu ^{(0)}\frac{4\psi _{20\mid 02}}{\kappa }\\ \nu ^{(0)}\frac{\kappa \psi _{02\mid 20}}{4}&{}\nu ^{(0)}\psi _{02\mid 02}-\frac{1}{2}\zeta ^{(0)} \end{bmatrix} \cdot \begin{bmatrix} \eta \\ \eta _\Omega \end{bmatrix} =\begin{bmatrix} n\tau _t T\\ 0 \end{bmatrix}, \end{aligned}$$

(42)

whose solution is

$$\begin{aligned} \eta _\Omega= & {} -\frac{\kappa }{4}\frac{\psi _{02\mid 20}}{\psi _{02\mid 02}-\frac{1}{2}\zeta ^{(0)}/\nu ^{(0)}}\eta , \end{aligned}$$

(43a)

$$\begin{aligned} \eta= & {} \frac{n\tau _t T}{\nu ^{(0)}}\left[ \psi _{20\mid 20}-\frac{1}{2}\frac{\zeta ^{(0)}}{\nu ^{(0)}}-\frac{\psi _{20\mid 02}\psi _{02\mid 20}}{\psi _{02\mid 02}-\frac{1}{2}\zeta ^{(0)}/\nu ^{(0)}}\right] ^{-1}. \end{aligned}$$

(43b)

3.3.3 Heat Flux and Torque-Vorticity Vectors

Let us now consider the first-order translational ($\textbf{q}_t$) and rotational ($\textbf{q}_r$) contributions to the heat flux ($\textbf{q}=\textbf{q}_t+\textbf{q}_r$), as well as the torque-vorticity vector $\textbf{Q}$.

The needed integrals are

$$\begin{aligned}{} & {} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} \frac{m}{2} V^2{V_i}\\ \frac{I}{2} \omega ^2{V_i}\\ \frac{I}{2} (\textbf{V}\cdot \varvec{\omega })\omega _i \end{array} \right\} {A_j}(\varvec{\Gamma })=-\frac{nT^2\tau _t}{2m}\left\{ \begin{array}{c} 5\tau _t\left( 1+2a_{20}\right) \\ 3\tau _r\left( 1+2a_{11}\right) \\ \tau _r\left( 1+2b_{00}\right) \end{array} \right\} \delta _{ij}, \end{aligned}$$

(44a)

$$\begin{aligned}{} & {} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} \frac{m}{2} V^2{V_i}\\ \frac{I}{2} \omega ^2{V_i}\\ \frac{I}{2} (\textbf{V}\cdot \varvec{\omega })\omega _i \end{array} \right\} {B_j}(\varvec{\Gamma })=-\frac{nT^2\tau _t}{2m}\left\{ \begin{array}{c} 5\tau _t a_{20}\\ 3\tau _r a_{11}\\ \tau _r b_{00} \end{array} \right\} \delta _{ij}, \end{aligned}$$

(44b)

$$\begin{aligned} \int d\varvec{\Gamma }\, \left\{ \begin{array}{c} \frac{m}{2} V^2\textbf{V}\\ \frac{I}{2} \omega ^2\textbf{V}\\ \frac{I}{2} (\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\end{array} \right\} {E}(\varvec{\Gamma }) =&\int d\varvec{\Gamma }\, \left\{ \begin{array}{c} \frac{m}{2} V^2\textbf{V}\\ \frac{I}{2} \omega ^2\textbf{V}\\ \frac{I}{2} (\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\end{array} \right\} T\partial _T f^{(0)}(\varvec{\Gamma })\nonumber \\ =&\int d\varvec{\Gamma }\, \left\{ \begin{array}{c} \frac{m}{2} V^2\textbf{V}\\ \frac{I}{2} \omega ^2\textbf{V}\\ \frac{I}{2} (\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\end{array} \right\} C_{k\ell }(\varvec{\Gamma })=0. \end{aligned}$$

(44c)

Here,

$$\begin{aligned} a_{20}= & {} \frac{3}{5}\frac{\left\langle V^4\right\rangle }{\left\langle V^2\right\rangle ^2}-1, \quad a_{11}=\frac{\left\langle V^2\omega ^2\right\rangle }{\left\langle V^2\right\rangle \left\langle \omega ^2\right\rangle }-1, \quad b_{00}=3\frac{\left\langle (\textbf{V}\cdot \varvec{\omega })^2\right\rangle }{\left\langle V^2\right\rangle \left\langle \omega ^2\right\rangle }-1, \end{aligned}$$

(45)

are cumulants associated with $f^{(0)}$. Their expressions within the IRMM can be found in Ref. [73].¹ The cumulants are finite only if the coefficient of tangential restitution $\beta$ is larger than a certain $\alpha$-dependent threshold value $\beta _0(\alpha )$ [73].

Now, we multiply both sides of Eq. (22) by $\{\frac{m}{2} V^2\textbf{V},\frac{I}{2} \omega ^2\textbf{V},\frac{I}{2} (\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\}$, and integrate over velocity. The result is

$$\begin{aligned}&{\mathcal {D}}_t^{(0)}\begin{Bmatrix} \textbf{q}_t^{(1)}\\ \textbf{q}_r^{(1)}\\ \textbf{Q}^{(1)} \end{Bmatrix} + \nu ^{(0)} \begin{Bmatrix} \varphi _{30\mid 30}\textbf{q}_t^{(1)}+\frac{4}{\kappa }\varphi _{30\mid 12} [2\textbf{q}_r^{(1)}-\textbf{Q}^{(1)} ]\\ \frac{\kappa }{4}\varphi _{12\mid 30}\textbf{q}_t^{(1)}+\varphi _{12\mid 12}^{(1)}\textbf{q}_r^{(1)}+\varphi _{12\mid 12}^{(2)}\textbf{Q}^{(1)}\\ \overline{\varphi }_{12\mid 12}^{(1)}\textbf{Q}^{(1)}+\overline{\varphi }_{12\mid 12}^{(2)}\textbf{q}_r^{(1)} \end{Bmatrix}\nonumber \\&\quad = -\frac{nT\tau _t}{2m}\left\{ \begin{array}{c} 5\tau _t\left( 1+2a_{20}\right) \\ 3\tau _r\left( 1+2a_{11}\right) \\ \tau _r\left( 1+2b_{00}\right) \end{array} \right\} \nabla T -\frac{T^2\tau _t}{2m}\left\{ \begin{array}{c} 5\tau _t a_{20}\\ 3\tau _r a_{11}\\ \tau _r b_{00} \end{array} \right\} \nabla n. \end{aligned}$$

(46)

The solution to Eq. (46) has the structure

$$\begin{aligned} \begin{Bmatrix} \textbf{q}_t^{(1)}\\ \textbf{q}_r^{(1)}\\ \textbf{Q}^{(1)} \end{Bmatrix}=-\begin{Bmatrix} \tau _t\lambda _t\\ \tau _r\lambda _r\\ \lambda _Q \end{Bmatrix} \nabla T -\begin{Bmatrix} \mu _t\\ \mu _r\\ \mu _Q \end{Bmatrix} \nabla n, \end{aligned}$$

(47)

where, by dimensional analysis, $\{\lambda _t,\lambda _r,\lambda _Q\}\propto nT/\nu ^{(0)}$ and $\{\mu _t,\mu _r,\mu _Q\}\propto T^2/\nu ^{(0)}$. As a consequence,

$$\begin{aligned}{} & {} {\mathcal {D}}_t^{(0)}\begin{Bmatrix} \lambda _t\\ \lambda _r\\ \lambda _Q \end{Bmatrix} \nabla T=-\begin{Bmatrix} \lambda _t\\ \lambda _r\\ \lambda _Q \end{Bmatrix}\zeta ^{(0)}\left( 2\nabla T+\frac{T}{n}\nabla n\right) , \end{aligned}$$

(48a)

$$\begin{aligned}{} & {} {\mathcal {D}}_t^{(0)}\begin{Bmatrix} \mu _t\\ \mu _r\\ \mu _Q \end{Bmatrix} \nabla n= -\frac{3}{2}\begin{Bmatrix} \mu _t\\ \mu _r\\ \mu _Q \end{Bmatrix} \zeta ^{(0)}\nabla n. \end{aligned}$$

(48b)

Equation (46) then yields the matrix equations

$$\begin{aligned}{} & {} \begin{bmatrix} \tau _t\lambda _t\\ \tau _r\lambda _r\\ \lambda _Q\\ \end{bmatrix}=\frac{nT\tau _t}{2m}\left( \nu ^{(0)}\mathsf {\Phi }-2\zeta ^{(0)}\textsf{I}\right) ^{-1} \cdot \begin{bmatrix} 5\tau _t\left( 1+2a_{20}\right) \\ 3\tau _r\left( 1+2a_{11}\right) \\ \tau _r\left( 1+2b_{00}\right) \end{bmatrix}, \end{aligned}$$

(49a)

$$\begin{aligned}{} & {} \begin{bmatrix} \mu _t\\ \mu _r\\ \mu _Q\\ \end{bmatrix}= \frac{T}{n}\left( \nu ^{(0)}\mathsf {\Phi }-\frac{3}{2}\zeta ^{(0)}\textsf{I}\right) ^{-1} \cdot \left( \frac{n\tau _t T}{2m} \begin{bmatrix} 5\tau _t a_{20}\\ 3\tau _r a_{11}\\ \tau _r b_{00} \end{bmatrix} +\zeta ^{(0)} \begin{bmatrix} \tau _t\lambda _t\\ \tau _r\lambda _r\\ \lambda _Q\\ \end{bmatrix} \right) , \end{aligned}$$

(49b)

where $\textsf{I}$ is the $3\times 3$ unit matrix and

$$\begin{aligned} \mathsf {\Phi }\equiv \begin{bmatrix} \varphi _{30\mid 30}&{}\frac{8}{\kappa }\varphi _{30\mid 12}&{}-\frac{4}{\kappa } \varphi _{30\mid 12}\\ \frac{\kappa }{4} \varphi _{12\mid 30}&{} \varphi _{12\mid 12}^{(1)}&{} \varphi _{12\mid 12}^{(2)}\\ 0&{} \overline{\varphi }_{12\mid 12}^{(2)}&{} \overline{\varphi }_{12\mid 12}^{(1)} \end{bmatrix}. \end{aligned}$$

(50)

Equation (49a) provides the thermal heat conductivities $\lambda _t$ and $\lambda _r$, as well as the thermal torque-vorticity conductivity $\lambda _Q$. Then, the diffusive (heat and torque-vorticity) conductivities $\mu _t$, $\mu _r$, and $\mu _Q$ are obtained from Eq. (49b).

4 Discussion

4.1 Summary of Final Expressions

To summarize, in the NSF constitutive equations, Eqs. (11), the zeroth-order translational temperature and cooling rate are given by Eqs. (18) and (20), respectively, the bulk and shear viscosities are given by Eqs. (38) and (43b), respectively, and the cooling-rate transport coefficient $\xi$ is obtained from Eq. (39). The thermal and diffusive heat conductivities are $\lambda =\tau _t\lambda _t+\tau _r\lambda _r$ and $\mu =\mu _t+\mu _r$, respectively, where $\lambda _t$ and $\lambda _r$ are given by Eq. (49a), whereas $\mu _t$ and $\mu _r$ are obtained from Eq. (49b). In Eqs. (49), the zeroth-order cumulants ($a_{20}$, $a_{11}$, and $b_{00}$) have exact explicit expressions [73]. In all those results, the relevant coefficients $\chi _{ij\mid k\ell }$, $\psi _{ij\mid k\ell }$, $\varphi _{ij\mid k\ell }$, $\varphi _{12\mid 12}^{(1)}$, $\varphi _{12\mid 12}^{(2)}$, $\overline{\varphi }_{12\mid 12}^{(1)}$, and $\overline{\varphi }_{12\mid 12}^{(2)}$ can be found in the Appendix.

It is noteworthy that the first-order stress tensor $\mathsf {\Pi }^{(1)}$ turns out to be intrinsically coupled to the first-order spin-spin tensor $\mathsf {\Omega }^{(1)}$, as Eq. (40) shows. Thus, the derivation of the shear viscosity $\eta$ requires the parallel derivation of the spin viscosity $\eta _\Omega$ introduced in Eq. (41), this new transport coefficient being given by Eq. (43a). This remark is important because the coupling between $\mathsf {\Pi }^{(1)}$ and $\mathsf {\Omega }^{(1)}$ is usually ignored in the approximate derivation of $\eta$ for the IRHSM, where a standard Sonine approximation is employed [17, 19].

Analogously, Eq. (46) shows that the first-order partial heat fluxes $\textbf{q}_t^{(1)}$ and $\textbf{q}_r^{(1)}$ are intrinsically coupled to the first-order torque-vorticity vector $\textbf{Q}^{(1)}$. The latter vector is characterized by the transport coefficients $\lambda _Q$ and $\mu _Q$ [see Eq. (47)], which are obtained from Eqs. (49a) and (49b), respectively. While the coupling between $\textbf{q}_t^{(1)}$ and $\textbf{q}_r^{(1)}$ is taken into account in the approximate derivation of $\lambda$ and $\mu$ for the IRHSM [17, 19], the coupling to $\textbf{Q}^{(1)}$ is not accounted for.

In order to nondimensionalize the coefficients, it is convenient to take as a reference the shear viscosity and thermal conductivity of a gas made of elastic and smooth spheres at the same translational temperature as that of the HCS of the granular gas. More specifically, we take

$$\begin{aligned} \eta _0= \frac{5}{2}\frac{n\tau _tT}{\nu ^{(0)}},\quad \lambda _0=\frac{15}{4}\frac{\tau _t\eta _0}{m}. \end{aligned}$$

(51)

Thus, the dimensionless transport coefficients are

$$\begin{aligned} \eta ^*= & {} \frac{\eta }{\eta _0},\quad \eta _b^*=\frac{\eta _b}{\eta _0},\quad \eta _\Omega ^*=\frac{\eta _\Omega }{\eta _0}, \end{aligned}$$

(52a)

$$\begin{aligned} \lambda ^*= & {} \frac{\lambda }{\lambda _0},\quad \mu ^*=\frac{\mu }{(T/n)\lambda _0},\quad \lambda _Q^*=\frac{\lambda _Q}{\lambda _0},\quad \mu _Q^*=\frac{\mu _Q}{(T/n)\lambda _0}. \end{aligned}$$

(52b)

Note that the cooling-rate transport coefficient $\xi$ is dimensionless by construction.

Table 1 Special limiting cases

4.2 Special Limits

The (reduced) transport coefficients $\eta ^*$, $\eta _b^*$, $\xi$, $\eta _\Omega ^*$, $\lambda ^*$, $\mu ^*$, $\lambda _Q^*$, and $\mu _Q^*$ are displayed in Table 1 for some special limits: pure smooth particles ($\beta =-1$), the quasi-smooth limit ($\beta \rightarrow -1$), and the conservative Pidduck case ($\alpha =\beta =1$). In the case of pure smooth particles, the rotational degrees of freedom are irrelevant and the results are obtained by formally setting $\theta \rightarrow 0$. This allows us to recover previous results [53]. In the quasismooth limit, the zeroth-order cumulants diverge, what produces the divergence of the transport coefficients associated with the heat flux and the torque-vorticity vector, whereas those associated with the pressure and spin-spin tensors take finite values. Finally, in the Pidduck model [84], the total kinetic energy is conserved by collisions, the model having been used for polyatomic gases [80, 84,85,86,87]. Even in that case, the couplings $\mathsf {\Pi }\leftrightarrow \mathsf {\Omega }$ and $\textbf{q}\leftrightarrow \textbf{Q}$ still hold, as made evident by the nonzero values of $\eta _\Omega ^*$ and $\lambda _Q^*$, respectively. Interestingly, $\lambda _Q^*$ is independent of the reduced moment of inertia $\kappa$ in the Pidduck model.

Fig. 1

Plot of a the reduced shear viscosity $\eta ^*$, b the reduced bulk viscosity $\eta _b^*$, c the cooling-rate transport coefficient $\xi$, and d the reduced spin viscosity $-\eta _\Omega ^*$ versus the coefficient of tangential restitution $\beta$ for $\alpha =0.5, 0.6,0.7,0.8,0.9,1$. The arrows indicate increasing values of $\alpha$. Here, the particles are assumed to have a uniform mass distribution ($\kappa =\frac{2}{5})$

4.3 Coefficients $\eta ^*$, $\eta _b^*$, $\xi$, and $\eta _\Omega ^*$

The NSF transport coefficients related to the pressure tensor, the cooling rate, and the spin-spin tensor are plotted in Fig. 1 as functions of the coefficient of tangential restitution for several values of the coefficient of normal restitution. Here, and in what follows, we assume a uniform mass distribution in each particle, i.e., $\kappa =\frac{2}{5}$.

Figure 1a and b show a dependence of $\eta ^*$ and $\eta _b^*$ on $\alpha$ and $\beta$ qualitatively similar to that observed in the case of the IRHSM within the standard Sonine approximation [17]. For a fixed $\alpha$, $\eta ^*$ exhibits a local maximum at a certain value of $\beta$. As $\alpha$ decreases, the maximum grows and moves toward smaller values of $\beta$. In the case of $\eta _b^*$, however, the magnitude of the maximum decreases with decreasing $\alpha$, with the exception of $\alpha =1$. Moreover, as roughness increases, $\eta _b^*$ becomes less and less dependent on $\alpha$. In what concerns the cooling-rate transport coefficient $\xi$, Fig. 1c presents a nonmonotonic dependence on both $\alpha$ and $\beta$ similar to that observed in the IRHSM with the standard Sonine approximation [17].

The new transport coefficient $\eta _\Omega ^*$, which is discarded in the standard Sonine approximation for the IRHSM, is displayed in Fig. 1d. The most remarkable feature is that, in contrast to $\eta ^*$, $\eta _\Omega ^*$ is negative. To better grasp the consequences of that, imagine the simple shear flow geometry $\nabla _i u_j={\dot{\gamma }}\delta _{iy}\delta _{jx}$, where the shear rate ${\dot{\gamma }}$ is assumed to be positive. Then, according to Eq. (11a), $\left\langle V_x V_y\right\rangle <0$, that is, particles moving with $V_y>0$ tend to have $V_x<0$, and vice versa (second and fourth quadrants on the xy plane). On the other hand, Eq. (41) implies that $\left\langle \omega _x\omega _y\right\rangle >0$, so that the projection on the xy plane of the angular velocity of the particles tends to lie on the first and third quadrants. This is consistent with the orthogonality condition $\left\langle \textbf{V}\cdot \varvec{\omega }\right\rangle =0$ implied by $\Upsilon _{ij}^{(1)}=0$ [cf. Eqs. (29b) and (34)].

Figure 1d shows that the magnitude of $\eta _\Omega ^*$ grows quasilinearly with $\beta$. Moreover, at fixed $\beta$, $\mid \!\eta _\Omega ^*\!\mid$ increases with decreasing $\alpha$. All of this represents a smoother and more regular dependence on $\alpha$ and $\beta$ than in the cases of the coefficients $\eta$, $\eta _b$, and $\xi$.

Fig. 2

Plot of a the reduced thermal heat conductivity $\lambda ^*$, b the reduced diffusive heat conductivity $\mu ^*$, c the reduced thermal torque-vorticity conductivity $\lambda _Q^*$, and d the reduced diffusive torque-vorticity conductivity $\mu _Q^*$ versus the coefficient of tangential restitution $\beta$ for $\alpha =0.5, 0.6,0.7,0.8,0.9,1$. The arrows indicate increasing values of $\alpha$. Note the vertical asymptotes signaling the threshold values $\beta _0(\alpha )$ below which the coefficients diverge. Note also that, in panels c and d, absolute values are taken and the dashed lines refer to negative values. Here, the particles are assumed to have a uniform mass distribution ($\kappa =\frac{2}{5})$

4.4 Coefficients $\lambda$, $\mu$, $\lambda _Q$, and $\mu _Q$

Now we turn to the coefficients related to the heat flux and the torque-vorticity vector, as displayed in Fig. 2. Except for $\alpha =1$, the shapes of the curves of $\lambda ^*$ and $\mu ^*$ (see Figs. 2a and b) differ qualitatively from those obtained from the IRHSM with the standard Sonine approximation because of the divergence of the IRMM coefficients at $\beta =\beta _0(\alpha )$. As said before, such a divergence is induced by that of the HCS fourth-degree cumulants [73]. This implies a breakdown of a hydrodynamic description, in analogy with what happens with the IMM [53, 59, 88].

Figure 2c and d show that the transport coefficients associated with the torque-vorticity vector $\textbf{Q}$ (which is ignored in the standard Sonine approximation of the IRHSM) have a dependence on $\alpha$ and $\beta$ much more intricate than the heat-flux coefficients. In particular, both $\lambda _Q^*$ and $\nu _Q^*$ may attain negative values for intermediate and small values of $\beta$. This effect takes place even at $\alpha =1$, in which case $\lambda ^*$ and $\mu ^*$ remain finite for all $\alpha$.

Fig. 3

a–c Plot of $k_\perp ^*$ (dashed lines) and $k_{\Vert }^*$ (solid lines) versus the coefficient of tangential restitution $\beta$ for a $\alpha =1$, b $\alpha =0.9$, and c $\alpha =0.8$. The vertical dotted lines in panels b and c signal the divergence of $k_{\Vert }^*$. d Loci $\beta =\beta _0(\alpha )$ (solid line), $k_\perp ^*=k_{\Vert }^*$ (dashed line), and $\lambda ^*=\mu ^*$ (dashed-dotted line) in the plane $\alpha$ versus $\beta$. Here, the particles are assumed to have a uniform mass distribution ($\kappa =\frac{2}{5})$

4.5 Instability of the Homogeneous Cooling State

Depending on the values of $\alpha$, $\beta$, and $\kappa$, the HCS can be unstable versus long-wavelength perturbations. A standard linear stability analysis of the NSF hydrodynamic equations for the IRHSM [18, 20] can be straightforwardly extended to the IRMM. According to it, the critical wave number below which the HCS becomes unstable is $k_c=\max \{k_\perp ,k_{\Vert }\}$, where

$$\begin{aligned} k_\perp =\sqrt{\frac{nm\zeta ^{(0)}}{2\eta }},\quad k_{\Vert }=\sqrt{\frac{3n\zeta ^{(0)}}{2(\lambda -\mu n/T)}}. \end{aligned}$$

(53)

While $k<k_\perp$ signals the appearance of vortices, a clustering phenomenon is present if $k<k_{\Vert }$. The wave number $k_{\Vert }$ is not well defined if $\beta <\beta _0(\alpha )$ since both $\lambda$ and $\mu$ diverge in that region because of the divergence of the HCS cumulants.

Figure 3a–c show $k_\perp ^*$ and $k_{\Vert }^*$, where $k^*=k\sqrt{\tau _t T/m}/\nu ^{(0)}$, versus $\beta$ for $\alpha =1$, 0.9, and 0.8, respectively. We observe that, in the case $\alpha =1$, $k_c=k_\perp$, whereas $k_c=k_{\Vert }$ in the interval $-0.203<\beta <0.355$ if $\alpha =0.9$. This behavior is qualitatively similar to that found in the IRHSM within the standard Sonine approximation. On the other hand, in the case $\alpha =0.8$, not only $k_c=k_{\Vert }$ in the interval $-0.494<\beta <0.517$, but $k_c\rightarrow \infty$ in the inner interval $-0.283<\beta <0.308$. This divergence reflects the fact that $\lambda <\mu n/T$ in that region, implying the absolute instability of the HCS for any perturbation. It turns out that the IRHSM with the standard Sonine approximation also predicts a region of absolute instability but confined to a much smaller region [20].

Figure 3d display the loci $\beta =\beta _0(\alpha )$, $k_{\perp }^*=k_{\Vert }^*$, and $\lambda ^*=\mu ^*$ in the plane $\alpha$ versus $\beta$. In the region below the locus $\beta =\beta _0(\alpha )$, the HCS cumulants and the NSF heat-flux coefficients diverge; in the region below the locus $k_\perp ^*=k_{\Vert }^*$, the critical wavenumber is the longitudinal one (i.e., $k_c^*=k_{\Vert }^*$); finally, in the region below the locus $\lambda ^*=\mu ^*$, $k_c^*\rightarrow \infty$, implying the absolute instability of the HCS.

Fig. 4

Plot of the relative deviations a $\Delta \eta /\eta$, b $\Delta \lambda /\lambda$, and c $\Delta \mu /\mu$ versus the coefficient of tangential restitution $\beta$ for $\alpha =0.5, 0.6,0.7,0.8,0.9,1$. The arrows indicate increasing values of $\alpha$. The dashed lines in panels b and c represent the values of $\Delta \lambda /\lambda$ and $\Delta \mu /\mu$, respectively, at $\beta =\beta _0(\alpha )$. Here, the particles are assumed to have a uniform mass distribution ($\kappa =\frac{2}{5})$

4.6 Impact of the Couplings $\mathsf {\Pi }\leftrightarrow \mathsf {\Omega }$ and $\textbf{q}\leftrightarrow \textbf{Q}$

As said before, one of the strong points of the IRMM is that it unveils the inherent couplings $\mathsf {\Pi }\leftrightarrow \mathsf {\Omega }$ and $\textbf{q}\leftrightarrow \textbf{Q}$. It is then in order to assess the impact of those couplings on the NSF transport coefficients $\eta$, $\lambda$, and $\mu$. On the other hand, the coefficients $\eta _b$ and $\xi$ are not affected by the couplings.

In the case of the shear viscosity, if the coupling $\mathsf {\Pi }\leftrightarrow \mathsf {\Omega }$ were ignored, then one would forget the term $\psi _{20\mid 02}\Omega _{ij}^{(1)}$ in the first row of Eq. (40). This is equivalent to formally setting $\psi _{20\mid 02}\rightarrow 0$ in Eq. (43b), with the result

$$\begin{aligned} \widetilde{\eta }=\frac{n\tau _t T}{\nu ^{(0)}}\left[ \psi _{20\mid 20}-\frac{1}{2}\frac{\zeta ^{(0)}}{\nu ^{(0)}}\right] ^{-1}, \end{aligned}$$

(54)

where we have used a tilde to distinguish the approximate shear viscosity ($\widetilde{\eta }$) from the true one ($\eta$).

In the case of the heat-flux coefficients, ignorance of the coupling is equivalent to formally setting $\textbf{Q}^{(1)}\rightarrow 0$ in the first two rows of Eq. (46). Thus, instead of Eqs. (49), one would have

$$\begin{aligned} \begin{bmatrix} \tau _t\widetilde{\lambda }_t\\ \tau _r\widetilde{\lambda }_r \end{bmatrix}= & {} \frac{nT\tau _t}{2m}\left( \nu ^{(0)}\widetilde{\mathsf {\Phi }}-2\zeta ^{(0)}\widetilde{\textsf{I}}\right) ^{-1} \cdot \begin{bmatrix} 5\tau _t\left( 1+2a_{20}\right) \\ 3\tau _r\left( 1+2a_{11}\right) \\ \end{bmatrix}, \end{aligned}$$

(55a)

$$\begin{aligned} \begin{bmatrix} \widetilde{\mu }_t\\ \widetilde{\mu }_r \end{bmatrix}= & {} \frac{T}{n}\left( \nu ^{(0)}\widetilde{\mathsf {\Phi }}-\frac{3}{2}\zeta ^{(0)}\widetilde{\textsf{I}}\right) ^{-1} \cdot \left( \frac{n\tau _t T}{2m} \begin{bmatrix} 5\tau _t a_{20}\\ 3\tau _r a_{11} \end{bmatrix} +\zeta ^{(0)} \begin{bmatrix} \tau _t\lambda _t\\ \tau _r\lambda _r \end{bmatrix} \right) , \end{aligned}$$

(55b)

where $\widetilde{\textsf{I}}$ is the $2\times 2$ unit matrix and

$$\begin{aligned} \widetilde{\mathsf {\Phi }}\equiv \begin{bmatrix} \varphi _{30\mid 30}&{}\frac{8}{\kappa }\varphi _{30\mid 12}\\ \frac{\kappa }{4} \varphi _{12\mid 30}&{} \varphi _{12\mid 12}^{(1)} \end{bmatrix}. \end{aligned}$$

(56)

The relative deviations $\Delta \eta /\eta$, $\Delta \lambda /\lambda$, and $\Delta \mu /\mu$, (where $\Delta \eta =\widetilde{\eta }-\eta$, $\Delta \lambda =\widetilde{\lambda }-\lambda$, and $\Delta \mu =\widetilde{\mu }-\mu$) are plotted in Fig. 4. As can be seen, the approximate shear viscosity ${\tilde{\eta }}$ underestimates the true value of $\eta$, this effect increasing monotonically with increasing inelasticity and, especially, increasing roughness. At $\beta =1$, the deviation of ${\tilde{\eta }}$ from $\eta$ ranges from about 9% for $\alpha =0.5$ to about 8% for $\alpha =1$ (Pidduck’s limit).

In what concerns $\Delta \lambda /\lambda$ and $\Delta \mu /\mu$, the combined influence of $\alpha$ and $\beta$ is much more intricate than in the case of $\Delta \eta /\eta$. In the region of high roughness (say $\beta >0.4$), both transport coefficients are overestimated if $\textbf{Q}$ is ignored, this effect increasing again monotonically with increasing inelasticity and roughness. However, at intermediate roughness (say $-0.4<\beta <0.4$), ${\tilde{\lambda }}$ and ${\tilde{\mu }}$ typically underestimate the true values. As the threshold values $\beta _0(\alpha )$ are approached, the relative deviations $\Delta \lambda /\lambda$ and $\Delta \mu /\mu$ reach finite values. At $\beta =1$, the deviations of the pair $({\tilde{\lambda }},{\tilde{\mu }})$ range from about (8%,3%) for $\alpha =0.5$ to about (6%,2%) for $\alpha =1$ (Pidduck’s limit). Note that, although at $\alpha =1$ both $\mu$ and $\widetilde{\mu }$ tend to 0 in the limit $\beta \rightarrow 1$, the relative deviation $\Delta \mu /\mu$ tends to a finite value in that limit.

Except for extremely low values of $\alpha$, the consequences of disregarding the couplings $\mathsf {\Pi }\leftrightarrow \mathsf {\Omega }$ and $\textbf{q}\leftrightarrow \textbf{Q}$ are generally not substantial. However, to maintain consistency, these couplings must be considered.

4.7 Proposal of an Augmented Sonine Approximation for the IRHSM

As said before, the derivation of the NSF transport coefficients in the IRHSM requires the use of approximations, since taking velocity moments in the kinetic equations for $f^{(0)}$ and $f^{(1)}$ generates an infinite hierarchy of moment equations. In the standard Sonine approximation [17, 19], the unknown vector functions $\varvec{{\mathcal {A}}}$ and $\varvec{{\mathcal {B}}}$ appearing in Eq. (25) are approximated by a Maxwellian times a linear combination of $\textbf{V}$, $V^2\textbf{V}$, and $\omega ^2\textbf{V}$. Likewise, the unknown tensor function ${\mathcal {C}}_{ij}$ is approximated by a Maxwellian times $V_iV_j-\frac{1}{3}V^2\delta _{ij}$. In the light of the results we have obtained in this paper for the IRMM, we propose here for the IRHSM the more consistent Sonine approximation

$$\begin{aligned}{} & {} \varvec{{\mathcal {A}}}\rightarrow -f_M \left[ \gamma _{A_t}\left( V^2-\frac{5\tau _t T}{m}\right) \textbf{V}+\gamma _{A_r}\left( \omega ^2-\frac{3\tau _r T}{I}\right) \textbf{V}+\gamma _{A_Q}(\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\right] , \end{aligned}$$

(57a)

$$\begin{aligned}{} & {} \varvec{{\mathcal {B}}}\rightarrow -f_M\left[ \gamma _{B_t}\left( V^2-\frac{5\tau _t T}{m}\right) \textbf{V}+\gamma _{B_r}\left( \omega ^2-\frac{3\tau _r T}{I}\right) \textbf{V}+\gamma _{B_Q}(\textbf{V}\cdot \varvec{\omega })\varvec{\omega }\right] , \end{aligned}$$

(57b)

$$\begin{aligned}{} & {} {\mathcal {C}}_{ij}\rightarrow -f_M\left[ {\gamma _{C_t}} \left( V_iV_j-\frac{1}{3}V^2\delta _{ij}\right) +{\gamma _{C_r}} \left( \omega _i\omega _j-\frac{1}{3}\omega ^2\delta _{ij}\right) \right] , \end{aligned}$$

(57c)

where

$$\begin{aligned} f_M=n\left( \frac{m I}{4\pi ^2\tau _t\tau _r T^2}\right) ^{\frac{3}{2}}\exp \left( -\frac{mV^2}{2\tau _t T}-\frac{I\omega ^2}{2\tau _r T}\right) . \end{aligned}$$

(58)

In the standard Sonine approximation [17, 19], one assumes $\gamma _{A_Q}=\gamma _{B_Q}=\gamma _{C_r}=0$. On the other hand, in the augmented version given by Eqs. (57), the coefficient $\gamma _{C_r}$ couples $\mathsf {\Pi }$ and $\mathsf {\Omega }$, whereas the coefficients $\gamma _{A_Q}$ and $\gamma _{B_Q}$ couple $\textbf{q}_t$ and $\textbf{q}_r$ to $\textbf{Q}$.

It is worthwhile noting that, in the case of hard disks on a plane [19], one has $\textbf{V}\perp \varvec{\omega }$ and $\varvec{\omega }\Vert \widehat{\textbf{z}}$, so that $(\textbf{V}\cdot \varvec{\omega })\varvec{\omega }=0$ and $\omega _i\omega _j-\frac{1}{3}\omega ^2\delta _{ij}=0$. Therefore, in that case the augmented Sonine approximation becomes identical to the standard one.

5 Conclusions

Quoting Ernst and Brito [34, 35], one can say that “what harmonic oscillators are for quantum mechanics, and dumb-bells for polymer physics, that is what elastic and inelastic Maxwell models are for kinetic theory.” The present paper, in conjunction with Ref. [73], aims to extend Ernst and Brito’s dictum into the domain of models featuring rough particles, specifically the IRMM.

Building upon the exact collisional production rates established in Ref. [73], in this work we have derived the NSF transport coefficients as explicit functions of the coefficients of normal ($\alpha$) and tangential ($\beta$) restitution, along with the reduced moment of inertia ($\kappa$). The resulting hydrodynamic constitutive equations, Eqs. (11), encompass the shear viscosity ($\eta$), the bulk viscosity ($\eta _b$), the thermal heat conductivity ($\lambda$), the diffusive heat conductivity ($\mu$), and the cooling-rate transport coefficient ($\xi$). The evaluation of $\eta$ requires the parallel evaluation of the spin viscosity $\eta _\Omega$ [see Eq. (41)], whereas the evaluation of $\lambda$ and $\mu$ requires the parallel evaluation of the torque-vorticity coefficients $\lambda _Q$ and $\mu _Q$ [see Eq. (47)]. Interestingly, Figs. 1 and 2 illustrate the complex dependence of these eight transport coefficients on $\alpha$ and $\beta$ for a uniform mass distribution ($\kappa =\frac{2}{5}$).

Our analysis reveals a noteworthy result as both the mean spin vector and the couple stress tensor vanish to the NSF order, contrary to the case of micropolar fluids [74, 75, 78, 79] and dense gases [80,81,82,83]. This discrepancy is attributed to the low-density regime and the absence of boundary effects in the NSF description.

While the coefficients tied to pressure and spin-spin tensors remain finite for any $\alpha$ and $\beta$, those linked to the heat flux and the torque-vorticity vector exhibit dependency on the HCS cumulants defined in Eqs. (45), diverging if $\beta$ is less than an $\alpha$-dependent value, $\beta _0(\alpha )$ [73]. This divergence, originating from an algebraic high-velocity tail in the HCS distribution function [34,35,36], manifests in the breakdown of hydrodynamics if $\beta <\beta _0(\alpha )$. Beyond this region, the NSF hydrodynamic framework aids in determining the instability of the HCS against weak inhomogeneous perturbations. That instability occurs if the wave number of the perturbations is smaller than $k_c=\max \{k_\perp ,k_{\Vert }\}$ (see Fig. 3). Remarkably, a dome-shaped region emerges in the $\alpha$ vs $\beta$ plane, within which $k_c\rightarrow \infty$, signifying an absolute instability of the HCS.

Additionally, considering the couplings $\mathsf {\Pi }\leftrightarrow \mathsf {\Omega }$ and $\textbf{q}\leftrightarrow \textbf{Q}$, their impact on the NSF coefficients is found to be below 10% if $\alpha >0.5$ and $\kappa =\frac{2}{5}$ (see Fig. 4). Nevertheless, a consistent treatment should incorporate these couplings, even within the IRHSM.

Although, for simplicity and conciseness, this paper presents graphs exclusively for a uniform mass distribution ($\kappa =\frac{2}{5}$), the reduced moment of inertia $\kappa$ serves as an additional control parameter influencing results. The nontrivial impact of $\kappa$ on the final outcomes can be exemplified by the Pidduck limit ($\alpha =\beta =1$). According to the fourth column of Table 1, the reduced shear viscosity, bulk viscosity, and thermal heat conductivity take the values $(\eta ^*,\eta _b^*,\lambda ^*)=(1,\infty ,1.31)$, (0.92, 0.61, 1.23), (0.91, 0.25, 1.29), and (0.98, 0.21, 1.30) for $\kappa =0$, $\frac{1}{10}$, $\frac{2}{5}$, and $\frac{2}{3}$, respectively.

In summary, the IRMM serves as a compelling mathematical model, providing exact results and offering insights into the intricacies anticipated in the more realistic IRHSM. Future work aims to explore the exact non-Newtonian properties of the IRMM under simple shear flow.

Explicit Expressions for the Coefficients Appearing in Eqs. (7)

As explained in Ref. [73], we have adopted a specific criterion for coefficient notation. Firstly, assume that $\Psi _{k_1k_2}(\varvec{\Gamma })$ represents a homogeneous velocity polynomial with degrees $k_1$ and $k_2$ relative to $\textbf{V}$ and $\varvec{\omega }$, respectively. Thus, a coefficient such as $Y_{k_1 k_2\mid \ell _1\ell _2}$ (with $Y=\chi$, $\varphi$, or $\psi$) corresponds to the collisional moment $\mathcal {J}[\Psi _{k_1k_2}]$, linked to a product like $\langle \Psi _{i_1i_2}(\varvec{\Gamma })\rangle \langle \Psi _{j_1j_2}(\varvec{\Gamma })\rangle$ with $i_1+j_1=\ell _1$ and $i_2+j_2=\ell _2$. If there exists more than one $Y_{k_1 k_2\mid \ell _1\ell _2}$ for a given $\mathcal {J}[\Psi _{k_1k_2}]$, each is identified by a superscript. Lastly, Greek letters $\chi$, $\varphi$, and $\psi$ signify coefficients in $\mathcal {J}[\Psi _{k_1k_2}]$ linked to scalar, vector, and tensor quantities $\Psi _{k_1k_2}(\varvec{\Gamma })$, respectively; additionally, an overline is used if $\Psi _{k_1k_2}(\varvec{\Gamma })$ involves the inner product $\textbf{V}\cdot \varvec{\omega }$.

The 17 coefficients appearing in Eqs. (7) are given by

$$\begin{aligned} \varphi _{01\mid 01}=&\frac{4{\widetilde{\beta }}}{3\kappa },\quad \psi _{11\mid 11}=\frac{1}{3}\left( {\widetilde{\alpha }}+2{\widetilde{\beta }}\frac{1+\kappa }{\kappa }\right) , \end{aligned}$$

(A1a)

$$\begin{aligned} \chi _{20\mid 20}=&\frac{2}{3}\left[ {\widetilde{\alpha }}\left( 1-{\widetilde{\alpha }}\right) +2{\widetilde{\beta }}\left( 1-{\widetilde{\beta }}\right) \right] ,\quad \chi _{20\mid 02}=-\frac{{\widetilde{\beta }}^2}{3},\end{aligned}$$

(A1b)

$$\begin{aligned} \chi _{02\mid 02}=&\frac{4{\widetilde{\beta }}}{3\kappa }\left( 1-\frac{{\widetilde{\beta }}}{\kappa }\right) ,\quad \chi _{02\mid 20}=-\frac{16{\widetilde{\beta }}^2}{3\kappa ^2},\end{aligned}$$

(A1c)

$$\begin{aligned} \psi _{20\mid 20}=&\frac{2}{15}\left( 5{\widetilde{\alpha }}-2{\widetilde{\alpha }}^2-6{\widetilde{\alpha }}{\widetilde{\beta }}+10{\widetilde{\beta }}-7{\widetilde{\beta }}^2\right) ,\end{aligned}$$

(A1d)

$$\begin{aligned} \psi _{20\mid 02}=&\frac{{\widetilde{\beta }}^2}{6}, \quad \psi _{02\mid 20}=\frac{8{\widetilde{\beta }}^2}{3\kappa ^2},\quad \psi _{02\mid 02}=\frac{2{\widetilde{\beta }}}{15\kappa }\left( 10-\frac{7{\widetilde{\beta }}}{\kappa }\right) , \end{aligned}$$

(A1e)

$$\begin{aligned} \varphi _{30\mid 30}=&\frac{1}{15}\left( 15{\widetilde{\alpha }}-11{\widetilde{\alpha }}^2-8{\widetilde{\alpha }}{\widetilde{\beta }}+30{\widetilde{\beta }}-26{\widetilde{\beta }}^2\right) ,\end{aligned}$$

(A1f)

$$\begin{aligned} \varphi _{30\mid 12}=&-\frac{{\widetilde{\beta }}^2}{6},\quad \varphi _{12\mid 30}=-\frac{8{\widetilde{\beta }}^2}{3\kappa ^2}\end{aligned}$$

(A1g)

$$\begin{aligned} \varphi _{12\mid 12}^{(1)}=&\frac{1}{15}\left[ 5{\widetilde{\alpha }}+10{\widetilde{\beta }}+\frac{2{\widetilde{\beta }}}{\kappa }\left( 10-4{\widetilde{\alpha }}-11{\widetilde{\beta }}-\frac{5{\widetilde{\beta }}}{\kappa }\right) \right] ,\end{aligned}$$

(A1h)

$$\begin{aligned} \overline{\varphi }_{12\mid 12}^{(1)}=&\frac{1}{15}\left[ 5{\widetilde{\alpha }}+10{\widetilde{\beta }}+\frac{{\widetilde{\beta }}}{\kappa }\left( 20-3{\widetilde{\alpha }}-7{\widetilde{\beta }}\frac{1+\kappa }{\kappa }\right) \right] ,\end{aligned}$$

(A1i)

$$\begin{aligned} \varphi _{12\mid 12}^{(2)}=&\frac{2{\widetilde{\beta }}}{15\kappa }\left( 2{\widetilde{\alpha }}+3{\widetilde{\beta }}\right) ,\quad \overline{\varphi }_{12\mid 12}^{(2)}=\frac{{\widetilde{\beta }}}{15\kappa }\left( {\widetilde{\alpha }}-{\widetilde{\beta }}\frac{1+\kappa }{\kappa }\right) , \end{aligned}$$

(A1j)

where

$$\begin{aligned} {\widetilde{\alpha }}\equiv {1+\alpha \over 2},\quad {\widetilde{\beta }}\equiv \frac{1+\beta }{2}\frac{\kappa }{1+\kappa }. \end{aligned}$$

(A2)