Kinetic Theory of Polydisperse Granular Mixtures: Influence of the Partial Temperatures on Transport Properties—A Review

It is well-recognized that granular media under rapid flow conditions can be modeled as a gas of hard spheres with inelastic collisions. At moderate densities, a fundamental basis for the determination of the granular hydrodynamics is provided by the Enskog kinetic equation conveniently adapted to account for inelastic collisions. A surprising result (compared to its molecular gas counterpart) for granular mixtures is the failure of the energy equipartition, even in homogeneous states. This means that the partial temperatures Ti (measuring the mean kinetic energy of each species) are different to the (total) granular temperature T. The goal of this paper is to provide an overview on the effect of different partial temperatures on the transport properties of the mixture. Our analysis addresses first the impact of energy nonequipartition on transport which is only due to the inelastic character of collisions. This effect (which is absent for elastic collisions) is shown to be significant in important problems in granular mixtures such as thermal diffusion segregation. Then, an independent source of energy nonequipartition due to the existence of a divergence of the flow velocity is studied. This effect (which was already analyzed in several pioneering works on dense hard-sphere molecular mixtures) affects to the bulk viscosity coefficient. Analytical (approximate) results are compared against Monte Carlo and molecular dynamics simulations, showing the reliability of kinetic theory for describing granular flows.


Introduction
It is well-known that when granular matter is subjected to a violent and sustained excitation, the motion of grains resembles to the random motion of atoms or molecules in an ordinary or molecular gas. In this situation (referred usually to as rapid flow conditions), the energy injected to the system compensates for the energy dissipated by collisions and the effects of gravity. A system of activated collisional grains is referred to as a granular gas; its study is the main objective of the present review.
Granular matter in nature is usually immersed in a fluid like water or air, so that a granular flow is a multiphase process. However, under some conditions (for instance, when the stress due to grains is larger than that exerted by the interstitial fluid), the effect of the fluid phase on grains can be neglected. Here, we will address our attention to the study of the so-called dry granular gases where the impact of the fluid phase on the dynamics of solid particles is not accounted for.
Since the grains which make up a granular material are of macroscopic size (their diameter is micrometers or larger), all the collisions among granular particles are inelastic. This is one of the main differences to molecular gases. Due to this fact, the conventional methods of equilibrium statistical mechanics and thermodynamics fail. However, kinetic theory (which essentially addresses the dynamics of grains) is still an appropriate tool since it applies to elastic or inelastic collisions [1,2]. As we are mainly interested in assessing the effect of inelasticity of collisions on the dynamical properties of the granular particles, it is quite usual to consider a relatively simple (idealized) model which isolates the collisional dissipation effect from other relevant properties of granular matter. The most popular model for granular gases is a system of identical smooth hard spheres with a constant (positive) coefficient of normal restitution α 1. This quantity measures the ratio between the magnitude of the normal component of the relative velocity (oriented along the line separating the centers of the two spheres at contact) before and after a collision. The case α = 1 corresponds to perfectly elastic collisions while when α < 1 part of the kinetic energy of the relative motion is lost.
Within the context of the inelastic hard sphere model, the Boltzmann and Enskog kinetic equations have been conveniently extended to account for the dissipative character of collisions [1][2][3][4][5][6][7][8][9][10]. While the Boltzmann equation applies to low-density gases, the Enskog equation holds for moderately dense gases. These kinetic equations have been employed in the last few years as the starting point to derive the corresponding granular hydrodynamic equations. In particular, in the case of monocomponent granular gases and assuming the existence of a normal (or hydrodynamic) solution for sufficiently long space and time scales, the Chapman-Enskog [11] and Grad's moment [12] methods have been applied to solve the Boltzmann and Enskog kinetic equations to the Navier-Stokes order and obtain explicit expressions for the transport coefficients [13][14][15][16][17][18][19][20].
On the other hand, since a real granular system is usually characterized by some degree of polydispersity in density and size, flows of granular mixtures are prevalent in both nature and industry. For instance, natural systems that are highly polidisperse and propagate as rapid granular flows are pyroplastic density currents [21], landslides and debris flows [22], and rock avalanches [23]. Examples of industrial systems include mixing of pharmaceutical powders and poultry feedstock.
Needless to say, in the context of kinetic theory, the determination of the Navier-Stokes transport coefficients of a granular mixture is more intricate than that of a monocomponent granular gas since not only is the number of transport coefficients larger than for a single gas but also they depend on many parameters (masses and diameters, concentrations, and coefficients of restitution). Thus, due to this type of technical difficulties, many of the early attempts [24][25][26][27] to obtain the transport coefficients of a granular mixture were carried out by assuming equipartition of energy: the partial temperatures T i of each species are equal to the (total) granular temperature T. A consequence of this assumption is that the Chapman-Enskog expansion was performed around Maxwellian distributions at the same temperature T for each species. The use of this Maxwellian distribution as the reference state in the Chapman-Enskog method can be only considered as reliable for nearly elastic spheres where the energy equipartition still holds. Moreover, within this level of approximation, the expressions of the transport coefficients are the same as those obtained for molecular (elastic) mixtures [11,28]; the inelasticity in collisions is only taken into account by the presence of a sink term in the energy balance equation.
However, many different works based on kinetic theory [29,30], computer simulations [31][32][33][34][35][36][37][38][39][40][41][42][43][44] and real experiments [45,46] have clearly shown the failure of the energy equipartition in granular mixtures. This failure occurs even in homogeneous situations (in the so-called homogeneous cooling state) and is a consequence of both the inelasticity in collisions and the mechanical differences of the particles (e.g., masses, diameters). In fact, nonequipartition disappears when collisions between the different species of the mixture are elastic or when they are mechanically equivalent. Although the possibility of energy nonequipartition in granular mixtures was already noted by Jenkins and Mancini [47], to the best of our knowledge the impact of nonequipartition on transport properties in granular mixtures was computed for the first time by Huilin et al. [48,49]. However, these authors do not attempt to solve the kinetic equation and they assume local Maxwellian distribution functions for each species even in inhomogeneous states. Although this procedure can be employed to get the collisional transfer contributions to the fluxes, it predicts vanishing Navier-Stokes transport coefficients for dilute granular mixtures which is of course a wrong result. A more rigorous way of incorporating energy nonequipartition in the Chapman-Enskog solution has been published in the past few years [50][51][52][53][54]. The results have clearly shown that in general the effect of temperatures differences on the Navier-Stokes transport coefficients are important, specially for disparate masses or sizes and/or strong inelasticity.
As the Chapman-Enskog procedure states [11], since the partial temperatures are kinetic quantities, they must be also expanded in terms of gradients of the hydrodynamic fields. The partial temperatures are scalars so that, their first-order contributions T (1) i must be proportional to the divergence of the flow velocity U. Thus, a different way of inducing a breakdown of the energy equipartition in granular mixtures is by the presence of the gradient ∇ · U. This effect is not generic of granular mixtures since it was already found in the pioneering works of dense hard-sphere mixtures with elastic collisions [55][56][57].
The non-vanishing divergence of the mean flow velocity ∇ · U causes that T (1) i is involved in the evaluation of the bulk viscosity (proportionality coefficient between the collisional part of the pressure tensor and ∇ · U) as well as in the first-order contribution to the cooling rate ζ (which accounts for the rate of kinetic energy dissipated by collisions).
The aim of this paper is to offer a short review on the influence of the energy nonequipartition on transport properties in granular mixtures. Since we will consider moderate densities, the one-particle velocity distribution functions of each species will obey the set of coupled Enskog kinetic equations. The review is structured as follows. The set of Enskog coupled kinetic equations for a multicomponent granular mixture and its associated macroscopic balance equations are introduced in Section 2. In particular, explicit forms for the collisional transfer contributions to the fluxes are given in terms of the one-particle velocity distribution function f i of each species. Section 3 deals with the solution to the Enskog equation in the homogeneous cooling state; a homogeneous state where the granular temperature decreases in time due to inelastic cooling. As for monocomponent granular gases, a scaling solution is proposed in which the time dependence of the distributions f i occurs entirely through the temperature of the mixture T. The temperature ratios T i /T = 1 (energy nonequipartition). In Section 3, the (approximate) theoretical results are compared against the results obtained from both the Direct Simulation Monte Carlo (DSMC) method and molecular dynamics (MD) simulations for conditions of practical interest. Comparison shows in general a good agreement between theory and simulations. The forms of the Navier-Stokes transport coefficients of the mixture in terms of the first-order distributions derived from the application of the Chapman-Enskog method around the local version of the homogeneous distributions obtained in Section 2 are displayed in Section 4. Section 5 addresses one of the main targets of the present paper: the study of the influence of the temperature ratios T (0) i /T on the transport coefficients. To show more clearly the impact of nonequipartition on transport, we focus here on our attention to the diffusion transport coefficients of a dilute granular binary mixture. As expected, we find that the effect of energy nonequipartition on transport is in general quite significant. This means that, beyond nearly elastic systems, any reliable theory devoted to granular mixtures must include this nonequipartition effect. The influence of the first-order contributions T (1) i to the partial temperatures on the bulk viscosity η b and the cooling rate ζ is widely analyzed in Section 6.
The contributions to η b coming from the coefficients T (1) i were implicitly neglected in several previous works [53,54,58] on dense granular mixtures. Our present results indicate that the impact of T (1) i on η b cannot be neglected for disparate masses and/or strong inelasticity. The paper is ended in Section 8 with a brief discussion of the results reported here.
Before ending this section, we want to remark that the present account is based on the authors' taste and perspective. In this sense, no attempt is made to include the extensive related work of many others in this field. The references given are selective and apologies are offered at the outset to the many other important contributions not recognized explicitly.

Enskog Kinetic Equation for Inelastic Hard Spheres
We consider a granular mixture of inelastic hard disks (d = 2) or spheres (d = 3) of masses m i and diameters σ i (i = 1, 2, . . . , s). The subscript i labels one of the s mechanically different species or components and d is the dimension of the system. For the sake of simplicity, we assume that the spheres are completely smooth; this means that inelasticity of collisions between particles of species i and j is only characterized by the constant (positive) coefficients of restitution α ij 1. The coefficient α ij measures the ratio between the magnitude of the normal component (along the line separating the centers of the two spheres at contact) of the relative velocity after and before the collision i-j. Figure 1 shows a schematic diagram of a ternary (s = 3) hard-sphere granular mixture. For moderate densities, the one-particle velocity distribution function f i (r, v, t) of species i verifies the set of s-coupled nonlinear integro-differential Enskog equations. In the absence of any external force, the Enskog kinetic equations are given by [10] where the Enskog collision operator In Equation (2), σ ij = σ ij σ, σ ij = (σ i + σ j )/2, σ is a unit vector directed along the line of centers from the sphere of species i to that of species j at contact, Θ is the Heaviside step function, and g 12 = v 1 − v 2 is the relative velocity of the colliding pair. Moreover, χ ij (r 1 , r 1 + σ ij ) is the equilibrium pair correlation function of two hard spheres, one of species i and the other of species j at contact, i.e., when the distance between their centers is σ ij .
As in the case of elastic hard spheres, the interactions between inelastic hard spheres are modeled by instantaneous collisions where momentum is transferred along the line joining the centers of the two colliding spheres. This reflection is inelastic, as illustrated in Figure 2 for colliding spheres of diameters σ 1 and σ 2 . The relationship between the pre-collisional velocities (v 1 , v 2 ) and the post-collisional velocities where µ ij = m i /(m i + m j ). Equations (3) give the so-called inverse or restituting collisions. Inversion of these collision rules provides the form of the so-called direct collisions, namely, collisions where the pre-collisional velocities (v 1 , v 2 ) lead to the post-collisional velocities (v 1 , v 2 ) [8]: From Equations (3) and (4), one gets the relations where For inelastic collisions, it is quite apparent from Equation (5) that the magnitude of the normal component of the pre-collisional relative velocity is larger than its post-collisional counterpart. In addition, comparison between Equations (3) and (4) shows that, except for molecular mixtures (elastic collisions), the direct and inverse collisions are not equivalent. This is essentially due to the lack of time reversal symmetry for inelastic collisions.
The change in kinetic energy of the colliding pair in a binary collision can be easily obtained from Equation (4): where m ij = m i m j /(m i + m j ) is the reduced mass. When α ij = 1 (elastic collisions), Equation (6) leads to ∆E ij = 0, as expected for molecular mixtures. When α ij < 1, ∆E ij < 0 so that, part of the kinetic energy is lost in a binary collision between a particle of species i and a particle of species j.

Macroscopic Balance Equations
The knowledge of the velocity distribution functions f i allows us to obtain the hydrodynamic fields of the multicomponent mixture. The quantities of interest in a macroscopic description of the granular mixture are the local number density n i of species i, the local mean flow velocity of the mixture U, and the granular temperature T. In terms of the distributions f i , they are defined, respectively, as In Equations (7)-(9), ρ = ∑ i m i n i is the total mass density, ρ i = m i n i is the mass density of the species i, n = ∑ i n i is the total number density, and V = v − U is the peculiar velocity. For the subsequent discussion, at a kinetic level, it is convenient to introduce the partial kinetic temperatures T i for each species. The temperature T i provides a measure of the mean kinetic energy of the species i. The partial temperatures are defined as From Equations (9) and (10), the granular temperature T of the mixture can be also written in terms of the partial temperatures T i as where x i = n i /n is the concentration or mole fraction of species i. Thus, due to the constraint (11), there are s − 1 independent partial temperatures in a mixture constituted by s components.
As occurs for molecular mixtures [28], the fields n i , U, and T are expected to be the slow variables that dominate the dynamics of the mixture for sufficiently long times through the set of hydrodynamic equations. For elastic collisions, the above fields are the densities of global conserved quantities and so, they persist at long times (in comparison with the mean free time) where the complex microscopic dynamics becomes negligible [59,60]. In the case of granular fluids, the energy is not conserved in collisions and the rate of energy dissipated by collisions is characterized (as we will see below) by a cooling rate. However, as confirmed by MD simulations (see, for instance, Ref. [33]), the cooling rate may be slow compared to the transient dynamics so that, the kinetic energy (or granular temperature) can be still considered as a slow variable.
The balance equations for n i , U, and T can be obtained by multiplying the set of Enskog equations (1) by 1, m i v, and m i 2 V 2 and summing over all the species in the momentum and energy equations. The result is (14) where D t = ∂ t + U · ∇ is the material derivative. In the above equations, is the mass flux for species i relative to the local flow U and the kinetic contributions P k and q k to the pressure tensor and heat flux are given, respectively, by A consequence of the definition (15) of the fluxes j i is that only s − 1 mass fluxes are independent since they have the constraint Needless to say, to end the derivation of the balance hydrodynamic equations one has to compute the right-hand side of Equations (13) and (14). These terms can be obtained by employing an important property of the integrals involving the Enskog collision operator J ij [r, v| f i , f j ] [10,53]: where ψ i (v 1 ) is an arbitrary function of v 1 , v 1 is defined by Equation (4) and The first term on the right hand side of Equation (19) represents a collisional effect due to scattering with a change in velocities. This term vanishes for elastic collisions. The second term on the right hand side of Equation (19) provides a pure collisional effect due to the spatial difference of the colliding pair. This term vanishes for low-density mixtures.
In the case ψ = m i v, the first term in the integrand (19) disappears since the momentum is conserved in all pair collisions, i.e., m i v 1 + m j v 2 = m i v 1 + m j v 2 . The second term in the integrand yields the result where the collision transfer contribution to the pressure tensor P c is [53] In the case ψ = 1 2 m i V 2 , the first term on the right hand side of Equation (19) does not vanish since the kinetic energy is not conserved in collisions. As before, the second term in the integrand gives the collisional transfer contribution to the heat flux q c . The result is [53] where and the (total) cooling rate ζ due to inelastic collisions among all species is given by The balance hydrodynamic equations for the densities of momentum and energy can be finally written when Equations (21)- (25) are substituted into the right hand sides of Equations (13) and (14). These balance equations can be written as where the pressure tensor P(r, t) and the heat flux q(r, t) have both kinetic and collisional transfer contributions, i.e., Equations (12), (26) and (27) are the balance equations for the hydrodynamic fields n i , U, and T, respectively, of a polydisperse granular mixture at moderate densities. This set of equations do not constitute a closed set of equations unless one expresses the fluxes and the cooling rate in terms of the above hydrodynamic fields and their spatial gradients. For small gradients, the corresponding constitutive equations for the fluxes and the cooling rate can be obtained by solving the set of Enskog kinetic equation (1) with the extension of the conventional Chapman-Enskog method [11] to dissipative dynamics.
Before closing this section, it is instructive to consider the case of dilute polydisperse granular mixtures. The corresponding balance equations can be obtained from Equations (12), (26) and (27) by taking χ ij → 1 and neglecting the different centers [(r, r ± σ ij )] of the colliding pair since the effective diameter σ ij is much smaller than that of the mean free path of this collision. This implies that the collision transfer contributions to the fluxes are much smaller than their corresponding kinetic counterparts (P → P k and q → q k ) and the cooling rate ζ is simply given by

Homogeneous Cooling State: Partial Temperatures
Let us consider a spatially homogeneous state of an isolated polydisperse granular mixture. In contrast to molecular mixtures, there is no longer an evolution toward the local Maxwellian distributions since those distributions are not a solution to the set of homogeneous (inelastic) Enskog equations. Instead, as we will show below, there is an special solution which is achieved after a few collision times by considering homogeneous initial conditions: the so-called homogeneous cooling state (HCS).
For spatially homogeneous isotropic states, the set of Enskog equations for the distri- where the Boltzmann collision operator J B ij is Upon writing Equations (30) and (31) we have taken into account that the dependence of the distributions f i on velocity v is only through its magnitude v. In Equation (30), note that χ ij refers to the pair correlation function for particles of species i and j when they are separated a distance σ ij . For homogeneous states, the balance Equations (12) and (26) trivially hold. On the other hand, the balance equation of the granular temperature (27) yields where the cooling rate ζ is defined in Equation (25) by making the replacement f ij (r 1 , v 1 , On the other hand, for homogeneous states, the integration in σ can be easily performed and ζ can be more explicitly written as Moreover, for symmetry reasons, the mass and heat fluxes vanish and the pressure tensor P k = pδ k , where the hydrostatic pressure p is [53] where denotes the partial temperature of species i in the homogeneous state (absence of spatial gradients).
To analyze the rate of change of the partial temperatures T (0) i , it is convenient to introduce the "partial cooling rates" ζ i . The definition of these quantities can be obtained by multiplying both sides of the Enskog equation (30) by m i 2 v 2 and integrating over velocity. The result is ∂T where From Equations (11), (32), and (35), one can express the total cooling rate ζ in terms of the partial cooling rates ζ i as The time evolution of the temperature ratios γ i (t) can be easily derived from Equations (32) and (35) as The term ζ ii gives the contribution to the partial cooling rate ζ i coming from the rate of energy loss from collisions between particles of the same species i. This term vanishes for elastic collisions but is different from zero when α ii < 1. The remaining contributions ζ ij (i = j) to ζ represent the transfer of energy between a particle of species i and particles of species j. In general, the term ζ ij = 0 (i = j) for both elastic and inelastic collisions. However, in the special state where the distribution functions f i are Maxwellian distributions at the same temperature (T (0) i = T for any species i), then ζ ij = 0 (i = j) for elastic collisions. This is a consequence of the detailed balance for which the energy transfer between different species is balanced by the energy conservation for this state [29].
The corresponding detailed balance state for inelastic collisions is the HCS. In this state, since the partial ζ i and total ζ cooling rates never vanish, the partial T (0) i and total T temperatures are always time dependent. As for monocomponent granular gases [8,61], whatever the initial uniform state considered is, we expect that the Enskog equation (30) tends toward the HCS solution where all the time dependence of the distributions f i (v; t) only occurs through the (total) temperature T(t). In this sense, the HCS solution qualifies as a normal or hydrodynamic solution since the granular temperature T(t) is in fact the relevant temperature at a hydrodynamic level. Thus, it follows from dimensional analysis that the distributions f i (v; t) have the form [29] where v th (t) = 2T(t)/m is a thermal velocity defined in terms of the global temperature T(t) of the mixture, m = (∑ i m i )/s, and ϕ i is a reduced distribution function whose dependence on the (global) granular temperature T(t) is through the dimensionless velocity v/v th (t).
Since the time dependence of the HCS solution (39) for f i only occurs through the (global) temperature T(t), then the temperature ratios γ i must be independent of time.
This means that all partial temperatures T and so, the temperatures of the species do not provide any new dynamical degree of freedom at the hydrodynamic stage. However, they still characterize the shape of the velocity distribution functions of each species and affect the quantitative averages (mass, momentum, and heat fluxes) calculated with these distributions.
As the temperature ratios do not depend on time, one possibility would be that However, the ratios T (0) i /T (i = 1, . . . , s) must be determined by solving the set of Enskog equations (30). As we will show latter, the above ratios are in general different from 1 and exhibit a complex dependence on the parameter space of the mixture.
Since γ i ≡ const, according to Equation (38), the partial cooling rates ζ i must be equal in the HCS: In addition, the right hand side of Equation (1) can be more explicitly written when one takes into account Equation (39): where use has been made of the identity Therefore, in dimensionless form, the Enskog equation (30) reads where is an effective collision frequency of the mixture and σ = (∑ i σ i )/s. The use of ζ * i instead of ζ * = ζ/ν on the left hand side of Equation (43) is allowed by the equality (40); this choice is more convenient since the first few velocity moments of Equation (43) are verified without any specification of the distributions ϕ i .
We are in front of a well-possed mathematical problem since we have to solve the set of s Enskog equation (30) for the velocity distribution functions f i (v; t) of the form (39) and subject to the s − 1 constraints (40). These 2s − 1 equations must be solved to determine the s distributions f i and the s − 1 temperature ratios γ i . As in the case of monocomponent granular gases [61], approximate expressions for the above quantities are obtained by considering the first few terms of the expansion of the distributions f i in a series of Sonine (or Laguerre) polynomials.
Before obtaining approximate expressions for the temperature ratios, it is important to remark that the failure of energy equipartition in granular fluids has been confirmed in computer simulation works [31][32][33][34][35][36][37][38][39][40][41][42][43][44] and even observed in real experiments of agitated mixtures [45,46]. All the studies conclude that the departure from energy equipartition depends on the mechanical differences between the particles of the mixture as well as the coefficients of restitution.

Approximate Solution
As usual, we expand the distributions ϕ i (c) in a complete set of orthogonal polynomials with a Gaussian measure. In practice, generalized Laguerre or Sonine polynomials   The leading Sonine approximation to the distribution ϕ i is [10] where Note that the parameters of the Gaussian prefactor in (44) are chosen such that ϕ i is normalized to 1 and its second moment ( d 2 θ −1 i ) is consistent with the exact moment (10). An advantage of this choice is that the leading Sonine polynomial is of degree 4. The coeffi- When ϕ i = ϕ i,M , Equation (46)   3 has been carried out for monocomponent granular gases [62]. The results show that the influence of the coefficient a Here, for the sake of simplicity, we will neglect the coefficients a (i)

.
The use of the leading Sonine approximation (44) to ϕ i permits to estimate the partial cooling rates through their definition (36). This involves the evaluation of some intricate collision integrals where nonlinear terms in a (i) 2 are usually neglected. Such approximation is based on the fact that the coefficients a (i) 2 are expected to be very small. In this case, ζ * i can be written as where and the expressions of the quantities ζ (1) ij are very large to be displayed here. They can be found for instance in Refs. [63,64]. The temperature ratios γ i can be already determined in the Maxwellian approximation (i.e., when a (i) 2 = 0) by using Equation (48) in the equality of the cooling rates (40). As will see later, the Maxwellian approximation to γ i leads to a quite accurate predictions.
On the other hand, beyond the Maxwellian approximation, it still remains to estimate the second Sonine coefficients (or kurtosis) a Equation (49) is still exact. However, as in the case of the evaluation of ζ * i , the computation of the collision integrals defining Λ i requires the use of the leading Sonine approximation (44) to achieve explicit results. Neglecting nonlinear terms in a (i) 2 , Λ i can be written as The forms of Λ (0) i and Λ (1) ij can be found in Refs. [63,64]. When Equations (47) and (50) are substituted into Equation (49) and only linear terms in a On the other hand, as noted in several papers on monocomponent granular gases [65][66][67], there is some ambiguity in considering the identity (49) and expands the right hand side as one gets the following system of linear algebraic equations: The solutions to the set of Equations (51) and (54) give the second Sonine coefficients a (i) 2 as functions of the temperature ratios γ i and the parameters of the mixture (masses and diameters, concentrations, coefficients of restitution, and volume fractions). The accuracy of these solutions will be assessed in Section 4 against Monte Carlo simulations in the case of binary mixtures (s = 2). When the expressions of the second Sonine coefficients a (i) 2 are substituted into the s − 1 conditions (40) one achieves the s − 1 temperature ratios γ i . The knowledge of a (i) 2 and γ i in terms of the parameter space of the system allows us to obtain the scaled distributions ϕ i (c) in the leading Sonine approximation (44). This approximate distribution is expected to describe fairly well the behavior of the true distribution in the region of thermal velocities (v ∼ v th , say). In the high velocity region (velocities much larger than that of the thermal one), the distributions ϕ i have an overpopulation [ϕ i (c) ∼ e −ac ] with respect to the Maxwell-Boltzmann tail e −c 2 [31,61,[68][69][70][71][72]. This exponential decay of the tails of the distribution function has been confirmed by computer simulations [31,73,74] and more recently, by means of a microgravity experiment [75].
To obtain the explicit dependence of the temperature ratios and the Sonine coefficients on the system parameters, the pair correlation functions χ ij must be given. Although some attempts have been made [76] for monocomponent granular fluids, we are not aware of any analytical expression of χ ij for granular mixtures. For this reason, we consider here the approximated expression of χ ij proposed for molecular mixtures. Thus, in the case of hard spheres (d = 3), a good approximation for the pair correlation function is [77][78][79] where φ = ∑ i n i πσ 3 i /6 is the solid volume fraction for spheres and M = ∑ i x i σ i .

Some Special Limits
Before illustrating the dependence of the temperature ratios and the second Sonine coefficients on the parameter space for binary (s = 2) and ternary (s = 3) mixtures, it is interesting to consider some simple limiting cases. For mechanically equivalent particles (m i = m, σ i = σ, and α ij = α), the solution to the conditions (40) if we solve Equation (51) or if we solve Equation (54). Equations (56) and (57) agree with the expressions obtained for a 2 for monocomponent granular gases [61,67], as expected.
Another interesting limit corresponds to the tracer limit, namely, a binary mixture where the concentration of one of the species (for example, species 1) is negligible (x 1 → 0). In this limit case, when the collisions between the particles of the excess gas 2 are elastic (α 22 = 1) then the solution to Equations (51) or (54) lead to a The expression (58) for the temperature ratio agrees with the one derived by Martin and Piasecki [30] who found that the Maxwellian distribution with the tracer temperature T (0) 1 defined by Equation (58) is an exact solution to the Boltzmann equation in the above conditions (α 22 = 1 and x 1 → 0).
We assume now that the tracer particles of the binary mixture are much heavier than particles of the excess gas (Brownian limit, i.e., m 2 /m 1 → 0). In this limit case, assuming that the temperature ratio 2 is finite, then the partial cooling rate ζ 1 can be written as where The temperature ratio is determined from the solution to the condition ζ 2 = ζ 1 . The expresion of T (0) where Here, depending on the approximation employed, a 2 is given by Equations (56) or (57). When α 22 = 1, ζ 2 = 0, and Equation (61) yields The expression (63) is consistent with the Brownian limit (m 2 /m 1 → 0) of Equation (58). The expressions (59) and (61) agree with the results obtained by Brey et al. [80]. It is important to remark that a "nonequilibrium" phase transition [81,82] has been found in the Brownian limit which corresponds to a extreme violation of energy equipartition. In other words, there is a region in the parameter space of the system where the temperature ratio goes to infinity and the mean square velocities of the excess gas and the tracer In this region, the Boltzmann-Lorentz collision operator can be well approximated by the Fokker-Planck operator [8,83].

Comparison between Theory and Computer Simulations
In the previous section, we have derived expressions for the temperature ratios γ i and the second Sonine coefficients a (i) 2 of an s-component granular mixture. These (approximate) expressions have been obtained (i) by considering the leading Sonine approximation (44) to the distribution functions and (ii) by retaining only linear terms in a (i) 2 in the algebraic equations defining the above coefficients. To asses the degree of accuracy of these theoretical results, in this section we will compare these predictions with those obtained by numerically solving the Enskog equation by means of the well-known DSMC method [84]. Although this computational method was originally devised for molecular (elastic) fluids, its extension to granular (inelastic) fluids is relatively simple. The simulations allow us to compute the velocity distribution functions over a quite wide range of velocities and obtain precise values of the temperature ratios and the fourth-degree velocity moments in the HCS.

DSMC
In this subsection we provide some details on the application of the DSMC method to a mixture of inelastic hard spheres. More specific details can be found for instance in Ref. [31]. The DSMC algorithm is composed in its basic form of a collision step that handles all particles collisions and a free drift step between particles collisions. As we are interested in solving the set of homogeneous Enskog equations, we take only care of the collisional stage. Thereby, we can consider a single cell wherein the positions of the particles need to be neither computed nor stored.
The velocity distribution function of each species i is represented by the velocities {v k } of N i simulated particles: where δ(x) is the Dirac delta distribution. The system is initialized by drawing the velocities of the particles from Maxwellian velocity distribution functions with temperatures T i,0 . Since the system is dilute enough, only binary collisions are considered. Collisions between particles of species i and j are simulated by choosing a sample of 1 2 max ∆t pairs at random with equiprobability. Here, ∆t is a time step, which is much smaller than mean free time, and ω (ij) max is an upper bound estimate of the probability that a particle of species i collides with a particle of species j per unit of time (typically ω where v th,0 = 2T(0)/m and T(0) is the initial granular temperature). For each pair of particles with velocities {v k , v } (being v k the velocity of a particle of the species i and v of the species j) a given direction σ k is chosen at random with equiprobability. Then, the collision between particles i and j is accepted with a probability equal to Θ(g k · σ k )ω If the collision is accepted, postcollisional velocities of each particle are assigned following the scattering rules (4). For the cases in which ω k . The former procedure is performed for i = 1, 2 and j = 1, 2 in binary and i = 1, 2, 3 and j = 1, 2, 3 in ternary mixtures.

Binary Mixtures
For illustrative purposes, we consider first a binary mixture (s = 2). The parameter space of this system is constituted by the coefficients of restitution (α 11 , α 22 , and α 12 ), the mass (m 1 /m 2 ) and diameter (σ 1 /σ 2 ) ratios, the concentration [x 1 = n 1 /(n 1 + n 2 )], and the solid volume fraction (φ). For the sake of simplicity, henceforth we will consider the case of common coefficients of restitution (α ≡ α ij ) and a three-dimensional system (d = 3). As discussed before, after an initial transient period, one expects that the scaled distribution functions ϕ i (c) reach stationary values independent of the initial preparation of the mixture. This hydrodynamic regime is identified as the HCS. In this regime, the temperature ratio 2 measure the deviations of the scaled distributions ϕ 1 and ϕ 2 , respectively, from their corresponding Maxwellian forms. The panels of Figure 4 show the α-dependence of the above Sonine coefficients for σ 1 /σ 2 = 1, x 1 = 1 2 , φ = 0, and three values of the mass ratio. As for monocomponent granular gases [61], we observe that the coefficients a (i) 2 exhibit a non-monotonic dependence on the coefficient of restitution since they decrease first as inelasticity increases until reaching a minimum value and then increase with decreasing α. We also find that the magnitude of a (i) 2 is in general very small for not quite strong inelasticity (for instance, α 0.5); this supports the assumption of a low-order truncation in the Sonine polynomial expansion of the distributions ϕ i . With respect to the comparison with computer simulations, it is quite apparent that both theoretical estimates for a (i) 2 display an excellent agreement with simulations for values of α 0.6. However, for large inelasticity (α 0.6), the best global agreement with simulations is provided by the approach (54), as Figure 4 clearly shows for values of α ≈ 0.1 (extreme dissipation). Regarding the dependence on the mass ratio, we find that the second Sonine coefficient of the heavier species is larger than that of the lighter species; this means that the departure of f 1 from its Maxwellian form accentuates when increasing the mass ratio m 1 /m 2 .
One of the most characteristic features of granular mixtures, as compared with molecular mixtures, is that the partial temperatures are different in homogenous states. The breakdown of energy equipartition is clearly illustrated in Figure 5 where T  Figure 5 highlights that, at a given value of α, the departure of the temperature ratio from unity increases with increasing the differences in the mass ratio. In general, the temperature of the heavier species is larger than that of the lighter species. Comparison with Monte Carlo simulations shows an excellent agreement in the complete range of values of α. In addition, although not shown here, the theoretical results obtained in the so-called Maxwellian approximation to ϕ i (i.e, when one takes a 2 are practically indistinguishable from those derived by considering the second Sonine coefficients. This means that the impact of these coefficients on the partial cooling rates is negligible and so, the Maxwellian approximation to ϕ i is sufficiently accurate to estimate the temperature ratio. After having studied the effect of the mass ratio on the temperature ratio, we now turn to further assessing the impact of inelastic collisions on this quantity. To analyze this influence, we consider the case m 1 = m 2 , σ 1 = σ 2 , but α 11 = α 22 = α 12 . The fact that the coefficients of restitution are different entails that there is breakdown of energy equipartition here either. In other words, the partial temperatures of both species are different when they differ only in their coefficients of restitution. This situation has been widely considered for analyzing segregation driven only by inelasticity [85][86][87][88]. The temperature ratio T is plotted versus α 12 in Figure 6 for the above case when x 1 = 1 2 , φ = 0, α 11 = 0.9, and α 22 = 0.5. The theoretical results have been obtained by solving Equation (54) for getting the coefficients a  2 is small but not negligible at all since the agreement with simulations improves when these coefficients are considered in the evaluation of the partial cooling rates. We also find that energy nonequipartition is still significant in this particular situation in spite of the fact that the species have the same mass and diameter. It is quite apparent that although the application of the DSMC method to dilute systems is more efficient than the MD method from a computational point of view, the latter method avoids a crucial assumption of the former method: molecular chaos hypothesis (e.g., it neglects the possible velocity correlations between the particles that are about to collide). The study of the HCS of a granular binary mixture from MD simulations allows to prove the existence of the scaled solution (39) in a broader context than the kinetic theory (which is based on molecular chaos assumption). The HCS solution (39) with different partial temperatures determined by equating the partial cooling rates [Equation (40)] has been clearly confirmed by MD simulations [33]. The occurrence of this sort of solution appears for a wide range of volume fractions, concentrations, and mass and diameter ratios as well as for weak and strong inelasticity. In addition, the comparison between the results obtained from kinetic theory (approximate theoretical results and DSMC results) and MD simulations for the temperature ratio in several conditions may be considered as an stringent assessment of the reliability of kinetic theory. Figure 7a Figure 7 shows that the agreement between Enskog theory and MD simulations is very good for α = 0.95 over the complete range of mass ratios considered. Agreement is also good at α = 0.8 and φ = 0.1, although significant discrepancies between the Enskog equation (theory and DSMC results) and MD appear for large mass ratios at α = 0.8 and φ = 0.2. Regarding the dependence of T  As widely discussed in molecular mixtures [28,89], the Enskog equation has some limitations for describing systems at high densities. In this range of densities, one has to take into account recollision events (ring collisions) which go beyond the Enskog description. The impact of these multiparticle collisions on dynamic properties seems to be more stronger for inelastic collisions due to the fact that colliding pairs tend to be more focused. Thus, one expects that the range of densities where the Enskog equation for granular systems provides reliable predictions diminishes as inelasticity increases [90,91]. This is the trend observed here for the temperature ratio and also in other type of problems [63,92]. Apart from this limitation, another approximation employed here is the use of Equation (55) for estimating the pair correlation functions χ ij . In particular, recent MD simulations [44] at high densities have shown that χ 11 = χ 22 = χ 12 even when x 1 = 1 2 and σ 1 = σ 2 [the approximation (55) yields χ 11 = χ 22 = χ 12 for this case]. Thus, MD simulations have shown a value of χ 11 about a 20% larger than that of Equation (55), while χ 22 is, however, about 15% smaller. These differences quantify the effect of the spatial correlations on the Enskog prediction of the temperature ratio.
In conclusion, the comparison carried out in Figure 7 gives again support to the use of the Enskog equation for the description of granular flows across a wide range of densities, length scales, and inelasticity. Despite this success, the observed discrepancies between Enskog equation and MD simulations opens the necessity of developing kinetic theories that go beyond the Enskog theory. In any case, as has been discussed in several previous works [10], no other theory with such generality exists yet.
. In spite of this simple approximation, Figure 8 highlights the excellent agreement found between theory and Monte Carlo simulations, even for quite extreme dissipation. As for binary mixtures, the mean kinetic energy of the heavier species is larger than that of the lighter species. Moreover, the departure of the energy equipartition increases with the disparity in the mass ratios. As a complement of Figure 8, a moderately dense ternary mixture is considered in Figure 9. Here, x 1 = x 2 = 1 3 , φ = 0.1, m 1 /m 3 = 5, m 2 /m 3 = 2, σ 1 /σ 3 = (m 1 /m 3 ) 1/3 , and σ 2 /σ 3 = (m 2 /m 3 ) 1/3 . We observe that the effect of volume fraction on the temperature ratios does not change the main trends observed for dilute ternary mixtures. However, given that the diameter ratios are disparate in this case, more discrepancies between theory and DSMC results are found for small values of the coefficient of restitution, specially when m 1 /m 3 = 5. The presence of the second Sonine coefficients in the evaluation of T (0) i /T could mitigate in part these differences.

Navier-Stokes Transport Coefficients
We assume that we slightly disturb the HCS by small spatial perturbations. These perturbations induce nonzero contributions to the mass, momentum, and heat fluxes. The corresponding constitutive equations for the irreversible fluxes allow us to identify the relevant Navier-Stokes transport coefficients of the mixture. As for molecular mixtures, a reliable way of determining the transport coefficients is by means of the Chapman-Enskog method [11]. This method solves the set of Enskog equations by expanding the distribution function f i (r, v; t) of each species around the local version of the HCS (namely, the state obtained from the HCS by replacing the density, flow velocity, and temperature by their local values). The HCS state plays the same role for granular mixtures as the local equilibrium distribution for a molecular mixture (elastic collisions).
Therefore, as in the HCS, after a transient period one assumes that the distributions f i adopt the form of a normal solution. In other words, we assume that all space and time dependence of f i (r, v; t) only occurs through a functional dependence on the hydrodynamic fields: Functional dependence here means that to know f i at the point r, we need to know the values of the fields and all their spatial derivatives at r. For small spatial gradients, the functional dependence (65) can be made local in space through an expansion of f i in powers of the gradients of the hydrodynamic fields n i , U, and T: where the distribution f i (r, v; t) are chosen in such a way that their first few velocity moments give the exact hydrodynamic fields: Thus, the remaining distributions f (k) i must obey the constraints: It is important to note that in the expansion (66) we have assumed that the spatial gradients are decoupled from the coefficients of restitution. As a consequence, the Navier-Stokes hydrodynamic equations hold for small spatial gradients but they are not limited in principle to weak inelasticity. This point is relevant in the case of granular mixtures since there are some situations (e.g., steady states such as the uniform shear flow problem [94][95][96][97][98]) where hydrodynamic gradients are coupled to inelasticity and so, the Navier-Stokes approximation is restricted to nearly elastic spheres. Thus, due to the possible lack of scale separation for strong inelasticity, Serero et al. [85,86] consider two different perturbation parameters in the Chapman-Enskog solution: the hydrodynamic gradients (or equivalently, the Knudsen number Kn = /L, where is the mean free path and L is a characteristic hydrodynamic length) and the degree of dissipation ij = 1 − α 2 ij . The results derived from this perturbation scheme [85,86] agree with those obtained here in the quasielastic limit ( ij → 0).
Another important issue in the Chapman-Enskog expansion of granular mixtures is the choice of the hydrodynamic fields. Here, as for molecular mixtures [11,55,[99][100][101], we use the conserved number densities n i , the flow velocity U associated with the conserved total momentum, and the granular temperature T associated with the total kinetic energy. On the other hand, due to energy nonequipartition, other authors [48,49,[102][103][104][105][106][107] employ the set consisting of the conserved number densities n i , the species flow velocities U i associated with the non-conserved species momenta, and the partial (or species) temperatures T i . However, this choice is potentially confusing since, although more detailed, has no predictive value on the relevant hydrodynamic large space and time scales [108]. In particular, the two-temperature Chapman-Enskog solution considered in these works [48,49,[102][103][104][105][106][107] is phenomenological and assumes local Maxwellian distributions even for non-homogeneous situations. Although this approach yields vanishing Navier-Stokes transport coefficients for low-density mixtures, it can be considered as reliable to estimate the collisional transfer contributions to the irreversible fluxes [106,109].
The Chapman-Enskog solution to the (inelastic) Enskog Equation (1) was obtained in Refs. [53,54] some years ago. In particular, to first order in spatial gradients, the first-order velocity distribution function f (72) where the unknowns A i (V), B ij (V), C i,λβ (V), and D i (V) are functions of the peculiar velocity V. These quantities are the solutions of a set of coupled linear integral equations [53]. Approximate solutions to this set of integral equations were obtained [54,58] by considering the leading terms in a Sonine polynomial expansion. This procedure allows us to get the explicit forms of the Navier-Stokes transport coefficients in terms of the mechanical parameters of the mixture (masses and sizes and the coefficients of restitution), the composition, and the solid volume fraction. The constitutive equations for the mass j (1) i , momentum P (1) λβ , and heat q (1) fluxes have the form In Equations (73)- (75), D ij are the mutual diffusion coefficients, D T i are the thermal diffusion coefficients, η is the shear viscosity coefficient, η b is the bulk viscosity, κ is the thermal conductivity coefficient, and D q,ij are the partial contributions to the Dufour coefficients D q,i = ∑ j D q,ji . The Navier-Stokes transport coefficients associated with the mass flux are defined as The Navier-Stokes transport coefficients associated with the pressure tensor and the heat flux have kinetic and collisional contributions. Their kinetic contributions are given by η k b = 0, The expressions of the collisional contributions to η, η b , D q,ij , and κ can be determined from Equations (22) and (24) by expanding the distribution functions f i to first order in gradients. Their explicit forms can be found in Ref. [10]. We will go back to this point in Section 7 when we analyze the impact of different partial temperatures on the bulk viscosity coefficient.

Influence of the Temperature Ratios T (0)
i /T on the Transport Coefficients As mentioned before, the determination of the set of Navier-Stokes transport coefficients requires to know the functions A i (V), B ij (V), C i,λβ (V), and D i (V). As in the study of the HCS, the usual approach is to expand these unknowns in a series expansion of Sonine polynomials and consider only the leading terms. This procedure involves a quite long and tedious task where several collision integrals must be computed.
As expected, all the transport coefficients depend explicitly on the temperature ratios i /T, which are defined in terms of the zeroth-order distributions f (0) i . As discussed in Section 3, given that the form of f In this approximation, the (reduced) partial cooling rates ζ * i → ζ i is given by Equation (48).
The forms of the Navier-Stokes transport coefficients can be found in Ref. [10] when one uses the Maxwellian approximation (81). Their explicit expressions are very large and so, they are omitted here for the sake of brevity. On the other hand, for the sake of concreteness and to show in a clean way the impact of γ i on transport, we focus on our attention in this section in the coefficients D ij and D T i of a binary mixture (s = 2) in the low-density regime (φ = 0). The diffusion coefficients play for instance a relevant role in one of the most important applications in granular mixtures: segregation by thermal diffusion [110]. Since j where Here, 2 is given by Equation (48) with χ ij = 1. It is quite apparent from Equations (84)-(86) that the coefficients D * ij and D T * 1 depend in a complex way on the temperature ratio γ 1 [recall that γ 2 = x −1 2 (1 − x 1 γ 1 ) in a binary mixture]. To show more clearly the influence of energy nonequipartition on diffusion coefficients, it is convenient to write the forms of the above dimensionless coefficients by assuming energy equipartition. In this approximation (γ 1 = 1), θ 1 = 2µ 12 , θ 2 = 2µ 21 , Here, we recall that σ 12 = (σ 1 + σ 2 )/2 and ρ = m 1 n 1 + m 2 n 2 . Thus, taking into account Equations (87)-(89), the forms of D * ij and D T * 1 by assuming energy equipartition read Figure 10 shows the scaled diffusion coefficients D 11 (α)/D 11 (1), D 12 (α)/D 12 (1), and D T 1 (α)/D T 1 (1) versus the (common) coefficient of restitution α 11 = α 22 = α 12 ≡ α for a threedimensional dilute binary mixture with σ 1 /σ 2 = 1, x 1 = 1 2 , and two values of the mass ratio: m 1 /m 2 = 0.5 and 2. Here, D ij (1) and D T 1 (1) refer to the values of D ij and D T 1 for elastic collisions (α = 1). The expressions of D T 1 , D 11 and D 12 are provided by Equations (82)- (85). We observe first that the deviation of the diffusion D ij and thermal diffusion D T 1 coefficients with respect to their forms for elastic collisions (molecular mixtures) is in general significant, as expected. The departure from unity appears even for relatively moderate dissipation (let's say, α 0.8). Figure 10 also shows that the coefficients D ij and D T 1 (scaled with their elastic values) exhibit a monotonic dependence on inelasticity, regardless the value of the mass ratio: they increase with increasing dissipation (or equivalently, decreasing α). Thus, inelasticity enhances the mass transport of species. This monotonic behavior found for dilute mixtures is not kept at finite densities (φ = 0) since, depending on the value of the mass ratio, the scaled coefficient D 11 (α)/D 11 (1) may exhibit a non-monotonic dependence on α (see Figure 5.5 of Ref. [10]). An important target of Figure 10 is to illustrate the impact of energy non-equipartition on the transport coefficients. The dashed (for m 1 /m 2 = 0.5) and dash-dotted (for m 1 /m 2 = 4) lines refer to the results obtained for the three coefficients by assuming the equality of the partial temperatures (T 2 ). The expressions of these coefficients in this approximation are given by Equations (90) and (91). As said in the Introduction section of this review, most of the previous studies reported in the granular literature on mixtures were based on this equipartition assumption [24][25][26][27]85,86]. Figure 10 highlights the significant effect of energy nonequipartition on mass transport, specially for strong inelasticity. The impact of different partial temperatures (T (0) 2 ) on diffusion coefficients is not only quantitative but also in some cases qualitative. Thus, for instance, while D 11 (α) > D 11 (1) for m 1 /m 2 = 4 when energy nonequipartition is accounted for, the opposite [D 11 (α) < D 11 (1)] occurs when energy equipartition is assumed. A similar behavior exhibits the coefficient D 12 in the case m 1 /m 2 = 0.5. As expected, the important differences found between both theories (with and without energy equipartition) clearly shows that the effect of different species' granular temperatures cannot be neglected in the study of transport properties in granular mixtures. This conclusion contrasts with the results derived by Yoon and Jenkins [111] who conclude that segregation is not greatly affected by the difference in temperatures of the two species, at least when the particles of both species are nearly elastic and their masses or sizes do not differ by too much. On the other hand, other studies [39,110,[112][113][114][115][116][117] have shown the important influence of energy nonequipartition on segregation.
As a complement of Figure 10, we consider now the tracer limit x 1 → 0. In this limit case, D 22 ∝ x 1 and D T 1 ∝ x 1 and so, both coefficients vanish when one of the species of the mixture is present in tracer concentration. The expression of the tracer diffusion coefficient D 11 simply reads where in the tracer limit Figure 11 shows the α-dependence of the (scaled) tracer diffusion coefficient D 11 (α)/D 11 (1) for d = 3, σ 1 /σ 2 = 2, and m 1 /m 2 = 8. As in Figure 10, the influence of energy nonequipartition on D 11 is quite relevant, specially at strong inelasticity. Moreover, the comparison with the simulation results obtained from the DSCM method shows an excellent agreement, showing again the accuracy of the first Sonine approximation to D 11 . To end this section, the shear viscosity coefficient η is considered. Figure 12 shows the shear viscosity coefficient (scaled with respect to its value for elastic collisions) for a three-dimensional moderately dense binary mixture (φ = 0.1) with the same parameters as in Figure 10. Although the qualitative behavior of η(α)/η(1) is quite similar with and without energy equipartition (there is a monotonic decrease in shear viscosity as inelasticity increases in all the cases), there are important quantitative discrepancies between both theories specially in the case m 1 /m 2 = 4. To complement Figure 12, we plot in Figure 13 η(α)/η(1) for a two-dimensional (d = 2) dilute (φ = 0) binary mixture with x 1 = 1 2 and m 1 /m 2 = (σ 1 /σ 2 ) 2 (particles of the same mass density). We observe good agreement with Monte Carlo simulations when energy nonequipartiton is accounted for in the theory. Thus, as in the case of the diffusion coefficients and based on the findings of Figures 12 and 13, we can conclude that a reliable kinetic theory for granular mixtures needs to take into account nonequipartition effects in momentum transport.

First-Order Contributions to the Partial Temperatures: Influence on the Bulk Viscosity
As mentioned in Section 1, the presence of a divergence ∇ · U of the flow velocity in a mixture induces nonzero first-order contributions T (1) i to the partial temperatures. This breakdown of the energy equipartition is additional to the one appearing in the HCS which is only due to the inelastic character of the binary collisions. In fact, T (1) i = 0 even in the case of molecular dense mixtures, namely, a dense hard-sphere mixture with elastic collisions [55][56][57].
The fact that the partial temperatures T (1) i are proportional to ∇ · U gives rise to a contribution to the bulk viscosity η b coming from these temperatures. In addition, for granular mixtures, the temperatures T (1) i are also involved in the evaluation of the first-order contribution ζ U (proportionality coefficient between ζ and ∇ · U) to the cooling rate. The coupling between η b and T (1) i was already recognized by the pioneering works of the Enskog equation for multicomponent molecular gases [55][56][57].
According to the definition (10) of T i , its first-order contribution is where f (1) i (r, V; t) is given by Equation (72). Since T (1) i is a scalar, it can be only coupled to ∇ · U because ∇n and ∇T are vectors and the tensor ∂ λ U β + ∂ β U λ − (2/d)δ λβ ∇ · U is a traceless tensor. As a consequence, T (1) i can be written as where the scalar quantities D i (V) obey the following set of coupled linear integral equations [53]: Here, 2 is obtained from Equation (36) (97) and the homogeneous term In Equation (98), p * ≡ p/(nT) is the (reduced) hydrostatic pressure [p is given by Equation (34)], and the collision operator It is worth mentioning that in this section, it is understood that χ ij is evaluated at the zeroth-order approximation.
As said before, as a byproduct the calculation of T (1) i allows us to compute the firstorder contribution to the cooling rate where the coefficients ζ (1,1) and ζ (1,0) are defined by Equations (97) and (99), respectively. As in the case of the Navier-Stokes transport coefficients, the evaluation of the firstorder contributions T (1) i requires to solve the integral Equation (96). These equations can be approximately solved by considering the leading Sonine approximation to D i (V). Before taking this sort of approximation, it is convenient to prove the solubility condition (71), or equivalently, Upon writing the condition (102) we have taken into account that and consequently, the granular temperature T is not affected by the spatial gradients, as expected in the Chapman-Enskog method [11]. According to Equation (103), only s − 1 partial temperatures T In the low-density regime ( (98) leads to D i (V) = 0. Thus, the homogeneous term D i vanishes in the integral Equation (96) and so, D i = 0. This implies that the first-order contributions i to the partial temperatures vanish for dilute granular mixtures [50][51][52]85,86].

Bulk Viscosity Coefficient
The bulk viscosity η b is defined through the constitutive equation (74). This transport coefficient plays a relevant role in problems where the gas density varies in the flow motion; it represents an additional resistance to contraction or expansion. Since η b has only collisional contributions, its form can be identified by expanding the collisional transfer contribution (22) to the pressure tensor to first order in the spatial gradients. The expression of η b can be written as [119] η where and While the first contribution η b to the bulk viscosity is given in terms of the zeroth-order dis- i , the second contribution η b is given in terms of the first-order contributions i to the partial temperatures. Although this second contribution has been in fact neglected in several previous works [10,53,54] on dense granular mixtures, as said before it was already computed in the pioneering studies on molecular hard-spheres mixtures [56,57]. The impact of η b on η b will be assessed in the next subsection when we estimate i by taking the corresponding leading Sonine approximation to D i . Note that the expression (105) for the bulk viscosity can be written as where the forms of the partial shear viscosity coefficients η i b can be easily obtained from Equations (106) and (107). These forms could provide some insight into a shear-induced segregation problem.
An accurate estimate of the first contribution η b to the bulk viscosity is obtained by replacing f (0) i (V) by the Maxwellian distribution (81). With this approximation, η b is [53] where θ i is given in Equation (45).

Leading Sonine Approximation to i
The coefficient i is defined by Equation (95). To estimate it, we take the following Sonine approximation to D i (V): The coefficients i can be determined by substituting (110) into the integral equation (96), multiplying them with the polynomial W i (V), and integrating over the velocity. The procedure is large but straightforward. Technical details for multicomponent mixtures can be found in Ref. [119]. Here, we focus on the case of a binary mixture (s = 2). In this case, 2 = −(x 1 /x 2 ) 1 and 1 = (T/(nσ d−1 In Equation (111), we have introduced the dimensionless quantities where we recall that m = (m 1 + m 2 )/2 for a binary mixture.
Another simple but interesting case corresponds to molecular mixtures of dense hardspheres. In this case (α 11 12 , θ 2 = 2µ 21 , and * 1 becomes * where the expressions of ω * 11,el and ω * 12,el are easily obtained from Equations (113) and (114), respectively, by considering elastic collisions. The expression (115) agrees with the one obtained many years ago by Karkheck and Stell [57] for a hard-sphere binary mixture (d = 3). On the other hand, for a two-dimensional system (d = 2), Equation (115) differs from the one derived by Jenkins and Mancini [47] for nearly elastic hard disks. As recognized by the authors of this paper, given that their prediction on * 1 was derived by assuming Maxwellian distributions for each species, a more accurate expression of * 1 is obtained when one evaluates this coefficient from the first-order distribution of the Chapman-Enskog solution. In particular, Equation (115) takes into account not only the different centers r and r ± σ of the colliding spheres in the Enskog collision operator (this is in fact the only ingredient accounted for in Ref. [47] for getting * 1 ) but also the form of the first-order distribution functions f (1) i given by Equation (72). Moreover, while * 1 → 0 for vanishing densities (φ → 0), the results found by Jenkins and Mancini [47] predict a nonvanishing * 1 for dilute binary mixtures if m 1 = m 2 . This result contrasts with those obtained for molecular mixtures [56,57].
To illustrate the differences between the results obtained in Ref. [47] and those derived here for disks, Figure 14 shows * 1 versus m 1 /m 2 when x 1 = 1 2 , φ = 0.25 and σ 1 /σ 2 = (m 1 /m 2 ) 1/2 (i.e., when the disks are made of the same material). In the case of disks [47], where φ = ∑ i n i πσ 2 i /4 is the solid volume fraction for disks and we recall that M = ∑ i x i σ i . It is quite apparent the differences found between both theories, specially for disparate masses.
For inelastic collisions, Figure 15 illustrates the dependence of * 1 on the (common) coefficient of restitution α for a binary mixture of hard spheres (d = 3) with x 1 = 1 2 , φ = 0.25, m 1 /m 2 = (σ 1 /σ 2 ) 3 and three different values of the mass ratio: m 1 /m 2 = 0.5, 2 and 5. We observe first that * 1 is significantly affected by inelasticity, specially for high mass ratios. With respect to the effect of the mass ratio on * 1 , we see that this coefficient decreases (increases) with increasing inelasticity when m 1 /m 2 > 1 (m 1 /m 2 < 1). As expected, Figure 15 also shows that the magnitude of * 1 is in general quite small in comparison with the remaining transport coefficients.

Influence of T
(1) i on the Bulk Viscosity and the Cooling Rate According to Equations (105)-(107), the coefficient * 1 is involved in the contribution η b to the bulk viscosity η b . To assess the impact of the first-order contributions to the partial temperatures on the bulk viscosity, we plot in Figure 16 the (reduced) bulk viscosity η b (α)/η b (1) as a function of the (common) coefficient of restitution α. As in Figure 15, x 1 = 1 2 , φ = 0.25 and m 1 /m 2 = (σ 1 /σ 2 ) 3 . Two different mass ratios are studied: m 1 /m 2 = 0.5 and 5. The value of the (reduced) bulk viscosity when the coefficient * 1 is neglected (dashed lines) is also plotted for the sake of comparison. Although both results (with and without the contribution coming from η b ) agree qualitatively, Figure 16 highlights that the impact of * 1 on the bulk viscosity cannot be neglected for high mass ratios and strong dissipation (let's say, for instance, α 0.5). Finally, Figure 17 shows the α-dependence of the first-order contribution ζ U to the cooling rate. This coefficient is defined by Equation (101) where ζ (1,1) is The coefficients ξ * 1 and ξ * 2 are given by Equation (112). As Figures 16 and 17 highlights that the influence of * 1 turns out to be relevant for strong inelasticities and high mass ratios.  Figure 17. Dependence of the magnitude of the (reduced) first-order contribution ζ U to the cooling rate on the (common) coefficient of restitution α for a granular binary mixture of hard spheres (d = 3) with x 1 = 1 2 , φ = 0.25, and σ 1 /σ 2 = 2, and two different values of the mass ratio: m 1 /m 2 = 0.5 (a) and m 1 /m 2 = 5 (b). The solid lines are the results obtained here while the dashed lines correspond to the results obtained for the (reduced) cooling rate when the contribution ζ (1,1) to ζ U is neglected.

Summary and Concluding Remarks
The primary objective of this review has been to analyze the influence of energy nonequipartition on the transport coefficients of an s-component granular mixture. Granular mixtures have been modeled here as a collection of inelastic hard spheres of masses m i and diameters σ i (i = 1, · · · , s). We have also assumed that spheres are completely smooth so that the inelasticity of collisions is only accounted for by the (positive) constant coefficients of normal restitution α ij ≤ 1. At a kinetic level, all the relevant information on the state of the mixture is given through the knowledge of the one-particle velocity distribution functions f i (r, v; t) of each species. At moderate densities, the distributions f i verify the set of s-coupled Enskog kinetic equations.
The study of the influence of different partial temperatures on transport has been carried out in two different steps. First, we have widely analyzed the failure of energy equipartition in granular mixtures in the HCS, namely, a homogeneous freely cooling state. The understanding of this simple situation is crucial because the HCS plays the role of the reference state in the Chapman-Enskog solution to the Enskog equation. Assuming the scaling solution (39) for the distributions f i , the temperature ratios γ i ≡ T To estimate the partial cooling rates ζ i , the leading Sonine approximation (44) to the scaled distributions ϕ i have been considered. This approximation also involves the calculation of the second Sonine coefficients a (i)

.
The temperature ratios γ i and the Sonine coefficients a These theoretical predictions have been tested via a comparison with Monte Carlo (DSMC) and MD simulations for conditions of practical interest. Comparison between DSMC simulations and theory shows in general an excellent agreement; more discrepancies are observed in the case of MD simulations, specially for high volume fractions and/or strong dissipation. This disagreement is a clear indication of the limitations of the Enskog theory in these ranges of values of volume fraction (or density) and/or inelasticity. On the other hand, the good agreement found with the DSMC results reinforces the reliability and accuracy of the approximate analytical predictions even for disparate mass and diameter ratios and/or small values of the coefficients of restitution. As expected, the deviations from the energy equipartition (T (0) i /T = 1) can be weak or strong depending on the mechanical differences between the different species of the mixture and the degree of inelasticity in collisions.
Once the dependence of the ratios T (0) i /T on the parameter space of the mixture has been characterized, the next step has been to study the influence of nonequipartition on the Navier-Stokes transport coefficients. This study is relevant since many of previous attempts reported in the literature [24][25][26][27] for obtaining the transport coefficients of granular mixtures had assumed the equality of the partial temperatures (T Thus, in contrast with the conclusions reached in previous works [111], our analysis shows that the impact of different partial temperatures on transport is in general quite significant, as has been clearly illustrated in Figures 10-13 for the diffusion and shear viscosity coefficients. As a second step in the paper, we have also analyzed the failure of energy equipartition due to the presence of spatial gradients in the system. More specifically, in a dense mixture a nonzero divergence ∇ · U of the flow velocity field induces a nonzero contribution to the partial temperatures T i . This first-order contribution T (1) i to T i is not generic of dense granular mixtures since it is also present for elastic collisions. In fact, previous pioneering works [55][56][57] on dense hard-sphere molecular mixtures (α ij = 1) determine these coefficients in terms of the parameters of the mixture. Here, we have extended those calculations to the case of granular mixtures (α ij < 1).
A careful analysis of the first-order Chapman-Enskog solution to the Enskog equation shows that the coefficients T (i) i are involved in the evaluation of the bulk viscosity η b (proportionality coefficient between the collisional part P c of the pressure tensor and ∇ · U) and the first-order contribution ζ U to the cooling rate ζ (proportionality coefficient between ζ and ∇ · U). Thus, although the coefficients T Our results indicate that the effect of T (1) i on both η b and ζ U is only relevant for high mass ratios and strong dissipation (see Figures 16 and 17). In this context, we can conclude that previous expressions of the bulk viscosity and cooling rate for dense granular mixtures [53,54,58] (which implicitly neglect the first-order contributions T (1) i ) must incorporate the contributions coming from T (1) i when the masses of the species are disparate and/or the degree of collisional dissipation turns out to be important. Granular hydrodynamics derived from hard-sphere models have been shown to be useful in the description of numerous industrial processes involving solid particles. Of particular relevance are high-speed, gas-solid flows, and fluidized beds. Such descriptions are now standard features of commercial and research codes. Since those codes rely upon accurate expressions of the Navier-Stokes transport coefficients, it is quite apparent that a first-order objective is to guarantee a reliable theoretical treatment. As shown in this paper and previous review works [10,123], the price of this accurate approach (in contrast to more phenomenological approaches) is an increasing complexity of the expressions derived for the transport coefficients.
As mentioned in Section 1, since grains in nature are generally surrounded by a fluid like water or air, a granular mixture is in fact a multiphase system. In this review, the influence of the interstitial fluid on the dynamical properties of the granular mixture has been neglected. A further step is to take into account the presence of the surrounding gas and develop a theory for moderately dense granular suspensions. This will provide a fundamental basis for the application of granular hydrodynamics under realistic conditions. Although some previous attempts based on the introduction of solid-fluid forces [120,[124][125][126] have been made in the past, it still remains to propose theories where the influence of collisions between the solid particles and molecules of the interestitial fluid is explicitly accounted for in the corresponding kinetic equation. The complexity for considering this type of collisions is not a real problem for implementation in a code.
As shown along this overview, granular mixtures exhibit a wide range of interesting phenomena for which the Navier-Stokes hydrodynamic equations can be considered as an accurate and practical tool. However, due to their complexity, many of their features are not fully understood. Kinetic theory and hydrodynamics (in the broader sense) can be expected to provide some insight into the understanding of such complex materials.
Author Contributions: Conceptualization, V.G.; software, M.G.C.; validation, M.G.C., R.G.G. and V.G.; formal analysis, R.G.G. and V.G.; writing-original draft preparation, V.G.; writing-review and editing, M.G.C., R.G.G. and V.G. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors acknowledge financial support from Grant PID2020-112936GB-I00 funded by MCIN/AEI/ 10.13039/501100011033, and from Grants IB20079 and GR18079 funded by Junta de Extremadura (Spain) and by ERDF A way of making Europe. The research of R.G.G. also has been supported by the predoctoral fellowship BES-2017-079725 from the Spanish Government.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors report no conflict of interest.