Mobility and diffusion of intruders in granular suspensions: Einstein relation

An interesting and challenging problem in statistical mechanics is the generalization of the fluctuation-response relation to non-equilibrium situations [1]. This problem has received considerable attention in the past few years by many researchers, who have tried to continue the research not only by means of theoretical tools but also by employing computer simulations. Among the different systems that are inherently out of equilibrium, granular matter can be considered as a good candidate to analyze this problem. It is well known that when granular matter is externally excited (rapid flow conditions), the motion of grains resembles the chaotic motion of atoms or molecules in a conventional molecular fluid. However, given that the size of the grains is mesoscopic (of the order of 1 µm, for instance), their interactions are inelastic and hence the total energy of the system decreases with time. To keep it in rapid conditions, one has to inject energy into the system to compensate for the energy dissipated by collisions, and hence a non-equilibrium steady state (NESS) is reached. For this sort of system, the fluctuation-response theorem has been proposed in terms of an effective temperature, which is clearly different from the environmental temperature [2–5].

It is quite apparent that an analysis of the validity of the fluctuation-response theorem requires knowledge of the complete dependence of the response and correlation functions on the frequency ω [6]. Since this is in fact quite a difficult problem, in order to gain some insight into the general problem, one usually considers the limit of small frequencies (ω → 0). In this limiting case, the classical relation between the diffusion coefficient D (autocorrelation function) and the mobility coefficient λ (linear response) is known as the Einstein relation [6].

In the case of granular gases, fluctuation response relations have been derived [7, 8] with respect to the so-called homogeneous cooling state (i.e. a state whose fate is a thermal death). In this situation, it has been proven that the response to an external force on an intruder (or impurity) particle violates the usual Einstein relation between the diffusion and mobility coefficients. There are three distinct origins for the violation of the Einstein relation: the deviation of the homogeneous cooling state from the Gibbs state (non-Gaussian distribution functions for the intruder and particles of the granular gas), the cooling of the reference state (yielding a different time dependence for D and λ), and energy non-equipartition (leading to different kinetic temperatures between the intruder and gas particles). A different approach widely employed in kinetic theory and computer simulations consists of considering driven granular gases where the system is heated by an external force (or thermostat) that compensates for the energy lost by collisions. This was the situation studied in [3] by computer simulations; it demonstrated the validity of the Einstein relation in NESS when the temperature of the bath was replaced by the temperature of the intruder T₀ ($\epsilon_0\equiv D/T_0 \lambda = 1$). This conclusion also agrees with the results derived from an exactly solvable model for driven dissipative systems [4].

Needless to say, thermostats are introduced to mimic the effects produced by bulk driving as in air-fluidized beds, for instance [9, 10]. Unfortunately, in most cases, the relationship between the results derived in driven (thermostated) granular gases and those obtained in real experiments is not clear. A more realistic example of thermostated granular systems consists of a set of solid particles surrounded by a gas of molecular particles. This provides a suitable starting point to model the behavior of granular suspensions. When the dynamics of grains are essentially ruled by their collisions, the tools of the classical kinetic theory (conveniently adapted to inelastic collisions) can be a reliable way to describe this type of granular flow [11]. However, due to the technical difficulties embodied in the study of two or more phases, a coarse-grained approach is generally adopted. In this description, the influence of the background gas on grains is usually incorporated in the kinetic equation through a fluid–solid interaction force. Usually, the gas-phase effects on the solid particles are described by the addition of a Fokker–Planck term (drag force term plus stochastic Langevin-like term) in the kinetic equation [12]. In fact, this way of driving the granular gas has been employed in computer simulations [3]. It is important to stress that this suspension model can also be derived in a more rigorous way by explicit consideration of the (elastic) collisions between grains and particles of the molecular gas. In this discrete description, the above collisions are accounted for via the Boltzmann–Lorentz collision operator [13]. In the limit where the grains are much heavier than the molecular gas particles, the Boltzmann–Lorentz operator reduces to the Fokker–Planck operator, and the results for the transport properties derived from the collision model [14] agree with those obtained from the coarse-grained approach [15].

The objective of this paper is to determine the diffusion D and mobility λ transport coefficients in a granular suspension. For moderate densities, our starting point is the (nonlinear) Enskog and the (linear) Enskog–Lorentz kinetic equations for the granular gas and the intruders, respectively, with the addition of Fokker–Planck operators to each one of these kinetic equations. The interaction between the grains and intruders with the interstitial gas is through two different drift coefficients γ and γ₀, respectively. To first order in both the concentration gradient and the external field, the Enskog–Lorentz equation is solved by means of the Chapman–Enskog method [16] adapted to dissipative dynamics. As for elastic collisions, the coefficients D and λ are given in terms of a set of coupled linear integral equations that are approximately solved by considering the second Sonine approximation (i.e. the second-order truncation of the Sonine polynomial expansion of the velocity distribution of intruders). As occurs in driven granular gases [17, 18], our results show that the deviations of the modified Einstein relation ε₀ from 1 are only due to the very small departure of the reference state (zeroth-order distribution of intruders) from the Maxwell–Boltzmann distribution. This departure is measured by the kurtosis (or fourth-degree cumulant) c₀. Since in general the magnitude of c₀ in granular suspensions is much smaller than the one obtained in freely cooling systems [7] and/or in driven granular gases [17, 18], one may conclude that the verification of the modified Einstein relation is much more accurate in gas–solid flows than in dry granular gases. This is likely one of the most relevant conclusions of the present work.

It is important to remark that our results are based on the Enskog equation. This equation is an extension of the Boltzmann equation (which holds for very dilute gases) to moderate densities. In this regime of densities, although spatial correlations are accounted for via the pair correlation function in this kinetic equation, velocity correlations between the particles that are about to collide are neglected (molecular chaos assumption) as in the Boltzmann description. This is the main limitation of the Enskog equation. In this particular context, it is worth highlighting the computer simulation results obtained by Puglisi et al [19]. These results demonstrate that the departure from the Einstein relation primarily arises from spatial and velocity correlations that emerge with increasing density, rather than from the non-Gaussian corrections to the distribution function. This conclusion has also been confirmed by experimental evidence [20] involving the Brownian motion of a rotating intruder immersed in a vibro-fluidized granular medium. It is shown that Einstein's relation holds in the dilute regime while it is violated for a high packing fraction; this violation cannot be explained in terms of effective temperatures. On the other hand, given that the spatial and velocity correlations are present for densities and inelasticities at which the Enskog equation does not presumably apply, the conclusions reached in [19] and [20] are not in conflict with those derived here.

The plan of the paper is as follows. In section 2, we describe the problem we are interested in. The steady homogeneous state of the intruders plus granular gas in contact with a thermal bath is studied in section 3. As expected, the intruder's temperature T₀ differs from that of the granular gas T, and so there is a breakdown of energy equipartition. In section 4, the Chapman–Enskog method is applied to solve the Enskog–Lorentz kinetic equation to first order in the concentration gradient and the external field. Some technical details concerning the calculations of the paper are provided in the appendices A and B. The theoretical results for D and λ are compared with Monte Carlo simulation results showing an excellent agreement, especially in the case of the second Sonine solution. The knowledge of T, T₀, D, and λ allows us to compute the conventional ε (defined in terms of the bath temperature T_b) and modified ε₀ (defined in terms of the intruder temperature T₀) Einstein relations. While $\epsilon_0\simeq 1$, ε clearly differs from 1, showing that the violation of the conventional Einstein relation in granular suspensions can be significant. We close the paper with some concluding remarks.

2.1. Granular gas

We consider a granular gas of inelastic hard spheres of mass m and diameter σ. The solid particles are immersed in a gas of viscosity η_g . Spheres (grains) are assumed to be completely smooth so that the inelasticity of collisions is only characterized by the constant (positive) coefficient of normal restitution $\alpha\unicode{x2A7D} 1$. When the suspensions are dominated by collisions (which are assumed to be nearly instantaneous) [11], a coarse-grained description can be adopted to account for the influence of the gas on the dynamics of solid particles. In this approach, the effect of the gas phase on grains is usually incorporated in the starting kinetic equation by means of a fluid–solid interaction force [21–23]. Some models for granular suspensions [24–33] only take into account the Stokes linear drag force law (which attempts to mimic the friction of grains with the interstitial gas) for gas–solid interactions. On the other hand, some works [34] have shown that the drag force term does not correctly capture the particle acceleration–velocity correlation observed in direct numerical simulations [35]. For this reason, an additional Langevin-like term is included in the effective fluid–solid force. This stochastic term models the additional effects of neighboring particles via the stochastic increment of a Wiener process [12]. In addition, this term (which randomly kicks the particles between collisions) also takes into account the energy gained by the solid particles due to their interaction with the background gas.

Thus, according to the above coarse-grained description, for moderate densities and assuming that the granular gas is in a steady homogeneous state, the one-particle velocity distribution function $f(\mathbf{v}, t)$ of the granular gas verifies the nonlinear Enskog equation [36]:

$$\begin{align} -\gamma\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}f-\frac{\gamma T_{\text{b}}}{m}\frac{\partial^2 f}{\partial v^2} = J\left[\mathbf{v}|f,f\right], \end{align} \tag{ 1 }$$

where the Enskog collision operator is

$$\begin{align} J\left[\mathbf{v}_1|f,f\right] & = \chi \sigma^{d-1}\int \text{d}\mathbf{v}_2\int \text{d}\widehat{\boldsymbol{\sigma}}\Theta\left(\widehat{\boldsymbol{\sigma}}\cdot\mathbf{g}_{12}\right) \left(\widehat{\boldsymbol{\sigma}}\cdot\mathbf{g}_{12}\right)\notag\\ & \quad\times\left[\alpha^{-2}f\left(\mathbf{v}_1^{^{\prime\prime}},t\right)f\left(\mathbf{v}_2^{^{\prime\prime}},t\right)-f\left(\mathbf{v}_1,t\right)f\left(\mathbf{v}_2,t\right)\right]. \end{align} \tag{ 2 }$$

Here, χ is the pair correlation function for grain–grain collisions at contact (i.e. when the distance between their centers is σ), $d$ is the dimensionality of the system, $\widehat{\boldsymbol{\sigma}}$ is a unit vector directed along the line of centers of the colliding spheres, Θ is the Heaviside step function [$\Theta(x) = 1$ for x > 0, $\Theta(x) = 0$ for $x\unicode{x2A7D} 0$], and $\mathbf{g}_{12} = \mathbf{v}_1-\mathbf{v}_2$ is the relative velocity of the two colliding spheres. The double primes on the velocities denote the initial values $(\mathbf{v}_1^{^{\prime\prime}},\mathbf{v}_2^{^{\prime\prime}})$ that yield $(\mathbf{v}_1,\mathbf{v}_2)$ following a binary collision:

$$\begin{align} \mathbf{v}_1^{^{\prime\prime}} = \mathbf{v}_1-\frac{1+\alpha^{-1}}{2}\left(\boldsymbol{\widehat{\sigma}} \cdot\mathbf{g}_{12}\right)\boldsymbol{\widehat{\sigma}}, \quad \mathbf{v}_2^{^{\prime}} = \mathbf{v}_2+\frac{1+\alpha^{-1}}{2}\left(\boldsymbol{\widehat{\sigma}} \cdot\mathbf{g}_{12}\right)\boldsymbol{\widehat{\sigma}} .\end{align} \tag{ 3 }$$

In equation (1), γ is the drift or friction coefficient (characterizing the interaction between particles of the granular gas and the background gas) and T_b is the bath temperature. As in previous works [15, 37], we assume here that γ is a scalar quantity proportional to the gas viscosity [38]. In the dilute limit every particle is only subjected to its respective Stokes drag and so, for hard spheres (d = 3), the drift coefficient γ is defined as

$$\begin{align} \gamma\equiv \gamma_\text{St} = \frac{3 \pi \sigma \eta_g}{m}. \end{align} \tag{ 4 }$$

Beyond the dilute limit, for moderate densities and low Reynolds numbers, one has the relationship

$$\begin{align} \gamma = \gamma_\text{St}R\left(\phi\right), \end{align} \tag{ 5 }$$

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

$$\begin{align} \phi = \frac{\pi^{d/2}}{2^{d-1}d\Gamma \left(\frac{d}{2}\right)}n\sigma^d. \end{align} \tag{ 6 }$$

The density dependence of the dimensionless function R can be inferred from computer simulations. Specific forms of R will be chosen later to assess the dependence of the dynamic properties of the system on the parameter space of the problem. On the other hand, it is worthwhile remarking that the results reported in this paper apply regardless of the specific choice of the function R.

In the homogeneous state, the only nontrivial balance equation is that of the granular temperature T, defined as

$$\begin{align} {d} n T = \int \textrm{d}\mathbf{v}\; m v^2 f\left(\mathbf{v}\right), \end{align} \tag{ 7 }$$

where

$$\begin{align} n = \int \textrm{d}\mathbf{v}\; f\left(\mathbf{v}\right) \end{align} \tag{ 8 }$$

is the number density of solid particles. The balance equation for T can be easily derived by multiplying both sides of equation (1) by $m v^2$ and integrating over velocity. This is given by

$$\begin{align} 2\gamma \left(T_\text{b}-T \right) = T \zeta, \end{align} \tag{ 9 }$$

where

$$\begin{align} \zeta = -\frac{1}{{d} n T}\int \textrm{d}\mathbf{v}\; m v^2\; J\left[f,f\right] \end{align} \tag{ 10 }$$

is the cooling rate. This quantity gives the rate of change of energy dissipated by collisions. When collisions are elastic (α = 1), ζ = 0. Since ζ is a functional of the distribution $f(\mathbf{v})$, it is quite obvious that one needs to know f to determine the cooling rate.

In the case of elastic collisions, equation (9) leads to the result $T = T_\text{b}$, and the Enskog equation (1) admits the simple Maxwell–Boltzmann solution

$$\begin{align} f\left(\mathbf{v}\right) = f_\text{b,M}\left(\mathbf{v}\right) = n\left(\frac{m}{2\pi T_\text{b}}\right)^{d/2}\text{exp}\left(-\frac{m v^2}{2 T_\text{b}}\right). \end{align} \tag{ 11 }$$

This result is nothing more than a consequence of the fluctuation-dissipation theorem [39]. On the other hand, for inelastic collisions (α ≠ 1), the exact solution of equation (1) is not known. However, in the region of thermal velocities, a good approximation can be obtained from an expansion in Sonine polynomials. In the leading order, the distribution f can be written as

$$\begin{align} f\left(\mathbf{v}\right) \to n \pi^{-d/2} v_\text{th}^{-d} e^{-\xi^2}\left\{1+\frac{c}{2}\left[\xi^4-\left(d+2\right)\xi^2+\frac{d\left(d+2\right)}{4}\right]\right\}, \end{align} \tag{ 12 }$$

where $\boldsymbol{\xi} = \mathbf{v}/v_\text{th}$ and $v_\text{th} = \sqrt{2T/m}$ is a thermal speed. The coefficient c (which measures the deviation of f from its Maxwellian form) is related to the kurtosis of the distribution. Its value has been estimated from the Enskog equation by considering linear terms in c [15]. Its explicit expression is

$$\begin{align} c = \frac{16\left(1-\alpha\right)\left(1-2\alpha^2\right)}{73+56d-3\alpha\left(35+8d\right)+30\left(1-\alpha\right)\alpha^2+\frac{64d\left(d+2\right)}{1+\alpha}\gamma^*} \ , \end{align} \tag{ 13 }$$

where

$$\begin{align} \gamma^* = \frac{\gamma}{\nu} = \frac{\sqrt{\pi}}{2^{d} d}\frac{R}{\phi\chi\sqrt{T^*}}. \end{align} \tag{ 14 }$$

Here, $T^* = T/\mathcal{T}$, $\mathcal{T} = m\sigma^2 \gamma_\text{St}^2$, and we have introduced the effective collision frequency

$$\begin{align} \nu = \frac{\sqrt{2}\pi^{\left(d-1\right)/2}}{\Gamma\left(\frac{d}{2}\right)}n\sigma^{d-1}\chi v_\text{th}. \end{align} \tag{ 15 }$$

The cooling rate ζ can be also determined from the Sonine approximation (12) with the result

$$\begin{align} \zeta = \frac{1-\alpha^2}{d}\left(1+\frac{3}{16}c\right)\nu. \end{align} \tag{ 16 }$$

Upon obtaining equation (16), nonlinear terms in c have been neglected.

For practical purposes, it is convenient to write equation (9) in dimensionless form. In this case, one achieves the equation

$$\begin{align} 2\delta \left(T_b^*-T^*\right) = \zeta^* T^{*3/2}, \end{align} \tag{ 17 }$$

where $T_\text{b}^* = T_\text{b}/\mathcal{T}$, $\zeta^* = \zeta/\nu$, and

$$\begin{align} \delta = \frac{\sqrt{\pi}}{2^{d} d}\frac{R}{\phi\chi} .\end{align} \tag{ 18 }$$

If one neglects the kurtosis c (which is in general very small [15]) in the expression (16) of $\zeta^*$, then equation (17) becomes a cubic equation for the (reduced) temperature $T^*$. In terms of the auxiliary parameter $\varepsilon\equiv \zeta^* \sqrt{T_\text{b}^*}/(2\delta)$, the physical (real) root of the cubic equation (17) can be written as

$$\begin{align} T^* = \frac{\left(\Xi^{1/3}+\Xi^{-1/3}-1\right)^2}{9\varepsilon^2}T_\text{b}^*, \end{align} \tag{ 19 }$$

where

$$\begin{align} \Xi = \frac{3\sqrt{3}\sqrt{27\varepsilon^4-4\varepsilon^2}+27\varepsilon^2-2}{2}. \end{align} \tag{ 20 }$$

As expected, for elastic collisions (α = 1), ε → 0 and so $T^* = T_\text{b}^*$ for any value of φ and $T_\text{b}^*$. When α < 1, $T^*\lt T_\text{b}^*$ since the granular temperature is smaller than that of the background gas. In the case that the coefficient c is not neglected, equation (17) is a quartic equation whose physical solution must be numerically determined.

The theoretical results for the (reduced) temperature $T^*$ are compared against direct simulation Monte Carlo (DSMC) simulations [40]. The DSMC simulations are performed following the same steps as described in [41] and [42]. Notably, modifications have been introduced to the collision stage of the simulation algorithm originally employed by Montanero and Garzó [43], with the primary objective of incorporating two key considerations: (i) the tracer concentration of intruder particles and (ii) the impact of the interstitial gas on the dynamic behavior of solid particles. The former adjustment entails the exclusion of intruder–intruder collisions and the preservation of grain velocities after grain–intruder collisions. In contrast, the latter modification exhibits a higher degree of complexity. For a three-dimensional system (d = 3), the influence of the interstitial fluid on grains is incorporated by iteratively updating the velocity vector v_k of each individual grain belonging to species i after every time increment δt, in accordance with the rule [44]:

$$\begin{align} \mathbf{v}_k\to e^{-\gamma_i\delta t}\mathbf{v}_k+\left(\frac{6\gamma_i T_\text{b} \delta t}{m_i}\right)^{1/2}\boldsymbol{\varpi}_k. \end{align} \tag{ 21 }$$

Here, $\boldsymbol{\varpi}_k$ is a random vector of zero mean and unit variance. Equation (21) converges to the Fokker–Plank operator when the time step δt is much shorter than the mean free time between collisions [44].

Figure 1 shows the dependence of the (reduced) granular temperature $T^*$ on the coefficient of restitution α for a three-dimensional (d = 3) granular gas with $T_\text{b}^* = 1$ and φ = 0.1. For d = 3, a good approximation to the pair correlation function χ is [45]

$$\begin{align} \chi = \frac{1-\frac{1}{2}\phi}{\left(1-\phi\right)^3}. \end{align} \tag{ 22 }$$

Moreover, for the sake of illustration, simulations for hard sphere systems [46–48] suggest the following form for the function $R(\phi)$:

$$\begin{align} R\left(\phi\right) = \frac{10 \phi}{\left(1-\phi\right)}+\left(1-\phi\right)^3\left(1+1.5 \sqrt{\phi}\right). \end{align} \tag{ 23 }$$

The line in figure 1 corresponds to the numerical solution to equation (17) while the symbols refer to the numerical results obtained from DSMC simulations [40]. As expected, the energy dissipated by collisions increases with increasing inelasticity, and so the kinetic energy of grains (or equivalently, their reduced temperature $T^*$) decreases. We also observe an excellent agreement between theory and simulations in figure 1 in the complete range of values α. Although not plotted, the curve given by the exact solution (19) (i.e. when one neglects c) is indistinguishable from the one represented in figure 1 taking into account the value of c. This feature can be easily explained by figure 2 where the α dependence of the fourth cumulant c is plotted for the same system as that of figure 1. We find that the magnitude of c is very small — much smaller than in the case of (dry) granular gases [49]. Moreover, the cumulant c exhibits a non-monotonic dependence on α since it decreases first as increasing inelasticity, reaches a minimum, and then increases with decreasing α. As in the case of $T^*$, an excellent agreement between theory and simulations is found.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.2. Intruders immersed in a granular gas

We now assume that a few intruders (of mass m₀ and diameter σ₀) are added to the system. Since the concentration of intruders is negligibly small, one can assume that the state of the granular gas is not disturbed by the presence of the intruders, and hence its distribution function $f(\mathbf{v})$ obeys the Enskog equation (1). Moreover, one can also neglect collisions among intruders themselves in the kinetic equation of the one-particle velocity distribution function $f_0(\mathbf{r}, \mathbf{v};t)$ of intruders. Thus, in this limiting tracer case, only the intruder–granular gas collisions (which are characterized by the coefficient of restitution $\alpha_0\neq \alpha$) will be considered in the above kinetic equation. Intruders also interact with the interstitial fluid through the friction coefficient γ₀, which is in general different from γ. Since we are also interested in obtaining the mobility of intruders, we will assume that intruder particles are also subjected to the action of a weak external field E (e.g. gravity or an electric field). This field only acts on intruders.

Note that, formally, the system (intruder plus granular gas) can be regarded as a binary granular suspension where one of the species is present in tracer concentration. For conciseness, in the remainder of the paper we will refer to intruders immersed in a granular suspension instead of a binary granular suspension with one tracer species.

Under the above conditions, the one-particle velocity distribution function $f_0(\mathbf{r}, \mathbf{v};t)$ of intruders verifies the Enskog–Lorentz kinetic equation

$$\begin{align} \frac{\partial f_0}{\partial t}+\mathbf{v}\cdot \nabla f_0+\frac{\mathbf{E}}{m_0}\cdot \frac{\partial}{\partial\mathbf{v}}\cdot f_0 -\gamma_0\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}f_0-\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2 f_0}{\partial v^2} = J_0\left[\mathbf{r}, \mathbf{v}|f_0,f\right], \end{align} \tag{ 24 }$$

where the Enskog–Lorentz collision operator $J_0[f_0,f]$ is [36]

$$\begin{align} J_{0}\left[\mathbf{r}_1, \mathbf{v}_1|f_0,f\right] & = \overline{\sigma}^{d-1}\int \text{d}\mathbf{v}_2\int \text{d}\widehat{\boldsymbol{\sigma}}\Theta\left(\widehat{\boldsymbol{\sigma}}\cdot\mathbf{g}_{12}\right) \left(\widehat{\boldsymbol{\sigma}}\cdot\mathbf{g}_{12}\right)\notag\\ & \quad \times \left[\alpha_{0}^{-2}\chi_{0} \left(\mathbf{r}_1,\mathbf{r}_1-\boldsymbol{\overline{\sigma}}\right)f_0\left(\mathbf{r}_1,\mathbf{v}_1^{^{\prime\prime}},t\right) f\left(\mathbf{v}_2^{^{\prime\prime}}\right)\right.\notag\\ & \quad \left.-\chi_{0}\left(\mathbf{r}_1,\mathbf{r}_1+\boldsymbol{\overline{\sigma}}\right)f_0\left(\mathbf{r}_1,\mathbf{v}_1,t\right) f\left(\mathbf{v}_2\right)\right]. \end{align} \tag{ 25 }$$

Here, χ₀ is the pair correlation function for intruder–granular gas collisions, $\boldsymbol{\overline{\sigma}} = \overline{\sigma} \widehat{\boldsymbol{\sigma}}$, $\overline{\sigma} = (\sigma+\sigma_0)/2$, and $\widehat{\boldsymbol{\sigma}}$ is the unit vector directed along the line of centers from the sphere of the intruder to the sphere of the granular gas at contact. The relationship between the velocities $(\mathbf{v}_1^{^{\prime\prime}},\mathbf{v}_2^{^{\prime\prime}})$ and $(\mathbf{v}_1,\mathbf{v}_2)$ is

$$\begin{align} \mathbf{v}_1^{^{\prime\prime}} = \mathbf{v}_1-\mu \left(1+\alpha_0^{-1}\right)\left(\boldsymbol{\widehat{\sigma}} \cdot\mathbf{g}_{12}\right)\boldsymbol{\widehat{\sigma}}, \quad \mathbf{v}_2^{^{\prime\prime}} = \mathbf{v}_2+\mu_0 \left(1+\alpha_0^{-1}\right)\left(\boldsymbol{\widehat{\sigma}} \cdot\mathbf{g}_{12}\right)\boldsymbol{\widehat{\sigma}}, \end{align} \tag{ 26 }$$

where

$$\begin{align} \mu = \frac{m}{m+m_0}, \quad \mu_0 = \frac{m_0}{m+m_0}. \end{align} \tag{ 27 }$$

Equation (26) gives the so-called inverse or restituting collisions. The so-called direct collisions are defined as collisions where the pre-collisional velocities $(\mathbf{v}_1,\mathbf{v}_2)$ yield $(\mathbf{v}_1^{^{\prime}},\mathbf{v}_2^{^{\prime}})$ as post-collisional velocities. Inversion of the collision rules (26) gives the forms

$$\begin{align} \mathbf{v}_1^{^{\prime}} = \mathbf{v}_1-\mu \left(1+\alpha_0\right)\left(\boldsymbol{\widehat{\sigma}} \cdot\mathbf{g}_{12}\right)\boldsymbol{\widehat{\sigma}}, \quad \mathbf{v}_2^{^{\prime}} = \mathbf{v}_2+\mu_0 \left(1+\alpha_0\right)\left(\boldsymbol{\widehat{\sigma}} \cdot\mathbf{g}_{12}\right)\boldsymbol{\widehat{\sigma}}. \end{align} \tag{ 28 }$$

Moreover, note that upon writing equation (25) we have accounted for the granular gas being in a steady homogeneous state.

In accordance with equation (5), the friction coefficient γ₀ for the intruder can be written as

$$\begin{align} \gamma_0 = \gamma_{0,\text{St}}R_0, \end{align} \tag{ 29 }$$

where for d = 3,

$$\begin{align} \gamma_{0,\text{St}} = \frac{3 \pi \sigma_0 \eta_g}{m_0} = \frac{\sigma_0 m}{\sigma m_0}\gamma_{\text{St}}. \end{align} \tag{ 30 }$$

As in the case of $R(\phi)$, the dependence of the function R₀ on the density φ and the remaining parameters of the system will be taken from the results obtained by computer simulations.

Apart from the granular temperature T, it is convenient at a kinetic level to introduce the partial temperature of intruders T₀. This quantity measures the mean kinetic energy of intruders. It is defined as

$$\begin{align} {d} n_0(\mathbf{r};t) T_0(\mathbf{r};t) = \int \textrm{d}\mathbf{v}\; m_0 v^2 f_0(\mathbf{r}, \mathbf{v}; t), \end{align} \tag{ 31 }$$

where

$$\begin{align} n_0\left(\mathbf{r};t\right) = \int \textrm{d}\mathbf{v}\; f_0\left(\mathbf{r}, \mathbf{v}; t\right) \end{align} \tag{ 32 }$$

is the number density of intruders. Upon writing equation (31) we have taken into account that the mean flow velocity of the granular gas vanishes in our problem. It must be recalled that n₀ is much smaller than its counterpart n for the particles of the granular gas.

Before considering the diffusion of intruders due to the presence of a weak concentration gradient $\nabla n_0$ and/or a weak external field E, it is convenient to first characterize the homogeneous steady state of intruders. This is a crucial point since the latter state plays the role of the reference state in the Chapman–Enskog solution to equation (24).

In the absence of diffusion (homogeneous steady state), equation (24) becomes

$$\begin{align} -\gamma_0\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}f_0-\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2 f_0}{\partial v^2} = \chi_0 J_0^\text{B}\left[f_0,f\right], \end{align} \tag{ 33 }$$

where the Boltzmann–Lorentz operator $J_0^\text{B}[f_0,f]$ is

$$\begin{align} J_0^\text{B}\left[f_0,f\right] & = \overline{\sigma}^{d-1}\int \textrm{d}\mathbf{v}_2\int \text{d}\widehat{\boldsymbol{\sigma}}\Theta \left(\widehat{\boldsymbol{\sigma}} \cdot\mathbf{g}_{12}\right) \left(\widehat{\boldsymbol{\sigma}}\cdot\mathbf{g}_{12}\right)\notag\\ & \quad\times \left[\alpha_0^{-2}f_0\left(\mathbf{v}_1^{^{\prime\prime}},t\right)f\left(\mathbf{v}_2^{^{\prime\prime}},t\right)-f_0\left(\mathbf{v}_1,t\right)f\left(,\mathbf{v}_2,t\right)\right]. \end{align} \tag{ 34 }$$

The equation for the (steady) partial temperature T₀ can be easily derived from equation (33) as

$$\begin{align} 2\gamma_0 \left(T_\text{b}-T_0\right) = T_0 \zeta_0, \end{align} \tag{ 35 }$$

where

$$\begin{align} \zeta_0 = -\frac{\chi_0}{\textrm{d} n_0 T_0}\int \textrm{d}\mathbf{v}\; m_0 v^2\; J_0^\text{B}\left[f_0,f\right] \end{align} \tag{ 36 }$$

is the partial cooling rate characterizing the rate of energy dissipated by intruder–grain collisions. As in the case of the granular gas, for elastic collisions ($\alpha_0 = \alpha = 1$), ζ₀ = 0, $T_\text{b} = T_0$, and equation (33) has the exact solution

$$\begin{align} f_0\left(\mathbf{v}\right) = n_0 \left(\frac{m_0}{2\pi T_\text{b}}\right)^{d/2} \text{exp} \left(-\frac{m_0 v^2}{2 T_\text{b}}\right). \end{align} \tag{ 37 }$$

As occurs for the granular gas, for inelastic collisions ($\alpha_0 \neq 1$) the solution to equation (33) is not known to date.

A good estimate for the partial temperature T₀ can be obtained by considering the leading Sonine approximation to $f_0(\mathbf{v})$ [36]:

$$\begin{align} f_0\left(\mathbf{v}\right) \to n_0 \pi^{-d/2} \beta^{d/2} v_\text{th}^{-d} e^{-\beta \xi^2}\left\{1+\frac{c_0}{2}\left[\beta^2 \xi^4-\left(d+2\right)\beta \xi^2+\frac{d\left(d+2\right)}{4}\right]\right\}. \end{align} \tag{ 38 }$$

Here,

$$\begin{align} \beta = \frac{m_0 T}{m T_0} \end{align} \tag{ 39 }$$

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

$$\begin{align} c_0 = \frac{1}{d\left(d+2\right)}\frac{m_0^2}{n_0 T_0^2}\int \textrm{d}\mathbf{v}\; v^4 f_0\left(\mathbf{v}\right)-1 \end{align} \tag{ 40 }$$

is the fourth-degree cumulant c₀. The use of the Sonine approximation (38) to f₀ allows us to compute the partial cooling rate ζ₀ by substituting (38) into equation (36) and retaining only linear terms in c and c₀. The expression of the (reduced) cooling rate $\zeta_0^* = \zeta_0/\nu$ can be written as

$$\begin{align} \zeta_0^* = \zeta_{00}+\zeta_{01}c_0+\zeta_{02}c, \end{align} \tag{ 41 }$$

where the explicit forms of ζ₀₀, ζ₀₁, and ζ₀₂ can be found in appendix A.

The cumulant c₀ can be determined by multiplying both sides of the Enskog equation (33) by v⁴ and integrating over v. In dimensionless form, the result is

$$\begin{align} \gamma_0^*\left(1+c_0-\frac{T_\text{b}^*}{T_0^*}\right) = \Sigma_0, \end{align} \tag{ 42 }$$

where $T_0^* = T_0/\mathcal{T}$,

$$\begin{align} \gamma_0^* = \frac{\gamma_0}{\nu} = \frac{\gamma_{0,\text{St}}}{\gamma_{\text{St}}} \frac{R_0}{R}\gamma^*, \end{align} \tag{ 43 }$$

and

$$\begin{align} \Sigma_0 = \frac{\chi_0}{4d\left(d+2\right)}\frac{m_0^2}{n_0 T_0^2 \nu}\int \textrm{d}\mathbf{v}\; v^4 J_0^\text{B}\left[f_0,f\right]. \end{align} \tag{ 44 }$$

Retaining only linear terms in c and c₀, one has the result

$$\begin{align} \Sigma_0 = \Sigma_{00}+\Sigma_{01}c_0+\Sigma_{02}c, \end{align} \tag{ 45 }$$

where the explicit forms of Σ₀₀, Σ₀₁, and Σ₀₂ are provided in appendix A. The expression of c₀ can be easily obtained when one takes into account equation (45) in equation (42). It is given by

$$\begin{align} c_0 = \frac{\gamma_0^*\left(1-\frac{T_\text{b}^*}{T_0^*}\right)-\Sigma_{00}-\Sigma_{02}c}{\Sigma_{01}-\gamma_0^*} .\end{align} \tag{ 46 }$$

Finally, in dimensionless form, equation (35) for $T_0^*$ can be written as

$$\begin{align} 2\gamma_0^* \left(T_\text{b}^*-T_0^*\right) = T_0^* \left(\zeta_{00}+\zeta_{01}c_0+\zeta_{02}c\right). \end{align} \tag{ 47 }$$

Substitution of equations (13) and (46) into equation (47) allows us to determine $T_0^*$ in terms of the parameter space of the system. When intruder and granular gas particles are mechanically equivalent ($m = m_0$, $\sigma = \sigma_0$, and $\alpha = \alpha_0$), then $\gamma^* = \gamma_0^*$, $\zeta^* = \zeta_0^*$, and equation (47) yields $T^* = T_0^*$. This means that energy equipartition applies in the self-diffusion problem. However, in the general case (namely, when collisions are inelastic and the intruder and grains are mechanically different), one has to numerically solve equation (47). As in the free cooling case [43, 50, 51], $T_0^*\neq T^*$ and so there is a breakdown of the energy equipartition, as expected.

The dependence of the temperature ratio $T_0/T$ on the (common) coefficient of restitution $\alpha = \alpha_0$ is plotted in figure 3 for d = 3, $T_\text{b}^* = 1$, and φ = 0.1. Two different mixtures have been considered. For d = 3, a good approximation for χ₀ is [52]

$$\begin{align} \chi_0 = \frac{1}{1-\phi}+3 \frac{\sigma_0}{\sigma+\sigma_0}\frac{\phi}{\left(1-\phi\right)^2}+2 \left(\frac{\sigma_0}{\sigma+\sigma_0}\right)^2 \frac{\phi^2}{\left(1-\phi\right)^3}. \end{align} \tag{ 48 }$$

In addition, in the case of an interstitial fluid with low Reynolds number and moderate densities, computer simulations for polydisperse gas–solid flows [46–48] estimate R₀ as

$$\begin{align} R_0 = 1+\left(R-1\right)\left[a \frac{\sigma_0}{\sigma}+\left(1-a\right)\frac{\sigma_0^2}{\sigma^2}\right], \end{align} \tag{ 49 }$$

where

$$\begin{align} a\left(\phi\right) = 1-2.660\phi+9.096\phi^2-11.338\phi^3. \end{align} \tag{ 50 }$$

Note that for mechanically equivalent particles, one has $R_0 = R$, as should be the case.

Zoom In Zoom Out Reset image size

In agreement with previous results [42] obtained by neglecting c and c₀, figure 3 shows a very tiny impact of the mass and diameter ratios on the temperature ratio $T_0/T$. In fact, this influence is amplified by the scale of the vertical axis. This means that the breakdown of energy equipartititon in granular suspensions is much more modest than in dry granular mixtures [43, 50] where the ratio $T_0/T$ clearly differs from 1 for both disparate mass and diameter ratios and/or strong inelasticity. At a more qualitative level, we see that $T_0\gt T$ ($T_0\lt T$) when the intruder is heavier (lighter) than the particles of the granular gas. This behavior is also present in dry granular mixtures. Again, we find an excellent agreement between theory and simulations. As a complement of figure 3, figure 4 shows c₀ versus α for the same systems as in figure 3. As in the case of c, the magnitude of the cumulant c₀ is very small, showing that the deviation of the homogeneous distribution f₀ from the Maxwell–Boltzmann distribution is imperceptible in granular suspensions. Good agreement between theory and DSMC results is observed, except when $m_0/m = 0.5$ for very small values of α ($\alpha\lesssim 0.3$). However, these discrepancies are in the order of 2$\%$, which is still lower than in the dry case.

Zoom In Zoom Out Reset image size

The objective of this section is to determine the diffusion and mobility transport coefficients of intruders immersed in a granular suspension. As stated before, the diffusion process is induced here by the presence of both a weak concentration gradient $\nabla n_0$ and a weak external field E. The corresponding transport coefficients are obtained by solving the Enskog–Lorentz kinetic equation (24) by means of the Chapman–Enskog method [16]. Since the granular gas is in a homogeneous state, χ₀ is constant in the tracer limit and the Enskog–Lorentz operator adopts the simple form $J_0[f_0,f] = \chi_0 J_0^\text{B}[f_0,f]$.

The intruders may freely exchange momentum and energy in their interaction with the particles of the granular gas; this means that these quantities are not invariants of the Enskog–Lorentz collision operator $J_0[f_0,f]$. Only the number density of intruders n₀ is conserved. Its continuity equation can be easily derived from equation (24) as

$$\begin{align} \frac{\partial n_0}{\partial t} = -\nabla \cdot \mathbf{j}_0, \end{align} \tag{ 51 }$$

where

$$\begin{align} \mathbf{j}_0\left(\mathbf{r}; t\right) = \int \textrm{d}\mathbf{v}\; \mathbf{v}\; f_0\left(\mathbf{r}, \mathbf{v};t\right) \end{align} \tag{ 52 }$$

is the intruder particle flux.

As usual in the Chapman–Enskog method, one assumes the existence of a normal solution where all the space and time dependence of f₀ only occurs through a functional dependence on the hydrodynamic fields. In this problem, the normal solution to f₀ is explicitly generated by expanding this distribution in powers of $\nabla n_0$ and E:

$$\begin{align} f_0 = f_0^{\left(0\right)}+\vartheta f_0^{\left(1\right)}+\vartheta^2 f_0^{\left(2\right)}+\cdots. \end{align} \tag{ 53 }$$

In equation (53), each factor ϑ corresponds to the implicit factors $\nabla n_0$ and E. Here, only terms to first-order in ϑ will be considered. The time derivative $\partial_t$ is also expanded as $\partial_t = \partial_t^{(0)}+\vartheta \partial_t^{(1)}+\cdots$, where

$$\begin{align} \partial_t^{\left(0\right)}n_0 & = 0, \quad \partial_t^{\left(0\right)}T = 2\gamma\left(T_\text{b}-T\right)-\zeta T, \end{align} \tag{ 54 }$$

$$\begin{align} \notag \\[-12pt] \partial_t^{\left(1\right)} n_0 & = -\nabla \cdot \mathbf{j}_0^{\left(0\right)}, \quad \partial_t^{\left(1\right)}T = 0, \end{align} \tag{ 55 }$$

and

$$\begin{align} \mathbf{j}_0^{\left(0\right)} = \int \textrm{d}\mathbf{v}\; \mathbf{v}\; f_0^{\left(0\right)}(\mathbf{v}) .\end{align} \tag{ 56 }$$

As noted in previous works [15, 37, 53], although we are interested in computing the diffusion coefficient under steady-state conditions, the presence of the interstitial fluid introduces the possibility of a local energy unbalance, and hence, the zeroth-order distribution $f_0^{(0)}$ is not in general a stationary distribution. This is because, for arbitrary small deviations from the homogeneous steady state, the energy gained by grains due to collisions with the background fluid cannot be locally compensated by the other cooling terms arising from the viscous friction and the collisional dissipation. Thus, in order to get the diffusion and mobility coefficients in the steady state, one has to first determine the unsteady integral equation obeying both coefficients and then solve it under the steady-state condition (9).

The zeroth-order approximation $f_0^{(0)}$ obeys the kinetic equation

$$\begin{align} \Delta T\frac{\partial f_0^{\left(0\right)}}{\partial T}-\gamma_0\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}f_0^{\left(0\right)}-\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2 f_0^{\left(0\right)}}{\partial v^2} = \chi_0 J_0^\text{B}\left[f_0^{\left(0\right)},f\right], \end{align} \tag{ 57 }$$

where $\Delta\equiv 2\gamma\left(\frac{T_\text{b}}{T}-1\right)-\zeta$. Upon deriving equation (57), we have accounted for the fact that $f_0^{(0)}$ depends on time through its dependence on temperature T. In the steady state ($\Delta = 0$), equation (57) has the same form as equation (33). This means that $f_0^{(0)}$ is the solution of equation (33) but taking into account the local dependence of the density n₀. An approximate form to $f_0^{(0)}$ is given by the Sonine approximation (38). Since f₀ is isotropic in velocity, then $\mathbf{j}_0^{(0)} = \mathbf{0}$ and hence $\partial_t^{(1)} n_0 = 0$.

To first order in ϑ, one achieves the kinetic equation

$$\begin{align} -\gamma_0\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}f_0^{\left(1\right)}-\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2 f_0^{\left(1\right)}}{\partial v^2}-\chi_0 J_0^\text{B}\left[f_0^{\left(1\right)},f\right] = -f_0^{\left(0\right)}\mathbf{v}\cdot \nabla \text{ln} n_0-\frac{\mathbf{E}}{m_0}\cdot \frac{\partial}{\partial \mathbf{v}}f_0^{\left(0\right)}. \end{align} \tag{ 58 }$$

Upon obtaining equation (58), we have considered steady conditions ($\Delta = 0$) and have taken into account that $\nabla f_0^{(0)} = f_0^{(0)}\nabla \text{ln} n_0$. The solution to equation (58) can be written as

$$\begin{align} f_0^{\left(1\right)}\left(\mathbf{v}\right) = \boldsymbol{\mathcal A}\left(\mathbf{v}\right) \cdot \nabla \text{ln} n_0+\boldsymbol{\mathcal B}\left(\mathbf{v}\right)\cdot \mathbf{E}, \end{align} \tag{ 59 }$$

where the coefficients $\boldsymbol{\mathcal A}$ and $\boldsymbol{\mathcal B}$ are functions of the velocity and the hydrodynamic fields. Substitution of equation (59) into equation (58) yields the following set of linear integral equations for the unknowns $\boldsymbol{\mathcal A}$ and $\boldsymbol{\mathcal B}$:

$$\begin{align} -\gamma_0\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}\boldsymbol{\mathcal A}-\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2 \boldsymbol{\mathcal A}}{\partial v^2}-\chi_0 J_0^\text{B}\left[\boldsymbol{\mathcal A},f\right] & = -\mathbf{v} f_0^{\left(0\right)}, \end{align} \tag{ 60 }$$

$$\begin{align} \notag \\[-12pt] -\gamma_0\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}\boldsymbol{\mathcal B}-\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2 \boldsymbol{\mathcal B}}{\partial v^2}-\chi_0 J_0^\text{B}\left[\boldsymbol{\mathcal B},f\right] & = -\frac{1}{m_0}\frac{\partial}{\partial \mathbf{v}}f_0^{\left(0\right)} .\end{align} \tag{ 61 }$$

In the first order of $\nabla n_0$ and E, the intruder particle flux has the form

$$\begin{align} \mathbf{j}_0^{\left(1\right)} = -D\nabla \text{ln} n_0+\lambda \mathbf{E}, \end{align} \tag{ 62 }$$

where D is the diffusion coefficient and λ is the mobility coefficient. Since

$$\begin{align} \mathbf{j}_0^{\left(1\right)} = \int \textrm{d}\mathbf{v}\; \mathbf{v}\; f_0^{\left(1\right)}\left(\mathbf{v}\right), \end{align} \tag{ 63 }$$

then, according to equation (59), D and λ are defined as

$$\begin{align} D = -\frac{1}{d}\int \textrm{d}\mathbf{v}\; \mathbf{v}\cdot \boldsymbol{\mathcal A}\left(\mathbf{v}\right), \quad \lambda = \frac{1}{d}\int \textrm{d}\mathbf{v}\; \mathbf{v}\cdot \boldsymbol{\mathcal B}\left(\mathbf{v}\right) .\end{align} \tag{ 64 }$$

For elastic collisions ($\alpha = \alpha_0 = 1$), $T = T_0 = T_\text{b}$ and $f_0^{(0)}(\mathbf{v})$ is the local equilibrium distribution (37). In this case, $\partial f_0^{(0)}/\partial \mathbf{v} = -(m_0 \mathbf{v}/T_\text{b})f_0^{(0)}$ and the integral equations (60) and (61) lead to the identity $\boldsymbol{\mathcal A} = -T_\text{b} \boldsymbol{\mathcal B}$. As a consequence, the conventional Einstein relation is verified; namely,

$$\begin{align} \epsilon = \frac{D}{T_\text{b} \lambda} = 1 .\end{align} \tag{ 65 }$$

On the other hand, for inelastic collisions, $T\neq T_0\neq T_\text{b}$, and hence the relationship between D and λ is no longer simple. There are in principle three different reasons for which the Einstein relation (65) is not verified for granular suspensions. First, when α < 1, the granular temperature T is different from the bath temperature T_b ($T\lt T_\text{b}$). Second, there is a breakdown of the energy equipartition ($T\neq T_0$) when intruders are mechanically different to the particles of the granular gas. Finally, as a third reason, since $f_0^{(0)}$ is not a Gaussian distribution, then $\partial f_0^{(0)}/\partial \mathbf{v}\neq -(m_0 \mathbf{v}/T_\text{b})f_0^{(0)}$ and hence D is not proportional to λ. The first two reasons of discrepancy can be avoided if one replaces the bath temperature T_b by the intruder particle T₀ in the Einstein relation (65). This change leads to the modified Einstein relation

$$\begin{align} \epsilon_0 = \frac{D}{T_{0} \lambda}. \end{align} \tag{ 66 }$$

The relation (66) was proposed by Barrat et al [3] to extend the Einstein relation (65) to granular gases.

Note that, in particular, if one takes the Maxwellian approximation (37) for $f_0^{(0)}$ with T₀ instead of T_b, then $\partial f_0^{(0)}/\partial \mathbf{v} = -(m_0 \mathbf{v}/T_{0})f_0^{(0)}$ and hence ε₀ = 1. Thus, it seems that the only reason for which $\epsilon_0\neq 1$ is due to the absence of the Gibbs state (non-Gaussian behavior of the distribution $f_0^{(0)}$). Since we have seen that the magnitude of the kurtosis c₀ is in general very small for granular suspensions (see for instance figure 4), one expects that the deviations of ε₀ from 1 can be quite difficult to detect in computer simulation experiments. In fact, molecular dynamics simulations [3] (for a similar sort of thermostat as the one employed in this paper) did not observe any deviation from the modified Einstein relation (ε₀ = 1) for a wide range of values of the coefficients of restitution and parameters of the mixture. Our objective here is to assess the departure of ε₀ from 1 in a granular suspension modeled by a stochastic bath with viscous friction.

4.1. Second Sonine approximation to D and λ

It is quite apparent that the transport coefficients D and λ are given in terms of the solution of the integral equations (60) and (61), respectively. These equations can be approximately solved by using a Sonine polynomial expansion. Here, as mentioned in section 1, we determine D and λ up to the second Sonine approximation. In this case, $\boldsymbol{\mathcal A}(\mathbf{v})$ and $\boldsymbol{\mathcal B}(\mathbf{v})$ are approximated by

$$\begin{align} \boldsymbol{\mathcal A}\left(\mathbf{v}\right)\to -f_{0,\text{M}}\left(\mathbf{v}\right)\left[a_1 \mathbf{v}+a_2 \mathbf{S}_0\left(\mathbf{v}\right)\right], \quad \boldsymbol{\mathcal B}\left(\mathbf{v}\right)\to -f_{0,\text{M}}\left(\mathbf{v}\right)\left[b_1 \mathbf{v}+b_2 \mathbf{S}_0\left(\mathbf{v}\right)\right], \end{align} \tag{ 67 }$$

where

$$\begin{align} f_{0,\text{M}}\left(\mathbf{v}\right) = n_0 \left(\frac{m_0}{2\pi T_0}\right)^{d/2} \text{exp} \left(-\frac{m_0 v^2}{2 T_0}\right), \end{align} \tag{ 68 }$$

and $\mathbf{S}_0(\mathbf{v})$ is the polynomial

$$\begin{align} \mathbf{S}_0\left(\mathbf{v}\right) = \left(\frac{1}{2}m_0 v^2-\frac{d+2}{2}T_0\right)\mathbf{v}. \end{align} \tag{ 69 }$$

The Sonine coefficients a₁, b₁, a₁, and a₂ are defined as

$$\begin{align} \left( \begin{array}{c} a_{1}\\ b_1 \end{array} \right) = -\frac{m_0}{{d}n_0T_0}\int \textrm{d}\textbf{v}\,\textbf{v} \cdot \left( \begin{array}{c} \boldsymbol{\mathcal A}\\ \boldsymbol{\mathcal B} \end{array} \right), \quad \left( \begin{array}{c} a_{2}\\ b_2 \end{array} \right) = -\frac{2}{d\left(d+2\right)}\frac{m_0}{n_0T_0^3}\int \textrm{d}\textbf{v}\,\textbf{S}_0 \cdot \left( \begin{array}{c} \boldsymbol{\mathcal A}\\ \boldsymbol{\mathcal B} \end{array} \right). \end{align} \tag{ 70 }$$

According to equation (64), $a_1 = m_0D/(n_0 T_0)$ and $b_1 = -m_0\lambda/(n_0 T_0)$. The evaluation of the coefficients a₁, b₁, a₁, and a₂ is carried out in appendix B.

Knowledge of the Sonine coefficients allows us to determine the first and second Sonine approximations to the diffusion coefficient D and the mobility coefficient λ. To write these expressions, it is convenient to introduce the dimensionless coefficients

$$\begin{align} D^* = \frac{m_0\nu}{T n_0}D, \quad \lambda^* = \frac{m_0\nu}{n_0}\lambda. \end{align} \tag{ 71 }$$

The second Sonine approximation $D^*[2]$ to $D^*$ can be written as

$$\begin{align} D^*\left[2\right] = \frac{\left(\nu_4^*+3\gamma_0^*-c_0 \nu_2^*\right)\tau_0}{\left(\nu_1^*+\gamma_0^*\right)\left(\nu_4^*+3\gamma_0^*\right)-\nu_2^*\left[\nu_3^*+2\gamma_0^* \left(1-\frac{T_\text{b}^*}{T_0^*}\right)\right]}, \end{align} \tag{ 72 }$$

where $\tau_0 = T_0/T$ is the temperature ratio. The expressions of the (reduced) collision frequencies $\nu_1^*$–$\nu_4^*$ can be found in appendix B. Equation (72) agrees with previous results derived in [42] when one takes $c = c_0 = 0$. The second Sonine approximation $\lambda^*[2]$ to $\lambda^*$ is given by

$$\begin{align} \lambda^*\left[2\right] = \frac{\nu_4^*+3\gamma_0^*}{\left(\nu_1^*+\gamma_0^*\right)\left(\nu_4^*+3\gamma_0^*\right)-\nu_2^*\left[\nu_3^*+2\gamma_0^* \left(1-\frac{T_\text{b}^*}{T_0^*}\right)\right]}.\end{align} \tag{ 73 }$$

4.2. DSMC simulations of D and λ

As in the case of the temperatures and the cumulants, to check the accuracy of the Sonine approximations we have solved the Enskog–Lorentz equation by means of the DSMC method described in section 2.1. The diffusion D and mobility λ transport coefficients have been computed separately.

Firstly, the calculation of the diffusion coefficient proceeds in the absence of any external field acting upon the intruder particles. In this scenario, equation (51) reads

$$\begin{align} \frac{\partial n_0}{\partial t} = -\frac{D}{n_0}\nabla^2 n_0, \end{align} \tag{ 74 }$$

where use has been made of the intruder particle flux equation (62) when $\mathbf{E} = \mathbf{0}$. In particular, the coefficient D can be ascertained by evaluating the mean square displacement of the intruders [6], as derived from the standard diffusion equation (74). Specifically, we have

$$\begin{align} \frac{\partial}{\partial t}\langle |\mathbf{r}\left(t\right)-\mathbf{r}\left(0\right)|^2\rangle = 2 d \frac{D_0}{n_0}. \end{align} \tag{ 75 }$$

In this context, $\langle|\mathbf{r}(t)-\mathbf{r}(0)|\rangle$ is the ensemble-average distance traveled by the intruder up to the time t.

On the other hand, the mobility of a tracer particle can be measured by applying a persistent yet small drag force $\mathbf{E} = E \mathbf{e}_x$ to the intruder particles. Over extended time intervals, the perturbed particles will reach a constant velocity λ, which is directly linked to the average distance traveled by the intruders by [54]

$$\begin{align} \langle\left(\mathbf{r}\left(t\right)-\mathbf{r}\left(0\right)\right)\cdot\mathbf{e}_x\rangle\approx\lambda E t. \end{align} \tag{ 76 }$$

The linearity of equation (76) has been checked in [54] by changing the amplitude of the perturbation E.

Figures 5 and 6 show the dependence of the (reduced) transport coefficients $D(\alpha)/D(1)$ and $\lambda(\alpha)/\lambda(1)$ for a three-dimensional (d = 3) system with $T_\text{b}^*$, φ = 0.1, and two different mixtures. Here, the diffusion D and mobility λ coefficients have been reduced with respect to their elastic limits D(1) and $\lambda(1)$, respectively. Theoretical predictions given by the first and second Sonine approximations are compared with DSMC simulations. Although we observe that the first Sonine approximation compares quite well with simulations, some small differences appear in the case of the diffusion coefficient for small mass ratios. These differences are mitigated by the second Sonine approximation since it yields an excellent agreement with the DSMC results. Moreover, while the (reduced) mobility coefficient always increases with decreasing α, a non-monotonic dependence of the (reduced) diffusion coefficient is present regardless of the mass ratio considered. Figures 5 and 6 also highlight that the effect of the mass ratio on $\lambda(\alpha)/\lambda(1)$ is much more significant than for $D(\alpha)/D(1)$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Einstein relation

Once the transport coefficients are known, the conventional and modified relations can be explicitly obtained in terms of the parameters of the system up to the second Sonine approximation. In the case of the conventional Einstein relation (65), one gets the result

$$\begin{align} \epsilon\left[2\right] = \frac{D\left[2\right]}{T_\text{b}\lambda\left[2\right]} = T^*\frac{D^*\left[2\right]}{\lambda^*\left[2\right]} = T_0^*\left(1-\frac{\nu_2^*}{\nu_4^*+3\gamma_0^*}c_0\right), \end{align} \tag{ 77 }$$

while in the case of the modified Einstein relation (66), one achieves the expression

$$\begin{align} \epsilon_0\left[2\right] = \frac{D\left[2\right]}{T_0\lambda\left[2\right]} = 1-\frac{\nu_2^*}{\nu_4^*+3\gamma_0^*}c_0. \end{align} \tag{ 78 }$$

It is quite apparent that while the conventional Einstein relation (77) fails due to both energy non-equipartition and non-Gaussian corrections to the distribution $f_0^{(0)}$, the departure of the modified Einstein relation (78) is only due to the latter feature ($c_0\neq 0$).

To illustrate the dependence of both Einstein relations on the (common) coefficient of restitution $\alpha = \alpha_0$, figure 7 shows ε and ε₀ for several mixtures. While $\epsilon_0 \simeq 1$ for all the mixtures (in fact, the three curves practically collapse in a common curve on the scale of the vertical axis of figure 7), there are significant deviations from 1 in the conventional Einstein relation. In fact, the deviations of ε₀ from 1 are much smaller than 1$\%$. This result contrasts with the ones previously derived for freely cooling [7] and driven (with a Gaussian thermostat) [17] granular gases.

Zoom In Zoom Out Reset image size

In summary, the results derived here show no new surprises relative to the earlier work for dry granular gases [17, 18]: the origin of the breakdown of the modified Einstein relation is only due to the departure of the reference state from the Maxwell–Boltzmann distribution. However, in contrast to the previous works [17, 18], the deviation of ε₀ from 1 in granular suspensions is much smaller than the one found in driven granular gases, even for moderate densities and/or strong inelasticity.

The main objective of this paper was to analyze the validity of the conventional $\epsilon = D/T \lambda = 1$ and modified $\epsilon_0 = D/T_0 \lambda = 1$ Einstein relations in a moderately dense granular suspension. The results were derived in the framework of the (inelastic) Enskog kinetic equation, which applies to moderate densities. As usual in granular suspensions and due to the difficulties embodied in the description of systems constituted by two or more phases, a coarse-grained approach was adopted. In this approach, the influence of the interstitial fluid on grains and intruders was modeled through two different forces. Each of the forces is composed of two terms: (i) a viscous drag term plus (ii) a stochastic Langevin-like term defined in terms of the background temperature T_b. Two different friction coefficients were introduced in the model; each of them accounts for the interaction between the grains and intruders with the external bath. Thus, the starting kinetic equations for grains and intruders are the Enskog and the Enskog–Lorentz equations, respectively, with the addition of Fokker–Planck terms to each one of the above master equations. The present work extends previous studies performed by one of the authors of this paper in the case of driven dry granular gases [17, 18].

To determine the explicit dependence of ε and ε₀ on the parameter space of the system, the corresponding Enskog–Lorentz kinetic equation for intruders was solved by means of the Chapman–Enskog method [16] up to the first order in both the density gradient and the external field. As for molecular mixtures, the diffusion D and mobility λ transport coefficients are given in terms of a set of coupled linear integral equations, which are approximately solved by expanding the unknowns in a series of Sonine polynomials. Here, the series was truncated by considering the two first relevant Sonine polynomials; this leads to the so-called first and second Sonine approximations to the coefficients D and λ. The reliability of these theoretical results was assessed via a comparison with computer simulations obtained by numerically solving the Enskog–Lorentz equation by means of the DSMC method [40]. The comparison shows in general a good agreement between theory and simulations, especially in the case of the second Sonine solution. This agreement shows again the reliability of kinetic theory to obtain the transport properties of granular gases.

As expected from previous works [17, 18], our results show that while the conventional Einstein relation is clearly violated, the deviations of the modified Einstein relation ε₀ from 1 are very tiny. In particular, the deviations of ε₀ from 1 are in general smaller than 1% in the range of inelasticities and densities studied. This means that these deviations are much smaller than the ones reported in [18] for moderate densities, especially when the gas is driven by the Gaussian thermostat (see figure 6 of [18]). The fact that $\epsilon_0\simeq 1$ for granular suspensions is essentially due to the small magnitude of the cumulant c₀, which is much smaller than in the driven case (compare, for instance, figure 2 of [18] with figure 4 of the present work).

On the other hand, the above conclusion disagrees with the computer simulation results obtained years ago by Puglisi et al [19], which were subsequently confirmed in an experiment [20]. In this experiment, the granular gas is driven by a shaker, and the tracer is a rotating wheel immersed in the gas. The main claim in both papers is that the violation of the modified Einstein relation mainly originates from the presence of spatial and velocity correlations, which are relevant as the density increases. Given that the Enskog equation takes into account the spatial correlations (through the pair correlation functions) but neglects the velocity correlation between the velocities of the particles that are about to collide (molecular chaos hypothesis), one could argue that the deviation of ε₀ from 1 in the Enskog theory could be more important as both the density and inelasticity increase. However, our results indicate that the violation of the modified Einstein relation is still very small (and, hence, undetectable in computer simulations) even when one considers high densities and/or strong inelasticity. In this context, and based on the Enskog results for granular suspensions at moderate densities, one can conclude that the origin of the deviation of ε₀ from 1 is mainly due to velocity correlations, which are absent in the Enskog theory. These velocity correlations are expected to have a significant impact on ε₀ for relatively high densities and/or high inelasticities.

In connection with the above point, one could include this sort of velocity correlation in the collision operator [55]. However, as mentioned in [18], the inclusion of this new ingredient in the Enskog collision operator makes analytical calculations intractable since higher-order correlations should be accounted for in the evaluation of the collision integrals. This type of calculation contrasts with the ones offered in this paper, where the diffusion and mobility transport coefficients have been explicitly determined in terms of masses, diameters, coefficients of restitution, density, and background temperature.

Although some simulation computer works [56–58] have clearly shown the failure of the molecular chaos assumption for inelastic collisions as the density increases, there is also some evidence in the granular literature on the usefulness of the Enskog theory for densities outside the dilute limit and inelasticities beyond the quasielastic limit. This evidence is supported by the agreement found at the level of the macroscopic properties between the Enskog results [59–62] and those obtained from computer simulations [60, 63–66] and real experiments [67–69].

One of the limitations of the theoretical results presented in this paper is that they are approximated since they have been obtained by considering the second Sonine approximation in the Chapman–Enskog solution. Exact results can be derived if one determines the coefficients D and λ by starting from the inelastic Maxwell model (IMM) for a dilute gas. As for molecular Maxwell gases [16], the collision rate of colliding particles in the IMM is assumed to be independent of the relative velocity. This simplification allows any moment of degree k of the Boltzmann collisional operator to be expressed in terms of velocity moments of degree k or less than k [36]. This feature of the Boltzmann collision operator of the IMM opens up the possibility of exactly determining the coefficients D and λ. These results are presented in the appendix C. According to these results, one concludes that the modified Einstein relation applies for IMM in any number of dimensions.

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.

In this appendix we display the explicit expressions of the (reduced) partial cooling rate $\zeta_0^*$ and the fourth-degree collisional moment Σ₀. Their forms are provided by equations (41) and (44), respectively, when nonlinear terms in c₀ and c are neglected. The expressions of ζ₀₀, ζ₀₁, and ζ₀₂ are given by [70]

$$\begin{align} \zeta_{00} & = \frac{2\sqrt{2}}{d}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_{0}}{\chi}\mu \left(\frac{1+\beta}{\beta}\right)^{1/2}\left(1+\alpha_{0}\right)\left[1-\frac{1}{2}\mu \left(1+\alpha_{0}\right)\left(1+\beta\right) \right], \end{align} \tag{ A1 }$$

$$\begin{align} \notag \\[-12pt] \zeta_{01} & = \frac{1}{2\sqrt{2}d}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_{0}}{\chi}\mu \frac{\left(1+\beta\right)^{-3/2}}{\beta^{1/2}}\left(1+\alpha_{0}\right) \left[3+4\beta-\frac{3}{2}\mu\left(1+\alpha_{0}\right)\left(1+\beta\right) \right], \end{align} \tag{ A2 }$$

$$\begin{align} \notag \\[-12pt] \zeta_{02} & = -\frac{1}{2\sqrt{2}d}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_{0}}{\chi}\mu \left(\frac{1+\beta}{\beta}\right)^{-3/2}\left(1+\alpha_{0}\right) \left[1+\frac{3}{2}\mu\left(1+\alpha_{0}\right)\left(1+\beta\right) \right]. \end{align} \tag{ A3 }$$

In the case of the fourth-degree collisional moment Σ₀, the expressions of Σ₀₀, Σ₀₁, and Σ₀₂ are [70]

$$\begin{align} \Sigma_{00} & = \frac{1}{\sqrt{2}d\left(d+2\right)} \left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_{0}}{\chi}\mu \left[ \beta\left(1+\beta\right)\right]^{-1/2}\left(1+\alpha_0\right)\Bigg\{-2\left[d+3+\left(d+2\right)\beta\right] \notag\\ & \quad +\mu\left( 1+\alpha_{0}\right) \left( 1+\beta \right)\left( 11+ d+\frac{d^2+5d+6}{d+3} \beta \right)-8\mu^{2}\left(1+\alpha_{0}\right)^{2}\left(1+\beta \right)^{2} \notag\\ & \quad +2\mu^{3}\left( 1+\alpha_{0}\right)^{3}\left( 1+\beta \right)^{3}\Bigg\}, \end{align} \tag{ A4 }$$

$$\begin{align} \Sigma_{01} & = \frac{1}{8\sqrt{2}d\left(d+2\right)} \left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_{0}}{\chi}\mu \beta^{-1/2}\left(1+\beta\right)^{-5/2} \left(1+\alpha_{0}\right) \notag\\ & \quad \times \Bigg\{-2\Big[ 45+15d+\left(114+39d\right)\beta+\left(88+32d\right)\beta^{2} +\left(16+8d\right)\beta^{3}\Big]\nonumber\\ & \quad +3\mu\left(1+\alpha_{0}\right) \left(1+\beta \right) \left[ 55+5d+9\left(10+d\right)\beta+4\left(8+d\right)\beta ^{2}\right]\nonumber\\ & \quad -24\mu^{2}\left(1+\alpha_{0}\right)^{2}\left(1+\beta \right)^{2}\left(5+4\beta \right)+30\mu^{3}\left(1+\alpha_{0}\right)^{3}\left(1+\beta \right)^{3}\Bigg\}, \end{align} \tag{ A5 }$$

$$\begin{align} \notag \\[-12pt] \Sigma_{02} & = \frac{1}{8\sqrt{2}d\left(d+2\right)} \left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_{0}}{\chi}\mu \beta^{3/2}\left(1+\beta\right)^{-5/2} \left( 1+\alpha_{0}\right) \Bigg\{2\left[d-1+\left(d+2\right)\beta\right]\nonumber\\ & \quad +3\mu\left(1+\alpha_{0}\right) \left(1+\beta\right) \left[d-1+\left(d+2\right)\beta\right] -24\mu^{2}\left( 1+\alpha_{0}\right)^{2}\left(1+\beta \right)^{2} \nonumber\\ & \quad +30\mu^{3}\left(1+\alpha_{0}\right)^{3}\left(1+\beta\right)^{3}\Bigg\} . \end{align} \tag{ A6 }$$

Some details are provided in this appendix in the calculation of the diffusion D and mobility λ coefficients up to the second Sonine approximation. The substitution of equation (67) into the integral equations (60) and (61), respectively, gives

$$\begin{align} \gamma_0\frac{\partial}{\partial\mathbf{v}}& \cdot\mathbf{v} (a_1 f_{0\text{M}}\mathbf{v} +a_2 f_{0\text{M}}\mathbf{S}_0) +\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2}{\partial v^2}\left(a_1 f_{0\text{M}}\mathbf{v}+a_2 f_{0\text{M}}\mathbf{S}_0\right)+a_1\chi_0 J_0^\text{B}\left[f_{0\text{M}}\mathbf{v},f\right] \notag\\ & \quad +a_2\chi_0 J_0^\text{B}\left[f_{0\text{M}}\mathbf{S}_0,f\right] = -\mathbf{v}f_0^{\left(0\right)}, \end{align} \tag{ B1 }$$

$$\begin{align} \notag \\[-12pt] \gamma_0\frac{\partial}{\partial\mathbf{v}}& \cdot\mathbf{v} (b_1 f_{0\text{M}}\mathbf{v} + b_2 f_{0\text{M}}\mathbf{S}_0 ) +\frac{\gamma_0 T_{\text{b}}}{m_0}\frac{\partial^2}{\partial v^2}\left(b_1 f_{0\text{M}}\mathbf{v}+b_2 f_{0\text{M}}\mathbf{S}_0\right)+b_1\chi_0 J_0^\text{B}\left[f_{0\text{M}}\mathbf{v},f\right]\notag\\ & \quad + b_2\chi_0 J_0^\text{B}\left[f_{0\text{M}}\mathbf{S}_0,f\right] = -\frac{1}{m_0} \frac{\partial}{\partial\mathbf{v}}f_0^{\left(0\right)}. \end{align} \tag{ B2 }$$

Next, equations (B1) and (B2) are multiplied by v and integrated over the velocity. The result is

$$\begin{align} \left(\gamma_0+\nu_1\right)D+\frac{n_0 T_0^2}{m_0} \nu_2 a_2 = \frac{n_0T_0}{m_0}, \end{align} \tag{ B3 }$$

$$\begin{align} \notag \\[-12pt] \left(\gamma_0+\nu_1\right)\lambda-\frac{n_0 T_0^2}{m_0} \nu_2 b_2 = \frac{n_0}{m_0}, \end{align} \tag{ B4 }$$

where use has been made of the identities $a_1 = (m_0 D/n_0 T_0)$, $b_1 = -(m_0 \lambda/n_0T_0)$, and we have introduced the quantities

$$\begin{align} \nu_1 = -\frac{m_0 \chi_0}{d n_0 T_0}\int \textrm{d}\mathbf{v}\; \mathbf{v}\cdot J_0^\text{B}\left[f_{0\text{M}}\mathbf{v},f\right], \quad \nu_2 = -\frac{m_0 \chi_0}{\textrm{d} n_0 T_0^2}\int \textrm{d}\mathbf{v}\; \mathbf{v}\cdot J_0^\text{B}\left[f_{0\text{M}}\mathbf{S}_0,f\right]. \end{align} \tag{ B5 }$$

If only the first Sonine correction is retained (i.e. $a_2 = b_2 = 0$), then $D[1] = T_0 \lambda[1]$ and the modified Einstein relation (66) is verified.

To get the second Sonine coefficients a₂ and b₂, one multiplies equations (B1) and (B2) by $\mathbf{S}_0(\mathbf{v})$ and integrates over the velocity. The result is

$$\begin{align} \frac{m_0}{n_0T_0^2}\left[2\gamma_0\left(1-\frac{T_\text{b}}{T_0}\right)+\nu_3\right]D+\left(3\gamma_0+\nu_4\right)a_2 = \frac{c_0}{T_0}, \end{align} \tag{ B6 }$$

$$\begin{align} \notag \\[-12pt] -\frac{m_0}{n_0 T_0^2}\left[2\gamma_0\left(1-\frac{T_\text{b}}{T_0}\right)+\nu_3\right]\lambda+\left(3\gamma_0+\nu_4\right)b_2 = 0, \end{align} \tag{ B7 }$$

where

$$\begin{align} \nu_3 & = -\frac{2}{d\left(d+2\right)}\frac{m_0 \chi_0}{n_0 T_0^2}\int \textrm{d}\mathbf{v}\; \mathbf{S}_0\cdot J_0^\text{B}\left[f_{0\text{M}}\mathbf{v},f\right], \notag\\ \nu_4 & = -\frac{2}{d\left(d+2\right)}\frac{m_0 \chi_0}{n_0 T_0^3}\int \textrm{d}\mathbf{v}\; \mathbf{S}_0\cdot J_0^\text{B}\left[f_{0\text{M}}\mathbf{S}_0,f\right] .\end{align} \tag{ B8 }$$

In reduced units and by using matrix notation, equations (B3) and (B6) along with equations (B4) and (B7) can be rewritten as

$$\begin{align} \left( \begin{array}{cc} \gamma_0^*+\nu_1^*&\tau_0^2 \nu_2^*\\ \frac{\nu_3^*+2\gamma_0^*\left(1-\frac{T_\text{b}^*}{T_0^*}\right)}{\tau_0^2}&3\gamma_0^*+\nu_4^* \end{array} \right) \left( \begin{array}{c} D^*\\ a_2^* \end{array} \right) & = \left( \begin{array}{c} \tau_0\\ \frac{c_0}{\tau_0} \end{array} \right), \end{align} \tag{ B9 }$$

$$\begin{align} \notag \\[-12pt] \left( \begin{array}{cc} \gamma_0^*+\nu_1^*&-\tau_0^2 \nu_2^*\\ \nu_3^*+2\gamma_0^*\left(1-\frac{T_\text{b}^*}{T_0^*}\right)&-\tau_0^2 \left(3\gamma_0^*+\nu_4^*\right) \end{array} \right) \left( \begin{array}{c} \lambda^*\\ b_2^* \end{array} \right) & = \left( \begin{array}{c} 1\\ 0 \end{array} \right). \end{align} \tag{ B10 }$$

Here, $\nu_i^* = \nu_i/\nu$ ($i = 1,\ldots,4$), $a_2^* = T\nu a_2$, and $b_2^* = T^2 \nu b_2$. From equations (B9) and (B10) one obtains the expressions (72) and (73) for the second Sonine approximations to $D^*$ and $\lambda^*$, respectively.

The integrals involving the (reduced) collision frequencies $\nu_i^*$ have been computed in previous works [62, 71, 72] for a d-dimensional system when f is replaced by the Maxwellian distribution

$$\begin{align} f_\text{M}\left(\mathbf{v}\right) & = n \left(\frac{m}{2\pi T}\right)^{d/2} \text{exp} \left(-\frac{m v^2}{2 T}\right). \end{align} \tag{ B11 }$$

In this case, the collision frequencies are given by

$$\begin{align} \nu_{1}^* & = \frac{\sqrt{2}}{d}\left(\frac{\overline{\sigma}}{\sigma} \right)^{d-1}\frac{\chi_0}{\chi}\mu \left(1+\alpha_0\right) \left(\frac{1+\beta}{\beta}\right)^{1/2}, \notag\\ \nu_{2}^* & = \frac{1} {\sqrt{2}d}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\frac{\chi_0}{\chi}\mu\left(1+\alpha_0\right)\left[\beta\left(1+\beta\right)\right]^{-1/2}, \end{align} \tag{ B12 }$$

$$\begin{align} \nu_{3}^* = \frac{\sqrt{2}} {d\left(d+2\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1} \frac{\chi_0}{\chi}\mu\left(1+\alpha_0\right)\left(\frac{\beta}{1+\beta}\right)^{1/2}A_c, \end{align} \tag{ B13 }$$

$$\begin{align} \nu_{4}^* = \frac{1} {\sqrt{2}d\left(d+2\right)}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1} \frac{\chi_0}{\chi}\mu\left(1+\alpha_0\right)\left(\frac{\beta}{1+\beta}\right)^{3/2} \left[A_d-\left(d+2\right)\frac{1+\beta}{\beta} A_c\right], \end{align} \tag{ B14 }$$

where

$$\begin{align} A_c & = \left(d+2\right)\left(1+2\varpi\right)+\mu\left(1+\beta\right)\left\{\left(d+2\right)\left(1-\alpha_0\right) -\left[\left(11+d\right)\alpha_0-5d-7\right]\varpi\beta^{-1}\right\} \notag\\ & \quad +3\left(d+3\right)\varpi^2\beta^{-1} +2\mu^2\left(2\alpha_0^{2}-\frac{d+3}{2}\alpha _{0}+d+1\right)\beta^{-1}\left(1+\beta\right)^2 - \left(d+2\right)\beta^{-1}\left(1+\beta\right), \end{align} \tag{ B15 }$$

$$\begin{align} A_d & = 2\mu^2\left(\frac{1+\beta}{\beta}\right)^{2} \left(2\alpha_0^{2}-\frac{d+3}{2}\alpha_0+d+1\right) \left[d+5+\left(d+2\right)\beta\right]-\mu\left(1+\beta\right) \nonumber\\ & \quad \times \Bigg\{\varpi\beta^{-2}\left[\left(d+5\right)+\left(d+2\right)\beta\right] \left[\left(11+d\right)\alpha_0 -5d-7\right]-\beta^{-1}[20+d\left(15-7\alpha_0\right)\nonumber\\ & \quad +d^2\left(1-\alpha_0\right)-28\alpha_0] -\left(d+2\right)^2\left(1-\alpha_0\right)\Bigg\} +3\left(d+3\right)\varpi^2\beta^{-2}\left[d+5+\left(d+2\right)\beta\right] \notag\\ & \quad + 2\varpi\beta^{-1}\left[24+11d+d^2+\left(d+2\right)^2\beta\right] +\left(d+2\right)\beta^{-1} \left[d+3+\left(d+8\right)\beta\right]\notag\\ & \quad -\left(d+2\right)\left(1+\beta\right)\beta^{-2} \left[d+3+\left(d+2\right)\beta\right].\nonumber \end{align} \tag{ B16 }$$

Here, $\varpi = (\mu_0/T_0)\left(T_0-T\right)$.

In this appendix we provide the exact results derived by considering the IMM for a dilute granular gas. The IMM is a further simplification of the inelastic hard sphere (IHS) model since it assumes that the collision rate of the colliding particles is independent of their relative velocity. In this model, the Boltzmann collision operator $J^{\text{IMM}}[f,f]$ of the granular gas reads [36]

$$\begin{align} J\left[\mathbf{v}_1|f,f\right] = \frac{\nu_\text{M}}{n S_d} \int \textrm{d}\mathbf{v}_2\int \textrm{d}\widehat{\boldsymbol{\sigma}} \left[\alpha^{-1}f\left(\mathbf{v}_1^{^{\prime\prime}}\right)f\left(\mathbf{v}_2^{^{\prime\prime}}\right)-f\left(\mathbf{v}_1\right)f\left(\mathbf{v}_2\right)\right], \end{align} \tag{ C1 }$$

where $S_d = 2\pi^{d/2}/\Gamma(d/2)$ is the total solid angle in d dimensions and the velocities $\mathbf{v}_{1,2}^{^{\prime\prime}}$ are related to $\mathbf{v}_{1,2}$ by equation (3). Moreover, $\nu_\text{M}$ is an effective collision frequency that is independent of velocity. This quantity can be seen as a free parameter of the model to be chosen to optimize the agreement with some proper quantity of interest obtained from the Boltzmann equation for IHS. In particular, if we chose $\nu_\text{M}$ to get the same expression of the cooling rate ζ as the one obtained in the Maxwellian approximation in the IHS model of diameter σ (equation (16) with c = 0), then one obtains the simple relationship $\nu_\text{M} = 2\nu$, where ν is defined by equation (15) with χ = 1.

In the context of IMM, the Boltzmann–Lorentz collision operator $J_0^{\text{IMM}}[f_0,f]$ reads [36]

$$\begin{align} J_0^{\text{IMM}}\left[f_0,f\right] = \frac{\nu_{\text{M},0}}{n S_d}\int \textrm{d}\mathbf{v}_2\int \textrm{d}\widehat{\boldsymbol{\sigma}} \left[\alpha_{0}^{-1}f_0\left(\mathbf{v}_1^{^{\prime\prime}}\right)f\left(\mathbf{v}_2^{^{\prime\prime}}\right)f_0\left(\mathbf{v}_1\right)f\left(\mathbf{v}_2\right)\right], \end{align} \tag{ C2 }$$

where $\nu_{\text{M},0}$ is an effective collision frequency for intruder–gas collisions and the relation between the velocities $\mathbf{v}_{1,2}^{^{\prime\prime}}$ and $\mathbf{v}_{1,2}$ is given by equation (26). When the form (C2) of the operator $J_0^{\text{IMM}}[f_0,f]$ is substituted into the definition (36), the partial cooling rate ζ₀ can be exactly determined for IMM. The result is [36]

$$\begin{align} \zeta_{0} = \frac{2\nu_{\text{M},0}}{d}\mu \left(1+\alpha_{0}\right)\left[1-\frac{1}{2}\mu \left(1+\alpha_{0}\right)\left(1+\beta\right) \right]. \end{align} \tag{ C3 }$$

Comparison of equation (C3) with equation (A1) (it gives ζ₀ for IHS in the Maxwellian approximation, i.e. when $c = c_0 = 0$) yields the relation

$$\begin{align} \nu_{\text{M},0} = \sqrt{2}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1}\left(\frac{1+\beta}{\beta}\right)^{1/2}\nu .\end{align} \tag{ C4 }$$

The determination of the transport coefficients D and λ follows similar mathematical steps to those made in the case of IHS, except that the corresponding collision integrals appearing in the evaluation of these coefficients can be exactly computed. They are given by

$$\begin{align} \int \textrm{d}\mathbf{v}\; \mathbf{v}\cdot J_0^{\text{IMM}}\left[\boldsymbol{\mathcal{A}},f\right] = \nu_{\text{M},0} \mu \left(1+\alpha_0\right)D, \quad \int \textrm{d}\mathbf{v}\; \mathbf{v}\cdot J_0^{\text{IMM}}\left[\boldsymbol{\mathcal{B}},f\right] = \nu_{\text{M},0} \mu \left(1+\alpha_0\right)\lambda. \end{align} \tag{ C5 }$$

The final expressions of D and λ can be easily derived when one takes into account equation (C5). The results are

$$\begin{align} D = \frac{n_0 T_0}{m_0}\left(\gamma_0+\nu_1\right)^{-1}, \quad \lambda = \frac{n_0}{m_0}\left(\gamma_0+\nu_1\right)^{-1}, \end{align} \tag{ C6 }$$

where

$$\begin{align} \nu_1 = \mu \left(1+\alpha_0\right)\frac{\nu_{\text{M},0}}{d} = \frac{\sqrt{2}}{d}\left(\frac{\overline{\sigma}}{\sigma}\right)^{d-1} \left(\frac{1+\beta}{\beta}\right)^{1/2}\mu \left(1+\alpha_0\right). \end{align} \tag{ C7 }$$

Equation (C6) shows that the expressions of D and λ derived for IMM coincide with the ones obtained from the Boltzmann equation for IHS in the first Sonine approximation when one neglects non-Gaussian corrections to the zeroth-order distributions ($c = c_0 = 0$). Thus, according to equation (C6), $\epsilon_0 = D/(T_0 \lambda) = 1$ for IMM.