\section{Introduction}
\label{sec1}


As is well known, the prototypical model of a granular gas is a system composed of smooth, frictionless hard spheres which collide inelastically with a constant coefficient of normal restitution $0<\alpha\leq 1$ \cite{BP04,C90,G03}. In the dilute and moderately dense regimes, the microscopic description of the gas is given by the one-particle velocity distribution function $f$ obeying the (inelastic) Boltzmann and Enskog kinetic equations \cite{BP04,GS95,BDS97}. At a more phenomenological level, the gas can also be described by the Navier--Stokes--Fourier (NSF) hydrodynamic equations for the densities of mass, momentum, and energy with appropriate constitutive equations for the stress tensor, heat flux, and cooling rate. The Chapman--Enskog method \cite{CC70,FK72} bridges the gap between the kinetic and hydrodynamic descriptions, thus providing explicit expressions for the NSF transport coefficients in terms of the coefficient of normal restitution. This task was first accomplished in the quasismooth limit \cite{LSJC84,JR85b,SG98}, the results being subsequently extended to finite degree of inelasticity for monocomponent \cite{BDKS98,GD99} and multicomponent \cite{GD02,GDH07,GHD07,MGH12} granular gases.

In spite of the interest and success of the smooth hard-sphere model of granular gases, grains in nature are typically frictional, and hence energy transfer between the translational and rotational degrees of freedom occurs upon particle collisions. The simplest model accounting for particle roughness (and thus including the particle angular velocity as an additional mechanical variable) neglects the effect of sliding collisions and is characterized by a constant coefficient of tangential restitution $\beta$ \cite{JR85}. This parameter ranges from $-1$ (perfectly smooth spheres) to $1$ (perfectly rough spheres). A more sophisticated model incorporates the Coulomb friction coefficient as a third collision constant, so that collisions become sliding beyond a certain impact parameter \cite{W93,FLCA94,GNB05,GNB05b}. On the other hand, this three-parameter model significantly complicates the kinetic description, while the simpler two-parameter $(\alpha,\beta)$ model captures the essential features of granular flows when particle rotations are relevant. This explains the wide use of the latter model in the literature
\cite{JR85,LS87,C89,L91,LB94,L95,GS95,L96,HZ97,ZTPSH98,ML98,LHMZ98,HHZ00,AHZ01,MHN02,CLH02,JZ02,PZMZ02,GNB05,%
GNB05b,Z06,BPKZ07,GA08,KBPZ09,SKG10,N11,SKS11,BB12,MDHEH13,VSK14}.

Needless to say, an important challenge is the derivation of the NSF hydrodynamic equations of a granular gas of inelastic \emph{rough} hard spheres, with explicit expressions for the transport coefficients as functions of $\alpha$ and $\beta$. Previous attempts have been restricted to nearly elastic collisions ($\alpha\lesssim 1$) and  either nearly smooth particles ($\beta\gtrsim -1$) \cite{JR85,GNB05,GNB05b} or nearly perfectly rough particles ($\beta\lesssim 1$) \cite{JR85,L91}. The goal of this paper is to  uncover the whole range of values of the two coefficients of restitution $\alpha$ and $\beta$ and derive explicit expressions for the NSF transport coefficients of a \emph{dilute} granular gas beyond the above limiting situations.

In the case of conventional gases (i.e., when energy is conserved upon collisions), the set of hydrodynamic variables is  related to densities of conserved quantities, namely, the particle density $n$ (conservation of mass), the flow velocity $\mathbf{u}$ (conservation of momentum), and the temperature $T$ (conservation of energy). If the particles are perfectly elastic and rough ($\alpha=\beta=1$), what is conserved is the sum of the translational and the rotational kinetic energies and thus the granular temperature $T$ has  translational ($T_t$) and  rotational ($T_r$) contributions \cite{P22,MSD66,CC70}. Moreover, since the angular velocity of the particles is not a collisional invariant, the mean spin $\boldsymbol{\Omega}$ is not included either in the set of hydrodynamic variables or in the definition of the rotational contribution ($T_r$) to the temperature \cite{P22,MSD66,CC70}.

For granular gases, although the total kinetic energy is dissipated by collisions,   the granular temperature is typically included as a hydrodynamic field  in most studies, and, consequently, a sink term appears in the corresponding balance equation. For nearly smooth spheres ($\beta\gtrsim -1$), some authors \cite{JR85,GNB05,GNB05b} have chosen (in addition to $n$ and $\mathbf{u}$) the two partial contributions $T_t$ and $T_r$ to the temperature, as well as the mean spin $\boldsymbol{\Omega}$,  as hydrodynamic variables. On the other hand, in this paper we will choose as hydrodynamic fields for dissipative gases the same as in conservative systems, i.e., $n$, $\mathbf{u}$, and $T$. In this way, the hydrodynamic description encompasses the conservative gases as special limits. The advantage of the choice of the set $\{n,\mathbf{u},T\}$ instead of the set $\{n,\mathbf{u},\boldsymbol{\Omega},T_t,T_r\}$ is analogous to the advantage of the set $\{n_1,n_2,\mathbf{u},T\}$ instead of $\{n_1,n_2,\mathbf{u}_1,\mathbf{u}_2,T_1,T_2\}$ in a binary mixture of inelastic smooth hard spheres, as discussed in Refs.\ \cite{GMD06,DB11}.

The plan of the paper is as follows. Section \ref{sec2} is devoted to the definition of the model of inelastic rough hard spheres and their description in the low-density regime by means of the Boltzmann equation. The  exact balance equations for the densities of mass, momentum, and energy are obtained from the Boltzmann equation, and the associated fluxes of momentum and energy, as well as the cooling rate, are identified.
The so-called homogeneous cooling state (HCS) is studied in Sec.\ \ref{sec3}. Special attention is paid to the time evolution of the mean spin $\boldsymbol{\Omega}$.  While the temperature ratio $T_r/T_t$ reaches a well-defined value for long times, the ratio $I\Omega^2/T_r$ ($I$ being the moment of inertia) decays to zero with a characteristic time typically smaller {(except near $\beta= -1$)} than the relaxation time of the temperature ratio (see Fig.\ \ref{fig1}). This clearly justifies the exclusion of $\boldsymbol{\Omega}$ as a hydrodynamic field. Next, the Chapman--Enskog method is applied in Sec.\ \ref{sec4} to derive the linear integral equations for the velocity-dependent functions characterizing the distribution function to first-order in the hydrodynamic gradients. In Sec.\ \ref{sec5},  the NSF transport coefficients (shear and bulk viscosities, thermal conductivity, Dufour-like coefficient, and cooling rate transport coefficient) are expressed in terms of integrals involving the solutions of the linear integral equations. Next, a first Sonine polynomial approximation is used to obtain practical results from this formulation, thus providing explicit forms for the transport coefficients as nonlinear functions of both $\alpha$ and $\beta$ (see Table \ref{table1}). Some technical details of the calculations are relegated to the Appendix.
The results are discussed in Sec.\ \ref{sec6}, where known expressions in the limiting cases of inelastic smooth spheres ($\alpha<1$, $\beta=-1$) and perfectly elastic and rough spheres ($\alpha=\beta=1$) are recovered. The intricate dependence of the transport coefficients on both coefficients of restitution is illustrated by some representative cases. Finally, the paper closes with some concluding remarks in Sec.\ \ref{sec7}.

{\renewcommand{\arraystretch}{2}
\begin{table}
   \caption{Summary of explicit expressions.}\label{table1}
\begin{ruledtabular}
\begin{tabular}{l}
$\displaystyle{\widetilde\alpha={1+\alpha\over2}}$, $\displaystyle{\widetilde\beta={1+\beta\over2}{\kappa\over\kappa+1}}$\\
$\displaystyle{\frac{T_t^{(0)}}{T}=\tau_t={2\over1+\theta_\hcs}}$, $\displaystyle{\frac{T_r^{(0)}}{T}=\tau_r={2\theta_\hcs\over1+\theta_\hcs}}$\\
$\displaystyle{\theta_\hcs=\sqrt{1+h^2}+h}$, $\displaystyle{h\equiv \frac{(1+\kappa)^2}{2\kappa(1+\beta)^2}\left[{1-\alpha^2}-(1-\beta^2)\frac{1-\kappa}{1+\kappa}\right]}$\\
$\displaystyle{\nu=\frac{16}{5}\sigma^2 n\sqrt{\pi \tau_tT/m}}$\\
$\displaystyle{\frac{\zeta^{(0)}}{\nu}=\zeta^*=\frac{5}{12}\frac{1}{1+\theta_\hcs}\left[1-\alpha^2+({1-\beta^2})
\frac{\theta_\hcs+\kappa}{1+\kappa}\right]}$\\[2mm]
\hline
 $\displaystyle{\eta=\frac{n\tau_t T}{\nu}\frac{1}{\nu_\eta^*-\frac{1}{2}\zeta^*}}$\\
$\displaystyle{\eta_b=\frac{n\tau_t\tau_r T}{\nu}\gamma_E}$\\
$\displaystyle{\lambda={\tau_t\lambda_t+\tau_r\lambda_r}}$, $\displaystyle{\lambda_t=\frac{5}{2}\frac{n\tau_t T}{m\nu}\gamma_{A_t}}$, $\displaystyle{\lambda_r=\frac{3}{2}\frac{n\tau_t T}{m\nu}\gamma_{A_r}}$\\
 $\displaystyle{\mu={\mu_t+\mu_r}}$, $\displaystyle{\mu_t=\frac{5}{2}\frac{\tau_t^2 T^2}{m\nu}\gamma_{B_t}}$, $\displaystyle{\mu_r=\frac{3}{2}\frac{\tau_t \tau_r T^2}{m\nu}\gamma_{B_r}}$\\
$\displaystyle{\xi={\frac{1}{2}\left(\tau_t\xi_t+\tau_r\xi_r\right)}=\gamma_E\Xi}$, $\displaystyle{\xi_t=\gamma_E\Xi_t}$, $\displaystyle{\xi_r=\gamma_E\Xi_r}$\\[1mm]
\hline
$\displaystyle{\nu_\eta^*=(\widetilde\alpha+\widetilde\beta)(2-\widetilde\alpha-\widetilde\beta)+\frac{\widetilde\beta^2\theta_\hcs}{6\kappa}}$\\
$\displaystyle{\gamma_E=\frac{2}{3}({\Xi_t-\Xi_r-\zeta^*})^{-1}}$\\
$\displaystyle{\Xi_t=\frac{5}{8}\tau_r\Big[1-{\alpha }^2+(1-{\beta }^2)\frac{\kappa}{1+\kappa}-\frac{\kappa}{3}({\theta_\hcs-5})\left(\frac{1+{\beta }}{1+\kappa}\right)^2
\Big]}$\\
$\displaystyle{\Xi_r=\frac{5}{8}\tau_t\frac{1+\beta}{1+\kappa}\left[\frac{\theta_\hcs-2}{3}(1-\beta)
+\frac{\kappa}{3}({\theta_\hcs-5})\frac{1+{\beta }}{1+\kappa}\right]}$\\
$\displaystyle{\Xi=\frac{5}{16}\tau_t\tau_r\left[1-{\alpha }^2+(1-{\beta }^2)\left(1+\frac{1}{3}\frac{\theta_\hcs-5}{1+\kappa}\right)\right]}$\\
$\displaystyle{\gamma_{A_t}=\frac{Z_r-Z_t-2\zeta^*}{\left(Y_t-2\zeta^*\right)\left(Z_r-2\zeta^*\right)-Y_r Z_t}}$\\
$\displaystyle{\gamma_{A_r}=\frac{Y_t-Y_r-2\zeta^*}{\left(Y_t-2\zeta^*\right)\left(Z_r-2\zeta^*\right)-Y_r Z_t}}$\\
$\displaystyle{\gamma_{B_t}=\zeta^*\frac{\gamma_{A_t}\left(Z_r-\frac{3}{2}\zeta^*\right)-\gamma_{A_r}Z_t}
{\left(Y_t-\frac{3}{2}\zeta^*\right)\left(Z_r-\frac{3}{2}\zeta^*\right)-Y_r Z_t}}$\\
$\displaystyle{\gamma_{B_r}=\zeta^*\frac{\gamma_{A_r}\left(Y_t-\frac{3}{2}\zeta^*\right)-\gamma_{A_t}Y_r}
{\left(Y_t-\frac{3}{2}\zeta^*\right)\left(Z_r-\frac{3}{2}\zeta^*\right)-Y_r Z_t}}$\\
$\displaystyle{Y_t=\frac{41}{12}\left(\widetilde\alpha+\widetilde\beta\right)-\frac{33}{12}\left(\widetilde\alpha^2+\widetilde\beta^2\right)-\frac{4}{3}\widetilde\alpha\widetilde\beta-\frac{7\theta_\hcs}{12}
\frac{\widetilde\beta^2}{\kappa}}$\\
$\displaystyle{Z_t=-\frac{5\theta_\hcs}{6}
\frac{\widetilde\beta^2}{\kappa}}$\\
$\displaystyle{Y_r=\frac{25}{36}\frac{\widetilde\beta}{\kappa}\left(1-3\frac{\widetilde\beta}{\theta_\hcs}-\frac{\widetilde\beta}{\kappa}\right)}$\\
$\displaystyle{Z_r=\frac{5}{6}\left(\widetilde\alpha+\widetilde\beta\right)+\frac{5}{18}\frac{\widetilde\beta}{\kappa}\left(7-3\frac{\widetilde\beta}{\kappa}-6{\widetilde\beta}-
 4{\widetilde\alpha}\right)}$\\
     \end{tabular}
 \end{ruledtabular}
 \end{table}
}


\section{Granular gas of inelastic rough hard spheres: Boltzmann description}
\label{sec2}


\subsection{Collision rules}


Let  $(\mathbf{v}, \mathbf{v}_1)$ and $(\boldsymbol{\omega}, \boldsymbol{\omega}_1)$ denote the linear and the angular  precollisional velocities, respectively, of two rough spherical particles with the same mass $m$, diameter $\sigma $, and moment of inertia $I$, while $(\mathbf{v}^\prime, \mathbf{v}_1^\prime)$ and $(\boldsymbol{\omega}^\prime,\boldsymbol{\omega}_1^\prime)$ correspond to their postcollisional velocities. The pre- and postcollisional velocities are related by
\begin{subequations}
\label{1,2}
\bal
\label{1}
 m \mathbf{v}^\prime=m\mathbf{v} -\mathbf{Q},\quad &I \boldsymbol{\omega}^\prime=I\boldsymbol{\omega} -{\sigma \over2}\,\hnn{\sigma}\times\mathbf{Q},\\
 \label{2}
 m \mathbf{v}_1^\prime=m\mathbf{v}_1 +\mathbf{Q},\quad
 &I \boldsymbol{\omega}_1^\prime=I\boldsymbol{\omega}_1 -{\sigma \over2}\,\hnn{\sigma}\times\mathbf{Q},
 \eal
\end{subequations}
where ${\mathbf{Q}}$ denotes the impulse exerted by the unlabeled particle on the labeled one and $\hnn{\sigma}$ is the unit collision vector  joining the centers of the two colliding spheres and pointing from the center of the {unlabeled} particle to the center of the {labeled one}. Furthermore, the relationship between the center-of-mass relative velocities $(\mathbf{g}=\mathbf{v}-\mathbf{v}_1,\mathbf{g}'=\mathbf{v}'-\mathbf{v}_1')$ and the  relative velocities $(\overline{\mathbf{g}},\overline{\mathbf{g}}')$ of the points of the spheres which are in contact during a binary encounter are
\begin{subequations}
\label{3}
\beq
 \overline{\mathbf{g}}=\mathbf{g}-{\sigma \over2}\hnn{\sigma}\times(\boldsymbol{\omega}+\boldsymbol{\omega}_1),
 \label{3a}
 \eeq
 \beq
 \label{3b}
 \overline{\mathbf{g}}^\prime=\mathbf{g}^\prime-{\sigma \over2}\hnn{\sigma}\times(\boldsymbol{\omega}^\prime+\boldsymbol{\omega}_1^\prime).
 \eeq
\end{subequations}

Combining Eqs.\ \eqref{1,2} and \eqref{3}, one obtains
\bal
\overline{\mathbf{g}}'=&\overline{\mathbf{g}}-\frac{2}{m}\mathbf{Q}+\frac{2}{m\kappa}\hnn{\sigma}\times(\hnn{\sigma}\times\mathbf{Q})\nonumber\\
=&\overline{\mathbf{g}}-\frac{2}{m}{\kappa+1\over \kappa}\mathbf{Q}+\frac{2}{m\kappa}(\hnn{\sigma}\cdot\mathbf{Q})\hnn{\sigma},
\label{3bg}
\eal
where in the second step use has been made of the mathematical property $\hnn{\sigma}\times(\hnn{\sigma}\times\mathbf{Q})=(\hnn{\sigma}\cdot\mathbf{Q})\hnn{\sigma}-\mathbf{Q}$ and
\beq
\kappa={4I\over m\sigma ^2}
\eeq
is a dimensionless moment of inertia, which may vary from zero to a maximum value of $2/3$, the former corresponding to a concentration of the mass at the center of the sphere, while the latter corresponds to a concentration of the mass on the surface of the sphere. The value $\kappa=2/5$ refers to a uniform distribution of the mass in the sphere.


The inelastic collisions of rough spherical particles are characterized by the relationships
 \begin{equation}
\hnn{\sigma}\cdot\overline{\mathbf{g}}^\prime=-\alpha(\hnn{\sigma}\cdot\overline{\mathbf{g}}),\quad
\hnn{\sigma}\times\overline{\mathbf{g}}^\prime=-\beta(\hnn{\sigma}\times\overline{\mathbf{g}}),
 \ee{4}
where $0<\alpha\leq1$ and $-1\leq\beta\leq1$ are the normal and tangential restitution coefficients, respectively. For an elastic collision of perfectly smooth spheres one has $\alpha=1$ and $\beta=-1$, while $\alpha=1$ and $\beta=1$ for an elastic encounter of perfectly rough spherical particles.

{}Insertion of Eq.\ \eqref{3bg} into Eq.\ \eqref{4} gives $\hnn{\sigma}\cdot\mathbf{Q}=m\widetilde\alpha(\hnn{\sigma}\cdot\overline{\mathbf{g}})$ and $\hnn{\sigma}\times\mathbf{Q}=m\widetilde\beta(\hnn{\sigma}\times\overline{\mathbf{g}})$, where the following abbreviations have been introduced:
 \beq\label{7}
 \widetilde\alpha\equiv{1+\alpha\over2},\quad \widetilde\beta\equiv{1+\beta\over2}{\kappa\over\kappa+1}.
 \eeq
Therefore, the impulse can be expressed as
 \bal
\mathbf{Q}=&{m\widetilde\alpha}(\hnn{\sigma}\cdot\overline{\mathbf{g}})\hnn{\sigma}-m\widetilde\beta\hnn{\sigma}\times(\hnn{\sigma}\times\overline{\mathbf{g}})\nonumber\\
=&{m\widetilde\alpha}(\hnn{\sigma}\cdot\mathbf{g})\hnn{\sigma}-m\widetilde\beta\hnn{\sigma}\times\left(\hnn{\sigma}\times\mathbf{g}+\sigma {\boldsymbol{\omega}+\boldsymbol{\omega}_1\over 2}\right).
\label{6}
 \eal


Equations \eqref{1,2} and (\ref{6})  express the  postcollisional velocities  in terms of the precollisional velocities and of the collision vector \cite{note_14_05_1}. {}From these results it is easy to obtain that the change of the total (translational plus rotational) kinetic energy reads
 \bal
 \Delta K=&{m\over2}\left(v^{\prime2}+v_1^{\prime2}-v^2-v_1^2\right)+{I\over2}\left(\omega^{\prime2}+\omega_1^{\prime2}-\omega^2-\omega_1^2\right)\nonumber\\
   =&-m{1-\beta^2\over4}{\kappa\over\kappa+1}\left[\hnn{\sigma}\times\left(\hnn{\sigma}\times\mathbf{g}+\sigma {\boldsymbol{\omega}+\boldsymbol{\omega}_1\over 2}\right)\right]^2\nonumber\\
   &-m{1-\alpha^2\over4}(\hnn{\sigma}\cdot\mathbf{g})^2.
  \label{Delta_E}
 \eal
The right-hand side  vanishes for elastic collisions of perfectly smooth spheres $(\alpha=1,\beta=-1)$ and for  elastic collisions of perfectly rough spherical particles  $(\alpha=1, \beta=1)$. In those cases the total energy is conserved in a collision.

Apart from the linear momentum, the angular momentum is conserved, namely,
\bal
m (\mathbf{r}\times \mathbf{v} +\mathbf{r}_1\times \mathbf{v}_1)+I(\boldsymbol{\omega}+\boldsymbol{\omega}_1)=&m (\mathbf{r}\times \mathbf{v}' +\mathbf{r}_1\times \mathbf{v}_1')\nonumber\\
&+I(\boldsymbol{\omega}'+\boldsymbol{\omega}_1'),
\label{bal_ang}
\eal
where $\mathbf{r}$ and $\mathbf{r}_1=\mathbf{r}+\sigma\hnn{\sigma}$ are the position vectors of the two colliding particles.

\subsection{Boltzmann equation}
A direct encounter is characterized by the precollisional velocities $(\mathbf{v}, \mathbf{v}_1; \boldsymbol{\omega},
 \boldsymbol{\omega}_1)$, by the postcollisional velocities $(\mathbf{v}^\prime, \mathbf{v}_1^\prime; \boldsymbol{\omega}^\prime,
 \boldsymbol{\omega}_1^\prime)$, and by the collision vector $\hnn{\sigma}$. For a restitution encounter
  the pre- and postcollisional velocities are denoted by $(\mathbf{v}^\ast, \mathbf{v}_1^\ast; \boldsymbol{\omega}^\ast,
 \boldsymbol{\omega}_1^\ast)$ and $(\mathbf{v}, \mathbf{v}_1; \boldsymbol{\omega}, \boldsymbol{\omega}_1)$, respectively,  and the collision vector by $\hnn{\sigma}^\ast=-\hnn{\sigma}$.
It is easy to verify  the relationship $\hnn{\sigma}^\ast\cdot\mathbf{g}=-\alpha(\hnn{\sigma}^\ast\cdot\mathbf{g}^\ast)=-\hnn{\sigma}\cdot\mathbf{g}$.
The modulus of the Jacobian of the  transformation $(\mathbf{v}^\ast, \mathbf{v}_1^\ast; \boldsymbol{\omega}^\ast,
 \boldsymbol{\omega}_1^\ast)\to(\mathbf{v}, \mathbf{v}_1; \boldsymbol{\omega}, \boldsymbol{\omega}_1)$ is given by
\beq
\left|{\partial(\mathbf{v}^\ast, \mathbf{v}_1^\ast; \boldsymbol{\omega}^\ast,
 \boldsymbol{\omega}_1^\ast)\over\partial(\mathbf{v}, \mathbf{v}_1; \boldsymbol{\omega}, \boldsymbol{\omega}_1)}\right|={1\over\alpha\beta^2}
 \label{13}
 \eeq
Thus,
 \begin{equation}
 (\hnn{\sigma}^\ast\cdot\mathbf{g}^\ast) d\mathbf{v}^\ast d\boldsymbol{\omega}^\ast d\mathbf{v}_1^\ast d\boldsymbol{\omega}_1^\ast={1\over\alpha^2\beta^2}
 (\hnn{\sigma}\cdot\mathbf{g})d\mathbf{v} d\boldsymbol{\omega} d\mathbf{v}_1 d\boldsymbol{\omega}_1.
 \ee{14}

From Eq.\ (\ref{14}) we may infer that the Boltzmann equation for granular gases of rough spherical particles without external forces and torques is given by
\begin{subequations}
  \label{15}
  \beq
{\partial f \over \partial t}+\mathbf{v}\cdot\boldsymbol{\nabla} f=J[f,f],
\label{15a}
\eeq
 \beq
J[f,f]=\sigma ^2\int d\mathbf{v}_1\int d\boldsymbol{\omega}_1\int_{+} d\hnn{\sigma}\,(\hnn{\sigma}\cdot\mathbf{g})\left({f_1^\ast f^\ast\over\alpha^2\beta^2} -f_1f\right),
\label{15b}
 \eeq
\end{subequations}
where  $f(\mathbf{r},\mathbf{v}, \boldsymbol{\omega}, t)$ is the one-particle distribution function in the phase space spanned by the positions and the linear  and angular
velocities of the particles.
As usual, in Eq.\ \eqref{15b} the notation $f_1=f(\mathbf{v}_1, \boldsymbol{\omega}_1)$, $f^*=f(\mathbf{v}^*, \boldsymbol{\omega}^*)$, $f_1^*=f(\mathbf{v}_1^*, \boldsymbol{\omega}_1^*)$ has been employed. Also, the subscript ($+$) in the integration over $\hnn{\sigma}$ denotes the constraint $\hnn{\sigma}\cdot\mathbf{g}>0$.

The so-called transfer equation  is obtained from the multiplication of the Boltzmann equation by an arbitrary function $\psi(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)$ and integration of the resulting equation over all values of the velocities $\mathbf{v}$ and $\boldsymbol{\omega}$, yielding
\beq
\label{16}
\partial_t \left(n\langle\psi\rangle\right) +\boldsymbol{\nabla}\cdot \left(n \langle \mathbf{v}\psi\rangle\right)-n\langle (\partial_t+\mathbf{v}\cdot\boldsymbol{\nabla})\psi\rangle=
\mathcal{J}[\psi|f,f],
\eeq
where
\beq
n(\mathbf{r},t)=\int d\mathbf{v}\int d\boldsymbol{\omega}\, f(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)
\label{n}
\eeq
is the local number density,
\beq
\langle\psi\rangle ={1\over n(\mathbf{r},t)}\int d\mathbf{v}\int d\boldsymbol{\omega}\, \psi(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)f(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)
\label{avpsi}
\eeq
is the local average value of $\psi$, and
\bal
\label{16J}
\mathcal{J}[\psi|f,f]\equiv&\int d\mathbf{v}\int d\boldsymbol{\omega}\,\psi(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)J[f,f]\nonumber\\
=& {\sigma ^2\over2}\int d\mathbf{v}\int d\boldsymbol{\omega}\int d\mathbf{v}_1\int d\boldsymbol{\omega}_1\int_+ d\hnn{\sigma}\,(\hnn{\sigma}\cdot\mathbf{g})\nonumber\\
&\times\left(\psi_1'+\psi'-\psi_1-\psi\right)f_1f
\eal
is the collisional production term of $\psi$.
In the second step  of Eq.\ \eqref{16J} we have used the relationship (\ref{14}) and the standard symmetry properties of the collision term.

 \subsection{Balance equations}

A macroscopic description of a rarefied granular gas of rough spheres can be characterized by the following basic fields: the particle number density
$n({\mathbf{r}},t)$ [see Eq.\ \eqref{n}], the hydrodynamic flow velocity
$\mathbf{u}({\mathbf{r}},t)$, and the granular temperature $T(\mathbf{r},t)$. The two latter quantities are defined as
\beq
\mathbf{u}=\langle \mathbf{v}\rangle, \quad T={1\over2}(T_t+T_r),
\label{u,T}
\eeq
where
\beq
T_t={m\over3}\langle V^2\rangle,\quad T_r={I\over3}\langle \omega^2\rangle
\label{Tt,Tr}
\eeq
are the (partial) translational and rotational temperatures, respectively, and the averages $\langle \cdots\rangle$ are defined by Eq.\ \eqref{avpsi}. In Eq.\ \eqref{Tt,Tr}, $\mathbf{V}=\mathbf{v}-\mathbf{u}$ is the (translational) peculiar velocity. On the other hand, in the definition of $T_r$ we have chosen not to refer the
angular velocities to the mean value
\beq
\boldsymbol{\Omega}=\langle \boldsymbol{\omega}\rangle
 \label{Omega}
 \eeq
because the latter is not a conserved quantity. Had we defined the granular temperature as $\overline{T}=(T_t+\overline{T}_r)/2$ with $\overline{T}_r={I\over3}\langle(\boldsymbol{\omega}- \boldsymbol{\Omega})^2\rangle$, then $\overline{T}$ would not be a conserved quantity in the case of completely rough and elastic collisions ($\alpha=\beta=1$), even though $\Delta K=0$ in that case [see Eq.\ \eqref{Delta_E}]

The balance equations for the basic fields are obtained from the transfer equation (\ref{16})  with the following choices for the arbitrary function $\psi(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)$:
\begin{enumerate}
\item
 Balance of particle number density ($\psi=1$),
\beq
\label{18}
{\cal D}_tn+n\boldsymbol{\nabla}\cdot \mathbf{u}=0.
\eeq
\item
Balance of momentum density ($\psi=m\mathbf{v}$),
\beq
\label{19}
\rho{\cal D}_t \mathbf{u}
+\boldsymbol{\nabla}\cdot \mathsf{P}
=0.
\eeq
\item
Balance of   temperature ($\psi=mV^2/2+I\omega^2/2$),
\beq
\label{20}
{\cal D}_t T +{1\over3n}\left(
\boldsymbol{\nabla}\cdot \mathbf{q}+\mathsf{P}:\boldsymbol{\nabla} \mathbf{u}\right)+ T  \zeta=0.
\eeq
\end{enumerate}
In the above balance equations, $\rho=mn$ is the mass density, $\mathcal{D}_t=\partial_t+\mathbf{u}\cdot\boldsymbol{\nabla}$ denotes the material time derivative, and the following quantities have been introduced: the pressure tensor
\beq
P_{ij}=\rho \langle V_iV_j\rangle,
\eeq
the heat flux vector
\beq
\mathbf{q}=\mathbf{q}_t+\mathbf{q}_r
\eeq
with
\beq
\mathbf{q}_t={\rho\over2}\langle V^2\mathbf{V}\rangle,\quad
\mathbf{q}_r={In\over2}\langle \omega^2\mathbf{V}\rangle,
\eeq
and the cooling rate
\beq
\label{zeta}
\zeta={T_t\over 2T}\zeta_t+{T_r\over 2T}\zeta_r
\eeq
with
\begin{subequations}
\label{ztrrot}
\beq
\label{ztr}
\zeta_t=-\frac{m}{3nT_t}\mathcal{J}[v^2|f,f],
\eeq
\beq
\label{zrot}
\zeta_r=-\frac{I}{3nT_r}\mathcal{J}[\omega^2|f,f].
\eeq
\end{subequations}


For further use, we note that the hydrostatic pressure $p$ is defined as one third of the trace of the pressure tensor, so that
\beq
p=nT_t.
\eeq
Also, the balance equations for the partial temperatures are
\beq
\label{20t}
{\cal D}_t T_t +{2\over3n}\left(
\boldsymbol{\nabla}\cdot \mathbf{q}_t+\mathsf{P}:\boldsymbol{\nabla} \mathbf{u}\right)+ T_t  \zeta_t=0,
\eeq
\beq
\label{20r}
{\cal D}_t T_r +{2\over3n}
\boldsymbol{\nabla}\cdot \mathbf{q}_r+ T_r  \zeta_r=0.
\eeq
Combination of Eqs.\ \eqref{20t} and \eqref{20r} yields Eq.\ \eqref{20}.


It is worth noting that the conservation of angular momentum, Eq.\ \eqref{bal_ang}, does not generate a balance equation with $\mathcal{J}[\psi|f,f]=0$ for the quantity $\psi= m \mathbf{r}\times\mathbf{v}+I\boldsymbol{\omega}$ because of the difference $\mathbf{r}_1-\mathbf{r}=\sigma\hnn{\sigma}$ between the centers of the two colliding spheres.

\section{Homogeneous cooling state}
\label{sec3}
Before considering the transport properties in inhomogeneous states, it is convenient to characterize the main properties of homogeneous states. In those cases, Eqs.\ \eqref{18} and \eqref{19} imply $n=\text{const}$ and $\mathbf{u}=\text{const}$, while Eqs.\ \eqref{20t} and \eqref{20r} become
\beq
\label{20tt}
\partial_tT_t + T_t  \zeta_t=0,\quad
\partial_t T_r + T_r  \zeta_r=0.
\eeq
The exact forms for the cooling rates $\zeta_t$ and $\zeta_r$ cannot be determined, unless the distribution function $f(\mathbf{v},\boldsymbol{\omega},t)$ is known. Good estimates of those quantities can be obtained by assuming the approximation
\beq
f(\mathbf{v},\boldsymbol{\omega})\to  \left({m \over 2\pi T_t}\right)^{3/2} e^{-mV^2/2T_t}f_r(\boldsymbol{\omega}),
\label{3.1}
\eeq
where
\beq
f_r(\boldsymbol{\omega})=\int d\mathbf{v}\, f(\mathbf{v},\boldsymbol{\omega})
\eeq
is the marginal distribution of angular velocities.
Equation \eqref{3.1} can be justified by maximum-entropy arguments, except that the explicit expression of $f_r(\boldsymbol{\omega})$ does not need to be specified. Substitution of \eqref{3.1} into Eqs.\ \eqref{ztrrot} and evaluation of the collision integrals given by Eq.\ \eqref{16J} yields \cite{SKG10}
\bal
\zeta_t=&\frac{5}{12}\Bigg[1-\alpha^2+\frac{\kappa}{1+\kappa}\left(1-\beta^2\right)-\frac{\kappa}{(1+\kappa)^2}\left(1+\beta\right)^2
\nonumber\\
&\times \theta\left(1+ X-{1\over\theta}\right)\Bigg]\nu,
\label{2.12}
\eal
\bal
\zeta_r=&\frac{5}{12}\frac{1+\beta}{1+\kappa}\Bigg[(1-\beta){(1+X)}+\frac{\kappa}{1+\kappa}\left(1+\beta\right)
\nonumber\\
&\times\left(1+X-{1\over\theta}\right)\Bigg]\nu,
\label{2.13}
\eal
where
\beq
\theta\equiv {T_r\over T_t},\quad X\equiv {I\Omega^2\over 3T_r},
\eeq
and
\beq
\nu\equiv\frac{16}{5}\sigma^2 n\sqrt{\pi T_t/m}
\label{nu}
\eeq
is an effective collision frequency. Note that $X=1-\overline{T}_r/T_r\leq 1$.
The total cooling rate is, according to Eq.\ \eqref{zeta},
\beq
\zeta=\frac{5}{12}\frac{1}{1+\theta}\left[1-\alpha^2+\frac{1-\beta^2}{1+\kappa}\theta\left(\frac{\kappa}{\theta}+1+X\right)\right]\nu.
\label{2.14}
\eeq




Since Eqs.\ \eqref{2.12} and \eqref{2.13} involve the norm $\Omega^2$ of the mean angular velocity, we need to complement Eq.\ \eqref{20tt} with the evolution equation for $\Omega^2$. By using again the approximation \eqref{3.1} one obtains \cite{SKG10}
\beq
\partial_t \Omega^2+2\zeta_\Omega \Omega^2=0,\quad \zeta_\Omega=\frac{5}{6}{1+\beta\over 1+\kappa}\nu.
\label{zetaO}
\eeq
By introducing the time variable $s(t)=\int_0^t dt'\,\nu(t')$, which measures the (nominal) number of collisions per particle from the initial time to time $t$, the solution to Eq.\ \eqref{zetaO} is
\beq
\Omega^2(s)=\Omega^2(0)e^{-2\zeta_\Omega^* s},\quad \zeta_\Omega^*\equiv \zeta_\Omega/\nu.
\label{Omega2}
\eeq
This allows us to solve the set of coupled equations \eqref{20tt} to obtain the time dependence of the temperatures $T_t(s)$ and $T_r(s)$. Both quantities asymptotically decrease in time due to energy dissipation. On the other hand, the relevant quantity is the temperature ratio $\theta(s)$. Analogously, rather than the time decay of $\Omega^2(s)$, and since $T_r(s)$ also tends to decay in time, the relevant quantity is the ratio $X(s)$. The evolution equations for both quantities are
\begin{subequations}
\beq
\partial_s \theta+(\zeta_r^*-\zeta_t^*)\theta=0,
\eeq
\beq
 \partial_s X+(2\zeta_\Omega^*-\zeta_r^*)X=0,
\eeq
\end{subequations}
where $\zeta_t^*\equiv\zeta_t/\nu$ and $\zeta_r^*\equiv\zeta_r/\nu$. Since
\bal
2\zeta_\Omega^*-\zeta_r^*=&\frac{5}{12}\frac{1+\beta}{(1+\kappa)^2}\Big[3-X+\beta(1+X)+2\kappa(1-X)\nonumber\\
&+\kappa\frac{1+\beta}{\theta}\Big]
\eal
is positive definite, it follows that $X\to 0$ monotonically, no matter the initial condition. On the other hand, the  evolution equation for $\theta$ admits a nonzero stationary solution given by the condition $\zeta_r^*=\zeta_t^*$. Such a solution is
\beq
\theta_\infty=\sqrt{1+h^2}+h
\label{1.8}
\eeq
with
\beq
h\equiv \frac{1+\kappa}{2\kappa(1+\beta)}\left[(1+\kappa)\frac{1-\alpha^2}{1+\beta}-(1-\kappa)(1-\beta)\right].
\label{1.9}
\eeq

\begin{figure}
\includegraphics[width=8cm]{nfig1.eps}
  \caption{{(Color online) Plot of the characteristic relaxation times $\ell_1^{-1}$ and $\ell_2^{-1}$ as functions of $\beta$ for $\kappa=\frac{2}{5}$ and $\alpha=0.8$}.}
  \label{nfig1}
\end{figure}

A standard linear stability analysis of the stationary solution $X=0$ and $\theta=\theta_\infty$ shows that the two associated eigenvalues are
\bal
\ell_{1}=&\theta_\infty\left.\frac{\partial (\zeta_r^*-\zeta_t^*)}{\partial\theta}\right|_{\theta=\theta_\infty,X=0}\nonumber\\
=&\frac{5}{12}\kappa\left(\frac{1+\beta}{1+\kappa}\right)^2\frac{1+\theta_\infty^2}{\theta_\infty},
\eal
\bal
\ell_{2}=&\left.{2\zeta_\Omega^*-\zeta_r^*}\right|_{\theta=\theta_\infty,X=0}\nonumber\\
=&\frac{5}{12}\frac{1+\beta}{(1+\kappa)^2}\left(3+\beta+2\kappa+\kappa\frac{1+\beta}{\theta_\infty}\right).
\eal
As expected, both eigenvalues are positive definite, what proves the stability of the stationary solution $(\theta,X)=(\theta_\infty,0)$  against homogeneous perturbations.
{The time evolution of $X(s)$ is governed by the eigenvalue $\ell_2$ only, so that the relaxation time (in units of number of collisions per particle) of $X(s)$ is $s\sim \ell_2^{-1}$. As for $\theta(s)$, its relaxation time is  $s\sim \ell_1^{-1}$ if $X(0)=0$ and $s\sim \max(\ell_1^{-1},\ell_2^{-1})$ if $X(0)\neq 0$. Figure \ref{nfig1} shows the dependence of both relaxation times on roughness for the representative case $\alpha=0.8$. Except in a  narrow roughness region adjacent to $\beta=-1$, one has $\ell_1^{-1}>\ell_2^{-1}$, so that $X(s)$ relaxes more rapidly than $\theta(s)$. }

\begin{figure}
\includegraphics[width=8cm]{fig1.eps}
  \caption{(Color online) Plot of (a) $X=I\Omega^2/3T_r$ and (b) $\theta=T_r/T_t$ as functions of the number of collisions per particle ($s$) scaled by $(1+\beta)^{-2}$ for $\kappa=\frac{2}{5}$, $\alpha=0.8$, and $\beta=-0.5$, $0$, $0.5$, and $1$, starting from the initial condition $\theta(0)=1$, $X(0)=\frac{1}{2}$. The horizontal lines in panel (b) correspond to the respective stationary values $\theta_\infty$.}
  \label{fig1}
\end{figure}

As an illustration of the evolution of $\theta(s)$ and $X(s)$, Fig.\ \ref{fig1} shows both quantities for $\alpha=0.8$ and $\beta=-0.5$, $0$, $0.5$, and $1$, starting from the initial condition $\theta(0)=1$, $X(0)=\frac{1}{2}$. It is quite apparent that $\theta$ reaches its stationary value \eqref{1.8} in about $30/(1+\beta)^2$ collisions per particle, while $X$ goes to zero in a significantly shorter period.


\begin{figure}
\includegraphics[width=8cm]{nfig2.eps}
  \caption{{(Color online) Plot of  $X=I\Omega^2/3T_r$ (dashed line) and  $(1+\beta)^2\theta=(1+\beta)^2T_r/T_t$ (solid line) as functions of the number of collisions per particle ($s$) scaled by $(1+\beta)^{-1}$ for $\kappa=\frac{2}{5}$, $\alpha=0.8$, and $\beta=-0.99$,  starting from the initial condition $\theta(0)=1$, $X(0)=\frac{1}{2}$. The dash-dotted line represents $(1+\beta)^2\theta$ if $X(0)=0$. The horizontal line corresponds to the  stationary value $(1+\beta)^2\theta_\infty$.}}
  \label{nfig2}
\end{figure}

{The special quasismooth limit $\beta\to -1$ deserves further comments \cite{S11b}. In that case, the asymptotic rotational-translational temperature ratio $\theta_\infty$ diverges as $\theta_\infty \to [(1+\kappa)^2/\kappa](1-\alpha^2)(1+\beta)^{-2}$, the two eigenvalues becoming $\ell_1\to \frac{5}{12}(1-\alpha^2)$, which is just the cooling rate of perfectly smooth spheres \cite{vNE98}, and $\ell_2\to \frac{5}{6}(1+\beta)/(1+\kappa)\to 0$. Therefore, if $X(0)=0$, $(1+\beta)^2\theta(s)$ relaxes to  $(1+\beta)^2\theta_\infty$ after a \emph{finite} number of collisions per particle on the order of $\ell_1^{-1}$. On the other hand, if $X(0)\neq 0$, $(1+\beta)^2\theta(s)$ evolves in two well-defined stages. The first stage lasts a characteristic time $\sim \ell_1^{-1}$ and is very similar to that of the case $X(0)=0$. This first stage is followed by a much slower relaxation (with $s\sim \ell_2^{-1}\sim (1+\beta)^{-1}\to\infty$ collisions per particle) toward the asymptotic value. This singular scenario in the quasismooth limit is illustrated by Fig.\ \ref{nfig2} for $\beta=-0.99$ and $\alpha=0.8$.}

\begin{figure}
\includegraphics[width=8cm]{nfig3a.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{nfig3b.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{nfig3c_lr.eps}
\caption{{(Color online) (a) Reduced cooling rate $\zeta^*$ as a function of $\beta$ for $\kappa=\frac{2}{5}$ and $\alpha=0.6$, $0.8$, and $1$. (b) Reduced cooling rate $\zeta^*$ as a function of $\alpha$ for $\kappa=\frac{2}{5}$ and $\beta=-0.5$, $0$, $0.5$, and $1$, as well as in the quasismooth limit $\beta\to -1$ and for purely smooth spheres ($\beta=-1$). (c) Density plot of the reduced cooling rate $\zeta^*$ for $\kappa=\frac{2}{5}$. The contour lines correspond to $\zeta^*=0.025,0.05,\ldots,0.35$.}}
\label{nfig3}
\end{figure}




The stationary solution \eqref{1.8} represents the value of the temperature ratio in the HCS \cite{HZ97,LHMZ98,HHZ00,AHZ01,Z06}. In such a state the whole time dependence of the distribution function only occurs through the granular temperature $T(t)$, so that the Boltzmann equation \eqref{15a} becomes
\beq
-\zeta T \partial_T f =J[f,f],
\label{B_HCS}
\eeq
{where, according to Eq.\ \eqref{2.14}, $\zeta=\zeta^*\nu$ with}
\beq
{\zeta^*=\frac{5}{12}\frac{1}{1+\theta_\infty}\left[1-\alpha^2+({1-\beta^2})
\frac{\theta_\infty+\kappa}{1+\kappa}\right].}
\label{zetaHCS}
\eeq
{The dependence of the HCS (reduced) cooling rate $\zeta^*$ on both $\alpha$ and $\beta$ is displayed in Fig.\ \ref{nfig3}. As can be observed, at a given value of $\alpha$ a maximum of $\zeta^*$ occurs around $\beta\approx 0$. In the region close to $\beta=-1$, where $\theta_\infty\gg 1$, Eq.\ \eqref{zetaHCS} becomes $\zeta^*\approx \frac{5}{12}(1-\beta^2)/(1+\kappa)$. Therefore, $\lim_{\beta\to -1}\zeta^*=0$.}

The HCS distribution function has the scaling form
\beq
f(\mathbf{v},\boldsymbol{\omega},t)=n\left(\frac{mI}{\tau_t\tau_r}\right)^{3/2} \left[T(t)\right]^{-3}\phi\left(\mathbf{c}(t),\mathbf{w}(t)\right),
\label{phi}
\eeq
where
\beq
\mathbf{c}(t)={\mathbf{V}\over\sqrt{2\tau_t T(t)/m}},\quad \mathbf{w}(t)={\boldsymbol{\omega}\over\sqrt{2\tau_r T(t)/m}},
\label{cw}
\eeq
with the time-independent temperature ratios
\beq
\tau_t\equiv {T_t(t)\over T(t)}={2\over1+\theta_\infty},\quad \tau_r\equiv {T_r(t)\over T(t)}={2\theta_\infty\over1+\theta_\infty}.
\eeq
Hence, according to Eq.\ \eqref{phi}, one has the relation
\beq
 T \partial_T f
=-\frac{1}{2}\left(\frac{\partial }{\partial\mathbf{V}}\cdot \mathbf{V} +\frac{\partial }{\partial\boldsymbol{\omega}}\cdot \boldsymbol{\omega}\right)f.
\eeq

{In the particular case of perfectly elastic and smooth particles ($\alpha=1$, $\beta= -1$), the rotational velocities must be ignored (i.e., $\tau_r=0$, $\tau_t=2$) and the solution of Eq.\ \eqref{B_HCS} is just the equilibrium distribution of translational velocities. In the other conservative case of perfectly elastic and rough spheres ($\alpha=\beta= 1$),  the solution of Eq.\ \eqref{B_HCS} is the equilibrium distribution with a common temperature (i.e., $\tau_r=\tau_t=1$).}
In the general dissipative case, however,  the  solution of Eq.\ \eqref{B_HCS} is not exactly known, although good approximations are given in the form of a two-temperature Maxwellian multiplied by  truncated Sonine polynomial expansions \cite{HZ97,ML98,LHMZ98,HHZ00,AHZ01,CLH02,Z06,BPKZ07,KBPZ09,BB12,SKS11,VSK14}.

\section{Chapman--Enskog method}
\label{sec4}

\subsection{Outline of the method}
To determine the distribution function from the Chapman--Enskog method \cite{CC70,FK72}, we write the Boltzmann equation (\ref{15a}) as
\begin{equation}
\mathcal{D}_tf+\epsilon \mathbf{V}\cdot\boldsymbol{\nabla} f =J[f,f],
\ee{22}
where  $\epsilon$ is a uniformity parameter (set equal to unity at the end of the calculations) measuring the strength of the spatial gradients.
According to the method, the distribution function and the material time derivative are expanded in terms of the parameter $\epsilon$ as follows,
\begin{equation}
f=f^{(0)}+\epsilon f^{(1)}+\epsilon^2 f^{(2)}+\cdots,
\ee{23a}
\begin{equation}
\mathcal{D}_t=\mathcal{D}_t^{(0)}+\epsilon \mathcal{D}_t^{(1)}+\epsilon^2 \mathcal{D}_t^{(2)}+\cdots.
\ee{23b}
Substitution of Eq.\ \eqref{23a} into the definitions of the fluxes and the cooling rate gives
\beq
P_{ij}=p^{(0)}\delta_{ij}+\epsilon P_{ij}^{(1)}+\epsilon^2 P_{ij}^{(2)}+\cdots,
\label{P0}
\eeq
\beq
\mathbf{q}=\epsilon \mathbf{q}^{(1)}+\epsilon^2 \mathbf{q}^{(2)}+\cdots,
\eeq
\begin{equation}
\zeta=\zeta^{(0)}+\epsilon \zeta^{(1)}+\epsilon^2 \zeta^{(2)}+\cdots,
\ee{23c}
where $p^{(0)}=n\tau_t  T$,
\beq
\zeta^{(0)}=-\frac{1}{6nT}\int d\mathbf{v}\int d\boldsymbol{\omega}\,\left({mV^2}+{I\omega^2}\right)J[f^{(0)},f^{(0)}],
\eeq
\bal
\label{zeta1}
\zeta^{(1)}=&\frac{1}{6nT}\int d\mathbf{v}\int d\boldsymbol{\omega}\,\left({mV^2}+{I\omega^2}\right)\mathcal{L}f^{(1)}\nonumber\\
=&\frac{\tau_t}{2}\zeta_t^{(1)}+\frac{\tau_r}{2}\zeta_r^{(1)}.
\eal
In Eq.\ \eqref{zeta1} $\mathcal{L}$ is the linearized collision operator defined as
\beq
\mathcal{L}\Phi=-J[\Phi,f^{(0)}]-J[f^{(0)},\Phi].
\eeq
Thus, $\zeta^{(1)}=\mathcal{J}_L[mV^2+I\omega^2|f^{(1)}]/6nT$, with the notation
\bal
\label{16L}
\mathcal{J}_L[\psi|\Phi]\equiv&\int d\mathbf{v}\int d\boldsymbol{\omega}\,\psi(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)\mathcal{L}\Phi\nonumber\\
=& {\sigma ^2\over2}\int d\mathbf{v}\int d\boldsymbol{\omega}\int d\mathbf{v}_1\int d\boldsymbol{\omega}_1\int_+ d\hnn{\sigma}\,(\hnn{\sigma}\cdot\mathbf{g})\nonumber\\
&\times\left(\psi_1+\psi-\psi_1'-\psi'\right)(f_1^{(0)}\Phi+f^{(0)}\Phi_1),
\eal
where $\psi_1\equiv \psi(\mathbf{v}_1,\boldsymbol{\omega}_1)$, $\psi'\equiv \psi(\mathbf{v}',\boldsymbol{\omega}')$, $\psi_1'\equiv \psi(\mathbf{v}_1',\boldsymbol{\omega}_1')$, $f_1^{(0)}\equiv f^{(0)}(\mathbf{v}_1,\boldsymbol{\omega}_1)$, $\Phi_1\equiv \Phi(\mathbf{v}_1,\boldsymbol{\omega}_1)$, and in the last step use has been made of Eq.\ \eqref{16J}.








Insertion of the expansions (\ref{23a}) and \eqref{23b} into the Boltzmann equation (\ref{22}) leads to the corresponding integro-differential equations for the different orders $f^{(k)}$. In particular, the two first equations are
\beq\label{24a}
  {\cal D}_t^{(0)} f^{(0)}=J[f^{(0)},f^{(0)}],
  \eeq
\beq
 \label{24b}
 \left({\cal D}_t^{(0)} +\mathcal{L}\right)f^{(1)}=-\left({\cal D}_t^{(1)}+\mathbf{V}\cdot\boldsymbol{\nabla} \right)f^{(0)}.
 \eeq

Since the distribution function in \eqref{23a} depends on time and space only through its dependence on the hydrodynamic fields $n$, $\mathbf{u}$, and $T$,
the action of the operator $\mathcal{D}_t^{(k)}$ can be written as
\beq
\mathcal{D}_t^{(k)}=\left(\mathcal{D}_t^{(k)} n\right)\frac{\partial}{\partial n}+ \left(\mathcal{D}_t^{(k)} \mathbf{u}\right)\cdot\frac{\partial}{\partial \mathbf{u}}+\left(\mathcal{D}_t^{(k)} T\right)\frac{\partial}{\partial T}.
\label{Dtk}
\eeq
{}From the balance equations \eqref{18}--\eqref{20} it follows that
\beq\label{26a}
{\cal D}_t^{(0)}n=0,\quad {\cal D}_t^{(1)}n=-n\boldsymbol{\nabla}\cdot\mathbf{u},
\eeq
\beq
\label{26b}
{\cal D}_t^{(0)}\mathbf{u}=0,\quad {\cal D}_t^{(1)}\mathbf{u}=-{\tau_t\over\rho}\left(n\boldsymbol{\nabla} T+T\boldsymbol{\nabla} n\right),
\eeq
\beq
\label{26c}
{\cal D}_t^{(0)} T =- T  \zeta^{(0)},\quad{\cal D}_t^{(1)} T =- T  \zeta^{(1)}-{ \tau_t\over3}T\boldsymbol{\nabla}\cdot\mathbf{u}.
\eeq

Taking into account that ${\cal D}_t^{(0)}\Phi=-\zeta^{(0)}T\partial_T\Phi$, it is obvious that Eq.\ \eqref{24a} is formally equivalent to the HCS equation \eqref{B_HCS}. This means that $f^{(0)}$ is the \emph{local} version of the HCS distribution. Moreover, in the approximation \eqref{3.1}, {$\zeta^*=\zeta^{(0)}/\nu$ is given by Eq.\ \eqref{zetaHCS}}.

\subsection{First-order distribution}

The first-order function $f^{(1)}$ obeys the linear equation \eqref{24b}. By using the properties \eqref{Dtk}--\eqref{26c}, the inhomogeneous term of Eq.\ \eqref{24b} becomes
\bal
-\left({\cal D}_t^{(1)}+\mathbf{V}\cdot\boldsymbol{\nabla}\right)f^{(0)}=&\mathbf{A}\cdot\boldsymbol{\nabla} \ln T+\mathbf{B}\cdot\boldsymbol{\nabla} \ln n{+C_{ij}\nabla_j u_i}\nonumber\\
&+E\boldsymbol{\nabla}\cdot\mathbf{u}+T\zeta^{(1)}\partial_T f^{(0)},
\label{ABCE}
\eal
where
\beq
\mathbf{A}=-T\left(\mathbf{V}\partial_T+\frac{\tau_t}{m}\partial_{\mathbf{V}}\right)f^{(0)},
\label{eqA}
\eeq
\beq
\mathbf{B}=-\left(\mathbf{V}+\frac{\tau_t T}{m}\partial_{\mathbf{V}}\right)f^{(0)},
\label{eqB}
\eeq
\beq
{C_{ij}=\left(\partial_{V_i}V_j-\frac{1}{3}\delta_{ij}\partial_{\mathbf{V}}\cdot\mathbf{V}\right)f^{(0)}},
\label{eqC}
\eeq
\beq
E=\frac{1}{3}\left(\partial_\mathbf{V}\cdot\mathbf{V}+\tau_tT\partial_T\right)f^{(0)}.
\label{eqE}
\eeq


In the case of pure smooth particles, $f^{(0)}$ is a function of the translational velocity only, there is no rotational energy  and hence $\tau_r=0$, $\tau_t=2$, and $T\partial_T f^{(0)}=-\frac{1}{2}\partial_{\mathbf{V}}\cdot\mathbf{V} f^{(0)}$. As a consequence, $E=0$ and one recovers the known results for inelastic smooth  particles \cite{BDKS98}. Thus, the presence of roughness induces a non-vanishing function $E$, even in the conservative case of perfectly rough particles ($\alpha=\beta=1$) \cite{CC70,K10a}. A subtler consequence of roughness is the symmetry breakdown of the traceless tensor $C_{ij}$. Isotropy implies that $f^{(0)}(\mathbf{V},\boldsymbol{\omega})$ is a function of the three scalars $V^2$, $\omega^2$ and $\chi^2\equiv(\mathbf{V}\cdot\boldsymbol{\omega})^2$. Therefore,
\beq
{C_{ij}- C_{ji}=2\frac{\partial f^{(0)}}{\partial \chi^2}(\mathbf{V}\cdot\boldsymbol{\omega})(V_j\omega_i-V_i\omega_j)}.
\eeq
If one neglects the orientational correlations between $\mathbf{V}$ and $\boldsymbol{\omega}$ in the HCS so that the dependence of $f^{(0)}$ on $\chi^2$ is ignored, then $C_{ij}=C_{ji}$, as happens in the pure smooth case.




Taking into account Eq.\ \eqref{ABCE}, the solution to Eq.\ \eqref{24b} has the form
\beq
f^{(1)}=
\boldsymbol{\mathcal{A}}\cdot\boldsymbol{\nabla} \ln T+\boldsymbol{\mathcal{B}}\cdot\boldsymbol{\nabla} \ln n+\mathcal{C}_{ij}{\nabla_j u_i}
+\mathcal{E}\boldsymbol{\nabla}\cdot\mathbf{u},
\label{f1}
\eeq
where the vectors $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$, the traceless tensor $\mathcal{C}_{ij}$, and the scalar $\mathcal{E}$ are unknown functions to be determined.
Combination of Eqs.\ \eqref{zeta1} and \eqref{f1} allows us to write the first-order contribution to the cooling rate as
\beq
\zeta^{(1)}={-\xi}\boldsymbol{\nabla}\cdot\mathbf{u}
\label{zeta1xi}
\eeq
with
\beq
\xi={\frac{1}{2}\left(\tau_t\xi_t+\tau_r\xi_r\right)},
\label{xi}
\eeq
where
\beq
{\xi_t=-\frac{m}{3n\tau_tT}\mathcal{J}_L[V^2|\mathcal{E}],\quad \xi_r=-\frac{I}{3n\tau_rT}\mathcal{J}_L[\omega^2|\mathcal{E}]}
\label{xitr}
\eeq
are the translational and rotational contributions to $\xi$.

Substitution of Eq.\ \eqref{f1} into Eq.\ \eqref{24b} gives the following set of linear integral equations:
\beq
\left(-\frac{\zeta^{(0)}}{2}-\zeta^{(0)}T\partial_T+\mathcal{L}\right)\boldsymbol{\mathcal{A}}=\mathbf{A},
\label{eqLA}
\eeq
\beq
\left(-\zeta^{(0)}T\partial_T+\mathcal{L}\right)\boldsymbol{\mathcal{B}}=\mathbf{B}+\zeta^{(0)}\boldsymbol{\mathcal{A}},
\label{eqLB}
\eeq
\beq
\left(-\zeta^{(0)}T\partial_T+\mathcal{L}\right)\mathcal{C}_{ij}=C_{ij},
\label{eqLC}
\eeq
\beq
\left(-\zeta^{(0)}T\partial_T+\mathcal{L}\right)\mathcal{E}{+\xi} T\partial_T f^{(0)}=E.
\label{eqLE}
\eeq
In Eqs.\ \eqref{eqLA} and \eqref{eqLB} use has been made of the property
\bal
\mathcal{D}_t^{(0)}\boldsymbol{\nabla} \ln T=&\boldsymbol{\nabla} \mathcal{D}_t^{(0)}\ln T=-\boldsymbol{\nabla}\zeta^{(0)}\nonumber\\
=&-\zeta^{(0)}\left(\boldsymbol{\nabla} \ln n+\frac{1}{2}\boldsymbol{\nabla} \ln T\right).
\eal

According to Eqs.\ \eqref{eqA}--\eqref{eqE}, the functions $\mathbf{A}$, $\mathbf{B}$, $C_{ij}$ and $E$ are orthogonal to $\{1,\mathbf{V}, mV^2+I\omega^2,\boldsymbol{\omega}\}$, i.e.,
\beq
\int d\mathbf{v}\int d\boldsymbol{\omega}\, \left(
\begin{array}{c}
  1\\
  \mathbf{V}\\
  mV^2+I\omega^2\\
  \boldsymbol{\omega}
  \end{array}
  \right)
  \cdot
  \left(\mathbf{A},\mathbf{B},C_{ij},E\right)=0.
  \eeq
Therefore, the Fredholm alternative \cite{DS67} implies that the necessary conditions for the existence of solutions (solubility conditions) to Eqs.\ \eqref{eqLA}--\eqref{eqLE}  are
\beq
\int d\mathbf{v}\int d\boldsymbol{\omega}\, \left(
\begin{array}{c}
  1\\
  \mathbf{V}\\
  mV^2+I\omega^2\\
  \boldsymbol{\omega}
  \end{array}
  \right)
  \cdot
  \left(\boldsymbol{\mathcal{A}},\boldsymbol{\mathcal{B}},\mathcal{C}_{ij},\mathcal{E}\right)=0.
  \label{solub}
  \eeq
The solubility conditions associated with $\{1,\mathbf{V}, mV^2+I\omega^2\}$ mean that, by construction, the hydrodynamic quantities $n$, $\mathbf{u}$, and $T$ are fully contained in $f^{(0)}$. As for the solubility condition associated with $\boldsymbol{\omega}$, it implies
\beq
\boldsymbol{\Omega}^{(1)}=\frac{1}{n}\int d\mathbf{v}\int d\boldsymbol{\omega}\, \boldsymbol{\omega} f^{(1)}=0.
\eeq
  \section{Navier--Stokes--Fourier transport coefficients}
  \label{sec5}
\subsection{Exact formal expressions}



{}From Eq.\ \eqref{f1} one can express the first-order pressure tensor and heat flux as
\beq
P_{ij}^{(1)}=-\eta \left(\nabla_i u_j+\nabla_j u_i-\frac{2}{3}\delta_{ij}\boldsymbol{\nabla}\cdot\mathbf{u} \right)-\eta_b\delta_{ij}\boldsymbol{\nabla}\cdot \mathbf{u},
\label{Pij1}
\eeq
\beq
\mathbf{q}^{(1)}=-\lambda \boldsymbol{\nabla} T-\mu\boldsymbol{\nabla} n,
\label{q1}
\eeq
where the transport coefficients are
\beq
\eta=-\frac{m}{10}\int d\mathbf{v}\int d\boldsymbol{\omega}\, \left(V_iV_j-\frac{V^2}{3}\delta_{ij}\right)\mathcal{C}_{ij},
\eeq
\beq
\eta_b=-\frac{m}{3}\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathcal{E},
\eeq
\beq
\lambda={\tau_t\lambda_t+\tau_r\lambda_r},\quad \mu={\mu_t+\mu_r},
\eeq
with
\beq
\lambda_t=-\frac{m}{6\tau_tT}\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathbf{V}\cdot\boldsymbol{\mathcal{A}},
\eeq
\beq
 \lambda_r=-\frac{I}{6\tau_rT}\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^2\mathbf{V}\cdot\boldsymbol{\mathcal{A}},
\eeq
\beq
\mu_t=-\frac{m}{6n}\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathbf{V}\cdot\boldsymbol{\mathcal{B}},
\eeq
\beq
\mu_r=-\frac{I}{6n}\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^2\mathbf{V}\cdot\boldsymbol{\mathcal{B}}.
\eeq
In the constitutive equations \eqref{Pij1} and \eqref{q1}, $\eta$ is the shear viscosity, $\eta_b$ is the bulk viscosity, $\lambda$ is the thermal conductivity, and $\mu$ is a Dufour-like coefficient. The two latter coefficients have translational ($\lambda_t$, $\mu_t$) and rotational ($\lambda_r$, $\mu_r$) contributions.

By multiplying Eq.\ \eqref{eqLA} by  $V^2\mathbf{V}$ and $\omega^2\mathbf{V}$, and integrating over velocity one obtains
\beq
{\lambda_t=\frac{5n\tau_t T}{2m}\frac{1+2 a_{20}^{(0)}}{\nu_{\lambda_t}-2\zeta^{(0)}}},
\label{lambdat}
\eeq
\beq
{\lambda_r=\frac{3n\tau_t T}{2m}\frac{1+2a_{11}^{(0)}}{\nu_{\lambda_r}-2\zeta^{(0)}}},
\label{lambdar}
\eeq
where
\beq
{\nu_{\lambda_t}=\frac{\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathbf{V}\cdot\mathcal{L}\boldsymbol{\mathcal{A}}}{\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathbf{V}\cdot\boldsymbol{\mathcal{A}}}},
\label{nulambdat}
\eeq
\beq
{\nu_{\lambda_r}=\frac{\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^2\mathbf{V}\cdot\mathcal{L}\boldsymbol{\mathcal{A}}}{\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^2\mathbf{V}\cdot\boldsymbol{\mathcal{A}}}}
\label{nulambdar}
\eeq
are the associated collision frequencies and
\beq
{a_{20}^{(0)}=\frac{m^2}{15n\tau_t^2 T^2}\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^4 f^{(0)}-1},
\eeq
\beq
{a_{11}^{(0)}=\frac{m I}{9n\tau_t\tau_r T^2}\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\omega^2 f^{(0)}-1},
\eeq
are cumulants of the HCS distribution $f^{(0)}$ \cite{VSK14}. Upon deriving Eqs.\ \eqref{lambdat} and \eqref{lambdar} we have taken into account that, by dimensional analysis, $T\partial_T\lambda_t=\frac{1}{2}\lambda_t$ and $T\partial_T\lambda_r=\frac{1}{2}\lambda_r$. Analogously, from Eq.\ \eqref{eqLB} one gets
\beq
{\mu_t=\frac{\tau_t T}{n}\frac{\lambda_t\zeta^{(0)}+(5n\tau_t T/2m)a_{20}^{(0)}}{\nu_{\mu_t}-\frac{3}{2}\zeta^{(0)}}},
\label{mut}
\eeq
\beq
{\mu_r=\frac{\tau_r T}{n}\frac{\lambda_r\zeta^{(0)}+(3n\tau_t T/2m)a_{11}^{(0)}}{\nu_{\mu_r}-\frac{3}{2}\zeta^{(0)}}},
\label{mur}
\eeq
where
\beq
{\nu_{\mu_t}=\frac{\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathbf{V}\cdot\mathcal{L}\boldsymbol{\mathcal{B}}}{\int d\mathbf{v}\int d\boldsymbol{\omega}\, V^2\mathbf{V}\cdot\boldsymbol{\mathcal{B}}}},
\label{numut}
\eeq
\beq
{\nu_{\mu_r}=\frac{\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^2\mathbf{V}\cdot\mathcal{L}\boldsymbol{\mathcal{B}}}{\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^2\mathbf{V}\cdot\boldsymbol{\mathcal{B}}}},
\label{numur}
\eeq
and use has been made of $T\partial_T\mu_t=\frac{3}{2}\mu_t$ and $T\partial_T\mu_r=\frac{3}{2}\mu_r$.

Next, multiplication of Eq.\ \eqref{eqC} by $V_iV_j-\frac{1}{3}\delta_{ij}V^2$ and integration over velocity yields
\beq
\eta=\frac{n\tau_t T}{\nu_\eta-\frac{1}{2}\zeta^{(0)}},
\label{eta}
\eeq
where
\beq
\nu_\eta=\frac{\int d\mathbf{v}\int d\boldsymbol{\omega}\, \left(V_iV_j-\frac{V^2}{3}\delta_{ij}\right)\mathcal{L}\mathcal{C}_{ij}}{\int d\mathbf{v}\int d\boldsymbol{\omega}\, \left(V_iV_j-\frac{V^2}{3}\delta_{ij}\right)\mathcal{C}_{ij}}.
\label{nueta}
\eeq
Finally, multiplying Eq.\ \eqref{eqLE} by $V^2$ allows us to obtain
\beq
\eta_b=\frac{\tau_t\tau_r nT}{\zeta^{(0)}}{\left(\xi_t-\xi_r-\frac{2}{3}\right)}.
\label{etab}
\eeq
In Eqs.\ \eqref{eta} and \eqref{etab} we have made use of the properties $T\partial_T\eta=\frac{1}{2}\eta$ and $T\partial_T\eta_b=\frac{1}{2}\eta_b$.



The existence of a nonzero bulk viscosity induces a breakdown of energy equipartition additional to the one already present in the HCS.
Taking the trace in Eqs.\ \eqref{P0} and \eqref{Pij1}, one has
\beq
T_t=\tau_t T-\frac{\eta_b}{n}\boldsymbol{\nabla}\cdot\mathbf{u}+\cdots,
\eeq
where the ellipses denote terms of at least second order in the hydrodynamic gradients. Since the total temperature $T$ is not affected by the gradients, then
\beq
 T_r=\tau_r T+\frac{\eta_b}{n}\boldsymbol{\nabla}\cdot\mathbf{u}+\cdots.
\eeq
Thus the temperature ratio becomes
\beq
\frac{T_r}{T_t}=\frac{\tau_r}{\tau_t}\left[1+2\frac{\xi_t-\xi_r-\frac{2}{3}}{\zeta^{(0)}}\boldsymbol{\nabla}\cdot\mathbf{u}+\cdots\right],
\eeq
where use has been made of Eq.\ \eqref{etab}.

Equations \eqref{lambdat}, \eqref{lambdar}, \eqref{mut}, \eqref{mur}, \eqref{eta}, and \eqref{etab} are formally exact but require the solution of the set of linear integral equations \eqref{eqLA}--\eqref{eqLE}. As happens in the conservative case \cite{CC70,FK72}, the exact solution of those equations is not known. Since $\boldsymbol{\mathcal{A}}$ is a vector, it can be expressed as a sum of projections along the three \emph{polar} vectors $\mathbf{V}$, $(\mathbf{V}\cdot\boldsymbol{\omega})\boldsymbol{\omega}$, and $\mathbf{V}\times\boldsymbol{\omega}$ \cite{CLD65}, namely
\beq
\boldsymbol{\mathcal{A}}=\mathcal{A}_1 \mathbf{V}+\mathcal{A}_2(\mathbf{V}\cdot\boldsymbol{\omega})\boldsymbol{\omega}+\mathcal{A}_3\mathbf{V}\times\boldsymbol{\omega},
\eeq
where $\mathcal{A}_i$ are unknown isotropic scalar functions, i.e., they depend on $V^2$, $\omega^2$, and $\chi^2=(\mathbf{V}\cdot\boldsymbol{\omega})^2$ only.
Of course, the vector function $\boldsymbol{\mathcal{B}}$ has a similar structure. The tensor $\mathcal{C}_{ij}$ can be expressed as a combination of traceless dyadic products of the three vectors $\mathbf{V}$, $(\mathbf{V}\cdot\boldsymbol{\omega})\boldsymbol{\omega}$, and $\mathbf{V}\times\boldsymbol{\omega}$ with unknown scalar coefficients. Finally, $\mathcal{E}$ is an unknown scalar function.

In Sec.\ \ref{sec5B} we derive explicit expressions for the transport coefficients by considering the leading terms in a Sonine polynomial expansion.

\subsection{Sonine approximation\label{sec5B}}

As said in Sec.\ \ref{sec3}, the HCS distribution function $f^{(0)}$ is not exactly known. In a recent paper \cite{VSK14}, the first four relevant cumulants (in particular, $a_{20}^{(0)}$ and $a_{11}^{(0)}$) have been studied theoretically by means of a Sonine polynomial expansion and also computationally by means of the direct simulation Monte Carlo (DSMC) method \cite{B94} and event-driven molecular dynamics. The results show that three of the cumulants are in general relatively small. On the other hand, the cumulant
\beq
{a_{02}^{(0)}=\frac{I^2}{15n\tau_r^2 T^2}\int d\mathbf{v}\int d\boldsymbol{\omega}\, \omega^4 f^{(0)}-1}
\eeq
measuring the kurtosis of the angular velocity distribution can take values larger than unity in the region of small roughness. Outside that region, $a_{02}^{(0)}$  decreases as roughness increases. As an example (see Fig.\ 3 of Ref.\ \cite{VSK14}), at $\alpha=0.9$, $|a_{20}^{(0)}|\lesssim 0.02$ and $|a_{11}^{(0)}|\lesssim 0.07$ for  the whole range $-1\leq \beta\leq 1$, and $a_{02}^{(0)}\lesssim 0.2$ in the range $-0.5\lesssim \beta\leq 1$.


{}From a practical point of view, it must be noted that most of the
materials are characterized by positive values of the roughness
parameter (typically, $\beta\sim 0.5$) \cite{L99} and in those cases the cumulants are small. As a consequence, the HCS distribution can be rather well approximated by the \emph{two-temperature} Maxwellian
\beq
f^{(0)}\to f^{(0)}_M=n\Bigg({m I\over 4\pi^2\tau_t\tau_r T^2}\Bigg)^{3\over2}e^{-c^2-w^2},
\label{Maxw}
\eeq
where we recall that the scaled translational and angular velocities $\mathbf{c}$ and $\mathbf{w}$ are defined by Eq.\ \eqref{cw}. In this Maxwellian approximation, the functions \eqref{eqA}--\eqref{eqE} reduce to
\beq
\mathbf{A}\to- f^{(0)}_Mv_0\left(c^2+w^2-4\right)\mathbf{c},
\label{A_M}
\eeq
\beq
\mathbf{B}\to 0,
\eeq
\beq
C_{ij}\to -2f^{(0)}_M\left(c_ic_j-\frac{1}{3}c^2\delta_{ij}\right),
\eeq
\beq
E\to\frac{1}{3}f^{(0)}_M\left[\tau_t\left(w^2-\frac{3}{2}\right)-\tau_r\left(c^2-\frac{3}{2}\right)\right],
\label{E_M}
\eeq
where
\beq
v_0\equiv \sqrt{2\tau_t T/m}
\eeq
is the (translational) thermal speed.

The Maxwellian forms \eqref{A_M}--\eqref{E_M} suggest to approximate the unknown functions $\boldsymbol{\mathcal{A}}$, $\boldsymbol{\mathcal{B}}$, $\mathcal{C}_{ij}$, and $\mathcal{E}$ by
\beq
\boldsymbol{\mathcal{A}}\to -f^{(0)}_M \frac{v_0}{\nu}\left[\gamma_{A_t}\left(c^2-\frac{5}{2}\right)+\gamma_{A_r}\left(w^2-\frac{3}{2}\right)\right]\mathbf{c},
\label{Aapprox}
\eeq
\beq
\boldsymbol{\mathcal{B}}\to  -f^{(0)}_M\frac{v_0}{\nu}\left[\gamma_{B_t}\left(c^2-\frac{5}{2}\right)+\gamma_{B_r}\left(w^2-\frac{3}{2}\right)\right]\mathbf{c},
\eeq
\beq
\mathcal{C}_{ij}\to -f^{(0)}_M\frac{\gamma_C}{\nu} \left(c_ic_j-\frac{1}{3}c^2\delta_{ij}\right),
\eeq
\beq
\mathcal{E}\to f^{(0)}_M\frac{\gamma_E}{\nu}\left[\tau_t\left(w^2-\frac{3}{2}\right)-\tau_r\left(c^2-\frac{3}{2}\right)\right],
\label{Eapprox}
\eeq
where $\nu$ is defined by Eq.\ \eqref{nu} (with $T_t=\tau_t T$) and the   $\gamma$ coefficients are directly related to the transport coefficients by
\beq
\eta=\frac{n\tau_t T}{\nu}\frac{\gamma_C}{2},\quad
\eta_b=\frac{n\tau_t\tau_r T}{\nu}\gamma_E,
\label{X1}
\eeq
\beq
\lambda_t=\frac{5}{2}\frac{n\tau_t T}{m\nu}\gamma_{A_t},\quad \lambda_r=\frac{3}{2}\frac{n\tau_t T}{m\nu}\gamma_{A_r},
\label{X2}
\eeq
\beq
\mu_t=\frac{5}{2}\frac{\tau_t^2 T^2}{m\nu}\gamma_{B_t},\quad \mu_r=\frac{3}{2}\frac{\tau_t \tau_r T^2}{m\nu}\gamma_{B_r}.
\label{X3}
\eeq
It can be checked that the forms \eqref{Aapprox}--\eqref{Eapprox} are consistent with the solubility conditions \eqref{solub}.




The basic Sonine approximations \eqref{Aapprox}--\eqref{Eapprox} allow us to evaluate explicitly the transport coefficients. The main steps are described in the Appendix and the final results are displayed in Table \ref{table1} \cite{note_14_05_2}.


\begin{figure}
\includegraphics[width=8cm]{fig2a.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig2b.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig2c_lr.eps}
\caption{(Color online) (a) Reduced shear viscosity $\eta^*$ as a function of $\beta$ for $\kappa=\frac{2}{5}$ and $\alpha=0.6$, $0.8$, and $1$. (b) Reduced shear viscosity $\eta^*$ as a function of $\alpha$ for $\kappa=\frac{2}{5}$ and $\beta=-0.5$, $0$, $0.5$, and $1$, as well as in the quasismooth limit $\beta\to -1$ and for purely smooth spheres ($\beta=-1$). (c) Density plot of the reduced shear viscosity $\eta^*$ for $\kappa=\frac{2}{5}$. The contour lines correspond to $\eta^*=1,1.05,\ldots,1.4$.}
\label{fig2}
\end{figure}


\begin{figure}
\includegraphics[width=8cm]{fig3a.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig3b.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig3c_lr.eps}
\caption{(Color online) Same as in Fig.\ \protect\ref{fig2} but for the reduced bulk viscosity $\eta_b^*$. The contour lines  in panel (c) correspond to $\eta_b^*=0.5,0.75,\ldots,3.5$.}
\label{fig3}
\end{figure}

\begin{figure}
\includegraphics[width=8cm]{fig4a.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig4b.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig4c_lr.eps}
\caption{(Color online) Same as in Fig.\ \protect\ref{fig2} but for the reduced thermal conductivity $\lambda^*$. The contour lines  in panel (c) correspond to $\lambda^*=1,1.25,\ldots,3.75$.}
\label{fig4}
\end{figure}

\begin{figure}
\includegraphics[width=8cm]{fig5a.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig5b.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig5c_lr.eps}
\caption{(Color online) Same as in Fig.\ \protect\ref{fig2} but for the reduced Dufour-like coefficient $\mu^*$. The contour lines  in panel (c) correspond to $\mu^*=0.25,0.5,\ldots,3.5$.}
\label{fig5}
\end{figure}

\begin{figure}
\includegraphics[width=8cm]{fig6a.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig6b.eps}\\
\vskip0.2cm
\includegraphics[width=8cm]{fig6c_lr.eps}
\caption{(Color online) Same as in Fig.\ \protect\ref{fig2} but for the cooling rate transport coefficient $\xi$. The contour lines  in panel (c) correspond to $\xi=0,0.05,\ldots,0.4$.}
\label{fig6}
\end{figure}


\section{Discussion}
\label{sec6}
Let us now analyze the dependence of the five transport coefficients on both $\alpha$ and $\beta$. In order to define dimensionless quantities, we will take as a reference the transport coefficients (shear viscosity and thermal conductivity) of a gas made of elastic and smooth spheres at the same translational temperature $\tau_t T$ as that of the HCS of the granular gas, i.e.,
\beq
\eta_0=\frac{n\tau_t T}{\nu},\quad \lambda_0=\frac{15}{4}\frac{\eta_0}{m}.
\eeq
More specifically, the dimensionless shear and bulk viscosities are
\beq
\eta^*=\frac{\eta}{\eta_0},\quad \eta_b^*=\frac{\eta_b}{\eta_0},
\eeq
while the dimensionless thermal conductivity and Dufour-like coefficients are
\beq
\lambda^*=\frac{\lambda}{\lambda_0},\quad \mu^*=\frac{\mu}{\lambda_0}\frac{n}{ T}.
\eeq


\subsection{Limiting cases}
Before considering the case of general $\alpha$ and $\beta$, it is worth considering some limiting cases. We start with the case of a granular gas constituted by purely smooth inelastic hard spheres ($\beta=-1$ with arbitrary $\alpha$). In such a gas, the rotational degrees of freedom are irrelevant and should not contribute to the transport properties. A convenient way of isolating the relevant translational properties consists in formally setting $\theta_\hcs=0$ (i.e., $\tau_r=0$ and $\tau_t=2$) in the general expressions of Table \ref{table1},  apart from taking $\beta=-1$.
Alternatively, one can take $\theta_\hcs=0$ and $\kappa=0$, the latter defining spheres with a vanishing moment of inertia.
In both cases one obtains the results presented in Table \ref{table2}, which are consistent  with those previously derived for smooth inelastic hard spheres \cite{BDKS98,note_14_04_1}.



{\renewcommand{\arraystretch}{1.2}
\begin{table*}
   \caption{Special limiting cases.}\label{table2}
\begin{ruledtabular}
\begin{tabular}{cccc}
Quantity&Pure smooth&Quasismooth limit&Perfectly rough and elastic  \\
&(${\beta }=-1$)&(${\beta }\to -1$)&($\alpha=\beta=1$) \\ \hline
$\eta^*$&$\displaystyle{\frac{24}{(1+{\alpha })(13-{\alpha })}}$&$\displaystyle{\frac{24}{(1+{\alpha })(19-7{\alpha })}}$&$\displaystyle{\frac{6(1+\kappa)^2}{6+13\kappa}}$\\
$\eta_b^*$&$0$&$\displaystyle{\frac{8}{5(1-{\alpha }^2)}}$&$\displaystyle{\frac{(1+\kappa)^2}{10\kappa}}$\\
$\lambda^*$&$\displaystyle{\frac{64}{(1+{\alpha })(9+7{\alpha })}}$&$\displaystyle{\frac{48}{25(1+{\alpha })}}$&$\displaystyle{\frac{12 (1+\kappa)^2 \left(37+151 \kappa +50 \kappa ^2\right)}{25 \left(12+75 \kappa +101 \kappa ^2+102 \kappa^3\right)}}$\\
$\mu^*$&$\displaystyle{\frac{1280(1-{\alpha })}{(1+{\alpha })(9+7{\alpha })(19-3{\alpha })}}$&$0$&$0$\\
$\xi$&$0$&$0$&$0$\\
 \end{tabular}
 \end{ruledtabular}
 \end{table*}
}




Next, we consider the quasismooth limit $\beta\to -1$. The results differ from those obtained before for pure smooth spheres because, {as shown in Sec.\ \ref{sec3}},  the HCS rotational-translational temperature ratio diverges as $\theta_\hcs\sim (1+\beta)^{-2}$, and thus the rotational contributions cannot be neglected. The final expressions, which are independent of the reduced moment of inertia $\kappa$, are also included in Table \ref{table2}. In contrast, it can be checked that the expressions in the limit of small moment of inertia ($\kappa\to 0$) depend on $\beta$.

Finally, let us analyze the physically important case of perfectly elastic and rough spheres ($\alpha=\beta=1$). This defines a conservative system (i.e., the total kinetic energy is conserved by collisions) that has been used for a long time to model polyatomic gases \cite{P22,MSD66,CC70,K10a}. The results obtained in this limit are displayed in the last column of Table \ref{table2} and  fully agree with those first derived by Pidduck \cite{P22}.

\subsection{General $\alpha$ and $\beta$}

Now we go back to a granular gas with general values of $\alpha$, $\beta$, and $\kappa$, in which case the expressions for the five transport coefficients are given in Table \ref{table1}. For the sake of concreteness, let us restrict ourselves to spheres with a uniform mass distribution, so that $\kappa=\frac{2}{5}$.

Figures \ref{fig2}--\ref{fig6} show the dependence of the reduced transport coefficients $\eta^*$, $\eta_b^*$, $\lambda^*$, $\mu^*$, and $\xi$, respectively, on both $\alpha$ and $\beta$. {As in Fig.\ \ref{nfig3},} in the top panels the quantities are plotted versus $\beta$ for three representative values of the coefficient of normal restitution ($\alpha=0.6$, $0.8$, and $1$). The middle panels present the dependence on $\alpha$ for a few representative values of the coefficient of tangential restitution, namely $\beta=-0.5$, $0$, $0.5$, and $1$. Additionally, the quasismooth limit ($\beta\to -1$) and the case of purely smooth spheres ($\beta=-1$) are also considered. Finally, the bottom panels represent density plots in the $\alpha$-$\beta$ plane.

Let us start analyzing the three transport coefficients that are also present in the case of purely smooth particles, i.e.,  the (reduced) shear viscosity, thermal conductivity, and Dufour-like coefficient (Figs.\ \ref{fig2}, \ref{fig4}, and \ref{fig5}, respectively).  We observe that, at fixed $\alpha$, those coefficients present a non-monotonic $\beta$ dependence with  maxima around $\beta\approx 0$. On the other hand,  the dependence on $\alpha$ is rather sensitive to the value of $\beta$, showing an intricate interplay between both coefficients of restitution. Typically, the transport coefficients increase with increasing inelasticity, although some exceptions are found (see, for instance, $\eta^*$ at $\beta=1$ and $\lambda^*$ and $\mu^*$ at $\beta=-0.5$). Furthermore, an interesting observation from Figs.\ \ref{fig2}, \ref{fig4}, and \ref{fig5} is that the impact of the coefficient of normal restitution on the transport coefficients $\eta^*$, $\lambda^*$, and $\mu^*$ is much milder in the case of rough spheres than for purely smooth spheres.
{Comparison between Fig.\ \ref{nfig3}, on the one hand, and Figs.\ \ref{fig2}, \ref{fig4}, and \ref{fig5}, on the other hand, shows that the general dependencies of the transport coefficients $\eta^*$, $\lambda^*$, and $\mu^*$ on both $\alpha$ and $\beta$ (in particular, the maxima around $\beta\sim 0$) are highly correlated to that of the cooling rate $\zeta^*$. This explains that  fair \emph{qualitative} estimates can be obtained by using the expressions for smooth spheres \cite{BDKS98,GD99} with the cooling rate replaced by the one for rough spheres \cite{MDHEH13}. }

Next, we consider the two transport coefficients that vanish for purely smooth particles. As observed from Fig.\ \ref{fig3}, the bulk viscosity $\eta_b^*$ exhibits a highly nontrivial behavior. It reaches especially high values in the quasielastic and quasismooth region, diverging in the limit $\alpha\to 1$, $\beta\to -1$. Outside that region, the bulk viscosity can be larger than the shear viscosity. For instance, $\eta_b^*/\eta^*\simeq 1.26$ at $\alpha=0.7$ and $\beta=-0.2$. As for the cooling rate transport coefficient $\xi$ (see Fig.\ \ref{fig6}), it also presents a complex behavior.  A remarkable feature is that it becomes negative in a certain region of the $\alpha$-$\beta$ plane near $\alpha=1$. {Of course, this does not mean that the cooling rate itself is negative or signals any breakdown of the Chapman--Enskog method, the Sonine approximation, or the friction model used. According to Eqs.\ \eqref{23c} and \eqref{zeta1xi}, a negative value of the transport coefficient $\xi$ simply implies that   the cooling rate $\zeta$ is larger (smaller) than its HCS value $\zeta^{(0)}$ if $\boldsymbol{\nabla}\cdot \mathbf{u}$ is positive (negative).}






\section{Concluding remarks}
\label{sec7}
In this work we have developed a hydrodynamic theory for a dilute granular gas modeled as a system of identical inelastic and rough hard spheres. Energy dissipation in collisions is characterized by two constant coefficients of restitution: the normal ($0<\alpha\leq 1$) and the tangential ($-1\leq\beta\leq 1$) coefficients. In this model both the translational velocity ($\mathbf{v}$) of the center of mass and the angular velocity ($\boldsymbol{\omega}$) of the particles are mutually influenced by collisions.

The methodology has been based on the Boltzmann kinetic equation for the one-particle velocity distribution function $f(\mathbf{r},\mathbf{v},\boldsymbol{\omega},t)$. The kinetic equation has been solved by means of the Chapman--Enskog method \cite{CC70,FK72} for states with small spatial gradients of the hydrodynamic fields: the number density $n(\mathbf{r},t)$, the flow velocity $\mathbf{u}(\mathbf{r},t)$, and the (total) granular temperature $T(\mathbf{r},t)$. The solution provides the constitutive equations for the pressure tensor $P_{ij}(\mathbf{r},t)$, the heat flux $\mathbf{q}(\mathbf{r},t)$, and the cooling rate $\zeta(\mathbf{r},t)$. The associated five transport coefficients (shear viscosity $\eta$, bulk viscosity $\eta_b$, thermal conductivity $\lambda$, Dufour-like coefficient $\mu$, and cooling rate transport coefficient $\xi$) are exactly expressed in terms of integrals involving the solutions of a set of linear integral equations. In particular, the existence of $\eta_b$ and $\xi$ imply that, in the case of a compressible flow (i.e., $\boldsymbol{\nabla}\cdot\mathbf{u}\neq 0$), the rotational-translational temperature ratio and the cooling rate differ from their forms in the reference HCS.

As happens in the conventional case of elastic particles \cite{CC70}, explicit expressions for the transport coefficients can be obtained by expanding the zeroth- and first-order distributions in Sonine polynomials and truncating the expansions at the simplest level (the so-called first Sonine approximation). This has allowed us to determine the five transport coefficients as nonlinear functions of $\alpha$, $\beta$, and the reduced moment of inertia $\kappa$. For easy reference, the final results are displayed in Table \ref{table1}.

Our results extend to arbitrary values of $\alpha$ and $\beta$ previous works for inelastic and purely smooth spheres ($\alpha<1$, $\beta=-1$) \cite{BDKS98} and elastic and perfectly rough spheres ($\alpha=\beta=1$) \cite{P22,CC70,K10a}, as shown in Table \ref{table2}. Due to the coupling between translational and rotational degrees of freedom in the HCS, the quasismooth limit $\beta\to -1$ yields results differing from those for purely smooth spheres.

As Figs.\ \ref{fig2}--\ref{fig6} clearly show, the dependence of the transport coefficients on both $\alpha$ and $\beta$ is rather intricate.
On the other hand, it must be noted that, since the deviations of the reference HCS from the two-temperature Maxwellian \eqref{Maxw} are important in the region $-1\leq\beta\lesssim -0.5$ only \cite{VSK14},  the expressions derived here are expected to be especially reliable in the region $-0.5\lesssim \beta\leq 1$, which is likely the one of practical interest from an experimental point of view \cite{L99}.
{A more rigorous hydrodynamic theory in the region of small roughness ($\beta\gtrsim -1$) would likely require, apart from accounting for strong non-Maxwellian features of the HCS distribution, the inclusion of the mean spin $\boldsymbol{\Omega}$ as an additional hydrodynamic variable.}

The present work opens new challenges to explore. First, we plan to carry out a linear stability analysis \cite{G05} of the NSF hydrodynamic equations  to determine the critical length $L_c$ beyond which the HCS becomes unstable and assess the impact of roughness on $L_c$ \cite{MDHEH13}. Given that most of the experimental setups consider granular systems confined in two dimensions, we  intend to determine the NSF transport coefficients for systems of inelastic rough hard disks by using a methodology similar to the one followed here. Moreover, the structure of the collisional frequencies derived in  the Appendix can be exploited to obtain the NSF transport coefficients of driven granular gases, in analogy to the case of smooth spheres \cite{GM02,GMT13,GCV13}. Finally, we will test the transport coefficients obtained from the first Sonine approximation against DSMC numerical solutions  of the Boltzmann equation  by methods similar to those employed for smooth spheres \cite{GM02,MG03b,GM03,BR04,MSG05,BRMG05,MSG07}.





\begin{acknowledgments}
The research of A.S. and V.G.  was supported by the Spanish Government through Grant No.\ FIS2013-42840-P and  by the Junta de Extremadura (Spain) through Grant No.\ GR10158, both partially financed by FEDER funds.
The work of G.M.K. has been supported by the Conselho Nacional de Desenvolvimento Cient\'ifico e Tecnol\'ogico (Brazil).
\end{acknowledgments}


