1 Introduction↩︎

When granular media are externally excited they behave like a fluid. In these conditions (rapid flow conditions), the tools of the classical kinetic theory of gases (conveniently adapted to dissipative dynamics) can be used to derive the corresponding hydrodynamic equations with explicit expressions for the transport coefficients [1]–[5]. A granular gas is typically modeled as a gas of hard spheres with inelastic collisions. In the simplest version, the spheres are assumed to be completely smooth, and inelasticity in collisions is accounted for only through a constant positive coefficient of normal restitution. In the dilute regime, the inelastic version of the Boltzmann equation [3], [5] has been solved by means of the Chapman–Enskog method [6] and the set of coupled linear integral equation verifying the Navier–Stokes transport coefficients has been also derived. As with elastic collisions [6], [7], for hard spheres usually these integral equations are approximately solved by considering the first few terms in a Sonine polynomial expansion of the distribution function. This procedure can be extended to the more realistic case of granular mixtures, namely when grains have different masses and/or diameters. However, the problem is more complex than that of a monocomponent gas since the number of transport coefficients is greater and they depend on more parameters (see, for example, Refs.[8]–[16]). Additionally, determining the transport properties from the Boltzmann equation for both elastic and/or inelastic hard spheres is a very difficult task for far from equilibrium states (i.e., beyond the Navier–Stokes domain).

Due to the aforementioned difficulties, alternative approaches are commonly used in kinetic theory to obtain exact results. One possibility is to maintain the intricate mathematical structure of the Boltzmann collision operator while assuming a different interaction model: the so-called inelastic Maxwell model (IMM). To the best of our knowledge, this interaction model was introduced independently in Refs.[17] and [18] at the beginning of the 21st century. The main reason for introducing IMM was to analyze in a clean way the overpopulation associated with the high energy tails of the distribution function in the homogeneous cooling state (HCS) [18]–[30].

As with the conventional Maxwell molecules [31] (which are defined by a repulsive potential inversely proportional to the fourth power of the distance for a three-dimensional gas), the collision rate for IMM is independent of the relative velocity of the colliding spheres. This contrasts with the inelastic hard sphere (IHS) model, in which the collision rate is proportional to the relative velocity. The main advantage of using IMM instead of IHS is that a collisional moment of degree \(k\) can be expressed in terms of velocity moments of degree \(k\) or smaller than \(k\), without knowing the velocity distribution function. This property allows one to obtain the exact forms of the Navier-Stokes transport coefficients, as well as the rheological properties of sheared granular gases [5]. However, since IMM’s scattering rules are the same as IHS’, these models do not describe real particles because they do not interact according to a given potential law. In any case, many researchers working in the field recognize that the cost of sacrificing physical realism can be offset by the number of exact analytical results obtained using IMM. Apart from this aspect, it should be noted that some experiments involving magnetic grains with dipolar interactions have been well described by IMM [32].

Another possible approach to achieving exact results is considering simpler kinetic models of the Boltzmann equation for IHS. These models are mathematically more tractable than the Boltzmann equation and typically maintain the primary physical properties of the true kinetic equation. In particular, for elastic collisions, the well-known Bhatnagar–Gross–Krook (BGK) kinetic model [33], [34] has proven to be an accurate tool to determine nonlinear transport properties, especially in shearing nonequilibrium states [34], [35]. Several models have been proposed for inelastic collisions in single gases [36]–[38]. However, the number of kinetic models proposed in the literature for multicomponent granular gases is much smaller. We are aware of only one kinetic model reported in the granular literature: the model proposed by Vega Reyes et al. (VGS model) [39]. This contrasts with the large number of kinetic models proposed in the case of molecular mixtures (see, for instance, Refs.[40]–[48]). The VGS model is based on the equivalence between a gas of elastic hard spheres subjected to a drag force proportional to the particle velocity and a gas of IHS [49]. The kinetic model is defined in terms of a relaxation term that can be selected from among the various models proposed for molecular mixtures. Recently [50], the VGS model has been solved in three different nonequilibrium problems. A comparison of the results derived in Ref.[50] with those obtained from the Boltzmann equation for IHS shows a general agreement, especially when dissipation is not strong.

The objective of this paper is to analyze the diffusion of tracer particles in sheared granular gases by starting from the Boltzmann equation for IMM and the VGS kinetic model. The study of tracer diffusion in granular shear flows has attracted the attention of engineers and physicists for years. Additionally, this is a problem of practical interest because granular materials must be mixed before processing can begin. Due to the complexity of the general problem, the limiting situation where the tracer particles are mechanically equivalent to the particles of the granular gas (self-diffusion problem) was widely studied in earlier works. Thus, experimental studies include both systems with macroscopic flows [51], [52] and vertical vibrated systems [53]. As a complement, computer simulation works [54]–[56] on dense systems have mainly analyzed the influence of the solid volume fraction on the elements of the self-diffusion tensor. Analytical studies on this problem are scarce. In the context of kinetic theory, two studies [57], [58] have evaluated the tracer diffusion coefficients in a granular gas under simple or uniform shear flow (USF) using the Boltzmann equation of IHS. Unlike previous studies, the analysis performed in Refs.[57], [58] considers tracer and granular gas particles as mechanically distinct, resulting in energy nonequipartition as the coefficients of restitution decrease. It is important to note that analyzing mass transport in a strongly shearing granular gas is difficult due to the anisotropy in velocity space induced by the shear flow. This gives rise to tensorial quantities (\(D_{k\ell}\), \(D_{p,k\ell}\), and \(D_{T,k\ell}\)) instead of the conventional scalar coefficients (the diffusion coefficient \(D\), the pressure diffusion coefficient \(D_p\), and the thermal diffusion coefficient \(D_T\)) for characterizing the mass transport.

Since the coefficients \(D_{k\ell}\), \(D_{p,k\ell}\), and \(D_{T,k\ell}\) were determined in Refs.[57] and [58] using approximate methods (the leading terms in a Sonine polynomial expansion), it is useful to revisit the study by considering both the Boltzmann equation for IMM and the VGS model, in which the exact forms of the diffusion coefficients can be obtained. Searching for exact solutions in kinetic theory is interesting from both a formal point of view and as a means of gauging the reliability of these types of solutions. Here, we compare these solutions with the (approximate) analytical results obtained form the Boltzmann equation for IHS and in some cases with computer simulations available in the granular literature. This type of comparison measures the degree to which both IMM and kinetic models accurately capture the influence of dissipation in realistic granular flows.

The plan of the paper is as follows. Section 2 introduces the Boltzmann equation and its balance hydrodynamic equations, and presents the form of the Boltzmann collision operators for IMM and the VGS kinetic model. Section 3 analyzes the rheology of a granular mixture under USF in the tracer limit within the context of IMM. Once the rheological properties are determined, the elements of the diffusion tensors are explicitly obtained in section 4 for IMM by solving the Boltzmann equation by means of a generalized Chapman–Enskog expansion [6] around the shear flow distribution. Section 5 briefly shows how to evaluate the tracer diffusion coefficients using the VGS model, while the comparison between the results obtained for IHS, IMM and the VGS model is carried out in section 6. The comparison generally shows reasonable agreement, especially for IMM. As an application, section 7 addresses the problem of segregation by thermal diffusion. For the sake of simplicity, we consider a situation in which the thermal gradient is perpendicular to the shear flow plane (\(xy\)-plane), so only segregation parallel to the thermal gradient occurs in the system. In this situation, the sign of the thermal diffusion factor \(\Lambda_z\) characterizes the tendency of the tracer particles to move towards the hot or cold plate. The paper ends with some concluding remarks in section 8.

2 Boltzmann kinetic equation for granular mixtures↩︎

We consider a granular binary mixture constituted by smooth hard disks (\(d=2\)) or spheres (\(d=3\)) of masses \(m_i\) and diameters \(\sigma_i\) (\(i=1,2\)). The collisions between grains are inelastic and characterized by the (constant) coefficients of normal restitution \(\alpha_{ij}\) (\(0\leq \alpha_{ij}\leq 1\)). The inelasticity in collisions only affects the translational degrees of freedom of grains. At a kinetic theory level, all the relevant information on the state of the mixture is given through the knowledge of the one-particle velocity distribution \(f_i(\mathbf{r}, \mathbf{v}; t)\) of species \(i\). In the low-density regime and in the absence of external forces, the distributions \(f_i\) obey the set of coupled nonlinear (inelastic) Boltzmann equations \[\label{2461} \left(\partial_t+{\boldsymbol{v}}\cdot \nabla \right)f_{i} ({\boldsymbol{r}},{\boldsymbol{v}};t) =\sum_{j}J_{ij}\left[{\boldsymbol{v}}|f_{i}(t),f_{j}(t)\right],\tag{1}\] where \(J_{ij}\left[{\boldsymbol{v}}|f_{i},f_{j}\right]\) are the Boltzmann collision operators [5]. The most relevant hydrodynamic fields in a binary mixture are the number densities \(n_i\), the mean flow velocity \(\mathbf{U}\), and the granular temperature \(T\). In terms of the distributions \(f_i\), those fields are defined as \[\label{2462} n_i=\int \mathrm{d}{\boldsymbol{v}} f_i({\boldsymbol{v}}),\tag{2}\] \[\label{2463} \rho{\boldsymbol{U}}=\sum_{i=1}^2\rho_i{\boldsymbol{U}}_i=\sum_{i=1}^2\int \mathrm{d}{\boldsymbol{v}}m_i{\boldsymbol{v}}f_i({\boldsymbol{v}}),\tag{3}\] \[\label{2464} nT=p=\sum_{i=1}^2 n_iT_i=\sum_{i=1}^2\frac{m_i}{d}\int \mathrm{d}{\boldsymbol{v}}V^2f_i({\boldsymbol{v}}),\tag{4}\] where \(\rho_i=m_in_i\) is the mass density of species \(i\), \(n=n_1+n_2\) is the total number density, \(\rho=\rho_1+\rho_2\) is the total mass density, \({\boldsymbol{V}}={\boldsymbol{v}}-{\boldsymbol{U}}\) is the peculiar velocity, and \(p\) is the hydrostatic pressure. Furthermore, the third equality of Eq.(4 ) defines the kinetic temperatures \(T_i\) of each species, which measure their mean kinetic energies. For inelastic collisions (\(\alpha_{ij}\neq 1\)), energy equipartition is in general broken and so \(T_i \neq T\) [59].

The collision operators conserve the particle number of each species and the total momentum, but the total energy is not conserved. These conditions lead to the following identities: \[\label{2465} \int\, d{\boldsymbol{v}} J_{ij}[{\boldsymbol{v}}|f_i,f_j]=0,\tag{5}\] \[\label{2466} \sum_{i=1}^2\sum_{j=1}^2\int \mathrm{d}{\boldsymbol{v}}m_{i}{\boldsymbol{v}} J_{ij}[{\boldsymbol{v}}|f_{i},f_{j}]=0,\tag{6}\] \[\label{2467} \sum_{i=1}^2\sum_{j=1}^2\int \mathrm{d}{\boldsymbol{v}}\frac{1}{2}m_{i}V^{2}J_{ij} [{\boldsymbol{v}}|f_{i},f_{j}]=-\frac{d}{2}nT\zeta.\tag{7}\] Here, \(\zeta\) is identified as the “cooling rate” due to inelastic collisions among all species. At a kinetic level, it is also convenient to introduce the “cooling rates” \(\zeta_i\) for the partial temperatures \(T_i\). They are defined as \[\label{2468} \zeta_i=\sum_{j=1}^2 \zeta_{ij}=-\sum_{j=1}^2 \frac{1}{dn_iT_i}\int \mathrm{d}{\boldsymbol{v}}m_iV^{2}J_{ij}[{\boldsymbol{v}}|f_{i},f_{j}],\tag{8}\] where the second equality defines the quantities \(\zeta_{ij}\). According to Eqs.@eq:2467 and 8 , the total cooling rate \(\zeta\) can be written in terms of the partial cooling rates \(\zeta_i\) as \[\label{2469} \zeta=T^{-1}\sum_{i=1}^2\; x_iT_i\zeta_i,\tag{9}\] where \(x_i=n_i/n\) is the concentration or mole fraction of species \(i\).

The macroscopic balance equations for the densities of mass, momentum and energy can be easily now derived from the constraints 5 –7 . They are given by \[D_{t}n_{i}+n_{i}\nabla \cdot {\boldsymbol{U}}+\frac{\nabla \cdot {\boldsymbol{j}}_{i}}{m_{i}} =0\;, \label{24610}\tag{10}\] \[D_{t}{\boldsymbol{U}}+\rho ^{-1}\nabla \cdot \mathsf{P}=0\;, \label{24611}\tag{11}\] \[D_{t}T-\frac{T}{n}\sum_{i=1}^2\frac{\nabla \cdot \mathbf{j}_{i}}{m_{i}}+\frac{2}{dn} \left( \nabla \cdot \mathbf{q}+\mathsf{P}:\nabla \mathbf{U}\right) =-\zeta T\;. \label{24612}\tag{12}\] In the above equations, \(D_{t}=\partial _{t}+{\boldsymbol{U}}\cdot \nabla\) is the material derivative, \[\mathbf{j}_{i}=m_{i}\int \mathrm{d}{\boldsymbol{v}}\,{\boldsymbol{V}}\,f_{i}({\boldsymbol{v}}) \label{24613}\tag{13}\] is the mass flux for species \(i\) relative to the local flow, \[\mathsf{P}=\sum_{i=1}^2\,\int \mathrm{d}{\boldsymbol{v}}\,m_{i}{\boldsymbol{V}}{\boldsymbol{V}}\,f_{i}({\boldsymbol{v}}) \label{24614}\tag{14}\] is the total pressure tensor, and \[\mathbf{q}=\sum_{i=1}^2\,\int \mathrm{d}{\boldsymbol{v}}\,\frac{1}{2}m_{i}V^{2}{\boldsymbol{V}} \,f_{i}({\boldsymbol{v}}) \label{24615}\tag{15}\] is the total heat flux. It must be remarked that the balance equations (10 )–(12 ) apply regardless of the details of the model for inelastic collisions considered. However, the influence of the collision model appears through the dependence of the cooling rate and the hydrodynamic fluxes on the coefficients of restitution.

2.1 Inelastic Maxwell models↩︎

On the other hand, the hydrodynamic equations 10 –12 do not constitute a closed set of differential equations for the hydrodynamic fields. To close them, one needs to solve the Boltzmann equation 1 to derive the corresponding constitutive equations for the fluxes and identify the explicit expressions of the transport coefficients. In the case of IHS, those expressions cannot be exactly determined and one has to resort to approximate solutions based on the use of Grad’s moment method [60]–[63] and/or the truncation of a Sonine polynomial expansion [3], [5]. As mentioned in the Introduction section, a possible way of obtaining exact forms for the transport coefficients is to consider IMM where the collision rate is independent of the relative velocity of the colliding spheres. For this interaction model, the Boltzmann collision operators \(J_{ij}^{\text{IMM}}[f_i,f_j]\) are [5], [26] \[J_{ij}^\text{IMM}[\mathbf{v}_1|f_i,f_j] =\frac{\nu_{\text{M},ij}}{n_jS_d} \int \mathrm{d}{\boldsymbol{v}}_{2}\int \mathrm{d}\widehat{\boldsymbol{\sigma }}\Big[ \alpha_{ij}^{-1}f_{i}({\boldsymbol{v}}_{1}'') f_{j}({\boldsymbol{v}}_{2}'')-f_{i}({\boldsymbol{v}}_{1})f_{j}({\boldsymbol{v}}_{2})\Big]. \label{24616}\tag{16}\] Here, \(\nu_{\text{M},ij}\) is an effective velocity-independent collision frequency for collisions of type \(i\)-\(j\) and \(S_d=2\pi^{d/2}/\Gamma(d/2)\) is the total solid angle in \(d\) dimensions. In Eq.@eq:24616 , in a binary collision the relationship between the pre-collisional velocities \((\mathbf{v}_1^{\prime\prime},\mathbf{v}_2^{\prime\prime})\) and the post-collisional velocities \((\mathbf{v}_1,\mathbf{v}_2)\) is \[\label{24617} {\boldsymbol{v}}_{1}^{\prime\prime}={\boldsymbol{v}}_{1}-\mu_{ji}\left( 1+\alpha_{ij}^{-1} \right)(\widehat{\boldsymbol{\sigma}}\cdot {\boldsymbol{g}}_{12})\widehat{\boldsymbol{\sigma}},\tag{17}\] \[\label{24617461} \quad {\boldsymbol{v}}_{2}^{\prime\prime}={\boldsymbol{v}}_{2}+\mu_{ij}\left( 1+\alpha_{ij}^{-1}\right) (\widehat{\boldsymbol{\sigma}}\cdot {\boldsymbol{g}}_{12})\widehat{\boldsymbol{\sigma}}\;,\tag{18}\] where \({\boldsymbol{g}}_{12}={\boldsymbol{v}}_1-{\boldsymbol{v}}_2\) is the relative velocity of the colliding pair, \(\widehat{\boldsymbol{\sigma}}\) is a unit vector directed along the centers of the two colliding spheres, and \(\mu_{ij}=m_i/(m_i+m_j)\). As for IHS [3], the scattering rules 17 and 18 yield \((\widehat{\boldsymbol{\sigma}}\cdot \mathbf{g}_{12}^{\prime\prime})=-\alpha_{ij}^{-1}(\widehat{\boldsymbol{\sigma}}\cdot \mathbf{g}_{12})\) where \(\mathbf{g}_{12}^{\prime\prime}={\boldsymbol{v}}_1^{\prime\prime}-{\boldsymbol{v}}_2^{\prime\prime}\). The collision frequencies \(\nu_{\text{M},ij}\) appearing in 16 can be seen as free parameters in the model. As usual, its dependence on the restitution coefficients \(\alpha_{ij}\) and the parameters of the mixture can be chosen to optimize the agreement with the results obtained from the Boltzmann equation for IHS. Of course, the choice is not unique and may depend on the property of interest.

As happens for elastic collisions [6], [34], the main advantage of using IMM instead of IHS is that a velocity moment of order \(k\) of the Boltzmann collision operator \(J_{ij}^\text{IMM}[f_i,f_j]\) only involves moments of order less than or equal to \(k\) of the distributions functions. This allows one to determine the Boltzmann collisional moments without the explicit knowledge of the velocity distribution functions. In particular, the first and second collisional moments of \(J_{ij}^\text{IMM}[f_i,f_j]\) are [64] \[\label{24619} \int \mathrm{d}{\boldsymbol{v}} m_i{\boldsymbol{V}}J_{ij}^\text{IMM}[f_i,f_j]=-\frac{\nu_{\text{M},ij}}{d\rho_j}\mu_{ji}(1+\alpha _{ij}) \left(\rho_j{\boldsymbol{j}}_i-\rho_i{\boldsymbol{j}}_j\right),\tag{19}\] \[\begin{align} \label{24620} \int \mathrm{d}{\boldsymbol{v}} m_i V_k V_\ell J_{ij}^\text{IMM}[f_i,f_j]&=& -\frac{\nu_{\text{M},ij}}{d\rho_j}\mu_{ji}(1+\alpha _{ij}) \Bigg\{2\rho_j P_{i,k \ell}-\left(j_{i,k}j_{j,\ell}+j_{j,k}j_{i,\ell}\right) \nonumber\\ & &-\frac{2}{d+2}\mu_{ji}(1+\alpha_{ij})\Bigg[\rho_jP_{i,k \ell}+\rho_iP_{j,k \ell}- \left(j_{i,k}j_{j,\ell}+j_{j,k}j_{i,\ell}\right) \nonumber\\ & &+ \Bigg(\frac{d}{2}\left(\rho_i p_j+\rho_jp_i\right)-{\boldsymbol{j}}_i\cdot {\boldsymbol{j}}_j\Bigg)\delta_{k \ell}\Bigg] \Bigg\}, \end{align}\tag{20}\] where \(p_i=n_i T_i=\text{Tr} (\mathsf{P}_i)/d\). The quantities \(\zeta_{ij}\) defined by Eq.@eq:2468 can be exactly obtained for IMM from Eq.@eq:24620 as \[\label{24621} \zeta_{ij}=\frac{2\nu_{\text{M},ij}}{d}\mu_{ji}(1+\alpha _{ij})\Bigg[1-\frac{\mu_{ji}}{2}(1+\alpha _{ij})\frac{\theta_i+\theta_j}{\theta_j}+\frac{\mu_{ji}(1+\alpha _{ij})-1}{d\rho_j p_i}{\boldsymbol{j}}_i\cdot {\boldsymbol{j}}_j\Bigg],\tag{21}\] where \[\label{24622} \theta_i=\frac{m_i T}{\overline{m} T_i}, \quad \overline{m}=\frac{m_1+m_2}{2}.\tag{22}\]

To optimize the agreement with the IHS results we adjust the collision frequencies \(\nu_{\text{M},ij}\) to obtain the same expression of \(\zeta_{ij}\) as the one found for IHS in the HCS [59]. However, given that the cooling rates are not exactly known for IHS, a good estimate for them can be achieved by considering the local equilibrium approximation for the velocity distribution functions \(f_i\), i.e., \[\label{24623} f_i({\boldsymbol{V}})\to f_{i,\text{M}}({\boldsymbol{V}})=n_i\left(\frac{m_i}{2\pi T_i}\right)^{d/2}\exp\left(-\frac{m_iV^2}{2T_i}\right).\tag{23}\] In this approximation, one has [5] \[\label{24624} \zeta_{ij}^{\text{IHS}}\to \frac{4\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)} n_j\mu_{ji}\sigma_{ij}^{d-1}\left(\frac{\theta_i+\theta_j}{\theta_i\theta_j}\right)^{1/2} (1+\alpha_{ij})\left[1-\frac{\mu_{ji}}{2}(1+\alpha_{ij}) \frac{\theta_i+\theta_j}{\theta_j}\right]\upsilon_\text{th},\tag{24}\] where \(\upsilon_\text{th}=\sqrt{2T/\overline{m}}\) is a thermal velocity defined in terms of the temperature \(T(t)\) of the mixture. Thus, according to Eqs.(21 ) and 24 , to get \(\zeta_{ij}=\zeta_{ij}^{\text{IHS}}\) one has to chose the collision frequencies \(\nu_{\text{M},ij}\) as \[\label{24625} \nu_{\text{M},ij}= \frac{2\pi^{(d-1)/2}}{\Gamma\left(\frac{d}{2}\right)}n_j \sigma_{ij}^{d-1}\left(\frac{2{T}_i}{m_i}+ \frac{2{T}_j}{m_j}\right)^{1/2}.\tag{25}\] Upon deriving Eq.@eq:24625 use has been made of the fact that the mass flux \({\boldsymbol{j}}_i\) vanishes in the HCS. In the remainder of this paper, we will take the choice (25 ) for \(\nu_{\text{M},ij}\).

2.2 Kinetic model for granular mixtures↩︎

Another different way of overcoming the mathematical intricacies of the Boltzmann collision operator for IHS is to consider a kinetic model. Here, as said in section 1, we consider the VGS kinetic model [39] for granular mixtures. To the best of our knowledge, this is the only kinetic model has been proposed so far in the granular literature for this sort of systems. The model is based on the equivalence between a system of elastic hard spheres subject to a nonconservative force proportional to the particle velocity with a gas of IHS [49]. This (approximate) mapping between a molecular hard sphere gas in the presence of a drag force with IHS allows to extend any kinetic model of molecular mixtures proposed in the literature to inelastic multicomponent gases. Here, the relaxation term appearing in the model reported in Ref.[39] has been chosen to be the well-known Gross-Krook (GK) kinetic model [40] proposed many years ago for studying transport properties of multicomponent molecular gases. In this sense, the kinetic model employed in this paper can be seen as a direct extension of the GK model [40] to granular mixtures.

In the tracer limit, within the VGS kinetic model for granular mixtures [39], the Boltzmann collision operators \(J_{22}[f_2,f_2]\) and \(J_{12}[f_1,f_2]\) are defined respectively, as \[\label{5461} J_{22}[f_2,f_2]\to K_{22}[f_2,f_2]=-\frac{1+\alpha_{22}}{2}\nu_{22}\left(f_2-f_{22}\right)+\frac{\epsilon_{22}}{2}\frac{\partial}{\partial \mathbf{V}}\cdot \Big[(\mathbf{v}-\mathbf{U}_2)\Big]f_2,\tag{26}\] \[\label{5462} J_{12}[f_1,f_2]\to K_{12}[f_1,f_2]=-\frac{1+\alpha_{12}}{2}\nu_{12}\left(f_1-f_{12}\right)+\frac{\epsilon_{12}}{2}\frac{\partial}{\partial \mathbf{V}}\cdot \Big[(\mathbf{v}-\mathbf{U}_1)\Big]f_1,\tag{27}\] where \[\label{5463} \nu_{ij}=\frac{4\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n_j\sigma_{ij}^{d-1}\left(\frac{2\widetilde{T}_i}{m_i}+ \frac{2\widetilde{T}_j}{m_j}\right)^{1/2}\tag{28}\] is an effective collision frequency for IHS, \(\widetilde{T}_i=T_i-(m_i/d)(\mathbf{U}_i-\mathbf{U})^2\), \[\label{5464} \epsilon_{ij}=\frac{1}{2}\nu_{ij}\mu_{ji}^2\left[1+\frac{m_i\widetilde{T}_j}{m_j \widetilde{T}_i}+\frac{3}{2d}\frac{m_i}{\widetilde{T}_i}\left({\boldsymbol{U}}_i-{\boldsymbol{U}}_j\right)^2\right](1-\alpha_{ij}^2),\tag{29}\] and \[\label{5465} f_{ij}(\mathbf{v})=n_i\left(\frac{m_i}{2\pi T_{ij}}\right)^{d/2} \exp\left[-\frac{m_i}{2T_{ij}}\left(\mathbf{v}-\mathbf{U}_{ij}\right)^2\right].\tag{30}\] In Eq.@eq:5465 , we have introduced the quantities \[\label{5466} \mathbf{U}_{ij}=\mu_{ij}\mathbf{U}_i+\mu_{ji}\mathbf{U}_j, \quad T_{ij}=\widetilde{T}_i+2\mu_{ij}\mu_{ji}\Bigg\{\widetilde{T}_j-\widetilde{T}_i+\frac{(\mathbf{U}_i-\mathbf{U}_j)^2}{2d} \left[m_j+ \frac{\widetilde{T}_j-\widetilde{T}_i}{\widetilde{T}_i/m_i+\widetilde{T}_j/m_j}\right]\Bigg\}.\tag{31}\] According to Eqs.@eq:24625 and 28 , note that \(\nu_{\text{M},ij}\neq \nu_{ij}\).

3 Rheology of a sheared granular mixture in the tracer limit. IMM↩︎

We consider the tracer limit (\(x_1\to 0\)) in a granular binary mixture. In this situation, since the concentration of the tracer species is negligibly small, the state of the excess granular gas is not perturbed by the presence of the tracer particles. Thus, the distribution function \(f_2\) of the excess granular gas obeys the nonlinear closed Boltzmann equation \[\label{24626} \partial_t f_2+{\boldsymbol{v}}\cdot \nabla f_{2}=J_{22}^{\text{IMM}}\left[f_{2},f_{2}\right] \;.\tag{32}\] Additionally, the collisions between tracer particles themselves can be also neglected in the kinetic equation of the distribution \(f_1\): \[\label{24627} \partial_t f_1+{\boldsymbol{v}}\cdot \nabla f_{1}=J_{12}^{\text{IMM}}\left[f_{1},f_{2}\right] \;.\tag{33}\]

We assume that the system (granular gas plus tracers) is under USF. At a macroscopic level, the USF state is characterized by constant densities \(n_2\simeq n\) and \(n_1\), a uniform granular temperature, and a linear velocity profile \[\label{3461} U_{1,k}=U_{2,k}=U_k=a_{k\ell}r_\ell, \quad a_{k\ell}=a \delta_{kx}\delta_{\ell y},\tag{34}\] \(a\) being the constant shear rate. This linear velocity profile assumes no boundary layer near the walls and is generated by the Lee-Edwards boundary conditions [65], which are simply periodic boundary conditions in the local Lagrangian frame moving with the mean flow velocity \(\mathbf{U}\). For elastic collisions, the temperature grows in time due to the viscous heating term (\(-aP_{xy}>0\)) and hence a steady state is not possible unless an external (artificial) force is introduced [66]. However, in the case of inelastic collisions, the temperature changes in time due to the competition between two (opposite) mechanisms: on the one hand, viscous (shear) heating (\(-aP_{xy}>0\)) and, on the other hand, energy dissipation in collisions (\(-\zeta T<0\)). A steady state is achieved when both mechanisms cancel each other and the fluid autonomously seeks the temperature at which the above balance occurs. Under stationary conditions, the balance equation (12 ) becomes \[\label{3462} aP_{2,xy}=-\frac{d}{2}\zeta_2 p_2,\tag{35}\] where we have taken into account that in the tracer limit, \(P_{xy}\simeq P_{2,xy}\), \(\zeta\simeq \zeta_2\) and \(p\simeq p_2=n_2 T_2\). According to Eq.@eq:3462 , it is quite apparent the intrinsic connection between the shear field and collisional dissipation in the system. Thus, the steady shear flow state characterized by Eq.@eq:3462 is inherently a non-Newtonian state [67] since the collisional cooling (which is fixed by the mechanical properties of the particles) sets the strength of the (reduced) velocity gradient in the steady state. This means that for given values of the shear rate \(a\) and the coefficient of restitution \(\alpha_{22}\), the steady state relation 35 gives the (reduced) shear rate \(a^*=a/\nu_{\mathrm{M},22}(T)\) as a unique function of the coefficient of restitution \(\alpha_{22}\).

At a microscopic level, the USF state becomes an homogeneous state when one refers the velocities of the particles to the local Lagrangian frame moving with the mean flow velocity \(\mathbf{U}\). In this frame, the distributions \(f_2\) and \(f_1\) adopt the forms \[\label{3462461} f_2(\mathbf{r}, \mathbf{v})=f_2(\mathbf{V}), \quad f_1(\mathbf{r}, \mathbf{v})=f_1(\mathbf{V}),\tag{36}\] where we recall that \(\mathbf{V}=\mathbf{v}-\mathbf{U}\) is the peculiar velocity. In the steady state, the corresponding set of Boltzmann equations 32 and 33 in the above Lagrangian frame become \[\label{3463} -a V_y\frac{\partial f_2}{\partial V_x}=J_{22}^{\text{IMM}}\left[f_{2},f_{2}\right], \quad -a V_y\frac{\partial f_1}{\partial V_x}=J_{12}^{\text{IMM}}\left[f_{1},f_{2}\right].\tag{37}\] Since the mass and heat fluxes vanish in the steady USF state, the relevant transport properties of the system are related to the pressure tensors \(P_{2,k \ell}\) and \(P_{1,k \ell}\) defined as \[\label{3464} P_{i,k \ell}=\int \mathrm{d}\mathbf{V}\; m_i V_k V_\ell f_i(\mathbf{V}).\tag{38}\]

Explicit expressions for the (reduced) nonzero elements of \(\mathsf{P}_2^*=\mathsf{P}_2/p_2\) and \(\mathsf{P}_1^*=\mathsf{P}_1/(x_1p_2)\) can be easily obtained from Eqs.@eq:3463 when one takes into account Eq.@eq:24620 with \(\mathbf{j}_i=\mathbf{0}\). In terms of the dimensionless quantities \[\label{3464461} \nu_\eta^*=\frac{(d+1-\alpha_{22})(1+\alpha_{22})}{d(d+2)},\tag{39}\] \[\zeta_2^*=\frac{\zeta_2}{\nu_{\text{M},22}}=\frac{1-\alpha_{22}^2}{2d},\] the nonzero elements of \(\mathsf{P}_2^*\) are [67] \[\label{3465} P_{2,yy}^*=P_{2,zz}^*=1-\frac{\zeta_2^*}{\nu_\eta^*}=\frac{d}{2}\frac{1+\alpha_{22}}{d+1-\alpha_{22}}, \quad P_{2,xx}^*=d-(d-1)P_{2,yy}^*=\frac{d}{2}\frac{d+3-(d+1)\alpha_{22}}{d+1-\alpha_{22}},\tag{40}\] \[\label{3466} P_{2,xy}^*=-\frac{P_{2,yy}^*}{\nu_\eta^*}a^*=-\frac{d^2(d+2)}{2(d+1-\alpha_{22})^2}a^*,\tag{41}\] \[a^*=\frac{a}{\nu_{\text{M},22}}=\sqrt{\frac{d}{2}\frac{\nu_\eta^*\zeta_2^*}{P_{2,yy}^*}}= \frac{d+1-\alpha_{22}}{d}\sqrt{\frac{1-\alpha_{22}^2}{2(d+2)}}.\] In the case of \(\mathsf{P}_1^*\), its nonzero elements are given by [64] \[\label{3467} P_{1,yy}^*=P_{1,zz}^*=-\frac{Y+X P_{2,yy}^*}{X_0}, \quad P_{1,xy}^*=\frac{a^*P_{1,yy}^*-XP_{2,xy}^*}{X_0}, \quad P_{1,xx}^*=d\gamma-(d-1)P_{1,yy}^*,\tag{42}\] where \(\gamma=T_1/T_2\) is the temperature ratio and \[\label{3468} Y=\frac{1}{d+2}\mu_{12}\mu_{21}\left(\frac{1+\theta}{\theta}\right)(1+\alpha_{12})^2 \nu_{\text{M},12}^*,\tag{43}\] \[\label{3468bis} X_0=-\frac{2}{d(d+2)}\mu_{21}(1+\alpha_{12})\left[d+2-\mu_{21}(1+\alpha_{12})\right] \nu_{\text{M},12}^*,\tag{44}\] \[\label{3469} X=\frac{2}{d(d+2)}\mu_{12}\mu_{21}(1+\alpha_{12})^2 \nu_{\text{M},12}^*.\tag{45}\] In Eqs.@eq:3468 and 45 , \(\theta=\theta_1/\theta_2=\mu/\gamma\), \(\mu=m_1/m_2\) is the mass ratio and \[\label{34610} \nu_{\text{M},12}^*=\frac{\nu_{\text{M},12}}{\nu_{\text{M},22}}=\left(\frac{\sigma_{12}}{\sigma_2}\right)^{d-1} \left(\frac{1+\theta}{2\theta}\right)^{1/2}.\tag{46}\] The temperature ratio is obtained by numerically solving the equation \[\label{34611} \frac{P_{1,xy}^*}{P_{2,xy}^*}=\frac{\gamma \zeta_1^*}{\zeta_2^*},\tag{47}\] where \[\label{34612} \zeta_1^*=\frac{2\nu_{\text{M},12}^*}{d}\mu_{21}(1+\alpha_{12})\Big[1-\frac{\mu_{21}}{2}(1+\alpha_{12})(1+\theta)\Big].\tag{48}\] Appendix 9 shows that the steady state solutions 40 –42 are indeed linearly stable solutions.

4 Tracer diffusion under USF: Results for IMM↩︎

We assume that the USF state characterized by the rheological properties 40 –42 is slightly perturbed by weak spatial gradients. Under these conditions, one expects that the mass flux of the intruder \(\mathbf{j}_1\) has nonzero contributions due to the existence of the gradients \(\nabla x_1\), \(\nabla p\) and \(\nabla T\). Additionally, the anisotropy induced by the strong shear flow in the velocity space gives rise to the presence of tensorial quantities (\(D_{k\ell}\), \(D_{p,k\ell}\) and \(D_{T,k\ell}\)) to describe the mass transport instead of the conventional scalar coefficients (\(D\), \(D_p\) and \(D_T\)) when the granular gas is in the HCS [11], [12]. The determination of the above tensors is the main goal of the present work.

As already made in Ref.[58], we have to start from the Boltzmann equation 33 with a general space and time dependence. Thus, we denote by \(U_{s,k}=a_{k\ell}r_\ell\) the mean flow velocity of the undisturbed USF state. However, when we disturb the USF state, the true mean flow velocity is \(U_k=U_{s,k}+\delta U_k\), where \(\delta U_k\) is a small perturbation to \(U_{s,k}\). In the frame moving with the (undisturbed) mean velocity \(\mathbf{U}_s\), Eq.@eq:24627 becomes \[\label{4461} \frac{\partial}{\partial t}f_1-aV_y\frac{\partial}{\partial V_x}f_1+\left({\boldsymbol{V}}+{\boldsymbol{U}}_s\right) \cdot \nabla f_1=J_{12}^\text{IMM}[f_1,f_2],\tag{49}\] where \(\mathbf{V}=\mathbf{v}-\mathbf{U}_s\) and the derivative \(\nabla f_1\) is taken at constant \({\boldsymbol{V}}\). The macroscopic balance equations associated with this disturbed USF state follows from the general equations (10 ), (11 ), and (12 ) when one takes into account that \(\mathbf{U}=\mathbf{U}_s+\delta \mathbf{U}\). The result is \[\label{4462} \partial_tn_1+\mathbf{U}_s\cdot \nabla n_1=-\nabla \cdot (n_1\delta \mathbf{U})- \frac{\nabla \cdot {\boldsymbol{j}}_1}{m_1},\tag{50}\] \[\label{4463} \partial_t\delta U_k+a_{k\ell}\delta U_\ell+(U_{s,\ell}+\delta U_\ell) \nabla_\ell \delta U_k=- \rho^{-1}\nabla_j P_{2,k\ell},\tag{51}\] \[\label{4464} \frac{d}{2}n\partial_tT+\frac{d}{2}n(\mathbf{U}_s+\delta \mathbf{U})\cdot \nabla T+aP_{xy}+\nabla \cdot {\boldsymbol{q}}+{\sf P}:\nabla \delta \mathbf{U}=-\frac{d}{2}p\zeta.\tag{52}\] Here, we recall that in the tracer limit \(\rho\simeq \rho_2=m_2 n_2\), \(n\simeq n_2\), \(T\simeq T_2\), \(\zeta\simeq \zeta_2\), and \(\mathsf{P}\simeq \mathsf{P}_2\). Additionally, the expressions of the cooling rate \(\zeta\), the mass flux \({\boldsymbol{j}}_1\), the pressure tensor \({\sf P}\), and the heat flux \({\boldsymbol{q}}\) are defined by Eqs.(7 ), (13 ), (14 ), and (15 ), respectively, with the replacement of \({\boldsymbol{V}}\) by the true peculiar velocity in the disturbed USF state \(\mathbf{c}=\mathbf{V}-\delta \mathbf{U}\).

Since the deviations from the USF are assumed to be small, our goal is to solve the Boltzmann kinetic equation 49 up to first order in the spatial gradients of the hydrodynamic fields \[\label{4465} A({\boldsymbol{r}},t)\equiv \{x_1({\boldsymbol{r}},t), p({\boldsymbol{r}}, t), T({\boldsymbol{r}}, t), \delta {\boldsymbol{U}}({\boldsymbol{r}},t)\}.\tag{53}\] In Eq.@eq:4465 , as in previous works on granular mixtures [11], [12], we represent the mass flux \(\mathbf{j}_1\) in terms of the spatial gradients of the fields \(x_1\), \(p=n T\), and \(T\). However, in contrast to previous studies [11], [12], [15], [16] where the granular gas is in the HCS, the system is strongly sheared and hence the conventional Chapman–Enskog method [6] cannot be applied. Thus, as in Ref.[58], we look for a solution to Eq.(49 ) by using a generalized Chapman–Enskog–like expansion where the velocity distribution function is expanded about a local shear flow reference state in terms of the small spatial gradients of the hydrodynamic fields relative to those of USF. This is the main new ingredient of the expansion. This type of generalized Chapman–Enskog expansion has been considered in the case of elastic gases to get the set of shear-rate dependent transport coefficients [34], [68] in a thermostatted shear flow problem and it has also been considered for monocomponent [69], [70] and multicomponent [58], [71] granular gases under USF.

Since the application of the generalized Chapman–Enskog method to Eq.@eq:4461 has been worked out in Ref.[58] for IHS, most of the mathematical steps involved in the determination of the first-order contribution \(\mathbf{j}_1^{(1)}\) to the mass flux for IMM will be omitted here. We refer the interested reader to Ref.[58] for more specific details. The first-order distribution \(f_1^{(1)}\) obeys the kinetic equation \[\begin{align} \label{4467461} \left(\partial_t^{(0)}-aV_y\frac{\partial}{\partial V_x}\right)f_1^{(1)}-J_{12}^\text{IMM}[f_1^{(1)},f_2^{(0)}]= {\boldsymbol{A}}_1\cdot \nabla x_1+{\boldsymbol{B}}_1\cdot \nabla p+{\boldsymbol{C}}_1\cdot \nabla T +\,\mathsf{D}_1:\nabla \delta {\boldsymbol{U}} +J_{12}^\text{IMM}[f_1^{(0)},f_2^{(1)}]. \end{align}\tag{54}\] In Eq.@eq:4467461 , the quantities \(\mathbf{A}_1(\mathbf{c})\), \(\mathbf{B}_1(\mathbf{c})\), \(\mathbf{C}_1(\mathbf{c})\), and \(\mathsf{D}_1(\mathbf{c})\) are defined by Eqs.(C.9)–(C.12), respectively, of Ref.[58] while the first-order distribution of the excess granular gas \(f_2^{(1)}(\mathbf{c})\) has the form [69], [70] \[\label{4468461} f_2^{(1)}(\mathbf{c})= {\boldsymbol{\mathcal{B}}}_{2}(\mathbf{c})\cdot \nabla p+{\boldsymbol{\mathcal{C}}}_{2}(\mathbf{c})\cdot \nabla T+{\sf {\mathcal{D}}}_{2}(\mathbf{c}):\nabla\delta {\boldsymbol{U}}.\tag{55}\] According to the right hand side of Eq.@eq:4467461 , the solution to this kinetic equation is \[f_1^{(1)}(\mathbf{c})={\boldsymbol{\mathcal{A}}}_{1}(\mathbf{c})\cdot \nabla x_1+ {\boldsymbol{\mathcal{B}}}_{1}(\mathbf{c})\cdot \nabla p+{\boldsymbol{\mathcal{C}}}_{1}(\mathbf{c})\cdot \nabla T+{\sf {\mathcal{D}}}_{1}(\mathbf{c}):\nabla\delta {\boldsymbol{U}}, \label{4468462}\tag{56}\] where the coefficients \({\boldsymbol{\mathcal{A}}}_{1}\), \({\boldsymbol{\mathcal{B}}}_{1}\), \({\boldsymbol{\mathcal{C}}}_{1}\), and \({\boldsymbol{\mathcal{D}}}_{1}\) are functions of the peculiar velocity and the hydrodynamic fields \(x_1\), \(p\), and \(T\).

To first-order, the mass flux \(\mathbf{j}_1^{(1)}\) is defined as \[\label{4466} \mathbf{j}_1^{(1)}=\int \mathrm{d}\mathbf{v}\; m_1 \mathbf{c} f_1^{(1)}(\mathbf{c}).\tag{57}\] Substitution of Eq.@eq:4468462 into Eq.@eq:4466 gives the expression \[\label{4467} j_{1,k}^{(1)}=-m_1 D_{k \ell}\nabla_\ell x_1-\frac{m_2}{T}D_{p,k \ell}\nabla_\ell p-\frac{\rho}{T}D_{T,k \ell}\nabla_\ell T,\tag{58}\] where \(\nabla_\ell\equiv \partial/\partial r_\ell\) and \[\label{4468} D_{k \ell}=-\int \mathrm{d}{\boldsymbol{c}}\,c_k\;{\mathcal{A}}_{1,\ell}({\boldsymbol{c}}),\tag{59}\] \[\quad D_{p,k \ell}=-\frac{Tm_1}{m_2}\int \mathrm{d}{\boldsymbol{c}}\,c_k\;{\mathcal{B}}_{1,\ell}({\boldsymbol{c}}),\] \[D_{T,k \ell}=-\frac{Tm_1}{\rho}\int \mathrm{d}{\boldsymbol{c}}\,c_k\;{\mathcal{C}}_{1,\ell}({\boldsymbol{c}}).\] In general, the set of generalized transport coefficients \(D_{k \ell}\), \(D_{p,k \ell}\), and \(D_{T,k \ell}\) are nonlinear functions of the shear rate and the coefficients of restitution \(\alpha_{22}\) and \(\alpha_{12}\). With respect to the conventional diffusion problem in the Navier–Stokes domain [11], there are important differences since the anisotropy induced by the presence of shear flow gives rise to new transport coefficients (such as the off-diagonal elements of the diffusion tensors) for the mass flux, reflecting broken symmetry. According to Eq.(58 ), the mass flux of the tracer particles is expressed in terms of a diffusion tensor \(D_{k \ell}\), a pressure diffusion tensor \(D_{p,k \ell}\), and a thermal diffusion tensor \(D_{T,k \ell}\).

The corresponding integral equations for the unknowns \({\boldsymbol{\mathcal{A}}}_{1}\), \({\boldsymbol{\mathcal{B}}}_{1}\), and \({\boldsymbol{\mathcal{C}}}_{1}\) defining the diffusion coefficients can be obtained by substituting Eqs.@eq:4468461 and 56 into Eq.@eq:4467461 and identifying coefficients of independent gradients. To achieve these equations, one needs to know the action of the operator \(\partial_t^{(0)}\) on the hydrodynamic fields. To lowest order in the expansion, the balance equations 50 –52 reduce to \[\label{4469} \partial_t^{(0)}x_1=0, \quad \partial_t^{(0)}\delta U_k=-a_{k \ell}\delta U_\ell,\tag{60}\] \[\label{44610} \partial_t^{(0)} \ln T=\partial_t^{(0)} \ln p=-\frac{2}{d p}a P_{xy}^{(0)}-\zeta.\tag{61}\] In the steady state (\(\partial_t^{(0)}T=\partial_t^{(0)}p=0\)), the expressions of \(\zeta\) and \(P_{xy}^{(0)}\) are given by Eqs.@eq:3464461 and 41 , respectively. To get the set of coupled linear integral equations for the unknowns, one has to take into account the contributions coming from the action of \(\partial_t^{(0)}\) on \(\nabla p\) and \(\nabla T\): \[\label{44611} \partial_t^{(0)} \nabla p=-\nabla \left(\frac{2}{d}a P_{xy}^{(0)}+p\zeta\right) =-\left(\frac{2a}{d}\frac{\partial P_{xy}^{(0)}}{\partial p}+ 2\zeta\right) \nabla p-\left(\frac{2a}{d}\frac{\partial P_{xy}^{(0)}}{\partial T}-\frac{1}{2} \frac{p\zeta}{T}\right) \nabla T,\tag{62}\] \[\begin{align} \label{44612} \partial_t^{(0)} \nabla T =-\nabla \left(\frac{2T}{d p}a P_{xy}^{(0)}+\zeta\right) =\left[\frac{2aT}{dp^2}\left(1-p\frac{\partial}{\partial p}\right)P_{xy}^{(0)}-\frac{T\zeta}{p}\right] \nabla p-\left[\frac{2a}{dp}\left(1+T\frac{\partial}{\partial T}\right)P_{xy}^{(0)}+\frac{1}{2}\zeta\right] \nabla T. \end{align}\tag{63}\] Additionally, within the Chapman–Enskog method [6], the distribution function \(f_1^{(1)}\) is qualified as a normal or hydrodynamic solution and hence, it depends on time through its dependence on the hydrodynamic fields \(A(\mathbf{r},t)\). As a consequence, \[\begin{align} \label{44613} \partial_t^{(0)}\mathcal{Z} =\frac{\partial \mathcal{Z}}{\partial p}\partial_t^{(0)} p+\frac{\partial \mathcal{Z}}{\partial T}\partial_t^{(0)} T+\frac{\partial \mathcal{Z}}{\partial \delta U_i}\partial_t^{(0)} \delta U_i = -\left(\frac{2}{d p}a P_{xy}^{(0)}+\zeta\right)\left(p\frac{\partial}{\partial p}+T\frac{\partial}{\partial T}\right) \mathcal{Z}+a_{k \ell}\delta U_\ell \frac{\partial \mathcal{Z}}{\partial c_k}, \end{align}\tag{64}\] where in the last step we have taken into account that \(\mathcal{Z}\) depends on \(\delta {\boldsymbol{U}}\) through \({\boldsymbol{c}}={\boldsymbol{V}}-\delta {\boldsymbol{U}}\).

As discussed in previous works [58], [71], the determination of the diffusion transport coefficients requires in general to numerically solve a set of differential equations. Thus, to get analytical forms for those coefficients, we restrict ourselves to linear deviations from the steady USF state. In this case, since the contributions to the mass flux are already of first order in the deviations from this steady state, we have to compute \(D_{ij}\), \(D_{p,ij}\), and \(D_{T,ij}\) to zero order in the deviations, namely, under steady state conditions. In this case, according to Eq.@eq:44610 , in the steady USF state \(\partial_t^{(0)}p=\partial_t^{(0)}T=0\) which implies that \((2/dp)a P_{2,xy}^{(0)}+\zeta=0\). Therefore, taking into account the results 62 –64 , the set of integral equations that the unknowns \(\mathcal{X}\equiv \{ {\boldsymbol{\mathcal{A}}}_{1}\), \({\boldsymbol{\mathcal{B}}}_{1}, {\boldsymbol{\mathcal{C}}}_{1} \}\) obey in the steady state is \[\label{44614} - a c_y\frac{\partial}{\partial c_x}{\boldsymbol{\mathcal{A}}}_{1}-J_{12}^\text{IMM}[{\boldsymbol{\mathcal{A}}}_{1},f_2^{(0)}]={\boldsymbol{A}}_{1},\tag{65}\] \[\begin{align} \label{44615} - a c_y\frac{\partial}{\partial c_x}{\boldsymbol{\mathcal{B}}}_{1} - \left(\frac{2a}{d}\partial_p P_{2,xy}^{(0)}+2\zeta\right){\boldsymbol{\mathcal{B}}}_{1} -J_{12}^\text{IMM}[{\boldsymbol{\mathcal{B}}}_{1},f_2^{(0)}] =\mathbf{B}_1-\Bigg[\frac{2aT}{d p^2}\left(1-p\partial_p\right)P_{2,xy}^{(0)}-\frac{T\zeta}{p}\Bigg] {\boldsymbol{\mathcal{C}}}_{1}+J_{12}^\text{IMM}[f_1^{(0)},{\boldsymbol{\mathcal{B}}}_{2}], \end{align}\tag{66}\] \[\begin{align} \label{44616} - a c_y\frac{\partial}{\partial c_x}{\boldsymbol{\mathcal{C}}}_{1} -\Bigg[\frac{2a}{d p}\left(1+T\partial_T\right)P_{2,xy}^{(0)}+\frac{1}{2}\zeta\Bigg]{\boldsymbol{\mathcal{C}}}_{1} -J_{12}^\text{IMM}[{\boldsymbol{\mathcal{C}}}_{1},f_2^{(0)}] =\mathbf{C}_1+\left(\frac{2a}{d}\partial_T P_{2,xy}^{(0)}-\frac{p\zeta}{2T}\right){\boldsymbol{\mathcal{B}}}_{1} +J_{12}^\text{IMM}[f_1^{(0)},{\boldsymbol{\mathcal{C}}}_{2}]. \end{align}\tag{67}\] It must be recalled that in Eqs.(65 )–(67 ) all the quantities are evaluated in the steady state. Henceforth, the calculations will be restricted to this particular condition. Moreover, as expected the structure of the integral equations (65 )–(66 ) is similar to the one derived for IHS [58], except by the replacement of the Boltzmann–Lorentz collision operators \(J_{12}^\text{IMM}[\mathcal{X},f_2^{(0)}]\) and \(J_{12}^\text{IMM}[f_1^{(0)},\mathcal{X}]\) by their corresponding IHS counterparts.

The main advantage of using IMM instead of IHS is that the collisional moments associated with the operator \(J_{12}^\text{IMM}[\mathcal{X},\mathcal{Y}]\) can be exactly computed. In particular, in the case of the mass flux, according to Eq.@eq:24619 one has the results \[\label{44617} \int \mathrm{d}{\boldsymbol{c}} \; m_1 c_k \left( \begin{array}{c} J_{12}^\text{IMM}[\mathcal{A}_{1,\ell},f_2^{(0)}]\\ J_{12}^\text{IMM}[\mathcal{B}_{1,\ell},f_2^{(0)}]\\ J_{12}^\text{IMM}[\mathcal{C}_{1,\ell},f_2^{(0)}] \end{array} \right) =\Omega \left( \begin{array}{c} m_1 D_{k \ell}\\ \frac{m_2}{T}D_{p,k \ell}\\ \frac{\rho}{T}D_{T,k \ell} \end{array} \right),\tag{68}\] \[\label{44618} \int \mathrm{d}{\boldsymbol{c}} \; m_1 c_k \left( \begin{array}{c} J_{12}^\text{IMM}[f_1^{(0)},\mathcal{B}_{2,\ell}]\\ J_{12}^\text{IMM}[f_1^{(0)},\mathcal{C}_{2,\ell}] \end{array} \right) = \left( \begin{array}{c} 0\\ 0 \end{array} \right),\tag{69}\] where we have taken into account that \(\mathbf{j}_1^{(0)}=\mathbf{j}_2^{(1)}=\mathbf{0}\) and \[\label{44622} \Omega=\frac{\nu_{\text{M},12}}{d}\mu_{21}(1+\alpha_{12}).\tag{70}\]

Thus, multiplying both sides of Eqs.@eq:44614 –67 by \(m_1 c_k\) and integrating over \(\mathbf{c}\), one gets the set of algebraic equations \[\label{44619} \left(a_{k\mu}+\Omega \delta_{k \mu}\right)D_{\mu \ell}=\frac{n}{\rho_1}P_{1,k \ell}^{(0)},\tag{71}\] \[\begin{align} \label{44620} \left(\frac{2a}{d}\partial_p P_{2,xy}^{(0)}+2\zeta\right)D_{p,k \ell}-\left(a_{k\mu}+\Omega \delta_{k \mu}\right)D_{p,\mu \ell}= \frac{T}{m_2}\left(\frac{\rho_1}{\rho}\partial_p P_{2,k\ell}^{(0)}-\partial_p P_{1,k\ell}^{(0)}\right) + \Bigg[\frac{2a}{d p}\left(1-p\partial_p\right)P_{2,xy}^{(0)}-\zeta\Bigg]D_{T,k\ell}, \end{align}\tag{72}\] \[\begin{align} \label{44621} \Bigg[\frac{2a}{d p}\left(1+T\partial_T\right)P_{2,xy}^{(0)}+\frac{1}{2}\zeta\Bigg]D_{T,k\ell}-\left(a_{k\mu}+\Omega \delta_{k \mu}\right)D_{T,\mu \ell}=\frac{T}{\rho}\left(\frac{\rho_1}{\rho}\partial_T P_{2,k\ell}^{(0)}-\partial_T P_{1,k\ell}^{(0)}\right) -\left(\frac{2a}{d n}\partial_T P_{2,xy}^{(0)}-\frac{\zeta}{2}\right)D_{p,k\ell}. \end{align}\tag{73}\] In Eqs.@eq:44620 and 73 , the derivatives of the pressure tensors with respect to the hydrostatic pressure \(p\) and/or the temperature \(T\) can be conveniently written in terms of the derivatives with respect to the (reduced) shear rate \(a^*\) when one takes into account that \(P_{2,k\ell}^{(0)}\) and \(P_{1,k\ell}^{(0)}\) depend on \(p\) and \(T\) explicitly and also through their dependence on the reduced shear rate \(a^*(p,T)\propto \sqrt{T}/p\). Thus, \[\label{44622461} \partial_p P_{2,k\ell}^{(0)}=\partial_p p P_{2,k\ell}^{*}(a^*)=\left(1-a^*\frac{\partial}{\partial a^*}\right)P_{2,k\ell}^{*}(a^*),\tag{74}\] \[\partial_T P_{2,k\ell}^{(0)}=\partial_T p P_{2,k\ell}^{*}(a^*)=\frac{p}{2T}a^*\frac{\partial}{\partial a^*}P_{2,k\ell}^{*}(a^*),\] \[\label{44623} \partial_p P_{1,k\ell}^{(0)}=\partial_p x_1 p P_{1,k\ell}^{*}(a^*)=x_1\left(1-a^*\frac{\partial}{\partial a^*}\right)P_{1,k\ell}^{*}(a^*),\tag{75}\] \[\partial_T P_{1,k\ell}^{(0)}=\partial_T x_1 p P_{1,k\ell}^{*}(a^*)=\frac{x_1 p}{2T}a^*\frac{\partial}{\partial a^*}P_{1,k\ell}^{*}(a^*).\] The derivatives of \(P_{2,k\ell}^{*}\) and \(P_{1,k\ell}^{*}\) with respect to \(a^*\) in the steady state for IMM are obtained in the Appendix 10.

As expected from the results derived for IHS [58], the coefficients \(D_{k\ell}\) obey an autonomous set of algebraic equations whose solution is \[\label{44624} \mathsf{D}=\frac{p}{m_1}\left(\mathsf{a}+\boldsymbol{\Omega}\right)^{-1}\cdot \mathsf{P}_1^*,\tag{76}\] where we recall that the tensors \(a_{k\ell}=a\delta_{kx}\delta_{\ell y}\) and \(\Omega_{k\ell}=\Omega \delta_{k\ell}\). The coefficients \(D_{p,k\ell}\) and \(D_{T,k\ell}\) obey a set of algebraic coupled equations. These equations can be easily solved.

For elastic collisions (\(\alpha_{22}=\alpha_{12}=1\)), Eqs.@eq:3465 –48 lead to \(a^*=0\), \(\zeta_2^*=0\), \(P_{2,k\ell}^{*}=P_{1,k\ell}^{*}=\delta_{k\ell}\), \(\gamma=1\), and \[\label{44625} \Omega\equiv \Omega^{\text{el}}=\frac{4\pi^{(d-1)/2}}{d\Gamma\left(\frac{d}{2}\right)}n \sigma_{12}^{d-1}\sqrt{\frac{2Tm_2}{m_1(m_1+m_2)}}.\tag{77}\] In this limiting case, the solution to Eqs.@eq:44619 –73 is \(D_{k\ell}=D^{\text{el}} \delta_{k\ell}\), \(D_{p,k\ell}=D_p^{\text{el}} \delta_{k\ell}\), and \(D_{T,k\ell}=D_T^{\text{el}} \delta_{k\ell}\), where \[\label{44626} D^{\text{el}}=\frac{p}{m_1 \Omega^{\text{el}}}, \quad D_p^{\text{el}}=x_1\left(1-\frac{m_1}{m_2}\right)\frac{T}{m_2 \Omega^{\text{el}}}, \quad D_T^{\text{el}}=0.\tag{78}\] As expected, the expressions 78 agree with those obtained in the Navier–Stokes domain for a molecular mixture of Maxwell molecules in the tracer limit [6]. Additionally, when the tracer particles are mechanically equivalent to the particles of the granular gas (i.e., \(\sigma_1=\sigma_2\), \(m_1=m_2\), and \(\alpha_{22}=\alpha_{12}\)), Eqs. 72 –73 lead to \(D_{p,k\ell}=D_{T,k\ell}=0\). In this limiting case, the mass flux obeys a generalized Fick’s law given by \[\label{44626461} j_{1,k}^{(1)}=-m_1 D_{k\ell} \frac{\partial x_1}{\partial r_\ell}.\tag{79}\] The dependence of the self-diffusion tensor \(D_{k\ell}\) on the coefficient of restitution has been analyzed by computer simulations for very dense granular fluids [56], [72].

5 Tracer diffusion under USF from a kinetic model of granular mixtures↩︎

To complement the results derived for IMM we consider here the kinetic model (referred to as the VGS model) defined by Eqs.@eq:5461 and 27 . First, the rheological properties in the USF state can be easily determined from Eqs.@eq:3463 when one makes the replacements \(J_{22}^\text{IMM}[f_2,f_2]\to K_{22}[f_2,f_2]\) and \(J_{12}^\text{IMM}[f_1,f_2]\to K_{12}[f_1,f_2]\). The nonzero elements of \(P_{2,ij}^*\) and \(P_{1,ij}^*\) are [50] \[\label{5468} P_{2,yy}^*=P_{2,zz}^*=\frac{1+\alpha_{22}}{1+\alpha_{22}+\epsilon_{22}^*}=\frac{2}{3-\alpha_{22}}, \quad P_{2,xx}^*=d-(d-1)P_{2,yy}^*=\frac{d}{2}\frac{d+3-(d+1)\alpha_{22}}{d+1-\alpha_{22}},\tag{80}\] \[\label{5469} P_{2,xy}^*=-\frac{2P_{2,yy}^*}{1+\alpha_{22}+2\epsilon_{22}^*}a^*=-\frac{8a^*}{(3-\alpha_{22})^2(1+\alpha_{22})}, \quad\tag{81}\] \[a^*=\frac{a}{\nu_{22}}=-\frac{d}{2}\frac{\zeta_2^*}{P_{2,xy}^*}= \frac{(3-\alpha_{22})(1+\alpha_{22})}{8}\sqrt{d(1-\alpha_{22})},\] \[\label{54610} P_{1,yy}^*=\frac{2}{(1+\alpha_{12})\nu_{12}^*+2\epsilon_{12}^*}, \quad P_{1,xy}^*=-\frac{2a^*}{(1+\alpha_{12})\nu_{12}^*+2\epsilon_{12}^*}P_{1,yy}^*, \quad P_{1,xx}^*=\gamma-(d-1)P_{1,yy}^*.\tag{82}\] In Eqs.@eq:5469 and 82 we have introduced the quantities \[\nu_{12}^*=\frac{\nu_{12}}{\nu_{22}}=\left(\frac{\sigma_{12}}{\sigma_2}\right)^{d-1}\sqrt{\frac{1+\theta}{2\theta}},\] \[\label{54610461} \epsilon_{22}^*=\frac{\epsilon_{22}}{\nu_{22}}=\frac{1-\alpha_{22}^2}{4}, \quad \epsilon_{12}^*=\frac{\epsilon_{12}}{\nu_{22}}=\frac{1}{2}\nu_{12}^*\mu_{21}^2(1+\theta)(1-\alpha_{12}^2).\tag{83}\] The temperature ratio \(\gamma=T_1/T_2\) can be obtained from the relationship 47 where \(\zeta_2^*\) must be replaced by \(\epsilon_{22}^*\) and \(\zeta_1^*\) is also given by Eq.@eq:34612 .

In the sheared diffusion problem, the first-order distribution function \(f_1^{(1)}\) obeys the kinetic equation 54 , except that \(J_{12}^\text{IMM}[f_1^{(0)},f_2^{(1)}]\to K_{12}[f_1^{(0)},f_2^{(1)}]=0\) and the Boltzmann–Lorentz collision operator \(J_{12}^\text{IMM}[f_1^{(1)},f_2^{(0)}]\) must be replaced by the operator \[\label{54611} K_{12}[f_1^{(1)},f_2^{(0)}]=-\frac{1+\alpha_{12}}{2}\nu_{12}\Big(f_1^{(1)}-f_{12}^{(1)}\Big)+\frac{\epsilon_{12}}{2}\Big( \frac{\partial}{\partial \mathbf{c}}\cdot \mathbf{c}f_1^{(1)}-\frac{\mathbf{j}_1^{(1)}}{\rho_1}\cdot \frac{\partial}{\partial \mathbf{c}} f_1^{(0)}\Big),\tag{84}\] where \[\label{54612} f_{12}^{(1)}(\mathbf{c})=\frac{\mu_{12}}{n_1T_{12}}\mathbf{c}\cdot \mathbf{j}_1^{(1)}f_{12}^{(0)}(\mathbf{c}), \quad f_{12}^{(0)}(\mathbf{c})=n_1\left(\frac{m_1}{2\pi T_{12}}\right)^{d/2} \exp\left(-\frac{m_1}{2T_{12}}c^2\right).\tag{85}\] With these replacements, it is easy to see that the diffusion transport coefficients of the BGK-type model are the solutions of the algebraic equations 71 –73 except that \(\Omega\) must be replaced by \[\label{54614} \Omega'=\frac{1+\alpha_{12}}{2}\mu_{21}\nu_{12}.\tag{86}\] Note that in Eqs. 71 –73 the derivatives of the pressure tensors \(P_{2,k \ell}^*\) and \(P_{1,k \ell}^*\) with respect to \(a^*\) in the steady state obtained from the kinetic model differ from those obtained from IMM and IHS. These derivatives are also displayed in the Appendix 10.

6 Comparison with the Boltzmann results for IHS↩︎

In this section, we want to compare the results obtained from both the Boltzmann equation for IMM and the VGS kinetic model with those previously derived from the IHS [57], [58]. In dimensionless form, the rheological properties and the diffusion coefficients of the system (granular gas plus tracer particles) depend on five quantities: the diameter \(\sigma_1/\sigma_2\) and mass \(m_1/m_2\) ratios, the coefficients of restitution \(\alpha_{22}\) and \(\alpha_{12}\), and the dimensionality \(d\) of the system. Thus, given that the parameter space of the system is large, henceforth we consider a three-dimensional gas (\(d=3\)) and a common coefficient of restitution \(\alpha_{22}=\alpha_{12}\equiv \alpha\). This reduces the parameter space to three quantities: \(\sigma_1/\sigma_2\), \(m_1/m_2\), and \(\alpha\).

First, we consider the rheological properties. In the case of the excess granular gas, figure [fig1] shows \(P_{2,k\ell}^*\) versus \(\alpha_{22}\). We have also included computer simulation results [73] obtained by numerically solving the (inelastic) Boltzmann equation for IHS by means of the direct simulation Monte Carlo (DSMC) method [74]. Comparing the three approaches, we see that the quantitative discrepancies between the VGS model’s theoretical predictions and the approximate IHS results are larger than those found for the IMM results. In fact, the theoretical results for IHS and IMM disagree very little; this difference increases slightly with inelasticity. Additionally, we observe excellent agreement between the Boltzmann theory for both interaction models and Monte Carlo simulations, even in the case of strong dissipation. Similar conclusions are reached for the reduced pressure tensor, \(P_{1,k\ell}^*\). Figure 2 illustrates this behavior, showing the dependence of \(P_{1,k\ell}^*\) on \(\alpha\) for the case \(\sigma_1/\sigma_2=1\) and \(m_1/m_2=0.5\). The good agreement between the IHS and IMM results is apparent again, especially in the case of the shear stress \(P_{1,xy}^*\), which is the most relevant rheological property in a shear flow problem as it defines the non-Newtonian shear viscosity. Although the VGS model predictions are good qualitatively, they exhibit larger discrepancies with the IHS results than the IMM predictions.

We analyze now the dependence of the diffusion coefficients \(T_{k\ell}\equiv \left\{D_{k\ell}, D_{p,k\ell},D_{T,k\ell}\right\}\) on the (common) coefficient of restitution \(\alpha\). According to the results derived in sections 4 and 5, \(T_{xz}=T_{zx}=T_{yz}=T_{zy}=0\) in agreement with the symmetry of the shearing field applied to the system. Thus, in a three-dimensional gas, there are five relevant elements of the tensors \(T_{k\ell}\): the three diagonal (\(T_{xx}\), \(T_{yy}\), and \(T_{zz}\)) and two off-diagonal elements (\(T_{xy}\) and \(T_{yx}\)). The results obtained here for IMM and from the kinetic model show that in general \(T_{xx}\neq T_{yy}=T_{zz}\) and \(T_{xy}\neq T_{yx}\). However, in the case of IHS, the approximate results derived in Ref.[58] show that \(T_{yy}\neq T_{zz}\) although the difference between both elements is in general very small.

Given that we are interested in this paper to assess the dependence of the coefficients \(T_{k\ell}\) on inelasticity, we have scaled them with respect to their elastic values except in the case of \(D_{T,k\ell}\) since this coefficient vanishes for elastic collisions. Thus, we define the dimensionless coefficients \(D_{k\ell}^*=D_{k\ell}/D^\text{el}\) and \(D_{p,k\ell}^*=D_{p,k\ell}/D_p^\text{el}\) while \(D_{T,k\ell}^*=D_{T,k\ell}/(x_1 T/m_2\nu)\). The effective collision frequency \(\nu=\nu_{\text{M},22}\) for IMM, \(\nu=\nu_{22}\) for the VGS model, and \(\nu=(2/(d+2))\nu_{\text{M},22}\) for IHS. The relevant elements of \(D_{k\ell}^*\), \(D_{p,k\ell}^*\), and \(D_{T,k\ell}^*\) are plotted in figures [fig3], [fig4], and [fig5], respectively, as functions of \(\alpha\) for the mixture \(\sigma_1/\sigma_2=1\) and \(m_1/m_2=0.5\). As expected, we observe that in general the influence of inelasticity on mass transport is quite significant regardless of the approximation used. Additionally, the anisotropy of the system (as measured by the differences \(|T_{xx}-T_{yy}|\) and \(|T_{zz}-T_{yy}|\)) is much important in the shear flow \(xy\)-plane. In fact, while \(T_{zz}-T_{yy}=0\) for IMM and kinetic model results, we observe that \(T_{zz}\simeq T_{yy}\) for IHS. With respect to the comparison between IHS, IMM, and kinetic model, it is quite apparent that although the predictions of the kinetic model reproduce qualitatively well the results obtained for IHS, significant discrepancies between both approaches are found for strong dissipation. As occurs for the rheology, the agreement between IMM and IHS is better than the one found between IHS and the kinetic model. In fact, the good agreement obtained here between IHS and IMM can justify the use of IMM as a reliable model to unveil in a clean way the effect of inelasticity on transport in real granular flows.

7 Thermal diffusion segregation↩︎

The knowledge of the diffusion transport coefficients allows us to apply the present theoretical results to the problem of segregation of tracer particles in a sheared granular gas. Segregation of dissimilar species in a granular mixture is likely one of the most relevant problems in granular flows not only from a fundamental point of view, but also from a more practical perspective. The problem was already studied in the context of IHS [75]. Our objective here is to revisit the problem by employing the results derived for IMM and from the VGS kinetic model.

Within the different mechanisms involved in the segregation, thermal diffusion segregation is perhaps one of the most studied in the literature [76], [77]. Thermal diffusion is produced by a thermal gradient which induces a relative motion within a mixture. This motion generates diffusion processes due to the concentration gradients created by the presence of thermal gradient. A steady state is achieved when thermal diffusion is balanced by the mixing effect coming from ordinary diffusion.

However, due to the anisotropy in velocity space induced by the shear flow, a complete description of the segregation problem requires to introduce a thermal diffusion tensor to characterize segregation in the different directions. Here, as in Ref.[75], for the sake of simplicity we consider a situation where the temperature gradient is perpendicular to the shear flow plane. Thus, for a three-dimensional system, \(\partial_x T=\partial_y T=0\) but \(\partial_z T\neq 0\). Additionally, since we are interested in a situation where the hydrodynamic equations 50 –52 admit a steady solution, we also assume that \(\delta \mathbf{U}=\mathbf{0}\). Under these conditions, the amount of segregation parallel to the \(z\)-axis can be measured by the thermal diffusion factor \(\Lambda_z\) defined as \[\label{7461} -\Lambda_z \frac{\partial \ln T}{\partial z}=\frac{\partial \ln x_1}{\partial z}.\tag{87}\] If we assume that the thermal gradient is directed downwards (\(\partial_z T<0\)), when \(\Lambda_z>0\) (\(\Lambda_z<0\)) the tracer particles tend to accumulate near the cold (hot) wall since \(\partial_z \ln x_1>0\) (\(\partial_z \ln x_1<0\)).

Let us write \(\Lambda_z\) in terms of the pressure tensors \(P_{zz}\) and \(P_{1,zz}\) as well as the diffusion coefficients \(D_{zz}\), \(D_{p,zz}\) and \(D_{T,zz}\). First, when only gradients along the \(z\)-axis exist, the momentum balance equation 51 yields \[\label{7462} \frac{\partial P_{zz}}{\partial z}=0.\tag{88}\] Since \(P_{zz}=pP_{zz}^*(a^*)\), then \[\label{7463} \frac{\partial P_{zz}}{\partial z}=\frac{\partial p}{\partial z}\left(1-a^*\partial_{a^*}\right)P_{zz}^*+ \frac{p}{2T}\frac{\partial T}{\partial z}a^*(\partial_{a^*}P_{zz}^*),\tag{89}\] where \(\partial_{a^*}P_{zz}^*=\partial_{a^*}P_{yy}^*\equiv \Delta_{yy}\) is given by Eq.@eq:b4 . Using Eqs.@eq:7462 and 89 , one gets a relationship between \(\partial_z p\) and \(\partial_z T\): \[\label{7464} \frac{\partial \ln p}{\partial z}=-\frac{1}{2}\frac{a^*(\partial_{a^*}P_{zz}^*)}{P_{zz}^*-a^*(\partial_{a^*}P_{zz}^*)}\frac{\partial \ln T}{\partial z}.\tag{90}\] In addition, according to the balance equation 50 , in the steady state with \(\delta \mathbf{U}=\mathbf{0}\) then \(j_{1,z}^{(1)}=0\). The constitutive equation for \(j_{1,z}^{(1)}\) is \[\label{7465} j_{1,z}^{(1)}=-m_1 D_{zz}\partial_z x_1-\frac{m_2}{T}D_{p,zz}\partial_z p-\frac{\rho}{T}D_{T,zz}\partial_z T.\tag{91}\] The condition \(j_{1,z}^{(1)}=0\) leads to \[\label{7466} \frac{\partial \ln x_1}{\partial z}=-\frac{\rho}{m_1 x_1}\left(\frac{D_{p,zz}}{D_{zz}}\frac{\partial \ln p}{\partial z} +\frac{D_{T,zz}}{D_{zz}}\frac{\partial \ln T}{\partial z}\right)= -\Bigg[(1-\mu)\frac{D_{p,zz}^*}{D_{zz}^*}\frac{\partial \ln p}{\partial z}+\Omega^{\text{el}*} \frac{D_{T,zz}^*}{D_{zz}^*}\frac{\partial \ln T}{\partial z}\Bigg].\tag{92}\] Here, \(\mu=m_1/m_2\) is the mass ratio and \(\Omega^{\text{el}*}=\Omega^{\text{el}}/\nu\), where the value of \(\nu\) depends on the approach followed.

From Eqs.@eq:7464 and 92 one easily gets the expression of the thermal diffusion factor \(\Lambda_z\) as \[\label{7468} \Lambda_z=\frac{\Omega^{\text{el}*}D_{T,zz}^*-\frac{1}{2}(1-\mu)a^*\Delta_{yy}\left(P_{zz}^*-a^*\Delta_{yy}\right)^{-1} D_{p,zz}^*}{D_{zz}^*}.\tag{93}\] The condition \(\Lambda_z=0\) gives the criterion for the upwards/downwards segregation transition.

In accordance to Eq.@eq:44624 the diffusion coefficient \(D_{zz}^*\) is positive. As a consequence, the marginal segregation curve (\(\Lambda_z=0\)) is obtained from the condition \[\label{7469} \Omega^{\text{el}*}D_{T,zz}^*-\frac{1}{2}(1-\mu)a^*\Delta_{yy}\left(P_{zz}^*-a^*\Delta_{yy}\right)^{-1} D_{p,zz}^*=0.\tag{94}\] For elastic collisions (\(\alpha_{22}=\alpha_{12}=1\)), \(D_{T,zz}^*=0\) and Eq.@eq:3466 for IMM or Eq.@eq:5469 for the VGS kinetic model yields \(a^*=0\). Thus, the condition 94 is trivially satisfied for any value of the mass and diameter ratios and so, no segregation occurs in this limiting case as expected. On the other hand, for inelastic collisions but the particles are mechanically equivalent, \(D_{p,k\ell}^*=D_{T,k\ell}^*=0\) and so, \(\Lambda_z=0\) for any value of the coefficient of restitution. This is the expected result since both species are indistinguishable.

For inelastic collisions, the zero contour of \(\Lambda_z\) exhibits a complex nonlinear dependence on the parameter space of the system. Thus, as did in section 6, for the sake of simplicity we take a three-dimensional granular gas (\(d=3\)) in the case of a common coefficient of restitution (\(\alpha_{22}=\alpha_{12}\equiv \alpha\)). The marginal segregation curve \(\Lambda_z=0\) separates regions of \(\Lambda_z>0\) (upwards segregation) and \(\Lambda_z<0\) (downwards segregation). At a fixed value of \(\alpha\), the points lying on the zero contour correspond to values of the diameter and mass ratios for which the intruder does not segregate in a sheared granular gas.

As an illustration, figure [fig6] shows the phase diagram in the \(\left(\sigma_1/\sigma_2, m_1/m_2\right)\)-plane for two different values of the (common) coefficient of restitution \(\alpha\). We compare the theoretical predictions for the marginal segregation curve \(\Lambda_z=0\) obtained previously for IHS in Ref.[75] with those derived here for IMM and from the VGS kinetic model. As previously mentioned, all curves pass through the point \((1,1)\) because it corresponds to the limiting case of mechanically equivalent particles. The three approaches show that, for \(\sigma_1<\sigma_2\), the main effect of inelasticity (or equivalently, the reduce shear rate \(a^*\)) is to enlarge the size of the downwards segregation region. The opposite occurs when the tracer particles are larger than the particles of the granular gas. In general, we see that the tracer particles tend to move toward hotter regions since upwards segregation occupies most of the system’s parameter space. Regarding the comparison of the three approaches, the results derived from IMM qualitatively agree well with the IHS results. Surprisingly, however, the segregation results obtained from the VGS model agree better with the IHS results than with the IMM results. This good agreement is especially noticeable when \(\sigma_1>\sigma_2\).

8 Concluding remarks↩︎

In this paper, we have analyzed the diffusion of tracer particles immersed in a sheared granular gas. Under these conditions, the mass flux \(\mathbf{j}_1^{(1)}\) is defined in terms of the tracer diffusion tensor \(D_{k\ell}\) (which couples \(\mathbf{j}_1^{(1)}\) with the concentration gradient \(\nabla x_1\)), the pressure diffusion tensor \(D_{p,k\ell}\) (which couples \(\mathbf{j}_1^{(1)}\) with the pressure gradient \(\nabla p\)), and the thermal diffusion tensor \(D_{T,k\ell}\) (which couples \(\mathbf{j}_1^{(1)}\) with the temperature gradient \(\nabla T\)). These tensorial quantities were evaluated years ago in the context of the Boltzmann equation for IHS [57], [58]. However, due to the intricate mathematical structure of the Boltzmann collision operator for IHS, the results obtained in Refs.[57], [58] involve several (uncontrolled) approximations at different stages of the derivation. Here, we revisit this problem by considering two different, complementary approaches that allow us to achieve exact results. First, we maintain the structure of the Boltzmann collision operators but consider a different interaction model: the so-called IMM, in which the collision rate of colliding spheres is independent of their relative velocity. This simplification enables us to obtain exact expressions for the rheological properties of the system (granular gas plus tracer particles), as well as the diffusion tensors. As a second approach, we keep the IHS interaction model but replace the true Boltzmann collision operators with simpler mathematical terms that retain their relevant physical properties. In this context, we consider the VGS kinetic model [39] proposed years ago for granular mixtures.

As in Ref.[58] for IHS, the diffusion tensors are obtained by solving the Boltzmann equation (or the VGS model) for tracer particles using a generalization of the Chapman–Enskog method [6] for far-from-equilibrium states. Since the granular gas is subjected to a strong shear rate, non Newtonian effects are relevant for finite inelasticity. Thus, the reference state (the zeroth-order distribution \(f_1^{(0)}\)) in the perturbation method is the shear flow distribution, not the local equilibrium distribution. Additionally, since collisional cooling cannot compensate for viscous heating locally, \(f_1^{(0)}\) is in general a time-dependent distribution even when the gas is slightly perturbed from the USF. Once the linear integral equations verifying the diffusion tensors are obtained, we restrict to steady state conditions and so, the reduced shear rate \(a^*\) is coupled to the coefficient of restitution \(\alpha_{22}\) which characterizes the inelasticity of grain-grain collisions. The consideration of the steady state allows us to achieve analytical, exact expressions for \(D_{k\ell}\), \(D_{p,k\ell}\), and \(D_{T,k\ell}\). These tensors depend nonlinearly on the diameter ratio, \(\sigma_1/\sigma_2\), the mass ratio, \(m_1/m_2\), and the coefficients of restitution \(\alpha_{22}\) and \(\alpha_{12}\) (which characterizes the inelasticity of tracer-grain collisions).

A comparison of the exact theoretical results of the IMM and VGS models with the approximate IHS results generally shows good qualitative agreement. At a more quantitative level, we observe excellent agreement in some cases (e.g., rheological properties) and reasonably good agreement in others (e.g., diffusion tensors), especially in the case of IMM. It is quite apparent that to confirm the reliability of the predictions offered by IMM and VGS model one should compare them with computer simulation results for IHS. At the level of rheology, the results reported here (see figure [fig1]) and in Ref.[2] for IMM clearly demonstrate the accuracy of this interaction model to capture the dependence of the pressure tensors on the coefficients of restitution. The lack of simulation data for the diffusion coefficients in the low-density regime prevents a comparison between theory and simulation. We hope that this paper will encourage simulators to perform simulations and confirm the results reported here.

As a nice application of the results exposed in this paper, the segregation of tracer particles in a sheared granular gas has been analyzed. The relative motion of the tracers with respect to the particles of the gas is caused by the presence of a temperature gradient. Here, for the sake of simplicity, we have assumed that the thermal gradient \(\partial_z T\) is perpendicular to the shear flow \(xy\)-plane. Under these conditions, the amount of segregation in the \(z\)-direction is measured by the thermal diffusion factor \(\Lambda_z\), defined in Eq.@eq:7461 . The condition of zero thermal diffusion (\(\Lambda_z=0\)) gives the segregation criterion for the transition from upwards segregation (regions where \(\Lambda_z>0\)) to downwards segregation (regions where \(\Lambda_z<0\)). A comparison with previous results obtained for IHS [75] (see figure [fig6]) shows reasonable agreement in general, especially for the VGS model when the tracer particles are larger than the granular gas particles.

In conclusion, the results found here give evidence of the accuracy of both IMM and the VGS kinetic model for studying far-from-equilibrium situations in granular flows, where using the original Boltzmann equation for IHS turns out to be quite intricate. Additionally, using a kinetic model instead of the Boltzmann equation for IHS and/or IMM allows us to obtain the explicit forms of the velocity distribution functions. This is likely one of the main advantages of starting from a kinetic model rather than the true Boltzmann kinetic equation.

V.G. acknowledges financial support from Grant No. PID2024-156352NB-I00 funded by MCIU/AEI/10.13039/501100011033/FEDER, UE and from Grant No. GR24022 funded by Junta de Extremadura (Spain) and by European Regional Development Fund (ERDF) “A way of making Europe”.

9 Linear stability analysis of the steady state solutions to USF↩︎

In this Appendix, we want to see whether the steady state solutions 40 –48 to the pressure tensors are (linearly) stable. We consider first the time evolution equation for the pressure tensor of the excess granular gas. The three relevant independent equations for \(P_{2,k\ell}\) are given by \[\label{a1} \partial_t p+\zeta_2 p+\frac{2a}{d}P_{2,xy}=0,\quad \partial_t P_{2,xy}+\nu_\eta P_{2,xy}+a P_{2,yy}=0,\tag{95}\] \[\label{a1bis} \partial_t P_{2,yy}+\nu_\eta P_{2,yy}-\left(\nu_\eta-\zeta_2\right) P_{2,yy}=0,\tag{96}\] where \(\nu_\eta=\nu_{\text{M},22}\nu_\eta^*\) and \(p\simeq p_2=n_2 T_2\). In terms of the dimensionless quantities, \(P_{2,k\ell}^*(t)=P_{2,k\ell}(t)/n_2 T(t)\), \(a^*(t)=a/\nu_{\text{M},22}(t)\) and \[\label{a2} \tau(t)=\int_0^{t} dt'\; \nu_{\text{M},22}(t'),\tag{97}\] Eqs.@eq:a1 and 96 become \[\label{a3} 2\partial_\tau \ln a^*=\zeta_2^*+\frac{2}{d}a^* P_{2,xy}^*,\tag{98}\] \[\label{a4} \partial_\tau P_{2,xy}^*=-a^*P_{2,yy}^*-P_{2,xy}^*\Big(\nu_\eta^*-\zeta_2^*-\frac{2}{d}a^* P_{2,xy}^*\Big),\tag{99}\] \[\label{a5} \partial_\tau P_{2,yy}^*=-P_{2,yy}^*\Big(\nu_\eta^*-\zeta_2^*-\frac{2}{d}a^* P_{2,xy}^*\Big)+\nu_\eta^*-\zeta_2^*,\tag{100}\] The variable \(\tau\) is the dimensionless time measured as the average collision number. A steady solution of Eqs.@eq:a3 –100 is given by Eqs.@eq:3465 and 41 . To carry out a linear stability analysis of these steady solutions, we look for solutions to the set 98 –100 given by \[\label{a6} a^*(\tau)=a_s^*+\delta a^*(\tau), \quad P_{2,xy}^*(\tau)=P_{2xy,s}^*+\delta P_{2,xy}^*(\tau), \quad P_{2,yy}^*(\tau)=P_{2yy,s}^*+\delta P_{2,yy}^*(\tau),\tag{101}\] where the subscript \(s\) means that the quantities are evaluated in the steady state. Substituting the identities 101 into Eqs.@eq:a3 –100 and neglecting nonlinear terms in the perturbations, one gets \[\label{a7} \partial_\tau \left( \begin{array}{c} \delta a_s^*\\ \delta P_{2,xy}^*\\ a_s^*\delta P_{2,yy}^* \end{array} \right)=-\mathsf{L}\cdot \left( \begin{array}{c} \delta a^*\\ \delta P_{2,xy}^*\\ a_s^*\delta P_{2,yy}^* \end{array} \right),\tag{102}\] where \(\mathsf{L}\) is the square matrix \[\mathsf{L}=\left( \begin{array}{ccc} \frac{\zeta_2^*}{2} & -\frac{\zeta_2^*\nu_\eta^{*2}}{2(\nu_\eta^*-\zeta_2^*)} & 0 \\ \left(\frac{\nu_\eta^*-\zeta_2^*}{\nu_\eta^*}\right)^2 & \nu_\eta^*+\zeta_2^*& 1\\ \frac{(\nu_\eta^*-\zeta_2^*)\zeta_2^*}{\nu_\eta^*} & -\zeta_2^*\nu_\eta^* & \nu_\eta^* \end{array} \right). \label{a8}\tag{103}\]

The time evolution of the deviations from the steady solution is governed by the three eigenvalues of the matrix \(\mathsf{L}\). If the real parts of those eigenvalues are positive the steady solution is linearly stable, while is unstable otherwise. The eigenvalues are determined from the solution of the secular equation \[\label{a9} \det \left(L_{k\ell}-\ell \delta_{k\ell}\right)=0.\tag{104}\] The solution to Eq.@eq:a9 leads to a real eigenvalue \(\ell_1\) and a pair of complex conjugate eigenvalues \(\ell_2\) and \(\ell_3\). As for IHS [67], the results show that \(\text{Re}(\lambda_i)>0\) (\(i=1,2,3\)) for any value of \(\alpha_{22}\). This means that the steady USF solution for the excess granular is linearly stable, and the characteristic relaxation time (measured by the number of collisions) is \(\ell_1^{-1}\). As an illustration, figure [appAfig1] shows the dependence of \(\ell_1\) and of the real part of \(\ell_{2,3}\) for a three-dimensional granular gas (\(d=3\)). Note that \(\ell_1\to 0\) in the elastic limit \(\alpha_{22}\to 1\). This is a consequence that for elastic collisions \(a^*\to 0\).

We consider now the (linear) stability of the steady solution to the set of time evolution equations associated with \(\gamma(t)\), \(P_{1,xy}(t)\) and \(P_{1,yy}(t)\). Since we have previously shown that the perturbations \((\delta a^*, \delta P_{2,xy}^*, \delta P_{2,yy}^*)\) tend to zero for sufficiently long times, we assume hence that \(a^*\equiv \text{const.}\), \(P_{2,xy}^*\equiv \text{const.}\) and \(P_{2,yy}^*\equiv \text{const.}\) in the evolution equations of \(\gamma(t)\), \(P_{1,xy}(t)\) and \(P_{1,yy}(t)\). In terms of the variable \(\tau\), the set of equations for \(\gamma\), \(P_{1,xy}\) and \(P_{1,yy}\) is \[\label{a10} \partial_\tau \gamma=-\frac{2}{d}a^* P_{1,xy}^*-\zeta_1^* \gamma,\tag{105}\] \[\label{a11} \partial_\tau P_{1,yy}^*=Y+ X_0 P_{1,yy}^*+X P_{2,yy}^*,\tag{106}\] \[\label{a12} \partial_\tau P_{1,xy}^*+a^* P_{1,yy}^*= X_0 P_{1,xy}^*+X P_{2,yy}^*,\tag{107}\] where the quantities \(Y\), \(X_0\) and \(X\) are defined by Eqs.@eq:3468 and 45 , respectively. As in the case of the excess granular gas, we want to solve the set of Eqs.@eq:a10 –107 by assuming small deviations from the steady state solution. Thus, we write \[\label{a13} \gamma(\tau)=\gamma_s+\delta \gamma (\tau), \quad P_{1,yy}^*(\tau)=P_{1yy,s}^*+\delta P_{1,yy}^*(\tau), \quad P_{1,xy}^*(\tau)=P_{1xy,s}^*+\delta P_{1,xy}^*(\tau).\tag{108}\] In the linear order in the perturbations, \[\label{a14} \zeta_1^*=\zeta_{1s}^*+\overline{\zeta}_1 \delta \gamma, \quad Y=Y_s+\overline{Y}\delta\gamma,\quad X_0=X_{0s}+\overline{X}_0\delta\gamma, \quad X=X_s+\overline{X}\delta\gamma,\tag{109}\] where \(\zeta_{1s}^*\), \(Y_s\), \(X_{0s}\), and \(X_s\) refer to the values of these quantities in the steady state and \[\begin{align} \label{a15} \overline{\zeta}_1= \frac{\sqrt{2}}{d}\left(\frac{\sigma_{12}}{\sigma_2}\right)^{d-1}\mu_{21}(1+\alpha_{12})\Bigg\{ \frac{(1+\theta_s^{-1})^{-1/2}}{2\mu}\left[1-\frac{\mu_{21}}{2}(1+\alpha_{12})(1+\theta_s)\right]+\frac{\mu_{21}^2}{2\mu_{12}}\theta_s^2(1+\theta_s^{-1})^{1/2}(1+\alpha_{12}) \Bigg\}, \end{align}\tag{110}\]

\[\label{a16} \overline{Y}=\frac{3}{2\sqrt{2}}\left(\frac{\sigma_{12}}{\sigma_2}\right)^{d-1}\frac{\mu_{21}^2}{(d+2)}(1+\alpha_{12})^2 (1+\theta_s^{-1})^{1/2},\tag{111}\] \[\label{a17} \overline{X}_0=-\frac{1}{\sqrt{2}d(d+2)}\left(\frac{\sigma_{12}}{\sigma_2}\right)^{d-1}\frac{\mu_{21}^2}{\mu_{12}} (1+\theta_s^{-1})^{-1/2}(1+\alpha_{12})\left[d+2-\mu_{21}(1+\alpha_{12})\right],\tag{112}\] \[\label{a18} \overline{X}=\frac{1}{\sqrt{2}d(d+2)}\left(\frac{\sigma_{12}}{\sigma_2}\right)^{d-1}\mu_{21}^2 (1+\theta_s^{-1})^{-1/2}(1+\alpha_{12})^2.\tag{113}\] Substitution of Eqs.@eq:a13 and 109 into Eqs.@eq:a10 –107 and neglecting nonlinear terms in the perturbations, after some algebra one gets the set of linear differential equations: \[\label{a19} \partial_\tau \left( \begin{array}{c} \delta \gamma\\ \delta P_{1,yy}^*\\ \delta P_{1,xy}^* \end{array} \right)=-\widetilde{\mathsf{L}}\cdot \left( \begin{array}{c} \delta \gamma\\ \delta P_{1,yy}^*\\ \delta P_{1,xy}^* \end{array} \right),\tag{114}\] where \(\widetilde{\mathsf{L}}\) is the square matrix \[\widetilde{\mathsf{L}}=\left( \begin{array}{ccc} \zeta_{1s}^*+\gamma_s \overline{\zeta}_1 & 0& \frac{2}{d}a_s^* \\ -\left(\overline{Y}+\overline{X}_0P_{1,yy,s}^*+\overline{X}P_{yy,s}^*\right) & -X_{0s}& 0\\ -\left(\overline{X}_0P_{1,xy,s}^*+\overline{X}P_{yy,s}^*\right) & a_s^* & -X_{0s} \end{array} \right). \label{a20}\tag{115}\] The eigenvalues of the matrix \(\widetilde{\mathsf{L}}\) are the roots of the secular equation \[\label{a21} \det \left(\widetilde{L}_{k\ell}-\widetilde{\ell} \delta_{k\ell}\right)=0.\tag{116}\] If the real parts of the eigenvalues \(\widetilde{\ell}\) are always positive for any value of the set \((\alpha_{22},\alpha_{12},m_1/m_2, \sigma_1/\sigma_2)\) then the steady USF solution for the pressure tensor \(P_{1,ij}^*\) is linearly stable.

A systematic analysis of the dependence of \(\text{Re}(\; \widetilde{\ell}_i)\) (\(i=1,2,3\)) on the parameter space shows that the real parts of the eigenvalues are always positive. As an illustration, figures [appAfig2] and [appAfig3] show the \(\alpha\)-dependence of \(\text{Re}(\; \widetilde{\ell}_i)\); being quite apparent that their real parts are positive.

10 Behavior of the zeroth-order pressure tensors near the steady state↩︎

In this Appendix we give the expressions of the derivatives of the zeroth-order pressure tensors \(P_{k\ell}^*\) and \(P_{1,k\ell}^*\) with respect to \(a^*\) near the steady state. We consider first IMM where the (dimensionless) elements of the pressure tensor \(P_{k\ell}^*\) obey the equation \[\label{b1} -\left(\frac{2}{d}a^* P_{xy}^{*}+\zeta^{*}\right)\left(1-\frac{1}{2}a^*\frac{\partial}{\partial a^*}\right) P_{k\ell}^{*}+ a_{k\mu}^*P_{\ell\mu}^{*}+a_{\ell\mu}^*P_{k\mu}^{*}= -\left[\nu_\eta^* P_{k\ell}^{*}+\left(\zeta^*-\nu_\eta^*\right)\delta_{k\ell}\right],\tag{117}\] where in the tracer limit \(P_{2,k\ell}^{*}\simeq P_{k\ell}^{*}\) and \(\zeta_2^*\simeq \zeta^*\). From Eq.@eq:b1 , one gets the set of equations \[\label{b2} \frac{\partial P_{yy}^*}{\partial a^*}=\frac{2(\nu_\eta^*-\zeta^*)-2P_{yy}^*\left(\nu_\eta^*-\zeta^*-\frac{2}{d}a^* P_{xy}^*\right)}{a^*\left(\frac{2}{d}a^* P_{xy}^{*}+\zeta^{*}\right)},\tag{118}\] \[\label{b3} \frac{\partial P_{xy}^*}{\partial a^*}=\frac{-2P_{yy}^*a^*-2P_{xy}^*\left(\nu_\eta^*-\zeta^*-\frac{2}{d}a^* P_{xy}^*\right)}{a^*\left(\frac{2}{d}a^* P_{xy}^{*}+\zeta^{*}\right)}.\tag{119}\] As expected, the numerators and denominators of Eqs.@eq:b2 and 119 vanish in the steady state \([(2/d)a^*P_{xy}^*+\zeta^*=0]\). As in the case of IHS [58], the steady-state limit of Eqs.@eq:b2 and 119 can be evaluated by means of l’Hopital’s rule. In this case, one achieves the results \[\label{b4} \Delta_{yy}=4 P_{yy}^*\frac{a^*\Delta_{xy}+P_{xy}^*}{2a^{*2}\Delta_{xy}+d\left(2\nu_\eta^*-\zeta^*\right)},\tag{120}\] where \(\Delta_{yy}\equiv \left(\partial P_{yy}^*/\partial a^*\right)_s\) and \(\Delta_{xy}\equiv \left(\partial P_{xy}^*/\partial a^*\right)_s\) is the real root of the cubic equation \[\begin{align} \label{b5} 2 a^{*4} \Delta_{xy}^3+4da^{*2}\nu_\eta^*\Delta_{xy}^2+\frac{d^2}{2}\left(4\nu_\eta^{*2}+6 \zeta^*\nu_\eta^*-3\zeta^{*2}\right)\Delta_{xy}+d^2\left(\nu_\eta^*-\zeta^*\right) \nu_\eta^{*-2}\left(\zeta^{*2}-5\zeta^*\nu_\eta^*+ 2\nu_\eta^{*2}\right)=0. \end{align}\tag{121}\] In the above equations, it is understood that all the quantities are computed in the steady state.

We consider now the derivatives of the elements \(P_{1,k\ell}^*\) with respect to \(a^*\). They verify the time dependent equation \[\label{b6} -\left(\frac{2}{d}a^* P_{xy}^{*}+\zeta^{*}\right)\left(1-\frac{1}{2}a^*\frac{\partial}{\partial a^*}\right) P_{1,k\ell}^{*}+ a_{k\mu}^*P_{1,\ell\mu}^{*}+a_{\ell\mu}^*P_{1,k\mu}^{*}= Y \delta_{k\ell}+X_0 P_{1,k\ell}^*+X P_{k\ell}^*,\tag{122}\] where we recall that the quantities \(Y\), \(X_0\) and \(X\) for IMM are defined by Eqs.@eq:3468 and 45 , respectively. The derivatives of \(\partial_{a^*}P_{1,yy}^{*}\), \(\partial_{a^*}P_{1,xy}^{*}\), and \(\partial_{a^*}\gamma\) can be easily obtained from the results derived for IHS in Ref.[58] by replacing the expressions of the quantities \(Y\), \(X_0\), and \(X\) of IHS by their corresponding counterparts for IMM given by Eqs.@eq:3468 and 45 , respectively. Thus, the expressions of the derivatives \(\partial_{a^*}P_{1,xy}^{*}\) and \(\partial_{a^*}P_{1,xy}^{*}\) are given by Eqs.(B.16) and (B.17), respectively, of Ref.[58] while the derivative \(\partial_{a^*}\gamma\) is ³ \[\label{b7} \left(\frac{\partial \gamma}{\partial a^*}\right)_s=\frac{\Lambda_1}{\Lambda_2},\tag{123}\] where \[\begin{align} \label{b8} \Lambda_1=d\left(\frac{1}{2}a^*\chi-X_0\right)\Bigg[\left(\frac{1}{2}a^*\chi-X_0\right) \left(\gamma\chi-\frac{2}{d}P_{1,xy}^*\right)-\frac{2}{d}a^*\Big(\chi P_{1,xy}^*-P_{1,yy}^*+X \Delta_{xy} \Big)\nonumber\\ +2a^{*2}\Big(\chi P_{1,yy}^*+X \Delta_{yy}\Big)\Bigg], \end{align}\tag{124}\] \[\begin{align} \label{b9} \Lambda_2=d\left(\frac{1}{2}a^*\chi-X_0\right)\Bigg[\left(\frac{1}{2}a^*\chi-X_0\right)\left(\zeta_1^*+\frac{1}{2}a^*\chi+\gamma \zeta_1^*\right)+\frac{2}{d}a^*\left(X_0'P_{1,xy}^*+X'P_{xy}^*\right)\Bigg]\nonumber\\ -2a^{*2}\left(Y'+X_0'P_{1,yy}^*+X'P_{yy}^*\right). \end{align}\tag{125}\] In Eqs.@eq:b8 and 125 , \(\chi=(2/d)(P_{xy}^*+a^* \Delta_{yy})\), \(Y'=\partial_\gamma Y\), \(X'=\partial_\gamma X\), and \(X_0'=\partial_\gamma X_0\). As in the case of the excess granular gas, all the quantities appearing in Eqs.@eq:b7 –125 are evaluated in the steady state.

In the case of the VGS kinetic model, the expressions of the derivatives \(\Delta_{yy}\) and \(\Delta_{xy}\) are given by Eqs. 120 and 121 , respectively, except that \(\nu_\eta^*=(1+\alpha_{22})/2+\epsilon_{22}^*\) and \(\zeta^*\) must be replaced by \(\epsilon_{22}^*\). With respect to the derivatives associated with the tracer particles, their forms are identical to those obtained for IMM except that \(X=0\), and the quantities \(Y\) and \(X_0\) are given by \[\label{b10} Y=\frac{1+\alpha_{12}}{2}\nu_{12}^* \Big[\gamma+2\mu_{12}\mu_{21}(1-\gamma) \Big],\tag{126}\] \[\label{b11} X_0=-\frac{1+\alpha_{12}}{2}\nu_{12}^*-\epsilon_{12}^*.\tag{127}\]

References↩︎