Rheology of a dilute binary mixture of inertial suspension under simple shear flow

The rheology of a dilute binary mixture of inertial suspension under simple shear flow is analyzed in the context of the Boltzmann kinetic equation. The effect of the surrounding viscous gas on the solid particles is accounted for by means of a deterministic viscous drag force plus a stochastic Langevin-like term defined in terms of the environmental temperature $T_\text{env}$. Grad's moment method is employed to determine the temperature ratio and the pressure tensor in terms of the coefficients of restitution, concentration, the masses and diameters of the components of the mixture, and the environmental temperature. Analytical results are compared against event-driven Langevin simulations for mixtures of hard spheres with the same mass density $m_1/m_2=(\sigma^{(1)}/\sigma^{(2)})^3$, $m_i$ and $\sigma^{(1)}$ being the mass and diameter, respectively, of the species $i$. It is confirmed that the theoretical predictions agree with simulations of various size ratios $\sigma^{(1)}/\sigma^{(2)}$ and for elastic and inelastic collisions in the wide range of parameters' space. It is remarkable that the temperature ratio $T_1/T_2$ and the viscosity ratio $\eta_1/\eta_2$ ($\eta_i$ being the partial contribution of the species $i$ to the total shear viscosity $\eta=\eta_1+\eta_2$) discontinuously change at a certain shear rate as the size ratio increases; this feature (which is expected to occur in the thermodynamic limit) cannot be completely captured by simulations due to small system size. In addition, a Bhatnagar--Gross--Krook (BGK)-type kinetic model adapted to mixtures of inelastic hard spheres is exactly solved when $T_\text{env}$ is much smaller than the kinetic temperature $T$. A comparison between the velocity distribution functions obtained from Grad's method, BGK model, and simulations is carried out.


Introduction
Rheology is the subject that studies the flow properties of materials.Although the viscosity of the Newtonian fluid is independent of the shear rate, there are many domestic substances (liquids containing microstructures such as suspensions and polymers) where the viscosity depends on the shear rate (non-Newtonian fluids).Within the class of non-Newtonian fluids, some of them exhibit shear thinning (namely, when the viscosity decreases with the shear rate) while others display shear thickening (namely, when the viscosity increases with the shear rate).The shear thickening is also categorized into two classes as the continuous shear thickening (CST) and discontinuous shear thickening (DST).The viscosity increases continuously in CST, while it abruptly changes discontinuously from a small value to a large value at a critical shear rate in DST.DST has attracted much attention among physicists in the last few years [1][2][3][4][5][6][7] as a typical nonequilibrium discontinuous phase transition between a liquid-like phase and a solid-like phase.In addition, the understanding of the origin of DST is also important for potential industrial applications such as protective vests and traction controls.
Although most of the previous studies on shear thickening have been oriented to dense suspensions, there are some other studies that analyze a DST-like process for the kinetic temperature of inertial suspensions.This type of inertial suspensions can be regarded as an idealistic model of aerosols [8].The DST-like process (or the ignited-quenched transition) of dilute inertial suspensions takes place as a result of a saddle-node bifurcation.On the other hand, the DST-like process for dilute suspensions becomes CST-like as the density of suspensions increases [9][10][11][12][13][14][15][16].
To gain some insight into the understanding of the generic features of rheological phase transitions, we use kinetic theory tools in this paper.This allows us to offer a quantitative theoretical analysis for the DST-like and CST-like processes in inertial suspensions.However, it should be noted that some previous kinetic theories for inertial suspensions have ignored thermal fluctuations in the dynamics of grains [9][10][11]14].A refined suspension model including a Langevin-like term has been more recently considered in Refs.[12,13,[15][16][17].The quantitative validity of these studies has already been verified by the event-driven Langevin simulation for hard spheres (EDLSHS) [15,16,18].
Most of the previous theoretical studies on the rheology of inertial suspensions have focused on monodisperse systems, namely, suspensions containing only identical spherical particles.In reality, suspended particles are not identical since the size of the particles is distributed and the shape and mechanical properties of the particles are also different.To quantify the impact of polydispersity on the rheological properties of inertial suspensions under simple or uniform shear flows (USF), we consider a binary mixture in this paper, namely, a suspension which contains two kinds of spherical particles having different sizes.We note that bidisperse systems are also studied in colloidal and blood suspensions [19][20][21][22].
A challenging and interesting problem in sheared granular binary mixtures is that of the diffusion.It is well established that in the absence of shear the mass flux is proportional to the density, pressure, and temperature gradients where the corresponding transport coefficients are scalar quantities [23].However, when the mixture is strongly sheared, due to the anisotropy induced by the shear flow tensorial quantities are required to characterize the mass transport instead of the conventional scalar diffusion coefficients.There have been some previous attempts for describing the self-diffusion problem in sheared granular mixtures [24,25].As expected, all previous studies indicate that the diffusion process in USF is highly anisotropic and the components of the diffusion can be observed in the directions parallel and perpendicular to the velocity gradient.To characterize such anisotropy of the diffusion tensor, there have been several theoretical studies based on kinetic theory [26,27], 2/39 simulation works of rapid granular shear flows [28,29], and experimental studies of dense, granular shear flows in a two-dimensional Couette geometry [30,31].
One of the key features of flows of polydisperse particles is segregation [32].This is likely one of the most relevant problems in granular mixtures, from practical and fundamental points of view.However, despite many industrial and scientific progresses made in the past few years, the mechanisms involved in the segregation phenomenon are still poorly understood.In particular, in the context of kinetic theory, many different papers have addressed the study of segregation [33][34][35][36][37][38][39][40][41][42].On the other hand, computer simulations of bidisperse granular mixtures under USF (and without any influence of gravity) [43] have not found any sign of large-scale size segregation.Another type of works has studied segregation in flows down inclined slopes in which approximate simple shear flows have been realized, at least, in the bulk regions away from the bottom boundaries and surfaces.It is remarkable that the trigger of the segregation is the deviation from the USF of the velocity profile as has been reported in Ref. [42].This suggests that segregation can be observed even for dilute mixtures without the influence of gravity, if we drive a shear flow through a boundary.In other words, the segregation is localized near the boundaries.
Previous studies of granular binary mixtures based on the kinetic theory have mainly focused on obtaining the Navier-Stokes transport coefficients of the mixture by considering states close to the homogeneous cooling state [23] and/or close to driven stationary homogeneous states [44][45][46].The results are more scarce in the study of the rheological properties of granular binary mixtures under USF [43,[47][48][49].As expected, the results show that the mixture is non-Newtonian and in some cases, the effect of bidispersity enhances the non-Newtonian character of the fluid.Since the USF is spatially homogeneous in the frame moving with the linear velocity field, no segregation appears in the system.However, when the USF state is slightly perturbed by small density and temperature gradients, a nonvanishing mass flux is present and the corresponding components of the diffusion tensors have been determined in the tracer limit in Refs.[26,27].The knowledge of the shear-rate dependence of the above diffusion tensors has allowed to analyze thermal diffusion segregation induced by the presence of a temperature gradient orthogonal to the shear flow plane [50].
Nevertheless, so far and to the best of our knowledge, there are few studies of binary mixtures of inertial suspensions including diffusion processes, in which the rheology of inertial suspensions drastically depends on the shear rate.Thus, as already did in the case of granular mixtures [26,27], one has to determine the rheological properties of sheared binary mixtures of inertial suspensions as a first step before considering the segregation problem.Once rheology is known, the components of the diffusion tensors can be determined by using a similar procedure as the one followed in (dry) granular mixtures.Therefore, the study of the rheology of a dilute binary mixture of inertial suspension is an important issue.
Beyond dilute granular flows, it is quite apparent that there are many exotic rheological processes in dense flows.These processes include glass transitions, shear jamming, jamming, and DST [1][2][3][4][5][6][7][51][52][53][54].Such exotic processes cannot be observed in monodisperse systems but they can be observed only in mixtures when the volume fraction exceeds the transition point for crystallization of identical spheres at the volume fraction ϕ = 0.49.

3/39
In this paper, we focus on the rheology of a dilute binary mixture under USF.As in our previous works [15,16], the influence of the interstitial gas on solid particles is accounted for in an effective way by means of (i) a deterministic drag force proportional to the particle velocity and (ii) a stochastic Langevin-like term.While the first contribution attempts to model the friction of grains on the viscous fluid (a collection of gas molecules), the second term mimics the energy gained by the solid particles due to their interactions with the particles of the surrounding gas.The corresponding set of two coupled Boltzmann kinetic equations is solved by two complementary and independent routes: (i) Grad's moment method and (ii) event-driven simulations for hard spheres (EDLSHS).The comparison between kinetic theory and EDLSHS allows us to verify the reliability of the theoretical predictions as the first step to tackle the behavior of sheared binary mixtures of inertial suspensions.Our (approximate) analytical results of the rheological properties of the mixture (the ratio T 1 /T 2 between the partial temperatures and the pressure tensor) agree well with simulations for conditions of practical interest.In particular, the temperature ratio T 1 /T 2 and the viscosity ratio η 1 /η 2 (where η i is the partial contribution of the component i to the total shear viscosity η = η 1 + η 2 ) exhibit a DST-like transition for sufficiently high values of the size ratio σ (1) /σ (2) .As a complement, we have also compared the velocity distribution function obtained by both Grad's moment method and a kinetic model with the one obtained by EDLSHS.
The contents of the paper are as follows.In Sect.2, we introduce the Langevin model and Boltzmann equation for a binary mixture of inertial suspensions under a simple shear.Section 3 deals with the theoretical procedure to derive the rheology of inertial suspensions in USF.Section 4 is the main part of this paper, in which we present the theoretical and numerical results and find a new rheological phase transition similar to DST.The velocity distribution function is also studied by comparing the results from Grad's approximation and simulations.In Sect.5, we discuss and conclude our results.Moreover, there are several appendices to explain the technical details of the paper.In Appendix A, the difference between P (i) yy and P (i) zz is discussed.In Appendix B, we provide some mathematical steps to compute the collisional moment needed to determine the components of the pressure tensor from Grad's method.The detailed rheological properties for the temperature ratio, temperature, and viscosity are discussed in Appendix C. Appendix D discusses how the discontinuous transition appears/disappears when we change the parameters of the mixture.The tracer limit of the theory is briefly presented in Appendix E while Appendix F gives the exact solution to a Bhatnagar-Gross-Krook (BGK)-like kinetic model for granular mixtures in the high shear rate regime.This solution provides a two-dimensional velocity distribution function.Finally, the one-dimensional velocity distribution function is displayed in Appendix G with a comparison with the one obtained from computer simulations.

Basic equations for a binary mixture of inertial suspension under uniform Shear Flows
In this section, we present the basic equations describing a dilute binary mixture of inertial suspensions under USF.In the first subsection, we introduce the Langevin equation characterizing the motion of each particle activated by the thermal noise caused by collisions with the environmental molecules.In the second subsection, we write the corresponding set of 4/39 two coupled nonlinear Boltzmann kinetic equations for the bidisperse inertial suspension in the low-density regime.

Langevin equation
We consider a three-dimensional binary mixture of inertial suspension modeled as a mixture of inelastic hard spheres of masses m i and diameters σ (i) (i = 1, 2).For the sake of simplicity, we assume that the spheres are completely smooth and hence, collisions among all pairs are characterized by (positive) constant coefficients of normal restitution e ij ≤ 1, where the subscripts ij denote the species i and j, respectively.Let us use the notations v 1 and v (j) 2 when the particle 1 (species i) collides with the particle 2 (species j).The post-collisional velocities v of particles 1 (species i) and v for 2 (species j) are expressed as where we have introduced the pre-collisional velocities of particles v 1 for 1 (species i) and 2 (species j), v 2 , the unit normal vector at contact σ, and the reduced mass Fig. 1 Setup of our system.Two species of particles are distributed in a fluidized inertial suspension characterized by the temperature T env .The shear is applied with the shear rate γ in the xy plane, where the x and y axes are the shear direction and the velocity gradient direction, respectively.Here, we use N = 30000 particles with the size and number ratio as σ (1) /σ (2) = 10.0 and N 1 /N 2 = 30/29970 = 1/999, respectively.
The inertial suspension we consider is subjected to a steady simple shear flow in the x direction as shown in Fig. 1.The equation of motion for the k-th particle of the species i is described by the Langevin equation where ζ i is the drag coefficient acting on the particle of species i from the environmental fluid, and p is the peculiar momentum of the k-th particle with velocity v Here, γ and e x are the shear rate and unit vector in the sheared (x) direction, respectively.If hard-core grains are subjected to the Stokes' drag, ζ i is simply proportional to σ (i) and √ T env , where T env is the environmental temperature.When we adopt the mean diameter σ ≡ (σ (1) + σ (2) )/2 and drag coefficient For denser flows, the dependence of ζ i on the parameters of the mixture is more complex [55,56].In Eq. ( 2), F imp k expresses the impulsive force accounting for the collisions while the noise term ξ k,α e α (the unit vector e α in the α−direction) satisfies the fluctuation-dissipation relation [57]:

Boltzmann equation
If the density of the solid particles is low enough, the Langevin equation ( 2) can be converted into the Boltzmann kinetic equation for the distribution function f i (r, v, t) for the species i of the dilute binary mixture of inertial suspensions.The set of coupled Boltzmann equations read with the collision integral [23] J where we have introduced σ 2 )/2.From the distribution f i , one can define the number density of species i as the flow velocity U i of species i as and the partial temperature T i of species i Here, V (r, t) ≡ v − U (r, t) is the peculiar velocity.The mean flow velocity U (r, t) and the kinetic temperature T (r, t) are defined, respectively, as where n ≡ n 1 + n 2 is the total number density, ρ i ≡ m i n i is the mass density of species i, ρ ≡ ρ 1 + ρ 2 is the total mass density, and ν i ≡ n i /n = N i /N is the fraction of species i.Here, N i is the number of particles of species i and N = N 1 + N 2 .

6/39
Let us consider the macroscopic velocity satisfying where γ is the constant shear rate.In terms of the peculiar velocity V, Eq. ( 4) can be rewritten as At a macroscopic level, the USF is characterized by uniform density and temperature and the linear velocity field (10).In addition, at a more fundamental level, the USF is defined as that which is spatially uniform in the Lagrangian frame moving with the velocity field (10).In this frame, f i (r, v, t) = f i (V , t) and Eq. ( 11) is reduced to the equation for the velocity distribution function: One of our theoretical goals is to determine the pressure tensor where the partial pressure tensor for species i is defined as We use the Greek and Latin characters for {x, y, z} and the species {1, 2}, respectively.The knowledge of the pressure tensor allows one to get the rheological properties of the inertial suspension.Needless to say, the flow under USF is independent of the spatial position by its definition.Therefore, we can start from Eq. ( 12) as the basic equation for the theoretical analysis.

Rheology under uniform shear flows
In this section, we present the results of rheology for a binary mixture of inertial suspension under USF obtained from the Boltzmann equation (12).There are three subsections in this Section.In the first subsection, we summarize a general framework to determine the rheology of inertial suspensions by deriving a set of equations for the partial pressure tensors by multiplying both sides of Eq. ( 12) by m i V V and integrating over velocity.No approximations are considered in this subsection, including the moment of the collision integral (5).In the second subsection, we focus on the steady rheology within the framework of the kinetic theory under Grad's moment method [58].In the third subsection, we explain the concrete procedure to determine the steady rheology.

Moment equation for the pressure tensor
As mentioned before, the set of equations for the pressure tensor of the species i is obtained by multiplying by m i V α V β both sides of the Boltzmann equation (12) and integrating over 7/39 V .The result is where Let us introduce the anisotropic temperatures It should be noted that there are some other anisotropic temperatures such as δT ≡ (P xx − P zz )/n, which differ from ∆T in general.Nevertheless, we ignore the difference between ∆T and δT in this paper, because (i) the detection of the difference between ∆T and δT is difficult [13], and (ii) the linear approximation of Grad's method used later yields P (i) and so, ∆T = δT .In general, δT differs from ∆T even for dilute systems (Refs.[10,11]), although the previous studies confirmed that the effect of δT = ∆T is small [11,13].We want also to indicate that the difference between ∆T and δT is almost imperceptible in the simulations (see Appendix A) despite the requirement of long and tedious calculations for evaluating this difference [13].Therefore, for simplicity, we ignore the difference between ∆T and δT in this paper.
If we adopt such a simplification, the evolution equations for T i , ∆T i , and xy are given by where we have introduced the environmental temperature T env of the suspension liquid.Note that the diagonal elements of the pressure tensor in dilute systems can be written as In this paper, we adopt Einstein's rule for the summation, i.e., P αα .Upon deriving Eqs.(19), we recall that we have made use of the identity P (i)

Kinetic theory of rheology for a dilute binary mixture of inertial suspension via
Grad's method 3.2.1.Grad's moment method for the velocity distribution function.To determine αβ , we need to know the explicit form of the collisional moments Λ (ij) αβ .This requires the knowledge of 8/39 the velocity distribution functions f i , which is an intricate problem even for the elastic case.As for monodisperse inertial suspensions, a useful way to estimate Λ (ij) is to adopt Grad's moment method [58] in which the true f i is approximated by the trial Grad's distribution [11,12,15,16,26,[58][59][60]: where and f i,M (V ) is the Maxwellian distribution at the temperature T i of the species i, i.e., Note that in Eq. ( 21) we have taken into account that the mass and heat fluxes of a binary mixture vanish in the USF state.
With the use of the distribution ( 21), the integrals appearing in the expression of the collisional moments Λ (ij) can be explicitly computed.After a lengthy calculation (see Appendix B for the derivation), one gets with Here, we have introduced ǫ i ≡ m i T /(mT i ), and the thermal velocity v T ≡ 2T /m with the mean mass m defined as m ≡ (m 1 + m 2 )/2.Upon deriving Eq. ( 24), nonlinear terms in the traceless stress tensor Π (i) αβ have been neglected (linear Grad's approximation).The expression (24) agrees with a previous derivation of the collision integral Λ (ij) αβ [47].Now, let us rewrite the set of equations ( 18) in dimensionless forms.We introduce the partial dimensionless temperatures θ i and the anisotropic temperatures ∆θ i for species i as We also introduce the global quantities θ ≡ 2 i=1 ν i θ i and ∆θ ≡ 2 i=1 ν i ∆θ i , where we recall that Then, the dimensionless collisional moment becomes with where we have introduced the packing fraction In addition, the dimensionless reduced mass characterizes the magnitude of the noise [13].In terms of the temperature ratio where we have introduced the scaled time τ ≡ t T env /m/σ, the dimensionless shear rate γ * ≡ γσ m/T env = γ/(ζξ env ), and the dimensionless drag coefficient For the sake of convenience, some explicit forms of Λ (ij) * αβ in Eq. ( 24) can be rewritten as where

Theoretical expressions in the steady rheology
Although the set of Eqs. ( 36) applies for time-dependent states, we are mainly interested in the rheology in the steady state.Hereafter, we focus on the steady rheology.
3.3.1.Theoretical flow curves in steady state.Let us solve the set of Eqs.(36) in the steady state.In this case (∂ τ = 0), the left hand side of the set (36) vanishes and one gets First, from Eq. (39a), one obtains Substituting Eq. ( 40) into Eq.(39b), a set of equations which determine ∆θ 1 and ∆θ 2 can be rewritten as for i = 1, 2. Here, we have introduced the quantities 11)  xy − Λ ′ (11)   xy 22)  xy − Λ ′ (22)   xy and Then, ∆θ 1 and ∆θ 2 can be expressed as Substituting Eqs. ( 40) and ( 44) into Eq.(39c) leads to the relationship where and Equation ( 45) determines the relationship between the (reduced) global kinetic temperature θ and the temperature ratio ϑ.For given values of the mixture and at a given value of θ, we determine the temperature ratio ϑ [defined in Eq. ( 32)] by numerically solving Eq. ( 45).As will be shown in the next section, we find a solution of θ by fixing ϑ in the intermediate shear regime where the size ratio becomes large by fixing the volume ratio.Once we determine this relationship, we can draw the flow curve with the aid of Eq. ( 45), where the shear rate is given by

Comparison between theory and simulation
In this section, we compare the theoretical results obtained in sec.3 with those of EDLSHS [18].We will consider binary mixtures constituted by species of the same mass density [m 1 /m 2 = (σ (1) /σ (2) ) 3 ] and a (common) coefficient of restitution [e 11 = e 12 = e 21 = e 22 ≡ e].
In the first subsection, we examine the case of an equimolar mixture (ν 1 = ν 2 = 1/2 or N 1 = N 2 ) while the general case of N 1 = N 2 will be analyzed in the second subsection.In particular, we find a new DST-like rheological phase transition for N 1 ≪ N 2 when we fix the volume ratio, i.e., a binary mixture in which the concentration of one of the species (the large tracer particles 1) is much more smaller than that of the other species (the small particles 2).

The rheology for
In this subsection, we present the results of EDLSHS to verify the validity of the predictions of the kinetic theory for N 1 = N 2 .In this case, we should note that the occupied volume is dominated by the large grains for a large size ratio σ (1) /σ (2) .For example, the ratio of occupied volumes V ≡ N 1 σ (1)3 /(N 2 σ (2)3 ) becomes 125 if we adopt σ (1) /σ (2) = 5.0.The results of EDLSHS under the control of N 1 /N 2 with fixing V will be discussed in the next subsection.For the comparison of the theoretical results with those of EDLSHS, we have used the steady solutions of Eqs.(39) for both elastic (e = 1) and inelastic (e = 0.5, 0.7, and 0.9) cases when we fix N = 1000, ϕ = 0.01, ξ env = 1.0, and 12/39 Figures 2 (for ϑ = T 1 /T 2 ) and 3 (for η 1 /η 2 ) show some characteristic rheological flow curves for binary mixtures for both elastic and inelastic cases.Here, we have introduced the viscosity η i for species i as Now, let us focus on the plot of the temperature ratio ϑ ≡ T 1 /T 2 against the reduced shear rate γ * in Fig. 2. In the low shear regime, the temperature ratio converges to unity as shown in Fig. 2.This is because the temperatures of both the larger and smaller particles are determined by the thermal noise of the background fluid.On the other hand, the temperature ratio converges to a constant in the high shear regime, which is determined by the interparticle inelastic collisions.Note that this converged value agrees with the one previously obtained for granular gases [47].Interestingly, our theory predicts the existence of a negative peak for ϑ in the intermediate shear regime.In particular, there exists a cusp for smaller values of σ (1) /σ (2) at a certain shear rate γ * c (at which |∂ϑ/∂ γ * | → ∞; see Fig. 2).Correspondingly, the ratios of the other quantities such as η 1 /η 2 exhibit cusps around γ * c (see Fig. 3).It is worth remarking that these observables exhibit common features since (i) they do not have sharp minima even near the DST-like transition point of one of two species, (ii) the deviations from unity become larger with increasing the size ratio σ (1) /σ (2) , and (iii) the ratios converge to values different from unity even in the low-shear limit.These singularities are inherently connected with the DSTlike transition observed (see Appendix C) for the global kinetic temperature θ and the shear viscosity with Eq. ( 49)] because the cusps vanish as the size ratio increases.Indeed, the flow curves for ϑ and η 1 /η 2 become smooth for σ (1) /σ (2) = 5.0 (see Figs. 2 and 3).Let us also discuss the existence of cusps in the flow curves observed in Figs. 2 and 3 when the size ratio is small.As shown in Fig. 4, the partial viscosities η * i also have discontinuous jumps when the mean viscosity η * = η * 1 + η * 2 has also this jump.At points ( γc ) where ∂η * i /∂ γ * → ±∞ (i = 1, 2) are satisfied, the viscosity ratio also diverges as This explains the reason for the existence of the cusps.
It is remarkable that the predictions of kinetic theory agree well with the simulation results of EDLSHS without any fitting parameter.Therefore, we conclude that our kinetic theory based on the Boltzmann equation with Grad's method is reliable to describe the rheology, at least, for N 1 = N 2 .
To close this subsection, we also note that the flow curves become discontinuous and continuous depending on the other parameters of the mixture.These behaviors are discussed in Appendix D.

The rheology for
In this subsection, we compare the simulation results for the rheology for N 1 = N 2 with those derived from the theoretical predictions by fixing the volume ratio V ≡ N 1 σ (1)3 /(N 2 σ (2)3 ) = 1.This means that the occupied volume by the large particles is the same as that by the small ones.From the definition of the volume ratio, the size ratio correspondingly becomes σ (1) /σ (2) = (N 2 /N 1 ) 1/3 = [(1 − ν 1 )/ν 1 ] 1/3 .Thus, as the size ratio increases, the number of collisions between large particles decreases.On the other hand, the impulse from the larger particle at each collision increases as compared with that from the smaller particle.
Figures 5(a ratios become (a) σ (1) /σ (2) = 10.0,(b) 6.93, (c) 4.63, and (d) 3.66, respectively.It should be noted that the number of particles of EDLSHS is fixed as N = 1000 in Fig. 5.It is quite apparent that the theory compares well with the simulation results in the wide range of the shear rate and without any fitting parameter when the size ratio is not large (or equivalently, ν 1 3.0 × 10 −3 in Fig. 5).On the other hand, some discrepancies between the theoretical prediction and the EDLSHS simulations are observed in the high shear regime when the size ratio becomes sufficiently large (see the data for ν 1 = 1.0 × 10 −3 in Fig. 5).
In particular, at γ * ≈ 30, the theory predicts a new DST-like transition in which the flow curve becomes an S-shape; in this region the temperature ratio versus the shear rate becomes a multivalued function (see Fig. 5(a)).Here, the upper branch becomes almost 100 times larger compared to the lower branch.This behavior is analogous to the ignited-quenched transitions for the shear-rate dependence of both the temperature and the viscosity for the monodisperse case [13,15].(See Appendix D for the minimum value of the size ratio at which this transition occurs.)The origin of the discrepancy between theory and simulations is essentially associated with the suppression of the collisions between the large (tracer) particles because the number of them becomes N 1 ∼ O(1) for σ (1) /σ (2) ≫ 1, as discussed in the following.
To verify our conjecture, we examine the simulation results obtained for different large system sizes: N = 10000 and 30000.We find that the disagreement between theory and simulation tends to decrease as the number of particles in the EDLSHS increases.As an illustration, Fig. 6 shows the dependence of both the temperature ratio ϑ and the viscosity ratio η 1 /η 2 on the number of particles N when we fix ν 1 = 1.0 × 10 −3 .Here, the relationships between the total number of particles and that of large particles correspond to (N, N 1 ) = (1000, 1), (10000, 10), and (30000, 30).As N 1 increases, the effect of collisions between large particles on rheology becomes non-negligible.The collisions between large 16/39  48) and ( 49) and the tracer limit explained in Appendix E, respectively.The dotted and dot-dashed lines represent the granular gas limits for Eqs. ( 48) and ( 49) and that under the tracer limit in Appendix E, respectively.
particles affect the flow curve in particular in the high shear regime.Correspondingly, the quantities discontinuously change at a certain shear rate; this shear rate depends on the number of particles.The above results suggest that (i) the discontinuous change predicted by the kinetic theory can be universally observed in the thermodynamic limit and (ii) the picture of an impurity enslaved to the host fluid (namely, when the tracer-tracer collisions are neglected) is insufficient to capture the above discontinuous transition.The fact that the seemingly natural "enslaved impurity" picture breaks down for large shear rates has been also shown the responsible for the extreme violation of energy equipartition in a sheared granular mixture in the tracer limit [61,62].We can also understand this finite size effect of EDLSHS when we observe the time evolution of the temperature ratio ϑ for a very large system.Figure 7 exhibits the typical evolution of ϑ for N = 30000.The solid line refers to the solution obtained for a binary mixture assisted by Eq. ( 45) with ν 1 = 1.0 × 10 −3 while the dashed line corresponds to the analytical result obtained in the Appendix E in the tracer limit (i.e., by neglecting collisions between tracer particles and by assuming that the excess species 2 is not affected by the presence of tracer particles).We observe the transient behavior in the result of EDLSHS from the tracer limit line (dashed line) to that of the (complete) solution including collisions between large particles (solid line).It is apparent that collisions between large tracer particles do not play any role in the early stage since the temperature ratio measured in the simulation agrees well with the tracer limit line (see the data for τ 2 in Fig. 7).As time goes on, however, those contributions become important for the rheology of the system.As a result, the temperature ratio measured in EDLSHS starts to increase abruptly (see the data for τ ≃ 3 in Fig. 7), and tends to converge to the asymptotic theoretical value (τ 5).
Let us check how the discontinuous changes of the temperature ratio and the viscosity ratio appear in the thermodynamic limit.According to Fig. 5, there must exist a critical 17/39  45) and ( 48) with ν 1 = 1.0 × 10 −3 while the dashed line corresponds to the analytical result obtained in the Appendix E in the tracer limit.

Velocity distribution function
In this subsection, let us compare the velocity distribution function (VDF) (21) of Grad's moment method with the one obtained by means of simulations.As a complement, we also include the exact VDF of a BGK-like kinetic model in the large shear limit (see the Appendix F).
For later analysis, let us introduce the dimensionless velocity c and the distribution function g i,G (c) as where f i,G stands for Grad's VDF (21) for species i.Now, we focus on the VDF of the larger particles 1.It is convenient to consider the marginal distribution g 1,G instead of using the full three-dimensional VDF.The distribution g 1,G is defined as where The VDF g (x,y) 1,G (c x , c y ) in Eq. ( 53) can characterize the anisotropy of the VDF induced by the shear flow.
Figures 9 and 10 present g 1,M (c x , c y ) for γ * = 0.32, 1.0, 3.2, 10, and 32.These values of the shear rate belong to the lower (0.32 and 1.0), intermediate (3.2 and 10), and higher (32) branches of the flow curve for e = 0.9 in Fig. 9 and e = 1 in Fig. 10, respectively.It is remarkable that Grad's distribution works well in the wide range of the shear rate as shown in Figs. 9 and 10.The corresponding trends are clearly observed when we consider the one-dimensional VDF in the x direction (see the Appendix G).Nevertheless, it seems that the enhancement of the VDF in the shear direction is underestimated in the theoretical prediction.It should be noted that this enhancement is suppressed for the one-dimensional VDF as shown in the Appendix G.
We also check whether the VDF obtained from the BGK-like model can be used (see the Appendix F for details).As expected, the deviation of the distribution of BGK-like model from that of the simulation is large for low shear regime.On the other hand, the agreement between BGK and simulations is reasonable in the high shear regime (see Fig. 11).In particular, the BGK distribution is more accurate than that of Grad's distribution in the  1,M (c x , c y ) in the (c x , c y )-plane for (a) γ * = 0.32, (b) 1.0, (c) 3.2, (d) 10, and (e) 32 when we fix ϕ = 0.01, e = 0.9, ξ env = 1.0, σ (1) /σ (2) = 2.0, and ν 1 = ν 2 = 1/2.The color plot corresponds to the simulation results.The solid, dashed, dotted, and dot-dashed lines represent the contours 0.2c max , 0.1c max , −0.1c max , and −0.2c max obtained from Grad's method ( 53) with (a) c max = 0.02, (b) 0.03, (c) 0.2, and (d, e) 0.1, respectively.small velocity region.In any case, it is important to recall that the BGK distribution obtained in the Appendix F only holds when T env = 0.This means that the possible discrepancies between the BGK distribution and simulations can be in part due to the fact that T env = 0 in the simulations.

Discussion and conclusion
In this paper, we have theoretically derived the rheology of a dilute binary mixture of inertial suspensions under USF.As in previous papers [15,16], two different but complementary approaches have been employed to solve the set of coupled Boltzmann kinetic equations.On the analytical side, Grad's moment method [58] has been used to approximately solve the Boltzmann equation.Since the mass and heat fluxes vanish in the USF, only the partial traceless stress tensors Π (i) αβ are retained in the trial distribution functions f i (V ).Then, the theoretical predictions for the temperature ratio T 1 /T 2 and the viscosity ratio η 1 /η 2 have been compared against computer simulations based on the event-driven Langevin simulation method.We have confirmed that the theoretical predictions agree with the results of simulation for hard spheres for various size ratios in most parameters' regions.We have found that the temperature ratio and viscosity ratio discontinuously change at a certain shear rate as 20/39 the size ratio increases.This feature cannot be captured by simulations when the size of the system is small.The above transition is similar to DST in dense suspensions or the first-order phase transition at equilibrium.Although the tracer limit of the theory is validated when the system size is small, the collisions between large tracer particles play dominant roles in the high shear regime.We have also compared the velocity distribution functions obtained by Grad's method and BGK-like model with those obtained from the simulations.
There are several future perspectives.First, we plan to analyze the mass transport of impurities in a sheared inertial suspension.As already did in Ref. [27], a Chapman-Enskoglike expansion around the local shear flow distribution obtained here will be considered to identify the shear-rate dependent diffusion D αβ , pressure diffusion D p,αβ , and thermal diffusion D T,αβ tensors.The determination of D αβ , D p,αβ , and D T,αβ will be discussed in a forthcoming paper.More importantly, the knowledge of the above diffusion tensors will allow us to analyze segregation by thermal diffusion [50].In the present paper, we have restricted to homogeneous systems, which makes the analysis easier than that for inhomogeneous systems.However, depending on the size or density of particles, the segregation is inevitable when one considers binary mixtures.In a sheared system, the segregation has been observed if there exists an inhomogeneous velocity profile [42].However, the velocity profile remains linear in 21/39 our simulations as far as we have checked.This linearity is violated if we consider systems under gravity or wall driven sheared systems.We believe that this scenario of segregation can be described by a dilute system described by the Boltzmann equation.We will analyze such systems in the near future.
Needless to say, we plan also to extend our analysis to moderately dense systems with the aid of the Enskog equation.The extension is tough but straightforward using a similar procedure as the one followed for monodisperse systems [16].This study will be also carried out in the future.and from Grant IB20079 funded by Junta de Extremadura (Spain) and by ERDF A way of making Europe.This work was initiated during a short stay of V.G. at the Yukawa Institute for Theoretical Physics (YITP, Kyoto University) supported by the YITP activity (YITP-T-18-03).V.G. appreciates the warm hospitality of the YITP at that time.S.T. and H.H. also acknowledge the warm hospitality of the Universidad de Extremadura during their stay there in 2020.
A. Difference between P (i) yy and P

(i) zz
In this Appendix, we show the difference between P (i) yy and P (i) zz .As mentioned in the main text, the second normal stress difference of species i, N zz is, in general, nonzero.However, the second difference 2 is treated as zero in the dilute limit of the kinetic theory.Figure A1 shows the plot of the ratio of the second to the first normal stress differences against the shear rate for σ (1) /σ (2) = 2.0 and 5.0 obtained from the simulations when we fix ϕ = 0.01, ξ env = 1.0, e = 0.9, and ν 1 = ν 2 = 1/2.Here, we have introduced the first normal stress difference of species i as 2 is values are much smaller than the values of 1 in the wide range of shear rates considered.Therefore, we do not consider the difference between them in this paper.It is noted that the second normal stress difference cannot be neglected when the volume fraction is finite.

B. Derivation of Λ (ij) αβ under the linear approximation of Grad's expansion
In this Appendix, we obtain the expression (24) for the collisional moment Λ (ij) αβ .For this purpose, we introduce the dimensionless velocities and equivalently ) in terms of G and g.Using the Grad's trial distribution (21), we can rewrite We note that nonlinear contributions of the stress tensor Π αβ are ignored in this Appendix.Let us rewrite the argument of the exponential part in Eq. (B3) as Introducing G ′ as one gets the identities Thus, we can rewrite Eq. (B3) as Let us rewrite Eq. ( 16).From Eqs. ( 1) and (B5), one gets Using Eqs.(B3) and (B8), we can rewrite the collisional moment Λ with the linear collisional moment For further calculation, let us first introduce with We also introduce with Using Eqs.(B11) and (B13), Eq. (B10) is rewritten as Let us go further by integrating over σ.Here, the following results are needed: with With the aid of Eqs.(B16)-(B17), one gets and there, αβ .(B20) Similarly, one achieves the result and then, In addition, one gets and then Substituting Eqs.(B20), (B22), and (B24) into Eq.(B15), one obtains Finally, the combination of Eqs.(B9) and (B25) yields Eq. ( 24).
When we focus on the reduced temperature θ (see Fig. C1), the effect of the bidispersity only appears around an intermediate shear regime ( γ * ≃ 5.0), where the discontinuous change corresponding to the DST is observed.Although this discontinuous change itself is reported even in monodisperse systems [13], the point at which the discontinuous change occurs depends on the size ratio.It is noteworthy that the change of the reduced temperature is drastic but continuous when the size ratio becomes large (see the data for σ (1) /σ (2) = 2.0 and 5.0 in Fig. C1).

27/39
As well as in Fig. C1, the viscosity η * is also plotted against the shear rate γ * in Fig. C2.If the size ratio σ (1) /σ (2) is close to unity, such as 1.4, the flow curves of θ and η * are similar to the corresponding ones for monodisperse gases, in which there are discontinuous changes of θ and η * around γ * ≈ 5.However, as the size ratio increases, the discontinuous changes of θ and η * become continuous.Moreover, these flow curves for inelastic inertial suspensions for large σ (1) /σ (2) are characteristic.Indeed, the slopes of θ and η * are oscillated with γ * before reaching their asymptotic values in the large shear rate limit.We also draw 3D-phase diagrams of the number of solutions obtained by the kinetic theory in the (ν 1 , γ * , e)-plane for ϕ = 0.01 in Fig. C3.The filled regions represent those whose number of solutions is three, while the empty regions represent only one solution.These plots show that the regions for the multiple solutions are localized in the narrow regimes in the (ν 1 , γ * , e)-plane.

D. Appearance/disappearance of the discontinuous transition
In this Appendix, let us show how the discontinuous transition appears/disappears when we change the parameters of the mixture.This appendix consists of three subsections.In the first part, we discuss how the results depend on the environmental temperature ξ env .In the second part, we distinguish the region of DST-like behavior from the CST-like behavior when we fix ν 1 = ν 2 = 1/2.In the last part, we also distinguish the region of DST-like behavior from the CST-like behavior if we fix the volume ratio V = 1, D.1.Effect of the environmental temperature ξ env

Fig. C2
The dimensionless viscosity η * against the dimensionless shear rate γ * for σ (1) /σ (2)   regime, but is independent in the high shear regime.The latter fact is understood because interparticle collisions are dominant in the latter regime.Figure D1 illustrates the above fact: The high shear regime is independent of the choice of the environmental temperature, but the low shear regime is determined by the value of the environmental temperature.It is interesting to note that the Newtonian regime becomes narrower as ξ env increases.More importantly, DST-like behavior for η * for low ξ env becomes CST-like as ξ env increases.

D.2. Effect of the size ratio for
Next, let us consider the size ratio dependence in the case of ν 1 = ν 2 = 1/2 based on the theoretical calculation.In this case, the discontinuous jumps are observed when the size ratio 29/39 no DST-like trans.
Fig. D2 Plot of the critical size ratio against the restitution coefficient when we fix ϕ = 0.01, ξ env = 1.0, and is not large such as σ (1) /σ (2) = 1.4 in Figs. 2 and 3 as shown in Fig. D2.On the other hand, the flow curves become continuous for larger size ratio.We can understand this behavior by considering first the discontinuous jump for the monodisperse system (ν 1 = 1, ν 2 = 0).Depending on the value of the (reduced) shear rate γ * ≡ γ/ζ 1 , there are two different regimes; high-shear and low-shear regimes.The former regime is known as Bagnold's expression, η * ∝ γ * /(ξ 2 env ϕ 2 ) for e < 1 [13].We note that, for the elastic case, a different expression is obtained as η * ∝ γ * 2 .However, the latter regime (low shear regime) is determined by the interaction between the particles and the solvent [13], and so η * ∼ 1.These two regimes switch to each other at γ * ≃ 1.Given that the difference between two regimes is proportional to the inverse of the volume fraction, the flow curve forms an S-shape connecting the two regimes.

30/39
Now, we consider binary systems.If the size ratio is not sufficiently large, such as σ (1) /σ (2) = 1.4 as shown in Figs.C1 and C2, the picture for the monodisperse system can also be used for a binary system.This means that the discontinuous jumps appear in this case.On the other hand, as the size ratio increases, collisions between smaller and larger particles compete with those between particles with the same size.This means that we need to discuss the mixing energy between smaller and larger particles in this case.Relating to this, we may use a discussion analogous to the phase coexistence and spinodal lines at equilibrium phase transitions, respectively, in the phase space of (θ, γ * , σ (1) /σ (2) ). Figure

D.3. Effect of the size ratio for
no DST-like trans.Let us consider the case of constant volume ratio V = 1.As shown in Fig. 6, the discontinuous transition occurs as the fraction ν 1 decreases, i.e., the size ratio increases.This transition 31/39 is different from the one found in Appendix D.2. Figure D4 plots the critical line between the discontinuous transition and continuous transition for ϕ = 0.01, ξ env = 1.0, and V = 1 based on the theoretical calculation.As the restitution coefficient e increases, the minimum size ratio also increases, which means that the fraction ν 1 decreases as ν 1 = 1/[1 + (σ (1) /σ (2) ) 3 ] from the definition of the volume ratio V. Unfortunately, it is a tough job to check this behavior in simulations.When the DST-like transition occurs, one needs to simulate a situation where multiple collisions between large particles occur.However, as the size ratio increases, the fraction of the larger particles, ν 1 , becomes small, then the collision frequency between them also decreases.This means that the time for multiple collisions exceeds the limit of realistic simulation time.
E. Detailed analysis in the tracer limit and the finite size effect of the simulation results In this Appendix, we display the explicit expressions of the partial pressure tensors of a binary mixture in the tracer limit.These expression are then compared with the simulation results when the number of particles is small.In the tracer limit (ν 1 → 0), the kinetic equation for the velocity distribution function f 2 of the excess granular gas 2 is the (closed) nonlinear Boltzmann equation since its state is not perturbed by the presence of the tracer particles 1.This means that collisions between tracer and gas particles in the kinetic equation for P  15) for i = 2.In addition, since the concentration of tracer particles is negligible, one can also neglect the tracer-tracer collisions in the kinetic equation for P (1) αβ .This implies that Λ (11) αβ + Λ (12) αβ → Λ (12) αβ in Eq. ( 15) for i = 1.The expressions of the (reduced) elements of the pressure tensor Π (2) αβ coincide with those obtained for a monodisperse granular suspension.The nontrivial components of Π (2) αβ are given by [13] Π (2) where γ = γ/ζ 2 and we have introduced with the partial volume fraction Here, it should be noted that the global temperature is approximately given by θ ≃ θ 2 in the tracer limit [26].Using the same procedure as in Ref. [13], the reduced shear rate γ * is Now, let us calculate the quantities for the tracer species 1.First, the quantities Λ (12) αα , Λ xy , and Λ ′ (12) xy are written as Λ (12)  αα ≡ 1 m * 3/2 2 (ϑ ′ + 1) where we have introduced ϑ ′ ≡ m 2 θ 1 /(m 1 θ 2 ).Then, the nonzero elements of Π (1) Π (1)  xx = −2Π (1)  yy , Π (1)  yy = Π Figure E1 presents the shear-rate dependence of both the temperature ratio ϑ and the viscosity ratio η 1 /η 2 in the tracer limit.It should be noted that the flow curves become smooth in the whole range of the shear rate even for a larger size ratio.The limitation of the tracer limit is also understood in Fig. E2, where the absolute values of the ratio Λ αα is given by Eq. ( 24).In the low shear regime, the values for (i, j) = (1, 1) and (2, 1) are smaller than unity, which means that the contributions coming from the collisions between the large tracer particles are negligible.This indicates that the tracer limit description is a reasonable approximation in this regime.In the high shear regime, on the other hand, the contributions from the collisions between tracer particles play an important role to the flow curve, though the number of collisions is small.Moreover, it is interesting that Λ (2,1) αα and Λ (1,2) αα become negative in the high and intermediate shear regimes, respectively, though their origins are not clear.As the number of particles used in the simulation increases, the simulation results recover the values of ϑ and η 1 /η 2 in the high shear regime.Then, we expect that the results of simulation for N → ∞ agree with the theoretical results.In other words, the results of EDLSHS containing a small number of particles is not reliable in this regime.for (i, j) = (1, 1) (solid line), (1, 2) (dashed line), and (2, 2) (dotted line) against the dimensionless shear rate γ * for ϕ = 0.01, ξ env = 1.0, e = 0.9, V = 1, and ν 1 = 1.0 × 10 −3 .

F. Two-dimensional velocity distribution function of BGK model
A possible way of overcoming the mathematical difficulties associated with the Boltzmann collision operators J ij [f i , f j ] is to use a kinetic model.As usual, the idea behind a kinetic model is to replace the true operator J ij by a simpler term that retains the main physical properties of the above operator.In the case of dilute granular mixtures, a BGK-like kinetic model was proposed in Ref. [63].In the case of the USF state (where U 1 = U 2 = U ), the BGK-like model is obtained by the replacement of the Boltzmann collision operator J ij [f i , f j ] by the diffusive term 34/39 where we have introduced the quantities , (F2) T ij = T i + 2m i m j (m i + m j ) 2 (T j − T i ). (F5) The corresponding BGK-like equation for the distribution f 1 in the steady USF is .
(F6) The kinetic equation for f 2 is obtained from Eq. (F6) by setting 1 ↔ 2. So far, we have not been able to obtain the explicit exact form of f i (V ) in Eq. ( F6).An exception corresponds to the simple limit case T env = 0 with keeping ζ i = const.It corresponds to a situation where the background temperature T env is much smaller than the kinetic temperature T under the high shear rate limit.Hence, the suspension model ignores the effects of thermal fluctuations on solid particles and the impact of the gas phase on grains is only accounted by the drag force term.Although ζ i should be proportional to √ T env for hard-core molecules, such a simplified model has been employed in some previous works [9,10,14].

Fig.
Fig. D4 Plot of the critical size ratio against the restitution coefficient e when we fix ϕ = 0.01, ξ env = 1.0, and V = 1.

Fig
Fig. E1 (a) Temperature ratio ϑ and (b) viscosity ratio η 1 /η 2 against the dimensionless shear rate γ * in the tracer limit for the same set of parameters of Fig. 5.The data of the simulation are obtained for N = 1000.