###### Abstract

The Boltzmann kinetic equation for dilute granular suspensions under simple
(or uniform) shear flow (USF) is considered to determine the non-Newtonian
transport properties of the system. In contrast to previous attempts based on
a coarse-grained description, our suspension model accounts for the real
collisions between grains and particles of the surrounding molecular gas. The
latter is modeled as a bath (or thermostat) of elastic hard spheres at a given
temperature. Two independent but complementary approaches are followed to
reach exact expressions for the rheological properties. First, the Boltzmann
equation for the so-called inelastic Maxwell models (IMM) is considered. The
fact that the collision rate of IMM is independent of the relative velocity of
the colliding spheres allows us to exactly compute the collisional moments of
the Boltzmann operator without the knowledge of the distribution function.
Thanks to this property the transport properties of the sheared granular
suspension can be _exactly_ determined. As a second approach, a
Bhatnagar–Gross–Krook (BGK)-type kinetic model adapted to granular suspensions
is solved to compute the velocity moments and the velocity distribution
function of the system. The theoretical results (which are given in terms of
the coefficient of restitution, the reduced shear rate, the reduced background
temperature, and the diameter and mass ratios) show in general a good
agreement with the approximate analytical results derived for inelastic hard
spheres (IHS) by means of Grad’s moment method and with computer simulations
performed in the Brownian limiting case ($m/m_{g}\to\infty$, where $m_{g}$ and
$m$ are the masses of the particles of the molecular and granular gases,
respectively). In addition, as expected the IMM and BGK results show that the
temperature and non-Newtonian viscosity exhibit and $S$ shape in a plane of
stress-strain rate (discontinuous shear thickening, DST). The DST effect
becomes more pronounced as the mass ratio $m/m_{g}$ increases.

###### keywords:

Granular suspensions; Boltzmann kinetic equation; inelastic Maxwell models;
BGK-type kinetic model; Non-Newtonian transport properties

1 1 0 Received: date; Accepted: date; Published: date Exact results for non-
Newtonian transport properties in sheared granular suspensions: inelastic
Maxwell models and BGK-type kinetic model Rubén Gómez González 1,†, Vicente
Garzó 2,† Rubén Gómez González and Vicente Garzó Correspondence:
<EMAIL_ADDRESS>URL: https://fisteor.cms.unex.es/investigadores/vicente-
garzo-puertos/ These authors contributed equally to this work.

##  1 Introduction

A very usual way of assessing the effect of the surrounding fluid on the
dynamics properties of solid particles is through an effective fluid-solid
force Koch (1990); Gidaspow (1994); Fullmer and Hrenya (2017); Jackson (2000).
In some models, this force is simply proportional to the velocity particle
(Stokes linear drag law) Tsao and Koch (1995); Sangani et al. (1996); Wylie et
al. (2009); Heussinger (2013); Chamorro et al. (2015); Saha and Alam (2017);
Alam et al. (2019); Saha and Alam (2020). This type of force attempts to mimic
the energy dissipated by grains due to their friction on the interstitial
viscous gas. A more sophisticated model Garzó et al. (2012); Hayakawa et al.
(2017); Hayakawa and Takada (2019); Gómez González and Garzó (2019a, b);
Hayakawa et al. (2020); Takada et al. (2020); Gómez González et al. (2020);
Takada et al. (2023) incorporates also a Langevin-like stochastic term that
accounts for the energy transferred to grains due to their “collisions” with
particles of the background gas. However, although this coarse-grained
approach has provided reliable results in the past, it would be desirable to
consider a suspension model that takes into account the real collisions
between grains and particles of the surrounding (molecular) gas. This sort of
suspension model (which has been inspired in a previous work of Biben _et al_.
Baldassarri et al. (2002)) has been recently proposed Gómez González and Garzó
(2022). In this model, granular particles are assumed to be sufficiently
rarefied so that, they do not disturb the state of the molecular (background)
gas. As a consequence, the interstitial gas may be treated as a _thermostat_
at the temperature $T_{g}$. Moreover, although the concentration (mole
fraction) of grains is quite small, apart from the elastic collisions between
solid and molecular gas particles one has to consider the _inelastic_
collisions between grains themselves. This model can be useful to analyze
transport properties in particle-laden suspensions Subramaniam (2020) where
very dilute particles (like aerosols) are immersed in a fluid (like air).

The rheological properties of a granular suspension under simple (or uniform)
shear flow (USF) has been recently determined Gómez González et al. (2023). In
contrast to previous attempts Tsao and Koch (1995); Sangani et al. (1996);
Wylie et al. (2009); Heussinger (2013); Chamorro et al. (2015); Saha and Alam
(2017); Hayakawa et al. (2017); Hayakawa and Takada (2019); Alam et al.
(2019); Saha and Alam (2020); Takada et al. (2020, 2023), the results obtained
in Ref. Gómez González et al. (2023) have been derived from the collisional
model proposed in Ref. Gómez González and Garzó (2022). On the other hand, a
limitation of these results is that they have been _approximately_ obtained by
employing Grad’s moment method Grad (1949), namely, a method based on the
truncation of a series expansion of the velocity distribution function in
(orthogonal) Sonine polynomials. The use of this approximate method is
essentially motivated by the form of the collision rate for inelastic hard
spheres (IHS) appearing inside the Boltzmann collision operator. The collision
rate for IHS is proportional to the magnitude of the normal component of the
relative velocity of the two spheres that are about to collide. This velocity
dependence of the collision rate for IHS prevents the possibility of deriving
exact expressions for the transport properties in the USF problem, even in the
case of elastic collisions.

A possible way of overcoming the technical difficulty of the hard-sphere
kernel is to consider the so-called inelastic Maxwell models (IMM). As for the
conventional Maxwell molecules Chapman and Cowling (1970); Garzó and Santos
(2003), the collision rate of IMM is independent of the relative velocity of
the two colliding spheres Ben-Naim and Krapivsky (2003). The use of IMM
instead of IHS opens up the possibility of getting exact analytical results of
the Boltzmann equation in some specific non-equilibrium situations, like the
USF. In particular, the knowledge of the collisional moments of the Boltzmann
equation for IMM Garzó and Santos (2007); Sánchez Romero and Garzó (2023)
enables a clear exploration of the impact of inelasticity on the non-Newtonian
transport properties of the granular suspension without introducing
uncontrolled approximations.

Another possible alternative for obtaining exact results is to consider a
kinetic model that retains the relevant physical properties of the Boltzmann
collision operator but turns out to be more tractable than the true kinetic
equation. This route has been widely employed in the past in the case of
molecular dilute gases Garzó and Santos (2003) where it has been shown that
several exact solutions in far from equilibrium states agree very well with
Monte Carlo simulations of the Boltzmann equation. Here, as in previous works
Gómez González and Garzó (2019b), we will consider a kinetic model for
granular suspensions Vega Reyes et al. (2007) to complement the theoretical
expressions obtained from the Boltzmann equation for IMM. Since this kinetic
model is based on the well known Bhatnagar–Gross–Krook (BGK) model Cercignani
(1988) for molecular gases, we will referred to it as a BGK-type kinetic
model.

The objective of this paper is to determine the rheological properties of
granular particles immersed in a bath of elastic hard spheres under USF. At a
macroscopic level, the USF is characterized by constant number densities for
solid and gas particles, a uniform temperature, and a (common) linear velocity
profile $U_{g,x}=U_{x}=ay$, where $a$ is the constant shear rate. Here,
$\mathbf{U}_{g}$ and $\mathbf{U}$ denote the mean flow velocities of the
molecular and granular gases, respectively. Since we are interested here in
the steady state where the system admits a non-Newtonian hydrodynamic
description, an external thermostat force (proportional to the peculiar
velocity) must be introduced to keep constant the temperature $T_{g}$ of the
molecular gas.

The use of IMM as well as a BGK-type kinetic model allows us to exactly
compute the rheological properties of the granular suspension. These
properties are expressed as nonlinear functions of the (reduced) shear rate
$a^{*}=a/\gamma$ (where $\gamma$ is a drift coefficient characterizing the
friction of solid particles on the viscous gas), the coefficient of
restitution, the (reduced) background temperature $T_{g}^{*}$, and the
diameter $\sigma/\sigma_{g}$ and mass $m/m_{g}$ ratios. Here, $\sigma_{g}$ and
$m_{g}$ are the diameter and mass of the particles of the molecular gas,
respectively, while $\sigma$ and $m$ are the diameter and mass of the solid
particles, respectively. As occurs for IHS Gómez González et al. (2023), our
results show that the kinetic granular temperature and the non-Newtonian
viscosity exhibit a discontinuous shear thickening (DST) effect for
sufficiently large values of the mass ratio $m/m_{g}$. In fact, in the
Brownian limiting case ($m/m_{g}\to\infty$) the expressions of the rheological
properties derived here reduce to those previously obtained Gómez González and
Garzó (2019b) from a coarse-grained description based on the Fokker–Planck
operator. This agreement justifies the use of this latter approach to analyze
the DST effect in dilute granular suspensions Pelargonio and Zaccone (2023).

Apart from the transport properties (which are related with the second-degree
velocity moments), the explicit forms of the higher degree velocity moments as
well as the velocity distribution function of the granular gas have been also
obtained from the BGK model. This is one of the main advantages of using a
kinetic model instead of the true Boltzmann equation. Our results show in
particular that the fourth-degree moments of the distribution function also
exhibit a DST effect. With respect to the velocity distribution function, as
expected we find that its distortion from equilibrium is more significant as
both the mass $m/m_{g}$ and diameter $\sigma/\sigma_{g}$ ratios depart from 1.
In addition, a comparison between the BGK results and numerical solutions of
the Boltzmann equation from the direct simulation Monte Carlo (DSMC) method
Bird (1994) for IHS shows a generally good qualitative agreement between both
approaches.

The plan of the paper is as follows. The Boltzmann kinetic equation for a
granular gas immersed in a bath of elastic hard spheres under USF is presented
in Section 2. The balance equations for the temperatures of the molecular and
granular gases are also displayed. Section 3 deals with the calculations
carried out for IMM of the rheological properties of the granular suspension.
While a shear thinning effect is always found for the nonlinear shear
viscosity of the molecular gas, the corresponding shear viscosity of the
granular gas exhibits a DST effect for sufficiently large values of the mass
ratio $m/m_{g}$. The results derived from the BGK-type kinetic model are
provided in Section 4 while a comparison between the theoretical results
obtained for IHS, IMM and BGK model is displayed in Section 5 for several
systems. Our results highlight a good agreement for the rheology between the
three different approaches. Moreover, theoretical results obtained from IMM
and BGK model are also compared against computer simulations in the Brownian
limit ($m/m_{g}\to\infty$) showing a good agreement. The paper is closed in
Section 6 with some concluding remarks.

##  2 Boltzmann kinetic equation for sheared granular suspensions

We consider a set of solid particles (granular gas) of mass $m$ and diameter
$\sigma$ which are immersed in a solvent (molecular gas) constituted by
particles of mass $m_{g}$ and diameter $\sigma_{g}$. As usual, the granular
gas is modeled as a gas of hard disks ($d=2$) or spheres ($d=3$) with
inelastic collisions. In the simplest model, the inelasticity of collisions is
characterized by a constant (positive) coefficient of normal restitution
$\alpha\leq 1$, where $\alpha=1$ refers to elastic collisions. On the other
hand, collisions between solid particles and particles of the molecular gas
are elastic. We also assume that the number density of grains is much smaller
than that of solvent so that, the state of the latter is not perturbed by the
presence of the former. In these conditions, we can treat the molecular gas as
a bath or _thermostat_ at the temperature $T_{g}$ (once the parameters of the
system, specifically the shear rate, have been set). Moreover, although the
granular gas is sufficiently rarefied, we take into account the collisions
among grains in its corresponding Boltzmann kinetic equation.

We assume that the system (granular particles plus solvent) is subjected to
USF. As said in the Introduction section, this state is characterized by
constant densities $n_{g}$ and $n$, uniform temperatures $T_{g}$ and $T$, and
by a (common) linear profile of the $x$ component of the flow velocities along
the $y$ axis:

$n_{g}\equiv\text{const},\quad n\equiv\text{const},$ (1) $\nabla T_{g}=\nabla
T=0,$ (2) $U_{g,i}=U_{i}=a_{ij}r_{j},\quad a_{ij}=a\delta_{ix}\delta_{jy},$
(3)

$a$ being the _constant_ shear rate. Here, $n_{g}$, $\mathbf{U}_{g}$, and
$T_{g}$ are the number density, the mean flow velocity, and the temperature,
respectively, of the molecular gas. In terms of its one-particle velocity
distribution function $f_{g}(\mathbf{r},\mathbf{v};t)$, these hydrodynamic
fields are defined as

$\left\\{n_{g},n_{g}\mathbf{U}_{g},dn_{g}T_{g}\right\\}=\int
d\mathbf{v}\left\\{1,\mathbf{v},m_{g}V^{2}\right\\}f_{g}(\mathbf{v}),$ (4)

where $\mathbf{V}=\mathbf{v}-\mathbf{U}$ is the peculiar velocity. Note that
in Eq. (4) the Boltzmann constant $k_{\text{B}}=1$. We will take this value
throughout the paper for the sake of simplicity. In addition, in Eqs. (1)–(3),
$n$, $\mathbf{U}$, and $T$ denote the number density, the mean flow velocity,
and the (granular) temperature, respectively, of the granular gas. They are
defined as

$\left\\{n,n\mathbf{U},dnT\right\\}=\int
d\mathbf{v}\left\\{1,\mathbf{v},mV^{2}\right\\}f(\mathbf{v}).$ (5)

Since the only spatial gradient present in the USF problem is the shear rate,
the pressure tensor

$\mathsf{P}_{g}=m_{g}\int d\mathbf{v}\;\mathbf{V}\mathbf{V}f_{g}(\mathbf{v})$
(6)

of the molecular gas, and the pressure tensor

$\mathsf{P}=m\int d\mathbf{v}\;\mathbf{V}\mathbf{V}f(\mathbf{v})$ (7)

of the granular gas are the relevant fluxes in the problem. They provide
information on the transport of momentum across the system. Our main target is
to determine $\mathsf{P}_{g}$ and $\mathsf{P}$ for arbitrary shear rates.

One of the main advantages of the USF at a microscopic level is that it
becomes a spatially homogeneous state when the velocities of the particles are
referred to a Lagrangian frame moving with the linear velocity
$U_{i}=a_{ij}r_{j}$. In this new frame and in the steady state, the
distribution functions of the molecular and granular gases adopt the form

$f_{g}(\mathbf{r},\mathbf{v})=f_{g}(\mathbf{V}),\quad
f(\mathbf{r},\mathbf{v})=f(\mathbf{V}).$ (8)

In addition, as the state of the solvent is not perturbed by the solid
particles, the temperature $T_{g}$ in the USF state increases in time due to
the viscous heating term $-aP_{xy}>0$. Thus, as usual in nonequilibrium
molecular dynamics simulations Evans and Morriss (1990), an external
nonconservative force (thermostat) is introduced in the molecular gas to
achieve a stationary state. Among the different possibilities, for simplicity,
a force proportional to the particle velocity (Gaussian thermostat) of the
form $\mathbf{F}_{g}=-m_{g}\xi\mathbf{V}$ is considered in this paper. The
parameter $\xi$ is chosen to be a function of the shear rate by the condition
that $T_{g}$ reaches a constant value in the long time limit. Analogously, the
granular gas is also subjected to this kind of Gaussian thermostat (i.e.,
$\mathbf{F}=-m\xi\mathbf{V}$), where $\xi$ is the same quantity for the
solvent and the solid particles.

Under the above conditions, in the low-density regime, the distribution
function $f_{g}(\mathbf{V})$ of the molecular gas obeys the nonlinear (closed)
Boltzmann equation

$-aV_{y}\frac{\partial f_{g}}{\partial
V_{x}}-\xi\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f_{g}=J_{g}[\mathbf{V}|f_{g},f_{g}],$
(9)

while the distribution function $f(\mathbf{V})$ of the granular gas obeys the
kinetic equation

$-aV_{y}\frac{\partial f}{\partial
V_{x}}-\xi\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f=J[\mathbf{V}|f,f]+J_{\text{BL}}[\mathbf{V}|f,f_{g}].$
(10)

Here, $J_{g}[f_{g},f_{g}]$ and $J[f,f]$ are the nonlinear Boltzmann collision
operators for the molecular and granular gases, respectively, and
$J_{\text{BL}}[f,f_{g}]$ is the linear Boltzmann–Lorentz collision operator
Brilliantov and Pöschel (2004); Garzó (2019). The balance equations for
$T_{g}$ and $T$ can easily obtained by multiplying both sides of Eqs. (9) and
(10) by $m_{g}V^{2}$ and $mV^{2}$, respectively, and integrating over
velocity. The results are

$-aP_{g,xy}=d\xi p_{g},$ (11) $-aP_{xy}-\frac{d}{2}p\zeta_{g}=d\xi
p+\frac{d}{2}p\zeta,$ (12)

where $p_{g}=n_{g}T_{g}$ and $p=nT$ are the hydrostatic pressures of the
molecular and granular gases, respectively, and the partial production rates
$\zeta$ and $\zeta_{g}$ are defined, respectively, as

$\zeta=-\frac{m}{dnT}\int
d\mathbf{v}\;V^{2}\;J[\mathbf{v}|f,f],\quad\zeta_{g}=-\frac{m}{dnT}\int
d\mathbf{v}\;V^{2}\;J_{\text{BL}}[\mathbf{v}|f,f_{g}].$ (13)

The cooling rate $\zeta$ gives the rate of kinetic energy loss due to
inelastic collisions between particles of the granular gas. It vanishes for
elastic collisions $(\alpha=1)$. The term $\zeta_{g}$ gives the transfer of
kinetic energy between the particles of the granular gas and the solvent. This
term vanishes when the granular and molecular gases are at the same
temperature ($T_{g}=T$). Equation (11) implies that in the steady state the
viscous heating term ($-aP_{g,xy}>0$) is exactly balanced by the heat
extracted in the gas by the external thermostat. On the other hand, since
$\zeta_{g}$ can be positive or negative, Eq. (12) implies that in the steady
state the term $-aP_{xy}-(d/2)p\zeta_{g}$ is exactly compensated for the
cooling terms arising from collisional dissipation ($\zeta p$) and the
thermostat term ($\xi p$).

The USF state is in general a non-Newtonian state characterized by shear-rate
dependent transport coefficients. In particular, one can define the non-
Newtonian shear viscosity of the molecular gas as

$\eta_{g}=-\frac{P_{g,xy}}{a}.$ (14)

Analogously, the non-Newtonian shear viscosity of the granular gas is given by

$\eta=-\frac{P_{xy}}{a}.$ (15)

In addition, beyond the Navier–Stokes domain, normal stress differences are
expected in the USF. This means that $P_{g,xx}\neq P_{g,yy}\neq P_{g,zz}$ and
$P_{xx}\neq P_{yy}\neq P_{zz}$.

It is quite evident that the evaluation of the rheological properties of the
molecular and granular gases requires the knowledge of the pressure tensors
$\mathsf{P}_{g}$ and $\mathsf{P}$. The non-zero elements of these tensors can
be obtained by multiplying by $m_{g}\mathbf{V}\mathbf{V}$ and
$m\mathbf{V}\mathbf{V}$ both sides of Eqs. (9) and (10), respectively, and
integrating over $\mathbf{V}$. However, to achieve explicit forms for
$\mathsf{P}_{g}$ and $\mathsf{P}$, one has to compute the collisional moments

$\mathsf{A}_{g}=m_{g}\int d{\bf V}{\bf V}{\bf V}J_{g}[f_{g},f_{g}],$ (16)
$\mathsf{B}=m\int d{\bf V}{\bf V}{\bf V}J[f,f],\quad\mathsf{C}=m\int d{\bf
V}{\bf V}{\bf V}J_{\text{BL}}[f,f_{g}].$ (17)

In the case of IHS, the collisional moments $\mathsf{A}_{g}$, $\mathsf{B}$,
and $\mathsf{C}$ cannot be exactly computed. As said in the Introduction
section, a good estimate of them for IHS has been made in Ref. Gómez González
et al. (2023) by means of Grad’s moment method Grad (1949). This method is
based on the expansion of the distributions $f_{g}(\mathbf{V})$ and
$f(\mathbf{V})$ in a complete set of orthogonal polynomials; the coefficients
being the corresponding velocity moments of those distributions. The above
expansion generates an infinite hierarchy of moment equations that must be
truncated at a given order. This truncation allows one to arrive to a closed
set of coupled equations for the velocity moments that can be recursively
solved. Thus, since the results derived in Ref. Gómez González et al. (2023)
for the rheological properties of molecular and granular gases are
approximated, it is convenient to revisit the problem and determine the
_exact_ expressions for the non-Newtonian transport properties of the granular
suspension. To achieve such exact forms, two independent approaches will be
considered in this paper: (i) the Boltzmann kinetic equation for IMM and (ii)
a BGK-type kinetic model for IHS. This task will be carried out in the next
two Sections.

##  3 Rheology from inelastic Maxwell models

In this Section, we will consider IMM, namely, a collisional model where the
collision rate of the two colliding spheres are independent of their relative
velocity. In this case, the Boltzmann collision operator $J_{g}[f_{g},f_{g}]$
111This is a simple version of the Boltzmann collision operator for Maxwell
molecules of the molecular gas can be written as Cercignani (1988); Garzó and
Santos (2003)

$J_{g}[\mathbf{v}_{1}|f_{g},f_{g}]=\frac{\nu_{g}^{\text{M}}}{n_{g}S_{d}}\int
d\mathbf{v}_{2}\int d\widehat{\boldsymbol{\sigma}}\,\left[f_{g}({\bf
V}_{1}^{\prime\prime})f_{g}({\bf V}_{2}^{\prime\prime})-f_{g}({\bf
V}_{1})f_{g}({\bf V}_{2})\right],$ (18)

where $S_{d}=2\pi^{d/2}/\Gamma(\frac{d}{2})$ is the total solid angle in $d$
dimensions and $\nu_{g}^{\text{M}}$ is an independent-velocity collision
frequency. In Eq. (18), the primes on the velocities denote the initial values
$\\{{\bf V}_{1}^{\prime\prime},{\bf V}_{2}^{\prime\prime}\\}$ that lead to
$\\{{\bf V}_{1},{\bf V}_{2}\\}$ following a binary collision:

$\mathbf{V}_{1}^{{}^{\prime\prime}}=\mathbf{V}_{1}-(\widehat{\boldsymbol{\sigma}}\cdot{\bf
g}_{12})\widehat{\boldsymbol{\sigma}},\quad\mathbf{V}_{2}^{{}^{\prime\prime}}=\mathbf{V}_{2}+(\widehat{\boldsymbol{\sigma}}\cdot{\bf
g}_{12})\widehat{\boldsymbol{\sigma}}.$ (19)

The effective collision frequency $\nu_{g}^{\text{M}}$ can be seen as a free
parameter of the model to be chosen to get agreement with the properties of
interest of the original Boltzmann equation for IHS. For instance, to
correctly capture the velocity dependence of the original IHS collision rate,
we can assume that the IMM collision rate is proportional to $\sqrt{T}_{g}$.

In the context of IMM, the inelastic Boltzmann collision operator $J[f,f]$ is
Ben-Naim and Krapivsky (2003); Garzó (2019)

$J[\mathbf{V}|f,f]=\frac{\nu^{\text{M}}}{nS_{d}}\int d{\bf V}_{2}\int
d\widehat{\boldsymbol{\sigma}}\,\left[\alpha^{-1}f({\bf
V}_{1}^{\prime\prime})f({\bf V}_{2}^{\prime\prime})-f({\bf V}_{1})f({\bf
V}_{2})\right],$ (20)

while the Boltzmann–Lorentz collision operator $J_{\text{BL}}[f,f_{g}]$ is
defined as Cercignani (1988); Garzó and Santos (2003)

$J_{\text{BL}}[\mathbf{V}|f,f_{g}]=\frac{\nu_{0}^{\text{M}}}{nS_{d}}\int d{\bf
V}_{2}\int d\widehat{\boldsymbol{\sigma}}\,\left[f({\bf
V}_{1}^{\prime\prime})f_{g}({\bf V}_{2}^{\prime\prime})-f({\bf
V}_{1})f_{g}({\bf V}_{2})\right].$ (21)

The relationship between
$(\mathbf{V}_{1}^{\prime\prime},\mathbf{V}_{2}^{\prime\prime})$ and
$(\mathbf{V}_{1},\mathbf{V}_{2})$ in Eq. (20) is

$\mathbf{V}_{1}^{\prime\prime}=\mathbf{V}_{1}-\frac{1+\alpha}{2\alpha}(\widehat{{\boldsymbol{\sigma}}}\cdot{\mathbf{g}}_{12})\widehat{\boldsymbol{\sigma}},\quad\mathbf{V}_{2}^{\prime\prime}=\mathbf{V}_{2}+\frac{1+\alpha}{2\alpha}(\widehat{{\boldsymbol{\sigma}}}\cdot{\mathbf{g}}_{12})\widehat{\boldsymbol{\sigma}},$
(22)

while in Eq. (21) is

$\mathbf{V}_{1}^{\prime\prime}=\mathbf{V}_{1}-2\mu_{g}(\widehat{{\boldsymbol{\sigma}}}\cdot{\mathbf{g}}_{12})\widehat{\boldsymbol{\sigma}},\quad\mathbf{V}_{2}^{\prime\prime}=\mathbf{V}_{2}+2\mu(\widehat{{\boldsymbol{\sigma}}}\cdot{\mathbf{g}}_{12})\widehat{\boldsymbol{\sigma}}$
(23)

where

$\mu_{g}=\frac{m_{g}}{m+m_{g}},\quad\mu=\frac{m}{m+m_{g}}.$ (24)

In addition, as in the case of the collision frequency $\nu_{g}^{\text{M}}$,
the collision frequencies $\nu^{\text{M}}$ and $\nu_{0}^{\text{M}}$ for
granular-granular and granular-molecular collisions, respectively, can be
chosen to optimize the agreement with the results derived from IHS. We will
chose them later.

As mentioned in previous works on IMM Santos and Garzó (2007); Garzó and
Santos (2007); Gómez González and Garzó (2019a), the main advantage of
computing the collisional moments of the Boltzmann operator for Maxwell models
(both elastic and inelastic models) is that they can be exactly provided in
terms of the velocity moments of the distribution functions without the
explicit knowledge of the latter. This property has been exploited to compute
the second, third, and fourth-degree collisional moments of IMM for
monocomponent Garzó and Santos (2007) and multicomponent Garzó (2003); Sánchez
Romero and Garzó (2023) granular gases. The exact knowledge of the second-
degree collisional moments allow us to get exact expressions for the
rheological properties of the molecular and granular gases. Let us evaluate
separately the rheology of both gases.

###  3.1 Rheological properties of the molecular gas

The pressure tensor of the molecular gas is defined by Eq. (6). To obtain the
non-zero elements of this tensor, one multiplies both sides of the Boltzmann
equation (9) by $m_{g}V_{i}V_{j}$ and integrates over velocity. The result is

$a_{ik}P_{g,kj}+a_{jk}{P}_{g,ki}+2\xi
P_{g,ij}=-\frac{2\nu_{g}^{\text{M}}}{d+2}\left({P}_{g,ij}-p_{g}\delta_{ij}\right),$
(25)

where $p_{g}=n_{g}T_{g}$ is the hydrostatic pressure of the molecular gas.
Upon obtaining Eq. (25), use has been made of the result Garzó and Santos
(2007)

${A}_{g,ij}=m_{g}\int d{\bf
V}V_{i}V_{j}J_{g}[f_{g},f_{g}]=-\frac{2\nu_{g}^{\text{M}}}{d+2}\left({P}_{g,ij}-p_{g}\delta_{ij}\right).$
(26)

The (reduced) elements of the pressure tensor
$\mathsf{P}_{g}^{*}=\mathsf{P}_{g}/(n_{g}T_{g})$ can be easily obtained from
Eq. (25). They are given by

$P_{g,xx}^{*}=\frac{1}{1+2\widetilde{\xi}}\left[1+\frac{2\widetilde{a}^{2}}{(1+2\widetilde{\xi})^{2}}\right],\quad
P_{g,yy}^{*}=P_{g,zz}^{*}=\frac{1}{1+2\widetilde{\xi}},\quad
P_{g,xy}^{*}=-\frac{\widetilde{a}}{(1+2\widetilde{\xi})^{2}}.$ (27)

Here, we have introduced the quantities

$\widetilde{a}=\frac{a}{\widetilde{\nu}_{g}^{\text{M}}},\quad\widetilde{\xi}=\frac{\xi}{\widetilde{\nu}_{g}^{\text{M}}},\quad\widetilde{\nu}_{g}^{\text{M}}=\frac{2}{d+2}\nu_{g}^{\text{M}}.$
(28)

The constraint $P_{g,xx}^{*}+(d-1)P_{g,yy}^{*}=d$ leads to a cubic equation
relating $\widetilde{\xi}$ and $\widetilde{a}$:

$\widetilde{a}^{2}=d\widetilde{\xi}(1+2\widetilde{\xi})^{2}.$ (29)

The real root of Eq. (29) gives the shear-rate dependence of
$\widetilde{\xi}(\widetilde{a})$. It is given by Ikenberry and Truesdell
(1956); Garzó and Santos (2003)

$\widetilde{\xi}(\widetilde{a})=\frac{2}{3}\sinh^{2}\left[\frac{1}{6}\cosh^{-1}\left(1+\frac{27}{d}\widetilde{a}^{2}\right)\right].$
(30)

Comparison between the results derived here for $P_{g,ij}^{*}$ with those
recently Gómez González et al. (2023) obtained for IHS by means of Grad’s
moment method Grad (1949) shows that both results are identical if the
effective collision frequency $\nu_{g}^{\text{M}}$ is given by

$\nu_{g}^{\text{M}}=\frac{4\pi^{(d-1)/2}}{\Gamma\left(\frac{d}{2}\right)}n_{g}\sigma_{g}^{d-1}\sqrt{\frac{T_{g}}{m_{g}}}.$
(31)

Henceforth, we will take the choice (31) for $\nu_{g}^{\text{M}}$.

From Eqs. (27), one can identify the (dimensionless) non-Newtonian shear
viscosity
$\eta_{g}^{*}=\nu_{g}^{\text{M}}\eta_{g}/p_{g}=-P_{g,xy}^{*}/\widetilde{a}$
and the (dimensionless) normal stress difference
$\Psi_{g}^{*}=P_{g,xx}^{*}-P_{g,yy}^{*}$ as

$\eta_{g}^{*}=\frac{1}{(1+2\widetilde{\xi})^{2}},\quad\Psi_{g}^{*}=\frac{2\widetilde{a}^{2}}{(1+2\widetilde{\xi})^{3}}.$
(32)

Note that the results derived here for Maxwell molecules yield
$P_{g,yy}=P_{g,zz}$. This result contrasts with the one obtained for hard
spheres by numerically solving the Boltzmann equation by means of the DSMC
method Bird (1994) where it has been shown that $P_{g,yy}\neq P_{g,zz}$.
However, the difference $P_{g,yy}-P_{g,zz}$ found in the simulations is in
general quite small Chamorro et al. (2015).

It is also important to remark that in the case of Maxwell molecules there is
an exact equivalence between the description with and without the drag force
$-m\xi\mathbf{V}$. Nevertheless, for non-Maxwell molecules this type of force
does not play a neutral role in the transport properties of the system Dufty
et al. (1986).

Figure 1: Plot of the nonlinear shear viscosity $\eta_{g}^{*}$ and the normal
stress difference $\Psi_{g}^{*}$ for hard spheres ($d=3$) as functions of the
(reduced) shear rate $\widetilde{a}$. Symbols refer to the DSMC results for
hard spheres.

The shear-rate dependence of $\eta_{g}^{*}$ and $\Psi_{g}^{*}$ is plotted in
Fig. 1 for a three-dimensional system ($d=3$). As expected, the nonlinear
viscosity $\eta_{g}^{*}$ decreases with increasing the (reduced) shear rate
$\widetilde{a}$ (shear thinning effect). The opposite effect is observed for
the normal stress difference function $\Psi_{g}^{*}$ since it increases with
the shear rate. Figure 1 also highlights the excellent agreement found between
the theoretical results for IMM with those obtained by numerically solving the
Boltzmann equation for hard spheres from the DSMC method Bird (1994).

###  3.2 Rheological properties of the granular gas

As in the case of the molecular gas, the rheology of the granular gas can be
also determined by multiplying both sides of Eq. (10) by $mV_{i}V_{j}$ and
integrating over $\mathbf{V}$. After some algebra, one achieves the result

$\displaystyle a_{ik}P_{kj}+a_{jk}{P}_{ki}+2\xi
P_{ij}=-\nu^{\text{M}}\left[\nu_{\eta}^{*}P_{ij}+p\left(\zeta^{*}-\nu_{\eta}^{*}\right)\delta_{ij}\right]-\frac{4}{d}\nu_{0}^{\text{M}}\mu_{g}$
$\displaystyle\times\Big{[}\Big{(}1-\frac{2}{d+2}\mu_{g}\Big{)}P_{ij}-\frac{2}{d+2}\mu\frac{n}{n_{g}}P_{g,ij}-\frac{d}{d+2}\mu_{g}nT_{g}\Big{(}\chi+\frac{m}{m_{g}}\Big{)}\delta_{ij}\Big{]},$
(33)

where use has been made of the results Garzó and Santos (2007); Garzó (2003);
Sánchez Romero and Garzó (2023)

$B_{ij}=m\int d{\bf
V}V_{i}V_{j}J[f,f]=-\nu^{\text{M}}\left[\nu_{\eta}^{*}P_{ij}+p\left(\zeta^{*}-\nu_{\eta}^{*}\right)\delta_{ij}\right],$
(34) $\displaystyle C_{ij}$ $\displaystyle=$ $\displaystyle m\int d{\bf
V}V_{i}V_{j}J_{\text{BL}}[f,f_{g}]=-\frac{4}{d}\nu_{0}^{\text{M}}\mu_{g}\Big{[}\Big{(}1-\frac{2}{d+2}\mu_{g}\Big{)}P_{ij}-\frac{2}{d+2}\mu\frac{n}{n_{g}}P_{g,ij}$
(35)
$\displaystyle-\frac{d}{d+2}\mu_{g}nT_{g}\Big{(}\chi+\frac{m}{m_{g}}\Big{)}\delta_{ij}\Big{]}.$

In Eq. (3.2),

$\zeta^{*}=\frac{\zeta}{\nu^{\text{M}}}=\frac{1-\alpha^{2}}{2d}$ (36)

is the (reduced) cooling rate for the granular gas, $\chi=T/T_{g}$ is the
temperature ratio, and

$\nu_{\eta}^{*}=\frac{\left(d+1-\alpha\right)(1+\alpha)}{d(d+2)}.$ (37)

The partial cooling rate $\zeta_{g}$ can be exactly obtained from Eq. (34) as

$\zeta_{g}=\frac{4\nu_{0}^{\text{M}}}{d}\mu_{g}\left[1-\mu_{g}\left(1+\theta\right)\right],$
(38)

where

$\theta=\frac{mT_{g}}{m_{g}T}$ (39)

is the ratio of the mean square velocities of granular and molecular gas
particles. The forms (36) and (38) can be employed to fix the values of the
free parameters $\nu^{\text{M}}$ and $\nu_{0}^{\text{M}}$. They are chosen
under the criterion that $\zeta$ and $\zeta_{g}$ of IMM are the same as that
of IHS of diameters $\sigma$ and $\sigma_{0}$. In this latter case, the above
cooling rates are estimated by using Grad’s approximation Gómez González et
al. (2023). In this approximation,

$\zeta^{\text{IHS}}=\frac{2\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n\sigma^{d-1}\sqrt{\frac{T}{m}}(1-\alpha^{2}),$
(40)
$\zeta_{g}^{\text{IHS}}=\frac{8\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n_{g}\overline{\sigma}^{d-1}\mu_{g}\left(\frac{1+\theta}{\theta}\right)^{1/2}\sqrt{\frac{2T_{g}}{m_{g}}}\left[1-\mu_{g}\left(1+\theta\right)\right],$
(41)

where $\overline{\sigma}=(\sigma+\sigma_{g})/2$. Equations (38), (40), and
(41) yield the identities

$\nu^{\text{M}}=\frac{4\pi^{(d-1)/2}}{\Gamma\left(\frac{d}{2}\right)}n\sigma^{d-1}\sqrt{\frac{T}{m}},\quad\nu_{0}^{\text{M}}=\frac{2\pi^{(d-1)/2}}{\Gamma\left(\frac{d}{2}\right)}n_{g}\overline{\sigma}^{d-1}\left(\frac{1+\theta}{\theta}\right)^{1/2}\sqrt{\frac{2T_{g}}{m_{g}}}.$
(42)

To compare with the rheological properties of IHS Gómez González et al.
(2023), it is convenient at this level of the description to identify the
friction (or drift) coefficient $\gamma$ appearing in the Brownian limiting
case ($m/m_{g}\to\infty$) when the molecular gas is at equilibrium. In fact,
this limiting case is the situation considered when one employs a coarse-
grained approach Koch (1990); Gidaspow (1994); Jackson (2000); Fullmer and
Hrenya (2017) to assess the impact of the interstitial gas on the dynamics
properties of grains. In this limiting case, the expression (35) of the
collisional moment $C_{ij}$ reduces to

$C_{ij}^{\text{Br}}=-\frac{8\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}nT_{g}n_{g}\overline{\sigma}^{d-1}\left(\frac{m_{g}}{m}\right)^{1/2}\sqrt{\frac{2T_{g}}{m}}\left(P_{k\ell}^{*}-\delta_{k\ell}\right),$
(43)

where $P_{ij}^{*}=P_{ij}/(nT_{g})$ and we have taken into account that in the
Brownian limit $\mu_{g}\to m_{g}/m$ and $(1+\theta)/\theta\to 1$ in the
expression (42) of $\nu_{0}^{\text{M}}$. The form of $C_{ij}$ derived in Ref.
Gómez González and Garzó (2019b) by replacing the Boltzmann–Lorentz
collisional operator (21) by the Fokker–Planck operator

$\gamma\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f+\gamma\frac{T_{g}}{m}\frac{\partial^{2}f}{\partial
V^{2}}$ (44)

is

$C_{ij}=-2\gamma nT_{g}\left(P_{k\ell}^{*}-\delta_{k\ell}\right).$ (45)

Comparison between Eqs. (43) and (45) allows us to identify $\gamma$ for IMM
as

$\gamma=\frac{4\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n_{g}\overline{\sigma}^{d-1}\left(\frac{m_{g}}{m}\right)^{1/2}\left(\frac{2T_{g}}{m}\right)^{1/2}.$
(46)

The expression (46) for the friction coefficient $\gamma$ for IMM is the same
as the one obtained for IHS Gómez González and Garzó (2022); Gómez González et
al. (2023).

We are now in conditions to obtain the nonzero elements of the (reduced)
pressure tensor $P_{ij}^{*}$. From Eqs. (3.2) and (46), one gets the equation

$\displaystyle a_{ik}^{*}P_{kj}^{*}+a_{jk}^{*}{P}_{ki}^{*}+2\xi^{*}P_{ij}^{*}$
$\displaystyle=$
$\displaystyle-\nu^{*\text{M}}\left[\nu_{\eta}^{*}P_{ij}^{*}+\chi(\zeta^{*}-\nu_{\eta}^{*})\delta_{ij}\right]-2\mu\left(\frac{1+\theta}{\theta}\right)^{1/2}$
(47) $\displaystyle\times\left(X\delta_{ij}+YP_{ij}^{*}+ZP_{g,ij}^{*}\right),$

where $a_{ij}^{*}=a_{ij}/\gamma$, $\xi^{*}=\xi/\gamma$,

$\nu^{*\text{M}}=\frac{\nu^{\text{M}}}{\gamma}=\frac{2^{d+1}d}{\sqrt{\pi}}\phi\sqrt{\chi
T_{g}^{*}},$ (48)

and

$X=-\frac{d}{d+2}\mu\left(\frac{1+\theta}{\theta}\right),\quad
Y=1-\frac{2\mu_{g}}{d+2},\quad Z=-\frac{2\mu}{d+2}.$ (49)

Here,

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

is the solid volume fraction of the granular gas,
$T_{g}^{*}=T_{g}/m\sigma^{2}\gamma^{2}$, and upon deriving Eq. (48) use has
been made of the identity

$\frac{n\sigma^{d-1}}{n_{g}\overline{\sigma}^{d-1}}=\frac{2^{d+\frac{3}{2}}}{\sqrt{\pi}}\left(\frac{m_{g}}{m}\right)^{1/2}\phi\sqrt{T_{g}^{*}}.$
(51)

As occurs for the rheology of the molecular gas, Eq. (47) shows that the
diagonal elements of the pressure tensor $P_{ij}^{*}$ orthogonal to the shear
plane $xy$ are equal to $P_{yy}^{*}$ (i.e.,
$P_{yy}^{*}=P_{zz}^{*}=\cdots=P_{dd}^{*}$). This implies that the $xx$ element
is given by $P_{xx}^{*}=d\chi-(d-1)P_{yy}^{*}$. The $yy$ and $xy$ elements of
the (reduced) pressure tensor can be written as

$P_{yy}^{*}=\frac{\Omega_{yy}}{\nu_{yy}},\quad
P_{xy}^{*}=\frac{\Omega_{xy}-a^{*}P_{yy}^{*}}{\nu_{yy}},$ (52)

where

$\nu_{yy}=2\xi^{*}+\nu^{*\text{M}}\nu_{\eta}^{*}+2\mu\left(\frac{1+\theta}{\theta}\right)^{1/2}Y,$
(53)
$\Omega_{yy}=-\nu^{*\text{M}}\chi\left(\zeta^{*}-\nu_{\eta}^{*}\right)-2\mu\left(\frac{1+\theta}{\theta}\right)^{1/2}\left(X+ZP_{g,yy}^{*}\right),$
(54)
$\Omega_{xy}=-2\mu\left(\frac{1+\theta}{\theta}\right)^{1/2}ZP_{g,xy}^{*}.$
(55)

Note that the elements of the pressure tensor $P_{g,yy}^{*}$ and
$P_{g,xy}^{*}$ of the molecular gas must be expressed in terms of the
(reduced) shear rate $a^{*}$ and the (reduced) thermostat parameter $\xi^{*}$.
For doing it, one has to take into account the relationships between
$\widetilde{a}$ and $\widetilde{\xi}$ with $a^{*}$ and $\xi^{*}$,
respectively. They are given by
$\widetilde{a}=(\gamma/\widetilde{\nu}_{g}^{\text{M}})a^{*}$ and
$\widetilde{\xi}=(\gamma/\widetilde{\nu}_{g}^{\text{M}})\xi^{*}$, where

$\frac{\gamma}{\widetilde{\nu}_{g}^{\text{M}}}=\frac{d+2}{\sqrt{2}d}\left(\frac{\overline{\sigma}}{\sigma_{g}}\right)^{d-1}\frac{m_{g}}{m}.$
(56)

The equation defining the temperature ratio $\chi$ can be easily derived from
Eq. (47) as

$\frac{2}{d}a^{*}P_{xy}^{*}+2\xi^{*}\chi=-\nu^{*\text{M}}\chi\zeta^{*}+2\mu^{2}\left(\frac{1+\theta}{\theta}\right)^{1/2}\left(1-\chi\right).$
(57)

From Eqs. (52) and (57), one gets finally $a^{*}$ in terms of the parameter
space of the system:

$a^{*}=\sqrt{\frac{d}{2}\frac{2\mu^{2}\left(1+\theta^{-1}\right)^{1/2}(1-\chi)-\left(\nu^{*\text{M}}\zeta^{*}+2\xi^{*}\right)\chi}{\frac{\Omega_{xy}/a^{*}}{\nu_{yy}}-\frac{\Omega_{yy}}{\nu_{yy}^{2}}}}.$
(58)

As happens in the case of IHS Gómez González et al. (2023), the temperature
ratio $\chi$ cannot be expressed in Eq. (58) as an explicit function of the
(reduced) shear rate and the remaining parameters of the system. On the other
hand, for given values of the parameter space
$\Xi\equiv\left(\chi,\alpha,\sigma/\sigma_{g},m/m_{g},\phi,T_{g}^{*}\right)$,
the temperature ratio can be implicitly determined from the physical solution
to Eq. (58).

###  3.3 Brownian limit

Before illustrating the shear-rate dependence of the rheological properties of
the molecular gas for arbitrary values of the mass ratio $m/m_{g}$, it is
convenient to check the consistency of the present results with those derived
in Ref. Gómez González and Garzó (2019b) for IMM by using the Fokker–Planck
operator (44). This consistency is expected to apply in the Brownian limit
$m/m_{g}\to\infty$. In this limiting case, at a given value of the (reduced)
shear rate $a^{*}$, $\theta\to\infty$, $\gamma/\nu_{g}^{\text{M}}\propto
m_{g}/m\to 0$, $\widetilde{a}\propto m_{g}/m\to 0$,
$\widetilde{\xi}\propto\widetilde{a}^{2}\propto(m_{g}/m)^{2}\to 0$, and
$\xi^{*}=\widetilde{\xi}(\widetilde{\nu}_{g}^{\text{M}}/\gamma)\propto
m_{g}/m\to 0$. Consequently, $P_{g,ij}^{*}=\delta_{ij}$ and the expressions of
$P_{yy}^{*}$, $P_{xy}^{*}$, and $a^{*}$ are

$P_{yy}^{*}=\frac{2-\nu^{*\text{M}}\chi\left(\zeta^{*}-\nu_{\eta}^{*}\right)}{2+\nu^{*\text{M}}\nu_{\eta}^{*}},\quad
P_{xy}^{*}=-\frac{2-\nu^{*\text{M}}\chi\left(\zeta^{*}-\nu_{\eta}^{*}\right)}{\left(2+\nu^{*\text{M}}\nu_{\eta}^{*}\right)^{2}}a^{*},$
(59)
$a^{*}=\sqrt{\frac{d}{2}\frac{\nu^{*\text{M}}\zeta^{*}+2\left(1-\chi^{-1}\right)}{\nu^{*\text{M}}\left(\nu_{\eta}^{*}-\zeta^{*}\right)+2\chi^{-1}}}\left(2+\nu^{*\text{M}}\nu_{\eta}^{*}\right).$
(60)

Equations (59) and (60) are consistent with Eqs. (32), (33), and (35) of Ref.
Gómez González and Garzó (2019b). It is important to note that, to assess
consistency with the Fokker–Planck results, the size ratio has been kept
constant or proportional to the mass ratio so that $\xi^{*}\to 0$.

##  4 Rheology from a BGK-type kinetic model of the Boltzmann equation

We consider in this Section the results derived for the USF from a BGK-type
kinetic model of the Boltzmann equation Brey et al. (1999); Vega Reyes et al.
(2007). In the problem for the granular suspension considered here, one has to
replace the true Boltzmann operators $J_{g}[f_{g},f_{g}]$, $J[f,f]$, and
$J_{\text{BL}}[f,f_{g}]$ by simpler relaxation terms that retain the relevant
physical properties of those operators but are more tractable than the true
kinetic equations. As in the case of IMM, let us determine separately the
rheological properties of the molecular and granular gases by starting from
these kinetic models.

###  4.1 Rheological properties of the molecular gas

In the case of the molecular gas, the Boltzmann collision operator
$J_{g}[f_{g},f_{g}]$ is replaced by the conventional BGK kinetic model
Bhatnagar et al. (1954); Cercignani (1988):

$J_{g}[f_{g},f_{g}]\to-\nu_{g}\left(f_{g}-f_{g}^{\text{M}}\right),$ (61)

where $\nu_{g}$ is an effective velocity-independent collision frequency and
$f_{g}^{\text{M}}$ is the Maxwellian distribution

$f_{g}^{\text{M}}(\mathbf{V})=n_{g}\left(\frac{m_{g}}{2\pi
T_{g}}\right)^{d/2}\exp\left(-\frac{m_{g}V^{2}}{2T_{g}}\right).$ (62)

Thus, according to Eq. (9), the velocity distribution function
$f_{g}(\mathbf{V})$ obeys the BGK kinetic equation

$-aV_{y}\frac{\partial f_{g}}{\partial
V_{x}}-\xi\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f_{g}=-\nu_{g}\left(f_{g}-f_{g}^{\text{M}}\right).$
(63)

The nonzero elements of the pressure tensor $P_{g,ij}$ can be easily obtained
from Eq. (63) by multiplying both sides of this equation by $m_{g}V_{i}V_{j}$
and integrating over $\mathbf{V}$. The BGK expressions of the (reduced)
elements of the pressure tensor $P_{g,ij}^{*}$ are given by Eqs. (27) with the
replacement $\widetilde{\nu}_{g}^{\text{M}}\to\nu_{g}$. As a consequence, the
results derived from the BGK equation for the rheological properties agree
with those obtained from the Boltzmann equation for IHS when

$\nu_{g}=\widetilde{\nu}_{g}^{\text{M}}=\frac{8\pi^{(d-1)/2}}{(d+2)\Gamma\left(\frac{d}{2}\right)}n_{g}\sigma_{g}^{d-1}\sqrt{\frac{T_{g}}{m_{g}}}.$
(64)

###  4.2 Rheological properties of the granular gas

In the case of the molecular gas, we consider the kinetic model proposed by
Vega Reyes _et al._ Vega Reyes et al. (2007) for granular mixtures where the
relaxation terms for the collision operators $J[f,f]$ and
$J_{\text{BL}}[f,f_{g}]$ are defined, respectively, as

$J[f,f]\to-\nu^{\prime}(\alpha)\left(f-f^{\text{M}}\right)+\frac{\epsilon}{2}\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f,$
(65) $J_{\text{BL}}[f,f_{g}]\to-\nu\left(f-\widetilde{f}_{g}\right).$ (66)

In Eqs. (65) and (66),

$f^{\text{M}}(\mathbf{V})=n\left(\frac{m}{2\pi
T}\right)^{d/2}\exp\left(-\frac{mV^{2}}{2T}\right)$ (67)

is the Maxwellian distribution of the granular gas,

$\epsilon=\frac{2\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n\sigma^{d-1}\sqrt{\frac{T}{m}}(1-\alpha^{2}),$
(68)

and the reference distribution function $\widetilde{f}_{g}(\mathbf{V})$ is
taken to be the same as in the well-known Gross and Krook (GK) model for
molecular (elastic) mixtures Gross and Krook (1956), i.e.,

$\widetilde{f}_{g}(\mathbf{V})=n\left(\frac{m}{2\pi\tilde{T}}\right)^{d/2}\exp\left(-\frac{mV^{2}}{2\widetilde{T}}\right).$
(69)

In Eqs. (65) and (66), the quantities $\nu^{\prime}$, $\nu$, and
$\widetilde{T}$ are chosen to optimize the agreement with some properties of
interest of the Boltzmann equation for IHS. The usual way of obtaining the
above parameters is to ensure that the kinetic model reproduces the
collisional transfer equations of momentum and energy for elastic collisions
($\alpha=1$). However, since $\mathbf{U}=\mathbf{U}_{g}$ in the USF, we only
have one constraint (the one associated with the transfer of energy) instead
of two, so that $\widetilde{T}$ and $\nu$ admit more than one form. Here, we
propose a choice (see Appendix A for more technical details) that leads to an
excellent agreement with the results obtained for IHS from Grad’s moment
method Gómez González et al. (2023). More specifically, we take the following
values of $\widetilde{T}$ and $\nu$:

$\widetilde{T}=T_{g},\quad\nu=\frac{8\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n_{g}\overline{\sigma}^{d-1}\frac{mm_{g}}{(m+m_{g})^{2}}\left(\frac{2T_{g}}{m_{g}}+\frac{2T}{m}\right)^{1/2}.$
(70)

It still remains to completely define the model to chose the effective
collision frequency $\nu^{\prime}$. It is defined here to reproduce the
collisional moment

$\int d\mathbf{v}\;mV_{i}V_{j}J[f,f]$ (71)

of the Boltzmann equation for IHS when one takes Grad’s trial distribution for
$f$ Gómez González et al. (2023). This leads to the expression

$\nu^{\prime}=\frac{2\pi^{(d-1)/2}}{d(d+2)\Gamma\left(\frac{d}{2}\right)}n\sigma^{d-1}\sqrt{\frac{T}{m}}\left(1+\alpha\right)\left[d+1+\left(d-1\right)\alpha\right].$
(72)

Therefore, the BGK kinetic equation for the sheared granular gas is given by

$-aV_{y}\frac{\partial f}{\partial
V_{x}}-\xi\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f=-\nu^{\prime}\left(f-f^{\text{M}}\right)+\frac{\epsilon}{2}\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f-\nu\left(f-\widetilde{f}_{g}\right),$
(73)

where $\nu$ and $\nu^{\prime}$ are defined by Eqs. (70) and (72),
respectively, and the Maxwellian distribution $\widetilde{f}_{g}$ is given by
Eq. (69) with $\widetilde{T}=T_{g}$.

The possibility of determining all the velocity moments of the distribution
function is likely one of the main advantages of employing a kinetic model
instead of the true Boltzmann equation. In the USF problem, it is convenient
to define the general velocity moments

$M_{k_{1},k_{2},k_{3}}=\int
d\mathbf{V}V^{k_{1}}_{x}V_{y}^{k_{2}}V_{z}^{k_{3}}f(\mathbf{V}).$ (74)

Although we are here essentially interested in the three-dimensional case, we
will compute the velocity moments for $d=3$ and $d=2$. Note that for hard
disks ($d=2$), $k_{3}=0$ in Eq. (74) since the $z$-axis is meaningless. The
hierarchy of moment equations can be obtained by multiplying Eq. (73) by
$V^{k_{1}}_{x}V_{y}^{k_{2}}V_{z}^{k_{3}}$ and integrating over $\mathbf{V}$.
The result is

$ak_{1}M_{k_{1}-1,k_{2}+1,k_{3}}+(\nu^{\prime}+\nu+k\lambda)M_{k_{1},k_{2},k_{3}}=N_{k_{1},k_{2},k_{3}},$
(75)

where $\lambda=\xi+\epsilon/2$, $k=k_{1}+k_{2}+k_{3}$, and

$N_{k_{1},k_{2},k_{3}}=n\left(\frac{2T_{g}}{m}\right)^{k/2}\Big{(}\nu+\chi^{k/2}\nu^{\prime}\Big{)}M_{k_{1},k_{2},k_{3}}^{L}.$
(76)

In Eq. (75), for hard spheres ($d=3$)

$M_{k_{1},k_{2},k_{3}}^{L}=\int
d\mathbf{c}\;c_{x}^{k_{1}}c_{y}^{k_{2}}c_{z}^{k_{3}}e^{-c^{2}}=\pi^{-3/2}\Gamma\left(\frac{k_{1}+1}{2}\right)\Gamma\left(\frac{k_{2}+1}{2}\right)\Gamma\left(\frac{k_{3}+1}{2}\right),$
(77)

if $k_{1}$, $k_{2}$, and $k_{3}$ are even, being zero otherwise. For hard
disks ($d=2$),

$M_{k_{1},k_{2},0}^{L}=\pi^{-1}\Gamma\left(\frac{k_{1}+1}{2}\right)\Gamma\left(\frac{k_{2}+1}{2}\right)$
(78)

if $k_{1}$ and $k_{2}$ are even, being zero otherwise.

The solution to Eq. (75) can be written as (see Appendix B of Ref. Gómez
González and Garzó (2019a) for some details)

$M_{k_{1},k_{2},k_{3}}=\sum_{q=0}^{k_{1}}\frac{k_{1}!}{(k_{1}-q)!}\frac{(-a)^{q}}{(\nu^{\prime}+\nu+k\lambda)^{1+q}}N_{k_{1}-q,k_{2}+q,k_{3}}.$
(79)

The nonzero elements of the pressure tensor $P_{ij}$ can be easily obtained
from Eq. (79). In dimensionless form, the BGK-expressions of the elements of
$P_{ij}^{*}=P_{ij}/nT_{g}$ are

$P_{yy}^{*}=P_{zz}^{*}=\frac{\nu+\nu^{\prime}\chi}{\nu^{\prime}+\nu+2\lambda},\quad
P_{xy}^{*}=-\frac{\nu+\nu^{\prime}\chi}{\left(\nu^{\prime}+\nu+2\lambda\right)^{2}}a,$
(80)
$P_{xx}^{*}=\frac{\nu+\nu^{\prime}\chi}{\nu^{\prime}+\nu+2\lambda}\Bigg{[}1+\frac{2a^{2}}{\left(\nu^{\prime}+\nu+2\lambda\right)^{2}}\Bigg{]}.$
(81)

The (steady) temperature ratio $\chi=T/T_{g}$ can be obtained from the
constraint $P_{xx}^{*}+(d-1)P_{yy}^{*}=d\chi$. This yields the implicit
equation

$a^{*}=\left(\nu^{{}^{\prime}*}+\nu^{*}+2\lambda^{*}\right)\sqrt{\frac{d}{2}\chi\frac{\nu^{*}\left(1-\chi^{-1}\right)+2\lambda^{*}}{\nu^{*}+\nu^{{}^{\prime}*}\chi}},$
(82)

where $a^{*}=a/\gamma$,

$\nu^{{}^{\prime}*}=\frac{\nu^{\prime}}{\gamma}=\frac{2^{d}}{d+2}\phi\sqrt{\frac{\chi
T_{g}^{*}}{\pi}}(1+\alpha)\left[d+1+\left(d-1\right)\alpha\right],$ (83)
$\nu^{*}=\frac{\nu}{\gamma}=2\mu^{2}\left(\frac{1+\theta}{\theta}\right)^{1/2},\quad\lambda^{*}=\frac{\lambda}{\gamma}=\xi^{*}+\frac{\epsilon^{*}}{2},$
(84)
$\xi^{*}=\frac{\sqrt{2}d}{d+2}\left(\frac{\sigma_{g}}{\overline{\sigma}}\right)^{d-1}\frac{m}{m_{g}}\widetilde{\xi},\quad\epsilon^{*}=2^{d}\phi\sqrt{\frac{\chi
T_{g}^{*}}{\pi}}(1-\alpha^{2}).$ (85)

Here, $\widetilde{\xi}$ is given by Eq. (30) and we recall that
$T_{g}^{*}=T_{g}/m\sigma^{2}\gamma^{2}$.

###  4.3 Brownian limit

As in the case of IMM, it is quite interesting to consider the Brownian
limiting case $m/m_{g}\to\infty$. In this case, $P_{g,ij}^{*}=\delta_{ij}$,
$\mu\to 1$, $\theta\to\infty$, and $\xi^{*}\to 0$. Thus, following similar
steps as those made for IMM, one gets the expressions

$P_{yy}^{*}=\frac{2+\nu^{{}^{\prime}*}\chi}{2+\nu^{{}^{\prime}*}+\epsilon^{*}},\quad
P_{xy}^{*}=-\frac{2+\nu^{{}^{\prime}*}\chi}{(2+\nu^{{}^{\prime}*}+\epsilon^{*})^{2}}a^{*},$
(86)
$a^{*}=\left(2+\nu^{{}^{\prime}*}+\epsilon^{*}\right)\sqrt{\frac{d}{2}\chi\frac{2\left(1-\chi^{-1}\right)+\epsilon^{*}}{2+\nu^{{}^{\prime}*}\chi}}.$
(87)

###  4.4 Velocity distribution of the granular gas

Apart from obtaining all the velocity moments, the use of kinetic models allow
us in some cases to explicitly determine the velocity distribution function
$f(\mathbf{V})$. The BGK-type equation (73) reads

$\left(1-d\widehat{\lambda}-\widehat{a}V_{y}\frac{\partial}{\partial
V_{x}}-\widehat{\lambda}\mathbf{V}\cdot\frac{\partial}{\partial\mathbf{V}}\right)f(\mathbf{V})=f_{\text{R}}(\mathbf{V}),$
(88)

where $\widehat{a}=a/(\nu^{\prime}+\nu)$,
$\widehat{\lambda}=\lambda/(\nu^{\prime}+\nu)$, and

$f_{\text{R}}(\mathbf{V})=\frac{\nu^{\prime}}{\nu^{\prime}+\nu}f^{\text{M}}(\mathbf{V})+\frac{\nu}{\nu^{\prime}+\nu}\widetilde{f}_{g}(\mathbf{V}).$
(89)

The hydrodynamic solution to Eq. (88) can be formally written as

$\displaystyle f(\mathbf{V})$ $\displaystyle=$
$\displaystyle\Big{(}1-d\widehat{\lambda}-\widehat{a}V_{y}\frac{\partial}{\partial
V_{x}}-\widehat{\lambda}\mathbf{V}\cdot\frac{\partial}{\partial\mathbf{V}}\Big{)}^{-1}f_{\text{R}}(\mathbf{V})$
(90) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}ds\;e^{-(1-d\widehat{\lambda})s}e^{\widehat{a}sV_{y}\frac{\partial}{\partial
V_{x}}}e^{\widehat{\lambda}s\mathbf{V}\cdot\frac{\partial}{\partial\mathbf{V}}}f_{\text{R}}(\mathbf{V}).$

The action of the velocity operators in Eq. (90) on an arbitrary function
$g(\mathbf{V})$ is Gómez González and Garzó (2019a)

$e^{\widehat{a}sV_{y}\frac{\partial}{\partial
V_{x}}}g(V_{x},V_{y},V_{z})=g(V_{x}+\widehat{a}sV_{y},V_{y},V_{z}),$ (91)
$e^{\widehat{\lambda}s\mathbf{V}\cdot\frac{\partial}{\partial\mathbf{V}}}g(V_{x},V_{y},V_{z})=g\left(e^{\widehat{\lambda}s}V_{x},e^{\widehat{\lambda}s}V_{y},e^{\widehat{\lambda}s}V_{z}\right).$
(92)

The explicit form of the one-particle velocity distribution function can be
explicitly obtained when one takes into account in Eq. (90) the action of the
velocity operators given by Eqs. (91) and (92). The result can be written as

$f(\mathbf{V})=n\left(\frac{m}{2T}\right)^{d/2}\varphi(\mathbf{c}),\quad\mathbf{c}=\left(\frac{m}{2T}\right)^{1/2}\mathbf{V},$
(93)

where

$\displaystyle\varphi(\mathbf{c})$ $\displaystyle=$
$\displaystyle\pi^{-d/2}\int_{0}^{\infty}ds\;e^{-(1-d\widehat{\lambda})s}\Bigg{\\{}\frac{\nu^{\prime}}{\nu^{\prime}+\nu}\exp\Big{[}-e^{2\widehat{\lambda}s}\Big{(}\mathbf{c}+s\widehat{\boldsymbol{a}}\cdot\mathbf{c}\Big{)}^{2}\Big{]}$
(94) $\displaystyle+\frac{\nu}{\nu^{\prime}+\nu}\chi^{d/2}\exp\Big{[}-\chi
e^{2\widehat{\lambda}s}\Big{(}\mathbf{c}+s\widehat{\boldsymbol{a}}\cdot\mathbf{c}\Big{)}^{2}\Big{]}\Bigg{\\}}.$

Here, we have introduced the tensor
$\widehat{a}_{ij}=\widehat{a}\delta_{ix}\delta_{jy}$.

In order to illustrate the shear-rate dependence of the distribution function,
it is convenient to consider the marginal distribution for $d=3$:

$\displaystyle\varphi_{x}(c_{x})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\;dc_{y}\int_{-\infty}^{\infty}\;dc_{z}\;\varphi(\mathbf{c})$
(95) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\pi}}\int_{0}^{\infty}\;ds\;\frac{e^{-(1-\widehat{\lambda})s}}{\sqrt{1+\widehat{a}^{2}s^{2}}}\Bigg{\\{}\frac{\nu^{\prime}}{\nu^{\prime}+\nu}\exp\Big{(}-e^{2\widehat{\lambda}s}\frac{c_{x}^{2}}{1+\widehat{a}^{2}s^{2}}\Big{)}$
$\displaystyle+\frac{\nu}{\nu^{\prime}+\nu}\chi^{1/2}\exp\Big{(}-\chi
e^{2\widehat{\lambda}s}\frac{c_{x}^{2}}{1+\widehat{a}^{2}s^{2}}\Big{)}\Bigg{\\}}.$

In the Brownian limit, $\xi^{*}\to 0$, $\mu\to 1$, $\theta\to\infty$, and so
$\nu^{*}\to 2$, $\lambda^{*}\to\epsilon^{*}/2$, and

$\widehat{\lambda}\to\frac{\epsilon^{*}/2}{2+\nu^{{}^{\prime}*}},\quad\widehat{a}\to\frac{a^{*}}{2+\nu^{{}^{\prime}*}}.$
(96)

Thus, when $m/m_{g}\to\infty$, Eq. (95) becomes

$\displaystyle\varphi_{x}(c_{x})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\pi}}\int_{0}^{\infty}\;ds\;\frac{e^{-(1-\widehat{\lambda})s}}{\sqrt{1+\widehat{a}^{2}s^{2}}}\Bigg{\\{}\frac{\nu^{{}^{\prime}*}}{2+\nu^{{}^{\prime}*}}\exp\Big{(}-e^{2\widehat{\lambda}s}\frac{c_{x}^{2}}{1+\widehat{a}^{2}s^{2}}\Big{)}$
(97) $\displaystyle+\frac{2}{2+\nu^{{}^{\prime}*}}\chi^{1/2}\exp\Big{(}-\chi
e^{2\widehat{\lambda}s}\frac{c_{x}^{2}}{1+\widehat{a}^{2}s^{2}}\Big{)}\Bigg{\\}},$

where $\widehat{\lambda}$ and $\widehat{a}$ are given by Eqs. (96).

##  5 Comparison between IMM and BGK results

In Sections 3 and 4, we made use of the Boltzmann equation for IMM and the
BGK-type kinetic model to investigate the shear-rate dependence of rheological
properties in a sheared granular suspension. These properties are expressed in
terms of the coefficient of restitution $\alpha$, the reduced background
temperature $T_{g}^{*}$, and the diameter $\sigma/\sigma_{g}$ and mass
$m/m_{g}$ ratios. Additionally, there exists a residual dependence on density
through the volume fraction $\phi$. To avoid that, one could for instance
reduce the shear rate using the effective collision frequencies
$\nu^{\text{M}}(T)$ or $\nu(T)$. However, for consistency with simulations and
considering the background temperature $T_{g}$ as a known quantity, we opted
to employ $\gamma(T_{g})$ as the reference frequency.

Given that in this Section the second-degree moments of the distribution
function are compared with molecular dynamics (MD) simulations for IHS in the
Brownian limiting case Hayakawa et al. (2017), we set fixed values of
$T_{g}^{*}=1$ and $\phi=0.0052$ for subsequent analysis. The selection of
$T_{g}^{*}$ as a free parameter imposes a constraint between the diameter
$\sigma/\sigma_{g}$ and mass $m/m_{g}$ ratios Gómez González et al. (2023):

$\frac{\sigma}{\sigma_{g}}=\left[\left(\frac{\sqrt{\pi}}{4\sqrt{2}}\frac{n}{n_{g}}\sqrt{\frac{m}{m_{g}}}\frac{1}{\phi\sqrt{T_{g}^{*}}}\right)^{1/(d-1)}-1\right]^{-1}.$
(98)

This relation ensures convergence of results to those obtained via the
Fokker–Planck equation as $m/m_{g}\to\infty$ since $\xi^{*}\to 0$.
Furthermore, since we want to recover the results obtained in Ref. Gómez
González et al. (2023) derived from Grad’s method, we take $n/n_{g}=10^{-3}$
(rarefied granular gas).

Figure 2: Plots of the (steady) granular temperature $\chi$ (a), the non-
Newtonian shear viscosity $\eta^{*}$ (b), and the normal stress difference
$\Psi^{*}$ (c) as functions of the (reduced) shear rate $a^{*}$ for two
different values of the coefficient of restitution $\alpha$: 0.9 (left panel)
and 1 (right panel). The graphs represent four distinct mass ratio values
$m/m_{g}$: $10^{3}$ (yellow lines), $10^{4}$ (blue lines), $10^{6}$ (red
lines), and the Brownian limit (black lines). Here, $T_{g}^{*}=1$, $d=3$, and
$\phi=0.0052$. The solid lines correspond to the IMM results, the dashed lines
are the BGK-like results, and the dotted lines refer to Grad’s solution for
IHS. Symbols denote computer simulation results performed in the Brownian
limit: circles refer to the DSMC data obtained in this paper for IHS while
squares are MD results obtained in Ref. Hayakawa et al. (2017) for IHS.

The second-degree moments expressed through the reduced temperature $\chi$,
non-Newtonian shear viscosity $\eta^{*}$, and the normal stress difference
$\Psi^{*}$ are plotted in Fig. 2 for $\alpha=0.9$ and $1$. Here,
$\Psi^{*}=P_{xx}^{*}-P^{*}_{yy}=d\chi-dP^{*}_{yy}$. Equations (52) and (80)
provide analytical expressions for IMM and BGK-type kinetic model,
respectively, from which rheological properties are illustrated. Notably,
there is nearly perfect agreement between Grad’s solution for IHS as obtained
in Ref. Gómez González et al. (2023) and both IMM and BGK-type results for any
mass ratio, highlighting the ability of relatively simple models to capture
essential properties of granular suspensions.

In particular, a DST transition characterized by $S$-shaped curves becomes
more pronounced as the mass ratio $m/m_{g}$ increases. Specifically, the non-
Newtonian shear viscosity $\eta^{*}$ exhibits a discontinuous transition (at a
certain value of $a^{*}$) which intensifies as the particles of the granular
gas become heavier. The theoretical results are validated with MD simulations
Hayakawa et al. (2017) in the Brownian limiting case ($m/m_{g}\to 0$), showing
generally good agreement despite slight discrepancies in the transition zone.
Simulations suggest a sharper transition, likely due to the absence of
molecular chaos in highly non-equilibrium situations. To address this, DSMC
simulations for IHS are performed in the same limit, showing good agreement
with theoretical results and further emphasizing a more pronounced transition.
This phenomenon is likely attributable to a sudden growth of higher-order
moments causing the distribution function to strongly deviate from the
reference Maxwellian approximation (67). Some technical details of the
application of the DSMC method are available in the supplementary material of
Ref. Gómez González and Garzó (2022) [https://doi.org/10.1017/jfm.2022.410].

The simplicity of the BGK and IMM models enables exploration beyond second-
degree moments. Accordingly, we utilize the BGK-type kinetic equation to
compute the fourth-degree moments. Although similar calculations could be
performed in the case of IMM, we opt to omit them due to their extensive
analytical effort. Additionally, drawing insights from the Fokker–Planck model
Gómez González and Garzó (2019b) and dry granular gases Santos and Garzó
(2007), we anticipate potential divergences of the moments derived from IMM
under certain shear rate conditions. We focus our efforts on calculating the
following fourth-degree moments

$M^{*}_{4|0}=\frac{m^{2}}{nT_{g}^{2}}\int d{\bf V}V^{4}f(\mathbf{V}),$ (99)
$M_{2|xy}^{*}=\frac{m^{2}}{nT_{g}^{2}}\int d{\bf
V}V^{2}V_{x}V_{y}f(\mathbf{V}).$ (100)

Thus, in terms of the generic moments $M_{k_{1},k_{2},k_{3}}$ and according to
the expression (79), the moments $M^{*}_{4|0}$ and $M_{2|xy}^{*}$ are given by

$M^{*}_{4|0}=\frac{m^{2}}{nT_{g}^{2}}\left(M_{4,0,0}+M_{0,4,0}+M_{0,0,4}+2M_{2,2,0}+2M_{2,0,2}+2M_{0,2,2}\right),$
(101)
$M^{*}_{2|xy}=\frac{m^{2}}{nT_{g}^{2}}\left(M_{3,1,0}+M_{1,3,0}+M_{1,1,2}\right).$
(102)

Figure 3: Plots of the (reduced) fourth-degree moments
$M_{4|0}^{*}(a^{*})/M_{4|0}^{*}(0)$ (a) and $-M_{2|xy}^{*}(a^{*})$ (b) as
functions of the (reduced) shear rate $a^{*}$ obtained from the BGK-type
equation for two different values of the coefficient of restitution $\alpha$:
0.9 (left panel) and 1 (right panel). The graphs represent four distinct mass
ratio values $m/m_{g}$: $10^{3}$ (yellow lines), $10^{4}$ (blue lines),
$10^{6}$ (red lines), and the Brownian limit (black lines). Here,
$T_{g}^{*}=1$, $d=3$, and $\phi=0.0052$. Symbols refer to the DSMC results
obtained for IHS in this paper in the Brownian limit.

Figure 3(a) illustrates the ratio $M^{*}_{4|0}(a^{*})/M^{*}_{4|0}(0)$ as a
function of $a^{*}$ for $\alpha=0.9$ and $1$ and three values of the mass
ratio. We observe that variations in the mass ratio do not significantly alter
the trends observed in the Brownian limiting case Gómez González and Garzó
(2019b). An abrupt transition in the higher-order moments is evident within a
small region of $a^{*}$. Specifically, the kurtosis $M^{*}_{4|0}$ increases
with the mass ratio $m/m_{g}$ until it converges to the value obtained in the
Brownian limit. Consistent with the conclusions drawn in Ref. Gómez González
et al. (2023), an increase in the mass of the particles of the granular gas
results in an elevation of the granular temperature. Consequently, non-
equipartition accentuates and moves the suspension away from equilibrium,
leading to an increase in kurtosis as the distribution function deviates from
its Maxwellian form. Regarding the influence of collisional dissipation, we
observe that the effect of $\alpha$ on $M^{*}_{4|0}$ remains relatively
discrete. Figure 3(b) illustrates the shear-rate dependence of the (reduced)
moment $M_{2|xy}^{*}$. This moment vanishes in the absence of shear rate
($a^{*}=0$). Similar conclusions to those made for the moment $M^{*}_{4|0}$
can be drawn.

Theoretical predictions for the fourth-degree moments are compared against
DSMC simulations for IHS conducted in this paper in the Brownian limiting
case. A qualitative agreement is observed, although simulations suggest a
sharper transition. Some quantitative discrepancies are noticeable, which are
mainly disguised by the scale. To assess the reliability of the BGK-type
results, we focus on the region $0<a^{*}<1$, where all the fourth-degree
velocity moments of IMM are well-defined functions of the shear rate. In
addition, non-Newtonian effects are still significant in the range of values
of the (reduced) shear rate $a^{*}\leq 1$. To this purpose, Fig. 4 shows the
(reduced) fourth-degree moments $M_{4|0}^{*}(a^{*})/M_{4|0}^{*}(0)$ and
$M_{2|xy}^{*}(a^{*})$ for $\alpha=0.7$ and $1$. These moments are also
illustrated as obtained for IMM in the Brownian limit Gómez González and Garzó
(2019b). It is worth noting that the results derived in Ref. Gómez González
and Garzó (2019b) stem from considering an effective force modeling the
interstitial gas, diverging from the limit of a Boltzmann–Lorentz operator
modeled by a BGK-type equation as $m/m_{g}\to\infty$. Consequently, since DSMC
simulations employ the exact Fokker–Planck operator, they perfectly align with
the IMM results, while discrepancies emerge when compared with the BGK-type
results. It is noteworthy that the BGK-type model slightly overestimates the
deviation from the Newtonian situation ($a^{*}=0$), a phenomenon also observed
for molecular gases Garzó and Santos (1995). Moreover, non-Newtonian effects
become apparent even at low values of $a^{*}$.

Figure 4: Plots of the (reduced) fourth-degree moments
$M_{4|0}^{*}(a^{*})/M_{4|0}^{*}(0)$ (a) and $-M_{2|xy}^{*}(a^{*})$ (b) as
functions of the (reduced) shear rate $a^{*}$ obtained from the BGK-type
equation for two different values of the coefficient of restitution $\alpha$:
0.7 (solid lines) and 1 (dashed lines). The graphs represent four distinct
mass ratio values $m/m_{g}$: $5$ (yellow lines), $10$ (blue lines), $50$ (red
lines), and the Brownian limit (black lines). Here, $T_{g}^{*}=1$, $d=3$, and
$\phi=0.0052$. Symbols refer to the DSMC results for IHS in the Brownian limit
(squares for $\alpha=1$ and circles for $\alpha=0.7$). The green lines are the
IMM results as obtained in Ref. Gómez González and Garzó (2019a) in the
Brownian limit.  Figure 5: Plot of the ratio
$R_{x}(c_{x})=\varphi_{x}(c_{x})/(\pi^{-1/2}e^{-c_{x}^{2}})$ for $a^{*}=1$ as
a function of the (reduced) velocity $c_{x}$ for $\alpha=0.9$ and five
different values of the mass ratio $m/m_{g}$: $1$ (yellow line), $5$ (blue
line), $10$ (red line), $50$ (green line), and the Brownian limit (black
lines). Here, $T_{g}^{*}=1$, $d=3$, and $\phi=0.0052$. The dashed line refer
to the DSMC results for IHS in the Brownian limit.

Finally, in Fig. 5, the ratio
$R_{x}(c_{x})=\varphi_{x}(c_{x})/(\pi^{-1/2}e^{-c_{x}^{2}})$ is plotted for
$a^{*}=1$ and four different values of the mass ratio. Here, the marginal
distribution $\varphi_{x}(c_{x})$ is given by Eq. (95). It is evident that the
deviation from equilibrium ($R_{x}\neq 1$) becomes more significant as the
mass ratio $m/m_{g}$ increases. Moreover, a comparison between theory and DSMC
simulations reveals some disagreement in the BGK-type solution. Although the
relative difference of these discrepancies is relatively small (it is about
8%), this contradicts what has been observed in Ref. Chamorro et al. (2015),
where good agreement between the BGK solution and DSMC data is shown in the
region of thermal velocities.

##  6 Concluding remarks

In our study, we have explored the non-Newtonian transport properties of a
dilute granular suspension subjected to USF using the Boltzmann kinetic
equation. The particles are represented as $d$-dimensional hard spheres with
mass $m$ and diameter $\sigma$, immersed in an interstitial gas acting as a
thermostat at temperature $T_{g}$. Various models for granular suspensions
incorporate a gas-solid force to represent the influence of the external
fluid. While some models consider only isolated body resistance via a linear
drag law Tsao and Koch (1995); Sangani et al. (1996); Parmentier and Simonin
(2012); Heussinger (2013); Wang et al. (2014); Saha and Alam (2017); Alam et
al. (2019); Saha and Alam (2020), others Garzó et al. (2012); Puglisi et al.
(1998) include an additional Langevin-type stochastic term. In this paper, we
consider a suspension model where the collisions between grains and particles
of the interstitial (molecular) gas are taken into account. Thus, based on
previous studies Biben et al. (2002); Santos (2003), we discretize the
surrounding molecular gas, assigning individual particles with mass $m_{g}$
and diameter $\sigma_{g}$, thereby accounting for _elastic_ collisions between
grains and background gas particles in the starting kinetic equation.

Under USF conditions, the system is characterized by constant density profiles
$n$ and $n_{g}$, uniform temperatures $T$ and $T_{g}$, and a (common) flow
velocity $U_{x}=U_{g,x}=ay$, where $a$ denotes the shear rate. In agreement
with previous investigations on uniform sheared suspensions, the mean flow
velocity $\mathbf{U}$ is coupled to that of the gas phase $\mathbf{U}_{g}$.
Consequently, the viscous heating term due to shear and the energy transferred
by the grains from collisions with the molecular gas are compensated by the
cooling terms derived from collisional dissipation, allowing the achievement
of a steady state. A distinctive feature of the USF is that the one-particle
velocity distribution function $f(\mathbf{r},\mathbf{v})$ depends on space
only through its dependence on the peculiar velocity
$\mathbf{V}=\mathbf{v}-\mathbf{U}$. Consequently, the velocity distribution
function becomes uniform in the Lagrangian reference frame moving with
$\mathbf{V}$. This means that $f(\mathbf{r},\mathbf{v})\equiv f(\mathbf{V})$.
Based on symmetry considerations, the heat flux $\mathbf{q}$ vanishes, making
the pressure tensor $\mathrm{P}$ the relevant flux. Therefore, to understand
the intricate dynamics of granular suspensions under shear flow, it is
imperative to focus on their non-Newtonian properties —derived from the
pressure tensor. These include the (reduced) temperature $\chi=T/T_{g}$, the
(reduced) nonlinear shear viscosity $\eta^{*}$, and the (reduced) normal
stress difference $\Psi^{*}$.

Given that the most challenging aspect of dealing with the Boltzmann equation
lies in the collision operator, it is reasonable to explore alternatives that
render this operator more analytically tractable than in the case of IHS.
Among the most sophisticated techniques in this regard is to consider the
Boltzmann equation for IMM. As in the case of elastic collisions Ikenberry and
Truesdell (1956); Garzó and Santos (2003), the collision rate for IMM is
independent of the relative velocity of the colliding particles. As a
consequence, the collisional moments of degree $k$ of the Boltzmann
collisional operator can be expressed as a linear combination of velocity
moments of degree less than or equal to $k$. To complement the results derived
for IMM, we have also considered in this paper the use of a BGK-type kinetic
model where the true Boltzmann operator is replaced by a simple relaxation
term. Here, we employed both approaches to compute the rheological properties
of the sheared granular suspension. Thus, our objective was twofold. Firstly,
we aimed to assess the reliability and compatibility of the proposed models
with previous results Gómez González et al. (2023) obtained for IHS using
Grad’s moment method. Additionally, DSMC simulations for IHS were performed as
an alternative method to validate any potential discrepancies identified.
Secondly, taking advantage of the capabilities provided by the BGK model, we
endeavored to calculate the velocity distribution function and the higher-
order moments that offer insights into its characteristics.

Before proceeding with the computation of rheological properties, it is
necessary to understand the response of the molecular fluid to shear stress.
This assessment was also conducted using both IMM and BGK-type kinetic model
that were later used to model the granular gas. A novelty here is the
application of a force (Gaussian thermostat) of the form
$\mathbf{F}=-m\xi\mathbf{V}$ to compensate for the energy gained through
viscous shear stresses. This force, by consistency, also applies to the
granular gas, maintaining convergence to a steady state. As anticipated, the
results agree well with those obtained through Grad’s moment method Gómez
González et al. (2023) [see Eqs. (32) and (64)]. Consequently, once the
problem conditions (including the shear rate $a$) are defined, the molecular
temperature $T_{g}$ is determined, effectively serving as a thermostat for the
granular gas.

After accurately describing the rheology of the molecular gas, we focused on
modeling the granular gas. Using both IMM and BGK-type model separately, we
have calculated the non-zero elements of the the pressure tensor. The
knowledge of these elements allows us to identify the relevant rheological
properties of the granular suspension. As shown in Fig. 2, these quantities
are represented as functions of the coefficient of restitution $\alpha$ and
the mass ratio $m/m_{g}$. In particular, we find that the theoretical results
obtained from the Grad’s method for IHS, IMM and BGK-type model show
remarkable agreement, with almost indistinguishable curves. This underlines
the effectiveness of _structurally simple_ models in capturing the
complexities of sheared granular suspensions. We observe a DST-type transition
starting at a certain value of $a^{*}$, which increases with the mass ratio
$m/m_{g}$. Interestingly, similar to the MD simulations performed for IHS
Hayakawa et al. (2017), the DSMC data suggest a more abrupt transition than
predicted by theory. Given that the main divergences between Grad’s (for IHS)
and DSMC results arise from the form of the distribution function, the
significance of investigating higher-order moments to assess the deviation of
the distribution function from its Maxwellian reference is then justified.

Based on previous literature where discrepancies in fourth-degree moments have
been observed Gómez González and Garzó (2019a); Santos and Garzó (2007), and
acknowledging the potential lengthiness of calculations, for the sake of
simplicity, we decided to employ only the BGK-type model to compute higher-
order moments. Specifically, we concentrated on the (symmetric) fourth-degree
moments $M_{4|0}$ and $M_{2|xy}$. The shear-rate dependence of these moments
is illustrated in Fig. 3 for the same parameter values of $\alpha$ and
$m/m_{g}$ as those employed for the rheological quantities. Initially, we note
that the fourth-degree moments also exhibit an abrupt transition at a value of
$a^{*}$ that increases with $m/m_{g}$ until reaching the Brownian limit.
Furthermore, DSMC simulations in the Brownian limit qualitatively capture the
profile of these moments, although some quantitative disparities are apparent.
To ascertain the extent of these discrepancies, we narrowed our focus to the
interval $0<a^{*}<1$ where non-Newtonian effects are apparent. Additionally,
we included IMM results directly as obtained from an effective Fokker-Planck-
type model Gómez González and Garzó (2019b). Figure 4 illustrates that BGK
results overestimate the deviation from the moments computed when no shear
stress is applied compared to DSMC simulations and the results obtained using
an effective force to model the interstitial gas. These disparities are also
observed in the marginal distribution function $\phi_{x}$.

The theoretical findings presented here motivate the comparison with computer
simulations. Although the observed agreement in the Brownian limit is
encouraging, there is scope to extend this agreement to scenarios with finite
mass ratios. Our plan is to carry out simulations of this type in the near
future, which we expect will further validate and improve our theoretical
framework. In addition, we plan to extend our current findings to finite
densities by exploring the Enskog kinetic equation, which will allow us to
evaluate the involvement of density in the occurrence of these phenomena. This
line of research will be one of the main objectives of our upcoming research.

“Conceptualization, R. G. G. and V. G.; software, validation, R. G. G. and V.
G.; formal analysis, R. G. G. and V. G.; writing—original draft preparation,
V. G.; writing—review and editing, R. G. G. and V. G. All authors have read
and agreed to the published version of the manuscript.”

The authors 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 authors report no conflict of interest.

no

##  Appendix A Some technical details in the BGK-type kinetic model

In this Appendix we give some details on the determination of the parameters
$\widetilde{T}$ and $\nu$ appearing in the relaxation term (66). To get them,
we require the collisional transfer of energy of grains due to their elastic
collisions with particles of the molecular gas to be the same as the one
obtained from the true Boltzmann–Lorentz collision operator. This implies that

$\int d\mathbf{v}V^{2}J_{\text{BL}}[\mathbf{V}|f,f_{g}]=-\nu\int
d\mathbf{v}V^{2}\left[f(\mathbf{V})-\widetilde{f}_{g}(\mathbf{V})\right].$
(103)

Given that the collisional moment involving the operator
$J_{\text{BL}}[f,f_{g}]$ cannot be exactly computer for IHS, one estimate this
moment by replacing $f$ and $f_{g}$ by their Grad’s solutions Grad (1949). In
this approximation, one achieves the results Gómez González et al. (2023)

$\int
d\mathbf{v}V^{2}J_{\text{BL}}[\mathbf{V}|f,f_{g}]=-\frac{8\pi^{(d-1)/2}}{\Gamma\left(\frac{d}{2}\right)}nn_{g}\overline{\sigma}^{d-1}\frac{m_{g}}{(m+m_{g})^{2}}\left(\frac{2T_{g}}{m_{g}}+\frac{2T}{m}\right)^{1/2}\left(T-{T}_{g}\right).$
(104)

Moreover,

$\int
d\mathbf{v}V^{2}\left[f(\mathbf{V})-\widetilde{f}_{g}(\mathbf{V})\right]=-d\frac{n\nu}{m}\left(T-\widetilde{T}\right).$
(105)

From Eqs. (104) and (105) one gets the identity

$\frac{8\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n_{g}\overline{\sigma}^{d-1}\frac{mm_{g}}{(m+m_{g})^{2}}\left(\frac{2T_{g}}{m_{g}}+\frac{2T}{m}\right)^{1/2}\left(T-{T}_{g}\right)=\nu\left(T-\widetilde{T}\right).$
(106)

Equation (106) allows to make the identifications (70).

Finally, from Grad’s moment method Grad (1949), the collisional moment (71)
can be written as Gómez González et al. (2023)

$\int
d\mathbf{v}\;mV_{i}V_{j}J[f,f]=-\nu^{\text{M}}nT_{g}\left[\nu_{\eta}^{*\text{IHS}}P_{ij}^{*}-\chi\left(\nu_{\eta}^{*\text{IHS}}-\zeta^{*\text{IHS}}\right)\delta_{ij}\right],$
(107)

where $\nu^{\text{M}}$ is defined in Eq. (42) and

$\nu_{\eta}^{*\text{IHS}}=\frac{(1+\alpha)\left[3\left(1-\alpha\right)+2d\right]}{2d(d+2)},\quad\zeta^{*\text{IHS}}=\frac{1-\alpha^{2}}{2d}.$
(108)

The BGK-type kinetic model (65) yields the result

$\int
d\mathbf{v}\;mV_{i}V_{j}\Big{[}-\nu^{\prime}\left(f-f^{\text{M}}\right)+\frac{\epsilon}{2}\frac{\partial}{\partial\mathbf{V}}\cdot\mathbf{V}f\Big{]}=-nT_{g}\Big{[}\left(\nu^{\prime}+\epsilon\right)P_{ij}^{*}-\nu^{\prime}\chi\delta_{ij}\Big{]},$
(109)

where we recall that $P_{ij}^{*}=P_{ij}/nT_{g}$. Comparison between Eqs. (107)
and (109) yields Eq. (72) for $\nu^{\prime}(\alpha)$.

## References

  * Koch (1990) Koch, D.L.  Kinetic theory for a monodisperse gas-solid suspension.  Phys. Fluids A 1990, 2, 1711–1722.
  * Gidaspow (1994) Gidaspow, D.  Multiphase Flow and Fluidization; Academic Press, 1994.
  * Fullmer and Hrenya (2017) Fullmer, W.D.; Hrenya, C.M.  The clustering instability in rapid granular and gas-solid flows.  Annu. Rev. Fluid Mech. 2017, 49, 485–510.
  * Jackson (2000) Jackson, R.  The Dynamics of Fluidized Particles; Cambridge University Press, New York, 2000.
  * Tsao and Koch (1995) Tsao, H.K.; Koch, D.L.  Simple shear flows of dilute gas-solid suspensions.  J. Fluid Mech. 1995, 296, 211–245.
  * Sangani et al. (1996) Sangani, A.S.; Mo, G.; Tsao, H.K.; Koch, D.L.  Simple shear flows of dense gas-solid suspensions at finite Stokes numbers.  J. Fluid Mech. 1996, 313, 309–341.
  * Wylie et al. (2009) Wylie, J.J.; Zhang, Q.; Li, Y.; Hengyi, X.  Driven inelastic-particle systems with drag.  Phys. Rev. E 2009, 79, 031301.
  * Heussinger (2013) Heussinger, C.  Shear thickening in granular suspensions: Interparticle friction and dynamically correlated clusters.  Phys. Rev. E 2013, 88, 050201 (R).
  * Chamorro et al. (2015) Chamorro, M.G.; Vega Reyes, F.; Garzó, V.  Non-Newtonian hydrodynamics for a dilute granular suspension under uniform shear flow.  Phys. Rev. E 2015, 92, 052205.
  * Saha and Alam (2017) Saha, S.; Alam, M.  Revisiting ignited-quenched transition and the non-Newtonian rheology of a sheared dilute gas-solid suspension.  J. Fluid Mech. 2017, 833, 206–246.
  * Alam et al. (2019) Alam, M.; Saha, S.; Gupta, R.  Unified theory for a sheared gas-solid suspension: from rapid granular suspension to its small-Stokes-number limit.  J. Fluid Mech. 2019, 870, 1175–1193.
  * Saha and Alam (2020) Saha, S.; Alam, M.  Burnett-order constitutive relations, second moment anisotropy and co-existing states in sheared dense gas-solid suspensions.  J. Fluid Mech. 2020, 887, A9.
  * Garzó et al. (2012) Garzó, V.; Tenneti, S.; Subramaniam, S.; Hrenya, C.M.  Enskog kinetic theory for monodisperse gas-solid flows.  J. Fluid Mech. 2012, 712, 129–168.
  * Hayakawa et al. (2017) Hayakawa, H.; Takada, S.; Garzó, V.  Kinetic theory of shear thickening for a moderately dense gas-solid suspension: From discontinuous thickening to continuous thickening.  Phys. Rev. E 2017, 96, 042903.
  * Hayakawa and Takada (2019) Hayakawa, H.; Takada, S.  Kinetic theory of discontinuous rheological phase transition for a dilute inertial suspension.  Prog. Theor. Exp. Phys 2019, 083J01.
  * Gómez González and Garzó (2019a) Gómez González, R.; Garzó, V.  Simple shear flow in granular suspensiones: Inelastic Maxwell models and BGK-type kinetic model.  J. Stat. Mech. 2019, 013206.
  * Gómez González and Garzó (2019b) Gómez González, R.; Garzó, V.  Transport coefficients for granular suspensions at moderate densities.  J. Stat. Mech. 2019, 093204.
  * Hayakawa et al. (2020) Hayakawa, H.; Takada, S.; Garzó, V.  Erratum: Kinetic theory of shear thickening for a moderately dense gas-solid suspension: from discontinuous thickening to continuous thickening.  Phys. Rev. E 2020, 101, 069904(E).
  * Takada et al. (2020) Takada, S.; Hayakawa, H.; Santos, A.; Garzó, V.  Enskog kinetic theory of rheology for a moderately dense inertial suspension.  Phys. Rev. E 2020, 102, 022907.
  * Gómez González et al. (2020) Gómez González, R.; Khalil, N.; Garzó, V.  Enskog kinetic theory for multicomponent granular suspensions.  Phys. Rev. E 2020, 101, 012904.
  * Takada et al. (2023) Takada, S.; Hayakawa, H.; Garzó, V.  Rheology of a dilute binary mixture of inertial suspension under simple shear flow.  Prog. Theor. Exp. Phys. 2023, 113J01.
  * Baldassarri et al. (2002) Baldassarri, A.; Marconi, U.M.B.; Puglisi, A.  Influence of correlations of the velocity statistics of scalar granular gases.  Europhys. Lett. 2002, 58, 14–20.
  * Gómez González and Garzó (2022) Gómez González, R.; Garzó, V.  Kinetic theory of granular particles immersed in a molecular gas.  J. Fluid Mech. 2022, 943, A9.
  * Subramaniam (2020) Subramaniam, S.  Multiphase flows: Rich physics, challenging theory, and big simulations.  Phys. Rev. Fluids 2020, 5, 110520.
  * Gómez González et al. (2023) Gómez González, R.; García Chamorro, M.; Garzó, V.  Rheology of granular particles immersed in a molecular gas under uniform shear flow.  arXiv:2311.16717 2023.
  * Grad (1949) Grad, H.  On the kinetic theory of rarefied gases.  Commun. Pure Appl. Math. 1949, 2, 331–407.
  * Chapman and Cowling (1970) Chapman, S.; Cowling, T.G.  The Mathematical Theory of Nonuniform Gases; Cambridge University Press, Cambridge, 1970.
  * Garzó and Santos (2003) Garzó, V.; Santos, A.  Kinetic Theory of Gases in Shear Flows. Nonlinear Transport; Springer, Netherlands, 2003.
  * Ben-Naim and Krapivsky (2003) Ben-Naim, E.; Krapivsky, P.L.  The Inelastic Maxwell Model. In Granular Gas Dynamics; Pöschel, T.; Luding, S., Eds.; Springer, 2003; Vol. 624, Lectures Notes in Physics, pp. 65–94.
  * Garzó and Santos (2007) Garzó, V.; Santos, A.  Third and fourth degree collisional moments for inelastic Maxwell models.  J. Phys. A: Math. Theor. 2007, 40, 14927–14943.
  * Sánchez Romero and Garzó (2023) Sánchez Romero, C.; Garzó, V.  High-degree collisional moments of inelastic Maxwell mixtures-Application to the homogeneous cooling and uniform shear flow states.  Entropy 2023, 25, 222.
  * Vega Reyes et al. (2007) Vega Reyes, F.; Garzó, V.; Santos, A.  Granular mixtures modeled as elastic hard spheres subject to a drag force.  Phys. Rev. E 2007, 75, 061306.
  * Cercignani (1988) Cercignani, C.  The Boltzmann Equation and Its Applications; Springer–Verlag, New York, 1988.
  * Pelargonio and Zaccone (2023) Pelargonio, S.; Zaccone, A.  Generalized Langevin equation with shear flow and its fluctuation-dissipation theorems derived from a Caldeira-Leggett Hamiltonian.  Phys. Rev. E 2023, 107, 064102.
  * Bird (1994) Bird, G.A.  Molecular Gas Dynamics and the Direct Simulation Monte Carlo of Gas Flows; Clarendon, Oxford, 1994.
  * Evans and Morriss (1990) Evans, D.J.; Morriss, G.P.  Statistical Mechanics of Nonequilibrium Liquids; Academic Press, London, 1990.
  * Brilliantov and Pöschel (2004) Brilliantov, N.; Pöschel, T.  Kinetic Theory of Granular Gases; Oxford University Press, Oxford, 2004.
  * Garzó (2019) Garzó, V.  Granular Gaseous Flows; Springer Nature, Cham, 2019.
  * Santos and Garzó (2007) Santos, A.; Garzó, V.  Simple shear flow in inelastic Maxwell models.  J. Stat. Mech. 2007, P08021.
  * Garzó (2003) Garzó, V.  Nonlinear transport in inelastic Maxwell mixtures under simple shear flow.  J. Stat. Phys. 2003, 112, 657–683.
  * Ikenberry and Truesdell (1956) Ikenberry, E.; Truesdell, C.  On the pressures and the flux of energy in a gas according to Maxwell’s kinetic theory.  J. Rat. Mech. Anal. 1956, 5, 1–54.
  * Dufty et al. (1986) Dufty, J.W.; Santos, A.; Brey, J.J.; Rodríguez, R.F.  Model for nonequilibrium computer simulation methods.  Phys. Rev. A 1986, 33, 459–466.
  * Brey et al. (1999) Brey, J.J.; Dufty, J.W.; Santos, A.  Kinetic models for granular flow.  J. Stat. Phys. 1999, 97, 281–322.
  * Bhatnagar et al. (1954) Bhatnagar, P.L.; Gross, E.P.; Krook, M.  A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component system.  Phys. Rev. 1954, 94, 511–525.
  * Gross and Krook (1956) Gross, E.P.; Krook, M.  Model for collision processes in gases. Small amplitude oscillations of charged two-component systems.  Phys. Rev. 1956, 102, 593–604.
  * Garzó and Santos (1995) Garzó, V.; Santos, A.  Comparison between the Boltzmann and BGK equations for uniform shear flow.  Physica A 1995, 213, 426–434.
  * Parmentier and Simonin (2012) Parmentier, J.F.; Simonin, O.  Transition models from the quenched to ignited states for flows of inertial particles suspended in a simple sheared viscous fluid.  J. Fluid Mech. 2012, 711, 147–160.
  * Wang et al. (2014) Wang, T.; Grob, M.; Zippelius, A.; Sperl, M.  Active microrhelogy of driven granular particles.  Phys. Rev. E 2014, 89, 042209.
  * Puglisi et al. (1998) Puglisi, A.; Loreto, V.; Marconi, U.M.B.; Petri, A.; Vulpiani, A.  Clustering and Non-Gaussian behavior in granular matter.  Phys. Rev. Lett. 1998, 81, 3848–3851.
  * Biben et al. (2002) Biben, T.; Martin, P.; Piasecki, J.  Stationary state of thermostated inelastic hard spheres.  Physica A 2002, 310, 308–324.
  * Santos (2003) Santos, A.  Granular fluid thermostated by a bath of elastic hard spheres.  Phys. Rev. E 2003, 67, 051101.