2. Granular gas in contact with a bath of elastic hard spheres. Boltzmann kinetic description

We consider a gas of inelastic hard disks ( $d=2$ ) or spheres ( $d=3$ ) of mass $m$ , diameter $\sigma$ and coefficient of normal restitution $\alpha$ . We assume that the spheres are perfectly smooth and therefore the collisions between particles are inelastic and characterised by a (positive) constant coefficient of normal restitution $\alpha \leqslant 1$ . For elastic collisions $\alpha =1$ , while $\alpha \lt 1$ for inelastic collisions. The granular gas is immersed in a molecular gas consisting of hard disks or spheres of mass $m_g$ and diameter $\sigma _g$ . The collisions between the granular particles and the gas molecules are assumed to be elastic. As mentioned in § 1, we also assume that the number density of the granular particles is much smaller than that of the molecular gas, so that the state of the latter is not significantly affected by the presence of grains. In this sense, the molecular gas can be treated as a thermostat or bath in equilibrium at temperature $T_g$ . Thus, its velocity distribution function $f_g(\boldsymbol{V}_{\!g})$ is

(2.1) \begin{equation} f_g(\boldsymbol{V}_{\!g})=n_g \Big (\frac {m_g}{2\pi T_g}\Big )^{d/2} \exp \Bigg (-\frac {m_g V_g^2}{2T_g}\Bigg ), \end{equation}

where $n_g$ is the number density of molecular gas and $\boldsymbol{V}_{\!g}=\boldsymbol{v}-\boldsymbol{U}_{\!g}$ . In principle, the mean flow velocity of molecular gas $\boldsymbol{U}_{\!g}$ is different from the mean flow velocity $\boldsymbol{U}$ of solid particles (see its definition in (2.8)). In addition, for the sake of simplicity, we take the Boltzmann constant $k_{{B}}=1$ throughout the paper.

In the low-density regime, the velocity distribution function $f(\boldsymbol{r}, \boldsymbol{v}, t)$ of granular particles verifies the Boltzmann kinetic equation. Moreover, although the granular gas is sufficiently rarefied and hence the properties of the molecular (interstitial) gas can be supposed to be constant, one has to take into account the collisions among grains themselves in the kinetic equation of $f(\boldsymbol{r}, \boldsymbol{v}, t)$ . Thus, in the presence of the gravitational field $\boldsymbol{g}$ , the distribution $f$ verifies the Boltzmann equation:

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

Here, the Boltzmann collision operator $J[f,f]$ gives the rate of change of $f$ due to inelastic collisions among granular particles. Its explicit form is (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004; Garzó Reference Garzó2019)

(2.3) \begin{equation} 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}) \left [\alpha ^{-2}f({\boldsymbol v}_{1}^{\prime\prime})f({\boldsymbol v}_{2}^{\prime\prime})-f({\boldsymbol v}_{1})f({\boldsymbol v}_{2})\right ]\!, \end{equation}

where $\boldsymbol{g}_{12}=\boldsymbol{v}_1-\boldsymbol{v}_2$ is the relative velocity, $\widehat {\boldsymbol{\sigma }}$ is a unit vector that joins the centres of the colliding particles and $\varTheta$ is the Heaviside step function. In (2.3), the double primes denote pre-collisional velocities. The relation between them and their corresponding post-collisional velocities $(\boldsymbol{v}_1,\boldsymbol{v}_2)$ is

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

In (2.2), the Boltzmann–Lorentz operator $J_g[f,f_g]$ accounts for the rate of change of $f$ due to elastic collisions between particles of the granular and molecular gas. Its form is (Résibois & de Leener Reference Résibois and de Leener1977)

(2.5) \begin{equation} 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}) \left [f({\boldsymbol v}_{1}^{\prime\prime})f_g({\boldsymbol v}_{2}^{\prime\prime})-f({\boldsymbol v}_{1})f_g({\boldsymbol v}_{2})\right ]\!, \end{equation}

where $\overline {\sigma }=(\sigma +\sigma _g)/2$ and the relationship between $(\boldsymbol{v}_{1}^{\prime\prime},\boldsymbol{v}_{2}^{\prime\prime})$ and $(\boldsymbol{v}_1,\boldsymbol{v}_2)$ is

(2.6) \begin{equation} \boldsymbol{v}_{1}^{\prime\prime}=\boldsymbol{v}_1-2 \mu _g(\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }}, \quad \boldsymbol{v}_{2}^{\prime\prime}=\boldsymbol{v}_2+2 \mu (\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }}. \end{equation}

Here,

(2.7) \begin{equation} \mu _g=\frac {m_g}{m+m_g}, \quad \mu =\frac {m}{m+m_g}. \end{equation}

As is customary, the effect of gravity on the properties of molecular gases under ordinary conditions is neglected in the present analysis. This approximation is justified by the fact that the influence of gravity on a molecule between successive collisions is negligible, i.e. $\ell \ll h$ , where $\ell$ denotes the mean free path for hard spheres and $h = v_{{\rm th}}^2/g$ represents the characteristic length scale over which gravitational effects become significant ( $v_{{\rm th}}$ being the thermal velocity). For instance, under terrestrial conditions at room temperature, this ratio is of the order of $\ell /h \sim 10^{-11}$ (Tij, Garzó & Santos Reference Tij, Garzó and Santos1999), thereby validating the omission of gravitational effects in the description of molecular gas behaviour.

The number density $n$ , mean flow velocity $\boldsymbol{U}$ and granular temperature $T$ of the granular gas are defined as the first few velocity moments of $f$ :

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

where $\boldsymbol{V}=\boldsymbol{v}-\boldsymbol{U}$ is the peculiar velocity. As said before, the difference $\Delta \boldsymbol{U}=\boldsymbol{U}-\boldsymbol{U}_{\!g}$ is in general different from zero (Gómez González & Garzó Reference Gómez González and Garzó2022b ). In fact, as we show later, $\Delta \boldsymbol{U}$ induces a non-vanishing contribution to the mass flux of intruders.

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 (Gómez González & Garzó Reference Gómez González and Garzó2022b )

(2.9) \begin{align}&\qquad\qquad D_t n+n\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}=0, \end{align}

(2.10) \begin{align}&\qquad \rho D_t\boldsymbol{U}=-\boldsymbol{\nabla }\boldsymbol{\cdot }\textsf{P}+\rho \boldsymbol{g}+\mathcal{F}[f], \end{align}

(2.11) \begin{align}& D_tT+\frac {2}{dn}\big (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}+\textsf{P}:\boldsymbol{\nabla }\boldsymbol{U}\big )=- T \zeta -T \zeta _g. \\[9pt] \nonumber \end{align}

In (2.9)–(2.11), $D_t=\partial _t+\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla}$ is the material derivative and the pressure tensor $\textsf{P}$ and the heat flux vector $\boldsymbol{q}$ are given, respectively, as

(2.12) \begin{equation} \textsf{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}

The production of momentum term $\mathcal{F}[f]$ appearing in (2.10) is defined as

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

This term is in general different from zero since the Boltzmann–Lorentz collision operator $J_g[f,f_g]$ does not conserve momentum. The form of $\mathcal{F}[f]$ can be made more explicit when one takes into account the property (Brilliantov & Pöschel Reference Brilliantov and Pöschel2004; Garzó Reference Garzó2019)

(2.14) \begin{eqnarray} \int {\rm d}\boldsymbol{v}_1 \varPsi (\boldsymbol{v}_1) J_g[\boldsymbol{v}_2|f,f_g]&=&\overline {\sigma }^{d-1}\int {\rm d}\boldsymbol{v}_{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)\nonumber \\ & & \times \left [\varPsi (\boldsymbol{v}_{1}^{\prime})-\varPsi (\boldsymbol{v}_1)\right ]\!, \end{eqnarray}

where $\boldsymbol{v}_{1}^{\prime}=\boldsymbol{v}_1-2 \mu _g(\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }}$ . Using (2.14), $\mathcal{F}[f]$ is

(2.15) \begin{equation} \mathcal{F}[f]=-\frac {2\pi ^{(d-1)/2}}{\varGamma \Big (\frac {d+3}{2}\Big )} \frac {m m_g}{m+m_g} \overline {\sigma }^{d-1}\int {\rm d}\boldsymbol{v}_1\int {\rm d}\boldsymbol{v}_2\; g_{12}\boldsymbol{g}_{12}\; f(\boldsymbol{v}_1) f_g(\boldsymbol{v}_2). \end{equation}

Finally, the partial production rates $\zeta$ and $\zeta _g$ appearing in the balance equation (2.11) are given, respectively, as

(2.16) \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}

While the cooling rate $\zeta$ provides the rate of change of kinetic energy of grains due to their inelastic collisions, the term $\zeta _g$ gives the transfer of kinetic energy in the collisions between the particles of the molecular and granular gases. The quantity $\zeta =0$ for elastic collisions ( $\alpha =1$ ) while $\zeta _g=0$ when the particles of the molecular and granular gases are mechanically equivalent.

It is interesting at this point to note the meaning of the granular temperature $T$ . To understand this, it is important to remember that our study is limited to the so-called rapid flow regime, namely a situation where grains are subjected to a strong external excitation (e.g. vibrating or shearing walls or air-fluidised beds). In this regime, the external energy supplied to the granular gas can compensate for the energy loss due to collisions and the effects of gravity. Since in this regime the motion of the grains is quite similar to the chaotic motion of atoms or molecules in an ordinary gas, as discussed in previous works (Gómez González & Garzó Reference Gómez González and Garzó2022b ), it is tempting to establish a relationship between the statistical motion of the grains and some kind of temperature. In this context, as usual in conventional kinetic theory (Chapman & Cowling Reference Chapman and Cowling1970), the granular temperature $T$ can be interpreted as a measure of the fluctuations of the velocities of grains with respect to its mean value $\boldsymbol{U}$ . Since granular gases are athermal systems (i.e. their thermal fluctuations have a negligible effect on the dynamics of grains), the granular temperature $T$ has no thermodynamic interpretation in contrast to the temperature $T_g$ of the molecular gas (see, for example, the review paper of Goldhirsch (Reference Goldhirsch2008) for a discussion of this issue). In any case, within the context of statistical thermodynamics, the thermodynamic temperature $T_g$ can be also understood as a statistical quantity measuring the deviations of molecular particle velocity $\boldsymbol{v}$ from its mean value $\boldsymbol{U}_{\!g}$ .

Before closing this section, it is worth analysing the limiting case $m\gg m_g$ . It is in fact a quite realistic case for granular suspensions (Subramaniam Reference Subramaniam2020) where the particles of the interstitial molecular gas are much lighter than the particles of the granular gas. In the limit $m/m_g\to \infty$ , a Kramers–Moyal expansion (Résibois & de Leener Reference Résibois and de Leener1977) allows us to approximate the Boltzmann–Lorentz operator $J_g[f,f_g]$ to the Fokker–Planck operator $J_g^{{FP}}[f,f_g]$ :

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

where the friction coefficient $\gamma$ is

(2.18) \begin{equation} \gamma =\frac {4\pi ^{(d-1)/2}}{{\rm d}\varGamma \Big (\frac {d}{2}\Big )}\Big (\frac {m_g}{m}\Big )^{1/2}\Bigg (\frac {2T_g}{m}\Bigg )^{1/2}n_g \overline {\sigma }^{d-1}. \end{equation}

Upon deriving (2.17), it has been assumed that $\boldsymbol{U}_{\!g}=\boldsymbol{0}$ and that the distribution function $f$ of the granular gas is a Maxwellian distribution.

Most of the theoretical works for suspension models reported in the granular literature are essentially based on the use of the Fokker–Planck operator (2.17) to account for in an effective way (coarse-grained approach) the influence of the surrounding fluid on the dynamics of grains (Koch & Hill Reference Koch and Hill2001; Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012; Chassagne et al. Reference Chassagne, Bonamy and Chauchat2023; Garzó Reference Garzó2023). These sorts of models have been considered to obtain the Navier–Stokes–Fourier transport coefficients of the suspension (Gómez González & Garzó Reference Gómez González and Garzó2019).

3. Intruders in granular suspensions

We assume now that a few intruders (or tracers) of mass $m_0$ and diameter $\sigma _0$ are added to the granular gas. In this situation, intruders and particles of the granular gas are surrounded by the molecular gas (bath of elastic hard spheres). The system can be seen as a ternary mixture where one of the components (intruders) is present in tracer concentration. Apart from the restitution coefficient $\alpha$ for inelastic grain–grain collisions, the coefficient of normal restitution $\alpha _0\leqslant 1$ characterises the inelastic collisions between the intruders and the particles of the granular gas. As in the case of the granular gas, collisions between intruders and particles of the surrounding molecular gas are elastic.

Since the concentration of intruders is much smaller than that of the granular gas (tracer limit), their presence does not affect the state of the granular gas. Under these conditions and in the presence of the gravitational field, the velocity distribution function $f_0(\boldsymbol{r}, \boldsymbol{v}, t)$ of the intruders obeys the kinetic equation:

(3.1) \begin{equation} \frac {\partial f_0}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }f_0+\boldsymbol{g}\boldsymbol{\cdot }\frac {\partial f_0}{\partial \boldsymbol{v}}=J_0[f_0,f]+J_{0g}[f_0,f_g], \end{equation}

where the (inelastic) version of the Boltzmann–Lorentz collision operator $J_0[f_0,f]$ gives the rate of change of $f_0$ due to the inelastic collisions between the intruders and particles of the granular gas. It is given by (Garzó Reference Garzó2019)

(3.2) \begin{equation} J_0[\boldsymbol{v}_1|f_0,f]=\sigma ^{\prime (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}) \left [\alpha _0^{-2}f_0({\boldsymbol v}_{1}^{\prime\prime})f({\boldsymbol v}_{2}^{\prime\prime})-\! f_0({\boldsymbol v}_{1})f({\boldsymbol v}_{2})\right ]\!. \end{equation}

As in (2.3), $\boldsymbol{g}_{12}=\boldsymbol{v}_1-\boldsymbol{v}_2$ is the relative velocity, $\widehat {\boldsymbol{\sigma }}$ is a unit vector and $\varTheta$ is the Heaviside step function. In addition, $\sigma '=(\sigma +\sigma _0)/2$ ,

(3.3) \begin{equation} \boldsymbol{v}_{1}^{\prime\prime}=\boldsymbol{v}_1-\frac {1+\alpha _0}{\alpha _0}\mu '(\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }}, \quad \boldsymbol{v}_{2}^{\prime\prime}=\boldsymbol{v}_2+\frac {1+\alpha _0}{2\alpha _0}\mu _{0}^{\prime}(\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }} \end{equation}

and

(3.4) \begin{equation} \mu '=\frac {m}{m+m_0}, \quad \mu _{0}^{\prime}=\frac {m_0}{m+m_0}. \end{equation}

In (3.1), the collision operator $J_{0g}[f_0,f_g]$ provides the rate of change of $f_0$ due to elastic collisions between intruders and particles of the molecular gas. Similarly to the operator $J_{0}[f,f_g]$ , it is given by

(3.5) \begin{equation} J_{0g}[\boldsymbol{v}_1|f_0,f_g]=\overline {\sigma }_0^{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}) \left [f_0({\boldsymbol v}_{1}^{\prime\prime})f_g({\boldsymbol v}_{2}^{\prime\prime})-f_0({\boldsymbol v}_{1})f_g({\boldsymbol v}_{2})\right ]. \end{equation}

Here, $\overline {\sigma }_0=(\sigma _0+\sigma _g)/2$ ,

(3.6) \begin{equation} \boldsymbol{v}_{1}^{\prime\prime}=\boldsymbol{v}_1-2 \mu _{g0}(\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }}, \quad \boldsymbol{v}_{2}^{\prime\prime}=\boldsymbol{v}_2+2 \mu _{0g} (\widehat {{\boldsymbol{\sigma }}} \boldsymbol{\cdot }{\boldsymbol g}_{12})\widehat {\boldsymbol{\sigma }} \end{equation}

and

(3.7) \begin{equation} \mu _{g0}=\frac {m_g}{m_0+m_g}, \quad \mu _{0g}=\frac {m_0}{m_0+m_g}. \end{equation}

Although the granular temperature $T$ is the relevant one at a hydrodynamic level, an interesting quantity at a kinetic level is the local temperature of the intruders $T_0$ . This quantity measures the mean kinetic energy of the intruders. It is defined as

(3.8) \begin{equation} T_0(\boldsymbol{r}, t)=\frac {m_0}{d n_0(\boldsymbol{r}, t)}\int \; {\rm d}\boldsymbol{v}\, V^2 f_0(\boldsymbol{r},\boldsymbol{v},t). \end{equation}

As confirmed by kinetic theory calculations (Garzó & Dufty Reference Garzó and Dufty1999b ), computer simulations (Garzó Reference Garzó2019) and experiments (Feitosa & Menon Reference Feitosa and Menon2002; Wildman & Parker Reference Wildman and Parker2002; Puzyrev et al. Reference Puzyrev, Trittel, Harth and Stannarius2024), the global temperature $T$ and the temperature of impurities $T_0$ are in general different.

Intruders may freely exchange momentum and energy with the particles of the granular and molecular gas. Thus, only the number density of intruders

(3.9) \begin{equation} n_0(\boldsymbol{r}, t)=\int {\rm d}\boldsymbol{v} f_0(\boldsymbol{r}, \boldsymbol{v}, t) \end{equation}

is conserved. This yields the balance equation

(3.10) \begin{equation} \frac {\partial \rho _0}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{j}_0=0, \end{equation}

where $\rho _0=m_0 n_0$ is the mass density of intruders and

(3.11) \begin{equation} \boldsymbol{j}_0(\boldsymbol{r}, t)= \int {\rm d}\boldsymbol{v}\; m_0 \boldsymbol{V}\; f_0(\boldsymbol{r}, \boldsymbol{v},t) \end{equation}

is the mass flux of intruders.

As in the case of the Boltzmann–Lorentz operator $J_g[f,f_g]$ , in the limiting case $m_0\gg m_g$ the operator $J_{0g}[f_0,f_g]$ reduces to the Fokker–Planck operator:

(3.12) \begin{equation} J_{0g}[f_0,f_g]\to J_{0g}^{{FP}}[f_0,f_g]=\gamma _0 \frac {\partial }{\partial \boldsymbol{v}}\boldsymbol{\cdot }\Bigg (\boldsymbol{v}+\frac {T_g}{m_0}\frac {\partial }{\partial \boldsymbol{v}}\Bigg )f_0(\boldsymbol{v}), \end{equation}

where the friction coefficient $\gamma _0$ is

(3.13) \begin{equation} \gamma _0=\frac {4\pi ^{(d-1)/2}}{d\varGamma \Big (\dfrac {d}{2}\Big )}\Big (\frac {m_g}{m_0}\Big )^{1/2}\Bigg (\dfrac {2T_g}{m_0}\Bigg )^{1/2}n_g \overline {\sigma }_0^{d-1}. \end{equation}

Note that the expression (3.13) for $\gamma _0$ differs from the macroscopic Stokes law describing the Brownian motion of a massive particle in an equilibrium host fluid (Bocquet, Hansen & Piaseck Reference Bocquet, Hansen and Piasecki1994a , Reference Bocquet, Piasecki and Hansenb ; Gómez González & Garzó Reference Gómez González and Garzó2022a ). This means that the results derived in this paper in the Brownian limit do not exactly reduce to those obtained by using a coarse-grained approach (Gómez González & Garzó Reference Gómez González and Garzó2022a ).

4. Homogeneous steady state

Before considering inhomogeneous situations, it is pertinent to analyse the HSS for the system (intruders and granular gas immersed in a molecular gas) in the absence of the gravitational field ( $\boldsymbol{g} = \boldsymbol{0}$ ). The study of this state for the intruders is crucial since the HSS plays the role of the reference state in the Chapman–Enskog method (Chapman & Cowling Reference Chapman and Cowling1970).

In the HSS, the densities $n$ and $n_0$ and the granular temperature $T$ are spatially uniform. Moreover, without loss of generality, the mean flow velocities vanish ( $\boldsymbol{U}=\boldsymbol{U}_{\!g}=\boldsymbol{0}$ ) with an appropriate choice of the frame of reference. Under these conditions, (2.2) for the granular gas reads

(4.1) \begin{equation} \partial _t f=J[f,f_0]+J_{g}[f,f_g], \end{equation}

while (3.1) for the intruders becomes

(4.2) \begin{equation} \partial _t f_0=J_0[f_0,f]+J_{0g}[f_0,f_g]. \end{equation}

4.1. Homogeneous steady state for the granular gas

Given that the HSS for the granular gas was already analysed by Gómez González & Garzó (Reference Gómez González and Garzó2022b ), only a few results are provided here. Since in the HSS the distribution $f(\boldsymbol{v})$ is isotropic in $\boldsymbol{v}$ , then $\mathcal{F}[f]=\boldsymbol{0}$ . On the other hand, the kinetic energy is not conserved by collisions and so $\zeta \neq 0$ and $\zeta _g \neq 0$ . The most interesting quantity in the HSS for the granular gas is the temperature ratio $\chi =T/T_g$ , which is in general different from 1. The only non-trivial balance equation in the homogeneous state for the granular gas is that of the temperature (2.11):

(4.3) \begin{equation} \frac {\partial T}{\partial t}=-T\left (\zeta +\zeta _g\right ). \end{equation}

After a transient period, it is expected that the granular gas achieves a steady state. Thus, according to (4.3), the steady-state condition is $\zeta +\zeta _g=0$ . As discussed by Gómez González & Garzó (Reference Gómez González and Garzó2022b ), as the molecular gas acts as a thermostat in the steady state, then the mean kinetic energy of the granular particles is smaller than that of the molecular gas ( $T\lt T_g$ ). This necessarily requires that $\zeta _g\lt 0$ so that, in the steady state, the production rates $\zeta$ and $|\zeta _g|$ exactly compensate each other and one achieves the condition $\zeta +\zeta _g=0$ .

However, to determine $\zeta$ and $\zeta _g$ one needs to know the velocity distribution function $f(\boldsymbol{v})$ for the granular gas. For inelastic collisions ( $\alpha \lt 1$ ) this distribution is not exactly known to date. On the other hand, the results obtained for the fourth cumulant or kurtosis $c$ of the distribution $f$ in Gómez González & Garzó (Reference Gómez González and Garzó2022b ) clearly show (see figure 3 of Gómez González & Garzó (Reference Gómez González and Garzó2022b )) that the magnitude of $c$ is in general very small for not quite strong inelasticity (e.g. $\alpha \gtrsim 0.5$ ). Thus, to estimate the production rates $\zeta$ and $\zeta _g$ one can replace the true distribution $f(\boldsymbol{v})$ by the Maxwellian distribution:

(4.4) \begin{equation} f_{{M}}(\boldsymbol{v})=n \Big (\frac {m}{2\pi T}\Big )^{d/2} \exp \left (-\frac {m v^2}{2T}\right )\!. \end{equation}

In the Maxwellian approximation, the dimensionless production rates $\zeta ^*=\zeta /\nu$ and $\zeta _g^*=\zeta _g/\nu$ are given by (Gómez González & Garzó Reference Gómez González and Garzó2022b )

(4.5) \begin{equation} \zeta ^*=\frac {\sqrt {2}\pi ^{(d-1)/2}}{d\varGamma \Big (\dfrac {d}{2}\Big )}(1-\alpha ^2), \quad \zeta _g^*=2x(1-x^2) \left (\frac {\mu T}{T_g}\right )^{1/2}\gamma ^*, \end{equation}

where $\nu =n\sigma ^{d-1}\sqrt {2T/m}$ is an effective collision frequency,

(4.6) \begin{equation} x=\Bigg (\mu _g+\mu \frac {T_g}{T}\Bigg )^{1/2} \end{equation}

and

(4.7) \begin{equation} \gamma ^*=\epsilon \; \chi ^{-1/2}, \quad \epsilon =\frac {\ell \gamma }{\sqrt {2T_g/m}}=\frac {\sqrt {2}\pi ^{d/2}}{2^d {d} \varGamma \left (\dfrac {d}{2}\right )}\frac {1}{\phi \sqrt {T_g^*}}. \end{equation}

Here, $\ell =1/n\sigma ^{d-1}$ is proportional to the mean free path of hard spheres,

(4.8) \begin{equation} \phi =\frac {\pi ^{d/2}}{2^{d-1}{d} \varGamma \left (\dfrac {d}{2}\right )} n\sigma ^d \end{equation}

is the solid volume fraction and

(4.9) \begin{equation} T_g^*=\frac {T_g}{m\sigma ^2 \gamma ^2}. \end{equation}

In the Maxwellian approximation (i.e. when one replaces $f$ by $f_{{M}}$ ), the steady temperature ratio $T/T_g$ can be obtained by inserting the expressions (4.5) of $\zeta ^*$ and $\zeta _g^*$ , respectively, into the (exact) steady-state condition $\zeta ^*+\zeta _g^*=0$ . This yields a cubic equation for the quantity $x$ whose physical solution is given by (4.14) of the supplementary material of Gómez González & Garzó (Reference Gómez González and Garzó2022b ). In terms of $x$ , the final expression of the temperature ratio $T/T_g$ is given by (4.15) of that supplementary material. In spite of considering the Maxwellian approximation for the distribution $f$ , an excellent agreement between theory and simulations for the temperature ratio is observed over the whole range of values of $\alpha$ studied (see figure 1 of Gómez González & Garzó (Reference Gómez González and Garzó2022b )).

4.2. Homogeneous steady state for the intruders

We analyse now the HSS for the intruders. The balance equation for the intruders’ temperature $T_0$ can be easily obtained from (4.2) as $\partial _t \ln T_0=- (\zeta _0+\zeta _{0g} )$ , where

(4.10) \begin{equation} \zeta _0=-\frac {m_0}{{d} n_0 T_0}\int {d}\boldsymbol{v}\; v^2\; J_0[\boldsymbol{v}|f_0,f], \quad \zeta _{0g}=-\frac {m_0}{{d} n_0 T_0}\int {\rm d}\boldsymbol{v}\; v^2\; J_{0g}[\boldsymbol{v}|f_0,f_g]. \end{equation}

In the HSS, $\partial _t T_0=0$ and the condition for obtaining $T_0$ is

(4.11) \begin{equation} \zeta _0+\zeta _{0g}=0. \end{equation}

As in the case of the granular gas, the exact form of the distribution function $f_0(\boldsymbol{v})$ for inelastic collisions is not known to date. The departure of $f_0(\boldsymbol{v})$ from its Maxwellian form

(4.12) \begin{equation} f_{0,{M}}(\boldsymbol{v})=n_0 \Big (\frac {m_0}{2\pi T_0}\Big )^{d/2} \exp \Big (-\frac {m_0 v^2}{2T_0}\Big ) \end{equation}

can be measured by the kurtosis $c_0$ . It is defined as (Garzó Reference Garzó2019)

(4.13) \begin{equation} c_0=\frac {1}{d(d+2)}\frac {m_0^2}{n_0T_0^2}\int {\rm d}\boldsymbol{v}\; v^4 f_0(\boldsymbol{v})-1. \end{equation}

Figure 1. Plot of the kurtosis $c_0$ associated with the distribution function of the intruders as a function of the coefficient of normal restitution $\alpha$ for $d=3$ , $\phi =0.0052$ , $T_g^*=1000$ and four different values of the mass ratio $m_0/m_g$ (from top to bottom, $m_0/m_g=20,\ 50,\ 100$ and 1000). Moreover, in all the curves $m_0/m=10$ , $\sigma _0/\sigma =5$ and $\sigma _0/\sigma _g=(m_0/m_g)^{1/3}$ . The solid lines are the theoretical results while the symbols are the DSMC simulation results. The dashed line is the result obtained from the Fokker–Planck approach (3.12) to the operator $J_{0g}[f_0,f_g]$ . Diamonds refer to DSMC simulations implemented using the time-driven approach (Gómez González & Garzó Reference Gómez González and Garzó2022b ).

Some technical details on the determination of $c_0$ are given in the supplementary material. Given that the expression of $c_0$ is very large and not very illuminating, its final form is not displayed here. In terms of dimensionless quantities, the parameter space of a $d$ -dimensional system is given by the set

(4.14) \begin{equation} \xi \equiv \left \{\frac {\sigma }{\sigma _g}, \frac {\sigma _0}{\sigma _g}, \frac {m}{m_g}, \frac {m_0}{m_g}, \alpha , \alpha _0, \phi , T_g^*\right \}. \end{equation}

In contrast to the monocomponent case (Gómez González & Garzó Reference Gómez González and Garzó2022b ), note that the diameter ratios $\sigma /\sigma _g$ and $\sigma _0/\sigma _g$ appear also as input parameters of the system.

Figure 1 shows $c_0$ versus the (common) coefficient of restitution $\alpha =\alpha _0$ for $d=3$ , $\phi =0.0052$ and $T_g^*=1000$ . Four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000) are considered keeping the ratio $m_0/m=10$ . In addition, $\sigma _0/\sigma =5$ and we have assumed that the intruders and molecular gas particles have the same mass density (i.e. $\sigma _0/\sigma =(m_0/m_g)^{1/3}$ ). As occurs for the kurtosis $c$ of the granular gas (see figure 3 of Gómez González & Garzó (Reference Gómez González and Garzó2022b )), it is quite apparent from figure 1 that the magnitude of $c_0$ is in general quite small. We observe that the agreement between theory and direct simulation Monte Carlo (DSMC) simulations is excellent even for quite extreme values of inelasticity.

The fact that $c_0$ is small allows us to guarantee that a good estimate of the production rates $\zeta _0$ and $\zeta _{0g}$ can be obtained by replacing $f_0(\boldsymbol{v})$ by the Maxwellian distribution $f_{0,{M}}(\boldsymbol{v})$ in (4.10). In this approximation, the dimensionless quantities $\zeta _0^*=\zeta _0/\nu$ and $\zeta _{0g}^*=\zeta _{0g}/\nu$ can be written as

(4.15) \begin{align}& \zeta _0^*=\frac {4\pi ^{(d-1)/2}}{{d}\varGamma \left (\frac {d}{2}\right )}\mu ' \left (\frac {\sigma '}{\sigma }\right )^{d-1}\left (1+\frac {m T_0}{m_0 T}\right )^{1/2}(1+\alpha _0)\left [1-\frac {\mu '}{2}(1+\alpha _0)\left (1+\frac {m_0 T}{m T_0}\right )\right ]\!, \end{align}

(4.16) \begin{align}&\qquad\qquad\qquad\qquad\qquad \zeta _{0g}^*=2x_0(1-x_0^2)\left (\mu _{0g}\frac {T_0}{T_g}\right )^{1/2}\gamma _0^*. \end{align}

In (4.15) and (4.16), we have introduced the quantities

(4.17) \begin{align}&\qquad x_0=\left (\mu _{g0}+\mu _{0g}\frac {T_g}{T_0}\right )^{1/2}\!, \end{align}

(4.18) \begin{align}& \gamma _0^*=\epsilon _0\chi ^{-1/2}, \quad \epsilon _0= \left (\frac {\overline {\sigma }_0}{\overline {\sigma }}\right )^{d-1}\frac {m}{m_0}\epsilon . \end{align}

The temperature ratio $\chi _0\equiv T_0/T_g$ can be finally determined by substituting (4.15) and (4.16) into the condition (4.11). It yields the nonlinear algebraic equation

(4.19) \begin{eqnarray} 2x_0(x_0^2-1)\left (\mu _{0g}\chi _0\right )^{1/2}\gamma _0^*&=&\frac {4\pi ^{(d-1)/2}}{{d}\varGamma \left (\dfrac {d}{2}\right )}\mu ' \left (\frac {\sigma '}{\sigma }\right )^{d-1}\left (1+\frac {m \chi _0}{m_0 \chi }\right )^{1/2}(1+\alpha _0)\nonumber \\ & & \times \left [1-\frac {\mu }{2}(1+\alpha _0)\left (1+\frac {m_0 \chi }{m \chi _0}\right )\right ]\!, \end{eqnarray}

where we recall that $\chi =T/T_g$ is given by (3.7) of the supplementary material of Gómez González & Garzó (Reference Gómez González and Garzó2022b ). The numerical solution to (4.17) provides the dependence of $T_0/T_g$ on the parameter space $\xi$ .

Figure 2. Temperature ratio $\chi _0\equiv T_0/T_g$ versus the (common) coefficient of normal restitution $\alpha _0=\alpha$ for $d=3$ , $\phi =0.0052$ , $T_g^*=1000$ and four different values of the mass ratio $m_0/m_g$ (from top to bottom, $m_0/m_g=20,\ 50,\ 100$ and 1000). Moreover, in all the curves $m_0/m=10$ , $\sigma _0/\sigma =5$ and $\sigma _0/\sigma _g=(m_0/m_g)^{1/3}$ . The solid lines are the theoretical results while the symbols are the DSMC simulation results. The dashed line is the result obtained by using the Fokker–Planck approach (3.12) to the operator $J_{0g}[f_0,f_g]$ . Diamonds refer to DSMC simulations implemented using the time-driven approach (Gómez González & Garzó Reference Gómez González and Garzó2022b ).

For the sake of illustration, figure 2 shows $\chi _0\equiv T_0/T_g$ as a function of the (common) coefficient of restitution $\alpha _0=\alpha$ for the same systems as in figure 1. As occurs for the ratio $T/T_g$ , due to the way of scaling the relevant quantities of the system, the deviation of $\chi _0$ from unity increases with decreasing the mass ratio $m_0/m_g$ . The agreement between the (approximate) theoretical results and simulations is again excellent; it clearly justifies the use of the Maxwellian distribution (4.12) to achieve accurate estimates of the cooling rates $\zeta _0$ and $\zeta _{0g}$ .

A point to consider here is that convergence to the results obtained using the Fokker–Planck model is achieved not only in the limit $m/m_g\to \infty$ and $m_0/m_g\to \infty$ but it is also necessary that $\sigma /\sigma _g\to \infty$ and $\sigma _0/\sigma _g\to \infty$ . For convenience, we assume in the rest of the work that $m/m_g\to \infty$ and $m_0/m_g\to \infty$ also imply $\sigma /\sigma _g\to \infty$ and $\sigma _0/\sigma _g\to \infty$ , and thus intruders and molecular gas particles have the same particle mass density (i.e. $m_0/\sigma _0^d=m_g/\sigma _g^d$ ).

5. Chapman–Enskog method. Diffusion transport coefficients

We assume that we perturb the homogeneous state by small spatial gradients. These perturbations induce non-vanishing contributions to the mass, momentum and heat fluxes. The determination of these fluxes allows us to identify the corresponding Navier–Stokes–Fourier transport coefficients of the granular suspension. As said in § 1, since in the tracer limit the pressure tensor and the heat flux vector of the binary mixture (intruders plus granular gas) are the same as that for the excess species (granular gas), the mass transport of intruders $\boldsymbol{j}_0$ is the relevant flux of the problem. The Navier–Stokes–Fourier transport coefficients of the granular gas were already determined by Gómez González & Garzó (Reference Gómez González and Garzó2022b ).

To get the mass flux $\boldsymbol{j}_0$ , the Boltzmann kinetic equation (3.1) is solved up to first order in spatial gradients by means of the Chapman–Enskog expansion Chapman & Cowling (Reference Chapman and Cowling1970) conveniently adapted to dissipative dynamics. As widely discussed in many textbooks (Chapman & Cowling Reference Chapman and Cowling1970; Ferziger & Kaper Reference Ferziger and Kaper1972), there are two different stages in the relaxation of a molecular gas towards equilibrium. For times of the order of the mean free time, one can identify a first stage (kinetic stage) where the main effect of collisions on the distribution function is to relax it towards the so-called local equilibrium state. Then, a hydrodynamic (slow) stage is achieved where the gas has completely forgotten its initial preparation. In this stage, the microscopic state of the gas is completely specified by the knowledge of the hydrodynamic fields (in the case of a binary mixture by $n_0$ , $n$ , $\boldsymbol{U}$ and $T$ ). These two stages are also expected in the case of granular gases except that in the kinetic stage the distribution function will generally relax towards a time-dependent non-equilibrium distribution (the homogeneous cooling state for freely cooling dry granular gases) instead of the local equilibrium distribution. A crucial point is that although the granular temperature $T$ is not a conserved field (due to the inelastic character of the collisions), it is assumed that $T$ can still be considered as a slow field. This assumption has been clearly supported by the good agreement found between granular hydrodynamics and computer simulations in several non-equilibrium situations (Lutsko et al. Reference Lutsko, Brey and Dufty2002; Dahl et al. Reference Dahl, Hrenya, Garzó and Dufty2002; Lois et al. Reference Lois, Lemaître and Carlson2007; Mitrano et al. Reference Mitrano, Dhal, Cromer, Pacella and Hrenya2011; Chialvo & Sundaresan Reference Chialvo and Sundaresan2013; Mitrano et al. Reference Mitrano, Garzó and Hrenya2014). More details on the application of the Chapman–Enskog method to dry (no gas phase) granular mixtures can be found, for example, in Garzó (Reference Garzó2019).

Based on the above arguments, in the hydrodynamic regime, the kinetic equation (3.1) admits a normal (or hydrodynamic) solution where all the space and time dependence of $f_0$ only occurs through a functional dependence on the hydrodynamic fields. As usual (Chapman & Cowling Reference Chapman and Cowling1970), this functional dependence can be made explicit by assuming small spatial gradients. In this case, $f_0$ can be written as a series expansion in powers of the spatial gradients of the hydrodynamic fields:

(5.1) \begin{equation} f_0=f_0^{(0)}+f_0^{(1)}+\cdots , \end{equation}

where the approximation $f_0^{(k)}$ is of order $k$ in the spatial gradients. In addition, in the presence of the gravitational force, it is necessary to characterise the magnitude of the force relative to that of the spatial gradients. Here, as in the case of conventional fluid mixtures (Chapman & Cowling Reference Chapman and Cowling1970), we assume that the magnitude of $\boldsymbol{g}$ is at least of first order in the perturbation expansion. The implementation of the Chapman–Enskog method to first order in the spatial gradients follows steps similar to those made in the conventional inelastic hard sphere model for dry granular mixtures (Garzó & Dufty Reference Garzó and Dufty2002; Garzó, Dufty & Hrenya et al. Reference Garzó, Dufty and Hrenya2007a ,Reference Garzó, Hrenya and Dufty b ).

In contrast to the application of the Chapman–Enskog method for dry granular mixtures, although we are interested here in calculating the diffusion transport coefficients in steady states, the presence of the surrounding molecular gas yields in inhomogeneous states a local energy unbalance between the energy supplied by the bath (or thermostat) and the energy lost by inelastic collisions. This means that we must first consider a time-dependent reference distribution $f^{(0)}(\boldsymbol{r}, \boldsymbol{v}, t)$ to obtain the time-dependent linear integral equations that verify the diffusion coefficients. One then assumes stationary conditions and solves (approximately) the above integral equations by considering the so-called first Sonine approximation. In addition, as discussed by Gómez González & Garzó (Reference Gómez González and Garzó2022b ), the term $\Delta \boldsymbol{U}=\boldsymbol{U}-\boldsymbol{U}_{\!g}$ must be considered to be at least of first order in spatial gradients. In this case, the Maxwellian distribution $f_g(\boldsymbol{r},\boldsymbol{v}, t)$ can be written as

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

where

(5.3) \begin{equation} f_g^{(0)}(\boldsymbol{V})=n_g \Big (\frac {m_g}{2\pi T_g}\Big )^{d/2} \exp \Bigg (-\frac {m_g V^2}{2T_g}\Bigg ), \quad f_g^{(1)}(\boldsymbol{V})=-\frac {m_g}{T_g}\boldsymbol{V}\boldsymbol{\cdot }\Delta \boldsymbol{U}f_g^{(0)}(\boldsymbol{V}). \end{equation}

The mathematical steps involved in the determination of the zeroth- and first-order distribution functions are quite similar to those made in previous works on granular mixtures (Garzó & Dufty Reference Garzó and Dufty2002; Garzó & Montanero Reference Garzó and Montanero2007). Technical details carried out in this derivation are provided in the supplementary material. In particular, the first-order distribution function $f^{(1)}(\boldsymbol{r}, \boldsymbol{v}, t)$ is given by

(5.4) \begin{eqnarray} f_0^{(1)}(\boldsymbol{V})&=&\boldsymbol{\mathcal{A}}_0(\boldsymbol{V})\boldsymbol{\cdot }\boldsymbol{\nabla }T+\boldsymbol{\mathcal{B}}_0(\boldsymbol{V})\boldsymbol{\cdot }\boldsymbol{\nabla }n+\boldsymbol{\mathcal{C}}_0(\boldsymbol{V})\boldsymbol{\cdot }\boldsymbol{\nabla }n_0+\mathcal{D}_{0}^{\prime}(\boldsymbol{V}) \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}\nonumber \\ & &+\mathcal{D}_{0,\textit{ij}}(\boldsymbol{V})\left (\partial _iU_{\!j}+\partial _{\!j} U_i-\frac {2}{d}\delta _{\textit{ij}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}\right ) +\boldsymbol{\varepsilon }_0(\boldsymbol{V})\boldsymbol{\cdot }\Delta \boldsymbol{U}, \end{eqnarray}

where the unknowns $(\boldsymbol{\mathcal{A}}_0, \boldsymbol{\mathcal{B}}_0, \boldsymbol{\mathcal{C}}_0, \mathcal{D}_{0,\textit{ij}},\mathcal{D}_{0}^{\prime},\boldsymbol{\varepsilon }_0)$ are the solutions of a set of coupled linear integral equations displayed in the supplementary material.

5.1. Diffusion transport coefficients

The constitutive equation for the mass flux $\boldsymbol{j}_0^{(1)}$ is given by (1.1). The diffusion transport coefficients are defined as

(5.5) \begin{align} D_T&=-\frac {m_0 T}{d\rho }\int{\rm d}\boldsymbol{v}\; \boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\mathcal{A}}_0(\boldsymbol{V}), \end{align}

(5.6) \begin{align} D&=-\frac {n}{d}\int {\rm d}\boldsymbol{v}\; \boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\mathcal{B}}_0(\boldsymbol{V}), \end{align}

(5.7) \begin{align} D_0&=-\frac {\rho }{dm_0}\int {\rm d}\boldsymbol{v}\; \boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\mathcal{C}}_0(\boldsymbol{V}), \end{align}

(5.8) \begin{align} D_0^U&=-\frac {m_0}{d}\int {\rm d}\boldsymbol{v}\; \boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\mathcal{\varepsilon }}_0(\boldsymbol{V}). \end{align}

The procedure for obtaining the expressions of the set of coefficients $(D_T, D, D_0, D_0^U)$ is described in the the supplementary material and only their final forms are provided in § 6. It is important to recall that the expressions of the diffusion coefficients cannot be analytically obtained for general unsteady conditions since it would require numerically solving a set of coupled differential equations for them. Thus, to reach analytical expressions for these diffusion coefficients, one assumes the validity of the steady-state constraint $\zeta +\zeta _g=0$ at each point of the system. This allows us to achieve explicit forms for the set $(D_T, D, D_0, D_0^U)$ .

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

As mentioned before, the diffusion transport coefficients are given in terms of the solutions of a set of coupled linear integral equations. As usual in kinetic theory of both molecular and granular gases, these integral equations can be approximately solved by considering the leading terms of a Sonine polynomial expansion of the unknowns $\boldsymbol{\mathcal{A}}_0$ , $\boldsymbol{\mathcal{B}}_0$ , $\boldsymbol{\mathcal{C}}_0$ and $\boldsymbol{\varepsilon }_0$ . In the case of the mass flux, these quantities are approximated by the polynomials

(6.1) \begin{align}& \boldsymbol{\mathcal{A}}_0(\boldsymbol{V})\to -f_{0,{M}}\boldsymbol{V}\frac {\rho }{Tn_0T_0}D_T, \quad \boldsymbol{\mathcal{B}}_0(\boldsymbol{V})\to -f_{0,{M}}\boldsymbol{V}\frac {m_0}{n n_0T_0}D, \end{align}

(6.2) \begin{align}&\quad \boldsymbol{\mathcal{C}}_0(\boldsymbol{V})\to -f_{0,{M}}\boldsymbol{V}\frac {m_0^2}{\rho n_0T_0}D_0, \quad \boldsymbol{\mathcal{\varepsilon }}_0(\boldsymbol{V})\to -f_{0,{M}}\boldsymbol{V}\frac {D_0^U}{n_0T_0}. \end{align}

To determine $D_0$ , $D$ , $D_T$ and $D_0^U$ , we substitute first $\boldsymbol{\mathcal{A}}_0$ , $\boldsymbol{\mathcal{B}}_0$ , $\boldsymbol{\mathcal{C}}_0$ and $\boldsymbol{\varepsilon }_0$ by their leading Sonine approximations (6.1) and (6.2) in the corresponding integral equations. Then we multiply these equations by $m_0\boldsymbol{V}$ and integrate over velocity. Technical details on these calculations are provided in the supplementary material.

6.1. Thermal diffusion coefficient $D_T$

The thermal diffusion coefficient $D_T$ is given by

(6.3) \begin{equation} D_T=\frac {n_0 T}{\rho \nu }D_T^*, \quad D_T^*=\frac {\tau _0-\frac {m_0}{m}+\chi \frac {\partial \tau _0}{\partial \chi }}{\beta \gamma ^*+\nu _D^*+\widetilde {\nu }_D\gamma _0^*}, \end{equation}

where $\gamma _0^*$ is defined by (4.18), $\tau _0=T_0/T$ and

(6.4) \begin{equation} \beta =\left (x^{-1}-3x\right )\mu ^{3/2}\chi ^{-1/2}. \end{equation}

Here, $x$ is given by (4.6) and in (6.3) we have introduced the (reduced) collision frequencies

(6.5) \begin{align}& \nu _D^*=\frac {2\pi ^{(d-1)/2}}{{d}\varGamma \left (\frac {d}{2}\right )}\left (\frac {\sigma '}{\sigma }\right )^{d-1}\mu ' \left (\frac {1+\theta _0}{\theta _0}\right )^{1/2}(1+\alpha _0), \end{align}

(6.6) \begin{align}&\qquad\quad \widetilde {\nu }_D=\left (\frac {m_0T_0}{m_g T_g}\right )^{1/2}\mu _{g0}\left (1+\theta _{0g}\right )^{1/2}, \end{align}

where

(6.7) \begin{equation} \theta _0=\frac {m_0T}{m T_0}, \quad \theta _{0g}=\frac {m_0T_g}{m_g T_0}. \end{equation}

While the quantity $\theta _0$ is the ratio of the mean square velocities of the intruders and granular gas particles, the quantity $\theta _{0g}$ gives the ratio of the mean square velocities of the intruders and molecular gas particles. Moreover, the explicit form of the derivative $\partial \tau _0/\partial \chi$ appearing in (6.3) is given in the supplementary material.

In the Brownian limiting case ( $m\gg m_g$ and $m_0\gg m_g$ ), $\widetilde {\nu }_D\to 1$ , $x\to \chi ^{-1/2}$ , $\beta \to 1-3\chi ^{-1}$ and (6.3) yields

(6.8) \begin{equation} D_T^*\to D_T^{*\text{B}}=\frac {\tau _0-\dfrac {m_0}{m}+\chi \dfrac {\partial \tau _0}{\partial \chi }}{\nu _D^*+\gamma _0^*-2\gamma ^*\chi ^{-1}-\dfrac {1}{2}\zeta ^*}. \end{equation}

Upon obtaining (6.8) use has been made of the steady-state condition

(6.9) \begin{equation} \chi ^{-1}\gamma ^*=\gamma ^*+\frac {1}{2}\zeta ^*. \end{equation}

The expression (6.8) agrees with previous results derived from the Langevin-like model based on the Fokker–Planck operators (2.17) and (3.12) (Gómez González & Garzó Reference Gómez González and Garzó2022a , Reference Gómez González and Garzó2023).

6.2. Mutual diffusion coefficient $D$

The mutual diffusion coefficient $D$ is

(6.10) \begin{equation} D=\frac {n_0T}{m_0\nu }D^*, \quad D^*=\frac {\zeta ^* D_T^*-\dfrac {m_0}{m}+\phi \dfrac {\partial \tau _0}{\partial \phi }}{\nu _D^*+\widetilde {\nu }_D\gamma _0^*}. \end{equation}

In (6.3) and (6.10), the derivatives $\partial \tau _0/\partial \chi$ and $\partial \tau _0/\partial \phi$ can be seen as a measure of the departure of the perturbed time-dependent state from the steady reference state. Their explicit forms and derivations are provided in the supplementary material.

In the Brownian limit ( $m\gg m_g$ and $m_0\gg m_g$ ), $\widetilde {\nu }_D\to 1$ and (6.10) reduces to

(6.11) \begin{equation} D^*\to D^{*{B}}=\frac {\zeta ^* D_T^*-\dfrac {m_0}{m}+\phi \dfrac {\partial \tau _0}{\partial \phi }}{\nu _D^*+\gamma _0^*}. \end{equation}

Equation (6.11) is consistent with the results obtained from the Langevin-like model (Gómez González & Garzó Reference Gómez González and Garzó2022a , Reference Gómez González and Garzó2023).

6.3. Tracer diffusion coefficient $D_0$

The tracer diffusion coefficient $D_0$ is given by

(6.12) \begin{equation} D_0=\frac {m n T}{m_0^2 \nu }D_0^*, \quad D_0^*=\frac {\tau _0}{\nu _D^*+\widetilde {\nu }_D\gamma _0^*}. \end{equation}

In the Brownian limit,

(6.13) \begin{equation} D_0^*\to D_0^{{B}*}=\frac {\tau _0}{\nu _D^*+\gamma _0^*}, \end{equation}

which agrees with previous results (Gómez González & Garzó Reference Gómez González and Garzó2022a , Reference Gómez González and Garzó2023).

6.4. Velocity diffusion coefficient $D_0^U$

The diffusion coefficient $D_0^{U}$ is given by

(6.14) \begin{equation} D_0^{U}=m_0 n_0 D_0^{U*}, \quad D_0^{U*}=\frac {\xi _0^*-\xi ^*}{\nu _D^*+\widetilde {\nu }_D\gamma _0^*}. \end{equation}

Here,

(6.15) \begin{align}& \xi _0^*=\frac {\xi _0}{\rho _0\nu }=\mu _{0g}\theta _{0g}^{-1/2}\left (1+\theta _{0g}\right )^{1/2}\gamma _0^*, \end{align}

(6.16) \begin{align}&\,\,\, \xi ^*=\frac {\xi }{\rho \nu }=\mu \theta ^{-1/2}\left (1+\theta \right )^{1/2}\gamma ^*, \end{align}

where $\rho _0=m_0n_0$ and $\theta =m T_g/(m_g T)$ .

In the Brownian limit, $\xi _0^*\to \gamma _0^*$ , $\xi ^*\to \gamma ^*$ and $\widetilde {\nu }_D\to 1$ , and (6.14) reduces to

(6.17) \begin{equation} D_0^{U{B}*}=\frac {\gamma _0^*-\gamma ^*}{\nu _D^*+\gamma _0^*}, \end{equation}

which is consistent with the previous results (Gómez González & Garzó Reference Gómez González and Garzó2022a , Reference Gómez González and Garzó2023) derived from the Langevin-like model.

6.5. Self-diffusion limiting case

Another interesting limiting case corresponds to the self-diffusion problem, namely when the intruders move in a granular gas whose particles are mechanically equivalent to it ( $\sigma =\sigma _0$ , $m=m_0$ , $\alpha =\alpha _0$ ). In this case, $\xi ^*=\xi _0^*$ , $\tau _0=1$ , $\partial _\chi \tau _0=\partial _\phi \tau _0=0$ , and so (6.3) and (6.14) yield $D_T=D_0^U=0$ as expected. Moreover, $D_0^*=-D^*$ and hence the constitutive equation (1.1) for the mass flux becomes

(6.18) \begin{equation} \boldsymbol{j}_0^{(1)}=-\frac {n T}{\nu }D_{0,{\textit{self}}}^* \boldsymbol{\nabla }x_0, \end{equation}

where $x_0=n_0/n$ is the concentration (or mole fraction) of the tracer species and the self-diffusion coefficient $D_{0,{\textit{self}}}^*$ is

(6.19) \begin{equation} D_{0,{\textit{self}}}^*=\frac {1}{\nu _{D,{\textit{self}}}^*+\gamma ^*\widetilde {\nu }_{D,{\textit{self}}}}. \end{equation}

In (6.16),

(6.20) \begin{equation} \nu _{D,{\textit{self}}}^*=\frac {\sqrt {2}\pi ^{(d-1)/2}}{d\varGamma \left (\dfrac {d}{2}\right )}(1+\alpha ),\quad \widetilde {\nu }_{D,{\textit{self}}}=\frac {m}{m+m_g}\left (1+\frac {m_g T}{m T_g}\right )^{1/2}. \end{equation}

Figure 3. Plot of the (scaled) thermal diffusion coefficient $D_T(\alpha )/D_T(1)$ versus the (common) coefficient of restitution $\alpha =\alpha _0$ for $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). In all the curves $m_0/m=8$ , $\sigma _0/\sigma =2$ and $\sigma _0/\sigma _g=(m_0/m_g)^{1/3}$ . The dashed line refers to the expression (6.8) derived in the Brownian limiting case for the ratio $D_T(\alpha )/D_T(1)$ .

Figure 4. Plot of the (scaled) mutual diffusion coefficient $D(\alpha )/D(1)$ versus the (common) coefficient of restitution $\alpha =\alpha _0$ for $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). In all the curves $m_0/m=8$ , $\sigma _0/\sigma =2$ and $\sigma _0/\sigma _g=(m_0/m_g)^{1/3}$ . The dashed line refers to the expression (6.11) derived in the Brownian limiting case for the ratio $D(\alpha )/D(1)$ .

Figure 5. Plot of the (scaled) tracer diffusion coefficient $D_0(\alpha )/D_0(1)$ versus the (common) coefficient of restitution $\alpha =\alpha _0$ for $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). In all the curves $m_0/m=8$ , $\sigma _0/\sigma =2$ and $\sigma _0/\sigma _g=(m_0/m_g)^{1/3}$ . The dashed line refers to the expression (6.13) derived in the Brownian limiting case for the ratio $D_0(\alpha )/D_0(1)$ .

Figure 6. Plot of the (scaled) velocity diffusion coefficient $D_0^{U}(\alpha )/D_0^{U}(1)$ versus the (common) coefficient of restitution $\alpha =\alpha _0$ for $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). In all the curves $m_0/m=8$ , $\sigma _0/\sigma =2$ and $\sigma _0/\sigma _g=(m_0/m_g)^{1/3}$ . The dashed line refers to the expression (6.17) derived in the Brownian limiting case for the ratio $D_0^{U}(\alpha )/D_0^{U}(1)$ .

7. An application: segregation of intruders in a granular suspension

As mentioned in § 1, an interesting application of the results displayed in § 6 is the the study of the segregation of intruders by thermal diffusion in a granular suspension. Needless to say, thermal diffusion segregation is likely one of the most common phenomena appearing in polydisperse systems. It occurs in a non-convective steady state ( $\boldsymbol{U}=\boldsymbol{U}_{\!g}=\boldsymbol{0}$ ) due to the existence of a temperature gradient, which causes the movement of the different species of the mixture. In the steady state, remixing of species generated by diffusion is balanced by segregation caused by temperature differences. To quantify the degree of segregation along the temperature gradient it is usual to introduce the thermal diffusion factor $\varLambda$ (Kincaid et al. Reference Kincaid, Cohen and López de Haro1987). In a steady state without convection and where the mass flux is zero ( $\boldsymbol{j}_0^{(1)}=\boldsymbol{0}$ ), the thermal diffusion factor is defined as

(7.1) \begin{equation} -\varLambda \frac {\partial \ln T}{\partial z}=\frac {\partial }{\partial z}\ln \Big (\frac {n_0}{n}\Big ), \end{equation}

where for the sake of simplicity we have assumed that gradients occur only along the $z$ axis. In addition, we also assume that the gravitational field is parallel to the thermal gradient, namely $\boldsymbol{g}=-g\widehat {e}_z$ , where $\widehat {e}_z$ is the unit vector in the positive direction of the $z$ axis.

We consider a scenario in which the intruders have a larger size than the granular gas particles ( $\sigma _0\gt \sigma$ ). Furthermore, as said before, since gravity and the thermal gradient point in the same direction, then the lower plate is hotter than the upper plate ( $\partial _z\ln T\lt 0$ ). Based on (7.1), when $\varLambda \gt 0$ , intruders rise relative to granular gas particles ( $\partial _z \ln (n_0/n)\gt 0$ ), leading to an accumulation of tracer particles near the cooler plate. This situation is commonly known as the Brazil nut effect (BNE). Conversely, for $\varLambda \lt 0$ , intruders fall relative to granular gas particles ( $\partial _z \ln (n_0/n)\lt 0$ ), resulting in an accumulation near the hotter plate. This situation is known as the reverse Brazil nut effect (RBNE). A representative diagram of segregation dynamics is shown in figure 7.

Figure 7. A representative diagram of the BNE ( $\varLambda \gt 0$ ) and RBNE ( $\varLambda \lt 0$ ) for a ternary system composed of molecular particles (blue), granular particles (green) and intruders (red).

We express the thermal diffusion factor in terms of the diffusion transport coefficients. In the steady state, to first order in spatial gradients, the momentum balance (2.10) reduces to

(7.2) \begin{equation} \left (1-\varepsilon \kappa ^*\right )+\left (1-\varepsilon \overline {\mu }^*\right )\frac {\partial _z \ln n}{\partial _z \ln T}=-g^*, \end{equation}

where

(7.3) \begin{equation} g^*=\frac {\rho g}{n\partial _z T}\lt 0 \end{equation}

is a dimensionless parameter measuring the gravity relative to the thermal gradient, $\kappa ^*$ and $\overline {\mu }^*$ are the (reduced) thermal conductivity and diffusive heat conductivity transport coefficients, respectively, and

(7.4) \begin{equation} \varepsilon =\frac {(d+2)\sqrt {\pi }}{2^{d+3}(d-1)}\frac {\mu \chi ^{-1/2}}{\phi \sqrt {T_g^*}}X(\theta ). \end{equation}

Here, $X(\theta )=\theta ^{-1/2}(1+\theta )^{-1/2}$ . The explicit forms of $\kappa ^*$ and $\overline {\mu }^*$ are displayed in the supplementary material. Upon obtaining (7.2), use has been made of the result in the leading Sonine approximation (Gómez González & Garzó Reference Gómez González and Garzó2022b ):

(7.5) \begin{equation} \mathcal{F}^{(1)}\big[f^{(1)}\big]\to \frac {1}{d+2}\frac {\rho \mu \gamma }{n}\kappa _0 X(\theta ) \big (\kappa ^* \partial _z \ln T+\overline {\mu }^* \partial _z \ln n\big ), \end{equation}

where

(7.6) \begin{equation} \kappa _0=\frac {d(d+2)^2\varGamma \left (\dfrac {d}{2}\right )}{16(d-1)\pi ^{\frac {d-1}{2}}}\sigma ^{1-d}\sqrt {\frac {T}{m}} \end{equation}

is the low-density value of the thermal conductivity for an ordinary gas of hard spheres.

According to (1.1), when $\Delta \boldsymbol{U}=\boldsymbol{0}$ , the condition $j_{0,z}^{(1)}=0$ yields the relation

(7.7) \begin{equation} -D_0^* \partial _z \ln n_0=D^* \partial _z \ln n+D_T^* \partial _z \ln T, \end{equation}

where $D_T^*$ , $D^*$ and $D_0^*$ are given by (6.3), (6.10) and (6.12), respectively. From (7.2) and (7.7), the thermal diffusion factor can be written as

(7.8) \begin{equation} \varLambda =\frac {\partial _z \ln n}{\partial _z \ln T}-\frac {\partial _z \ln n_0}{\partial _z \ln T} =\frac {D_T^*+\left (D_0^*+D^*\right )\left (\varepsilon \kappa ^*-1-g^*\right )\left (1-\varepsilon \overline {\mu }^*\right )^{-1}}{D_0^*}. \end{equation}

In the Brownian limit ( $m/m_g\to \infty$ ), $\mu \to 1$ and $\theta \to \infty$ , so that $\varepsilon \to 0$ . In this limiting case, $\varLambda$ reduces to

(7.9) \begin{equation} \varLambda =\frac {D_T^*-\left (D_0^*+D^*\right )\left (1+g^*\right ) }{D_0^*}. \end{equation}

Equation (7.9) agrees with previous results derived for the segregation of massive intruders in a granular suspension (Gómez González & Garzó Reference Gómez González and Garzó2022a ).

Since (6.12) clearly shows that $D_0^*\gt 0$ , then the curves delineating the regimes between the segregation towards the cold and the hot wall (BNE/RBNE transition) are determined from the condition

(7.10) \begin{equation} \left (1-\varepsilon \overline {\mu }^*\right )D_T^*=-\left (D_0^*+D^*\right )\left (\varepsilon \kappa ^*-1-g^*\right ). \end{equation}

7.2. Different mechanical properties

7.2.1. Absence of gravity ( $|{g}^*|\to 0$ )

Let us first consider a scenario where gravitational effects are negligible ( $|{g}^*|\to 0$ ). Under this assumption, the condition $\varLambda =0$ reads

(7.11) \begin{equation} \left (1-\varepsilon \overline {\mu }^*\right )D_T^*=-\left (D_0^*+D^*\right )\left (\varepsilon \kappa ^*-1\right ). \end{equation}

Figure 8. Plot of the marginal segregation curve ( $\varLambda =0$ ) for a (common) coefficient of restitution $\alpha =\alpha _0=0.7$ and $|{g}^*|\to 0$ . The parameters used are: $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ (20, 50, 100 and 1000). The points below the curve correspond to $\varLambda \gt 0$ (BNE), while the points above the curve correspond to $\varLambda \lt 0$ (RBNE). The dashed line is the result obtained in the Brownian limiting case.

Figure 8 shows the BNE/RBNE phase diagram for a three-dimensional system ( $d=3$ ) with a common coefficient of restitution $\alpha =\alpha _0=0.7$ , $T_g^*=10$ and five distinct values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100,\ 1000$ , and the Brownian limit $m_0/m_g\to \infty$ ). At first glance, we notice a quantitative discrepancy with the diagrams shown in Gómez González & Garzó (Reference Gómez González and Garzó2023) for $|{g}^*|\to 0$ . These discrepancies were expected since, in that work, segregation was calculated for moderate densities, unlike here ( $\phi =0.0052$ ). Also, when using the Fokker–Planck approach (3.12) for the operator $J_{0g}[f_0,f_g]$ , we see that the expression (4.18) of $\gamma _0^*$ employed here slightly differs from the expression of $\gamma ^*_{0,{FP}}$ defined by Gómez González & Garzó (Reference Gómez González and Garzó2023). Specifically, $\gamma ^*_{0,{FP}}=(\sigma _0/\sigma )(\overline {\sigma }/\overline {\sigma }_0)^{d-1}\gamma _0^*$ . This discrepancy is because $\gamma ^*_{0,\text{FP}}$ was obtained from the granular literature using phenomenological arguments.

However, what is surprising is the complete reversal we see in the diagram when we observe the RBNE effect as we increase $m_0/m$ , contrary to what is observed in Gómez González & Garzó (Reference Gómez González and Garzó2023). This requires a more subtle explanation. When calculating the diagrams in the case of an effective model where the thermostat is modelled by a Fokker–Planck equation, the thermostat intensity is regulated only by the (dimensionless) bath temperature $T_{g}^*$ . If this parameter is kept constant, as in the figures 11 and 12 reported by Gómez González & Garzó (Reference Gómez González and Garzó2023), the thermostat does not change when the mass ( $m_0/m$ ) or size ( $\sigma _0/\sigma$ ) ratios are modified. However, in our case, we modify $m_0/m$ (or $\sigma _0/\sigma$ ) keeping $m_0/m_g$ constant for each particular curve. Thus, by modifying the mass ratio $m_0/m$ for a particular value of $m_0/m_g$ , the relative mass between the granular and molecular gas changes. Therefore, the effect of the collisions between the grains and the particles of the molecular gas will be different. Concretely, if we increase $m_0/m$ without changing $m_0/m_g$ , the molecular gas will have a mass increasingly similar to that of the grains. Consequently, the thermalising effect of the molecular gas that compensates for the effect of inelasticity will be more effective, causing the temperature of the granular gas to be higher and thus tend to increase with respect to that of the intruders (RBNE). We observe the same when we increase the size ratio ( $\sigma _0/\sigma$ ). In this case, for a given $m_0/m_g$ (or equivalently $\sigma _0/\sigma _g$ ), as the size of the intruders increases, the grains will have a size increasingly similar to the particles of the surrounding molecular gas; thus the effective area of the grain in a collision decreases, and, with it, the number of collisions. Therefore, the grain will have less effective thermalisation and will move to the cold zone (BNE).

On the other hand, as we increase $m_0/m_g$ , we see that the transition to the RBNE phase occurs earlier. This is because, by increasing the mass more rapidly than the size ( $m_0/m_g=(\sigma _0/\sigma _g)^d$ ), the effect of intruder–molecular gas collisions in the motion of intruders becomes less effective. In the curves where the mass ratios $m_0/m_g$ are close to each other, the relative size also plays an important role, and crossings can be observed, as was the case in figure 4.

Moreover, for elastic collisions ( $\alpha = \alpha _0 = 1$ ), the segregation criterion notably deviates from the classical result obtained for molecular mixtures of hard spheres by Kincaid et al. (Reference Kincaid, Cohen and López de Haro1987), where, in the first Sonine approximation, the condition $\varLambda =0$ yields the simple segregation criterion $m_0/m = 1$ . For the present system (intruders moving in a granular gas immersed into a molecular gas), our analysis reveals that there is no segregation for the remaining parameters considered in figure 8.

7.2.2. Thermalised systems ( $\partial _z T\to 0$ )

Let us explore a scenario where gravity is the main factor influencing segregation dynamics. In this situation, $|{g}^*|\to \infty$ , which makes the temperature gradient negligible ( $\partial _y T\to 0$ ), and the condition $\varLambda =0$ leads to the relation

(7.12) \begin{equation} D_{0}^*+D^*=0. \end{equation}

Figure 9. Plot of the marginal segregation curve ( $\varLambda =0$ ) for a (common) coefficient of restitution $\alpha =\alpha _0=0.7$ and and $|{g}^*|\to \infty$ . The parameters used are: $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). The points below the curve correspond to $\varLambda \gt 0$ (BNE), while the points above the curve correspond to $\varLambda \lt 0$ (RBNE). The dashed line is the result obtained by using the Fokker–Planck approach (2.17) to the operator $J_g[f,f_g]$ .

Figure 9 shows the case $|{g}^*|\to \infty$ with the same parameters as figure 8. In this strong-gravity case, the explanation is simpler. Gravity is much stronger than both the energy coming from the molecular gas and the energy lost in inelastic collisions, as seen in granular suspensions and driven granular gases (Garzó Reference Garzó2008; Gómez González & Garzó Reference Gómez González and Garzó2023; Gómez González et al. Reference Gómez González, Garzó, Brito and Soto2024). Because of this, all particles fall to the bottom plate. Heavier intruders are harder to lift, so they move downward (RBNE). However, if the intruders are larger but keep the same mass, smaller particles hit more often. These collisions push the smaller particles below the intruders, lifting them and creating a buoyancy effect (BNE)

For elastic collisions ( $\alpha =\alpha _0=1$ ), the cooling rate vanishes ( $\zeta ^*=0$ ), and there is equipartition of energy ( $\tau _0=1$ ). Thus, it is straightforward to verify from (6.10) and (6.12) that the segregation criterion simplifies to $m_0/m=1$ . This result aligns with previous findings for dry granular gases (i.e. in the absence of an interstitial gas phase) (Brey et al. Reference Brey, Ruiz-Montero and Moreno2005; Garzó Reference Garzó2006, Reference Garzó2011), as well as with those obtained in the coarse-grained approach by using the Fokker–Planck equation (Gómez González & Garzó Reference Gómez González and Garzó2023).

7.2.3. General case

As a final situation, we consider the general case where the effects of the temperature gradient and gravity are comparable. To exemplify this, figure 10 shows the marginal segregation curve for a reduced gravity $|{g}^*|=1$ and for the same systems as depicted in figures 8 and 9.

Figure 10. Plot of the marginal segregation curve ( $\varLambda =0$ ) for a (common) coefficient of restitution $\alpha =\alpha _0=0.7$ and and $|{g}^*|=1$ . The parameters used are: $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). The points below the curve correspond to $\varLambda \gt 0$ (BNE), while the points above the curve correspond to $\varLambda \lt 0$ (RBNE). The dashed line is the result obtained by using the Fokker–Planck approach (2.17) to the operator $J_g[f,f_g]$ .

Figure 11. Plot of the marginal segregation curve ( $\varLambda =0$ ) for a (common) coefficient of restitution $\alpha =\alpha _0=1$ and and $|{g}^*|=1$ . The parameters used are: $d=3$ , $\phi =0.0052$ , $T_g^*=10$ and four different values of the mass ratio $m_0/m_g$ ( $20,\ 50,\ 100$ and 1000). The points below the curve correspond to $\varLambda \gt 0$ (BNE), while the points above the curve correspond to $\varLambda \lt 0$ (RBNE). The dashed line is the result obtained by using the Fokker–Planck approach (2.17) to the operator $J_g[f,f_g]$ .

The main point is that, unlike in dry granular mixtures and granular suspensions (Garzó Reference Garzó2019; Gómez González & Garzó Reference Gómez González and Garzó2023; Gómez González et al. Reference Gómez González, Garzó, Brito and Soto2024), gravity has a weaker effect than the thermal gradient, keeping the BNE/RBNE transition quite similar to that shown in figure 8. The explanation can be the same as in the case without gravity since, when increasing $m_0/m$ , the relative mass of the grains compared with molecular particles decreases. As a result, the thermalisation due to the surrounding fluid does not remain constant, making thermal effects more pronounced even when gravity is present. To complement the phase diagram shown in figure 10 for inelastic collisions, in figure 11 we plot the marginal segregation curve ( $\varLambda =0$ ) for the same values of the mass ratio $m_0/m_g$ as in figure 10 but for elastic collisions. It is quite apparent from the comparison between figures 10 and 11 that the inelasticity of collisions plays a secondary role in the segregation behaviour in this case since we observe practically no differences between both figures.

8. Summary and concluding remarks

This paper aims to determine the diffusion transport coefficients for tracer (or intruder) particles within a granular gas modelled as a gas of inelastic hard spheres. Intruders and granular gas are immersed in a bath of elastic hard spheres (molecular gas). We examine scenarios where the granular particles are sufficiently dilute, ensuring the molecular gas remains unaffected and serves as a thermostat at a given temperature $T_g$ . Unlike other suspension models, which consider an effective fluid–solid force, our model accounts for both inelastic collisions between the tracer and granular particles, as well as among the granular particles themselves. Additionally, it takes into account the elastic collisions between the grains and molecular gas particles, as well as between the intruders and molecular gas. We consider the low-density regime for the suspended solid particles and, hence, the velocity distribution functions $f(\boldsymbol{r}, \boldsymbol{v}, t)$ and $f_0(\boldsymbol{r}, \boldsymbol{v}, t)$ for grains and intruders, respectively, obey the (inelastic) Boltzmann equations.

To ensure a consistent theoretical framework, we first analyse homogeneous reference states as a basis for applying the Chapman–Enskog method and obtaining then the diffusion transport coefficients. In the homogeneous state, we present new results for the temperature ratio $\chi _0 = T_0/T_g$ and the kurtosis $c_0$ associated with the intruders. As in the case of the granular gas (Gómez González & Garzó Reference Gómez González and Garzó2022b ), we find that the tracer temperature ratio $\chi _0$ shows larger deviations from unity as the mass ratio $m_0/m_g$ decreases. Regarding the kurtosis $c_0$ (which measures the departure of the tracer distribution from its Maxwellian form), as for the granular gas (Gómez González & Garzó Reference Gómez González and Garzó2022b ), our results clearly show that this quantity remains small enough to validate the use of the Maxwellian approximation for the velocity distributions.

We compared the theoretical results derived in the HSS for $\chi _0$ and $c_0$ with DSMC data, and the agreement is remarkable. This justifies the use of the Maxwellian distribution (4.12) to achieve accurate estimates of the cooling rates $\zeta _0$ and $\zeta _{0g}$ . We also observe the convergence of the results in the Brownian limit ( $m/m_g \to \infty$ and $m_0/m_g \to \infty$ ) with those obtained in previous works (Gómez González & Garzó Reference Gómez González and Garzó2019, Reference Gómez González and Garzó2023) by using the Fokker–Planck approach.

Once the homogeneous state is characterised, the corresponding set of kinetic equations for the mixture were addressed using the Chapman–Enskog method (Chapman & Cowling Reference Chapman and Cowling1970), up to first order in spatial gradients. From this solution, the four diffusion transport coefficients are derived by considering the leading terms in a Sonine polynomial expansion of the first-order distribution function. In dimensionless form, these diffusion coefficients are given in terms of eight dimensionless parameters: the diameter ratios $\sigma /\sigma _g$ and $\sigma _0/\sigma _g$ , the mass ratios $m/m_g$ and $m_0/m_g$ , the coefficients of restitution $\alpha$ and $\alpha _0$ , the (reduced) bath temperature $T_g^*$ and the volume fraction $\phi$ (which is considered to be quite small since we are considering a very dilute granular gas). Compared with previous attempts to obtain tracer diffusion coefficients in granular suspensions (see e.g. Gómez González & Garzó Reference Gómez González and Garzó2023) by starting from a coarse-grained approach, our results provide a more general description of the dependence of diffusion on mass ( $m/m_g$ and $m_0/m_g$ ) and diameter ( $\sigma /\sigma _g$ and $\sigma _0/\sigma _g$ ) ratios. As expected, the expressions for the tracer diffusion coefficients agree with previous results derived in the Brownian limit from the Fokker–Planck operator (Gómez González & Garzó Reference Gómez González and Garzó2023). Thus, with respect to previous attempts for obtaining tracer diffusion coefficients in granular suspensions (Gómez González & Garzó Reference Gómez González and Garzó2023) by starting from a coarse-grained approach, our results provide a more general description of the dependence of the diffusion on the mass ( $m/m_g$ and $m_0/m_g$ ) and diameter ( $\sigma /\sigma _g$ and $\sigma _0/\sigma _g$ ) ratios. In particular, as expected, the expressions for the tracer diffusion coefficients agree with previous results derived in the Brownian limit from the Fokker–Planck operator (Gómez González & Garzó Reference Gómez González and Garzó2023). In this context, the present results encompass the previous ones as a special case, thus covering interesting physical situations that had not previously been analysed from a theoretical perspective.

In general, the diffusion coefficients exhibit significant deviations from their elastic counterparts, particularly under strong inelasticity. A key result is that convergence to the predictions of the effective model is achieved only at large mass ratios $m_0/m_g \approx 1000$ , which contrasts with the convergence observed in a monocomponent granular suspension where agreement was found at $m/m_g= 50$ (Gómez González & Garzó Reference Gómez González and Garzó2022a ). This emphasises the ability of the model to capture more realistic scenarios beyond the scope of Fokker–Planck-based models.

As an application, we investigate segregation driven by both thermal gradients and gravity. We find that increasing the mass ratio $m_0/m_g$ tends to push the intruders towards the bottom plate (RBNE), regardless of the gravitational strength. This behaviour contrasts with the results derived from the Fokker–Planck model (Gómez González & Garzó Reference Gómez González and Garzó2023). In that model, the thermostat effect remains constant (set by $T_g^*$ ), whereas in our model it depends on $m_0/m_g$ , resulting in different thermalisation dynamics. The main novelty of the study reported here, compared with previous works on dry granular mixtures (see e.g. Garzó Reference Garzó2008, Reference Garzó2011), is the analysis of the role of $m_0/m_g$ at fixed $\alpha$ . Unlike the dry case, where inelasticity strongly affects BNE/RBNE transitions due to energy non-equipartition, here its impact is reduced because the molecular gas injects energy homogeneously, compensating collisional dissipation and effectively suppressing partial temperature differences through thermalisation.

The results reported in this work suggest a framework to study experimentally the influence of an interstitial gas on segregation dynamics. In this context, a feasible experiment can be designed following the approach discussed in Gómez González et al. (Reference Gómez González, Abad, Bravo Yuste and Garzó2023). A suitable system to represent the low-Reynolds-number regime, in which the effect of granular collisions is comparable to the effect of the thermal bath, can be realised by immersing gold grains in hydrogen gas at low pressure. In this set-up, the molecular gas provides a homogeneous thermal bath with appropriate viscosity $\eta _g$ , such that the Reynolds numbers remain very small ( ${Re} \sim 10^{-4}$ ) and the Stokes number, $T_g^*$ , as well as $\gamma /\nu$ , are all close to unity. This ensures that the molecular gas effectively thermalises the granular particles, while the effect of granular collisions remains significant.

Once the particles composing the system are selected, an experimental set-up to investigate segregation of intruders in a granular gas immersed in an interstitial fluid can be designed based on previous studies (Möbius et al. Reference Möbius, Lauderdale, Nagel and Jaeger2001; Naylor, Swift & King Reference Naylor, Swift and King2003; Sánchez, Swift, & King Reference Sánchez, Swift and King2004; Wylie et al. Reference Wylie, Zhang, Xu and Sun2008; Clement et al. Reference Clement, Pacheco-Martínez, Swift and King2010; Pastenes, Géminard & Melo Reference Pastenes, Géminard and Melo2014). The set-up can consist of a transparent container with a porous base to allow fluid flow while maintaining particle collisions, filled with small granular particles and a single intruder whose mass and size can be adjusted. The interstitial fluid, either a gas or a liquid, occupies the voids between particles, enabling the study of drag and thermalisation effects. The container is subjected to controlled vertical or horizontal vibrations with adjustable amplitude, and the position of the intruder is tracked using high-speed imaging. By varying the intruder’s mass and size and the properties of the fluid, this set-up allows for systematic measurement of segregation phenomena, including the influence of mass ratio and interstitial fluid on the BNE/RBNE transition.

Finally, it is important to remark that to facilitate the solution of the integral equations verifying the mixture, we have considered the tracer limit where one of the species is present at a negligible concentration. A possible extension of this work is to generalise diffusion to a binary mixture with arbitrary concentrations. In addition, it may be interesting to study the role of density in the diffusion and segregation diagrams using the Enskog kinetic equation which applies to moderate densities. Moreover, revisiting the problem within a simplified random-walk framework could shed light on the underlying mechanisms responsible for the observed non-monotonicities and the crossing of curves in the tracer diffusion coefficient. These studies will be developed in future projects.