Rheology of Dense Sheared Granular Liquids
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec Rheology of Dense Sheared Granular Liquids
Koshiro Suzuki ∗ and Hisao Hayakawa † ∗ Canon Inc., 30-2 Shimomaruko 3-chome, Ohta-ku, Tokyo 146-8501, Japan † Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwake-cho, Kyoto 606-8502, Japan
Abstract.
The rheology of dense sheared granular liquids is investigated based on the mode-coupling theory (MCT). Thisextended MCT includes correlations for the density-current mode as well as the density-density correlation mode, and a self-consistent coupling equation for the energy balance condition. The extended MCT exhibits disappearance of the two-steprelaxation of the density-density correlation function, and also successfully reproduces the density dependence of the shearviscosity for volume fractions between 0.50 and 0.60, if we shift the density. However, it predicts unphysical tendency for thegranular temperature. The cause of this drawback and the possibilities of its amendment are discussed.
Keywords:
Granular matter, Liquid theory, Mode coupling theory
PACS:
INTRODUCTION
Establishing a macroscopic description of granular materials has been a long-term challenge for both science andengineering. The problem extends to a vast range, from creep motion or force chain dynamics of frictional particles,two-phase flow of air-fluidized beds, to nonequilibrium transport of sheared granular flows [1]. Similar to solid-liquidtransitions, critical features of the jamming transition in the vicinity of the random close packing and their relation tothe glass transition have attracted much interest in the last decade [2, 3]. Even when we focus only on the classicalproblem of the flow properties well below the jamming transition density j J , which can be traced back to Bagnold’swork [4], we have not yet understood the rheological properties of dense granular flows. One of the remarkableachievements is the extension of the Boltzmann-Enskog kinetic theory to inelastic hard disks and spheres [5], whichstimulated the following works [6, 7]. However, it has been recognized that the kinetic theory breaks down at highdensities with volume fraction j > . . < j < j J .On the other hand, a continuum description with long-time correlations has been constructed for thermal glassy liq-uids. The mode-coupling theory (MCT) exhibits the two-step relaxation of density correlation functions characteristicof glasses [10]. Although MCT has presented remarkable success, it is marred with problems; for instance, it predictsa non-ergodic transition which is not observed in experiments [11], or its critical density j MCT = .
516 [12] is farbelow the glass transition density, which is observed to be j g ≃ . − .
60 in numerical simulations [13].Despite these problems, it is still tempting to extend MCT to granular materials, because it might be the simplestmethod to include correlations. In MCT, these correlations are encoded in a memory kernel, which is absent in thekinetic theory. The extension of MCT to randomly driven granular systems has been proposed and analyzed [14].However, the physics of sheared flows is completely distinct from random driving cases, and its study has to beaddressed independently. Indeed, the characteristic plateau in the density correlation function for both glassy andrandomly driven granular systems does not exist in sheared granular liquids [15]. The extension of MCT to shearedgranular liquids can be traced back to Ref. [16], and an explicit calculation of time correlation functions has beenreported in Ref. [17] without considering the energy balance condition. However, the problem of rheology and thecharacterization of the steady state have not been addressed so far. In this article, we demonstrate that the extendedMCT reproduces the results of molecular dynamics (MD) simulations for the relaxation of time correlation functionsand the density dependence of the shear viscosity.The organization of this paper is as follows. We first explain the microscopic set up of the theory, in particular theLiouville equation and the steady-state condition. Then we derive the mode-coupling equations, without specifyingthe initial distribution function. We perform a concrete calculation for the case of a canonical initial distribution andevaluate the validity of the results. Finally, we discuss the problem of the formulation and the possibilities of itsamendment.
ICROSCOPIC SETUPEquations of Motion
The starting point of our theory is an assembly of N identical soft, smooth, inelastic spheres with mass m anddiameter d , contained in a box of volume V . Thus, the number density is given by n = N / V . The interaction betweenthe spheres are elastic as well as inelastic (dissipative), both of which are repulsive and emerge at contact. The systemis subjected to a bulk shearing, which is uniform and constant with shear rate ˙ g . At first the system is equilibratedwithout shear and dissipation. Then, at t = t ( < ) , both of them are switched on, which eventually leads the systemto a nonequilibrium steady state. The Newtonian equation of motion for the i th sphere ( i = , · · · , N ) is given by thefollowing set of Sllod equations [18]: ˙ rrr i ( t ) = ppp i ( t ) m + ˙ ggg · rrr i ( t ) , (1)˙ ppp i ( t ) = FFF ( el ) i ( t ) + FFF ( vis ) i ( t ) − ˙ ggg · ppp i ( t ) . (2)Here, GGG ( t ) ≡ { rrr i ( t ) , ppp i ( t ) } Ni = is a set of the positions and the momenta of all the grains at time t , ˙ ggg is the shearrate tensor whose components are assumed to be given by ˙ g mn = ˙ gd m x d n y , and FFF ( el ) i ( t ) , FFF ( vis ) i ( t ) are the elasticand dissipative interactions, respectively. The Greek indices m , n , l , · · · denote spatial components { x , y , z } , and therule for the summation over repeated indices is implied. Here, the elastic force is given by FFF ( el ) i ( t ) ≡ − (cid:229) j = i Q ( d − r i j ) ¶ u ( r i j ( t )) / ¶ rrr i j ( t ) , where rrr i j ( t ) ≡ rrr i ( t ) − rrr j ( t ) and r i j ( t ) ≡ | rrr i j ( t ) | are, respectively, the relative position and thedistance between the i th and j th spheres, u ( r ) is the two-body potential, and Q ( x ) is a step function which is 1 for x > FFF ( vis ) i ( t ) ≡ − z (cid:229) j = i Q ( d − r i j ) ˆ rrr i j ( ˙ rrr i j · ˆ rrr i j ) , where ˆ rrr ≡ rrr / | rrr | is a unit vector and z is aviscous constant corresponding to the harmonic potential. Liouville Equation
To formulate a theory, we rewrite the equations of motion, Eqs. (1) and (2), to the form of the Liouville equation[18]. In this formulation, the equation of motion for an arbitrary phase-space variable A ( GGG ) casts into the form ddt A ( GGG ( t )) = i L ( GGG ) A ( GGG ( t )) , (3)where the operator i L ( GGG ) is referred to as the Liouvillian. The explicit form of i L ( GGG ) for Eqs. (1) and (2) is given by i L ( GGG ) = i L ( el ) ( GGG ) + i L ˙ g ( GGG ) + i L ( vis ) ( GGG ) , (4)where the elastic part i L ( el ) ( GGG ) , the shear part i L ˙ g ( GGG ) , and the viscous part i L ( vis ) ( GGG ) are given by i L ( el ) ( GGG ) = (cid:229) Ni = h ppp i m · ¶¶ rrr i + FFF ( el ) i · ¶¶ ppp i i , i L ˙ g ( GGG ) = (cid:229) Ni = h ( ˙ ggg · rrr i ) · ¶¶ rrr i − ( ˙ ggg · ppp i ) · ¶¶ ppp i i , i L ( vis ) ( GGG ) = (cid:229) Ni = FFF ( vis ) i · ¶¶ ppp i . On the otherhand, the equation for the phase-space distribution function r ( GGG , t ) reads ¶¶ t r ( GGG , t ) = − i L † ( GGG ) , (5)where the adjoint Liouvillian i L † ( GGG ) is defined as i L † ( GGG ) = i L ( GGG ) + L ( GGG ) , together with the phase-space contrac-tion factor, L ( GGG ) = ¶ / ¶ GGG · ˙ GGG = − ( z / m ) (cid:229) h i , j i Q ( d − r i j ) < . Here, we denote (cid:229) Ni = (cid:229) j = i · · · as (cid:229) h i , j i · · · . The formalsolution of Eqs. (3) and (5) are given by A ( GGG ( t )) = e i L t A ( GGG ( )) , r ( GGG , t ) = e − i L † t r ( GGG , ) . (6)To proceed, it is convenient to rewrite r ( GGG , t ) in Eq. (6), by use of an identity [18], as r ( GGG , t ) = r ini ( GGG ) + Z t ds e − i L † s (cid:0) − i L † (cid:1) r ini ( GGG ) , (7)where r ini ( GGG ) ≡ r ( GGG , ) is the initial distribution function. Then, we can introduce the nonequilibrium work function, W ( GGG ) , as the eigenvalue of i L † : L † r ini ( GGG ) = − r ini ( GGG ) W ( GGG ) , (8)from which we obtain r ( GGG , t ) = r ini ( GGG ) + Z t ds e − i L † s [ r ini ( GGG ) W ( GGG )] . (9)Note that W ( GGG ) includes the nonequilibrium information of the dynamics of the system, which is caused by shearingand dissipation. In the remainder, we denote the ensemble average with respect to r ini ( GGG ) by h· · ·i = Z d GGG r ini ( GGG ) · · · . (10)The nonequilibrium ensemble average of a phase-space variable A ( GGG ( t )) , which we denote h A ( t ) i , can be expressedas h A ( t ) i = Z d GGG r ini ( GGG ) A ( GGG ( t )) = Z d GGG r ( GGG , t ) A ( GGG ( )) . (11)The equivalence of the two expressions in Eq. (11) can be shown by the adjoint relation of the Liouvillians [18].From Eqs. (9), (11), and the adjoint relation, we obtain h A ( t ) i = h A ( ) i + R t ds h A ( s ) W ( ) i , which is referred to as thegeneralized Green-Kubo formula [18]. The steady-state average, A ¥ , is obtained from this formula as A ¥ = h A ( ) i + Z ¥ ds h A ( s ) W ( ) i . (12) Steady-State Condition
The steady state of a sheared granular flow is characterized by the balance of energy. The internal energy H ( GGG ) isgiven H ( GGG ) = (cid:229) Ni = (cid:2) ppp i / m + (cid:229) j = i u ( r i j ) (cid:3) . The following equality can be shown from Eqs. (1) and (2),˙ H ( GGG ) = − ˙ gs xy ( GGG ) − R ( GGG ) , (13)where s xy ( GGG ) = V (cid:229) Ni = (cid:2) p xi p yi / m + y i F xi (cid:3) is the shear stress and R ( GGG ) = − (cid:229) h i , j i ˙ rrr i j · FFF ( vis ) i j = z (cid:229) h i , j i Q ( d − r i j ) ( ˙ rrr i j · ˆ rrr i j ) is Rayleigh’s dissipation function. From Eq. (13) the steady state is characterized by the condition,˙ g ( s xy ) ¥ + R ¥ = , (14)where the steady-state values ( s xy ) ¥ and R ¥ can be evaluated by Eq. (12). Equation (14) implies the balance betweenthe heating due to shear and the cooling due to dissipation. As will be discussed later, the time evolution of thetemperature is not included in Eq. (14) in our framework, and hence it will be imposed as a condition which determinesthe steady-state temperature. SHEARED GRANULAR MCT
So far the formulation is exact but not solvable. To obtain a calculable theory for the low-frequency modes, Moriequations are approximated by MCT to derive a closed set of equations. MCT is also applied to the generalized Green-Kubo formula to calculate the shear stress and the energy dissipation rate. First we briefly review the essence of thisprocedure [17].
MCT Equations
We identify the density and the current density fluctuations, which are n qqq ≡ (cid:229) Ni = e iqqq · rrr i − N d qqq , and j l qqq ≡ (cid:229) Ni = (cid:0) p l i / m (cid:1) e iqqq · rrr i in Fourier space, as the slow modes. In standard granular hydrodynamics, the tempera-ture fluctuation should also be considered on the same grounds as n qqq and j l qqq . However, we attempt to formulate atheory without the temperature fluctuation, since otherwise it will be too complicated and almost untractable. How tocompensate the dynamics of the temperature will be discussed below.he Mori equations for n qqq and j l qqq are obtained by applying the projection operator P ( t ) X = (cid:229) kkk h Xn ∗ kkk ( t ) i NS k ( t ) n kkk ( t ) + (cid:229) kkk h X j m ∗ kkk ( t ) i Nv T j m kkk ( t ) , (15)where X ( GGG ) is an arbitrary phase-space variable. Here, S k = h n kkk n ∗ kkk i / N is the static structure factor, v T = h j l kkk j l ∗ kkk i / ( N ) is the equal-time correlation of the current, and kkk ( t ) ≡ ( k x , k y − ˙ g tk x , k z ) is the wave vector in the sheared frame. Notethat h· · ·i denotes the ensemble average with respect to the initial distribution function, Eq. (10). The orthogonalityof n kkk and j l kkk , i.e. h n kkk j l ∗ kkk i =
0, is necessary for the idempotency, P ( t ) = P ( t ) . We assume that r ini ( GGG ) is even withrespect to the momentum, which is sufficient for h n kkk j l ∗ kkk i =
0. Then, to obtain a closure for these equations, the secondprojection operator P mc ( t ) = P nn ( t ) + P n j ( t ) , (16)where P nn ( t ) = (cid:229) kkk > ppp h Xn ∗ kkk ( t ) n ∗ ppp ( t ) i N S k ( t ) S p ( t ) n kkk ( t ) n ppp ( t ) , P n j ( t ) = (cid:229) kkk , ppp h Xn ∗ kkk ( t ) j m ∗ ppp ( t ) i N S k ( t ) v T n kkk ( t ) j m ppp ( t ) , (17)is applied to the memory kernels, together with the factorization approximation. Then we obtain a set of coupledequations for four time correlation functions, F qqq ( t ) ≡ (cid:10) n qqq ( t ) ( t ) n qqq ( ) ∗ (cid:11) / N , H l qqq ( t ) ≡ i D j l qqq ( t ) ( t ) n qqq ( ) ∗ E / N , ¯ H l qqq ( t ) ≡ i (cid:10) n qqq ( t ) ( t ) j l qqq ( ) ∗ (cid:11) / N , and C mn qqq ( t ) ≡ D j m qqq ( t ) ( t ) j n qqq ( ) ∗ E / N . However, we might expect that anisotropic effects are sub-dominant, considering that sheared colloidal systems are almost isotropic [19, 20]. Hence, we resort to the isotropicapproximation , and reduce the time correlation functions to two scalar functions, F q ( t ) and Y q ( t ) , which depend onlyon the modulus of the wave vector, q ≡ | qqq | . Here, Y q ( t ) is introduced by the definition ¯ H l qqq ( t ) ≡ − q l ( t ) Y q ( t ) , andhence is essentially the density-current correlation. These two functions are related to the remaining two by the rela-tions H l qqq ( t ) ≃ [ q l ( t ) / q ( t ) ] ddt F q ( t ) and C mn qqq ( t ) ≃ d mn (cid:8) ddt Y q ( t ) + (cid:2) Y q ( t ) / q ( t ) (cid:3) ddt q ( t ) (cid:9) . It should be noted that,while H l qqq ( t ) is obtained by the time derivative of F qqq ( t ) , ¯ H qqq ( t ) cannot be derived from F qqq ( t ) ; i.e., the minimal set oftime correlation functions is { F q , Y q } . Then, the MCT equations for F q ( t ) and Y q ( t ) read d dt F q ( t ) = − Z q ( t ) F q ( t ) − A q ( t ) ddt F q ( t ) − Z t ds M qqq ( s ) ( t − s ) dds F q ( s ) , (18) d dt Y q ( t ) = − Z q ( t ) Y q ( t ) − A ll q ( t ) ddt Y q ( t ) − Z t ds M ll qqq ( s ) ( t − s ) dds Y q ( s ) , (19)where the memory kernels M qqq ( t ) and M ll qqq ( t ) are quadratic forms in terms of the time correlation functions. Thecoefficients Z q , A q , A ll q and the coefficients of the quadratic forms in the memory kernels are equal-time correlations,which originate from the projection operators. Hence, their explicit forms depend on the choice of r ini ( GGG ) , which wewill discuss later. We simply note here that the memory kernels consist of elastic and dissipative terms. For instance, M qqq ( t ) has the following form, M qqq ( t ) = M ( el ) qqq ( t ) + ( z / m ˙ g ) M ( vis ) qqq ( t ) , where M ( el ) qqq ( t ) is the conventional term whichinduces the glass transition for sheared non-dissipative particles [21], and M ( vis ) qqq ( t ) represents the dissipative terms.These features also hold for the kernel M ll qqq ( t ) . It is significant that the dissipative kernels have negative contributionsand are possible to suppress the glassy elastic term. The intuitive picture of this effect will be discussed below.Accordingly, the second projection operator P mc ( t ) is applied to the time correlation h A ( s ) W ( ) i in Eq. (12),together with the factorization approximation. Then we obtain a quadratic form of the time correlation functions.Again, the explicit forms of the steady-state formulas for s xy and R depend on r ini ( GGG ) , but exhibit the following formsin general, ( s xy ) SS = ( s ( el ) xy ) SS + z m ˙ g ( s ( vis ) xy ) SS , R SS = z m ˙ g h R ( loc ) SS + D R SS i . (20)Here, ( s ( el ) xy ) SS and z ( s ( vis ) xy ) SS / ( m ˙ g ) are the elastic and the dissipative components of the shear stress, and R ( loc ) SS and D R SS are the local and the non-local contributions, which correspond to the first and the second terms of Eq. (12),respectively. anonical Initial Distribution To proceed to explicit calculations, it is necessary to specify the initial distribution function, r ini ( GGG ) . We considerthe following form, r ini ( GGG ) = e − I ( GGG ) / R d GGG e − I ( GGG ) , where I ( GGG ) is the effective potential. Then, from Eq. (8), the workfunction W ( GGG ) is expressed as W ( GGG ) = i L ( GGG ) I ( GGG ) − L ( GGG ) = ˙ I ( GGG ) − L ( GGG ) . (21)The most naive choice of I ( GGG ) is I ( GGG ) = H ( GGG ) T , (22)which corresponds to the canonical initial distribution. In this case, the system is equilibrated at temperature T at t < t =
0, and shear and dissipation are switched on at t = t =
0. In this section, we show the results for the case ofcanonical initial distribution and discuss its validity.
Interpretation and Determination of the Temperature
A simple feature of r ini ( GGG ) under I ( GGG ) of Eq. (22) is that the momentum integral in Eq. (10) is reduced to aGaussian which can be performed. Then, the canonical temperature T in Eq. (22) appears in the coefficients of theMCT equations and the memory kernels. This poses a problem, since the dynamics governed by the MCT equations,and hence the steady-state averages, apparently depend on T , while it is expected to be independent from physicalgrounds. The origin of this problem resides in our treatment where the dynamics is projected onto the density and thecurrent density fluctuations, but not on the temperature fluctuation, and hence the time evolution of the temperature isnot included.To remedy this problem, we identify T with the steady-state temperature, and determine it from the energy balancecondition, Eq. (14). As shown below, T appears in the coefficients of the time correlation functions via the projectionoperators; that is, it appears in MCT equations, Eqs. (18) and (19), as well as in the energy balance condition, Eq. (14).We solve MCT equations iteratively until the temperature T is self-consistently determined by the energy balancecondition. The initial conditions for Eqs. (18) and (19) are F q ( ) = S q , ddt F q ( ) = Y q ( ) = ddt Y q ( ) = v T . Explicit Expressions
For the case of the canonical initial distribution, we can explicitly calculate the coefficients Z q , A q , and A ll q inEqs. (18), (19), and the coefficients which appear in the memory kernels M qqq ( t ) , M ll qqq ( t ) as well as in the steady-stateformulas for s xy and R , i.e. Eqs. (20). Coefficients–
The coefficients Z q , A q , and A ll q are given by Z q = v T q / S q , A q = p ( z H / m ) [ − j ( qd ) + j ( qd )] , and A ll q = p ( z H / m ) [ − j ( qd )] , where v T = p T / m is the thermal velocity, j l ( qd ) is the spherical Bessel functionof the l -th order, and z H m = √ p (cid:0) − e (cid:1) g ( d ) nd r Tm (23)is the effective frequency of the dissipation in the hard-core limit. Here, e is the normal restitution coefficient, and g ( d ) is the contact value of the equilibrium radial distribution function. We take this limit because we address the rheologybelow the jamming transition density. Memory Kernels–
The memory kernel for F q ( t ) , i.e. M qqq ( t ) , is given up to linear order in z H / ( m ˙ g ) by M qqq ( t ) = M ( el ) qqq ( t ) + ( z H / m ˙ g ) (cid:229) i = M ( vis i ) qqq ( t ) , where M ( el ) qqq ( t ) , M ( vis1 ) qqq ( t ) , and M ( vis2 ) qqq ( t ) are given by M ( el ) qqq ( t ) = nv T p q Z ¥ dk k Z q + k | q − k | d p pV ( el ) M ( ¯ k ( t ) , ¯ p ( t )) V ( el ) M ( k , p ) F k ( t ) F p ( t ) , (24) M ( vis1 ) qqq ( t ) = ˙ g p mq Z ¥ dk k Z q + k | q − k | d pV ( vis ) M ( ¯ k ( t ) , ¯ p ( t ) , a ( t ) , b ( t )) V ( el ) M ( k , p ) F k ( t ) ˙ F p ( t ) , (25) M ( vis2 ) qqq ( t ) = − ˙ g p mq Z ¥ dk k Z q + k | q − k | d p p V ( el ) M ( ¯ k ( t ) , ¯ p ( t )) V ( vis ) M ( k , p , a , b ) F k ( t ) Y p ( t ) . (26)ere, ¯ k ( t ) = kh ( ˙ g t ) with h ( x ) = p + x / V ( el ) M ( k , p ) = c k (cid:0) q + k − p (cid:1) + c p (cid:0) q + p − k (cid:1) , (27) V ( vis ) M ( k , p , a , b ) = sin a sin ( a + b ) · mI ( ) ( k ) (cid:2) cos a cos ( a + b ) − sin a sin ( a + b ) (cid:3) · mI ( ) ( k ) − cos b · mI ( ) ( p ) , (28)where I ( ) ( q ) = [ − j ( qd )] and I ( ) ( q ) = [ − j ( qd ) + j ( qd )] , and the angles a , b are defined by cos a =( q + k − p ) / ( qk ) , cos b = ( q + p − k ) / ( qp ) . Note that V ( el ) M ( k , p ) is the vertex function of the equilibrium MCT[10]. Note also that Y q ( t ) couples to the MCT equation of F q ( t ) , Eq. (18), only via M ( vis2 ) qqq ( t ) , so Y q ( t ) decouplesfrom F q ( t ) in the limit z →
0. Hence, our theory reduces to the sheared MCT of thermal glassy systems in the limit z →
0, together with appropriate thermostats, e.g. Gaussian Isokinetic thermostat, which controls the temperature.Similarly, the memory kernel for Y q ( t ) , i.e. M ll qqq ( t ) , is given up to linear order in z H / ( m ˙ g ) by M ll qqq ( t ) = M ( el ) ll qqq ( t )+( z H / m ˙ g ) (cid:229) i = M ( vis i ) ll qqq ( t ) , where M ( el ) ll qqq ( t ) , M ( vis1 ) ll qqq ( t ) , and M ( vis2 ) ll qqq ( t ) are given by M ( el ) ll qqq ( t ) = nv T p q Z ¥ dk k Z q + k | q − k | d p pW ( el ) ( ¯ k ( t ) , ¯ p ( t ) , k , p , a , b ) F k ( t ) F p ( t ) , (29) M ( vis1 ) ll qqq ( t ) = ˙ g p mq Z ¥ dk k Z q + k | q − k | d p h kc k W ( vis1 ) ( ¯ k ( t ) , ¯ p ( t ) , a ( t ) , b ( t )) + pc p W ( vis2 ) ( ¯ k ( t ) , ¯ p ( t ) , a ( t ) , b ( t )) i F k ( t ) ˙ F p ( t ) , (30) M ( vis2 ) ll qqq ( t ) = − ˙ g p mq Z ¥ dk k Z q + k | q − k | d p p (cid:26) kc ¯ k ( t ) W ( vis1 ) ( k , p , a , b ) + pc ¯ p ( t ) (cid:20) W ( vis2 ) ( k , p , a , b ) + ( ˙ g t ) W ( vis3 ) ( k , a , b ) (cid:21)(cid:27) × F k ( t ) Y p ( t ) . (31)The detailed expressions for the coefficients W ( el ) , W ( vis1 ) , W ( vis2 ) , and W ( vis3 ) are omitted here. Note that M ll qqq ( t ) isthe trace of the tensor M mn qqq ( t ) , and hence its coefficients cannot be expressed as a product of two vertex functions, asis the case for M qqq ( t ) . Work Function and the Steady-State Formula–
The explicit form of the work function can also be calculated fromEqs. (21) and (22) as W ( GGG ) = − b ˙ g V s xy ( GGG ) − b R ( GGG ) − L ( GGG ) , where b = / T is the inverse temperature. Thisexpression can be identically rewritten as W ( GGG ) = − b ˙ g V s ( el ) xy ( GGG ) + b ˙ g V s ( vis )( ) xy ( GGG ) − b D R ( GGG ) , (32)where s ( el ) xy ( GGG ) , s ( vis )( ) xy ( GGG ) , and D R ( GGG ) are given by s ( el ) xy ( GGG ) = V (cid:229) Ni = h p xi p yi / m + y i F ( el ) xi i , s ( vis )( ) xy ( GGG ) = − z H V (cid:229) h i , j i (cid:16) ppp ij m · ˆ rrr i j (cid:17) ˆ x i j ˆ y i j r i j Q ( d − r i j ) , and D R ( GGG ) = R ( ) ( GGG ) + T L ( GGG ) , (33)with R ( ) ( GGG ) = z H (cid:229) h i , j i (cid:16) ppp ij m · ˆ rrr i j (cid:17) Q ( d − r i j ) . Note that D R ( GGG ) is the purely dissipative contribution which exists inthe absense of shear, among which R ( ) ( GGG ) is the contribution from the energy dissipation due to inelastic collisions.The steady-state shear stress reads ( s xy ) SS = ( s ( el ) xy ) SS + ( z H / m ˙ g ) (cid:229) i = ( s ( vis i ) xy ) SS , where ( s ( el ) xy ) SS , ( s ( vis1 ) xy ) SS , and ( s ( vis2 ) xy ) SS are given by ( s ( el ) xy ) SS = ˙ g T p Z ¥ dth ( ˙ g t ) Z ¥ dk k V ( el ) s ( ¯ k ( t )) V ( el ) s ( k ) F k ( t ) , (34) ( s ( vis1 ) xy ) SS = − ˙ g p Z ¥ dth ( ˙ g t ) Z ¥ dk k V ( vis ) s ( ¯ k ( t )) V ( el ) s ( k ) F k ( t ) ˙ F k ( t ) , (35) ( s ( vis2 ) xy ) SS = ˙ g p Z ¥ dth ( ˙ g t ) Z ¥ dk k V ( el ) s ( ¯ k ( t )) V ( vis ) s ( k ) F k ( t ) Y k ( t ) . (36)The vertex functions V ( el ) s ( q ) and V ( vis ) s ( q ) are given by V ( el ) s ( q ) = ( ¶ ln S q / ¶ q ) / S q and V ( vis ) s ( q ) = − m h I ( ) s ( q ) − I ( ) s ( q ) i , with I ( ) s ( q ) ≡ d j ( qd ) and I ( ) s ( q ) ≡ d (cid:2) j ( qd ) − j ( qd ) (cid:3) . Note that ( s ( el ) xy ) SS , Eq. (34), is coincident to the steady-state shear stress in the sheared MCT of thermal glassy systems [22, 19].he steady-state energy dissipation rate is given by R SS = z H m ˙ g " R ( loc ) SS + (cid:229) i = R ( i ) SS , (37)where R ( loc ) SS is the local contribution, R ( loc ) SS =
12 ˙ g g ( d ) Nnd T , (38)and R ( i ) SS ( i = , ,
3) are the non-local contributions, R ( ) SS = ˙ g T p V Z ¥ dt ˙ g th ( ˙ g t ) Z ¥ dk k I ( ) s ( k ) V ( el ) s ( ¯ k ( t )) F ¯ k ( t ) ( t ) , (39) R ( ) SS = − ˙ g p V Z ¥ dth ( ˙ g t ) Z ¥ dk k V ( vis ) s ( ¯ k ( t )) V ( el ) s ( k ) F k ( t ) ˙ F k ( t ) , (40) R ( ) SS = ˙ g p V Z ¥ dt ˙ g th ( ˙ g t ) Z ¥ dk k I ( ) s ( k ) V ( vis ) s ( ¯ k ( t )) F ¯ k ( t ) ( t ) Y ¯ k ( t ) ( t ) . (41)Note that R ( loc ) SS , Eq. (38), is formally coinicident to the kinetic theory expression of the energy dissipation rate [6]. Numerical Calculations
The inputs necessary for the calculation are S q and g ( d ) , both at equilibrium. We adopt for S q the Percus-Yevicksolution for hard spheres [23] ( j ≤ .
52) and the numerical results from the MD simulation ( j > . g ( d ) the interpolation formula valid in the range 0 . < j < .
64 [24].In order to verify the validity of our theory, we also perform MD simulations of the Sllod equations for softspheres, Eqs. (1) and (2), under the Lees-Edwards boundary condition [18]. The equations are integrated by the Verletalgorithm. The conditions are: the time step D t = . p m / k , where k is the spring constant, the number of particles N = g p m / k = . × − . The values of z are determined by k and e . Theresults are averaged from 10 independent samples. Time correlation functions–
The left (right) panel of Fig. 1 is the relaxation of the density-density (density-current)time correlation function F q ( t ) ( Y q ( t ) ). The reader is referred to the caption for the conditions. The result of MCTfor F q ( t ) shows a clear two-step relaxation in the nearly elastic cases ( e = 0.99, 0.98), while it cannot be seen forthe dissipative cases ( e = 0.94, 0.92). This is qualitatively consistent with the result of the MD simulations, thoughthe plateau is suppressed even for e = .
999 in MD. Note that our results for MCT and MD which exhibit no plateauin F q ( t ) are consistent with the results of MD reported previously [15]. The result of MCT for Y q ( t ) shows a clearpeak and a shoulder for the dissipative cases, while it has only a peak in short time scales for the nearly elastic cases.This is again qualitatively consistent with MD, though the amplitude of MD is much smaller than that of MCT. Itis interesting to note that a plateau or a shoulder appears complementarily in F q ( t ) and Y q ( t ) . This is an evidenceof the significant role of Y q ( t ) in sheared granular liquids. From this result, we can derive the following intuitivepicture when the dissipation is significant. The spheres lose density correlations because they cease to collide with thesurrounding caging spheres at the time scale ∼ m / z H due to inelasticity. On the other hand, the spheres "memorize"the information of the current at the onset of shearing and inelastic collisions, since they lose the kinetic energy morerapidly when they interact with currents which have larger relative velocity. The spheres become crowded where thedissipation rate is large, which results in a density fluctuation; i.e., the density fluctuation is correlated with the initialcurrent, even at the time scales of slow motions. This correlation is eventually destroyed by shearing at the time scale ∼ ˙ g − > m / z H . In contrast, for elastic or driven cases, the information of the initial current is lost and no densityfluctuation is generated. Shear viscosity–
The density dependence of the shear viscosity h = s xy / ˙ g is shown in the left panel of Fig. 2. Inthis figure, we also plot the result of MD for hard spheres [8]. We find a remarkable agreement between MCT andMD for j < .
60 by using the shift of the density in MCT as j → j eff ≡ j + D j with
D j = .
11. This leads tothe MCT transition at the effective density, j MCT , eff = j MCT + D j ≃ . j g ≃ . − .
60, it is close to 0.62, which is used to explain the result of the simulations for
IGURE 1.
The relaxation of the time correlation functions. The left (right) panel is for the density-density (density-current)time correlation function F q ( t )( Y q ( t )) . The results of MCT for j = .
52 and e =0.99, 0.98, 0.94, 0.92 are shown in thick solid,broken, dashed, and thin solid lines. The results of MD for j = .
63 and e =0.9999, 0.999 are shown in triangles and diamonds. Thewave number qd = . Y q ( t ) is multiplied by a factor 30 to be comparable with the results ofMCT in a single figure. colloidal hard-sphere dispersions by MCT, for both the long-time self-diffusion coefficient [25] and the shear viscosity[26]. The results of MCT approximately scales as h (cid:181) ( j J − j eff ) − a for j eff < .
60, where a cross-over from a = j J . However, for j eff > . a becomessmaller and MCT fails to reproduce the divergence observed in MD. In fact, the result of MD shows a cross-over at j ≃ .
60 to a stronger divergence as j → j J . Comparison with the related works on the cross-over of the exponentfrom 1 to 4 for 2D disks [28] and understanding this cross-over is a future task. The divergent feature for j → j J isexpected to originate in the contact network of the spheres, which is not considered in the present framework. Hence,it is reasonable to find a discrepancy between MCT and MD in the vicinity of the jamming transition density. Temperature–
The density dependence of the temperature is shown in the right panel of Fig. 2, together withthe granular temperature measured in MD. The results of MCT show that the temperature monotonically increasesup to j ≃ .
60 in accordance with MD, although the magnitudes are systematically larger. However, qualitativediscrepancies are found for j > .
60, where T decreases. This indicates that the balance of the shearing and theenergy dissipation, which determines T , is invalid for j > .
60 in this framework .
DISCUSSIONSProblem of Canonical Initial Distribution
First we discuss the problem of the canonical initial distribution. As discussed in the previous section, we caninfer that the correct balance between the shearing and the dissipation is not realized in MCT. To inspect in detail,let us observe the steady-state energy dissipation rate R SS , Eq. (37). While the local contribution z H R ( loc ) SS / ( m ˙ g ) isindependent of ˙ g , the non-local contributions with time integrations, z H R ( i ) SS / ( m ˙ g ) ( i = , , g or ˙ g . This indicates that the purely dissipative contribution is absent in z H R ( i ) SS / ( m ˙ g ) ( i = , , h A ( t ) W ( ) i in the steady-state formula, Eq. (12), is approximated as h A ( t ) W ( ) i ≈ (cid:10)(cid:2) U ( t , ) P ( t ) A ( t ) (cid:3) (cid:2) P ( ) W ( ) (cid:3)(cid:11) , (42)where U ( t , ) is the time evolution operator in the space orthogonal to the projected space [29], and P ( t ) isthe second projection operator, Eq. (16), restricted to zero-wavevector. The specific feature of the canonical initial IGURE 2.
The density dependence of the shear viscosity (left panel) and the temperature (right panel). The results of MCTfor e =0.98, 0.96, 0.94, 0.92 are shown in thick solid, dashed, thin solid, and dotted lines, where the density is shifted as j → j eff = j + D j with
D j = .
11. The results of MD for e =0.98, 0.96, 0.94, 0.92 are shown in diamonds, squares, circles,and triangles. The filled symbols are from the MD we have performed, and the open symbols are cited from Ref. [8]. distribution is that the purely dissipative contribution in the work function, i.e. D R of Eq. (33), is projected out by P , i.e. P ( ) D R ( ) = , (43)which results in h A ( t ) D R ( ) i = . This indicates that the purely dissipative effect is missing in the steady-stateformula. This observation suggests us to adopt a non-canonical initial distribution.In general, it is expected that the steady-state averages are insensitive to the choice of the initial distribution. In fact,it is shown explicitly that this is true for some special cases [30]. However, this feature is violated in the approximateformulation of MCT. This is because MCT is intended to describe long-time dynamics, and sacrifices the accuracy forthe description of the initial relaxation. One evidence we have already seen is its dependence on the initial temperature,which appears in the coefficients. Thus, the choice of the initial distribution is crucial in MCT for sheared granularliquids. A basic idea to remedy this problem is to choose the initial distribution which is much closer to the steady-state distribution. To obtain an exact steady-state distribution is of course an extremely difficult task which is out ofour scope at present. However, we might be able to speculate a valid distribution.
Relation to Previous Works
Next we discuss the relation of this work to the previous related works. In the former MCTs for thermal shearedglassy systems [22] or randomly driven granular systems [14], only the projection to the pair-density modes, i.e. P nn of Eq. (17), has been included. In fact, we have already shown in Ref. [31] that the effect of the projection to thedensity-current modes, i.e. P n j of Eq. (17), is negligible in thermal sheared underdamped systems [29]. However, wehave figured out that including the correlations to the density-current modes, or, equivalently, considering Y q ( t ) as wellis crucial to take into account the dissipation effect for sheared granular liquids. This can be convinced by observingthat the dissipative force is proportional to the relative velocity of spheres, and hence dissipation is correlated withcurrent fluctuations. Scaling of the Shear Stress
As for the scaling of the shear viscosity with the shear rate, we only obtain the Bagnold scaling h (cid:181) ˙ g . This seemsto originate from the hard-core limit we have adopted. To reproduce a departure from the Bagnold scaling, which issignificant near the jamming transition point, we should take into account the soft-core nature in our theory [27]. Morespecifically, it might be necessary to further incorporate the information of the contact networks formed between theparticles. However, constructing a theory based on the generalized Green-Kubo formula is still promising, because itholds even above the jamming transition point [32]. ONCLUSION
We have succeeded in predicting the time correlations and the shear viscosity of dense sheared granular liquids byextending MCT, although the prediction of time correlations requires refinement to improve quantitative accuracy. Ithas been demonstrated that, in contrast to thermal glassy systems, the density-current correlation plays an essentialrole. Its validity for the prediction of the shear viscosity has been verified for j < .
60 by comparing it with MD,with the aid of the shift of the density. This scheme is expected to be further extended to higher densities j > .
60 byincorporating the soft-core nature of the spheres. On the other hand, our theory fails to predict the appropriate tendencyof the granular temperature for j > .
60. We have discussed that this problem might reside in the choice of the initialdistribution. To remedy this problem, it seems to be crucial to adopt a distribution which is expected to be close to thesteady-state distribution.
ACKNOWLEDGMENTS
The authors are grateful to S.-H. Chong and M. Otsuki for their collaboration in the initial stage of this project andextensive discussions. They are also grateful to M. Sperl, W. T. Kranz, and M. Fuchs for fruitful discussions, and T.G. Sano and S. Takada for providing us the prototype of the program for the MD simulation. This work is partiallysupported by the Grant-in-Aid of MEXT (Grant No. 25287098). The numerical calculations in this work were carriedout at the computer facilities at the Yukawa Institute.
REFERENCES
1. H. M. Jaeger, S. R. Nagel, and R. P. Behringer,
Rev. Mod. Phys , 1259–1273 (1996).2. A. J. Liu, and S. R. Nagel, Nature (London) , 21–22 (1998).3. A. Ikeda, L. Berthier, and P. Sollich,
Phys. Rev. Lett. , 018301 (2012).4. R. A. Bagnold,
Proc. R. Soc. Lond. A , 49–63 (1954).5. J. T. Jenkins, and M. W. Richman,
Arch. Rat. Mech. Anal. , 355–377 (1985).6. V. Garzó, and J. W. Dufty, Phys. Rev. E , 5895–5911 (1999).7. K. Saitoh, and H. Hayakawa, Phys. Rev. E , 021302 (2007).8. N. Mitarai, and H. Nakanishi, Phys. Rev. E , 031305 (2007).9. J. T. Jenkins, and D. Berzi, Gran. Matt. , 151–158 (2010).10. W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory , Oxford University Press, Oxford, 2009.11. W. Kob, and H. C. Andersen,
Phys. Rev. Lett. , 1376–1379 (1994).12. U. Bengzelius, W. Götze, and A. Sjölander, J. Phys. C , 5915–5934 (1984).13. J. M. Gordon, J. H. Gibbs, and P. D. Fleming, J. Chem. Phys. , 2771–2778 (1976).14. W. T. Kranz, M. Sperl, and A. Zippelius, Phys. Rev. E , 022207 (2013).15. M. P. Ciamarra, and A. Coniglio, Phys. Rev. Lett. , 235701 (2009).16. H. Hayakawa, and M. Otsuki,
Prog. Theor. Phys. , 381–396 (2008).17. K. Suzuki, and H. Hayakawa,
AIP Conf. Proc. , 670–673 (2013).18. D. J. Evans, and G. P. Morriss,
Statistical Mechanics of Nonequilibrium Liquids, 2nd ed. , Cambridge University Press,Cambridge, 2008.19. K. Miyazaki, D. R. Reichman, and R. Yamamoto,
Phys. Rev. E , 011501 (2004).20. O. Henrich, F. Weysser, M. E. Cates, and M. Fuchs, Phil. Trans. R. Soc. A , 5033–5050 (2009).21. M. Fuchs, and M. E. Cates,
J. Rheol. , 957–1000 (2009).22. M. Fuchs, and M. E. Cates, Phys. Rev. Lett. , 248304 (2002).23. J.-P. Hansen, and I. R. McDonald, Theory of Simple Liquids, 3rd ed. , Academic Press, Oxford, 2006.24. S. Torquato,
Phys. Rev. E , 3170–3182 (1995).25. A. J. Banchio, J. Bergenholtz, and G. Nägele, Phys. Rev. Lett. , 1792–1795 (1999).26. Z. Cheng, J. Zhu, P. M. Chaikin, S.-E. Phan, and W. B. Russel, Phys. Rev. E , 041405 (2002).27. A. Ikeda, and L. Berthier, Phys. Rev. E , 052305 (2013).28. M. Otsuki, H. Hayakawa, and S. Luding, Prog. Theor. Phys. Suppl. , 110–133 (2010).29. K. Suzuki, and H. Hayakawa,
Phys. Rev. E , 012304 (2013).30. H. Hayakawa, S.-H. Chong, and M. Otsuki, AIP Conf. Proc. , 19–30 (2010).31. K. Suzuki, and H. Hayakawa,
AIP Conf. Proc. , 750–757 (2013).32. H. Hayakawa, and M. Otsuki,
Phys. Rev. E88