# Tracer diffusion coefficients in a moderately dense granular suspension.
Stability analysis and thermal diffusion segregation

Rubén Gómez González<EMAIL_ADDRESS>Departamento de Física, Universidad de
Extremadura, E-06006 Badajoz, Spain  Vicente Garzó<EMAIL_ADDRESS>https://fisteor.cms.unex.es/investigadores/vicente-garzo-puertos Departamento
de Física and Instituto de Computación Científica Avanzada (ICCAEx),
Universidad de Extremadura, E-06006 Badajoz, Spain

###### Abstract

The diffusion transport coefficients of a binary granular suspension where one
of the components is present in tracer concentration are determined from the
(inelastic) Enskog kinetic equation. The effect of the interstitial gas on the
solid particles is accounted for in the kinetic equation through two different
terms: (i) a viscous drag force proportional to the particle velocity and (ii)
stochastic Langevin-like term defined in terms of the background temperature.
The transport coefficients are obtained as the solutions of a set of coupled
linear integral equations recently derived for binary granular suspensions
with arbitrary concentration [Gómez González _et al._ , Phys. Rev. E 101,
012904 (2020)]. To achieve analytical expressions for the diffusion
coefficients which can be sufficiently accurate for highly inelastic
collisions and/or disparate values of the mass and diameter rations, the above
integral equations are approximately solved by considering the so-called
second Sonine approximation (two terms in the Sonine polynomial expansion of
the distribution function). The theoretical results for the tracer diffusion
coefficient $D_{0}$ (coefficient connecting the mass flux with the gradient of
density of tracer particles) are compared with those obtained by numerically
solving the Enskog equation by means of the direct simulation Monte Carlo
method. Although the first-Sonine approximation to $D_{0}$ yields in general a
good agreement with simulation results, we show that the second-Sonine
approximation leads to an improvement over the first-Sonine correction,
especially when the tracer particles are much lighter than the granular gas.
The expressions derived here for the diffusion coefficients are also used for
two different applications. First, the stability of the homogeneous steady
state is discussed. Second, segregation induced by a thermal gradient is
studied. As expected, the results show that the corresponding phase diagrams
for segregation clearly differ from those found in previous works when the
effect of gas phase on grains is neglected.

##  I Introduction

The understanding of transport processes occurring in a granular suspension
(granular particles immersed in a fluid, like air or water) is still a
challenging open problem. This is essentially due to the fact that the
corresponding granular flow is a multiphase process involving quite different
spatial and temporal scales. The problem is much more difficult when one
considers _multicomponent_ gas-solid flows (i.e ., when grains have different
masses and sizes), specially when one intends to describe the granular
suspension in terms of a set of coupled kinetic equations for each one of the
one-particle velocity distributions of the different phases.

Thus, because of the complexity embodied in the description of multicomponent
granular suspensions, a coarse-grained approach is usually adopted where the
influence of the background (interstitial) fluid on grains is accounted for
through an effective fluid-solid force. Koch (1990); Gidaspow (1994); Jackson
(2000); Koch and Hill (2001) In some works, Gidaspow (1994); Zhao and Wang
(2021) this force is simply proportional to the relative mean flow velocity
between the solid and gas phases. In another studies, Tsao and Koch (1995);
Sangani  _et al._ (1996); Wylie  _et al._ (2009); Parmentier and Simonin
(2012); Heussinger (2013); Wang  _et al._ (2014); Saha and Alam (2017); Alam,
Saha, and Gupta (2019); Saha and Alam (2020) each solid particle is only
subjected to the Stokes drag force which takes into account the friction of
grains on the viscous surrounding fluid. A more sophisticated model Garzó  _et
al._ (2012) accounts for the thermal fluctuations arising from the fluid and
incorporates a stochastic Langevin-like term defined in terms of the (known)
temperature of the interstitial gas. The stochastic term takes into account
the energy transferred to grains due to their “interactions” with the
background particles.

The Navier–Stokes transport coefficients of a binary granular suspension have
been recently obtained Gómez González, Khalil, and Garzó (2020) in the context
of the (inelastic) Enskog kinetic equation. As in Refs. Garzó  _et al._ ,
2012; Gómez González and Garzó, 2019, a coarse-grained description has been
also adopted and hence, the influence of the interstitial gas on solid
particles is via a viscous drag force plus a stochastic Langevin-like term.
The Enskog equation has been solved by means of the Chapman–Enskog method
Chapman and Cowling (1970) where, in contrast to other perturbation schemes,
Sela and Goldhirsch (1998); Serero  _et al._ (2006); Zhao and Wang (2021) the
Knudsen number (which measures the strength of the spatial variations of the
hydrodynamic fields on the scale of the mean free path) and not the degree of
inelasticity has been used as the only relevant perturbation parameter. As a
consequence, although the results reported in Ref. Gómez González, Khalil, and
Garzó, 2020 are limited to the Navier–Stokes order (first order in spatial
gradients), they go beyond the weak dissipation limit and hence can be in
principle applicable to a wide range of coefficients of restitution. In
addition, they take into account the nonequipartition of granular energy and
can be easily extended to any number of components. In this context, the
kinetic theory exposed in Ref. Gómez González, Khalil, and Garzó, 2020 extends
the theory developed by Garzó _et al._ Garzó, Dufty, and Hrenya (2007); Garzó,
Hrenya, and Dufty (2007) for _dry_ (no gas phase) granular mixtures to gas-
solid flows. It is important to note that the theoretical results for the
shear viscosity coefficient of a binary granular mixture obtained in Refs.
Garzó, Dufty, and Hrenya, 2007; Garzó, Hrenya, and Dufty, 2007 from the
Chapman–Enskog solution has been recently shown Garzó (2021); García Chamorro
and Garzó (2023) to compare much better with computer simulations than those
obtained by assuming Maxwellian distributions for each species. Solsvik and
Manger (2021)

As occurs for elastic Chapman and Cowling (1970) and inelastic Garzó, Dufty,
and Hrenya (2007); Garzó, Hrenya, and Dufty (2007) collisions, the explicit
forms of the kinetic contributions to the Navier–Stokes transport coefficients
displayed in Ref. Gómez González, Khalil, and Garzó, 2020 requires to solve a
set of coupled linear integral equations. The usual way of obtaining these
kinetic contributions consists of approximating the solutions to these
integral equations by Maxwellian distributions (defined at different partial
temperatures) corrected by truncated Sonine polynomial expansions. For the
sake of simplicity, the leading Sonine polynomial (the so-called first Sonine
approximation) is retained in most of the works. However, although the results
obtained from this approximation for _dry_ (no gas phase) granular mixtures
compare in general quite well with simulations for not too large values of the
coefficients of restitution, Garzó and Montanero (2004, 2007); García Chamorro
and Garzó (2023) they exhibit significant discrepancies with computer
simulations for strong inelasticity and/or disparate values of the mass and
diameter ratios. Garzó and Vega Reyes (2009, 2012) These differences can be
mitigated in part if one retains second-order terms in the Sonine polynomial
expansion (second Sonine approximation). Thus, it is important to assess the
degree of accuracy of at least the two first Sonine solutions for dense binary
granular suspensions via a comparison of the corresponding theoretical results
with appropriate computer simulations.

Needless to say, the determination of the Navier–Stokes transport coefficients
for gas-solid flows beyond the first Sonine approximation is a quite
cumbersome problem. These technical difficulties increase considerably when
one studies multicomponent systems due specially to the coupling among the
different integral equations obeying the transport coefficients. Thus, in
order to gain some insight, one considers simple situations where explicit
results can be in principle attainable. Here, we consider the so-called
_tracer_ limit, namely, a granular binary suspension where the concentration
of one of the species (of mass $m_{0}$ and diameter $\sigma_{0}$) is much
smaller than the other one (of mass $m$ and diameter $\sigma$). As usual, in
the tracer limit, one can assume that (i) the state of the excess component is
not perturbed by the presence of the tracer particles and (ii) one can also
neglect the collisions among tracer particles in its kinetic equation. In
these conditions, while the distribution function $f$ of the granular fluid
obeys a nonlinear (closed) Enskog equation, the distribution $f_{0}$ of the
tracer species verifies the (linear) Enskog–Lorentz equation. The influence of
the interstitial gas in both kinetic equations on the distributions $f$ and
$f_{0}$ is accounted for by means of Fokker–Planck terms with different
coefficients of friction $\gamma$ and $\gamma_{0}$, respectively.

In the tracer limit, the stress tensor and the heat flux of the binary
granular suspension (granular gas plus tracers) are the same as those for the
monodisperse granular suspension. The expressions of these fluxes were already
determined in Ref. Gómez González and Garzó, 2019. Consequently, the mass
transport of tracer particles $\mathbf{j}_{0}$ is the relevant flux of the
problem. The constitutive equation of $\mathbf{j}_{0}$ to Navier–Stokes order
is given byGómez González, Khalil, and Garzó (2020)

$\mathbf{j}_{0}^{(1)}=-\frac{m_{0}^{2}}{\rho}D_{0}\nabla
n_{0}-\frac{mm_{0}}{\rho}D\nabla n-\frac{\rho}{T}D^{T}\nabla
T-D^{U}\Delta\mathbf{U},$ (1)

where $\rho=mn$ is the total mass density, $n$ is the number density of the
gas particles, $n_{0}$ is the number density of tracer particles, $T$ is the
granular temperature, and $\Delta\mathbf{U}=\mathbf{U}-\mathbf{U}_{g}$;
$\mathbf{U}$ and $\mathbf{U}_{g}$ being the mean flow velocities of the
granular and the interstitial gas, respectively. In Eq. (1), $D_{0}$ is the
tracer diffusion coefficient, $D$ is the mutual diffusion coefficient, $D^{T}$
is the thermal diffusion coefficient, and $D^{U}$ is the velocity diffusion
coefficient. One of the main goals of the present paper is to obtain the
diffusion transport coefficients $D_{0}$, $D$, $D^{T}$, and $D^{U}$ up to the
second Sonine approximation in terms of the coefficients of restitution for
the tracer-grain ($\alpha_{0}$) and grain-grain ($\alpha$) collisions, the
masses and diameters, the solid volume fraction $\phi$ occupied by the
granular gas, and the background temperature $T_{\text{ex}}$. The expression
of the coefficient $D_{0}$ has been obtained in a previous work Gómez González
_et al._ (2023) where a complete study of the mean square displacement of an
intruder in a granular suspension has been carried out.

To check the accuracy of both Sonine approximations, the theoretical
predictions for the tracer diffusion coefficient $D_{0}$ are compared with
numerical solutions of the Enskog equation obtained from the direct simulation
Monte Carlo (DSMC) method. Bird (1994) As in previous papers, Garzó and
Montanero (2004); Garzó and Vega Reyes (2009, 2012); Gómez González  _et al._
(2023) the coefficient $D_{0}$ is computed from the Einstein formula
connecting the tracer diffusion coefficient with the mean square displacement
of intruders in a granular suspension in a homogeneous steady state. Here, to
complement the simulations performed in Ref. Gómez González  _et al._ , 2023
(which were carried out for mixtures with $\alpha=\alpha_{0}$) we consider
situations where both coefficients of restitution are different
($\alpha\neq\alpha_{0}$). Comparison between theory and simulations shows an
excellent agreement, especially in the case of the second Sonine approximation
for systems with high degree of inelasticity and/or very disparate mass and
diameter ratios.

The knowledge of the complete set of diffusion coefficients opens up the
possibility of applying our results to specific problems. As a first
application, we perform a linear stability analysis of the homogeneous steady
state of the mixture. This study is important by itself and also because this
state plays the role of the reference state in the Chapman–Enskog expansion.
Since the linearized hydrodynamic equations associated with the hydrodynamic
fields of the granular gas ($n$, $\mathbf{U}$, and $T$) are decoupled from the
tracer density $n_{0}$, the stability analysis (involving $d+3$ differential
equations) is quite simple since a previous study Gómez González and Garzó
(2019) has shown that the $d+2$ equations corresponding to the fields $n$,
$\mathbf{U}$, and $T$ are linearly stable. The fact that the coefficient
$D_{0}$ is positive implies automatically the stability of the homogeneous
steady state of the system (tracer plus granular gas). As a second
application, thermal diffusion segregation of an intruder immersed in a
granular suspension is analyzed. A segregation criterion showing the
transition between the well-known Brazil-nut effect (BNE) and the reverse
Brazil-nut effect (RBNE) by varying the parameters of the system is derived.
The present study complements a previous analysis recently carried in the low-
density regime. Gómez González and Garzó (2022a) As expected, our results show
that the form of the phase-diagrams for the BNE/RBNE transition depends
sensitively on the presence of the interstitial gas.

The plan of the paper is as follows. In Sec. II, the Enskog kinetic equation
for the granular suspension is introduced and the homogeneous steady state is
analyzed. Section III deals with the Enskog–Lorentz kinetic equation for the
tracer particles immersed in a granular suspension. We study first the
corresponding homogeneous steady state for intruders where it is shown that
their partial temperature is different from that of the granular gas. The set
of integral equations obeying the diffusion transport coefficients is
displayed in Sec. IV while the the first and second Sonine approximations to
these coefficients are obtained in Sec. V. These coefficients are explicitly
determined in terms of both the granular and the background temperatures, the
volume fraction, the coefficients of restitution, and the masses and diameters
of the mixture. Some technical details of the calculations are offered in four
appendices. The dependence of the diffusion coefficients on the parameter
space is illustrated in Sec. VI where the theoretical predictions of the
coefficient $D_{0}$ are also compared with Monte Carlo simulations. The
stability analysis of the homogeneous steady state is carried out in Sec. VI
while the thermal diffusion segregation is studied in Sec. VIII. We close the
paper in Sec. VII with a brief discussion of the results reported in this
work.

##  II Granular suspension

Let us consider a set of solid particles of mass $m$ and diameter $\sigma$
with collision rules according to the smooth hard sphere model. In this case,
collisions among solid particles are _inelastic_ and are characterized by a
positive constant coefficient of normal restitution $\alpha\leq 1$, where
$\alpha=1$ corresponds to elastic collisions (ordinary gases). We assume that
grains are immersed in a molecular gas of viscosity $\eta_{g}$. The system is
subjected to the action of the gravitational field and hence, each particle
feels the action of the force $m\mathbf{g}$ where $\mathbf{g}$ is the gravity
acceleration.

As said in Sec. I, we adopt here a coarse-grained description and the
influence of the interstitial gas on solid particles is via an instantaneous
fluid force. For low Reynolds numbers, this force is composed by two
independent terms. The first term accounts for the friction of grains on the
surrounding molecular gas and hence, it is proportional to the particle
velocity $\mathbf{v}$ (viscous drag force). The second term tries to model the
energy gained by the solid particles due to their “collisions” with the more
energetic particles of the gas phase. Williams and MacKintosh (1996) Both
terms are represented in the Enskog equation by a Fokker–Planck collision
operator. van Noije and Ernst (1998)

According to this way of modeling a granular suspension, at moderate
densities, the one-particle velocity distribution function of solid particles
$f(\mathbf{r},\mathbf{v};t)$ obeys the Enskog kinetic equation Gómez González
and Garzó (2019)

$\displaystyle\frac{\partial f}{\partial t}+\mathbf{v}\cdot\nabla
f-\gamma\Delta\mathbf{U}\cdot\frac{\partial
f}{\partial\mathbf{v}}+\mathbf{g}\cdot\frac{\partial
f}{\partial\mathbf{v}}-\gamma\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{V}f$
(2) $\displaystyle-\gamma\frac{T_{\text{ex}}}{m}\frac{\partial^{2}f}{\partial
v^{2}}=J[\mathbf{r},\mathbf{v}|f(t),f(t)],$

where $T_{\text{ex}}$ is the background temperature,
$\Delta\mathbf{U}=\mathbf{U}-\mathbf{U}_{g}$, $\mathbf{U}_{g}$ being the mean
fluid velocity of the gas phase and

$\mathbf{U}(\mathbf{r},t)=\frac{1}{n(\mathbf{r},t)}\int
d\mathbf{v}\;\mathbf{v}f(\mathbf{r},\mathbf{v},t)$ (3)

is the mean particle velocity. In addition, the Enskog collision operator
$J[\mathbf{r},\mathbf{v}|f,f]$ is

$\displaystyle J\left[\mathbf{r}_{1},\mathbf{v}_{1}\mid f(t),f(t)\right]$ (4)
$\displaystyle=$ $\displaystyle\sigma^{d-1}\int d\mathbf{v}_{2}\int
d\widehat{{\bm{\sigma}}}\Theta(\widehat{{\bm{\sigma}}}\cdot\mathbf{g}_{12})(\widehat{{\bm{\sigma}}}\cdot\mathbf{g}_{12})$
$\displaystyle\times\Big{[}\alpha^{-2}\chi\left(\mathbf{r}_{1},\mathbf{r}_{1}-{\bm{\sigma}}\right)f(\mathbf{r}_{1},\mathbf{v}_{1}^{\prime\prime};t)f(\mathbf{r}_{1}-{\bm{\sigma}},\mathbf{v}_{2}^{\prime\prime};t)$
$\displaystyle-\chi\left(\mathbf{r}_{1},\mathbf{r}_{1}+{\bm{\sigma}}\right)f(\mathbf{r}_{1},\mathbf{v}_{1};t)f(\mathbf{r}_{1}+{\bm{\sigma}},\mathbf{v}_{2};t)\Big{]},$

where $d$ is the dimensionality of the system ($d=2$ for disks and $d=3$ for
spheres), $\bm{\sigma}=\sigma\widehat{\bm{\sigma}}$, $\widehat{\bm{\sigma}}$
being a unit vector, $\Theta$ is the Heaviside step function, and ${\bf
g}_{12}={\bf v}_{1}-{\bf v}_{2}$. The double primes on the velocities in the
Enskog operator (4) denote the initial values
$\\{\mathbf{v}_{1}^{\prime\prime},\mathbf{v}_{2}^{\prime\prime}\\}$ that lead
to $\\{\mathbf{v}_{1},\mathbf{v}_{2}\\}$ following a binary collision:

${\bf v}_{1}^{\prime\prime}={\bf
v}_{1}-\frac{1}{2}\left(1+\alpha^{-1}\right)(\widehat{{\bm{\sigma}}}\cdot{\bf
g}_{12})\widehat{\bm{\sigma}},$ (5) ${\bf v}_{2}^{\prime\prime}={\bf
v}_{2}+\frac{1}{2}\left(1+\alpha^{-1}\right)(\widehat{{\bm{\sigma}}}\cdot{\bf
g}_{12})\widehat{\bm{\sigma}}.$ (6)

In addition, $\chi[{\bf r},{\bf r}\pm\bm{\sigma}|\\{n(t)]$ is the equilibrium
pair correlation function at contact as a functional of the nonequilibrium
density field $n({\bf r},t)$ defined by

$n({\bf r},t)=\int\;d{\bf v}f({\bf r},{\bf v},t).$ (7)

Apart from $n$ and $\mathbf{U}$, the other relevant hydrodynamic field is the
_granular_ temperature $T(\mathbf{r},t)$ defined as

$T(\mathbf{r},t)=\frac{m}{dn(\mathbf{r},t)}\int
d\mathbf{v}\;V^{2}f(\mathbf{r},\mathbf{v},t),$ (8)

where $\mathbf{V}=\mathbf{v}-\mathbf{U}$ is the peculiar velocity.

As in previous works, Gómez González, Khalil, and Garzó (2020); Gómez González
and Garzó (2019) the drift (or friction) coefficient $\gamma$ is assumed to be
a scalar quantity for the sake of simplicity. It is proportional to the gas
viscosity $\eta_{g}$. Given that in the dilute regime every particle is only
subjected to its respective Stokes drag force, then for hard spheres ($d=3$)
$\gamma$ is defined as

$\gamma\equiv\gamma_{\text{St}}=\frac{3\pi\sigma\eta_{g}}{m}.$ (9)

For moderate densities and low Reynolds numbers, $\gamma$ can be written as

$\gamma=\gamma_{\text{St}}R(\phi),$ (10)

where $R(\phi)$ is a function of the solid volume fraction

$\phi=\frac{\pi^{d/2}}{2^{d-1}d\Gamma\left(\frac{d}{2}\right)}n\sigma^{d}.$
(11)

The density dependence of the dimensionless function $R$ can be inferred from
computer simulations. For instance, for hard spheres ($d=3$), a good
approximation for $R$ is given by Van der Hoef, Beetstra, and Kuipers (2005);
Beetstra, Van der Hoef, and Kuipers (2007); Yin and Sundaresan (2009)

$R(\phi)=\frac{10\phi}{(1-\phi)}+\left(1-\phi\right)^{3}\left(1+1.5\sqrt{\phi}\right).$
(12)

The macroscopic balance equations for the granular suspension are obtained
when one multiplies the Enskog equation (11) by $\\{1,m{\bf v},mv^{2}\\}$ and
integrates over velocity. After some algebra, one achieves the balance
equations

$D_{t}n+n\nabla\cdot{\bf U}=0\;,$ (13) $D_{t}{\bf
U}+\rho^{-1}\nabla\cdot\mathsf{P}=\mathbf{g}-\gamma\Delta\mathbf{U}\;,$ (14)
$D_{t}T+\frac{2}{dn}\left(\nabla\cdot{\bf q}+\mathsf{P}:\nabla{\bf
U}\right)=2\gamma\left(T_{\text{ex}}-T\right)-\zeta\,T.$ (15)

In the above equations, $D_{t}=\partial_{t}+{\bf U}\cdot\nabla$ is the
material derivative. The cooling rate $\zeta$ is proportional to
$1-\alpha^{2}$ and is due to dissipative collisions. The pressure tensor ${\sf
P}({\bf r},t)$ and the heat flux ${\bf q}({\bf r},t)$ have both kinetic and
collisional transfer contributions, i.e., ${\sf P}={\sf P}^{\text{k}}+{\sf
P}^{\text{c}}$ and ${\bf q}={\bf q}^{\text{k}}+{\bf q}^{\text{c}}$. The
microscopic expressions of $\zeta$, $\mathsf{P}$, and $\mathbf{q}$ in terms of
$f$ can be found in Ref. Gómez González, Khalil, and Garzó, 2020. They will be
omitted here for brevity.

###  II.1 Homogeneous steady state

Let us assume that the granular suspension is in a homogeneous steady state.
In this state, $\mathbf{U}=\mathbf{U}_{g}$, $\nabla n=\nabla T=0$,
$\partial_{t}n=\partial_{t}T=0$ and hence, there is an exact balance between
the cooling effects (coming from the viscous friction and the inelasticity of
collisions) and the energy gained by grains due to their interaction with
particles of the bath. The balance equation (15) for the granular temperature
$T$ reads

$2\gamma\left(T_{\text{ex}}-T\right)=T\zeta.$ (16)

For elastic collisions ($\alpha=1$), $\zeta=0$, and $T=T_{\text{ex}}$, as
expected. Moreover, it is quite simple to show that the Maxwellian
distribution at the temperature $T_{\text{ex}}$ is a solution of the Enskog
equation (2). This is a consequence of the fluctuation-dissipation theorem.
van Kampen (2007) For inelastic collisions ($\alpha\neq 1$), it follows from
dimensional analysis that $f(\mathbf{v})$ has the scaling form

$f(\mathbf{v})=n\pi^{-d/2}v_{\text{th}}^{-d}\varphi(\mathbf{c}),$ (17)

where $v_{\text{th}}=\sqrt{2T/m}$ is the thermal velocity of the granular gas
and $\mathbf{c}\equiv\mathbf{v}/v_{\text{th}}$. The exact form of the scaling
distribution $\varphi$ is not known to date. On the other hand, previous
calculations Gómez González and Garzó (2019) has shown that the fourth-
cumulant (measuring the deviation of $f$ from its Maxwellian form) is in
general very small. Thus, a good estimate for the temperature $T$ can be made
by approximating $f$ by the Maxwellian distribution

$f_{\text{M}}(\mathbf{v})=n\left(\frac{m}{2\pi
T}\right)^{d/2}\exp\left(-\frac{mv^{2}}{2T}\right).$ (18)

In this case, the cooling rate $\zeta$ is

$\zeta=\frac{\sqrt{2}\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}\chi(1-\alpha^{2})\nu,$
(19)

where

$\nu=n\sigma^{d-1}\sqrt{\frac{2T}{m}}$ (20)

is an effective collision frequency. Equation (16) can be more explicitly
written when one takes into account Eq. (19) for $\zeta$. In dimensionless
form, Eq. (16) can be rewritten as

$2\lambda\theta^{-1/2}\left(\theta^{-1}-1\right)=\zeta^{*},$ (21)

where $\theta\equiv T/T_{\text{ex}}$, $\zeta^{*}=\zeta/\nu$, and

$\lambda=\frac{\sqrt{2}\pi^{d/2}}{2^{d}d\Gamma\left(\frac{d}{2}\right)\phi}\frac{R(\phi)}{\sqrt{T_{\text{ex}}^{*}}}.$
(22)

In Eq. (22), the (reduced) bath temperature $T_{\text{ex}}^{*}\equiv
T_{\text{ex}}/(m\sigma^{2}\gamma_{\text{St}}^{2})$. Equation (21) is a cubic
equation for the (reduced) temperature $\theta$. It physical solution gives
$\theta$ in terms of the density $\phi$, the (reduced) background temperature
$T_{\text{ex}}^{*}$, and the coefficient of restitution $\alpha$. It must
satisfy the requirement $\theta=1$ for $\alpha=1$ and any value of $\phi$ and
$T_{\text{ex}}^{*}$. In contrast, for inelastic collisions, $\theta<1$ since
the granular temperature is smaller than that of the interstitial gas. More
explicitly, the physical solution to Eq. (21) is given by Gómez González  _et
al._ (2023)

$\theta=\frac{\left(\Xi^{1/3}+\Xi^{-1/3}-1\right)^{2}}{9\varepsilon^{2}},$
(23)

where $\varepsilon=\zeta^{*}/2\lambda$ and

$\Xi=\frac{3\sqrt{3}\sqrt{27\varepsilon^{4}-4\varepsilon^{2}}+27\varepsilon^{2}-2}{2}.$
(24)

For elastic collisions ($\alpha=1$), $\varepsilon\to 0$, and Eq. (22) yields
$\theta\to 1$ as expected.

##  III Tracer particles immersed in a granular suspension

Let us suppose now that a few tracer or intruder particles of mass $m_{0}$ and
diameter $\sigma_{0}$ are added to the granular suspension. Since the
concentration of the tracer particles is negligible, their presence does not
have any effect on the state of the granular suspension. This implies that the
velocity distribution function $f(\mathbf{v})$ of the granular particles still
obeys the (nonlinear) Enskog equation (2). In addition, the balance equations
for the macroscopic flow velocity $\mathbf{U}$ and temperature $T$ for the
multicomponent granular suspension composed by the (monocomponent) granular
suspension plus the tracer particles are the same as those for the granular
suspension, namely they are given by Eqs. (14) and (15), respectively. Under
these conditions, the velocity distribution function $f_{0}({\bf r},{\bf
v},t)$ of the tracer species satisfies the linear Enskog-Lorentz kinetic
equation

$\displaystyle\frac{\partial f_{0}}{\partial t}+\mathbf{v}\cdot\nabla
f_{0}-\gamma_{0}\Delta\mathbf{U}\cdot\frac{\partial
f_{0}}{\partial\mathbf{v}}+\mathbf{g}\cdot\frac{\partial
f_{0}}{\partial\mathbf{v}}-\gamma_{0}\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{V}f_{0}$
$\displaystyle-\gamma_{0}\frac{T_{\text{ex}}}{m_{0}}\frac{\partial^{2}f_{0}}{\partial
v^{2}}=J_{0}[\mathbf{r},\mathbf{v}|f_{0}(t),f(t)],$ (25)

where the Enskog–Lorentz collision operator $J_{0}[{\bf v}|f_{0}(t),f(t)]$ is

$\displaystyle J_{0}[{\bf r}_{1},{\bf v}_{1}|f_{0}(t),f(t)]$ $\displaystyle=$
$\displaystyle\overline{\sigma}^{d-1}\int d\mathbf{v}_{2}\int
d\widehat{{\bm{\sigma}}}\Theta(\widehat{{\bm{\sigma}}}\cdot\mathbf{g}_{12})(\widehat{{\bm{\sigma}}}\cdot\mathbf{g}_{12})$
$\displaystyle\times\Big{[}\alpha_{0}^{-2}\chi_{0}\left(\mathbf{r}_{1},\mathbf{r}_{1}-\overline{{\bm{\sigma}}}\right)f_{0}(\mathbf{r}_{1},\mathbf{v}_{1}^{\prime\prime};t)f(\mathbf{r}_{1}-\overline{{\bm{\sigma}}},\mathbf{v}_{2}^{\prime\prime};t)$
$\displaystyle-\chi_{0}\left(\mathbf{r}_{1},\mathbf{r}_{1}+\overline{{\bm{\sigma}}}\right)f_{0}\left(\mathbf{r}_{1},\mathbf{v}_{1};t\right)f\left(\mathbf{r}_{1}+\overline{{\bm{\sigma}}},\mathbf{v}_{2};t\right)\Big{]}.$

Here, $\overline{{\bm{\sigma}}}=\overline{\sigma}\widehat{{\bm{\sigma}}}$,
$\overline{\sigma}=(\sigma_{0}+\sigma)/2$, $\alpha_{0}\leq 1$ is the
(positive) coefficient of restitution for tracer-gas collisions, and
$\chi_{0}$ is the pair correlation function for tracer-gas pairs at contact.
The relationship between the precollisional velocities
$\left\\{\mathbf{v}_{1}^{\prime\prime},\mathbf{v}_{2}^{\prime\prime}\right\\}$
and the postcollisional velocities
$\left\\{\mathbf{v}_{1},\mathbf{v}_{2}\right\\}$ is

${\bf v}_{1}^{\prime\prime}={\bf
v}_{1}-M\left(1+\alpha_{0}^{-1}\right)(\widehat{{\bm{\sigma}}}\cdot{\bf
g}_{12})\widehat{{\bm{\sigma}}},$ (27) ${\bf v}_{2}^{\prime\prime}={\bf
v}_{2}+M_{0}\left(1+\alpha_{0}^{-1}\right)(\widehat{{\bm{\sigma}}}\cdot{\bf
g}_{12})\widehat{\bm{\sigma}},$ (28)

where $M=m/(m+m_{0})$ and $M_{0}=m_{0}/(m+m_{0})$. To summarize and clarify, a
collision diagram along with the forces acting on the grains can be seen in
Figure 1.

The number density for the tracer particles is defined as

$n_{0}({\bf r},t)=\int\;d{\bf v}f_{0}({\bf r},{\bf v},t).$ (29)

The tracer or intruder may freely lose or gain momentum and energy in its
interactions with the particles of the granular gas. This means that these
quantities are not invariants of the collision operator $J_{0}[{\bf
v}|f_{0},f]$. Only the number density $n_{0}$ is conserved; its continuity
equation is directly obtained from the Enskog–Lorentz equation (15):

$D_{t}n_{0}+n_{0}\nabla\cdot{\bf U}+\frac{\nabla\cdot{\bf j}_{0}}{m_{0}}=0\;,$
(30)

where

${\bf j}_{0}(\mathbf{r},t)=m_{0}\int d{\bf v}\,{\bf V}\,f_{0}({\bf r},{\bf
v},t)$ (31)

is the mass flux for the tracer particles relative to the mean local flow
${\bf U}$ of the granular gas. Apart from $n_{0}$, an interesting quantity at
a kinetic level is the local temperature of the intruder $T_{0}$. It is
defined as

$T_{0}({\bf r},t)=\frac{m_{0}}{dn_{0}({\bf r},t)}\int\;d{\bf
v}\,V^{2}f_{0}({\bf r},{\bf v},t).$ (32)

The partial temperature $T_{0}$ measures the mean kinetic energy of the
intruders; it is in general different from that of the granular gas ($T\neq
T_{0}$).

The friction coefficient $\gamma_{0}$ takes into account the “interaction”
between the tracer particles and the background molecular gas. As the
coefficient $\gamma$, it can be written as
$\gamma_{0}=\gamma_{0,\text{St}}R_{0}(\phi)$ where for hard spheres ($d=3$)

$\gamma_{0,\text{St}}=\frac{3\pi\sigma_{0}\eta_{g}}{m_{0}}=\frac{\sigma_{0}m}{\sigma
m_{0}}\gamma_{\text{St}}.$ (33)

As in the case of the function $R$, one must chose an specific form of $R_{0}$
to get explicit results. Computer simulations Van der Hoef, Beetstra, and
Kuipers (2005); Beetstra, Van der Hoef, and Kuipers (2007); Yin and Sundaresan
(2009) indicate that a reliable choice for the function $R_{0}$ is

$R_{0}=1+\left(R-1\right)\left[a\frac{\sigma_{0}}{\sigma}+(1-a)\frac{\sigma_{0}^{2}}{\sigma^{2}}\right],$
(34)

where

$a(\phi)=1-2.660\phi+9.096\phi^{2}-11.338\phi^{3}.$ (35)

Note that for mechanically equivalent particles, $R_{0}=R$ as expected. It
must remarked that our results apply for any specific choice of the functions
$R$ and $R_{0}$.

Figure 1: (color online) A sketch representation of a collision between a
particle of species $i$ and a particle of species $j$ together with the drag
and stochastic gas-solid forces acting on these particles.

###  III.1 Homogeneous steady state

As in the case of the granular gas, the study of the homogeneous steady state
for the tracer particles is a crucial step since this state plays the role of
the reference state in the Chapman–Enskog perturbative method. Chapman and
Cowling (1970) In the absence of spatial gradients, the Enskog–Lorentz
equation (III) becomes

$-\gamma_{0}\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{V}f_{0}-\gamma_{0}\frac{T_{\text{ex}}}{m_{0}}\frac{\partial^{2}f_{0}}{\partial
v^{2}}=J_{0}[f_{0}(t),f(t)],$ (36)

where use has been made of the identity $\mathbf{U}=\mathbf{U}_{g}$ for
homogeneous states. The equation for the (steady) partial temperature $T_{0}$
can be easily obtained from Eq. (36) as

$2\gamma_{0}\left(T_{\text{ex}}-T_{0}\right)=T_{0}\zeta_{0},$ (37)

where

$\zeta_{0}=-\frac{1}{dn_{0}T_{0}}\int d\mathbf{v}\;m_{0}v^{2}\;J_{0}[f_{0},f]$
(38)

is the partial cooling rate characterizing the rate of energy dissipated by
tracer-grain collisions. As in the case of the granular suspension, for
elastic collisions ($\alpha_{0}=\alpha=1$), $\zeta_{0}=0$,
$T_{\text{ex}}=T_{0}$, and the Maxwellian distribution is an exact solution of
Eq. (36). Again, for inelastic collisions ($\alpha\neq 1$) the solution is not
known to date. However, since the corresponding cumulant associated with the
distribution $f_{0}(\mathbf{v})$ is in general very small, Gómez González and
Garzó (2021) an accurate estimate of $\zeta_{0}$ is obtained from the
Maxwellian distribution defined at the partial temperature $T_{0}$, i.e.,

$f_{0,\text{M}}(\mathbf{v})=n_{0}\left(\frac{m_{0}}{2\pi
T_{0}}\right)^{d/2}\exp\left(-\frac{mv^{2}}{2T_{0}}\right).$ (39)

In this approximation,

$\displaystyle\zeta_{0}$ $\displaystyle=$
$\displaystyle\frac{4\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}\chi_{0}^{(0)}M\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\left(\frac{1+\beta}{\beta}\right)^{1/2}(1+\alpha_{0})$
(40)
$\displaystyle\times\left[1-\frac{1}{2}M(1+\beta)(1+\alpha_{0})\right]\nu,$

where $\chi_{0}^{(0)}$ is the tracer-gas pair correlation at zeroth-order and

$\beta=\frac{m_{0}T}{mT_{0}}$ (41)

is the ratio between the mean square velocities of intruders and grains.

In dimensionless form, Eq. (37) reads

$2\lambda_{0}\theta^{-1/2}\left(\theta^{-1}-\tau_{0}\right)=\zeta_{0}^{*}\tau_{0},$
(42)

where $\tau_{0}=T_{0}/T$ is the temperature ratio,
$\zeta_{0}^{*}=\zeta_{0}/\nu$, and

$\lambda_{0}=\frac{\sigma_{0}m}{\sigma m_{0}}\frac{R_{0}}{R}\lambda.$ (43)

We recall that $\lambda$ is defined by Eq. (22). When intruder and granular
gas particles are mechanically equivalent ($m=m_{0}$, $\sigma=\sigma_{0}$, and
$\alpha=\alpha_{0}$), $\lambda=\lambda_{0}$, $\zeta^{*}=\zeta_{0}^{*}$,
$T=T_{0}$, and hence energy equipartition applies. However, in the general
case (namely, when collisions are inelastic and intruder and grains are
mechanically different), the solution to the cubic equation (42) provides the
temperature ratio $\tau_{0}$ in terms of the parameter space of the system. As
in the free cooling case, there is a breakdown of the energy equipartition
($\tau_{0}\neq 1$) as expected.

##  IV Diffusion transport coefficients

The diffusion transport coefficients associated with the mass flux
$\mathbf{j}_{0}$ of intruders can be obtained by solving the Enskog–Lorentz
equation (III) by means of the Chapman–Enskog method Chapman and Cowling
(1970) conveniently adapted to account for the inelastic character of
collisions. Garzó (2019) This long and hard task has been recently carried out
in Ref. Gómez González, Khalil, and Garzó, 2020 up to first-order in spatial
gradients and arbitrary concentration. Using those results, we consider here
the tracer limit ($x_{0}\equiv n_{0}/n\to 0$) of the linear integral equations
obeying the corresponding diffusion transport coefficients.

The first-order contribution $\mathbf{j}_{0}^{(1)}$ to the mass flux is given
by Eq. (1) where the diffusion transport coefficients $D^{T}$, $D_{0}$, $D$,
and $D^{U}$ are defined, respectively, as

$D^{T}=-\frac{m_{0}}{\rho d}\int
d\mathbf{v}\mathbf{V}\cdot\bm{\mathcal{A}}_{0}\left(\mathbf{V}\right),$ (44)
$D_{0}=-\frac{\rho}{m_{0}n_{0}d}\int
d\mathbf{v}\mathbf{V}\cdot\bm{\mathcal{B}}_{0}\left(\mathbf{V}\right),$ (45)
$D=-\frac{1}{d}\int
d\mathbf{v}\mathbf{V}\cdot\bm{\mathcal{C}}_{0}\left(\mathbf{V}\right),$ (46)

and

$D^{U}=-\frac{m_{0}}{d}\int
d\mathbf{v}\mathbf{V}\cdot\bm{\mathcal{E}}_{0}\left(\mathbf{V}\right).$ (47)

For arbitrary concentration $x_{0}$, the quantities
$\bm{\mathcal{A}}_{0}(\mathbf{V})$,
$\bm{\mathcal{B}}_{0}\left(\mathbf{V}\right)$,
$\bm{\mathcal{C}}_{0}\left(\mathbf{V}\right)$, and
$\bm{\mathcal{E}}_{0}(\mathbf{V})$ are the solutions of the set of coupled
linear integral equations given by Eqs. (73)–(74) and (77), respectively, of
Ref. Gómez González, Khalil, and Garzó, 2020. In order to write this set of
integral equations in the limit $x_{0}\to 0$, one has to take into account
that while in the tracer limit $D_{0}$ is independent of $x_{0}$, the
coefficients $D$, $D^{T}$, and $D^{U}$ are proportional to $x_{0}$. This
dependence must be then self-consistently confirmed. This means that
$\bm{\mathcal{A}}_{0}\propto x_{0}$, $\bm{\mathcal{C}}_{0}\propto x_{0}$, and
$\bm{\mathcal{E}}_{0}\propto x_{0}$. Thus, in the tracer limit ($x_{0}\to 0$),
the set of coupled linear integral equations for the unknowns
$\bm{\mathcal{A}}_{0}(\mathbf{V})$,
$\bm{\mathcal{B}}_{0}\left(\mathbf{V}\right)$,
$\bm{\mathcal{C}}_{0}\left(\mathbf{V}\right)$, and
$\bm{\mathcal{E}}_{0}(\mathbf{V})$ are:

$-\Bigg{(}2\gamma\theta^{-1}+\frac{1}{2}\zeta^{(0)}\Bigg{)}\bm{\mathcal{A}}_{0}-\gamma_{0}\frac{\partial}{\partial{\bf
V}}\cdot\left({\bf
V}\bm{\mathcal{A}}_{0}\right)-\gamma_{0}\frac{T_{\text{ex}}}{m_{0}}\frac{\partial^{2}\bm{\mathcal{A}}_{0}}{\partial
v^{2}}-J_{0}^{(0)}[\bm{\mathcal{A}}_{0},f^{(0)}]=\mathbf{A}_{0}+J_{0}^{(0)}[f_{0}^{(0)},\bm{\mathcal{A}}],$
(48) $-\gamma_{0}\frac{\partial}{\partial{\bf V}}\cdot\left({\bf
V}\bm{\mathcal{B}}_{0}\right)-\gamma_{0}\frac{T_{\text{ex}}}{m_{0}}\frac{\partial^{2}\bm{\mathcal{B}}_{0}}{\partial
v^{2}}-J_{0}^{(0)}[\bm{\mathcal{B}}_{0},f^{(0)}]=\mathbf{B}_{0},$ (49)
$-\gamma_{0}\frac{\partial}{\partial{\bf V}}\cdot\left({\bf
V}\bm{\mathcal{C}}_{0}\right)-\gamma_{0}\frac{T_{\text{ex}}}{m_{0}}\frac{\partial^{2}\bm{\mathcal{C}}_{0}}{\partial
v^{2}}-J_{0}^{(0)}[\bm{\mathcal{C}}_{0},f^{(0)}]=\mathbf{C}_{0}+J_{0}^{(0)}[f_{0}^{(0)},\bm{\mathcal{C}}]+\Bigg{[}n\frac{\partial\zeta^{(0)}}{\partial
n}-2n\gamma\left(\theta^{-1}-1\right)\frac{\partial\ln\gamma}{\partial
n}\Bigg{]}\bm{\mathcal{A}}_{0},$ (50) $-\gamma_{0}\frac{\partial}{\partial{\bf
V}}\cdot\left({\bf
V}\bm{\mathcal{E}}_{0}\right)-\gamma_{0}\frac{T_{\text{ex}}}{m_{0}}\frac{\partial^{2}\bm{\mathcal{E}}_{0}}{\partial
v^{2}}-J_{0}^{(0)}[{\bm{\mathcal{E}}}_{0},f^{(0)}]=\mathbf{E}_{0}.$ (51)

In the integral equations (48)–(51), $\zeta^{(0)}$ is the zeroth-order
approximation to the cooling rate $\zeta$, and $J_{0}^{(0)}[X,Y]$ is the
operator

$\displaystyle J_{0}^{(0)}\left[{\bf v}_{1}|X,Y\right]$ $\displaystyle=$
$\displaystyle\chi_{0}^{(0)}\overline{\sigma}^{d-1}\int d\mathbf{v}_{2}\int
d\widehat{{\bm{\sigma}}}\Theta(\widehat{{\bm{\sigma}}}\cdot\mathbf{g}_{12})(\widehat{{\bm{\sigma}}}\cdot\mathbf{g}_{12})$
$\displaystyle\times\left[\alpha_{0}^{-2}X(\mathbf{V}_{1}^{\prime\prime})Y(\mathbf{V}_{2}^{\prime\prime})-X(\mathbf{V}_{1})Y(\mathbf{V}_{2})\right],$

where $\chi_{0}^{(0)}$ is the intruder-gas pair correlation function at zeroth
order. The inhomogeneous terms appearing in Eqs. (48)–(51) are given by

$\mathbf{A}_{0}\left(\mathbf{V}\right)=-\mathbf{V}T\frac{\partial
f_{0}^{(0)}}{\partial T}-\frac{p}{\rho}\frac{\partial
f_{0}^{(0)}}{\partial\mathbf{V}}-\bm{\mathcal{K}}_{0}\Bigg{[}T\frac{\partial
f^{(0)}}{\partial T}\Bigg{]},$ (53) ${\bf B}_{0}\left(\mathbf{V}\right)=-{\bf
V}n_{0}\frac{\partial f_{0}^{(0)}}{\partial n_{0}},$ (54)
$\displaystyle\mathbf{C}_{0}\left(\mathbf{V}\right)$ $\displaystyle=$
$\displaystyle-\mathbf{V}n\frac{\partial f_{0}^{(0)}}{\partial
n}-m^{-1}\frac{\partial p}{\partial n}\frac{\partial
f_{0}^{(0)}}{\partial\mathbf{V}}-\frac{(1+\omega)^{-d}}{\chi_{0}^{(0)}T}$ (55)
$\displaystyle\times\left(\frac{\partial\mu_{0}}{\partial\phi}\right)_{T,n_{0}}\bm{\mathcal{K}}_{0}\left[f^{(0)}\right],$
$\mathbf{E}_{0}(\mathbf{V})=\left(\gamma_{0}-\gamma\right)\frac{\partial
f_{0}^{(0)}}{\partial\mathbf{V}}.$ (56)

In Eqs. (53)–(56),

$p=nT\left[1+2^{d-2}\chi^{(0)}\phi(1+\alpha)\right]$ (57)

is the pressure, $\omega\equiv\sigma_{0}/\sigma$ is the size ratio, $\mu_{0}$
is the chemical potential of the intruder and the operator
$\bm{\mathcal{K}}_{0}[X]$ is defined as

$\displaystyle\bm{\mathcal{K}}_{0}[X]$ $\displaystyle=$
$\displaystyle\overline{\sigma}^{d}\chi_{0}^{(0)}\int d\mathbf{v}_{2}\int
d\widehat{\bm{\sigma}}\Theta(\widehat{\bm{\sigma}}\cdot\mathbf{g}_{12})(\widehat{\bm{\sigma}}\cdot\mathbf{g}_{12})\bm{\widehat{\sigma}}$
$\displaystyle\times\left[\alpha_{0}^{-2}f_{0}^{(0)}(\mathbf{V}_{1}^{\prime\prime})X(\mathbf{V}_{2}^{\prime\prime})+f_{0}^{(0)}(\mathbf{V}_{1})X(\mathbf{V}_{2})\right].$

In addition, upon writing Eq. (55), we have neglected the non-Gaussian
corrections to $f^{(0)}$ and $f_{0}^{(0)}$, have taken into account that
$I_{121}\propto x_{0}\to 0$ and have used the explicit form of $I_{122}$ given
by Garzó, Dufty, and Hrenya (2007); Garzó (2019)

$I_{122}=\frac{d\Gamma\left(\frac{d}{2}\right)}{\pi^{d/2}T\overline{\sigma}^{d}\chi_{0}^{(0)}}\left(\frac{\partial\mu_{0}}{\partial
n}\right)_{T,n_{0}}-2-n\frac{\partial\ln\chi_{0}^{(0)}}{\partial n}.$ (59)

The explicit expressions of the quantities $I_{ij\ell}$ for a binary mixture
can be found in the Appendix A of chapter 5 of Ref. Garzó, 2019.

Note that in Eqs. (48), (50), and (51), the quantities
$\bm{\mathcal{A}}_{0}(\mathbf{V})$ and $\bm{\mathcal{C}}_{0}(\mathbf{V})$
associated with the tracer particles are coupled with their corresponding
counterparts $\bm{\mathcal{A}}(\mathbf{V})$ and $\bm{\mathcal{C}}(\mathbf{V})$
of the host granular gas. A direct consequence of this coupling is that the
mass flux $\mathbf{j}_{0}^{(1)}$ inherits terms coming from the autonomous
host integral equations (those proportional to the spatial gradients $\nabla
n$ and $\nabla T$). The expressions of $\bm{\mathcal{A}}$ and
$\bm{\mathcal{C}}$ up to the second Sonine approximation are displayed in the
Appendix A.

##  V Second Sonine polynomial approximation

As for ordinary (elastic) mixtures, Chapman and Cowling (1970) the integral
equations (48)–(51) are solved by expanding the unknowns
$\bm{\mathcal{A}}_{0}$, $\bm{\mathcal{B}}_{0}$, $\bm{\mathcal{C}}_{0}$, and
$\bm{\mathcal{E}}_{0}$ in a complete set of orthogonal polynomials with a
Gaussian measure. In kinetic theory, one employs the Sonine (or generalized
Laguerre) polynomials. On the other hand, for practical purposes, to achieve
explicit results one has to resort to some kind of truncation in the series
expansion. Usually, the leading term in the above series expansion is
retained. In this paper, as said in Sec. I, to provide accurate results for
the diffusion transport coefficients for quite inelastic systems and/or
mixtures with quite disparate masses and diameters, we will consider two terms
in the series expansion. This approximation is referred to as the second
Sonine approximation. In this case, the quantities $\bm{\mathcal{A}}_{0}$,
$\bm{\mathcal{B}}_{0}$, $\bm{\mathcal{C}}_{0}$, and $\bm{\mathcal{E}}_{0}$ are
given by

$\bm{\mathcal{A}}_{0}({\bf V})\to-f_{0,M}({\bf
V})\left[\frac{\rho}{n_{0}T_{0}}{\bf V}D^{T}+a_{0}{\bf S}_{0}({\bf
V})\right],$ (60) $\bm{\mathcal{B}}_{0}({\bf V})\to-f_{0,M}({\bf
V})\left[\frac{m_{0}^{2}}{\rho T_{0}}{\bf V}D_{0}+b_{0}{\bf S}_{0}({\bf
V})\right],$ (61) $\bm{\mathcal{C}}_{0}({\bf V})\to-f_{0,M}({\bf
V})\left[\frac{m_{0}}{n_{0}T_{0}}{\bf V}D+c_{0}{\bf S}_{0}({\bf V})\right],$
(62) $\bm{\mathcal{E}}_{0}\to-f_{0,M}({\bf V})\left[\frac{1}{n_{0}T_{0}}{\bf
V}D^{U}+e_{0}{\bf S}_{0}({\bf V})\right],$ (63)

where

${\bf S}_{0}({\bf
V})=\left(\frac{1}{2}m_{0}V^{2}-\frac{d+2}{2}T_{0}\right){\bf V},$ (64)

and $f_{0,M}({\bf V})$ is the Maxwellian distribution (39) with the change
$\mathbf{v}\to\mathbf{V}$. While the diffusion coefficients are defined by
Eqs. (44)–(47), the second Sonine coefficients $a_{0}$, $b_{0}$, $c_{0}$, and
$e_{0}$ are defined as

$a_{0}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{3}}\int\;d{\bf v}\;{\bf
S}_{0}({\bf V})\cdot\bm{\mathcal{A}}_{0}({\bf V}),$ (65)
$b_{0}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{3}}\int\;d{\bf v}\;{\bf
S}_{0}({\bf V})\cdot\bm{\mathcal{B}}_{0}({\bf V}),$ (66)
$c_{0}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{3}}\int\;d{\bf v}\;{\bf
S}_{0}({\bf V})\cdot\bm{\mathcal{C}}_{0}({\bf V}),$ (67)
$e_{0}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{3}}\int\;d{\bf v}\;{\bf
S}_{0}({\bf V})\cdot\bm{\mathcal{E}}_{0}({\bf V}).$ (68)

The transport coefficients $D_{0}$, $D$, $D^{T}$ and $D^{U}$ as well as the
second Sonine coefficients $a_{0}$, $b_{0}$, $c_{0}$ and $e_{0}$ are
determined by substitution of Eqs. (60)–(63) into the integral equations
(48)–(51), multiplication of these equations by $m_{0}{\bf V}$ and by ${\bf
S}_{0}({\bf V})$, and integration over velocity. Some technical details on
these calculations are given in the Appendix B.

To provide the final expressions, it is convenient to write them in
dimensionless form. Let us introduce the reduced transport coefficients

$D^{T*}=\frac{\rho\nu}{n_{0}T}D^{T},\quad D_{0}^{*}=\frac{m_{0}^{2}\nu}{\rho
T}D_{0},$ (69) $D^{*}=\frac{m_{0}\nu}{n_{0}T}D,\quad
D_{U}^{*}=\rho_{0}^{-1}D^{U},$ (70)

where $\rho_{0}=m_{0}n_{0}$ is the mass density of the intruders and $\nu$ is
defined by Eq. (20).

According to the results exposed in the Appendix B, the first Sonine
approximations $D^{T*}[1]$, $D_{0}^{*}[1]$, $D^{*}[1]$, and $D_{U}^{*}[1]$ can
be written as

$D^{T*}[1]=\frac{X_{1}^{*}}{\nu_{1}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}+\gamma_{0}^{*}},$
(71) $D_{0}^{*}[1]=\frac{\tau_{0}}{\nu_{1}^{*}+\gamma_{0}^{*}},$ (72)
$D^{*}[1]=\frac{Y_{1}^{*}}{\nu_{1}^{*}+\gamma_{0}^{*}},$ (73)
$D_{U}^{*}[1]=\frac{\gamma_{0}^{*}-\gamma^{*}}{\nu_{1}^{*}+\gamma_{0}^{*}}.$
(74)

In Eqs. (71)–(74), $X_{1}^{*}$, $a^{*}=\nu T^{2}a$, $Y_{1}^{*}$, and
$c^{*}=\nu nTc$ are given by Eqs. (128), (109), (139), and (110),
respectively. The dimensionless cooling rate $\zeta^{*}=\zeta^{(0)}/\nu$ is
defined by Eq. (19) while the expressions of the (reduced) collision
frequencies $\nu_{i}^{*}=\nu_{i}/\nu$ are displayed in the Appendix C. The
forms of the second Sonine approximations $D^{T*}[2]$, $D_{0}^{*}[2]$,
$D^{*}[2]$, and $D_{U}^{*}[2]$ are given by Eqs. (130), (133), (135), and
(141), respectively.

In the low-density limit ($\phi\to 0$), the first Sonine approximations to the
diffusion coefficients reduce to

$D^{T*}[1]=\frac{\tau_{0}+\theta\Delta_{\theta,0}-m_{0}/m}{\nu_{1}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}+\gamma_{0}^{*}},$
(75) $D_{0}^{*}[1]=\frac{\tau_{0}}{\nu_{1}^{*}+\gamma_{0}^{*}},$ (76)
$D^{*}[1]=\frac{\zeta^{*}D^{T^{*}}[1]-m_{0}/m}{\nu_{1}^{*}+\gamma_{0}^{*}},$
(77)

where $\nu_{1}^{*}$ is given by Eq. (142) with $\chi_{0}^{(0)}=1$ while Eq.
(154) gives the derivative
$\Delta_{\theta,0}\equiv\partial\tau_{0}/\partial_{\theta}$. This derivative
arises from the fact that the zeroth-order distribution $f_{0}^{(0)}$ depends
not only on hydrodynamic quantities but, in addition, on the partial
temperature $\tau_{0}$. Therefore, when performing the Chapman–Enskog
expansion around $f_{0}^{(0)}$, there are contributions to the transport
coefficients coming from the derivatives of the kinetic quantities in the
vicinity of the steady state. These contributions can be seen as a measure of
the departure of the perturbed time-dependent state from the steady reference
state. The expressions (75)–(77) are consistent with those obtained in a
previous work Gómez González and Garzó (2022a) from the Boltzmann equation.
111Note that there is a mistake in Eq. (35) for the coefficient $D^{*}$.
Equation (77) corrects this mistake.

It is quite apparent that, in general, the first and second Sonine
approximations to the diffusion transport coefficients have a complex
dependence on the coefficients of restitution, the mass and diameter ratios,
the solid volume fraction, and the background temperature. Thus, before
studying this dependence, it is quite instructive to consider some special
limits.

In the limiting case of mechanically equivalent particles ($m_{0}=m$,
$\sigma_{0}=\sigma$, $\alpha_{0}=\alpha$, $\gamma=\gamma_{0}$), as expected
one gets $D^{T*}[2]=D_{U}^{*}[2]=0$, $D_{0}^{*}[2]=-D^{*}[2]$ and so,

$\mathbf{j}_{0}^{(1)}=-\frac{nT}{\nu}D_{0}^{*}[2]\nabla x_{0}.$ (78)

The expression of $D_{0}^{*}[2]$ for arbitrary $d$ and finite $\phi$ is still
a very long expression. In the particular case of a dilute ($\phi=0$)
suspension of hard spheres ($d=3$), $D_{0}^{*}[2]$ is

$D_{0}^{*}[2]=\frac{98+\alpha\left[47+3\alpha_{0}(4\alpha-13)\right]+\frac{90\sqrt{2}}{\sqrt{\pi}}\gamma^{*}}{\frac{4\sqrt{2\pi}}{3}(1+\alpha)^{2}\left[50+\alpha(4\alpha-25)\right]+\gamma^{*}\left\\{(1+\alpha)\left[218+3\alpha(4\alpha-17)\right]+\frac{90\sqrt{2}}{\sqrt{\pi}}\gamma^{*}\right\\}}.$
(79)

Apart from $D_{0}^{*}[2]$, the first Sonine approximation $D_{0}^{*}[1]$
reaches the simple form

$D_{0}^{*}[1]=\frac{1}{\frac{2\sqrt{2\pi}}{3}(1+\alpha)+\gamma^{*}}.$ (80)

In the case of elastic collisions ($\alpha=\alpha_{0}=1$) but mechanically
different particles, $\tau_{0}=\theta=1$ and the tracer diffusion coefficient
for a three-dimensional dilute system reads

$D^{*}_{0}[2]=\frac{D^{*}_{0}[1]}{1-\left[\left(1+\frac{3}{2}\mu
M_{0}^{-1/2}(1+\omega)^{-2}\gamma_{0}^{*}\right)\left(30\mu^{2}+(16+\frac{45}{\sqrt{\pi}}M^{-2}M_{0}^{-1/2}(1+\omega)^{-2}\gamma_{0}^{*})\mu+13\right)\right]^{-1}},$
(81)

where $\mu\equiv m_{0}/m$ is the mass ratio and

$D^{*}_{0}[1]=\left[\frac{2\sqrt{\pi}}{3}\mu^{-1}M_{0}^{1/2}(1+\omega)^{2}+\gamma_{0}^{*}\right]^{-1}.$
(82)

Here, $\omega\equiv\sigma_{0}/\sigma$ is the diameter ratio. When
$\gamma^{*}=\gamma_{0}^{*}=0$, the expression (81) agrees with the one derived
many years ago for a gas mixture of elastic hard spheres. Mason (1954); López
de Haro and Cohen (1984) All these results show the self-consistency of the
present results with those reported in previous works in some limiting cases.

##  VI Some illustrative systems and comparison with computer simulations

Figure 2: (color online) Plot of reduced tracer diffusion coefficient
$D_{0}(\alpha)/D_{0}(1)$ as a function of the (common) coefficient of
restitution $\alpha$ for the systems $(m_{0}/m=4,\sigma_{0}/\sigma=2)$ (a) and
$(m_{0}/m=0.5,\sigma_{0}/\sigma=0.8)$ (b) in the case of a three-dimensional
gas with $\phi=0.1$ and $T_{\text{ex}}^{*}=0.1$. The solid lines correspond to
the second Sonine approximation while the dashed lines refer to the first
Sonine approximation. Here, $D_{0}(1)$ is the elastic value of the tracer
diffusion coefficient consistently obtained in each approximation.  Figure 3:
(color online) Plot of reduced mutual diffusion coefficient $D(\alpha)/D(1)$
as a function of the (common) coefficient of restitution $\alpha$ for the
systems $(m_{0}/m=4,\sigma_{0}/\sigma=2)$ (a) and
$(m_{0}/m=0.5,\sigma_{0}/\sigma=0.8)$ (b) in the case of a three-dimensional
gas with $\phi=0.1$ and $T_{\text{ex}}^{*}=0.1$. The solid lines correspond to
the second Sonine approximation while the dashed lines refer to the first
Sonine approximation. Here, $D(1)$ is the elastic value of the tracer
diffusion coefficient consistently obtained in each approximation. Figure 4:
(color online) Plot of reduced thermal diffusion coefficient
$D^{T}(\alpha)/D^{T}(1)$ as a function of the (common) coefficient of
restitution $\alpha$ for the systems $(\mu=4,\omega=2)$ (a) and
$(\mu=0.5,\omega=0.8)$ (b) in the case of a three-dimensional gas with
$\phi=0.1$ and $T_{\text{ex}}^{*}=0.1$. The solid lines correspond to the
second Sonine approximation while the dashed lines refer to the first Sonine
approximation. Here, $D^{T}(1)$ is the elastic value of the thermal diffusion
coefficient consistently obtained in each approximation.  Figure 5: (color
online) Plot of reduced velocity diffusion coefficient
$D^{U}(\alpha)/D^{U}(1)$ as a function of the (common) coefficient of
restitution $\alpha$ for the systems $(\mu=4,\omega=2)$ (a) and
$(\mu=0.5,\omega=0.8)$ (b) in the case of a three-dimensional gas with
$\phi=0.1$ and $T_{\text{ex}}^{*}=0.1$. The solid lines correspond to the
second Sonine approximation while the dashed lines refer to the first Sonine
approximation. Here, $D^{U}(1)$ is the elastic value of the velocity diffusion
coefficient consistently obtained in each approximation.

The expressions derived in Sec. IV for the (reduced) diffusion coefficients
depend on many parameters:
$\left\\{d,m_{0}/m,\sigma_{0}/\sigma,\alpha,\alpha_{0},\phi,T_{\text{ex}}^{*}\right\\}$.
This complexity appears in the elastic limit as well, except for the
dependence on the coefficients of restitution $\alpha$ and $\alpha_{0}$. Thus,
as we have done in previous papers, Garzó and Vega Reyes (2009, 2012) we
normalize the transport coefficients with respect to their values for elastic
collisions ($\alpha=\alpha_{0}=1$). This will show us more clearly the impact
of inelasticity in collisions on the diffusion transport coefficients. Also,
for the sake of simplicity, we consider the physical case of hard spheres
($d=3$) and in some plots we take a common coefficient of restitution
$\alpha=\alpha_{0}$. This reduces the parameter space of the system to five
independent quantities:
$\left\\{m_{0}/m,\sigma_{0}/\sigma,\alpha,\phi,T_{\text{ex}}^{*}\right\\}$.

In order to get the explicit dependence of the transport coefficients on the
above parameters one has to give the forms of the pair correlation functions
$\chi^{(0)}$ and $\chi_{0}^{(0)}$ for a three-dimensional system. In this
case, a good approximation for $\chi^{(0)}$ is provided by the Carnahan-
Starling form Carnahan and Starling (1969)

$\chi^{(0)}=\frac{1-\frac{1}{2}\phi}{(1-\phi)^{3}},$ (83)

while the intruder-gas pair correlation function is given by Boublik (1970)

$\chi_{0}^{(0)}=\frac{1}{1-\phi}+3\frac{\omega}{1+\omega}\frac{\phi}{(1-\phi)^{2}}+2\frac{\omega^{2}}{(1+\omega)^{2}}\frac{\phi^{2}}{(1-\phi)^{3}}.$
(84)

The expression for the chemical potential of the intruder consistent with the
approximation (84) is Reed and Gubbins (1973)

$\displaystyle\frac{\mu_{0}}{T}$ $\displaystyle=$ $\displaystyle C_{3}+\ln
n_{0}-\ln(1-\phi)+3\omega\frac{\phi}{1-\phi}$ (85)
$\displaystyle+3\omega^{2}\left[\ln(1-\phi)+\frac{\phi(2-\phi)}{(1-\phi)^{2}}\right]$
$\displaystyle-\omega^{3}\left[2\ln(1-\phi)+\frac{\phi(1-6\phi+3\phi^{2})}{(1-\phi)^{3}}\right],$

where $C_{3}$ is a constant.

Figures 2–5 show the dependence of the diffusion transport coefficients on
inelasticity for two different systems. As said before, each transport
coefficient has been reduced with respect to its elastic value consistently
obtained in each approximation. The dashed and solid lines refer to the first
and second Sonine approximations, respectively. We observe first that in
general both Sonine approximations yield quite similar results even for strong
inelasticity and quite disparate values of the mass and/or diameter ratios.
These results clearly differ from those reported in the case of _dry_ (no gas
phase) granular mixtures [see, for instance, Figs. 1-3 of Ref. Garzó and Vega
Reyes, 2009 for the same mixture parameters] where both Sonine approximations
lead to quite different results in the range of large inelasticity when the
tracer particles are lighter than the particles of the granular gas. This
means that the Sonine polynomial expansion exhibits a better convergence for
granular suspensions than in the dry granular case.

With respect to the influence of inelasticity on diffusion coefficients, we
observe that the impact of $\alpha$ on the diffusion transport coefficients is
in general relevant (since their forms differ clearly from their elastic
counterparts), although less important than in the dry granular case. The
comparison of Figs. 2–4 with Figs. 1–3 of Ref. Garzó and Vega Reyes, 2009
reveals important differences for diffusion between the results obtained with
and without the gas-phase. Thus, while the (scaled) tracer diffusion
coefficient $D_{0}(\alpha)/D_{0}(1)$ exhibits in some cases a non-monotonic
dependence on inelasticity for gas-solid flows, this reduced coefficient
increases with decreasing $\alpha$ in the dry granular limit. In addition,
Fig. 2 also shows the tiny dependence of $D_{0}(\alpha)/D_{0}(1)$ on both the
mass and diameter ratios. In the case of the (scaled) coefficients
$D(\alpha)/D(1)$ and $D^{T}(\alpha)/D^{T}(1)$, qualitative differences between
the predictions with and without the gas phase are observed for large mass
and/or diameter ratios. At a more quantitative level, we find that the effect
of $\alpha$ on those coefficients is more significant in dry granular
mixtures. The reduced velocity diffusion coefficient $D^{U}(\alpha)/D^{U}(1)$
is not present in the dry case. Figure 5 shows in general a weak influence of
inelasticity on this coefficient.

There are several competing effects involved in the mechanisms in the
diffusion of tracer particles in a granular suspension. On the one hand, the
collisional effects are accounted for in the cooling rate $\zeta^{*}$ and the
(reduced) collision frequencies $\nu_{i}^{*}$. On the other hand, the
influence of the interstitial gas on mass transport appears through the
(reduced) friction coefficients $\gamma^{*}$ and $\gamma_{0}^{*}$, which forms
also depend on the volume fraction $\phi$, the mass and diameter ratios, and
the (reduced) background temperature $T_{\text{ex}}^{*}$. Given the size of
the parameter space of the problem (six independent parameters), it is
generally difficult to predict which effect (collisions between grains or
interstitial gas) dominates.

In summary, as already noted in previous papers on granular mixtures, Gómez
González, Khalil, and Garzó (2020) the diffusion transport coefficients for a
binary granular suspension where one of the species is present in tracer
concentration differ in general significantly from those obtained for dry
granular mixtures. Garzó, Dufty, and Hrenya (2007); Garzó, Hrenya, and Dufty
(2007); Garzó and Vega Reyes (2009) These differences between both systems
increase with increasing inelasticity in some cases, and depending on the
cases, the influence of both mass and diameter ratios can be important.

###  VI.1 Tracer diffusion coefficient. DSMC results

It is quite apparent that the accuracy of the Sonine approximations must be
assessed via a comparison with computer simulation results. Here, as in a
previous work, Gómez González  _et al._ (2023) we have numerically solved the
Enskog equation by means of the DSMC method Bird (1994) for the tracer
particles immersed in a granular suspension in the homogeneous steady state.
In this situation, $\Delta\mathbf{U}=\mathbf{0}$, $\nabla n=\nabla T=0$, and
hence Eq. (30) yields

$\frac{\partial x_{0}}{\partial t}=2d\frac{m_{0}D_{0}}{\rho}\nabla^{2}x_{0},$
(86)

where use has been made of the constitutive equation (1). Equation (86) is a
standard diffusion equation with the time-independent diffusion coefficient
$m_{0}D_{0}/\rho$. It follows that the mean square displacement of the tracer
particles after a time interval $t$ is McLennan (1989)

$\frac{\partial}{\partial
t}\langle|\mathbf{r}(t)-\mathbf{r}(0)|^{2}\rangle=2d\frac{m_{0}D_{0}}{\rho}.$
(87)

Here, $|\mathbf{r}(t)-\mathbf{r}(0)|$ is the distance traveled by the tracer
particles from $t=0$ until time $t$. Equation (87) is the Einstein formula
relating the tracer diffusion coefficient $D_{0}$ with the mean square
displacement. The relation (87) is also employed in Monte Carlo simulations to
obtain the coefficient $D_{0}$.

Since we want to numerically solve the Enskog–Lorentz kinetic equation (III)
under steady conditions, we slightly adapt the method proposed by Montanero
and Garzó Montanero and Garzó (2002) oriented to model freely cooling granular
mixtures. The main changes to be taken into account are three: (i) the
intruders are in tracer concentration, and so we do not need to account for
collisions among intruder particles themselves; (ii) if a granular-intruder
collision is accepted, only the velocities of the intruders are updated; and
(iii) to account for the action of the interstitial fluid on grains we update
the velocity $\mathbf{v}_{k}$ of every single grain of species $i$ (intruder
or granular particles) after each time step $\delta t$ according to the rule:
Khalil and Garzó (2014)

$\mathbf{v}_{k}\to e^{-\gamma_{i}\delta
t}\mathbf{v}_{k}+\left(\frac{6\gamma_{i}T_{\text{ex}}\delta
t}{m_{i}}\right)^{1/2}\bm{\varpi}_{k}.$ (88)

Here, $\bm{\varpi}_{k}$ is a random vector of zero mean and unit variance. The
Fokker-Planck equation is recovered when the time step $\delta t$ is much
shorter than the mean free time between collisions. Khalil and Garzó (2014)

Figure 6: (color online) Plot of the temperature ratio $T_{0}/T$ as a function
of the intruder-grain coefficient of restitution $\alpha_{0}$ for the systems
($\mu=8$, $\omega=2$) (a) and ($\mu=0.5$, $\omega=0.5^{1/3}$) (b) in the case
of a three-dimensional gas with $\alpha=0.9$, $\phi=0.2$, and
$T^{*}_{\text{ex}}=1$. The solid lines correspond to the theoretical results
whereas the symbols refer to DSMC results. Figure 7: (color online) Plot of
the temperature ratio $T_{0}/T$ as a function of the intruder-grain
coefficient of restitution $\alpha_{0}$ for the systems ($\mu=10$, $\omega=5$)
(a) and ($\mu=0.5$, $\omega=0.5$) (b) in the case of a three-dimensional gas
with $\alpha=0.9$, $\phi=0.2$, and $T^{*}_{\text{ex}}=1$. The solid lines
correspond to the theoretical results whereas the symbols refer to DSMC
results.

To complement the simulation results displayed in Ref. Gómez González  _et
al._ , 2023 that were obtained considering a common coefficient of restitution
$\alpha=\alpha_{0}$, we simulate here systems where $\alpha\neq\alpha_{0}$.
Moreover, to assess the influence of $\alpha_{0}$ in the tracer diffusion
coefficient $D_{0}$, the same mixtures as in Sec. V.C of Ref. Gómez González
_et al._ , 2023 are studied. However, we consider denser suspensions
($\phi=0.2$) to enlarge the difference between the first and second Sonine
approximations to $D_{0}$.

Figure 6 depicts the temperature ratio $T_{0}/T$ versus the coefficient of
restitution $\alpha_{0}$ for $d=3$, $\alpha=0.9$, $T^{*}_{\text{ex}}=1$, and
$\phi=0.2$. Two different mixtures with the same mass density
($m/\sigma^{3}=m_{0}/\sigma_{0}^{3}$) are considered. As in the case of a
common coefficient of restitution, we find a small influence of the mechanical
properties (masses and sizes) on $T_{0}/T$ when the mass densities are equal.
However, we observe an enhance breakdown of energy equipartition when
$\alpha\neq\alpha_{0}$, as expected. This occurs because $T_{0}\sim T$ when
$\alpha_{0}\sim\alpha$ since the loss of energy is comparable in both species.
Using similar arguments to those presented in Sec. VI.B of Ref. Gómez González
_et al._ , 2023, we conclude that for moderately dense suspensions, the action
of the external bath is accentuated when $\omega$ increases and/or $\mu$
decreases. Although a combined effect of both masses and sizes is observed, we
can conclude that the influence of the diameter ratio $\omega$ on the
breakdown of energy equipartition is more relevant than the one caused by the
mass ratio $\mu$. This is due to the action of the interstitial gas resulting
from gas-solid collisions, whose dependence on particle size
($\gamma_{0,\text{St}}\propto\sigma_{0}/\sigma$) can be inferred from the
change of an effective surface, which overcomes the inertial forces determined
by particle masses. For this reason, the bath thermalizes the larger species
more effectively causing its temperature to rise and therefore, the curve (a)
of Fig. 6 ($\omega>1$) appears above the curve (b) of Fig. 6 ($\omega<1$).
Furthermore, we note that a departure from energy equipartition occurs when
$\alpha_{0}=1$. In contrast to Fig. 4 of Ref. Gómez González  _et al._ , 2023,
where $T_{\text{ex}}=T_{0}=T$ for $\alpha_{0}=1$, here we still observe energy
dissipation since $\alpha\neq 1$. Specifically, the loss of kinetic energy
when $\alpha_{0}=1$ ($\alpha=0.9$) is greater for grains than for intruders,
resulting in $T_{0}>T$. Figure 6 also highlights an excellent agreement
between theory and DSMC results, ensuring the reliability of the Maxwellian
approximations (18) and (39) taken for the distribution functions.

We now consider two different mixtures in Fig. 7. In contrast to the prior
case depicted in Fig. 6, the systems under consideration possess different
mass densities, accentuating the disparity between the temperature ratios.
Similar to the observations made in Fig. 6, we note that for larger intruders,
the thermalization of the temperature by the bath occurs more efficiently.
Thus, the action of the bath facilitates the classification of diverse
mixtures based on their mechanical properties and deviation from energy
equipartition. For moderate mass densities ($\mu\sim\omega$), the curves
exhibited in both Figs. 6 and 7 display a bottom-to-top ordering as the value
of $\omega$ increases. We can also analyze the temperature ratio $T_{0}/T$ in
the specific case in which $\alpha=\alpha_{0}=0.9$. While energy equipartition
practically holds for the system considered in Fig. 6, a more noticeable
breakdown appears in Fig. 7. This result has been previously discussed in Ref.
Gómez González  _et al._ , 2023 and is due in part to the mechanical
differences between the mixtures. We can conclude again that the theoretical
results perfectly reproduce the DSMC simulations, as demonstrated in Fig. 7.

Figure 8: (color online) Plot of the reduced tracer diffusion coefficient
$D_{0}(\alpha_{0})/D_{0}(1)$ as a function of the intruder-grain coefficient
of restitution $\alpha_{0}$ for the systems ($\mu=8$, $\omega=2$) (a) and
($\mu=0.5$, $\omega=0.5^{1/3}$) (b) in the case of a three-dimensional gas
with $\alpha=0.9$, $\phi=0.2$, and $T^{*}_{\text{ex}}=1$. Solid and dashed
lines are for the second and first Sonine approximations, respectively,
whereas the symbols refer to DSMC results. Here, $D_{0}(1)$ is the value of
the tracer diffusion coefficient consistently obtained in each approximation
when $\alpha_{0}=1$. Figure 9: (color online) Plot of the reduced tracer
diffusion coefficient $D_{0}(\alpha_{0})/D_{0}(1)$ as a function of the
intruder-grain coefficient of restitution $\alpha_{0}$ for the systems
($\mu=10$, $\omega=5$) (a) and ($\mu=0.5$, $\omega=0.5$) (b) in the case of a
three-dimensional gas with $\alpha=0.9$, $\phi=0.2$, and
$T^{*}_{\text{ex}}=1$. Solid and dashed lines are for the second and first
Sonine approximations, respectively, whereas the symbols refer to DSMC
results. Here, $D_{0}(1)$ is the value of the tracer diffusion coefficient
consistently obtained in each approximation when $\alpha_{0}=1$.

Finally, the accuracy of the first and second Sonine approximations to the
tracer diffusion coefficient $D_{0}$ is evaluated by performing DSMC
simulations. The ratio $D_{0}(\alpha_{0})/D_{0}(1)$ is plotted in Figs. 8 and
9 as a function of the intruder-grain coefficient of restitution $\alpha_{0}$
($\alpha\neq\alpha_{0}$). Here, $D_{0}(1)$ denotes the tracer diffusion
coefficient value when $\alpha_{0}=1$. Initially, it is observed that the
tracer diffusion coefficient $D_{0}$ exhibits minimal sensitivity to the
mechanical properties of the particles when the mass densities are equal.
However, in Fig. 9, significant qualitative differences in the behavior of
$D_{0}$ for the two mixtures are evident. It is also quite apparent that a
non-monotonic trend in the $\alpha_{0}$-dependence of $D_{0}$ is consistently
observed in both cases. Furthermore, a deliberate selection of a density value
of $\phi=0.2$ and a reduced background temperature of $T_{\text{ex}}^{*}=1$
allows for a larger discrepancy between both Sonine approximations, enabling a
clearer assessment of their accuracy through simulations. As seen in the dry
granular case when $\alpha=\alpha_{0}$, the convergence between both Sonine
approximations is practically achieved when the intruders possess both greater
mass and size compared to the grains. Although noticeable deviations are
observed as $\omega$ and $\mu$ decrease, these discrepancies are smaller than
those observed in the dry granular case. Garzó (2019) The distinctions between
the Sonine approximations are used to evaluate their accuracy in comparison to
simulations. Thus, Figs. 8 and 9 provide clear evidence that the second Sonine
approximation for the ratio $D_{0}(\alpha_{0})/D_{0}(1)$ improves the
theoretical predictions of the first Sonine approximation, yielding excellent
agreement with the simulation results.

##  VII Stability analysis of the homogeneous steady state

The knowledge of the diffusion transport coefficients along with the
Navier–Stokes transport coefficients of the granular suspension opens up the
possibility of solving the hydrodynamic equations for the densities $n_{0}$
and $n$, the flow velocity $\mathbf{U}$, and the granular temperature $T$ for
states near to the homogeneous steady state. This solution will give us
information on the linear stability of the above homogeneous state. The
linearized hydrodynamic equations are obtained when one substitutes the
constitutive equations for the mass flux, pressure tensor, and heat flux into
Eqs. (13)–(15) and (30). These equations are given by

$D_{t}n_{0}+n_{0}\nabla\cdot\mathbf{U}=\nabla\cdot\Bigg{(}\frac{m_{0}}{\rho}D_{0}\nabla
n_{0}+\frac{m_{0}}{\rho}D\nabla n+\frac{\rho}{m_{0}T}D^{T}\nabla T\Bigg{)},$
(89) $D_{t}n+n\nabla\cdot\mathbf{U}=0,$ (90)
$D_{t}U_{i}+\rho^{-1}\partial_{i}p=\rho^{-1}\partial_{j}\left[\eta\left(\partial_{i}U_{j}+\partial_{j}U_{i}-\frac{2}{d}\delta_{ij}\nabla\cdot\mathbf{U}\right)+\eta_{\text{b}}\delta_{ij}\nabla\cdot\mathbf{U}\right]-\gamma\Delta\mathbf{U},$
(91)
$\displaystyle\Big{(}D_{t}+2\gamma\left(1-\theta^{-1}\right)+\zeta^{(0)}\Big{)}T$
$\displaystyle=$ $\displaystyle\frac{2}{dn}\nabla\cdot\left(\kappa\nabla
T+\mu\nabla
n\right)+\frac{2}{dn}\bigg{[}\eta\left(\partial_{i}U_{j}+\partial_{j}U_{i}-\frac{2}{d}\delta_{ij}\nabla\cdot\mathbf{U}\right)$
(92)
$\displaystyle+\eta_{\text{b}}\delta_{ij}\nabla\cdot\mathbf{U}\bigg{]}\partial_{i}U_{j}-T\zeta_{U}\nabla\cdot\mathbf{U}-\frac{2}{dn}p\nabla\cdot\mathbf{U}.$

Here, $\eta$ and $\eta_{\text{b}}$ are the shear and bulk viscosities,
respectively, $\kappa$ is the thermal conductivity, $\mu$ is the diffusive
heat conductivity, and $\zeta_{U}$ is the first-order contribution to the
cooling rate. All these quantities have been determined in Ref. Gómez González
and Garzó, 2019 as functions of the coefficient of restitution $\alpha$, the
volume fraction $\phi$, and the background temperature $T_{\text{ex}}$.

In the homogeneous steady state, the hydrodynamic fields take the homogeneous
steady values $n_{0,\text{s}}\equiv\text{const.}$,
$n_{\text{s}}\equiv\text{const.}$, $T_{\text{s}}\equiv\text{const.}$, and
$\mathbf{U}_{g}=\mathbf{U}\equiv\mathbf{0}$. The subscript $s$ means that the
quantities are evaluated in the steady state. For small spatial gradients, we
assume that the deviations $\delta
y_{\beta}(\mathbf{r},t)=y_{\beta}(\mathbf{r},t)-y_{\beta,\text{s}}$ are small,
where $\delta y_{\beta}(\mathbf{r},t)$ denotes the deviations of the
hydrodynamic fields
$\left\\{y_{\beta};\beta=1,\cdots,d+3\right\\}=\left\\{n_{0},n,\mathbf{U},T\right\\}$
from their values in the homogeneous steady state. Moreover, as usual we also
suppose that the interstitial fluid is not perturbed and hence,
$\mathbf{U}_{g}=\mathbf{U}=\mathbf{0}$.

It is quite apparent that the hydrodynamic equations (90)–(92) for $n$,
$\mathbf{U}$, and $T$ are decoupled from Eq. (89) for $n_{0}$. This is a
consequence of the tracer limit. Equations (90)–(92) were already solved in
the Fourier space in Ref. Gómez González and Garzó, 2019 showing that this set
of equations is linearly stable. This means that the corresponding
perturbations ($\delta n$, $\delta\mathbf{U}$, $\delta T$) tend to zero for
sufficiently long times. In this case, Eq. (89) becomes an autonomous
differential equation for $\delta n_{0}$. To solve it, as in Ref. Gómez
González and Garzó, 2019, we introduce the following space and time variables:

$t^{\prime}=\frac{1}{2}n_{\text{s}}\sigma^{d-1}\sqrt{\frac{T_{\text{s}}}{m}}t,\quad\mathbf{r}^{\prime}=\frac{1}{2}n_{\text{s}}\sigma^{d-1}\mathbf{r}.$
(93)

The dimensionless time scale $t^{\prime}$ measures the average number of
collisions per particle in the time interval between 0 and $t$. The unit
length $\mathbf{r}^{\prime}$ is proportional to the mean free path of solid
particles. Moreover,

$\rho_{0,\mathbf{k}}(\tau)=\frac{\delta
n_{0,\mathbf{k}}(\tau)}{n_{0,\text{s}}}$ (94)

denotes the Fourier transformed dimensionless tracer density where

$\delta
n_{0,\mathbf{k}}(\tau)=\int\text{d}\mathbf{r}^{\prime}\text{e}^{-i\mathbf{k}\cdot\mathbf{r}^{\prime}}\delta
n_{0}(\mathbf{r}^{\prime},\tau).$ (95)

Note that here the wave vector $\mathbf{k}$ is dimensionless. In terms of the
above dimensionless variables, $\rho_{0,\mathbf{k}}(\tau)$ verifies the
equation

$\frac{\partial}{\partial
t^{\prime}}\rho_{0,\mathbf{k}}=-\frac{1}{2\sqrt{2}}\frac{m}{m_{0}}D_{0}^{*}k^{2}\rho_{0,\mathbf{k}},$
(96)

whose solution is

$\rho_{0,\mathbf{k}}(t^{\prime})=\rho_{0,\mathbf{k}}(0)\exp\left(-\frac{1}{2\sqrt{2}}\frac{m}{m_{0}}D_{0}^{*}k^{2}t^{\prime}\right).$
(97)

The first Sonine approximation to $D_{0}^{*}$ given by Eq. (72) clearly shows
that $D_{0}^{*}[1]$ is always positive. As expected, a careful analysis of the
dependence of the second Sonine approximation $D_{0}^{*}[2]$ on the parameter
space of the system also shows that $D_{0}^{*}[2]>0$. As a consequence,
$\rho_{0,\mathbf{k}}$ tends to zero for $t^{\prime}\to\infty$ and hence, the
homogeneous steady state for the system (tracer plus granular suspension) is
_linearly_ stable.

##  VIII Thermal diffusion segregation

Another nice and interesting application of the results derived in sections IV
and V is the study of thermal diffusion segregation of an intruder in a
granular suspension. The understanding of the physical mechanisms involved in
segregation and mixing of dissimilar grains is likely one of the most
challenging problems in granular mixtures. In addition, the problem is
relevant not only from a fundamental point of view but also from a more
practical perspective. Needless to say, segregation has been widely studied in
the past few years for dry granular mixtures. Among the different segregation
problems, the so-called Brazil-nut effect (BNE) is perhaps one of the most
famous examples of (size) segregation in vertically vibrated mixtures. In the
Brazil-nut effect, a relatively large particle (intruder) tends to climb to
the top of the sample against gravity. Rosato  _et al._ (1987); Duran,
Rajchenbach, and Clément (1993); Cooke  _et al._ (1996) On the other hand,
another different experimental works Shinbrot and Muzzio (1998); Hong, Quinn,
and Luding (2001) have observed the reverse buoyancy effect: the so-called
reverse Brazil-nut effect (RBNE) where the intruder may also sink to the
bottom of the sample under certain conditions. Several theoretical Arnarson
and Willits (1998); Jenkins and Yoon (2002); Trujillo, Alam, and Herrmann
(2003); Alam, Trujillo, and Herrmann (2006); Yoon and Jenkins (2006); Garzó
(2006, 2008); Garzó and Vega Reyes (2009); Garzó (2009) and computational
Galvin, Dahl, and Hrenya (2005); Brey, Ruiz-Montero, and Moreno (2005, 2006);
Brito  _et al._ (2008); Brito and Soto (2009) works have attempted to explain
the segregation induced by the presence of a thermal gradient (thermal
diffusion segregation) and/or the gravitational force. All these results refer
to situations where the effect of the interstitial gas (bath) on the dynamics
of tracer particles has been neglected. Our objective here is to assess the
influence of the background gas on the segregation criterion.

Thermal diffusion is the transport of matter caused by the presence of a
thermal gradient. As a result of the motion of the species of the mixture, a
steady state is finally achieved where the remixing effect of ordinary
diffusion is balanced by the segregation effect arising from thermal
diffusion. Kincaid, Cohen, and López de Haro (1987) The amount of segregation
parallel to the thermal gradient may be characterized by the so-called thermal
diffusion factor $\Lambda$. This quantity is defined in an inhomogeneous non-
convecting ($\mathbf{U}=\mathbf{U}_{g}=\mathbf{0}$) steady state with zero
mass flux ($\mathbf{j}_{0}^{(1)}=\mathbf{0}$) through the relation

$-\Lambda\frac{\partial\ln T}{\partial z}=\frac{\partial}{\partial
z}\ln\Big{(}\frac{n_{0}}{n}\Big{)},$ (98)

where only gradients along the $z$ axis have been assumed for the sake of
simplicity. Let us assume that the tracer is larger than the particles of the
granular gas ($\sigma_{0}>\sigma$) and also that gravity and the thermal
gradient point in parallel directions (i.e., the bottom plate is hotter than
the top plate, $\partial_{z}\ln T<0$, see Fig. 10). If the temperature
gradient is sufficiently small, one can assume that $\Lambda$ is practically
constant over the relevant ranges of composition and temperature of the
system. Thus, according to Eq. (98), when $\Lambda>0$ the intruders (or tracer
particles) tend to rise with respect to the particles of the granular gas
(i.e, $\partial_{z}\ln(n_{0}/n)>0$) and hence, the intruder particles tend to
accumulate near the cold plate (BNE). On the other hand, when $\Lambda<0$ the
intruders tend to fall with respect to the particles of the granular gas (i.e,
$\partial_{z}\ln(n_{0}/n)<0$) and so, the intruder particles tend to
accumulate near the hot plate (RBNE).

Figure 10: (color online) Illustration of the behavior of the gas-solid
system. The small circles in the diagram represent granular particles, while
the large circles depict intruders. All of them are immersed in a sea of
molecular particles (very small particles). Specifically, the BNE (RBNE)
effect occurs when the intruder rises (falls) to the top (bottom) plate of the
system.

Let us determine the thermal diffusion factor. The mass flux $j_{0,z}^{(1)}$
of tracer particles is given by Eq. (1) with $\Delta\mathbf{U}=\mathbf{0}$.
Since $j_{0,z}^{(1)}=0$ in the steady state and
$\mathbf{U}=\mathbf{U}_{g}=\mathbf{0}$, then the momentum balance equation
(14) yields $\partial_{z}p=-\rho g$. In dimensionless form, this relation can
be written as

$p^{*}+\left(p^{*}+\phi\frac{\partial
p^{*}}{\partial\phi}\right)\frac{\partial_{z}\ln n}{\partial_{z}\ln
T}=-g^{*},$ (99)

where $p^{*}=p/(nT)$ and $g^{*}=\rho g/n\partial_{z}T<0$ is a dimensionless
parameter measuring the gravity relative to the thermal gradient. The
condition $j_{0,z}^{(1)}=0$ along with Eq. (99) allows one to express the
factor $\Lambda$ in terms of the (dimensionless) diffusion coefficients
$D_{0}^{*}$, $D^{*}$, and $D^{T*}$ as

$\Lambda=\frac{\xi
D^{T^{*}}-\left(D_{0}^{*}+D^{*}\right)\left(g^{*}+p^{*}\right)}{\xi
D_{0}^{*}},$ (100)

where $\xi=p^{*}+\phi\partial_{\phi}p^{*}$. The dependence of $\Lambda$ on the
parameters of the system can be obtained when one substitutes the expressions
of the (reduced) pressure $p^{*}$ [Eq. (57)] and the (reduced) diffusion
transports coefficients into Eq. (100). Since we have seen that the first
Sonine approximation to these coefficients is quite accurate in the complete
range of values of the coefficients of restitution, we will consider
henceforth those expressions for the sake of simplicity. In this
approximation, the coefficients $D^{T*}$, $D_{0}^{*}$, and $D^{*}$ are given
by Eqs. (120), (132), and (134), respectively.

As expected, it is quite apparent from Eq. (100) that the influence of the
parameter space of the system on the sign of $\Lambda$ is quite complex, given
the large number of parameters involved. On the other hand, since Eq. (132)
clearly shows that $D_{0}^{*}>0$, then the condition $\Lambda=0$ yields the
relation

$\xi D^{T^{*}}=\left(D_{0}^{*}+D^{*}\right)\left(g^{*}+p^{*}\right).$ (101)

Equation (101) gives the curves delineating the regimes between the
segregation towards the cold and the hot wall (BNE/RBNE transition). However,
the explicit dependence of the segregation criterion (101) on the parameter
space is still cumbersome even if one employs the first Sonine solution to the
diffusion transport coefficients. For this reason, to disentangle the
different competing mechanisms appearing in the segregation problem, it is
first convenient to consider some simple situations where a more simple
criterion may be offered.

###  VIII.1 Mechanically equivalent particles

In this limit case ($m_{0}=m$, $\sigma_{0}=\sigma$, $\alpha=\alpha_{0}$,
$\gamma=\gamma_{0}$), $\tau_{0}=1$ and according to Eqs. (120), (132), and
(134), $D^{T^{*}}=0$ and $D_{0}^{*}=-D^{*}=(\nu_{1}^{*}+\gamma^{*})^{-1}$.
Thus, Eq. (101) holds for any value of the coefficients of restitution,
masses, diameters, and volume fraction. In this case, as expected, no
segregation is possible.

###  VIII.2 Low density limit

Let us consider a binary granular suspension in the low-density regime
($\phi\to 0$) and in the absence of gravity ($g^{*}=0$). In this regime,
$p^{*}=1$, $\xi=1$, and the diffusion coefficients are given by Eqs.
(75)–(77). According to these expressions, the segregation criterion (101)
yields

$\displaystyle\left(\tau_{0}-\mu\right)\left(\nu_{1}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}+\gamma_{0}^{*}\right)=\left(\nu_{1}^{*}+\gamma_{0}^{*}-\zeta^{*}\right)$
$\displaystyle\times\left(\tau_{0}+\theta\Delta_{\theta,0}-\mu\right).$ (102)

The segregation criterion (VIII.2) is still quite a complex relation in
comparison with the one obtained in the absence of gas phase (dry limit) where
$\Lambda=0$ when $\mu=\tau_{0}$. Brey, Ruiz-Montero, and Moreno (2005); Garzó
(2006) For elastic collisions, $\zeta^{*}=0$, $\theta=\tau_{0}=1$ and hence,
Eq. (VIII.2) becomes

$\left(1-\mu\right)\left(\nu_{1}^{*}-2\gamma^{*}+\gamma_{0}^{*}\right)=\left(\nu_{1}^{*}+\gamma_{0}^{*}\right)\Delta_{\theta,0}.$
(103)

As for inelastic collisions, Equation (103) clearly differs from the one
derived many years ago for molecular mixtures of hard spheres Kincaid, Cohen,
and López de Haro (1987) where the line separating the segregation towards the
cold or hot wall is simply $\mu=1$ in the first Sonine approximation.

In the case that the inhomogeneities in both the temperature and the mixture
volume fraction are neglected ($\partial_{z}T\to 0$) but gravity is different
from zero, then $|g^{*}|\to\infty$ (thermalized systems). In this case, the
segregation dynamics of the intruders is essentially driven by the
gravitational force. This situation (gravity dominates the temperature
gradient) can be achieved in the shaken or sheared systems employed in
numerical simulations and physical experiments. Hong, Quinn, and Luding
(2001); Wildman, Huntley, and Parker (2001); Breu  _et al._ (2003); Schautz
_et al._ (2005) Under these conditions ($|g^{*}|\to\infty$), the criterion
(VIII.2) leads to $D_{0}^{*}+D^{*}=0$. In the low-density limit, one gets

$\tau_{0}-\mu=\frac{\tau_{0}+\theta\Delta_{\theta,0}-\mu}{\nu_{1}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}+\gamma_{0}^{*}}\zeta^{*}.$
(104)

For elastic collisions, $\zeta^{*}=0$ and Eq. (104) reduces to $1-\mu=0$, as
expected.

###  VIII.3 Dense granular suspensions

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

Let us now study the segregation dynamics for moderately dense systems. We
analyze first situations where the influence of gravity can be neglected
($|g^{*}|\to 0$). In this particular scenario, the determination of the
segregation criterion is exclusively governed by the temperature gradient.

Let us consider first a three dimensional dry _molecular_ mixture
($\alpha=\alpha_{0}=1$) at moderate densities ($\phi=0.2$); this is
illustrated by the dashed line of Fig. 11. In this case, the particles of the
molecular fluid (“grains”) that are close to the hot plate move faster. When
tracer particles are introduced, they undergo substantial velocity changes
upon colliding with the grains as they approach to the hot plate. As the size
of the intruder particles increases (keeping the mass ratio $m_{0}/m$ fixed),
collisions between grains and intruders become more frequent due to the larger
effective surface area of the intruders. Consequently, these collisions cause
the intruders to rebound and migrate back to the colder side of the system
(RBNE). On the other hand, if we increase the mass ratio $m_{0}/m$ (keeping
the size ratio $\sigma_{0}/\sigma=$ fixed), the inertia of the massive
particle leads it to retain its initial trajectory upon collision.
Consequently, a greater number of collisions is required to modify the
trajectory of the intruder. The application of a thermal gradient implicitly
induces a gradient in particle density ($\partial_{z}n>0$) because the grains
experience a higher frequency of collisions in the hot plate, resulting in a
lower density compared to the cold plate. As a consequence of this density
gradient, collisions between the intruder and grains increase in the cold
plate, impeding the progress of the massive particle and causing it to migrate
toward the hotter zone (BNE).

When the mixture is immersed in an interstitial gas, we observe that the RBNE
zone is expanded, as depicted in the solid line of Fig. 11. This phenomenon
can be understood by considering the arguments previously discussed in Sec.
VI. First, the presence of the interstitial gas facilitates the thermalization
of the temperatures $T$ and $T_{0}$, leading to an equipartition of the
kinetic energy ($T=T_{0}=T_{\text{ex}}$). Thus, when grains are located in the
hotter (colder) region, the interstitial gas tends to decrease (increase) its
mean kinetic energy to establish thermal equilibrium. This effect is more
pronounced when larger particles are present, as their increased surface
amplifies the action of the interstitial gas
($\gamma_{0,\text{St}}\propto\sigma_{0}/\sigma$). When larger intruders are
located in the hotter region, the interstitial gas undergoes a faster
temperature decrease compared to that of the gas particles. This leads to an
increased disparity in the mean velocities between grains and intruders. This
discrepancy enables the grains to more effectively modify the trajectory of
the intruders, causing them to move back towards the colder zone. As a result,
the cooling effect caused by the presence of the bath combines with the size
effect observed in a dry gas and the region where the RBNE effect occurs
expands. Conversely, when the intruders are heavier than the grains, the
influence of the interstitial gas on intruders is weakened
($\gamma_{0,\text{St}}\propto(m_{0}/m)^{-1}$), thereby accentuating the
inertial effects observed in the dry case. This explains the reduction of the
BNE region for suspensions when smaller size ratios are considered.

The effects of the inelasticity in collisions on segregation are illustrated
in Fig. 12 in the case of the common coefficient of restitution
$\alpha=\alpha_{0}=0.8$. Notably, in contrast to the elastic case (depicted in
the dashed line of Fig. 11), we observe a pronounced enhancement of the RBNE
for dry granular mixtures; this enhancement is more significant within the
region where the intruders exhibit larger sizes. The inelastic nature of
collisions leads to a consequent loss of kinetic energy and a decrease in
temperature. Given the higher collision frequency in the hot plate, the
dissipative nature of these collisions becomes more pronounced. This
phenomenon bears similarity to the thermalizing influence exerted by the
interstitial gas on particles located in the hotter zone. Consequently, when
considering larger intruders, the BNE/RBNE diagram in a dry granular gas
exhibits a closer resemblance to the corresponding diagram observed in a
granular suspension. Interestingly, we find that the RBNE zone exhibits a
reduction in size as the degree of inelasticity increases in the case of
granular suspensions (solid line of Fig. 12). This is attributed to the bath’s
ability to lower the temperature of the intruder particles near the hot plate,
which in turn reduces their collision frequency. As a result, it becomes more
challenging for the tracer particles to undergo collisions that result in
rebound and subsequent return to the colder plate.

Figure 11: (color online) Plot of the marginal segregation curve ($\Lambda=0$)
for $d=3$, $\alpha=\alpha_{0}=1$, $\phi=0.2$, $T^{*}_{\text{ex}}=0.1$, and
$|g^{*}|\to 0$. The solid line corresponds to the segregation criterion for
molecular suspensions while the dashed line refers to the one derived for
molecular mixtures (no gas phase). The points below the curve correspond to
$\Lambda<0$ (RBNE), while the points above the curve correspond to $\Lambda>0$
(BNE). Figure 12: (color online) Plot of the marginal segregation curve
($\Lambda=0$) for $d=3$, $\alpha=\alpha_{0}=0.8$, $\phi=0.2$,
$T^{*}_{\text{ex}}=0.1$, and $|g^{*}|\to 0$. The solid line corresponds to the
segregation criterion for granular suspensions while the dashed line refers to
the one derived for dry granular mixtures. The points below the curve
correspond to $\Lambda<0$ (RBNE), while the points above the curve correspond
to $\Lambda>0$ (BNE).

####  VIII.3.2 Thermalized systems ($\partial_{z}T\to 0$)

We now consider the situation where the segregation dynamics is essentially
driven by the gravitational force. In this case, $|g^{*}|\to\infty$ and hence,
the temperature gradient can be neglected ($\partial_{z}T\to 0$). This
condition can be achieved in experimental setups and numerical simulations
involving shaken or sheared systems. Hong, Quinn, and Luding (2001); Breu  _et
al._ (2003); Schautz  _et al._ (2005); Qiao  _et al._ (2021); Shao  _et al._
(2023) Consequently, it is worthwhile to explore the influence of the
interstitial gas on the segregation criterion in this particular case. Figures
13 and 14 show the marginal segregation curve for the same systems as depicted
in Figs. 11 and 12. As apparent from Figs. 13 and 14, we observe that the
impact of the gravitational force on segregation is much more important than
that of the interstitial gas. In fact, both curves (with and without gas
phase) practically collapse in a common curve. Furthermore, similar to
granular mixtures driven by a stochastic external force (see Figs. 4 and 6 of
Ref. Garzó, 2009), we see that the presence of gravity completely reverses the
RBNE/BNE transition found in the phase diagrams in absence of gravity.
Regarding the effect of inelasticity in collisions, our results clearly
indicate that this effect is also negligible, even less significant than that
observed in a driven granular mixtures. Garzó (2009)

Therefore, to provide an explanation for BNE/RBNE transition, we focus on the
influence of gravity rather than dissipative dynamics or the force exerted by
the interstitial gas. Evidently, when the intruder particles are heavier, they
accumulate at the bottom of the container, leading to the observation of the
RBNE effect. This explains the left region of Figs. 13 and 14. Conversely, as
the size ratio $\sigma_{0}/\sigma$ increases, the tracer particles suffer more
collisions per time. The main consequence of these collisions is the existence
of a “buoyancy” effect on the intruder, resulting from the pressure exerted by
the sea of granular particles in the colder region. As a consequence, the
intruder is lifted against gravity, and the BNE effect is observed.

Figure 13: (color online) Plot of the marginal segregation curve ($\Lambda=0$)
for $d=3$, $\alpha=\alpha_{0}=1$, $\phi=0.2$, $T^{*}_{\text{ex}}=0.1$, and
$|g^{*}|\to\infty$. The solid line corresponds to the segregation criterion
for granular suspensions while the dashed line refers to the one derived for
dry granular mixtures. The points below the curve correspond to $\Lambda>0$
(BNE), while the points above the curve correspond to $\Lambda<0$ (RBNE).
Figure 14: (color online) Plot of the marginal segregation curve ($\Lambda=0$)
for $d=3$, $\alpha=\alpha_{0}=0.8$, $\phi=0.2$, $T^{*}_{\text{ex}}=0.1$, and
$|g^{*}|\to\infty$. The solid line corresponds to the segregation criterion
for granular suspensions while the dashed line refers to the one derived for
dry granular mixtures. The points below the curve correspond to $\Lambda>0$
(BNE), while the points above the curve correspond to $\Lambda<0$ (RBNE).

####  VIII.3.3 General case

Finally, we consider the general case where the effect of the temperature
gradient is comparable to that of gravity. To illustrate it, Fig. 15 shows the
marginal segregation curve for $|g^{*}|=1$, a (common) coefficient of
restitution $\alpha=\alpha_{0}=0.7$, $T_{\text{ex}}^{*}=0.1$, and two values
of the volume fraction $\phi$: $0.1$ and $0.2$. We observe that the effect of
gravity in the segregation dynamics overcomes that of the temperature
gradient, as seen from the reversion in the RBNE/BNE transition (compare Figs.
11 and 12 with Fig. 15). This finding is consistent with the previous results
obtained for dry granular mixtures.Garzó (2009) We also see that the effect of
the interstitial gas is to enhance the BNE effect, since the buoyancy
increases with the diameter ratios $\sigma_{0}/\sigma$. With regards to
density, an increase of shaking strength mimics the effect of decreasing the
volume fraction $\phi$ in the BNE/RBNE transition. Trujillo, Alam, and
Herrmann (2003); Alam, Trujillo, and Herrmann (2006) Figure 15 also
illustrates that the BNE effect is intensified as density increases (or
shaking strength is reduced). This behavior was observed in driven granular
mixtures Garzó (2009) and we see that is also applicable to both dry granular
mixtures and suspensions.

Figure 15: (color online) Plot of the marginal segregation curve ($\Lambda=0$)
for $d=3$, $\alpha=\alpha_{0}=0.7$, $T^{*}_{\text{ex}}=0.1$, $|g^{*}|=1$, and
two values of $\phi$. The solid line corresponds to the segregation criterion
for granular suspensions while the dashed line refers to the one derived for
dry granular mixtures. The points below the curve correspond to $\Lambda>0$
(BNE), while the points above the curve correspond to $\Lambda<0$ (RBNE).

##  IX Discussion

In this paper, the Enskog-Lorentz kinetic equation has been considered as the
starting point to obtain the diffusion transport coefficients of a binary
granular suspension where one of the species is present in tracer
concentration. As in previous works, Gómez González, Khalil, and Garzó (2020)
the influence of the interstitial gas on grains has been accounted for via a
drag force plus a stochastic Langevin term defined in terms of the (known)
background temperature $T_{\text{ex}}$. For arbitrary concentration, the
corresponding set of kinetic equations for the binary mixture were solved in a
previous work Gómez González, Khalil, and Garzó (2020) by means of the
Chapman–Enskog method Chapman and Cowling (1970) up to the first order in
spatial gradients. From this solution, the Navier–Stokes transport
coefficients of the mixture are given in terms of the solutions of a set of
coupled linear integral equations. In particular, in the case of the mass
flux, there are four relevant diffusion transport coefficients obeying
integral equations that are usually solved by considering only the leading
term in a Sonine polynomial expansion. Due to the technical difficulties
embodied in the general problem, we have considered here the tracer limit for
the sake of simplicity. In this limit, to solve the corresponding integral
equations, we have retained in the Sonine polynomial expansion terms up to the
second order and hence, we have explicitly determined the first and second
Sonine approximations to the diffusion coefficients as functions of the
parameter space of the system.

It should be noted that the accuracy of the Sonine approximation may be
questionable under the conditions of strong dissipation and/or quite disparate
size or mass ratios. However, our results have shown a relatively quite rapid
convergence between both Sonine approximations in comparison with the dry (no
gas phase) granular limit. In this context, we can conclude that the presence
of the interstitial gas tends to enhance the convergence of the Sonine
polynomial series expansion.

To assess the reliability of the Sonine approximations considered in this
paper, theoretical results have been compared against DSMC simulations. First,
a remarkable agreement is found between the theoretical predictions for the
temperature ratio $T_{0}/T$ (which has been obtained by approximating the
zeroth-order distributions by their Maxwellian forms) and the DSMC
simulations. Additionally, by employing the Einstein formula for the mean
square displacement of intruders moving in a granular suspension, the tracer
diffusion coefficient $D_{0}$ has been computationally determined. In order to
extend the study recently done in Ref. Gómez González  _et al._ , 2023 for the
evaluation of $D_{0}$ when $\alpha=\alpha_{0}$, we have performed in this
paper simulations when both coefficients of restitution are different
($\alpha\neq\alpha_{0}$). It is confirmed first that the (reduced) coefficient
$D_{0}(\alpha_{0})/D_{0}(1)$ [$D_{0}(1)$ being the tracer diffusion
coefficient for elastic collisions] exhibits a non-monotonic dependence on
$\alpha_{0}$. Furthermore, intentional selection of mixture parameters reveals
stronger disparities between the results obtained from the second and first
Sonine approximations for significantly small mass and size ratios.
Nevertheless, these discrepancies remain smoother than those have been
observed in the dry granular case. Garzó and Vega Reyes (2009, 2012) As
expected from previous works, Gómez González  _et al._ (2023) although the
theoretical predictions of the first Sonine approximation agree relatively
well with simulations, we find that the results of the second Sonine
approximation agree perfectly with the DSMC simulations. Thus, the second
Sonine approximation corrects the small discrepancies observed between the
first Sonine solution and simulations, indicating an improved convergence of
the solution.

The knowledge of the diffusion coefficients allow us to perform a linear
stability analysis of the homogeneous steady state of the suspension. This
analysis is relatively simple since the hydrodynamic equations for $n$,
$\mathbf{U}$, and $T$ are decoupled from that for $n_{0}$. Moreover, given
that a previous Gómez González and Garzó (2019) linear stability analysis had
shown that the hydrodynamic equations obeying the fields $n$, $\mathbf{U}$,
and $T$ are linearly stable, the hydrodynamic equation for $n_{0}$ becomes an
autonomous differential equation. We can at this point obtain the time
evolution of the Fourier transform to the (dimensionless) tracer density.
Since both the first ($D^{*}_{0}[1]$) and second ($D^{*}_{0}[2]$) (reduced)
Sonine approximations are always positive, any perturbation in the tracer
density tends to vanish over sufficiently long times. As a consequence, the
homogeneous steady state of the system (tracer plus granular suspension) has
been shown to be _linearly_ stable.

As a second application, segregation induced by a thermal gradient and/or
gravity has been also studied. The main target has been to assess the impact
of the interstitial gas on the phase diagrams previously obtained in the
absence of the gas phase. Garzó (2011) In the absence of gravity, our results
(see Figs. 11 and 12) show that the effect of the gas phase on segregation is
in general significant for large mass and/or diameter ratios. The main effect
of the surrounding gas is to increase the size of the RBNE region (intruders
attempt to accumulate near the hot plate) with respect to the one observed in
the absence of the gas phase. On the other hand, in the situations where
segregation is mainly driven by gravity (namely, when thermal gradient can be
neglected) the phase diagrams with and without gas phase are practically
identical (see Figs. 13 and 14). In addition, as for dry granular mixtures,
the influence of inelasticity of collisions on segregation is very tiny since
the corresponding marginal segregation curve for elastic collisions collapses
practically with the one obtained for inelastic collisions.

The results reported in this paper have been derived from a suspension model
where the effect of the background gas on the dynamics of grains has been
accounted for via a fluid-solid force (coarse-grained description). It would
be interesting to revisit the tracer diffusion problem by considering a
collisional model that explicitly takes into account the collisions between
grains and particles of the surrounding molecular gas. This sort of suspension
model has been recently proposed Gómez González and Garzó (2022b) and the
results derived from this collisional model has been shown to reduce to those
obtained from the coarse-grained approach Gómez González and Garzó (2019) when
the grains are much heavier than the particles of the background gas. We
expect a similar conclusion in the tracer diffusion problem analyzed in the
present paper. Last but not least, it would also be desirable to carry on
computer simulations to assess the reliability of the phase diagrams obtained
here for thermal diffusion segregation of intruders. We plan to make some
progresses along these lines in the near future.

###### Acknowledgements.

We acknowledge financial support from Grant PID2020-112936GB-I00 funded by
MCIN/AEI/10.13039/501100011033, and from Grant IB20079 funded by Junta de
Extremadura (Spain) and by ERDF A way of making Europe. The research of R.G.G.
also has been supported by Plan Propio de Iniciación a la Investigación,
Desarrollo Tecnológico e Innovación de la Universidad de Extremadura (ACCIÓN
III).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Rubén Gómez González: Formal analysis (equal); Investigation (equal); Software
(equal); Writing – review & editing (equal).
Vicente Garzó: Conceptualization (equal); Investigation (equal); Supervision
(equal); Writing – original draft (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A The second Sonine coefficients $a$ and $c$ of the granular gas

The first-order velocity distribution function $f^{(1)}$ of the granular gas
has the form Gómez González and Garzó (2019)

$f^{(1)}=\bm{\mathcal{A}}\cdot\nabla T+\bm{\mathcal{C}}\cdot\nabla
n+\bm{\mathcal{D}}:\nabla\mathbf{U}+E\nabla\cdot\mathbf{U}.$ (105)

To evaluate the diffusion transport coefficients of the tracer particles, only
the quantities $\bm{\mathcal{A}}$ and $\bm{\mathcal{C}}$ are needed. Both
quantities vanish in the first Sonine approximation. In the second Sonine
approximation, $\bm{\mathcal{A}}$ and $\bm{\mathcal{C}}$ are given,
respectively, by

$\bm{\mathcal{A}}({\bf V})\to-f_{\text{M}}({\bf V})a\;{\bf S}({\bf
V}),\quad\bm{\mathcal{C}}({\bf V})\to-f_{\text{M}}({\bf V})c\;{\bf S}({\bf
V}),$ (106)

where the Maxwellian distribution $f_{\text{M}}({\bf V})$ is defined by Eq.
(18) with the replacement $\mathbf{v}\to\mathbf{V}$ and
$\mathbf{S}(\mathbf{V})$ is

${\bf S}({\bf V})=\left(\frac{1}{2}mV^{2}-\frac{d+2}{2}T\right){\bf V}.$ (107)

The coefficients $a$ and $c$ are related to the kinetic contributions to the
thermal conductivity $\kappa_{k}$ and the diffusive heat conductivity
$\mu_{k}$. They are defined as

$\left\\{a,c\right\\}=\frac{2}{d+2}\frac{m}{nT^{3}}\left\\{\kappa_{\text{k}},\mu_{\text{k}}\right\\}.$
(108)

Up to the second Sonine approximation, the coefficients $a$ and $c$ are Gómez
González and Garzó (2019)

$a=\frac{1}{\nu
T^{2}}\frac{1+3\frac{2^{d-3}}{d+2}\chi^{(0)}\phi\left(1+\alpha\right)^{2}\left(2\alpha-1\right)}{\nu_{\kappa}^{*}+\gamma^{*}-\frac{3}{2}\zeta^{*}},$
(109) $\displaystyle c$ $\displaystyle=$ $\displaystyle\frac{1}{\nu
nT}\left(\nu_{\kappa}^{*}+3\gamma^{*}\right)^{-1}\Bigg{\\{}a^{*}\Bigg{[}\zeta^{*}\left(1+\phi\frac{\partial\ln\chi^{(0)}}{\partial\phi}\right)$
(110)
$\displaystyle-2\left(\theta^{-1}-1\right)\gamma^{*}\phi\frac{\partial\ln
R}{\partial\phi}\Bigg{]}+3\frac{2^{d-2}}{\left(d+2\right)}\chi^{(0)}\phi\left(1+\alpha\right)$
$\displaystyle\times\left(1+\frac{1}{2}\phi\frac{\partial\ln\chi^{(0)}}{\partial\phi}\right)\alpha(\alpha-1)\Bigg{\\}}.$

In Eqs. (109) and (110), $\chi^{(0)}$ is the pair correlation at zeroth order,
$\gamma^{*}=\gamma/\nu=\lambda\theta^{-1/2}$, $a^{*}=\nu T^{2}a$, and

$\displaystyle\nu_{\kappa}^{*}$ $\displaystyle=$
$\displaystyle\frac{8}{\sqrt{2}(d+2)}\frac{\pi^{(d-1)/2}}{\Gamma\left(\frac{d}{2}\right)}\frac{1+\alpha}{d}\chi^{(0)}$
(111)
$\displaystyle\times\left[\frac{d-1}{2}+\frac{3}{16}(d+8)(1-\alpha)\right].$

It must be recalled again that Eqs. (109)–(111) have been displayed neglecting
the non-Gaussian corrections to the zeroth-order distribution function
$f^{(0)}$. As said along the paper, the effect of these contributions on the
Sonine coefficients $a$ and $c$ are significant only for quite extreme values
of inelasticity. Gómez González and Garzó (2019)

##  Appendix B Sonine approximations to the diffusion transport coefficients

In this appendix we determine the first and second Sonine approximations to
the diffusion transport coefficients $D_{0}$, $D$, $D^{T}$, and $D^{U}$. These
coefficients are given in terms of the solutions of the set of integral
equations (48)–(51), where the inhomogeneous terms $\mathbf{A}_{0}$,
$\mathbf{B_{0}}$, $\mathbf{C}_{0}$, and $\mathbf{E}_{0}$ are defined by Eqs.
(53)–(56), respectively. Up to the second Sonine approximation, the functions
$\bm{\mathcal{A}}_{0}$, $\bm{\mathcal{B}}_{0}$, $\bm{\mathcal{C}}_{0}$, and
$\bm{\mathcal{E}}_{0}$ are approximated by Eqs. (60)–(63), respectively, while
$\bm{\mathcal{A}}$ and $\bm{\mathcal{C}}$ are defined by Eq. (106). Let us
evaluate each transport coefficient separately.

###  B.1 Thermal diffusion coefficient $D^{T}$

We consider first the thermal diffusion coefficient $D^{T}$. To get it, we
substitute $\bm{\mathcal{A}}$ and $\bm{\mathcal{A}}_{0}$ by their Sonine
approximations and then, multiply the integral equation (48) by
$m_{0}\mathbf{V}$ and integrate over velocity. After some algebra, one gets

$\Bigg{(}\nu_{1}-2\gamma\theta^{-1}-\frac{1}{2}\zeta^{(0)}+\gamma_{0}\Bigg{)}D^{T}+\frac{n_{0}T_{0}^{(0)2}}{\rho}\nu_{2}a_{0}=Z_{1},$
(112)

where

$Z_{1}=-\frac{n_{0}T_{0}^{(0)2}}{\rho}\nu_{3}a+\frac{n_{0}T_{0}^{(0)}}{\rho}+\frac{n_{0}T}{\rho}\theta\Delta_{\theta,0}-\frac{p\rho_{0}}{\rho^{2}}+\frac{m_{0}}{d\rho}\int
d\mathbf{V}\;\mathbf{V}\cdot\bm{\mathcal{K}}_{0}\Bigg{[}T\frac{\partial
f^{(0)}}{\partial T}\Bigg{]}.$ (113)

Here, $\Delta_{\theta,0}=\partial\tau_{0}/\partial\theta$,
$\tau_{0}=T_{0}^{(0)}/T$, and we have introduced the collision frequencies

$\nu_{1}=-\frac{1}{dn_{0}T_{0}^{(0)}}\int d{\bf v}\,m_{0}{\bf V}\cdot
J_{0}^{(0)}[f_{0,M}{\bf V},f^{(0)}],$ (114)
$\nu_{2}=-\frac{1}{dn_{0}T_{0}^{(0)2}}\int d{\bf v}\,m_{0}{\bf V}\cdot
J_{0}^{(0)}[f_{0,M}{\bf S}_{0},f^{(0)}],$ (115)
$\nu_{3}=-\frac{1}{dn_{0}T_{0}^{(0)2}}\int d{\bf v}\,m_{0}{\bf V}\cdot
J_{0}^{(0)}[f_{0}^{(0)},f_{M}{\bf S}].$ (116)

In Eqs. (112) and (113), $\zeta^{(0)}$ and $T_{0}^{(0)}$ are the zeroth-order
contributions to the cooling rate $\zeta$ and the partial temperature $T_{0}$,
respectively. In the Maxwellian approximation, the cooling rate $\zeta^{(0)}$
is defined by Eq. (19) and $T_{0}^{(0)}$ is given by the physical solution to
Eq. (42). Upon obtaining (113), use has been made of the result

$\int
d\mathbf{V}m_{0}\mathbf{V}\cdot\mathbf{A}_{0}=-dn_{0}T_{0}^{(0)}-dn_{0}T\theta\Delta_{\theta,0}+\frac{dp}{\rho}\rho_{0}-\int
d{\bf v}m_{0}\mathbf{V}\cdot\bm{\mathcal{K}}_{0}\Bigg{[}T\frac{\partial
f^{(0)}}{\partial T}\Bigg{]}.$ (117)

The collision integral of Eq. (117) involving the operator
$\bm{\mathcal{K}}_{0}$ can be computed when one replaces the distribution
function $f^{(0)}$ by its Maxwellian form (39). In this approximation, one has
the identity

$T\frac{\partial f^{(0)}}{\partial
T}=-\frac{1}{2}\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f^{(0)}$
(118)

and hence, one obtains the result Garzó, Hrenya, and Dufty (2007)

$\int d{\bf v}m_{0}\mathbf{V}\cdot\bm{\mathcal{K}}_{0}\Bigg{[}T\frac{\partial
f^{(0)}}{\partial
T}\Bigg{]}=2^{d-1}dn_{0}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d}\chi_{0}^{(0)}\phi
M_{0}(1+\alpha_{0}).$ (119)

If only the first Sonine correction is retained in Eq. (112) (which means
$a_{0}=a=0$), then the (reduced) first Sonine approximation
$D^{*T}[1]=(m\nu/x_{0}T)D^{T}[1]$ to the thermal diffusion coefficient is

$D^{*T}[1]=-\Bigg{(}\nu_{1}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}+\gamma_{0}^{*}\Bigg{)}^{-1}\Bigg{[}\frac{m_{0}}{m}p^{*}-\tau_{0}-\theta\Delta_{\theta,0}-2^{d-1}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d}\chi_{0}^{(0)}M_{0}\phi(1+\alpha_{0})\Bigg{]},$
(120)

where $\nu_{i}^{*}=\nu_{i}/\nu$ and $p^{*}=p/(nT)$ is given by Eq. (57).

To close the evaluation of $D^{T}$ up to the second Sonine approximation, we
multiply now Eq. (48) by $\mathbf{S}_{0}(\mathbf{V})$ and integrates over
$\mathbf{v}$. After a tedious algebra, one achieves the result

$\Bigg{(}\nu_{4}+3\gamma_{0}-2\gamma\theta^{-1}-\frac{1}{2}\zeta^{(0)}\Bigg{)}a_{0}+\frac{\rho}{n_{0}T_{0}^{(0)2}}\Bigg{[}\nu_{5}+2\gamma_{0}\Bigg{(}1-\frac{1}{d+2}\frac{T_{\text{ex}}}{T_{0}^{(0)}}\Bigg{)}\Bigg{]}D^{T}=Z_{2},$
(121)

where

$Z_{2}=-\nu_{6}a+\frac{1}{T_{0}^{(0)}}+\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{(0)3}}\int
d\mathbf{V}\;\mathbf{S}_{0}\cdot\bm{\mathcal{K}}_{0}\Bigg{[}T\frac{\partial
f^{(0)}}{\partial T}\Bigg{]},$ (122)
$\nu_{4}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{(0)3}}\int d{\bf v}\,{\bf
S}_{0}\cdot J_{0}^{(0)}[f_{0,M}{\bf S}_{0},f^{(0)}],$ (123)
$\nu_{5}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{(0)2}}\int d{\bf v}\,{\bf
S}_{0}\cdot J_{0}^{(0)}[f_{0,M}{\bf V},f^{(0)}],$ (124)
$\nu_{6}=-\frac{2}{d(d+2)}\frac{m_{0}}{n_{0}T_{0}^{(0)3}}\int d{\bf v}\,{\bf
S}_{0}\cdot J_{0}^{(0)}[f_{0}^{(0)},f_{M}{\bf S}].$ (125)

As before, if one replaces $f_{0}^{(0)}$ by $f_{0,\text{M}}$, the collision
integral involving the operator $\bm{\mathcal{K}}_{0}$ in Eq. (122) is given
by Garzó, Hrenya, and Dufty (2007)

$\displaystyle\int d{\bf v}\mathbf{S}_{0}({\bf V})\cdot\bm{{\cal
K}}_{0}\left[T\frac{\partial f^{(0)})}{\partial T}\right]$ $\displaystyle=$
$\displaystyle\frac{d}{4}n_{0}\frac{M_{0}T^{2}}{m}(1+\omega)^{d}\chi_{0}^{(0)}\phi(1+\alpha_{0})\Bigg{\\{}\frac{M\tau_{0}}{M_{0}}\Big{[}(d+2)(M_{0}^{2}-1)$
(126)
$\displaystyle+(2d-5-9\alpha_{0})M_{0}M+(d-1+3\alpha_{0}+6\alpha_{0}^{2})M^{2}\Big{]}+6M^{2}(1+\alpha_{0})^{2}\Bigg{\\}}.$

The (reduced) thermal diffusion coefficient $D^{T*}$ is defined in Eq. (69).
Thus, in reduced units and by using matrix notation, Eqs. (112) and (121) can
be rewritten as

$\left(\begin{array}[]{cc}\nu_{1}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}+\gamma_{0}^{*}&\tau_{0}^{2}\nu_{2}^{*}\\\
\frac{\nu_{5}^{*}+2\gamma_{0}^{*}\Big{(}1-(\theta\tau_{0})^{-1}\Big{)}}{\tau_{0}^{2}}&\nu_{4}^{*}+3\gamma_{0}^{*}-2\gamma^{*}\theta^{-1}-\frac{1}{2}\zeta^{*}\end{array}\right)\left(\begin{array}[]{c}D^{T*}\\\
a_{0}^{*}\end{array}\right)=\left(\begin{array}[]{c}X_{1}^{*}-\tau_{0}^{2}\nu_{3}^{*}a^{*}\\\
X_{2}^{*}-\nu_{6}^{*}a^{*}\end{array}\right).$ (127)

Here, $a_{0}^{*}=T^{2}\nu a_{0}$, $a^{*}=T^{2}\nu a$, and

$X_{1}^{*}=\tau_{0}+\theta\Delta_{\theta,0}-\frac{pm_{0}}{\rho
T}+\frac{1}{2}(1+\omega)^{d}M_{0}\chi_{0}^{(0)}\phi(1+\alpha_{0}),$ (128)
$X_{2}^{*}=\frac{1}{\tau_{0}}+\frac{2}{d(d+2)}\frac{m_{0}T}{n_{0}T_{0}^{(0)3}}\int
d\mathbf{V}\;\mathbf{S}_{0}\cdot\bm{\mathcal{K}}_{0}\Bigg{[}T\frac{\partial
f^{(0)}}{\partial T}\Bigg{]}.$ (129)

The second Sonine approximation $D^{T*}[2]$ to $D^{T*}$ can be obtained from
Eq. (127) and the result is

$D^{T*}[2]=\frac{(3\gamma_{0}^{*}+\nu_{4}^{*}-\frac{1}{2}\zeta^{*}-2\gamma^{*}\theta^{-1})(a^{*}\nu_{3}\tau_{0}^{2}-X_{1}^{*})-\nu_{2}^{*}\tau_{0}^{2}(a^{*}\nu_{6}^{*}-X_{2}^{*})}{\nu_{2}^{*}(\nu_{5}^{*}+2\gamma_{0}^{*}(1-(\tau_{0}\theta)^{-1})-(\gamma_{0}^{*}+\nu_{1}^{*}-\frac{1}{2}\zeta^{*}-2\gamma^{*}\theta^{-1})(3\gamma_{0}^{*}+\nu_{4}^{*}-\frac{1}{2}\zeta^{*}-2\gamma^{*}\theta^{-1})}.$
(130)

To obtain the dependence of $D^{T*}[1]$ and $D^{T*}[2]$ on the parameter space
of the system, it still remains to determine the (reduced) collision
frequencies $\nu_{i}^{*}$. Explicit expressions of them have been achieved in
previous works Garzó, Hrenya, and Dufty (2007); Garzó and Vega Reyes (2012);
Garzó (2019) when the distributions $f^{(0)}$ and $f_{0}^{(0)}$ are replaced
by their Maxwellian approximations $f_{M}$ and $f_{0,M}$, respectively. These
expressions will be provided in the Appendix C for the sake of completeness.

The determination of the first and second Sonine approximations of the
diffusion coefficients $D_{0}$ and $D$ follows similar steps as those made for
$D^{T}$. Below, only the final expressions will be provided.

###  B.2 Tracer diffusion coefficient $D_{0}$

The (reduced) tracer diffusion coefficient $D_{0}^{*}$ is defined by Eq. (69).
It obeys the matrix equation

$\left(\begin{array}[]{cc}\nu_{1}^{*}+\gamma_{0}^{*}&\tau_{0}^{2}\nu_{2}^{*}\\\
\frac{\nu_{5}^{*}+2\gamma_{0}^{*}\Big{(}1-(\theta\tau_{0})^{-1}\Big{)}}{\tau_{0}^{2}}&\nu_{4}^{*}+3\gamma_{0}^{*}\end{array}\right)\left(\begin{array}[]{c}D_{0}^{*}\\\
b_{0}^{*}\end{array}\right)=\left(\begin{array}[]{c}\tau_{0}\\\
0\end{array}\right),$ (131)

where $b_{0}^{*}=T\nu b_{0}$. Thus, the first $D_{0}^{*}[1]$ and second
$D_{0}^{*}[2]$ Sonine approximations to $D_{0}^{*}$ read, respectively, as

$D_{0}^{*}[1]=\frac{\tau_{0}}{\nu_{1}^{*}+\gamma_{0}^{*}},$ (132)
$D_{0}^{*}[2]=\frac{(3\gamma_{0}^{*}+\nu_{4}^{*})\tau_{0}}{(\gamma_{0}^{*}+\nu_{1}^{*})(3\gamma_{0}^{*}+\nu_{4}^{*})-\nu_{2}^{*}(\nu_{5}^{*}+2\gamma_{0}^{*}(1-(\tau_{0}\theta)^{-1}))}.$
(133)

###  B.3 Mutual diffusion coefficient $D$

The (reduced) diffusion coefficient $D^{*}$ is defined in Eq. (70). After some
calculations, the first $D^{*}[1]$ and second $D^{*}[2]$ Sonine approximations
to $D^{*}$ can be written, respectively, as

$\displaystyle D^{*}[1]$ $\displaystyle=$
$\displaystyle(\nu_{1}^{*}+\gamma_{0}^{*})^{-1}\Bigg{\\{}\Bigg{[}\zeta^{*}\left(1+\phi\frac{\partial\ln\chi^{(0)}}{\partial\phi}\right)-2\left(\theta^{-1}-1\right)\gamma^{*}\phi\frac{\partial\ln
R}{\partial\phi}\Bigg{]}D^{T*}[1]+\phi\Delta_{\phi,0}-\frac{m_{0}}{m}\left(p^{*}+\phi\frac{\partial
p^{*}}{\partial\phi}\right)$ (134)
$\displaystyle+\frac{1}{2}M_{0}\phi(1+\alpha_{0})\left(\frac{1+\beta}{\beta}\right)\frac{\partial}{\partial\phi}\left(\frac{\mu_{0}}{T}\right)_{T,n_{0}}\Bigg{\\}},$
$D^{*}[2]=\frac{(3\gamma_{0}^{*}+\nu_{4}^{*})(c^{*}\nu_{3}^{*}\tau_{0}^{2}-Y_{1}^{*})-\nu^{*}_{2}\tau_{0}^{2}(c^{*}\nu_{6}^{*}-Y_{2}^{*})}{\nu_{2}^{*}(\nu_{5}^{*}+2\gamma_{0}^{*}(1-(\tau_{0}\theta)^{-1}))-(\gamma_{0}^{*}+\nu_{1}^{*})(3\gamma_{0}^{*}+\nu_{4}^{*})}.$
(135)

In Eq. (135), the dimensionless quantities $Y_{1}^{*}$ and $Y_{2}^{*}$ are
given by

$\displaystyle Y_{1}^{*}$ $\displaystyle=$
$\displaystyle\Bigg{[}\zeta^{*}\left(1+\phi\frac{\partial\ln\chi^{(0)}}{\partial\phi}\right)-2\left(\theta^{-1}-1\right)\gamma^{*}\phi\frac{\partial\ln
R}{\partial\phi}\Bigg{]}D^{T*}[2]$
$\displaystyle+\phi\Delta_{\phi,0}-\frac{m_{0}}{m}\left(p^{*}+\phi\frac{\partial
p^{*}}{\partial\phi}\right)+\frac{1}{2}M_{0}\phi(1+\alpha_{0})\left(\frac{1+\beta}{\beta}\right)\frac{\partial}{\partial\phi}\left(\frac{\mu_{0}}{T}\right)_{T,n_{0}},$
$\displaystyle Y_{2}^{*}$ $\displaystyle=$
$\displaystyle\Bigg{[}\zeta^{*}\left(1+\phi\frac{\partial\ln\chi^{(0)}}{\partial\phi}\right)-2\left(\theta^{-1}-1\right)\gamma^{*}\phi\frac{\partial\ln
R}{\partial\phi}\Bigg{]}a_{0}^{*}+\tau_{0}^{-2}\phi\Delta_{\phi,0}+\frac{1}{2(d+2)}\frac{M^{2}}{M_{0}}\phi(1+\alpha_{0})\frac{\partial}{\partial\phi}\left(\frac{\mu_{0}}{T}\right)_{T,n_{0}}$
(137)
$\displaystyle\times\Bigg{\\{}\Big{[}(d+8)M_{0}^{2}+(7+2d-9\alpha_{0})M_{0}M+(2+d+3\alpha_{0}^{2}-3\alpha_{0})M^{2}\Big{]}\beta$
$\displaystyle+3M^{2}(1+\alpha_{0})^{2}\beta^{3}+\Big{[}(d+2)M_{0}^{2}+(2d-5-9\alpha_{0})M_{0}M+(d-1+3\alpha_{0}+6\alpha_{0}^{2})M^{2}\Big{]}\beta^{2}$
$\displaystyle-(d+2)\beta(1+\beta)\Bigg{\\}},$

where

$a_{0}^{*}=\frac{\tau_{0}^{-2}\left[\nu_{5}^{*}+2\gamma_{0}^{*}(1-(\tau_{0}\theta)^{-1})\right]\left(X_{1}^{*}-a^{*}\nu_{3}\tau_{0}^{2}\right)-\left(\gamma_{0}^{*}+\nu_{1}^{*}-\frac{\zeta^{*}}{2}-2\gamma^{*}\theta^{-1}\right)\left(X_{2}^{*}-a^{*}\nu_{6}^{*}\right)}{\nu_{2}^{*}\left[\nu_{5}^{*}+2\gamma_{0}^{*}(1-(\tau_{0}\theta)^{-1})\right]-(\gamma_{0}^{*}+\nu_{1}^{*}-\frac{1}{2}\zeta^{*}-2\gamma^{*}\theta^{-1})(3\gamma_{0}^{*}+\nu_{4}^{*}-\frac{1}{2}\zeta^{*}-2\gamma^{*}\theta^{-1})}.$
(138)

###  B.4 Velocity diffusion coefficient $D^{U}$

The (reduced) velocity diffusion coefficient $D_{U}^{*}$ is defined by Eq.
(70). It obeys the matrix equation

$\left(\begin{array}[]{cc}\nu_{1}^{*}+\gamma_{0}^{*}&\tau_{0}^{2}\nu_{2}^{*}\\\
\frac{\nu_{5}^{*}+2\gamma_{0}^{*}\Big{(}1-(\theta\tau_{0})^{-1}\Big{)}}{\tau_{0}^{2}}&\nu_{4}^{*}+3\gamma_{0}^{*}\end{array}\right)\left(\begin{array}[]{c}D_{U}^{*}\\\
e_{0}^{*}\end{array}\right)=\left(\begin{array}[]{c}\gamma_{0}^{*}-\gamma^{*}\\\
0\end{array}\right),$ (139)

where $e_{0}^{*}=(T^{2}/m_{0})e_{0}$. The first $D_{U}^{*}[1]$ and second
$D_{U}^{*}[2]$ Sonine approximations to $D_{U}^{*}$ are given, respectively,
by

$D_{U}^{*}[1]=\frac{\gamma_{0}^{*}-\gamma^{*}}{\nu_{1}^{*}+\gamma_{0}^{*}},$
(140)
$D_{U}^{*}[2]=\frac{(3\gamma_{0}^{*}+\nu_{4}^{*})(\gamma^{*}-\gamma_{0}^{*})}{\nu_{2}^{*}(\nu_{5}^{*}+2\gamma_{0}^{*}(1-(\tau_{0}\theta)^{-1}))-(\gamma_{0}^{*}+\nu_{1}^{*})(3\gamma_{0}^{*}+\nu_{4}^{*})}.$
(141)

##  Appendix C Collision frequencies

In this appendix we display the expressions of the collision frequencies
$\nu_{1}^{*}$–$\nu_{6}^{*}$. They are given by Garzó, Hrenya, and Dufty
(2007); Garzó and Vega Reyes (2012)

$\nu_{1}^{*}=\frac{2\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\chi_{0}^{(0)}M(1+\alpha_{0})\left(\frac{1+\beta}{\beta}\right)^{1/2},$
(142)
$\nu_{2}^{*}=\frac{\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\chi_{0}^{(0)}M(1+\alpha_{0})[\beta(1+\beta)]^{-1/2},$
(143)
$\nu_{3}^{*}=-\frac{\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\chi_{0}^{(0)}\frac{M^{2}}{M_{0}}(1+\alpha_{0})\beta^{5/2}(1+\beta)^{-1/2},$
(144)
$\nu_{4}^{*}=\frac{\pi^{(d-1)/2}}{d(d+2)\Gamma\left(\frac{d}{2}\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\chi_{0}^{(0)}M(1+\alpha_{0})\left(\frac{\beta}{1+\beta}\right)^{3/2}\left[A-(d+2)\frac{1+\beta}{\beta}B\right],$
(145)
$\nu_{5}^{*}=\frac{2\pi^{(d-1)/2}}{d(d+2)\Gamma\left(\frac{d}{2}\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\chi_{0}^{(0)}M(1+\alpha_{0})\left(\frac{\beta}{1+\beta}\right)^{1/2}B,$
(146)
$\nu_{6}^{*}=-\frac{\pi^{(d-1)/2}}{d(d+2)\Gamma\left(\frac{d}{2}\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\chi_{0}^{(0)}\frac{M^{2}}{M_{0}}(1+\alpha_{0})\left(\frac{\beta}{1+\beta}\right)^{3/2}\left[C+(d+2)(1+\beta)D\right].$
(147)

In Eqs. (142)–(147), we have introduced the quantities

$\displaystyle A$ $\displaystyle=$ $\displaystyle
2M^{2}\left(\frac{1+\beta}{\beta}\right)^{2}\left(2\alpha_{0}^{2}-\frac{d+3}{2}\alpha_{0}+d+1\right)\big{[}d+5+(d+2)\beta\big{]}-M(1+\beta)\Big{\\{}\xi\beta^{-2}[(d+5)+(d+2)\beta]$
$\displaystyle\times[(11+d)\alpha_{0}-5d-7]-\beta^{-1}[20+d(15-7\alpha_{0})+d^{2}(1-\alpha_{0})-28\alpha_{0}]-(d+2)^{2}(1-\alpha_{0})\Big{\\}}$
$\displaystyle+3(d+3)\xi^{2}\beta^{-2}[d+5+(d+2)\beta]+2\xi\beta^{-1}[24+11d+d^{2}+(d+2)^{2}\beta]$
$\displaystyle+(d+2)\beta^{-1}[d+3+(d+8)\beta]-(d+2)(1+\beta)\beta^{-2}[d+3+(d+2)\beta],$
$\displaystyle B$ $\displaystyle=$
$\displaystyle(d+2)(1+2\xi)+M(1+\beta)\Big{\\{}(d+2)(1-\alpha_{0})-[(11+d)\alpha_{0}-5d-7]\xi\beta^{-1}\Big{\\}}+3(d+3)\xi^{2}\beta^{-1}$
(149)
$\displaystyle+2M^{2}\left(2\alpha_{0}^{2}-\frac{d+3}{2}\alpha_{0}+d+1\right)\beta^{-1}(1+\beta)^{2}-(d+2)\beta^{-1}(1+\beta),$
$\displaystyle C$ $\displaystyle=$ $\displaystyle
2M^{2}(1+\beta)^{2}\left(2\alpha_{0}^{2}-\frac{d+3}{2}\alpha_{0}+d+1\right)\left[d+2+(d+5)\beta\right]-M(1+\beta)\Big{\\{}\xi[d+2+(d+5)\beta][(11+d)\alpha_{0}$
(150)
$\displaystyle-5d-7]+\beta[20+d(15-7\alpha_{0})+d^{2}(1-\alpha_{0})-28\alpha_{0}]+(d+2)^{2}(1-\alpha_{0})\Big{\\}}+3(d+3)\xi^{2}[d+2+(d+5)\beta]$
$\displaystyle-2\xi[(d+2)^{2}+(24+11d+d^{2})\beta]+(d+2)\beta[d+8+(d+3)\beta]-(d+2)(1+\beta)[d+2+(d+3)\beta],$
$\displaystyle D$ $\displaystyle=$
$\displaystyle(d+2)(2\xi-\beta)+M(1+\beta)\Big{\\{}(d+2)(1-\alpha_{0})+[(11+d)\alpha_{0}-5d-7]\lambda\Big{\\}}-3(d+3)\xi^{2}$
(151)
$\displaystyle-2M^{2}\left(2\alpha_{0}^{2}-\frac{d+3}{2}\alpha_{0}+d+1\right)(1+\beta)^{2}+(d+2)(1+\beta).$

In Eqs. (C)–(151), $\xi=M_{0}(1-\tau_{0}^{-1})$.

##  Appendix D Derivatives of $\tau_{0}$ with respect to $\theta$ and $\phi$

It is quite apparent that to determine the transport coefficients $D^{T*}[2]$
and $D^{*}[2]$ one needs to evaluate the derivatives
$\Delta_{\theta,0}=\partial\tau_{0}/\partial\theta$ and
$\Delta_{\phi,0}=\partial\tau_{0}/\partial\phi$ in the steady state. To
evaluate them, we start from Eq. (69) of Ref. Gómez González, Khalil, and
Garzó, 2020 in the tracer limit ($x_{0}\to 0$):

$\Omega^{*}\theta\frac{\partial\tau_{0}}{\partial\theta}=-\tau_{0}\Omega^{*}+\Omega_{0}^{*},$
(152)

where

$\Omega^{*}=2\gamma^{*}(\theta^{-1}-1)-\zeta^{*},\quad\Omega_{0}^{*}=2\gamma_{0}^{*}(\theta^{-1}-\tau_{0})-\tau_{0}\zeta_{0}^{*},$
(153)

and $\zeta_{0}^{*}=\zeta_{0}/\nu$ is given by Eq. (40). Since in the steady
state, $\Omega^{*}=\Omega_{0}^{*}=0$, then according to Eq. (152) the
derivative $\partial\tau_{0}/\partial\theta$ becomes indeterminate. As in
previous works, Gómez González, Khalil, and Garzó (2020); Khalil and Garzó
(2013) the above problem can be solved by employing l’Hôpital’s rule. Thus, we
take first the derivative with respect to $\theta$ in both sides of Eq. (152)
and then take the steady-state condition ($\Omega^{*}=\Omega_{0}^{*}=0$).
After some algebra, we achieve the following expression for
$\Delta_{\theta,0}$:

$\Delta_{\theta,0}=\frac{\gamma^{*}_{0}\theta^{-1}(3\theta^{-1}-\tau_{0})-\gamma^{*}\tau_{0}\theta^{-1}(3\theta^{-1}-1)}{\gamma^{*}(3\theta^{-1}-1)-2\gamma_{0}^{*}-\zeta^{*}_{0}+\beta\frac{\partial\zeta_{0}^{*}}{\partial\beta}}.$
(154)

Once the derivative is known, the derivative
$\Delta_{\phi,0}\equiv\partial\tau_{0}/\partial\phi$ can be computed in a
similar way. The expression is

$\Delta_{\phi,0}=\frac{\left[\zeta^{*}\frac{\partial\ln\chi^{(0)}}{\partial\phi}-2\frac{\partial\lambda}{\partial\phi}\theta^{-1/2}(\theta^{-1}-1)\right](\theta\Delta_{\theta,0}+\tau_{0})+2\frac{\partial\lambda_{0}}{\partial\phi}\theta^{-1/2}(\theta^{-1}-\tau_{0})-\tau_{0}\zeta_{0}^{*}\frac{\partial\ln\chi_{0}^{(0)}}{\partial\phi}}{2\gamma_{0}^{*}+\zeta^{*}_{0}-\beta\frac{\partial\zeta_{0}^{*}}{\partial\beta}}.$
(155)

Note that in the limit of mechanically equivalent particles,
$\gamma^{*}=\gamma_{0}^{*}$, $\tau_{0}=1$ and hence, according to Eq. (154)
$\Delta_{\theta,0}=0$. As a result, since $\chi^{(0)}=\chi_{0}^{(0)}$ and
$\zeta^{*}=\zeta^{*}_{0}$, Eq (154) yields $\Delta_{\phi,0}=0$.

## References

  * Koch (1990) D. L. Koch, “Kinetic theory for a monodisperse gas-solid suspension,” Phys. Fluids A 2, 1711–1722 (1990).
  * Gidaspow (1994) D. Gidaspow,  _Multiphase Flow and Fluidization_ (Academic Press, 1994).
  * Jackson (2000) R. Jackson,  _The Dynamics of Fluidized Particles_ (Cambridge University Press, New York, 2000).
  * Koch and Hill (2001) D. L. Koch and R. J. Hill, “Inertial effects in suspensions and porous-media flows,” Annu. Rev. Fluid Mech. 33, 619–647 (2001).
  * Zhao and Wang (2021) B. Zhao and J. Wang, “Kinetic theory of polydisperse gas-solid flow: Navier-stokes transport coefficients,” Phys. Fluids 33, 103322 (2021).
  * Tsao and Koch (1995) H.-K. Tsao and D. L. Koch, “Simple shear flows of dilute gas–solid suspensions,” J. Fluid Mech. 296, 211–245 (1995).
  * Sangani  _et al._ (1996) A. S. Sangani, G. Mo, H.-K. Tsao, and D. L. Koch, “Simple shear flows of dense gas-solid suspensions at finite Stokes numbers,” J. Fluid Mech. 313, 309–341 (1996).
  * Wylie  _et al._ (2009) J. J. Wylie, Q. Zhang, Y. Li, and X. Hengyi, “Driven inelastic-particle systems with drag,” Phys. Rev. E 79, 031301 (2009).
  * Parmentier and Simonin (2012) J.-F. Parmentier and O. Simonin, “Transition models from the quenched to ignited states for flows of inertial particles suspended in a simple sheared viscous fluid,” J. Fluid Mech. 711, 147–160 (2012).
  * Heussinger (2013) C. Heussinger, “Shear thickening in granular suspensions: Interparticle friction and dynamically correlated clusters,” Phys. Rev. E 88, 050201 (R) (2013).
  * Wang  _et al._ (2014) T. Wang, M. Grob, A. Zippelius, and M. Sperl, “Active microrheology of driven granular particles,” Phys. Rev. E 89, 042209 (2014).
  * Saha and Alam (2017) S. Saha and M. Alam, “Revisiting ignited-quenched transition and the non-Newtonian rheology of a sheared dilute gas-solid suspension,” J. Fluid Mech. 833, 206–246 (2017).
  * Alam, Saha, and Gupta (2019) M. Alam, S. Saha, and R. Gupta, “Unified theory for a sheared gas-solid suspension: from rapid granular suspension to its small-Stokes-number limit,” J. Fluid Mech. 870, 1175–1193 (2019).
  * Saha and Alam (2020) S. Saha and M. Alam, “Burnett-order constitutive relations, second moment anisotropy and co-existing states in sheared dense gas-solid suspensions,” J. Fluid Mech. 887, A9 (2020).
  * Garzó  _et al._ (2012) V. Garzó, S. Tenneti, S. Subramaniam, and C. M. Hrenya, “Enskog kinetic theory for monodisperse gas-solid flows,” J. Fluid Mech. 712, 129–168 (2012).
  * Gómez González, Khalil, and Garzó (2020) R. Gómez González, N. Khalil, and V. Garzó, “Enskog kinetic theory for multicomponent granular suspensions,” Phys. Rev. E 101, 012904 (2020).
  * Gómez González and Garzó (2019) R. Gómez González and V. Garzó, “Transport coefficients for granular suspensions at moderate densities,” J. Stat. Mech. 093204 (2019).
  * Chapman and Cowling (1970) S. Chapman and T. G. Cowling,  _The Mathematical Theory of Nonuniform Gases_ (Cambridge University Press, Cambridge, 1970).
  * Sela and Goldhirsch (1998) N. Sela and I. Goldhirsch, “Hydrodynamic equations for rapid flows of smooth inelastic spheres to Burnett order,” J. Fluid Mech. 361, 41–74 (1998).
  * Serero  _et al._ (2006) D. Serero, I. Goldhirsch, S. H. Noskowicz, and M. L. Tan, “Hydrodynamics of granular gases and granular gas mixtures,” J. Fluid Mech. 554, 237–258 (2006).
  * Garzó, Dufty, and Hrenya (2007) V. Garzó, J. W. Dufty, and C. M. Hrenya, “Enskog theory for polydisperse granular mixtures. I. Navier–Stokes order transport,” Phys. Rev. E 76, 031303 (2007).
  * Garzó, Hrenya, and Dufty (2007) V. Garzó, C. M. Hrenya, and J. W. Dufty, “Enskog theory for polydisperse granular mixtures. II. Sonine polynomial approximation,” Phys. Rev. E 76, 031304 (2007).
  * Garzó (2021) V. Garzó, “Comment on “Kinetic theory models for granular mixtures with unequal granular temperature: Hydrodynamic velocity” [Phys. Fluids 33, 043321 (2021)],” Phys. Fluids 33, 089101 (2021).
  * García Chamorro and Garzó (2023) M. García Chamorro and V. Garzó, “Assessment of kinetic theories for moderately dense granular binary mixtures: Shear viscosity coefficient,” Phys. Fluids 35, 027121 (2023).
  * Solsvik and Manger (2021) J. Solsvik and E. Manger, “Kinetic theory models for granular mixtures with unequal granular temperature: Hydrodynamic velocity,” Phys. Fluids 33, 043321 (2021).
  * Garzó and Montanero (2004) V. Garzó and J. M. Montanero, “Diffusion of impurities in a granular gas,” Phys. Rev. E 69, 021301 (2004).
  * Garzó and Montanero (2007) V. Garzó and J. M. Montanero, “Navier–Stokes transport coefficients of $d$-dimensional granular binary mixtures at low-density,” J. Stat. Phys. 129, 27–58 (2007).
  * Garzó and Vega Reyes (2009) V. Garzó and F. Vega Reyes, “Mass transport of impurities in a moderately dense granular gas,” Phys. Rev. E 79, 041303 (2009).
  * Garzó and Vega Reyes (2012) V. Garzó and F. Vega Reyes, “Segregation of an intruder in a heated granular gas,” Phys. Rev. E 85, 021308 (2012).