2. Boltzmann kinetic equation for a granular gas surrounded by a molecular gas

We consider a gas of inelastic hard disks ($d=2$) or spheres ($d=3$) of mass $m$ and diameter $\sigma$. The spheres are assumed to be perfectly smooth, so that collisions between any two particles of the granular gas are characterised by a (positive) constant coefficient of normal restitution $\alpha \leq 1$. When $\alpha =1$ ($\alpha <1$), the collisions are elastic (inelastic). The granular gas is immersed in a gas of elastic hard disks or spheres of mass $m_g$ and diameter $\sigma _g$ (‘molecular gas’). A collision between a granular particle and a particle of the molecular gas is considered to be elastic. As discussed in § 1, we are interested here in describing a situation where the granular gas is sufficiently rarefied (the number density of granular particles is much smaller than that of the molecular gas) so that the state of the molecular gas is not affected by the presence of solid (grains) particles. In this sense, the background (molecular) gas may be treated as a thermostat, which is at equilibrium at a temperature $T_g$. Thus, the velocity distribution function $f_g$ of the molecular gas is the Maxwell–Boltzmann distribution:

(2.1)\begin{equation} f_g(\boldsymbol{V}_g)=n_g \left(\frac{m_g}{2{\rm \pi} T_g}\right)^{d/2} \exp \bigg(-\frac{m_g V_g^{2}}{2T_g}\bigg), \end{equation}

where $n_g$ is the number density of the molecular gas and $\boldsymbol {V}_g=\boldsymbol {v}-\boldsymbol {U}_g$, in which $\boldsymbol {U}_g$ is the mean flow velocity of the molecular gas. Note that here, for the sake of generality, we have assumed that $\boldsymbol {U}_g \neq \boldsymbol {U}$ ($\boldsymbol {U}$ being the mean flow velocity of the granular gas; see its definition in (2.7)). In addition, for the sake of simplicity, the Boltzmann constant $k_{B}=1$ throughout the paper.

In the low-density regime, the time evolution of the one-particle velocity distribution function $f(\boldsymbol {r}, \boldsymbol {v}, t)$ of the granular gas is given by the Boltzmann kinetic equation. Since the granular particles collide among themselves and with the particles of the molecular gas, in the absence of external forces the velocity distribution $f(\boldsymbol {r}, \boldsymbol {v}, t)$ verifies the kinetic equation:

(2.2)\begin{equation} \frac{\partial f}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f=J[\boldsymbol{v}\,|\,f,f]+J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

Here, the Boltzmann collision operator $J[\boldsymbol {v}\,|\,f,f]$ gives the rate of change of the distribution $f(\boldsymbol {r}, \boldsymbol {v}, t)$ due to binary inelastic collisions between granular particles. On the other hand, the Boltzmann–Lorentz operator $J_g[\boldsymbol {v}\,|\,f,f_g]$ accounts for the rate of change of the distribution $f(\boldsymbol {r}, \boldsymbol {v}, t)$ due to elastic collisions between granular and molecular gas particles.

The explicit form of the nonlinear Boltzmann collision operator $J[\boldsymbol {v}\,|\,f,f]$ is (Garzó Reference Garzó2019)

(2.3)\begin{align} J[\boldsymbol{v}_1\,|\,f,f]=\sigma^{d-1}\int {\rm d}\boldsymbol{v}_{2} \int {\rm d}\widehat{\boldsymbol {\sigma }}\varTheta (\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12}) (\widehat{\boldsymbol {\sigma }} \boldsymbol{\cdot} \boldsymbol{g}_{12}) [\alpha^{{-}2}f(\boldsymbol{v}_{1}'')f(\boldsymbol{v}_{2}'')-f(\boldsymbol{v}_{1})f(\boldsymbol{v}_{2})], \end{align}

where $\boldsymbol {g}_{12}=\boldsymbol {v}_1-\boldsymbol {v}_2$ is the relative velocity, $\widehat {\boldsymbol {\sigma }}$ is a unit vector along the line of centres of the two spheres at contact and $\varTheta$ is the Heaviside step function. In (2.3), the double primes denote pre-collisional velocities. The relationship between pre-collisional $(\boldsymbol {v}_1'', \boldsymbol {v}_2'')$ and post-collisional $(\boldsymbol {v}_1, \boldsymbol {v}_2)$ velocities is

(2.4a,b)\begin{equation} \boldsymbol{v}_1''=\boldsymbol{v}_1-\frac{1+\alpha}{2\alpha}(\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }},\quad \boldsymbol{v}_2''=\boldsymbol{v}_2+\frac{1+\alpha}{2\alpha}(\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }}. \end{equation}

The form of the linear Boltzmann–Lorentz collision operator $J_g[\boldsymbol {v}\,|\,f,f_g]$ is (Résibois & de Leener Reference Résibois and de Leener1977; Garzó Reference Garzó2019)

(2.5)\begin{align} J_g[\boldsymbol{v}_1\,|\,f,f_g]=\overline{\sigma}^{d-1}\int {\rm d}\boldsymbol{v}_{2} \int {\rm d}\widehat{\boldsymbol {\sigma }}\varTheta (\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12}) (\widehat{\boldsymbol {\sigma }}\boldsymbol{\cdot} \boldsymbol{g}_{12})[f(\boldsymbol{v}_{1}'')f_g(\boldsymbol{v}_{2}'')-f(\boldsymbol{v}_{1})f_g(\boldsymbol{v}_{2})], \end{align}

where $\overline{\sigma }=(\sigma +\sigma _g)/2$. In (2.5), the relationship between $(\boldsymbol {v}_1'',\boldsymbol {v}_2'')$ and $(\boldsymbol {v}_1,\boldsymbol {v}_2)$ is

(2.6a,b)\begin{equation} \boldsymbol{v}_1''=\boldsymbol{v}_1-2 \mu_g(\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }},\quad \boldsymbol{v}_2''=\boldsymbol{v}_2+2 \mu (\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }}, \end{equation}

where $\mu _g=m_g/(m+m_g)$ and $\mu =m/(m+m_g)$.

The relevant hydrodynamic fields of the granular gas are the number density $n(\boldsymbol {r}; t)$, the mean flow velocity $\boldsymbol {U}(\boldsymbol {r}; t)$ and the granular temperature $T(\boldsymbol {r}; t)$. They are defined, respectively, as

(2.7)\begin{equation} \left\{n, n \boldsymbol{U}, dn T\right\}=\int {\rm d}\boldsymbol{v} \{1, \boldsymbol{v}, m V^{2}\} f(\boldsymbol{v}), \end{equation}

where $\boldsymbol {V}=\boldsymbol {v}-\boldsymbol {U}$ is the peculiar velocity. As said before, in general the mean flow velocity $\boldsymbol {U}$ of solid particles is different from the mean flow velocity $\boldsymbol {U}_g$ of molecular gas particles. As we show later, the difference $\boldsymbol {U}-\boldsymbol {U}_g$ induces a non-vanishing contribution to the heat flux.

The macroscopic balance equations for the granular gas are obtained by multiplying (2.2) by $\{1, \boldsymbol {v}, m V^{2}\}$ and integrating over velocity. The result is

(2.8)$$\begin{gather} D_t n+n\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0, \end{gather}$$

(2.9)$$\begin{gather}\rho D_t\boldsymbol{U}={-}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}+\boldsymbol{\mathcal{F}}[f], \end{gather}$$

(2.10)$$\begin{gather}D_tT+\frac{2}{dn}\left(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q}+\boldsymbol{\mathsf{P}}:\boldsymbol{\nabla}\boldsymbol{U}\right)={-} T \zeta-T \zeta_g. \end{gather}$$

Here, $D_t=\partial _t+\boldsymbol {U}\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the material derivative, $\rho =m n$ is the mass density of solid particles and the pressure tensor $\boldsymbol{\mathsf{P}}$ and the heat flux vector $\boldsymbol {q}$ are given, respectively, as

(2.11a,b)\begin{equation} \boldsymbol{\mathsf{P}}=\int {\rm d}\boldsymbol{v} m \boldsymbol{V}\boldsymbol{V} f(\boldsymbol{v}),\quad \boldsymbol{q}=\int {\rm d}\boldsymbol{v} \frac{m}{2} V^{2} \boldsymbol{V} f(\boldsymbol{v}). \end{equation}

Since the Boltzmann–Lorentz collision term $J_g[\boldsymbol {v}\,|\,f,f_g]$ does not conserve momentum, then the production of momentum $\boldsymbol {\mathcal {F}}[f]$ is in general different from zero. It is defined as

(2.12)\begin{equation} \boldsymbol{\mathcal{F}}[f]=\int {\rm d}\boldsymbol{v} m \boldsymbol{V} J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

In addition, the partial production rates $\zeta$ and $\zeta _g$ are given, respectively, as

(2.13a,b)\begin{equation} \zeta={-}\frac{m}{d n T}\int {\rm d}\boldsymbol{v} V^{2} J[\boldsymbol{v}\,|\,f,f], \quad \zeta_g={-}\frac{m}{d n T}\int {\rm d}\boldsymbol{v} V^{2} J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

The cooling rate $\zeta$ gives the rate of kinetic energy loss due to inelastic collisions between particles of the granular gas. It vanishes for elastic collisions. The term $\zeta _g$ gives the transfer of kinetic energy between the particles of the granular and molecular gases. It vanishes when the granular and molecular gases are at the same temperature ($T_g=T$).

The macroscopic size of grains entails that their gravitational potential energy is much larger than the usual thermal energy scale $k_{B}T$. For example, a grain of common sand at room temperature ($T=300$ K) would require a energy of the order of $10^{5}k_{B}T$ to rise a distance equal to its diameter when subjected to the action of gravity (Heinrich, Nagel & Behringer Reference Heinrich, Nagel and Behringer1996). Therefore, thermal fluctuations have a negligible effect on the dynamics of grains, and so they are considered athermal systems. On the other hand, when particles are subjected to a strong excitation (e.g. vibrating walls or air-fluidised beds), the external energy supplied to the system can compensate for the energy dissipated by collisions and the effects of gravity. In this situation (rapid-flow conditions), particles’ velocities acquire some kind of random motion that looks much like the motion of the atoms or molecules in an ordinary or molecular gas. The rapid-flow regime opens up the possibility of establishing a relation between particles’ response to the external supply of energy and some kind of temperature. In this context, the granular temperature $T$ can be interpreted as a statistical quantity measuring the deviations (or fluctuations) of the velocities of grains with respect to its mean value $\boldsymbol {U}$. Just as in the ordinary case, the velocity fluctuations approach is the basic assumption in the construction of out-of-equilibrium theories such as the kinetic theory. Since granular gases are athermal, the granular temperature $T$ has no relation to the conventional thermodynamic temperature associated with the second thermodynamic law. On the other hand, although the classical ensemble averages provide a thermodynamic interpretation for the temperature $T_g$ of the bath (or molecular gas) through the definition of entropy, statistical mechanics shows that the thermodynamic temperature is the same as that obtained using the velocity fluctuations approach (up to a factor including the mass of the particles and the Boltzmann constant) (see for instance the review paper of Goldhirsch (Reference Goldhirsch2008) for a more detailed discussion of this issue). Moreover, coming back to the suspensions framework, the fluctuation–dissipation theorem supports the athermal statistical interpretation of $T_g$ since its value has been demonstrated to coincide with an effective temperature derived from the Einstein relation in both experiments and simulations (Puglisi, Baldassarri & Loreto Reference Puglisi, Baldassarri and Loreto2002; Garzó Reference Garzó2004; Chen & Hou Reference Chen and Hou2014). Thus, $T_g$ is treated also here in the same way as $T$. Namely, $T_g$ (defined through the distribution $f_g(\boldsymbol {v})$) is used to assess the deviations of the velocities of the molecular gas with respect to the mean value $\boldsymbol {U}_g$. In addition, we also assume that the values of $m$, $n$, $T$ and $T_g$ are such that a description based on classical mechanics is appropriate.

2.1. Brownian limit ($m/m_g \to \infty$)

The Boltzmann equation (2.2) applies in principle for arbitrary values of the mass ratio $m/m_g$. On the other hand, a physically interesting situation arises in the so-called Brownian limit, namely when the granular particles are much heavier than the particles of the surrounding molecular gas ($m/m_g \to \infty$). In this case, a Kramers–Moyal expansion (Résibois & de Leener Reference Résibois and de Leener1977; Rodríguez, Salinas-Rodríguez & Dufty Reference Rodríguez, Salinas-Rodríguez and Dufty1983; McLennan Reference McLennan1989) in the velocity jumps $\delta \boldsymbol {v}=(2/(1+m/m_g)) (\widehat {\boldsymbol {\sigma }}\boldsymbol {\cdot } \boldsymbol {g}_{12}) \boldsymbol {g}_{12}$ allows us to approximate the Boltzmann–Lorentz operator $J_g[\boldsymbol {v}\,|\,f,f_g]$ by the Fokker–Planck operator $J_g^{{FP}}[\boldsymbol {v}\,|\,f,f_g]$ (Résibois & de Leener Reference Résibois and de Leener1977; Rodríguez et al. Reference Rodríguez, Salinas-Rodríguez and Dufty1983; McLennan Reference McLennan1989; Brey, Dufty & Santos Reference Brey, Dufty and Santos1999a; Sarracino et al. Reference Sarracino, Villamaina, Costantini and Puglisi2010):

(2.14)\begin{equation} J_g[\boldsymbol{v}\,|\,f,f_g]\to J_g^{{FP}}[\boldsymbol{v}\,|\,f,f_g]=\gamma \frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot} \left(\boldsymbol{v}+\frac{T_g}{m} \frac{\partial}{\partial \boldsymbol{v}}\right)f(\boldsymbol{v}), \end{equation}

where the drift or friction coefficient $\gamma$ is defined as

(2.15)\begin{equation} \gamma=\frac{4{\rm \pi}^{(d-1)/2}}{d\varGamma\left(\dfrac{d}{2}\right)} \biggl(\frac{m_g}{m}\biggr)^{1/2}\left(\frac{2T_g}{m}\right)^{1/2}n_g \overline{\sigma}^{d-1}. \end{equation}

While obtaining (2.14) and (2.15), it has been assumed that $\boldsymbol {U}_g=\boldsymbol {0}$ and that the distribution $f(\boldsymbol {v})$ of the granular gas is a Maxwellian distribution.

Most of the suspension models employed in the granular literature to fully account for the influence of an interstitial molecular fluid on the dynamics of grains are based on the replacement of $J_g[\boldsymbol{v}\,|\,f,f_g]$ by the Fokker–Planck operator (2.14) (Koch & Hill Reference Koch and Hill2001). More specifically, for general inhomogeneous states, the impact of the background molecular gas on solid particles is through an effective force composed of three different terms: (i) a term proportional to the difference ${\rm \Delta} \boldsymbol {U}=\boldsymbol {U}-\boldsymbol {U}_g$, (ii) a drag force term mimicking the friction of grains on the viscous interstitial gas and (iii) a stochastic Langevin-like term accounting for the energy gained by grains due to their interactions with particles of the molecular gas (neighbouring particles effect) (Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012). This yields the following kinetic equation for gas–solid suspensions:

(2.16)\begin{equation} \frac{\partial f}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f-\gamma{\rm \Delta} \boldsymbol{U}\boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{v}}-\gamma \frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot} \boldsymbol{V} f-\gamma \frac{T_{g}}{m}\frac{\partial^{2} f}{\partial v^{2}}=J[\boldsymbol{v}\,|\,f,f]. \end{equation}

Note that there are three different scalars $(\beta,\gamma,\xi )$ in the suspension model proposed by Garzó et al. (Reference Garzó, Tenneti, Subramaniam and Hrenya2012); each one of the coefficients is associated with the different terms of the fluid–solid force. For the sake of simplicity, the results derived by Gómez González & Garzó (Reference Gómez González and Garzó2019) were obtained by assuming that $\beta =\gamma =\xi$.

Since the model attempts to mimic gas–solid flows where $\boldsymbol {U}\neq \boldsymbol {U}_g$, note that one has to make the replacement $\boldsymbol {v}\to \boldsymbol {v}-\boldsymbol {U}_g$ in (2.14) to obtain (2.16). The Boltzmann equation (2.16) has been solved by means of the Chapman–Enskog method (Chapman & Cowling Reference Chapman and Cowling1970) to first order in spatial gradients. Explicit forms for the Navier–Stokes–Fourier transport coefficients have been obtained in steady-state conditions, namely when the cooling terms are compensated for by the energy gained by the solid particles due to their collisions with the bath particles (Garzó, Chamorro & Vega Reyes Reference Garzó, Chamorro and Vega Reyes2013; Gómez González & Garzó Reference Gómez González and Garzó2019). Thus, the results derived in the present paper must be consistent with those previously obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) when the limit $m/m_g \to \infty$ is considered in our general results.

3. Homogeneous steady state

As a first step and before studying inhomogeneous states, we consider the HSS. The HSS is the reference base state (zeroth-order approximation) used in the Chapman–Enskog perturbation method (Chapman & Cowling Reference Chapman and Cowling1970). Therefore, its investigation is of great importance. The HSS was widely analysed by Santos (Reference Santos2003) for a three-dimensional granular gas. Here, we extend these calculations to a general dimension $d$.

In the HSS, the density $n$ and temperature $T$ are spatially uniform, and with an appropriate selection of the frame reference, the mean flow velocities vanish ($\boldsymbol {U}=\boldsymbol {U}_g=\boldsymbol {0}$). Consequently, the Boltzmann equation (2.2) reads

(3.1)\begin{equation} \frac{\partial f(\boldsymbol{v};t)}{\partial t}=J[\boldsymbol{v}\,|\,f,f]+J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

Moreover, the velocity distribution $f(\boldsymbol {v};t)$ of the granular gas is isotropic in $\boldsymbol {v}$ so that the production of momentum $\boldsymbol {\mathcal {F}}[f]=\boldsymbol {0}$, according to (2.12). Thus, the only non-trivial balance equation is that of the temperature (2.10), namely $\partial _t \ln T=-(\zeta +\zeta _g)$. As mentioned in § 2, since collisions among granular particles are inelastic, the cooling rate $\zeta >0$. A collision between a particle of the granular gas and a particle of the molecular gas is elastic, and so the total kinetic energy of two colliding particles in such a collision is conserved. On the other hand, since in the steady state the background gas acts as a thermostat, the mean kinetic energy of granular particles is smaller than that of the molecular gas, and so $T< T_g$. This necessarily implies that $\zeta _g<0$. Therefore, in the steady state, the terms $\zeta$ and $|\zeta _g|$ exactly compensate each other and one gets the steady-state condition $\zeta +\zeta _g=0$. This condition allows one to get the steady granular temperature $T$. However, according to the definitions (2.13a,b), the determination of $\zeta$ and $\zeta _g$ requires knowledge of the velocity distribution $f(\boldsymbol {v})$. For inelastic collisions ($\alpha \neq 1$), to date the solution of the Boltzmann equation (3.1) has not been found. On the other hand, a good estimate of $\zeta$ and $\zeta _g$ can be obtained when the first Sonine approximation to $f$ is considered (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004). In this approximation, $f(\boldsymbol {v})$ is given by

(3.2)\begin{equation} f(\boldsymbol{v})\simeq f_{{MB}}(\boldsymbol{v})\left\{1+\frac{a_2}{2} \left[\left(\frac{m v^{2}}{2T}\right)^{2}-(d+2)\frac{mv^{2}}{2T}+\frac{d(d+2)}{4}\right]\right\}, \end{equation}

where

(3.3)\begin{equation} f_{{MB}}(\boldsymbol{v})=n \biggl(\frac{m}{2{\rm \pi} T}\biggr)^{d/2} \exp \left(-\frac{m v^{2}}{2T}\right) \end{equation}

is the Maxwell–Boltzmann distribution and

(3.4)\begin{equation} a_2=\frac{1}{d(d+2)}\frac{m^{2}}{n T^{2}}\int {\rm d}\boldsymbol{v}v^{4} f(\boldsymbol{v})-1 \end{equation}

is the kurtosis or fourth cumulant. This quantity measures the departure of the distribution $f(\boldsymbol {v})$ from its Maxwellian form $f_{{MB}}(\boldsymbol {v})$. From experience with the dry granular case (van Noije & Ernst Reference van Noije and Ernst1998; Garzó & Dufty Reference Garzó and Dufty1999; Montanero & Santos Reference Montanero and Santos2000; Santos & Montanero Reference Santos and Montanero2009), the magnitude of the cumulant $a_2$ is expected to be very small, and so the Sonine approximation (3.2) to the distribution $f$ turns out to be reliable. In the case that $|a_2|$ does not remain small for high inelasticity, one should include cumulants of higher order in the Sonine polynomial expansion of $f$. However, the possible lack of convergence of the Sonine polynomial expansion for very small values of the coefficient of restitution (Brilliantov & Pöschel Reference Brilliantov and Pöschel2006a,Reference Brilliantov and Pöschelb) puts in doubt the reliability of the Sonine expansion in the high-inelasticity region. Here, we restrict ourselves to values of $\alpha$ where $|a_2|$ remains relatively small.

The expressions of $\zeta$ and $\zeta _g$ can now be obtained by replacing in (2.13a,b) $f$ by its Sonine approximation (3.2). Retaining only linear terms in $a_2$, the forms of the dimensionless production rates $\zeta ^{*}=(\ell \zeta /v_{th})$ and $\zeta _g^{*}=(\ell \zeta /v_{th})$ can be written as (van Noije & Ernst Reference van Noije and Ernst1998; Brilliantov & Pöschel Reference Brilliantov and Pöschel2006a)

(3.5a,b)\begin{equation} \zeta^{*}=\tilde{\zeta}^{(0)}+\tilde{\zeta}^{(1)} a_2, \quad \zeta_g^{*}=\tilde{\zeta}_g^{(0)}+\tilde{\zeta}_g^{(1)} a_2, \end{equation}

where

(3.6a,b)$$\begin{gather} \tilde{\zeta}^{(0)}=\frac{\sqrt{2}{\rm \pi}^{(d-1)/2}}{d\varGamma\left(\dfrac{d}{2}\right)}(1-\alpha^{2}),\quad \tilde{\zeta}^{(1)}=\frac{3}{16} \tilde{\zeta}^{(0)}, \end{gather}$$

(3.7a,b)$$\begin{gather}\tilde{\zeta}_g^{(0)}=2x(1-x^{2}) \left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*},\quad \tilde{\zeta}_g^{(1)}=\frac{\mu_g}{8}x^{{-}3}[x^{2}\left(4-3\mu_g\right)-\mu_g] \left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*}. \end{gather}$$

Here, $\ell =1/(n\sigma ^{d-1})$ is proportional to the mean free path of hard spheres, $v_{th}=\sqrt {2T/m}$ is the thermal velocity and we have introduced the auxiliary parameters

(3.8)\begin{equation} x=\left(\mu_g+\mu \frac{T_g}{T}\right)^{1/2} \end{equation}

and

(3.9a,b)\begin{equation} \gamma^{*}=\varepsilon \left(\frac{T_g}{T}\right)^{1/2},\quad \varepsilon=\frac{\ell \gamma}{\sqrt{2T_g/m}}= \frac{\sqrt{2}{\rm \pi}^{d/2}}{2^{d} d \varGamma\left(\dfrac{d}{2}\right)}\frac{1}{\phi \sqrt{T_g^{*}}}. \end{equation}

Here, $\phi =[{\rm \pi} ^{d/2}/2^{d-1}d\varGamma (d/2)]n\sigma ^{d}$ is the solid volume fraction and $T_g^{*}=T_g/(m\sigma ^{2} \gamma ^{2})$ is the (reduced) bath temperature. The dimensionless coefficient $\gamma ^{*}$ characterises the rate at which the collisions between grains and molecular particles occur. Equations (3.6a,b) and (3.7a,b) agree with those obtained by Santos (Reference Santos2003) for $d=3$.

To close the problem, we have to determine the kurtosis $a_2$. In this case, one has to compute the collisional moments

(3.10a,b)\begin{equation} \varLambda\equiv \int {\rm d}\boldsymbol{v}v^{4} J[\boldsymbol{v}\,|\,f,f],\quad \varLambda_g\equiv \int {\rm d}\boldsymbol{v} v^{4} J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

In the steady state, one has the additional condition $\varLambda +\varLambda _g=0$. The moments $\varLambda$ and $\varLambda _g$ have been obtained in previous works (van Noije & Ernst Reference van Noije and Ernst1998; Brilliantov & Pöschel Reference Brilliantov and Pöschel2006a; Garzó, Vega Reyes & Montanero Reference Garzó, Vega Reyes and Montanero2009; Garzó Reference Garzó2019) by replacing $f$ by its first Sonine form (3.2) and neglecting nonlinear terms in $a_2$. In terms of $\gamma ^{*}$, the expressions of $\{\varLambda ^{*},\varLambda _g^{*}\}=(\ell /(n v_{th}^{5}))\{\varLambda,\varLambda _g\}$ are given by

(3.11a,b)\begin{equation} \varLambda^{*}=\varLambda^{(0)}+\varLambda^{(1)} a_2,\quad \varLambda_g^{*}=\varLambda_g^{(0)}+\varLambda_g^{(1)} a_2, \end{equation}

where

(3.12)$$\begin{gather} \varLambda^{(0)}={-}\frac{{\rm \pi}^{(d-1)/2}}{\sqrt{2}\varGamma\left(\dfrac{d}{2}\right)} \left(d+\frac{3}{2}+\alpha^{2}\right)(1-\alpha^{2}), \end{gather}$$

(3.13)$$\begin{gather}\varLambda^{(1)}={-}\frac{{\rm \pi}^{(d-1)/2}}{\sqrt{2}\varGamma\left(\dfrac{d}{2}\right)} \left[\frac{3}{32}\left(10d+39+10\alpha^{2}\right)+\frac{d-1}{1-\alpha}\right](1-\alpha^{2}), \end{gather}$$

(3.14)$$\begin{gather}\varLambda_g^{(0)}={{d}x}^{{-}1}\left(x^{2}-1\right) [8 \mu_g x^{4}+x^{2}\left(d+2-8\mu_g\right)+\mu_g]\left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*}, \end{gather}$$

(3.15)\begin{gather} \hspace{-2pc}\varLambda_g^{(1)} = \frac{d}{8} x^{{-}5}\{4x^{6}[30\mu_g^{3}-48\mu_g^{2}+3(d+8)\mu_g-2(d+2)]+ \mu_g x^{4} [{-}48\mu_g^{2} \nonumber\\ \hspace{4pc} +\,3(d+26)\mu_g-8(d+5)]+\mu_g^{2} x^{2}(d+14-9\mu_g)-3\mu_g^{3}\} \left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*}. \end{gather}

For hard spheres ($d=3$), equations (3.12)–(3.15) are consistent with those previously obtained by Santos (Reference Santos2003).

Inserting (3.5a,b)–(3.7a,b) and (3.11a,b)–(3.15) into the steady-state conditions ($\zeta +\zeta _g=0$, $\varLambda +\varLambda _g=0$), one gets a set of coupled equations:

(3.16)$$\begin{gather} \tilde{\zeta}^{(0)}+\tilde{\zeta}_g^{(0)}+(\tilde{\zeta}^{(1)}+\tilde{\zeta}_g^{(1)})a_2=0, \end{gather}$$

(3.17)$$\begin{gather}\varLambda^{(0)}+\varLambda_g^{(0)}+(\varLambda^{(1)}+\varLambda_g^{(1)})a_2=0. \end{gather}$$

Eliminating $a_2$ in (3.16) and (3.17), one achieves the following closed equation for the temperature ratio $T/T_g$:

(3.18)\begin{equation} (\tilde{\zeta}^{(1)}+\tilde{\zeta}_g^{(1)})(\varLambda^{(0)}+\varLambda_g^{(0)})= (\tilde{\zeta}^{(0)}+\tilde{\zeta}_g^{(0)})(\varLambda^{(1)}+\varLambda_g^{(1)}). \end{equation}

For given values of $\alpha$, $\phi$ and $T_g^{*}$, the solution of (3.18) gives $T/T_g$. Once the temperature ratio is determined, the cumulant $a_2$ is simply given by

(3.19)\begin{equation} a_2={-}\frac{\tilde{\zeta}^{(0)}+\tilde{\zeta}_g^{(0)}}{\tilde{\zeta}^{(1)}+\tilde{\zeta}_g^{(1)}}={-}\frac{\varLambda^{(0)}+\varLambda_g^{(0)}}{\varLambda^{(1)}+\varLambda_g^{(1)}}. \end{equation}

The set of dimensionless control parameters of the problem considered here ($m/m_g$, $\alpha$, $\phi$, $T_g^{*}$) has been essentially chosen to perform a close and clean comparison with the previous results obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) by means of the suspension model (2.16). On the other hand, another possible set of parameters are $m/m_g$, $\alpha$ and $\omega \equiv (n\sigma ^{d-1})/(n_g\overline {\sigma }^{d-1})$. The parameter $\omega$ represents the mean free path associated with the grain–gas collisions relative to that associated with grain–grain collisions. In fact, $\omega$ was considered independently by Biben et al. (Reference Biben, Martin and Piasecki2002) and Santos (Reference Santos2003) in their study of the HSS. The relationship between $\omega$ and $T_g^{*}$ is

(3.20)\begin{equation} \omega=\frac{2^{d+{3}/{2}}}{\sqrt{\rm \pi}}\left(\frac{m_g}{m}\right)^{1/2}\phi \sqrt{T_g^{*}}. \end{equation}

Note that $\omega$ encompasses the dependence on the volume fraction $\phi$ and the (reduced) bath temperature $T_g^{*}$ through the combination $\phi \sqrt {T_g^{*}}$. For this reason, according to (3.9a,b), the number of independent parameters in this set is $\alpha$, $m/m_g$ and $\omega$.

3.1. Brownian limit

Before illustrating the dependence of $T/T_g$ and $a_2$ on $\alpha$ for given values of $m/m_g$, $\phi$ and $T_g^{*}$, it is interesting to consider the Brownian limit $m/m_g\to \infty$. In this limiting case, $\mu _g\to 0$, $\mu \to 1$, $x\to \sqrt {T_g/T}$ and so

(3.21a,b)$$\begin{gather} \tilde{\zeta}_g^{(0)}\to 2\left(1-\frac{T_g}{T}\right)\gamma^{*}, \quad \tilde{\zeta}_g^{(1)}\to 0, \end{gather}$$

(3.22a,b)$$\begin{gather}\varLambda_g^{(0)}\to d(d+2)\left(\frac{T_g}{T}-1\right)\gamma^{*}, \quad \varLambda_g^{(1)}\to -d(d+2)\gamma^{*}. \end{gather}$$

Taking into account these results, the set of (3.16) and (3.17) can be written in the Brownian limit as

(3.23a,b)\begin{equation} 2\gamma^{*}\left(\frac{T_g}{T}-1\right)=\zeta^{*},\quad d(d+2)\left(\gamma^{*} a_2-\frac{1}{2}\zeta^{*}\right)=\varLambda^{*}. \end{equation}

These equations are the same as those derived by Gómez González & Garzó (Reference Gómez González and Garzó2019) (see (29) and (34) of that paper) by using the suspension model (2.16). This shows the consistency of the present results in the HSS with those obtained in the Brownian limit.

3.2. Direct simulation Monte Carlo simulations

The previous analytical results have been obtained by using the first Sonine approximation (3.2) to $f$. Thus, it is worth solving the Boltzmann kinetic equation by means of an alternative method to test the reliability of the theoretical predictions for $T/T_g$ (3.18) and $a_2$ (3.19). The direct simulation Monte Carlo (DSMC) method developed by Bird (Reference Bird1994) is considered here to numerically solve the Boltzmann equation in the homogeneous state. Some technical details of the application of the DSMC method to the system studied in this paper are provided as supplementary material available at https://doi.org/10.1017/jfm.2022.410.

The dependence of the temperature ratio $\chi \equiv T/T_g$ on the coefficient of restitution $\alpha$ is plotted in figure 1 for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and several values of the mass ratio $m/m_g$. The value $T_g^{*}=1000$ has been chosen to guarantee that the grain–grain collisions play a relevant role in the dynamics of granular gas. Namely, that the value of the friction coefficient $\gamma$ is comparable to the value of the grain–grain collision frequency ($\nu =v_{th}/\ell$), so that the dimensionless coefficient $\gamma ^{*}\equiv \gamma /\nu$ is of the order of unity. This is fulfilled for the system studied here ($d=3$, $\phi =0.001$ and $T_g^{*}=1000$) since $\gamma ^{*}\simeq 10\sqrt {T_g/T}$, $T_g/T$ being of the order of unity (see figure 1). Thus, although the results displayed in this paper significantly differ from those found in the dry case (no gas phase), the effects of inelastic collisions still have importance for the dynamics of grains. The value $T_g^{*}=1000$ will therefore be maintained throughout this work.

Figure 1. Temperature ratio $\chi \equiv T/T_g$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=50, 10, 5$ and 1). The solid lines are the theoretical results obtained by numerically solving (3.18) and the symbols are the Monte Carlo simulation results. The dotted line is the result obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the Langevin-like suspension model (2.16) while black circles refer to DSMC simulations implemented using the time-driven approach (see the supplementary material).

Theoretical results are compared against DSMC simulations in figure 1, which ensures the reliability of the results derived in this section for two different reasons: (i) a good agreement between theory and simulation is found and (ii) the convergence towards the Brownian limit can be clearly observed. Surprisingly, this convergence is fully reached for relatively small values of the mass ratio ($m/m_g\approx 50$). We also find that the departure of $\chi$ from unity increases as the masses of the granular and gas particles are comparable. However, this unexpected result is only due to the way of scaling the variables. This is illustrated in figure 2 where we take $\omega$ instead of $T_g^{*}$ as input as in figure 2 of Santos (Reference Santos2003). In contrast to figure 1, as expected (Barrat & Trizac Reference Barrat and Trizac2002; Dahl et al. Reference Dahl, Hrenya, Garzó and Dufty2002), it is quite apparent that the lack of energy equipartition is more noticeable as $m\gg m_g$. In addition, we also see that the impact of the mass ratio on the temperature ratio is apparently more significant when one fixes $\omega$ instead of $T_g^{*}$. The fact that the difference $1-\chi$ increases with decreasing mass ratio when $T_g^{*}$ is fixed (see figure 1) can be easily understood. According to (2.6a,b), the transmission of energy per individual collision from a molecular particle to a grain is greater when their masses are similar. Nonetheless, the constraint imposed by the way of scaling $\gamma$ leads to a dependence of $N/N_g$ on the mass ratio $m/m_g$ for fixed $\sigma _g$. Thus, $N_g/N\propto m/m_g$, and so the number density of the molecular gas increases with increasing mass ratio. In this way, the mean force exerted by the molecular particles on the grains is greater, and therefore the thermalisation caused by the presence of the interstitial fluid is much more effective. The steady temperature ratio $\chi$ is reached when the energy lost by collisions is compensated for by the energy provided by the bath. Hence, the non-equipartition of energy turns out then to be remarkable to small values of $m/m_g$ and $\alpha$. Finally, figure 3 shows the $\alpha$ dependence of the cumulant $a_2$ for the same parameters as in figure 1. As expected, we find that the magnitude of $a_2$ is in general small for not quite large inelasticity (for instance, $\alpha \gtrsim 0.5$); this result supports the assumption of a low-order truncation (first Sonine approximation) in the polynomial expansion of the distribution function. Figure 3 highlights the excellent agreement between theory and simulations, except for $m/m_g=1$ where small differences are present for very strong inelasticities. However, these discrepancies are of the same order as those found for dry granular gases (Montanero & Garzó Reference Montanero and Garzó2002).

Figure 2. Temperature ratio $\chi \equiv T/T_g$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $\omega =0.1$ and four different values of the mass ratio $m/m_g$: $m/m_g=1$ (solid line), $m/m_g=10$ (dashed line), $m/m_g=100$ (dotted line) and $m/m_g=1000$ (dash-dotted line). The (reduced) bath temperature $T_g^{*}=61.36 (m/m_g)$.

Figure 3. Plot of the fourth cumulant $a_2$ as a function of the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). The solid lines are the theoretical results obtained from (3.18) and the symbols are the Monte Carlo simulation results. The dotted line is the result obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the Langevin-like suspension model (2.16) while black circles refer to DSMC simulations implemented using the time-driven approach (see the supplementary material).

4. Chapman–Enskog expansion. First-order approximation

We perturb now the homogeneous state by small spatial gradients. These perturbations will give non-zero contributions to the pressure tensor and the heat flux vector. The determination of these fluxes will allow us to identify the Navier–Stokes–Fourier transport coefficients of the granular gas. For times longer than the mean free time, we assume that the system evolves towards a hydrodynamic regime where the distribution function $f(\boldsymbol {r}, \boldsymbol {v};t)$ adopts the form of a normal or hydrodynamic solution. This means that all space and time dependence of $f$ only occurs through the hydrodynamic fields $n$, $\boldsymbol {U}$ and $T$:

(4.1)\begin{equation} f(\boldsymbol{r},\boldsymbol{v},t)=f\left[\boldsymbol{v}\,|\,n (t), \boldsymbol{U}(t), T(t) \right]. \end{equation}

The notation on the right-hand side indicates a functional dependence on the density, flow velocity and temperature. For low Knudsen numbers (i.e. small spatial variations), the functional dependence (4.1) can be made local in space by means of an expansion in powers of the gradients $\boldsymbol {\nabla } n$, $\boldsymbol {\nabla } \boldsymbol {U}$ and $\boldsymbol {\nabla } T$ (Chapman & Cowling Reference Chapman and Cowling1970). In this case, $f$ can be expressed in the form $f=f^{(0)}+f^{(1)}+f^{(2)}+\cdots$, where the approximation $f^{(k)}$ is of order $k$ in spatial gradients. Here, since we are interested in the Navier–Stokes hydrodynamic equations, only terms up to first order in gradients are considered in the constitutive equations for the momentum and heat fluxes.

As has been noted in previous works on granular suspensions (Garzó et al. Reference Garzó, Chamorro and Vega Reyes2013; Gómez González & Garzó Reference Gómez González and Garzó2019), although one is interested in computing transport in steady conditions, the presence of the background molecular gas may induce a local energy unbalance between the energy lost due to inelastic collisions and the energy transfer via elastic collisions. Thus, we have to consider first the time-dependent distribution $f^{(0)}(\boldsymbol {r},\boldsymbol {v}; t)$ in order to arrive at the linear integral equations obeying the Navier–Stokes–Fourier transport coefficients. Then, to get explicit forms for the transport coefficients, the above integral equations are (approximately) solved under steady-state conditions.

To collect the different approximations in (2.2), one has to characterise the magnitude of the velocity difference ${\rm \Delta} \boldsymbol {U}$ relative to the gradients as well. In the absence of spatial gradients, the production of momentum term $\mathcal {F}[f]\propto {\rm \Delta} \boldsymbol {U}$ (see (6.5)) and hence, according to the momentum balance equation (2.9), the mean flow velocity $\boldsymbol {U}$ of the granular gas relaxes towards that of the molecular gas $\boldsymbol {U}_g$ after a transient period. Thus, the term ${\rm \Delta} \boldsymbol {U}$ must be considered to be at least of first order in the spatial gradients. In this case, the Maxwellian distribution $f_g(\boldsymbol {v})$ must also be expanded as

(4.2)\begin{equation} f_g(\boldsymbol{v})=f_g^{(0)}(\boldsymbol{V})+f_g^{(1)}(\boldsymbol{V})+\cdots, \end{equation}

where

(4.3)\begin{equation} f_g^{(0)}(\boldsymbol{V})=n_g \left(\frac{m_g}{2{\rm \pi} T_g}\right)^{d/2} \exp \bigg(-\frac{m_g V^{2}}{2T_g}\bigg) \end{equation}

and

(4.4)\begin{equation} f_g^{(1)}(\boldsymbol{V})={-}\frac{m_g}{T_g}\boldsymbol{V}\boldsymbol{\cdot} {\rm \Delta} \boldsymbol{U} f_g^{(0)}(\boldsymbol{V}). \end{equation}

According to the expansion (4.1), the pressure tensor $P_{ij}$, the heat flux $\boldsymbol {q}$ and the partial production rates $\zeta$ and $\zeta _g$ must also be expressed according to the perturbation scheme in the forms

(4.5a–d)\begin{equation} \left.\begin{gathered} P_{ij}=P_{ij}^{(0)}+P_{ij}^{(1)}+\cdots,\quad \boldsymbol{q}=\boldsymbol{q}^{(0)}+\boldsymbol{q}^{(1)}+\cdots,\\ \zeta=\zeta^{(0)}+\zeta^{(1)}+\cdots,\quad \zeta_g=\zeta_g^{(0)}+\zeta_g^{(1)}+\cdots . \end{gathered}\right\} \end{equation}

In addition, the time derivative $\partial _t$ is also given as $\partial _t=\partial _t^{(0)}+\partial _t^{(1)}+\cdots$. The action of the operators $\partial _t^{(k)}$ on the hydrodynamic fields can be identified when the expansions (4.5a–d) of the fluxes and the production rates are considered in the macroscopic balance equations (2.8)–(2.10). This is the conventional Chapman–Enskog method (Chapman & Cowling Reference Chapman and Cowling1970; Garzó Reference Garzó2019) for solving the Boltzmann kinetic equation.

As usual in the Chapman–Enskog method (Chapman & Cowling Reference Chapman and Cowling1970), the zeroth-order distribution function $f^{(0)}$ defines the hydrodynamic fields $n$, $\boldsymbol {U}$ and $T$:

(4.6)\begin{equation} \left\{n, n \boldsymbol{U},d n T\right\}=\int {\rm d}\boldsymbol{v}\{1,\boldsymbol{v}, m V^{2}\}f^{(0)}(\boldsymbol{V}). \end{equation}

The requirements (4.6) must be fulfilled at any order in the expansion, and so the distributions $f^{(k)}$ ($k\geq 1$) must thus obey the orthogonality conditions

(4.7)\begin{equation} \int {\rm d}\boldsymbol{v}\{1,\boldsymbol{v}, m V^{2}\}f^{(k)}(\boldsymbol{V})=\{0,\boldsymbol{0},0\}. \end{equation}

These are the usual solubility conditions of the Chapman–Enskog scheme.

The mathematical steps involved in the determination of the zeroth- and first-order distributions are quite similar to those made in previous works (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Garzó & Dufty Reference Garzó and Dufty1999; Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012, Reference Garzó, Chamorro and Vega Reyes2013; Gómez González & Garzó Reference Gómez González and Garzó2019). Some of the technical details involved in this derivation are given in the supplementary material.

4.1. Navier–Stokes transport coefficients

To first order in spatial gradients and based on symmetry considerations, the pressure tensor $P_{ij}^{(1)}$ and the heat flux $\boldsymbol {q}^{(1)}$ are given, respectively, by

(4.8)$$\begin{gather} P_{ij}^{(1)}={-}\eta \left(\frac{\partial U_i}{\partial r_j}+ \frac{\partial U_j}{\partial r_i}-\frac{2}{d}\delta_{ij}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{U}\right), \end{gather}$$

(4.9)$$\begin{gather}\boldsymbol{q}^{(1)}={-}\kappa \boldsymbol{\nabla} T-\bar{\mu} \boldsymbol{\nabla} n-\kappa_U {\rm \Delta} \boldsymbol{U}. \end{gather}$$

Here, $\eta$ is the shear viscosity, $\kappa$ is the thermal conductivity, $\bar {\mu }$ is the diffusive heat conductivity and $\kappa _U$ is the velocity conductivity. While $\eta$, $\kappa$ and $\mu$ are the coefficients of proportionality between fluxes and hydrodynamic gradients, the coefficient $\kappa _U$ connects the heat flux with the velocity difference ${\rm \Delta} \boldsymbol {U}$ (‘convection’ current). Although this contribution is not present in dry granular gases (Garzó Reference Garzó2019), it also appears in the case of driven granular mixtures (Khalil & Garzó Reference Khalil and Garzó2013, Reference Khalil and Garzó2018). The coefficient $\kappa _U$ can be seen as a measure of the contribution to the heat flow due to ‘diffusion’ (in the sense that we have a ‘binary mixture’ of granular and gas particles where both species have different mean velocities). In this context, $\kappa _U$ can be regarded as an effect inverse to thermal diffusion (diffusion thermo-effect) (Chapman & Cowling Reference Chapman and Cowling1970). It is important to recall that ${\rm \Delta} \boldsymbol {U}$ vanishes for HSSs, and hence one expects that ${\rm \Delta} \boldsymbol {U}$ can be expressed in terms of $\boldsymbol {\nabla } T$ and $\boldsymbol {\nabla } n$ in particular inhomogeneous situations.

The Navier–Stokes–Fourier transport coefficients $\eta$, $\kappa$ and $\mu$ are defined, respectively, as

(4.10)$$\begin{gather} \eta={-}\frac{1}{(d-1)(d+2)}\int {\rm d}\boldsymbol{v} R_{ij}(\boldsymbol{V})\mathcal{C}_{ij}(\boldsymbol{V}), \end{gather}$$

(4.11a,b)$$\begin{gather}\kappa={-}\frac{1}{dT}\int {\rm d}\boldsymbol{v}\boldsymbol{S}(\boldsymbol{V})\boldsymbol{\cdot} \boldsymbol{\mathcal{A}}(\boldsymbol{V}),\quad \bar{\mu}={-}\frac{1}{dn}\int {\rm d}\boldsymbol{v}\boldsymbol{S}(\boldsymbol{V})\boldsymbol{\cdot} \boldsymbol{\mathcal{B}}(\boldsymbol{V}), \end{gather}$$

while $\kappa _U$ is

(4.12)\begin{equation} \kappa_U={-}\frac{1}{d}\int {\rm d}\boldsymbol{v}\boldsymbol{S}(\boldsymbol{V})\boldsymbol{\cdot} \boldsymbol{\mathcal{E}}(\boldsymbol{V}). \end{equation}

In (4.10)–(4.12), we have introduced the quantities

(4.13a,b)\begin{equation} R_{ij}(\boldsymbol{V})=m\left(V_i V_j-\frac{1}{d}V^{2} \delta_{ij}\right),\quad \boldsymbol{S}(\boldsymbol{V})=\left(\frac{m}{2}V^{2}-\frac{d+2}{2}T\right)\boldsymbol{V}. \end{equation}

5. Sonine polynomial approximation to the transport coefficients in steady-state conditions

So far, all the results displayed in § 4 for the transport coefficients $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$ are exact. More specifically, their expressions are given by (4.10)–(4.12), where the unknowns $\boldsymbol {\mathcal {A}}$, $\boldsymbol {\mathcal {B}}$, $\mathcal {C}_{ij}$, $\mathcal {D}$ and $\boldsymbol {\mathcal {E}}$ are the solutions of a set of coupled linear integral equations displayed in the supplementary material. However, it is easy to see that the solution for general unsteady conditions requires one to solve numerically a set of coupled differential equations for $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$. Thus, in a desire of achieving analytical expressions of the transport coefficients, we consider steady-state conditions. In this case, the constraint $\zeta ^{(0)}+\zeta _g^{(0)}=0$ applies locally, and so the transport coefficients can be explicitly obtained. The procedure for deriving the expressions of the transport coefficients is described in the supplementary material and only their final forms are provided here.

5.1. Shear viscosity

The shear viscosity coefficient $\eta$ is given by

(5.1)\begin{equation} \eta=\frac{\eta_0}{\nu_\eta^{*}+K'\tilde{\nu}_\eta\gamma^{*}}, \end{equation}

where $\gamma ^{*}$ is defined in (3.9a,b), $\eta _0=[(d+2)\varGamma (d/2)/(8{\rm \pi} ^{(d-1)/2})]\sigma ^{1-d}\sqrt {m T}$ is the low-density value of the shear viscosity of an ordinary gas of hard spheres ($\alpha =1$) and $K'= \sqrt {2}(d+2)\varGamma (d/2)/(8{\rm \pi} ^{(d-1)/2})$. Moreover, we have introduced the (reduced) collision frequencies

(5.2)\begin{equation} \nu_\eta^{*}=\frac{3}{4d}\left(1-\alpha+\frac{2}{3}d\right)(1+\alpha), \end{equation}

(5.3)\begin{align} \tilde{\nu}_\eta &= \frac{1}{(d-1)(d+2)}\left(\frac{m}{m_g}\right)^{3} \mu_g \left(\frac{T_g}{T}\right)^{2}\theta^{{-}1/2} \nonumber\\ &\quad \times\{2(d+3)(d-1)\left(\mu-\mu_g \theta\right)\theta^{{-}2}(1+\theta)^{{-}1/2} \nonumber\\ &\quad +2d(d-1)\mu_g \theta^{{-}2}(1+\theta)^{1/2}+2(d+2)(d-1)\theta^{{-}1}(1+\theta)^{{-}1/2}\}, \end{align}

where $\theta =m T_g/(m_g T)$ is the ratio of the mean square velocities of granular and molecular gas particles. It is important to recall that all the quantities appearing in (5.1) are evaluated at the steady-state conditions.

5.2. Thermal conductivity, diffusive heat conductivity and velocity conductivity

We consider here the transport coefficients associated with the heat flux. The thermal conductivity coefficient $\kappa$ is

(5.4)\begin{equation} \kappa=\frac{d-1}{d}\frac{\kappa_0}{\nu_\kappa^{*}+K'\left(\tilde{\nu}_\kappa+\beta\right)\gamma^{*}}, \end{equation}

where $\kappa _0=[d(d+2)/2(d-1)](\eta _0/m)$ is the low-density value of the thermal conductivity for an ordinary gas of hard spheres and

(5.5)\begin{equation} \beta=\left(x^{{-}1}-3x\right)\mu^{3/2}\left(\frac{T_g}{T}\right)^{1/2}. \end{equation}

In (5.4), we have introduced the (reduced) collision frequencies

(5.6)$$\begin{gather} \nu_\kappa^{*}=\frac{1+\alpha}{d}\left[\frac{d-1}{2}+\frac{3}{16}(d+8)(1-\alpha)\right](1+\alpha), \end{gather}$$

(5.7)$$\begin{gather}\tilde{\nu}_\kappa=\frac{1}{2(d+2)}\mu \frac{\theta}{1+\theta} \left[G-(d+2)\frac{1+\theta}{\theta}F\right], \end{gather}$$

where

(5.8)\begin{align} F &= (d+2)(2 \delta +1)+4(d-1)\mu_g \delta \theta^{{-}1}(1+\theta)+3(d+3)\delta^{2}\theta^{{-}1} \nonumber\\ &\quad +(d+3)\mu_g^{2} \theta^{{-}1}(1+\theta)^{2} -(d+2)\theta^{{-}1}(1+\theta), \end{align}

(5.9)\begin{align} G &= (d+3)\mu_g^{2} \theta^{{-}2}(1+\theta)^{2}\left[d+5+(d+2)\theta\right] \nonumber\\ &\quad -\mu_g(1+\theta)\{4(1-d)\delta \theta^{{-}2}\left[d+5+(d+2)\theta\right]-8(d-1)\theta^{{-}1}\} \nonumber\\ &\quad +3(d+3)\delta^{2}\theta^{{-}2}\left[d+5+(d+2)\theta\right]+2\delta\theta^{{-}1} [24+11d+d^{2}+(d+2)^{2}\theta] \nonumber\\ &\quad +(d+2)\theta^{{-}1}\left[d+3+(d+8)\theta\right]-(d+2)\theta^{{-}2}(1+\theta)\left[d+3+(d+2)\theta\right] \end{align}

and $\delta \equiv \mu -\mu _g \theta$.

The diffusive heat conductivity $\bar {\mu }$ can be written as

(5.10)\begin{equation} \bar{\mu}=\frac{K' T}{n}\frac{\kappa \zeta^{*}}{\nu_\kappa^{*}+K'\tilde{\nu}_\kappa \gamma^{*}}. \end{equation}

Finally, the velocity conductivity $\kappa _U$ is given by

(5.11)\begin{equation} \kappa_U={-}\frac{n T}{2}\frac{K' \mu (1+\theta)^{{-}1/2}\theta^{{-}1/2}H}{\nu_\kappa^{*}+K'\tilde{\nu}_\kappa \gamma^{*}}\gamma^{*}, \end{equation}

where

(5.12)\begin{equation} H=(d+2)(1+2\delta)+4(1-d)\mu_g (1+\theta) \delta-3(d+3)\delta^{2}-(d+3)\mu_g^{2} (1+\theta)^{2}. \end{equation}

5.3. Brownian limit

Equations (5.1), (5.4), (5.10) and (5.11) provide the expressions of the transport coefficients $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$, respectively, for arbitrary values of the mass ratio $m/m_g$. As before in the homogeneous state, it is quite interesting to consider the limiting case $m/m_g \to \infty$ (Brownian limit). In this limit case, $\mu _g\to 0$, $\mu \to 1$, $\chi \equiv \textrm {finite}$, and so $\theta \to 0$, $x\to \chi ^{-1/2}$, $\delta \to 1-\chi ^{-1}$ and $\beta \to 1-3\chi ^{-1}$. This yields the results $\tilde {\nu }_\eta \to 2$ and $\tilde {\nu }_\kappa \to 3$, so that in the Brownian limit (5.1), (5.4), (5.10) and (5.11) reduce to (note that there is a typographical error in expression (78) of Gómez González & Garzó (Reference Gómez González and Garzó2019) for the coefficient $\bar {\mu }$ since the denominator should be $\nu _\kappa ^{*}+3K'\gamma ^{*}$)

(5.13a,b)$$\begin{gather} \eta\to \frac{\eta_0}{\nu_\eta^{*}+2K'\gamma^{*}},\quad \kappa\to \frac{d-1}{d}\frac{\kappa_0}{\nu_\kappa^{*}+K'\left(\gamma^{*}- \frac{3}{2}\zeta^{*}\right)}, \end{gather}$$

(5.14a,b)$$\begin{gather}\bar{\mu}\to \frac{\kappa T}{n}\frac{K'\zeta^{*}}{\nu_\kappa^{*}+3K' \gamma^{*}},\quad \kappa_U \to 0. \end{gather}$$

Equations (5.13a,b) and (5.14a,b) agree with the results obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) by using the suspension model (2.16). This confirms the self-consistency of the results obtained in this paper for general values of the mass ratio.

5.4. Some illustrative systems

In the steady state, the expressions of the Navier–Stokes–Fourier transport coefficients $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$ are provided by (5.1), (5.4), (5.10) and (5.11), respectively. As in previous works on transport in granular gases (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Garzó & Dufty Reference Garzó and Dufty1999; Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012; Gómez González & Garzó Reference Gómez González and Garzó2019), to highlight the $\alpha$ dependence of the transport coefficients, they are scaled with respect to their values for elastic collisions. This scaling cannot be made in the case of the diffusive heat conductivity $\bar {\mu }$ since this coefficient vanishes for $\alpha =1$. In this case, we consider the scaled coefficient $n\bar {\mu }/T\kappa (1)$, where $\kappa (1)$ refers to the value of the thermal conductivity (5.4) for elastic collisions. All these scaled coefficients exhibit a complex dependence on the coefficient of restitution $\alpha$, the mass ratio $m/m_g$, the volume fraction $\phi$ (through the parameter $\varepsilon$ defined by (3.9a,b)) and the reduced temperature $T_g^{*}$ of the molecular gas.

Figures 4–7 show $\eta (\alpha )/\eta (1)$, $\kappa (\alpha )/\kappa (1)$, $n\bar {\mu }/T\kappa (1)$ and $\kappa _U(\alpha )/\kappa _U(1)$, respectively, as functions of the coefficient of restitution $\alpha$. Here, $\eta(1)$, $\kappa(1)$ and $\kappa_U(1)$ correspond to the values of $\eta$, $\kappa$ and $\kappa_U$ for elastic collisions. Moreover, in those plots we consider a three-dimensional system ($d=3$) with $\phi =0.001$ (very dilute granular gas), $T_g^{*}=1000$ and four different values of the mass ratio: $m/m_g=1, 5, 10$ and 50. We have also plotted the results obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the suspension model (2.16). The results obtained from this model are expected to apply when $m \gg m_g$ (Brownian limit).

Figure 4. Plot of the (scaled) shear viscosity coefficient $\eta (\alpha )/\eta (1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=50, 10, 5$ and 1). The solid lines are the results derived in this paper while the dotted line is the result obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the suspension model (2.16). Here, $\eta (1)$ refers to the shear viscosity coefficient when collisions between grains are elastic ($\alpha =1$).

Figure 5. Plot of the (scaled) thermal conductivity coefficient $\kappa (\alpha )/\kappa (1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). The solid lines are the results derived in this paper while the dotted line is the result obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the suspension model (2.16). Here, $\kappa (1)$ refers to the thermal conductivity coefficient when collisions between grains are elastic ($\alpha =1$).

Figure 6. Plot of the (scaled) diffusive heat conductivity coefficient $n\bar {\mu }(\alpha )/T\kappa (1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). The solid lines are the results derived in this paper while the dotted line is the result obtained by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the suspension model (2.16). Here, $\kappa (1)$ refers to the thermal conductivity coefficient when collisions between grains are elastic ($\alpha =1$).

Figure 7. Plot of the (scaled) velocity conductivity coefficient $\kappa _U(\alpha )/\kappa _U(1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). Here, $\kappa _U(1)$ refers to the velocity conductivity coefficient when collisions between grains are elastic ($\alpha =1$).

We observe that the deviations of the transport coefficients from their elastic forms are in general significant, especially when $m=m_g$. While the (scaled) shear viscosity and thermal conductivity coefficients exhibit a non-monotonic dependence on inelasticity, the (scaled) heat diffusive and velocity conductivity coefficients increase with increasing inelasticity, regardless of the value of the mass ratio considered. In addition, while $\eta (\alpha )<\eta (1)$, the opposite happens for the thermal conductivity since $\kappa (\alpha )>\kappa (1)$. With respect to the dependence on the mass ratio $m/m_g$, at a fixed value of the coefficient of restitution, it is quite apparent that while the (scaled) shear viscosity increases with increasing mass ratio, the (scaled) thermal conductivity decreases with increasing mass ratio. The same happens for the (scaled) coefficients $n\bar {\mu }/T\kappa (1)$ and $\kappa _U(\alpha )/\kappa _U(1)$ since both scaled coefficients decrease as the mass ratio increases. We also see that in the case $m/m_g=50$, the results derived here for $\eta (\alpha )/\eta (1)$, $\kappa (\alpha )/\kappa (1)$ and $n\bar {\mu }/T\kappa (1)$ practically coincide with those obtained in the Brownian limit by Gómez González & Garzó (Reference Gómez González and Garzó2019). However, in the case $m/m_g=50$, the (scaled) velocity conductivity coefficient $\kappa _U(\alpha )/\kappa _U(1)$ ($\kappa _U=0$ for any value of $\alpha$ in the Brownian limit) is still clearly different from zero.

Although the results obtained here for the (scaled) transport coefficients depend on the values of the mass ratio and the (reduced) temperature of the molecular gas, it is worthwhile comparing the present results with those obtained for dry granular gases (i.e. in the absence of the molecular gas). In the case of the shear viscosity, a comparison between both systems (with and without the gas phase) shows significant discrepancies (see, for instance, figure 3.1 of Garzó Reference Garzó2019) even at a qualitative level since while $\eta$ increases with inelasticity for dry granular gases, the opposite happens here whatever the mass ratio considered. On the other hand, a more qualitative agreement is found for the the thermal conductivity (see, for instance, figure 3.2 of Garzó Reference Garzó2019) since $\kappa$ increases with decreasing $\alpha$ in both systems. In any case, important quantitative differences appear at strong dissipation since the influence of inelasticity on $\kappa$ is more relevant in the dry case than in the presence of the molecular gas. A similar conclusion is reached for the heat diffusive coefficient $\bar {\mu }$ (see, for instance, figure 3.3 of Garzó Reference Garzó2019) where the magnitude of this (scaled) coefficient for dry granular gases is much larger than that found here for granular suspensions.

6. Linear stability analysis of the HSS

Once the transport coefficients of the granular gas are known, the corresponding Navier–Stokes hydrodynamic equations can be explicitly displayed. To derive them, one has to take into account first that the the production of momentum $\boldsymbol {\mathcal {F}}[f]$ (defined in (2.12)) to first order in spatial gradients can be written as

(6.1)\begin{align} \boldsymbol{\mathcal{F}}^{(1)}[f^{(1)}] &= \int {\rm d}\boldsymbol{v} m\boldsymbol{V} \{J_g[f^{(0)},f_g^{(1)}]+J_g[f^{(1)},f_g^{(0)}]\} \nonumber\\ &={-}\xi {\rm \Delta} \boldsymbol{U}+\int {\rm d}\boldsymbol{v} m\boldsymbol{V}J_g[f^{(1)},f_g^{(0)}], \end{align}

where

(6.2)\begin{equation} \xi=\rho \mu \theta^{{-}1/2}(1+\theta)^{1/2}\gamma. \end{equation}

By symmetry reasons, the contributions of the first-order distribution $f^{(1)}(\boldsymbol {V})$ to the second integral in (6.1) come from the terms $\boldsymbol {\mathcal {A}}(\boldsymbol {V})$, $\boldsymbol {\mathcal {B}}(\boldsymbol {V})$ and $\boldsymbol {\mathcal {E}}(\boldsymbol {V})$. In the leading Sonine approximation (see the supplementary material), one gets

(6.3)\begin{align} \int {\rm d}\boldsymbol{v} m\boldsymbol{V}J_g[f^{(1)},f_g^{(0)}] &={-}\frac{2}{d(d+2)} \left(\frac{m}{nT^{2}}\kappa \boldsymbol{\nabla} \ln T+\frac{m}{T^{3}}\bar{\mu} \boldsymbol{\nabla} \ln n +\frac{m}{nT^{3}}\kappa_U {\rm \Delta} \boldsymbol{U}\right) \nonumber\\ &\quad \times \int {\rm d}\boldsymbol{v} m \boldsymbol{V} \boldsymbol{\cdot} J_g[\boldsymbol{S}f^{(0)},f_g^{(0)}]. \end{align}

The collision integral appearing in (6.3) can be computed with the result

(6.4)\begin{equation} \int {\rm d}\boldsymbol{v} m \boldsymbol{V}\boldsymbol{\cdot} J_g[\boldsymbol{S}f^{(0)},f_g^{(0)}]={-}\frac{d}{2}n T^{2} \mu \gamma \theta^{{-}1/2}\left(1+\theta\right)^{{-}1/2}. \end{equation}

Substitution of (6.4) into (6.3) leads to the final expression of $\boldsymbol {\mathcal {F}}^{(1)}[f^{(1)}]$:

(6.5)\begin{equation} \boldsymbol{\mathcal{F}}^{(1)}[f^{(1)}]={-}\xi {\rm \Delta} \boldsymbol{U}+\frac{\rho}{d+2}\mu\gamma\left(\frac{\kappa}{n}\boldsymbol{\nabla} \ln T+\frac{\bar{\mu}}{T}\boldsymbol{\nabla} \ln n+\frac{\kappa_U}{nT}{\rm \Delta} \boldsymbol{U}\right)X, \end{equation}

where

(6.6)\begin{equation} X(\theta)=\theta^{{-}1/2}\left(1+\theta\right)^{{-}1/2}. \end{equation}

Thus, when the constitutive equations (4.8)–(4.9) and (6.5) are substituted into the (exact) balance equations (2.8)–(2.10), one gets the Navier–Stokes hydrodynamic equations for a granular gas immersed in a molecular gas:

(6.7)\begin{gather} D_t n+n\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0, \end{gather}

(6.8)\begin{gather} \hspace{-3pc}\rho D_t U_i+\frac{\partial p}{\partial r_i} = \frac{\partial}{\partial r_j} \left[\eta\left(\frac{\partial U_j}{\partial r_i}+\frac{\partial U_i}{\partial r_j}-\frac{2}{d}\delta_{ij} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}\right)\right]-\xi {\rm \Delta} U_i \nonumber\\ \hspace{5pc}+\,\frac{\rho}{d+2}\mu\gamma X\left(\frac{\kappa}{n}\frac{\partial \ln T}{\partial r_i}+ \frac{\bar{\mu}}{T}\frac{\partial \ln n}{\partial r_i}+\frac{\kappa_U}{nT}{\rm \Delta} U_i\right), \end{gather}

(6.9)\begin{gather} \hspace{-8.9pc}D_tT+T(\zeta^{(0)}+\zeta_g^{(0)}) = \frac{2}{dn}\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\kappa\boldsymbol{\nabla} T+\bar{\mu}\boldsymbol{\nabla} n+\kappa_U {\rm \Delta} \boldsymbol{U}\right) \nonumber\\ \hspace{10pc} +\,\frac{2}{dn}\left[\eta\left(\frac{\partial U_j}{\partial r_i}+ \frac{\partial U_i}{\partial r_j}-\frac{2}{d}\delta_{ij} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}\right) \frac{\partial U_i}{\partial r_j}-\frac{2}{d}T\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{U}\right]. \end{gather}

As said in § 5, we have not considered in (6.9) the first-order contributions to $\zeta$ and $\zeta _g$ since they vanish when non-Gaussian corrections to the distribution function $f^{(0)}$ are neglected. In addition, as already mentioned in several previous works (Garzó Reference Garzó2005; Garzó, Montanero & Dufty Reference Garzó, Montanero and Dufty2006), the above production rates should also include second-order contributions in spatial gradients. However, in the case of a dry dilute granular gas (Brey et al. Reference Brey, Dufty, Kim and Santos1998), it has been shown that these contributions are very small, and hence they can be neglected in the hydrodynamic equations. We expect here that the same happens for a granular suspension. Apart from the above approximations, the Navier–Stokes hydrodynamic equations (6.7)–(6.9) are exact to second order in the spatial gradients of $n$, $\boldsymbol {U}$ and $T$.

A simple solution of (6.7)–(6.9) corresponds to the HSS studied in § 2. A natural question is whether actually the HSS may be unstable with respect to long enough wavelength perturbations, as occurs for dry granular fluids in freely cooling flows (Goldhirsch & Zanetti Reference Goldhirsch and Zanetti1993; McNamara Reference McNamara1993). This is one of the most characteristic features of granular gases; its origin is associated with the inelasticity of collisions. On the other hand, the stability of the HSS was also analysed by Gómez González & Garzó (Reference Gómez González and Garzó2019) in the Brownian limit case showing that the HSS is always linearly stable. The objective of this section is to check whether the HSS is still linearly stable for arbitrary values of the mass ratio $m/m_g$.

Given that the present analysis is quite similar to that previously made by Gómez González & Garzó (Reference Gómez González and Garzó2019), only the final results are provided here. Some mathematical steps are displayed in the supplementary material. In Fourier space, as expected (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Garzó Reference Garzó2005), the $d-1$ transverse velocity components $\boldsymbol {w}_{\boldsymbol {k}\perp }=\boldsymbol {w}_{\boldsymbol {k}}-(\boldsymbol {w}_{\boldsymbol {k}} \boldsymbol {\cdot }\hat {\boldsymbol {k}})\hat {\boldsymbol {k}}$ (orthogonal to the wave vector $\boldsymbol {k}$) decouple from the other three modes. In terms of the dimensionless time $\tau$, their evolution equation is $\partial _\tau \boldsymbol {w}_{\boldsymbol {k}\perp }- \lambda _\perp (k)\boldsymbol {w}_{\boldsymbol {k}\perp }=0$. The explicit form of $\lambda _\perp (k)$ can be found in the supplementary material. A systematic analysis of the dependence of $\lambda _\perp (k)$ on the parameter space of the system shows that $\lambda _\perp (k)$ is always negative, and hence the transversal shear modes $\boldsymbol {w}_{\boldsymbol {k}\perp }(\boldsymbol {k},\tau )$ are linearly stable.

A careful analysis of the corresponding eigenvalues $\lambda _\ell (k)$ of the matrix $\boldsymbol{\mathsf{M}}$ obeying the remaining three longitudinal modes ($\rho _{\boldsymbol {k}}$, $\theta _{\boldsymbol {k}}$ and the longitudinal velocity component of the velocity field, $w_{\boldsymbol {k}\parallel }=\boldsymbol {w}_{\boldsymbol {k}}\boldsymbol {\cdot }\hat {\boldsymbol {k}}$) shows that one of the longitudinal modes could be unstable for values of the wavenumber $k< k_{\textrm {h}}$. The expression of $k_{\textrm {h}}$ is

(6.10)\begin{equation} k_{{h}}^{2}=\frac{\sqrt{2}(2\overline{\zeta}_g\gamma^{*}-\zeta^{*})+ \dfrac{2d}{d+2}\mu X \gamma^{*}[2\left(\zeta^{*}D_T^{*}-\overline{\zeta}_g \overline{\mu}^{*}\gamma^{*}\right)-\zeta^{*}\overline{\mu}^{*}]}{\overline{\mu}^{*}-D_T^{*}}, \end{equation}

where all the quantities appearing in (6.10) are defined in the supplementary material. As in the case of $\lambda _\perp (k)$, a study of the dependence of $k_{{h}}$ on the parameters of the system shows that $k_{{h}}^{2}$ is always negative. Consequently, there are no physical values of the wavenumber for which the longitudinal modes become unstable, and hence the longitudinal modes are also linearly stable.

In summary, the linear stability analysis of the HSS carried out here for a dilute granular gas surrounded by a molecular gas shows no surprises relative to the earlier study performed in the Brownian limit: the HSS is linearly stable for arbitrary values of the mass ratio $m/m_g$. However, the dispersion relations defining the dependence of the eigenvalues $\lambda _\perp (k)$ and $\lambda _{||}(k)$ on the parameter space are very different from those previously obtained when $m/m_g\to \infty$ (Gómez González & Garzó Reference Gómez González and Garzó2019). As an illustration, figure 8 shows the real parts of the eigenvalues $\lambda _{i,||}$ ($i=1,2,3$) and $\lambda _{\perp }$ as a function of the wavenumber $k$ for $\phi =0.001$, $T_g^{*}=1000$, $m/m_g=1$ and $\alpha =0.8$. It is quite apparent that all the eigenvalues are negative, as expected. In particular, although the longitudinal mode $\lambda _{1,||}$ is quite close to 0, the inset clearly shows that it is always negative. In addition, we also observe that in general the eigenvalues exhibit a very weak dependence on $k$; this contrasts with the results obtained for dry granular fluids (see, for instance, figure 4.7 of Garzó Reference Garzó2019).

Figure 8. Dispersion relations for a three-dimensional granular gas with $\phi =0.001$, $T_g^{*}=1000$, $m/m_g=1$ and $\alpha =0.8$. From top to bottom, the curves correspond to the longitudinal mode $\lambda _{1,||}$, the two degenerate shear (transversal) modes $\lambda _\perp$ (dotted line) and the two remaining longitudinal modes $\lambda _{3,||}$ and $\lambda _{2,||}$. The dependence of $\lambda _{1,||}$, $\lambda _{3,||}$ and $\lambda _\perp$ on $k$ is shown more clearly in the insets. Only the real parts of the eigenvalues are plotted.

7. Summary and concluding remarks

The main goal of this paper has been to determine the Navier–Stokes–Fourier transport coefficients of a granular gas (modelled as a gas of inelastic hard spheres) immersed in a bath of elastic hard spheres (molecular gas). We have been interested in a situation where the solid particles are sufficiently dilute, and hence one can assume that the state of the bath (molecular gas) is not affected by the presence of the granular gas. Under these conditions, the molecular gas can be considered as a thermostat kept at equilibrium at a temperature $T_g$. This system (granular gas thermostatted by a gas of elastic hard spheres) was originally proposed years ago by Biben et al. (Reference Biben, Martin and Piasecki2002) and it can be considered as an idealised collisional model for particle-laden suspensions. Thus, in contrast to the previous suspension models employed in the granular literature (Tsao & Koch Reference Tsao and Koch1995; Sangani et al. Reference Sangani, Mo, Tsao and Koch1996; Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012; Saha & Alam Reference Saha and Alam2017) where the effect of the interstitial fluid on grains is accounted for via an effective fluid–solid force, the model considered here takes into account not only the inelastic collisions among grains themselves but also the elastic collisions between particles of the granular and molecular gases. Moreover, we have also assumed that the volume fraction occupied by the suspended solid particles is very small (low-density regime). In this case, the one-particle velocity distribution function $f(\boldsymbol {r}, \boldsymbol {v};t)$ of grains verifies the Boltzmann kinetic equation.

Before analysing inhomogeneous states, we have considered first homogeneous situations. The study of this state is important because it plays the role of the reference base state in the Chapman–Enskog solution to the Boltzmann equation. In this simple situation, the granular ‘temperature’ $T$ tends to thermalise to the bath temperature $T_g$. In the steady state, both competing effects (characterised by the cooling rates $\zeta$ and $\zeta _g$) cancel each other and a breakdown of energy equipartition appears ($T< T_g$). In the HSS, the temperature ratio $T/T_g$ and the kurtosis $a_2$ (measuring the deviation of the distribution function from its Maxwellian form) have been here estimated by considering the so-called first Sonine approximation (3.2) to the distribution function $f(\boldsymbol {v})$. In this approximation, the temperature ratio is obtained by solving the steady-state condition $\zeta +\zeta _g=0$ while the kurtosis is given by (3.19). These equations provide the dependence of $T/T_g$ and $a_2$ on the parameter space of the system. Our theoretical results extend to arbitrary dimensions the results obtained by Santos (Reference Santos2003) for hard spheres ($d=3$). To assess the accuracy of the (approximate) analytical results, a suite of Monte Carlo simulations have been also performed. Comparison between theory and simulations shows in general a very good agreement, especially in the case of the temperature ratio.

Once the homogeneous state is characterised, the next step has been to solve the Boltzmann equation by means of the Chapman–Enskog-like expansion to first order in spatial gradients (Chapman & Cowling Reference Chapman and Cowling1970; Brilliantov & Pöschel Reference Brilliantov and Pöschel2004; Garzó Reference Garzó2019). The Navier–Stokes–Fourier transport coefficients have been explicitly determined by considering the leading terms in a Sonine polynomial expansion when steady-state conditions apply. Their forms are given by (5.1) for the shear viscosity $\eta$, (5.4) for the thermal conductivity $\kappa$, (5.10) for the diffusive heat conductivity $\bar {\mu }$ and (5.11) for the so-called velocity conductivity coefficient $\kappa _U$. This latter coefficient takes into account a contribution to the heat flux coming from the velocity difference ${\rm \Delta} \boldsymbol {U}$. It is quite apparent that the expressions for the transport coefficients show a complex dependence on $\alpha$, $m/m_g$, $T_g^{*}$ and $\phi$ (see figures 4–7). The dependence on the latter two parameters appears via the dimensionless drift coefficient $\gamma ^{*}$; this coefficient provides a characteristic rate for the elastic collisions between granular and bath particles. In general, we find that significant quantitative differences between dry granular gases (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Garzó & Dufty Reference Garzó and Dufty1999; Garzó Reference Garzó2013, Reference Garzó2019; Gupta Reference Gupta2020) and granular suspensions appear as the inelasticity in collisions increases. Thus, the impact of the gas phase on the transport coefficients of the granular gas cannot be in general neglected.

Interestingly, in the Brownian limit ($m/m_g \to \infty$), a careful analysis shows that the expressions of $\eta$, $\kappa$ and $\bar {\mu }$ reduce to those previously derived by Gómez González & Garzó (Reference Gómez González and Garzó2019) using the Langevin-like model (2.16) for the instantaneous gas–solid force. In this limiting case, the present results show that the coefficient $\kappa _U$ vanishes, in agreement with Gómez González & Garzó (Reference Gómez González and Garzó2019). Therefore, the results reported in this paper extend to arbitrary values of the mass ratio $m/m_g$ the results derived in previous works (Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012; Gómez González & Garzó Reference Gómez González and Garzó2019). Nonetheless, the convergence to the Brownian limit is reached for relatively small values of the mass ratio ($m/m_g\simeq 50$ for $T_g^{*}=1000$). Thus, from a practical point of view, it seems that in most cases the Langevin-like model is able to assess the effect of the interstitial gas on transport properties of granular gas. As discussed in § 3.2, the convergence of the present results to the Brownian limit depends on the way of scaling the variables. Here, since we have been interested in recovering the results provided by Gómez González & Garzó (Reference Gómez González and Garzó2019) in the Brownian limit, $T_g^{*}$ has been used as the control parameter regarding the force exerted by the molecular gas on grains. This fact implies that the influence of the mass ratio $m/m_g$ appearing in $\gamma$ (see (2.15)) on the thermalisation process is absorbed in the selection of $T_g^{*}$. On the other hand, following the scaling proposed by Biben et al. (Reference Biben, Martin and Piasecki2002) and Santos (Reference Santos2003), another possibility could have been to choose $\omega$ (see (3.20)) as the bath parameter. In such a way, an extra dependence on the mass ratio in the coefficient $\gamma ^{*}$ emerges and leads to a different dependence of both the temperature ratio (see figure 2) and the transport coefficients on the mass ratio.

Knowledge of the forms of the transport coefficients opens up the possibility of performing a linear stability analysis on the resulting continuum hydrodynamic equations. As in the Brownian limit case (Gómez González & Garzó Reference Gómez González and Garzó2019), the analysis shows that the HSS is always linearly stable whatever the mass ratio considered.

Although in real experiments the energy input for keeping rapid-flow conditions is usually done either by driving through the boundaries (Yang et al. Reference Yang, Huan, Candela, Mair and Walsworth2002) or by bulk driving (as in air-fluidised beds; Schröter, Goldman & Swinney Reference Schröter, Goldman and Swinney2005; Abate & Durian Reference Abate and Durian2006), these ways of supplying energy produce in many cases strong spatial gradients in the bulk domain, and so the Navier–Stokes hydrodynamics does not hold. To overcome this difficulty, it is quite common in computer simulations (see e.g. Puglisi et al. Reference Puglisi, Loreto, Marconi, Petri and Vulpiani1998; Paganobarraga et al. Reference Paganobarraga, Trizac, van Noije and Ernst2002; Prevost, Egolf & Urbach Reference Prevost, Egolf and Urbach2002; Fiege, Aspelmeier & Zippelius Reference Fiege, Aspelmeier and Zippelius2009; Shaebani, Sarabadani & Wolf Reference Shaebani, Sarabadani and Wolf2013) to heat the system homogeneously by the presence of an external driving force or thermostat (Evans & Morriss Reference Evans and Morriss1990). A different (and likely more realistic) way of thermostatting a granular gas is by means of a sea of elastic hard spheres (Biben et al. Reference Biben, Martin and Piasecki2002). In this context, given that in most of the computer simulation works the effects of the thermostat on the properties of the granular gas are ignored, the results derived in this paper may be useful for simulators when studying problems in granular fluids thermostatted by a bath of elastic hard spheres. Regardless of practical applications, needless to say a complete comprehension of gas–solid flows is still missing (Fullmer & Hrenya Reference Fullmer and Hrenya2017; Morris Reference Morris2020; Han et al. Reference Han, Yue, Geng, Hu, Liu, Suleiman, Cui, Bai and Xu2021). For this reason, the development of theoretical models for granular suspensions where collisions play a significant role becomes of great relevance to understand from a more fundamental view the results derived by means of effective models such as the Langevin-like model (Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012).

One of the main limitations of the results derived in this paper is its restriction to the low-density regime. The extension of the present theory to a moderately dense granular suspension described by the Enskog kinetic equation is an interesting project for the future. These results could stimulate the performance of MD simulations to assess the reliability of the theory for finite densities. Furthermore, in an attempt to model gas–solid flows in fluidised beds where the fluid–solid interactions are governed not only by the drag force (Stokes drag force) but also by the Archimedes force (gas pressure gradient), the present model could be also extended to account for these new ‘additional forces’. Another challenging work could be the determination of the non-Newtonian rheological properties of a granular suspension under simple shear flow. This study would allow the extension of previous studies (Tsao & Koch Reference Tsao and Koch1995; Sangani et al. Reference Sangani, Mo, Tsao and Koch1996; Chamorro et al. Reference Chamorro, Vega Reyes and Garzó2015; Saha & Alam Reference Saha and Alam2017; Alam et al. Reference Alam, Saha and Gupta2019; Takada et al. Reference Takada, Hayakawa, Santos and Garzó2020) to arbitrary values of the mass ratio $m/m_g$. The accuracy of the results for the shear viscosity can be tested by means of DSMC simulations. For instance, by performing local deviations of the velocity field from its value in the homogeneous state and analysing the relaxation process (Brey, Ruiz-Montero & Cubero Reference Brey, Ruiz-Montero and Cubero1999b; Brey & Cubero Reference Brey and Cubero2001). Another possible project could be to revisit the results obtained in this paper by considering the charge transport equation recently considered by Ceresiat, Kolehmainen & Ozel (Reference Ceresiat, Kolehmainen and Ozel2021). Work along these lines will be carried out in the future.