Kinetic Theory of Response Functions for the Hard Sphere Granular Fluid
aa r X i v : . [ c ond - m a t . s o f t ] A ug Kinetic Theory of Response Functions for the Hard SphereGranular Fluid
Aparna Baskaran
Physics Department, Syracuse University, Syracuse, NY 13244
James W. Dufty
Department of Physics, University of Florida, Gainesville, FL 32611
J. Javier Brey
F´ısica Te´orica, Universidad de Sevilla,Apartado de Correos 1065, E-41080, Sevilla, Spain (Dated: November 15, 2018)
Abstract
The response functions for small spatial perturbations of a homogeneous granular fluid have beendescribed recently. In appropriate dimensionless variables, they have the form of stationary statetime correlation functions. Here, these functions are expressed in terms of reduced single particlefunctions that are expected to obey a linear kinetic equation. The functional assumption requiredfor such a kinetic equation, and a Markov approximation for its implementation are discussed.If, in addition, static velocity correlations are neglected, a granular fluid version of the linearizedEnskog kinetic theory is obtained. The derivation makes no a priori limitation on the density, spaceand time scale, nor degree of inelasticity. As an illustration, recently derived Helfand and Green-Kubo expressions for the Navier-Stokes order transport coefficients are evaluated with this kinetictheory. The results are in agreement with those obtained from the Chapman-Enskog solution tothe nonlinear Enskog kinetic equation.
PACS numbers: 05.20.Dd,45.70.-n,05.60.-k,47.10.ab . INTRODUCTION One of the most productive methods to study transport in normal fluids is through themeasurement, simulation, and theory of linear response functions [1–4]. Of particular inter-est are those that describe the linear response of the “hydrodynamic fields” (mass, energy,and momentum densities) to small spatial perturbations of the homogeneous equilibriumstate. The terminology, hydrodynamic fields, is due to the fact that these are the variablesexpected to obey the phenomenological hydrodynamic equations on large space and timescales. The response functions provide the means to study such hydrodynamic excitationsstarting from their fundamental basis in non-equilibrium statistical mechanics. For exam-ple, they provide the formally exact Helfand [5] and Green-Kubo [4, 6] representations fortransport coefficients. More generally, the response functions describe the broader range ofexcitations on shorter space and time scales as well. One of the most instructive theoreticalapproaches to their evaluation has been kinetic theory, with the greatest progress made forthe idealized fluid of hard spheres [4, 7].Recently, this linear response approach has been extended to granular fluids [8, 9]. Theobjective here is to demonstrate the application of kinetic theory methods for the evalua-tion of the granular response functions. Only the case of smooth, inelastic hard spheres isconsidered both for simplicity and to parallel closely the corresponding developments fornormal fluids. Such an idealized model still captures the most important features of manygranular fluids [10]. The usual notion of kinetic theory is a nonlinear equation for the prob-ability density in single particle phase space. An advantage of the linear response functionsis that their kinetic equation is inherently linear. For practical purposes, a Markovian ap-proximation to this linear kinetic equation is described, based on the neglect of dynamicalcorrelations. The approximations leading to this equation, and the differences between itsimplications for normal and granular fluids, are discussed. The nature of the approximationdoes not a priori assume weak dissipation, low density, or large length and time scales. Inthe elastic limit, it becomes the linear Enskog kinetic equation for the response functions ofa normal equilibrium fluid. A granular Enskog limit is described here as well, by the furtherneglect of all velocity correlations in the Markov approximation.A related set of time correlation functions involving the fluxes of the hydrodynamic fields,instead of the fields themselves, determine the Helfand and Green-Kubo representations2or the transport coefficients of the Navier-Stokes order hydrodynamic equations [8, 9].These can be evaluated by the same kinetic theory developed for the response functions.This is described in detail for the shear viscosity in the Markov approximation, indicatinghow the time dependence of these flux correlation functions can be determined, as well asthe associated transport coefficient. The remaining Navier-Stokes transport coefficients areevaluated in the Appendix D. It is confirmed that, in the linear Enskog approximation, theresults agree in detail with those obtained from the Chapman-Enskog method to solve thenonlinear Enskog equation for the distribution function [11].A similar program has already been carried out in the limited context of dynamics foran impurity particle in a granular fluid [12, 13]. In that case, the only hydrodynamic fieldis the impurity particle probability density, and the response function is its autocorrelationfunction. The hydrodynamic equation is a diffusion equation, and the Green-Kubo expres-sion for the diffusion coefficient is given by the time integral of the velocity autocorrelationfunction. This was evaluated by kinetic theory in the Enskog approximation, and both thecorrelation function and the transport coefficient were compared with molecular dynam-ics simulation data over a wide range of densities and degrees of inelasticity. The resultsprovide an instructive characterization of the domain of validity for the Enskog (Enskog-Lorentz in this case) kinetic equation, and expose important differences between impuritydynamics in normal and granular fluids. The presentation here constitutes an extension ofthat theoretical analysis to the full range of multi-particle mass, energy, and momentumtransport.The origin of a kinetic theory is an exact hierarchy (the BBGKY hierarchy [4, 7]) ofequations for the reduced few particle representations for a property of interest. A kineticequation is comprised of the first hierarchy equation together with a “closure”, expressingthe solution to the second hierarchy equation as a functional of that for the first. This leadsto a closed, deterministic equation for the latter which is the kinetic equation. Practicalmethods have been developed for normal fluids based on inversion of cluster expansions andpartial resummations, as well as more phenomenological estimates. This is the point atwhich kinetic theory confronts the difficult many body problem, and one objective of thecurrent work is to motivate a corresponding attention to such details for the granular fluid.Only in this way can the qualitative speculations about differences between normal andgranular fluid be made more precise. An example of such uncertainties is the role of velocity3orrelations in the construction of the closure. It is well known that velocity correlationsgenerated by collective many particle collisions (e.g., ring collisions) are responsible forthe dominant density dependence of transport coefficients at very high densities, but theyare relatively unimportant at low to moderate densities where simple spatial correlations(e.g., excluded volume effects) are dominant. The latter are incorporated in the Enskogapproximation for accurate corrections to the predictions of the low density Boltzmannequation. Granular fluids introduce a complication to this separation of dynamical velocitycorrelations and static structural correlations, according to the density considered. In thesesystems, there are inherent static velocity correlations, not directly associated with manyparticle dynamics, that are present even at low to moderate densities. The kinetic theoryfor response functions provides an appropriate setting for the study of the quantitativeimportance of these correlations on properties of interest.The response functions in the Markovian approximation are expressed in terms of the lin-ear generator for dynamics. Questions about the existence and dominance of hydrodynamicscan be made precise at this point, by asking if the hydrodynamic modes (eigenvalues) appearin the spectrum of this operator and if they are the slowest modes. The first part, existenceof hydrodynamic modes, can be demonstrated at long wavelengths if expected conditions ofanalyticity are satisfied (see Sec. VI below). The second issue of dominance at long timescan be addressed practically using simplified kinetic models for this generator. Such resultssupport the primary assumptions in references [8] and [9], for a derivation of formally exactexpressions for the hydrodynamic transport coefficients in terms of time correlation func-tions. However, the linear kinetic equation applies as well to short space and time scales,that are important for non-hydrodynamic response.The main points of this analysis are summarized in the last section.
II. LINEAR RESPONSE FUNCTIONS
An idealized granular fluid of N smooth, inelastic hard spheres ( d = 3) or disks ( d = 2)of mass m and diameter σ is considered. The inelasticity is characterized by a constantcoefficient of normal restitution α . The properties of interest are the average number density n ( r , t ), the granular temperature T ( r , t ), and the local flow velocity U ( r , t ), collectivelydenoted by y ( r , t ) ≡ { y β ( r , t ) } . The response to be studied here is the variation of these4elds at time t due to a variation in their initial values at time t = 0. The initial conditionsare spatial variations of a homogeneous reference state, y β ( r ,
0) = y β,h (0) + δy β ( r , δy β [ r , t | y (0)] depends nonlinearly on the perturbations δy ( r , δy β [ r , t | y (0)] = X γ Z d r ′ C βγ ( r − r ′ , t ) δy γ ( r ′ , , (1)with the linear response functions defined as C βγ ( r − r ′ , t ) = (cid:20) δy β [ r , t | y (0)] δy γ ( r ′ , (cid:21) δy (0)=0 . (2)Since the reference state is homogeneous, the linear response functions depend on r and r ′ only through their difference.The difference between granular and normal fluids occurs already at the level of thehomogeneous reference state. For normal fluids, this is the equilibrium stationary state andall the time dependence of the response functions is due to the spatial perturbations. Forgranular fluids, the homogeneous reference state is inherently time-dependent, even withoutperturbation, due to the “cooling” of inelastic collisions. As a consequence, the temperature T h of the homogeneous granular fluid decreases in time according with a cooling law ∂T h ( t ) ∂t = − ζ [ T h ( t )] T h ( t ) , (3)where the cooling rate ζ ( T h ) is a characteristic function of the homogeneous cooling ref-erence state. Thus the relevant response at time t is measured relative to the referencehomogeneous state at the same time rather than to the initial state. Then, dimensionlessfields δy ∗ β are introduced by (cid:8) δy ∗ β (cid:9) ≡ (cid:26) δy β y β,h ( t ) (cid:27) ≡ (cid:26) δnn h , δTT h ( t ) , δ U v ( t ) (cid:27) , (4)where the definition of the reference fields y β,h ( t ) follows from the second identity and v ( t ) ≡ (cid:20) T h ( t ) m (cid:21) / (5)is a thermal velocity. The dependence of the linear response functions on r − r ′ suggeststhe utility of a Fourier representation of Eq. (1) that is expressed in the form δ e y ∗ β ( k ∗ , s ) = X γ e C ∗ βγ ( k ∗ , s ) δ e y ∗ γ ( k ∗ , , (6)5here a tilde over a function denotes its dimensionless Fourier transform, defined by e f ( k ∗ ) ≡ ℓ − d Z d r e i k ∗ · r /ℓ f ( r )= Z d r ∗ e i k ∗ · r ∗ f ( r ) , (7)with r ∗ = r /ℓ , and dimensionless response functions have been identified as e C ∗ βγ ( k ∗ , s ) = y − β,h ( t ) ℓ d e C βγ ( k ∗ , t ) y γ,h (0)= y − β,h ( t ) Z d r e i k ∗ · r /ℓ C βγ ( r , t ) y γ,h (0) . (8)Here k ∗ = k ℓ is a dimensionless wavevector, with ℓ being a characteristic length of thesystem. Moreover, a dimensionless time scale s = Z t dt ′ v ( t ′ ) ℓ (9)has been introduced. This time s has the interpretation of an average collision number perparticle up to time t , if ℓ is chosen to be the mean free path of the particles. The remainderof this presentation focuses on the dimensionless linear response functions e C ∗ βγ .For small k ∗ and large s , the functions δy ∗ β ( k ∗ , s ) are expected to obey the linearizedhydrodynamic Navier-Stokes equations, and the response functions correspond to the Greenfunctions for the solution to the initial value problem associated with those equations [9].That description is only phenomenological, since it is parameterized by the unknown trans-port coefficients. A more complete and exact description for all k ∗ and s is provided bynon-equilibrium statistical mechanics, as described in references [8] and [9]. Briefly, theconstruction is as follows.The hydrodynamic fields y β ( r , t ) are defined from averages of the microscopic numberdensity, energy density, and momentum density over the phase space density ρ (Γ , t ), repre-senting the probability that the positions q r and velocities v r of the particles have specifiedvalues denoted by Γ ≡ { x , . . . , x N } . The latter is a point in the 2 N d dimensional phasespace with the notation x r ≡ { q r , v r } . For any specified initial state, ρ (Γ , (cid:18) ∂∂t + L (cid:19) ρ (Γ , t ) = 0 , (10)with L = N X r =1 v r · ∂∂ q r − N X r =1 N X s = r T ( x r , x s ) . (11)6he time independent operator L is the generator for the hard sphere dynamics, where thesingular binary collisions are described by T ( x r , x s ) = δ ( q rs − σ ) | b q rs · g rs | (cid:2) Θ ( b q rs · g rs ) α − b − rs − Θ ( − b q rs · g rs ) (cid:3) . (12)In this expression, q rs = q r − q s , g rs = v r − v s , Θ is the Heaviside step function, b q rs ≡ q rs /q rs ,and b − rs is the substitution operator that replaces the velocities v r , v s by their “precollisional”values v ′′ r , v ′′ s , b − rs F ( v r , v s ) = F ( v ′′ r , v ′′ s ) , (13) v ′′ r = v r − α α ( b q rs · g rs ) b q rs , (14) v ∗′′ s = v s + 1 + α α ( b q rs · g rs ) b q rs . (15)For an isolated system, instead of the equilibrium state for a normal fluid, there is aspacial solution to the Liouville equation of the form ρ h (Γ , t ) = [ ℓv ( t )] − Nd ρ ∗ h ( { q rs /ℓ, v r /v ( t ) } ) ≡ [ ℓv ( t )] − Nd ρ ∗ h (Γ ∗ ) . (16)The dimensionless phase point Γ ∗ ≡ { x ∗ , . . . , x ∗ N } is now expressed in terms of the scaledpositions and velocities x ∗ r = { q ∗ r , v ∗ r } ≡ { q r /ℓ, v r /v ( t ) } . This solution depends on the po-sitions only through the relative variables q rs and, therefore, it has translational invariance,representing a homogeneous state. All of the time dependence occurs through the thermalvelocity defined by Eqs. (5) and (3), so the Liouville equation becomes for ρ ∗ h (Γ ∗ ) L ∗ ρ ∗ h (Γ ∗ ) = 0 . (17)The operator L ∗ is the sum of the original generator for trajectories, now in the dimensionlessvariables, plus a scaling operator representing the time dependenc of v ( t ) using Eq. (3), L ∗ = L ∗ + ζ ∗ N X r =1 ∂∂ v ∗ r · v ∗ r , (18)where ζ ∗ is the dimensionless cooling rate ζ ∗ = ℓv ( t ) ζ [ T h ( t )] . (19)The solution to Eq. (17) will be referred to as the homogeneous cooling state (HCS).7or more general solutions to the Liouville equation, it is useful to introduce the samescaled velocities to account for this inherent cooling. Then, the Liouville equation in dimen-sionless form becomes (cid:18) ∂∂s + L ∗ (cid:19) ρ ∗ (Γ ∗ , s ) = 0 . (20)In this form, it is seen that the HCS is a stationary solution to the dimensionless Liouvilleequation. This is an important result for the representation of granular response functions. Itshows that, in the appropriate dimensionless form, the reference state is again stationary, justas for equilibrium fluids. However, the introduction of this stationary representation comesat the price of changing the generator for the dynamics from L ∗ to L ∗ . For the purposes ofthe discussion here, it is assumed that all properties of this homogeneous reference state areknown. Further comments on the HCS are given in the Appendix A.The deviations of the relative hydrodynamic fields form their values in the HCS, δ e y ∗ β ( k ∗ , s ), are the averages of associated phase functions e a ∗ β (Γ ∗ , s ), δ e y ∗ β ( k ∗ , s ) = Z d Γ ∗ e a ∗ β (Γ ∗ ; k ∗ ) [ ρ ∗ (Γ ∗ , s ) − ρ ∗ h (Γ ∗ )] , (21)with e a ∗ β (Γ ∗ ; k ∗ ) = 1 n ∗ h N X r =1 e i k ∗ · q ∗ r a ∗ β ( v ∗ r ) , (22)where n ∗ h ≡ n h ℓ d is the number density of the HCS in the reduced units and the singleparticle functions a ∗ β ( v ∗ r ) are defined by (cid:8) a ∗ β ( v ) (cid:9) ≡ (cid:26) , v ∗ d − , v ∗k , v ∗⊥ (cid:27) . (23)For later convenience, the components of the flow field and the velocity of the particles havebeen chosen to be a longitudinal component along k ∗ , and d − v ∗k ≡ b k · v ∗ and v ∗⊥ i ≡ b e i · v ∗ , with nb k ≡ k ∗ /k ∗ , b e i ; i = 1 , . . . , d − o forming a setof d pairwise perpendicular unit vectors. To formulate the linear response problem, theinitial state ρ ∗ (Γ ∗ ,
0) for the solution to the Lioville equation in (20) is chosen to be closeto the HCS, in the sense that it is a functional of the initial fields δy ∗ β ( r ∗ ,
0) and becomesthe HCS for δy ∗ β →
0. More specifically, the system is viewed as partitioned into small cellssuch that the distribution function is the HCS in each cell, but with different values for thehydrodynamic fields. This is the analogue of the local equilibrium distribution for normalfluids, and will be referred to as the local
HCS, ρ ∗ lh [Γ ∗ | δy ∗ ]. Its construction is discussed in8ppendix A. For small δy ∗ β , the solution to the Liouville equation for the initial local HCSexpanded to first order is ρ ∗ (Γ ∗ , s ) − ρ ∗ h (Γ ∗ ) = e − s L ∗ ρ ∗ (Γ ∗ , → e − s L ∗ d +2 X γ =1 Z d r ∗ (cid:20) δρ ∗ ℓh [Γ ∗ | δy ∗ ] δy ∗ γ ( r ∗ ) (cid:21) δy ∗ =0 δy ∗ γ ( r ∗ , . (24)Substitution into Eq. (21) allows identification of the response functions C ∗ βγ ( r ∗ , s ) and theirequivalent Fourier representation (using translational invariance of the functional derivativeat δy ∗ = 0), e C ∗ βγ ( k , s ) = Z d Γ ∗ e a ∗ β (Γ ∗ ; k ∗ ) ψ ∗ γ (Γ ∗ ; ∗ , s ) , (25)with ψ ∗ γ (Γ ∗ ; r ∗ , s ) = e − s L ∗ ψ ∗ γ (Γ ∗ ; r ∗ ) (26)and ψ ∗ γ (Γ ∗ ; r ∗ ) = (cid:20) δρ ∗ ℓh [Γ ∗ | δy ∗ ] δy ∗ γ ( r ∗ ) (cid:21) δy =0 . (27)These results are closely analogous to those for a normal fluid, where ρ ∗ ℓh [Γ ∗ | δy ∗ ] is an initiallocal equilibrium ensemble and ρ ∗ h (Γ ∗ ) is the reference equilibrium Gibbs ensemble, ρ ∗ e (Γ ∗ ).For example, if the grand canonical ensemble were used, the e ψ ∗ γ ’s would become linearcombinations of the set (cid:8)e a ∗ γ (cid:9) times the equilibrium ensemble, and the response functionswould be equilibrium time correlation functions for the local conserved quantities. III. KINETIC THEORY
All of the following analysis is carried out in terms of the dimensionless variables, so theasterisk will be left implicit for simplicity. The response functions of Eq. (25) are expressedin terms of the full N -particle phase space. A reduced description in terms of the singleparticle phase space is possible because the functions e a β (Γ , k ) defined in Eq. (22) are sumsof single particle functions. Consequently, integrating over the positions and momenta forall except one particle leads to the exact alternative representation: e C βγ ( k , s ) = 1 n h Z d v a β ( v ) e ψ (1) γ ( v , k , s ) , (28)where e ψ (1) γ ( v , k , s ) = Z d q e i k · q ψ (1) γ ( x , s ) . (29)9he function ψ (1) γ ( x , s ) is the first element of a hierarchy of N functions ψ ( m ) γ ( x , x , . . . , x m , s ), m = 1 , . . . , N , defined through ψ ( m ) γ ( x , . . . , x m , s ) ≡ (cid:20) δf ( m ) ( x , . . . , x m , s ) δy γ ( , (cid:21) δy =0 , (30) f ( m ) ( x , . . . , x m , s ) ≡ N !( N − m )! Z dx m +1 ..dx N ρ (Γ , s ) , (31)1 ≤ m ≤ N . The functions f ( m ) ( x , . . . , x m , s ) are the reduced distribution functions associ-ated with the solution to the Liouville equation (20). They obey the corresponding BBGKYhierarchy of equations [4]. The first equation of this hierarchy is (cid:18) ∂∂s + v · ∂∂ q (cid:19) f (1) ( x , s ) + ζ ∂∂ v · (cid:2) v f (1) ( x , s ) (cid:3) = Z dx T ( x , x ) f (2) ( x , x , s ) . (32)Then, it follows directly from the definition (30) that the function ψ (1) γ ( x , s ) obeys theanalogous equation (cid:18) ∂∂s + v · ∂∂ q (cid:19) ψ (1) γ ( x , s ) + ζ ∂∂ v · (cid:2) v ψ (1) γ ( x , s ) (cid:3) = Z dx T ( x , x ) ψ (2) γ ( x , x , s ) . (33)The representation given by Eq. (28) is very appealing, since the N -particle problem hasbeen expressed without approximation in terms of the effective dynamics in the single par-ticle phase space. The fundamental difficulty, however, is that the first hierarchy equation(33) does not determine this effective dynamics without specifying ψ (2) γ ( x , x , s ). An equa-tion similar to (33) can be written for ψ (2) γ ( x , x , s ), the second BBGKY hierarchy equation,but it in turn requires specification of ψ (3) γ ( x , x , x , s ). In this way, a coupling to the full N particle problem recurs. This coupling is broken if ψ (2) γ ( x , x , s ) can be specified as anexplicit functional of ψ (1) γ ( x , s ) . It is argued in Appendix B that this functional, when itexists, is independent of δy γ ( r , ψ (1) γ ( x, s ), and has the general form ψ (2) γ ( x , x , s ) = Z dx K ( x , x , s ; x ) ψ (1) γ ( x, s ) . (34)The kernel defining the functional is K ( x , x , s ; x ) = (cid:20) δf (2) ( x , x , s ) δf (1) ( x, s ) (cid:21) δy =0 . (35)Once K ( x , x , s ; x ) is known, Eq. (33) becomes a closed, deterministic, linear kinetic equa-tion for ψ (1) γ in its most general form, (cid:20) ∂∂s + v · ∂∂ q + M ( s ) (cid:21) ψ (1) γ ( x , s ) + ζ ∂∂ v · (cid:2) v ψ (1) γ ( x , s ) (cid:3) = 0 , (36)10ith the formal “collision operator” M given by M ( x , s ) ψ (1) γ ( x , s ) ≡ − Z dx (cid:20)Z dx T ( x , x ) K ( x , x , s ; x ) (cid:21) ψ (1) γ ( x, s ) . (37)This notion of a kinetic equation makes no a priori assumptions regarding the density ordegree of dissipation, as is sometimes assumed. In fact, it is formally exact so any restrictionsarise only when specific approximations are introduced to construct the functional.The utility of the kinetic theory representation in any specific application depends onan accurate construction of the kernel K ( x , x , s ; x ) defined in Eq. (35). This is the pointat which the full many-body problem must be confronted. The idea originated with Bo-goliubov [14], who conceived that all f ( m ) ( x , . . . , x m , s ) with m > f ( m ) [ x , .., x m | f (1) ( s )], after a brief “synchronization” time, where all time dependence occursonly through f (1) ( s ). For normal fluids at low density, construction of this functional canbe accomplished by formal density expansions, leading to a sequence of contributions fromclusters of particles of increasing size [15]. At lowest order, the closure is that associated withthe Boltzmann equation, while at next order three particle scattering is described. Even fornormal fluids, these sophisticated cluster expansions have limited direct use due to many-particle recollisions events that contribute secular terms in the formal expansions, violatingthe notion of a short synchronization time [16]. These recollisions (“rings”) in turn signalnon-analytic density dependence and slow algebraic decay of correlations in time [17]. Theexperience gained in such studies over the past forty years has provided important insight forthe construction of more phenomenological closures, such as the Enskog approximation formoderately dense gases and mode coupling models for very dense and metastable (glassy)fluids [18]. The development of a kinetic theory for granular fluids provides an opportunityto revisit many of these issues in an even more challenging context [19]. The analysis of thenext sections illustrates this for the simplest “mean field” approximation to describe densityeffects beyond the Boltzmann limit. IV. THE MARKOVIAN APPROXIMATION
The approximation described in this section is based on the assumption that the formof the correlations and their effect on collisional properties are essentially the same at alltimes. If so, the kernel K ( x , x , s ; x ) that determines these properties at the two-particle11evel, can be represented approximately by its form at the initial time s = 0, K ( x , x , s ; x ) ≃ K ( x , x , x ) . (38)The collision operator term in Eq. (37) then becomes M ( x ) ψ (1) γ ( x , s ) = − Z dx Z dx T ( x , x ) " δf (2) ℓh [ x , x | δy ] δf (1) ℓh [ x | δy ] δy =0 ψ (1) γ ( x, s ) , (39)where it has been used that the initial conditions f ( m ) ( x , .., x m ,
0) are the correspondingreduced distribution functions associated with the local HCS, ρ ℓh [Γ | δy ], given in Eq. (A8), f ( m ) ( x , .., x m ,
0) = f ( m ) ℓh [ x , .., x m | δy ]) = N !( N − m )! Z dx m +1 ..dx N ρ ℓh [Γ | δy ] . (40)Since the approximate operator M ( x ) is time-independent, the entire generator for thedynamics of ψ (1) γ ( x , s ) in Eq. (36) also is independent of time. The equation is written incompact form as (cid:20) ∂∂s + Λ( x ) (cid:21) ψ (1) γ ( x , s ) = 0 , (41)with the generator Λ( x ) given byΛ( x ) X ( x ) ≡ (cid:20) v · ∂∂ q + M ( x ) (cid:21) X ( x ) + ζ ∂∂ v · [ v X ( x )] , (42)for arbitrary X ( x ). This is a necessary condition for a Markovian description and, conse-quently, Eq. (38) will be referred to as the Markovian approximation. The resulting kinetictheory is exact at asymptotically short times, and the nature of the approximation makes noexplicit limitation on the density or the degree of dissipation. Of course, two particle corre-lations that develop over time are neglected. For a normal fluid, this Markov approximationleads to the Enskog approximation, where only time independent two-particle correlationsare taken into account. The neglected time-dependent correlations are found in that caseto be important only at high densities, and the Enskog approximation provides relevantcorrections to the Boltzmann results up to moderate densities. It is reasonable to expect asimilar context for granular fluids, although conditioned by the additional parameter spaceof the coefficient of restitution α .The construction of K ( x , x , x ) now involves only analysis of the initial reduced dis-tribution functions associated with the local HCS, and the corresponding collision operator12 is entirely characterized by properties of this state. It is useful to make this more explicitin terms of the pair correlation function g (2) ℓh [ x , x | δy ] for the HCS, defined by f (2) ℓh [ x , x | δy ] ≡ f (1) ℓh [ x | δy ] f (1) ℓh [ x | δy ] g (2) ℓh [ x , x | δy ] . (43)Then Eq. (35) in the Markovian approximation yields K ( x , x , x ) = g (2) h ( x , x ) h f (1) h ( v ) δ ( x − x ) + f (1) h ( v ) δ ( x − x ) i + f (1) h ( v ) f (1) h ( v ) " δg (2) ℓh [ x , x | δy ] δf (1) ℓh [ x | δy ] δy =0 . (44)The first term on the right side is given explicitly in terms of the HCS correlations, whilethe second term requires further analysis of the dependence of local HCS correlations on f (1) ℓh . For a normal fluid , these correlations are independent of the velocities and depend onthe local equilibrium distribution f (1) ℓe only through the local density, " δg (2) ℓh [ x , x | δy ] δf (1) ℓh [ x | δy ] δy =0 → n h " δg (2) ℓe [ q , q | δy ] δn ( q ) δy =0 . (45)The functional form of the local equilibrium pair correlation functional g (2) ℓe [ q , q | δy ] isknown, so explicit construction of the collision operator M is possible in this case. Theresult is the linear revised Enskog kinetic equation [20]. For granular fluids, g (2) ℓh [ x , x | δy ]has a more general functional functional dependence on f (1) ℓh through its explicit functionaldependence on δy, since δy γ ( q ) = 1 n h Z dx a γ ( x ) δ ( q − q ) n f (1) ℓh [ x | δy ] − f (1) h ( v ) o , (46)where the a γ ( x )’s are the single particle functions given in Eq. (23). Thus " δg (2) ℓh [ x , x | δy ] δf (1) ℓ [ x | δy ] δy =0 = 1 n h X λ " δg (2) ℓh [ x , x | δy ] δy λ ( q ) δy =0 a λ ( v ) . (47)This provides the practical route for constructing the Markovian kinetic theory for granularfluids. The collision operator becomes M ( x ) ψ (1) γ ( x , s ) = − Z dx T ( x , x ) g (2) h ( x , x ) h f (1) h ( v ) ψ (1) γ ( x , s ) + f (1) h ( v ) ψ (1) γ ( x , s ) i − X λ Z d q c λ ( v , q ) 1 n h Z d v a λ ( v ) ψ (1) γ ( x , s ) , (48)13ith c λ ( v , q ) = Z dxT ( x , x ) f (1) h ( v ) f (1) h ( v ) " δg (2) ℓh [ x , x | δy ] δy λ ( q ) δy =0 . (49)Further discussion and simplification of c λ ( v , q ) is given in Appendix C.The response functions of Eq. (28) are given in terms of e ψ (1) γ ( v , − k , s ), that is propor-tional to the Fourier transform of ψ (1) γ ( x , s ) as indicated in Eq. (29). An equation for thelatter can be easily derived from Eq. (36), (cid:20) ∂∂s + e Λ( k ) (cid:21) e ψ (1) γ ( v , − k , s ) = 0 , (50)where the generator for the dynamics is the Fourier transform of Λ defined in Eq. (42), e Λ( k ) = − i k · v + f M ( k ) + ζ (cid:18) d + v · ∂∂ v (cid:19) . (51)with f M ( v , k ) defined by f M ( v , k ) X ( v ) ≡ − Z dx T ( x , x ) g (2) h ( x , x ) h e i k · q f (1) h ( v ) X ( v ) + f (1) h ( v ) X ( v ) i − X λ e c λ ( v , k ) 1 n h Z d v a λ ( v ) X ( v ) , (52)for arbitrary X ( v ), e c λ ( v , k ) being the Fourier transform of c λ ( v , q ).The response functions in Eq. (28) become e C βγ ( k , s ) = 1 n h Z d v a β ( v ) e − s e Λ( k ) e ψ (1) γ ( v , k ) , (53) e ψ (1) γ ( v , − k ) = Z d q e i k · q " δf (1) ℓh [ x | δy ] δ e y γ ( ) δ e y =0 . (54)Equation (53) is the primary practical result of our analysis here. It provides a realistickinetic theory description of the most fundamental time dependent fluctuations in a granularfluid: those induced by perturbations of the hydrodynamic fields. To appreciate the scopeand generality of this result, note that the Markovian approximation is exact at short timesfor all k , densities, and degrees of restitution. In this short time limit, Eq. (53) yieldslim s → e C βγ ( k , s ) = lim s → (cid:0) e − sN ( k ) (cid:1) βγ , (55)with N βγ ( k ) = 1 n h Z d v a β ( v ) e Λ( k ) e ψ (1) γ ( v , k ) . (56)14t longer times, the Markovian approximation is expected to continue to provide a goodapproximation across this parameter space, since the derivation does not explicitly requireany limitations on k , densities, or degrees of restitution. This expectation is borne out inthe elastic limit, where the linear Enskog theory for equilibrium time correlation functionsis recovered, as discussed in the next section. In that case, comparisons with moleculardynamics simulations (wavevector dependent transport [21]) and neutron scattering experi-ments [22], confirm the accuracy and practical utility of this kinetic theory over a wide rangeof wavevectors and densities. A determination of the corresponding domain of accuracy forgranular fluids, awaits similar comparisons with simulation and experiments, but there is nosimple reason to expect qualitative rather than quantitative differences from normal fluids.In Sec. VI, it is shown that the operator e Λ( k ) contains the hydrodynamic modes in itsspectrum for small k . Thus the response functions in Eq. (53) provide a means to study thetransition from short time dynamics to a presumed dominant hydrodynamics at long timesfor granular fluids, in a manner similar to that done for normal fluids [4]. In addition, sinceit is valid for all k , the nature of hydrodynamics beyond the Navier-Stokes approximationcan be studied.An important application of general linear response methods, is the derivation of Helfandand Green-Kubo expressions for the transport coefficients [8, 9]. These are formally exactresults given in terms of time correlation functions. The above analysis for the responsefunctions applies to these as well, and their dynamics in the Markovian approximation isgenerated by the same operator e Λ( k ). The evaluation of the Helfand and Green-Kuboexpressions for the shear viscosity is illustrated in Sec. VII. The results provide a general-ization of those from Enskog kinetic theory [11], to include pair velocity correlations. Whensuch correlations are neglected, the results of ref. [11] are recovered in detail. This is verifiedfor the other Navier-Stokes transport coefficients as well in Appendix E. V. GRANULAR ENSKOG APPROXIMATION
The Markovian approximation discussed in the previous section, requires specification of g (2) ℓh ( x , x ) and c λ ( x , q ), or equivalently δg (2) ℓh [ x , x | δy ] /δy λ ( q ) for δy = 0. While theseare well defined in terms of the local HCS distribution, little is know about their detailedforms as yet, except in the elastic limit where they are accurately determined from liquid15tate theory of the pair correlation function. The important simplification in that case isthe absence of velocity correlations. It is plausible to assume that such correlations remainweak for the granular fluid as well, and to make the approximation g (2) ℓh [ x , x | δy ] ≃ g (2) ℓh [ q , q | δy ] = g (2) ℓh [ q , q | δn ] . (57)The last equality recognizes that the neglect of velocity correlations leads to a functionalthat is independent of δT and δ U (for hard spheres or disks) and hence is a functional onlyof the density. Furthermore, g (2) ℓh [ q , q |
0] = g (2) h ( q ) as a consequence of fluid symmetry.Then, Eq. (49) reduces to c ( E ) λ ( v , q ) = δ λ c ( E ) ( v , , q ) = δ λ Z dx T ( x , x ) f (1) h ( v ) f (1) h ( v ) " δg (2) ℓh ( q , q | δn ) δn ( q ) δn =0 . (58)The collision operator (52) now simplifies to f M E X ( v ) = − g (2) h ( σ ) Z dx T ( x , x ) h e − i k · q f (1) h ( v ) X ( v ) + f (1) h ( v ) X ( v ) i − n h e c ( E ) ( v , k ) Z d v X ( v ) . (59)It remains to give the explicit density dependence for g (2) ℓh [ r , r | δn ] [23]. As a practicalmatter, it can be chosen to be the pair distribution for a nonuniform normal fluid forwhich a well-developed theory exists. In that case, g (2) h ( | r − r | ) = g (2) eq ( | r − r | ), the radialdistribution function for a uniform hard sphere fluid. Also, for this choice the functionalderivative appearing in the expression of e c ( E ) ( v , k ) can be evaluated for the first few termsof a k expansion, as is required for evaluation of transport coefficients. The context of sucha choice would be that the static spatial correlations of a hard sphere system are due toexcluded volume effects, and these can be captured using the pair correlation function of afluid of elastic hard spheres. The generator Λ( k ) obtained in this approximation from Eq.(59) gives the granular Enskog kinetic theory, and is the linearized version of the one studiedin ref. [11]. VI. HYDRODYNAMIC MODES
An important feature of the response functions considered here is their relationship tohydrodynamic response. At small k and large s , these response functions should correspond16o those from the phenomenological hydrodynamic equations. For example, it is this rela-tionship that allows the derivation of Helfand and Green-Kubo expressions for the transportcoefficients. Any acceptable approximate kinetic theory for the response functions shouldpreserve this relationship to hydrodynamics. More specifically, the hydrodynamic excita-tions should appear in the spectrum of the linear operator e Λ( k ), e Λ( k , v ) φ ( β ) ( k , v ) = λ ( β ) ( k ) φ ( β ) ( k , v ) , β = 1 , . . . , d + 2 , (60)where the set (cid:8) λ ( β ) ( k ) (cid:9) are the eigenvalues of the d + 2 linearized hydrodynamic equations.The above eigenfunctions and eigenvalues are determined in the limit k = in Appendix Dwith the results (cid:8) λ ( β ) ( ) (cid:9) = (cid:26) , ζ , − ζ (cid:27) , (61) (cid:8) φ ( β ) ( v , ) (cid:9) = ( e ψ (1)1 ( v , ) − (cid:18) ∂ ln ζ ln n h (cid:19) T h e ψ (1)2 ( v , ) , e ψ (1)2 ( v , ) , e ψ (1)3 ( v , ) ) . (62)The eigenvalue − ζ / d -fold degenerate, and the associated eigenfunctions are the com-ponents of the vector e ψ (1)3 ( v , ) ≡ − ∂f (1) h ( v ) /∂ v .In the elastic limit, these eigenvalues are all zero, corresponding to the d + 2 conservationlaws, and the eigenfunctions become Maxwellians times linear combinations of the summa-tional invariants (1 , v , v ). For inelastic collisions, the nonzero eigenvalues describe responseof the cooling temperature to linear perturbations and growth of a constant velocity per-turbation relative to the characteristic cooling thermal velocity. In both cases, these arealso the eigenvalues of the phenomenological linearized hydrodynamic equations in the longwavelength limit.With the eigenfunctions and eigenvalues known at k = 0, their values for finite butsmall k can be obtained by perturbation theory. In this way the Navier-Stokes transportcoefficients can be determined directly from the coefficients up through order k . This directcalculation of the spectrum for the generator of a linear kinetic theory has been described indetail recently for the granular Boltzmann equation [24], and its extension to the Markoviankinetic theory given here is straightforward. Instead, the remainder of this presentationaddresses the calculation of the transport coefficients from an approximate evaluation oftheir Helfand and Green-Kubo representations that have been obtained from linear response.17 II. HELFAND AND GREEN-KUBO EXPRESSIONS
As noted above, the exact response functions defined in Sec. II must agree with those ofthe linearized phenomenological hydrodynamic equations in the long wavelength and longtime limit. This relationship allows identification of the parameters of those phenomenolog-ical equations in terms of the response functions in this limit. The results of this analysisfor granular fluids has been given recently, leading to expressions for the transport coeffi-cients in terms of certain time correlation functions derived from the response functions [9].These correlation functions can be evaluated approximately by the Markov kinetic theory, toobtain explicit results for all transport coefficients appearing in the Navier-Stokes hydrody-namic equations. Further, neglecting the velocity correlations in the Markov theory, allowsthe evaluation of these quantities in the granular Enskog theory, reproducing the resultsreported in [11]. In this section, only the shear viscosity is considered as an example, whileall remaining transport coefficients are analyzed in Appendix E.The exact Helfand and Green-Kubo expressions for the shear viscosity are (dimensionlessunits are still assumed) η = lim Ω ηH ( s ) = Ω ηH (0) + lim Z s ds ′ Ω ηG ( s ′ ) , (63)respectively, where the symbol lim denotes the hydrodynamic limit of V → ∞ , followed by s → ∞ . The correlation functions in the above equation are defined byΩ ηH ( s ) = − V − d + d − d X i =1 d X j =1 Z d Γ H ij (Γ) e − s ( L + ζ ) M η,ij (Γ) , (64)Ω ηG ( s ) = ∂∂s Ω ηH ( s ) = − V − d + d − d X i =1 d X j =1 Z d Γ H ij (Γ) e − s ( L + ζ )Υ η,ij (Γ) . (65)Here, H ij (Γ) is the volume integrated momentum flux, H ij (Γ) = N X r =1 v r,i v r,j + N X r =1 N X s = r H (2) ij ( x r , x r ) , (66) H (2) ij ( x r , x s ) = (1 + α ) σ δ ( q rs − σ ) Θ ( − b q rs · g rs ) ( b q rs · g lm ) b q rs,i b q rs,j , (67) M η,ij is the traceless tensor M η,ij = − N X r =1 (cid:18) q ri ∂∂v r,j + q r,j ∂∂v r,i − d δ ij q r · ∂∂ v r (cid:19) ρ h (Γ) , (68)18nd Υ η,ij (Γ) is the associated Green-Kubo conjugate flux,Υ η,ij (Γ) = − (cid:18) L + ζ (cid:19) M η,ij (Γ) . (69)It is seen in Eqs. (64) and (65) that the generator of dynamics is L + ζ . This reflects thefact that the k = 0 mode λ = − ζ / , of Eq. (61), has been subtracted out.The time independent contribution Ω ηH (0) in the Green-Kubo expression, can be evalu-ated exactly from the definitions (66) and (68) with the result:Ω ηH (0) = V − d + d − d X i =1 d X j =1 Z dx h ij ( v ) 12 (cid:18) q ,i ∂∂v ,j + q ,j ∂∂v ,y − d δ ij q · ∂∂ v (cid:19) f (1) h ( v ) , (70)where h ij ( v ) = v i v j + Z dx Z dx H (2) ij ( x , x ) K ( x , x , x ) . (71)The first term in the above expression of h ij gives no contribution to Ω ηH (0), from fluidsymmetry. The second term can be recognized as being proportional to the average collisionfrequency, ν av , as determined by the loss part of the right hand side of the hard sphereBBGKY hierarchy (32) specialized for the HCS,Ω ηH (0) = (1 + α ) σ d + 2 d ) ν av (72) ν av = 2 σ d − Z d b σ Z d v Z d v Θ ( − b σ · g ) ( b σ · g ) f (2) h ( σ , v , v ) . (73) A. Evaluation in the Markov Approximation
A complete evaluation of the correlation function Ω ηH ( s ) is possible using the Markovkinetic theory. As in Sec. III, the correlation function can be given a representation interms of one and two particle functionsΩ ηH ( s ) = − V − d + d − d X i =1 d X j =1 Z dx v ,i v ,j M (1) η,ij ( x , s ) − V − d + d − d X i =1 d X j =1 Z dx Z dx H (2) ij ( x , x ) M (2) η,ij ( x , x , s ) , (74)with the reduced functions M ( m ) η,ij ( x , .., x m , s ) defined by M ( m ) η,ij ( x , . . . , x m , s ) ≡ N !( N − m )! Z dx m +1 . . . Z dx N e − s ( L + ζ ) M η,ij (Γ) . (75)19imilarly to the functions ψ γ considered in Sec. III, the above functions obey a BBGKYhierarchy, the first equation of which is (cid:18) ∂∂s + ζ v · ∂∂ q (cid:19) M (1) η,ij ( x , s ) + ζ ∂∂ v · h v M (1) η,ij ( x , s ) i = Z dx T ( x , x ) M (2) η,ij ( x , x , s ) . (76)The Markovian approximation in the present case is the same as that defined by Eqs. (34)and (38) M (2) η,ij ( x , x , s ) ≃ Z dx K ( x , x , x ) M (1) η,ij ( x, s ) . (77)Then, Eq. (76) becomes the Markovian kinetic equation (cid:18) ∂∂s + ζ (cid:19) M (1) η,ij = 0 , (78)where the linear operator Λ is the same as defined in Eq. (42).In this approximation, the Helfand expression of the shear viscosity of a hard sphere ordisk granular fluid becomes η = lim Ω ηH ( s ) ≃ V − d + d − d X i =1 d X j =1 Z dx h ij ( v ) e − s ( Λ+ ζ ) × (cid:18) q ,i ∂∂v ,j + q ,j ∂∂v ,y − d δ ij q · ∂∂ v (cid:19) f (1) h ( v ) , (79)where h ij ( v ) is defined in Eq. (71).Next, the Green-Kubo expression for η in the Markov approximation can be identifiedfrom Eqs. (63) and (65),Ω ηH (0) = V − d + d + 2 d X i =1 d X j =1 Z dx h ij ( v ) 12 (cid:18) q ,i ∂∂v ,j + q ,j ∂∂v ,y − d δ ij q · ∂∂ v (cid:19) f (1) h ( v )(80)and Ω ηG ( s ) = V − d + d − d X i =1 d X j =1 Z dx h ij ( v ) e − s ( Λ+ ζ ) γ ij ( v ) . (81)The reduced conjugate flux γ ij is γ ij ( v ) = − (cid:18) Λ + ζ (cid:19) (cid:18) q ,i ∂∂v ,j + q ,j ∂∂v ,y − d δ ij q · ∂∂ v (cid:19) f (1) h ( v ) . (82)Comparison of Eqs. (70) and (80) shows that Ω ηH (0) is given exactly in the Markov approx-imation. The Green-Kubo representation for the shear viscosity requires the large s limit of20he integral over s in Eq. (63). It can be verified that Ω ηG ( s ) has no invariant part, so thatthis limit is expected to exist. This issue is discussed in some detail in Appendix E. Thenthe Green-Kubo expression for shear viscosity can be written as η = Ω ηH (0) + 1 d + d − d X i =1 d X j =1 Z d v h ij ( v ) D ij ( v ) , (83)where D ij ( v ) is a solution to the integral equation (cid:18) ζ ∂∂ v · v + M + ζ (cid:19) D ij ( v ) = γ ij ( v ) . (84)Upon writing the above equation, it has been taken into account that terms involving spatialderivatives give a vanishing contribution to the expression of the shear viscosity. This is thetraditional form in which expressions for transport coefficients are obtained from a ChapmanEnskog expansion of a normal solution to the kinetic equation governing the dynamics ofthe system. B. Evaluation in the Granular Enskog Approximation
The further neglect of velocity correlations in the collision operator M , leads to thegranular Enskog approximation, i.e., the results given by Eqs. (79) and (81) apply with onlythe replacement Λ by Λ E , withΛ E ( x ) ≡ v · ∂∂ q + M E ( x ) + ζ ∂∂ v · v , (85)and the operator M E given by, M E X ( x ) ≡ − g h ( σ ) Z dx T ( x , x ) h f (1) h ( v ) X ( x ) + X ( x ) f (1) h ( v ) i − n h Z d q c ( x , , q ) Z d v X ( x ) . (86)The function c is given in Eq. (58). In the Enskog approximation, the conjugate flux γ ij inEq. (82) becomes γ ij ( v ) = − (cid:18) Λ E + ζ (cid:19) (cid:18) q ,i ∂∂v ,j + q ,j ∂∂v ,y − d δ ij q · ∂∂ v (cid:19) f (1) h ( v ) . (87)One final simplification occurs for the shear viscosity and some other transport coeffi-cients. The mean field term in Eq. (86) vanishes when acting on γ ij ( v ) and, therefore, (cid:18) Λ E + ζ (cid:19) γ ij ( v ) = (cid:18) J ( x ) + ζ (cid:19) γ ij ( v ) , (88)21here the operator J ( x ) has been introduced, J ( x ) ≡ ζ ∂∂ v · v − g (2) h ( σ ) I . (89)Here, I is the linearized Boltzmann collision operator for inelastic hard spheres or disks, I X ( x ) ≡ Z dx T ( x , x ) h f (1) h ( v ) X ( x ) + X ( x ) f (1) h ( v ) i . (90)Therefore, the correlation function in Eq. (81) and the expression for the shear viscosity inEq. (83) take the final formsΩ ηEG ( s ) = 1 d + d − " α ) σ d n h π d/ g (2) h ( σ )4Γ (cid:0) d +42 (cid:1) × d X i =1 d X j =1 Z d v v i v j exp (cid:20) − s (cid:18) J + ζ (cid:19)(cid:21) γ ij ( v ) (91)and η = Ω ηEH (0) + 1 d + d − d X i =1 d X j =1 Z d v h ij ( v ) D Eij ( v ) , (92)respectively. The Enskog approximation for Ω ηH (0) isΩ ηEH (0) = π ( d − / (1 + α ) σ d +1 g (2) h ( σ )2( d + 2 d )Γ (cid:0) d +12 (cid:1) Z d v Z d v | v − v | f (1) h ( v ) f (1) h ( v ) (93)and D Eij is a solution to the integral equation (cid:18) J + ζ (cid:19) D Eij ( v ) = γ ij ( v ) , (94)The above result for the shear viscosity agrees in detail with that obtained in ref. [11],through a Chapman-Enskog procedure applied to the non-linear granular Enskog equation.The results in Eqs. (81) and (91) for the Green-Kubo integrand are new. At the formallyexact level, the integrand is given by the correlation between the flux and the conjugate flux.In detail, the contributions from H (2) ij and the T ( i, j ) terms of L , appear to yield singularitiesat t = 0, signaling a possible nonanalytic dependence on t . This occurs even in the elasticlimit, and is a peculiarity of hard particle dynamics. Consequently, previous theoretical andsimulation studies have avoided this by studying the Helfand forms for transport coefficients.The approximate kinetic theory described here gives an explicit analytic estimate for Ω ηG ( s ),whose integral yields a good estimate for the transport coefficients. This suggests that Ω ηG ( s )may have a dominant analytic part with a relatively small non-analytic correction.22 III. DISCUSSION
Kinetic theory has been used extensively as a formal tool for approximate evaluationof response functions, in the study of hard spheres as a model for normal fluids. Theobjective of this work is to take a first step in the development and application of thistool in the analogous field of granular fluids. The advantage of developing kinetic theoryin the context of linear response functions, lies in the fact that the resulting theories areinherently linear, and provide a more tractable setting to explore questions such as agingto hydrodynamics and short wavelength behavior of the exact hydrodynamic response. Thetwo primary contributions here are: 1) the development of a practical kinetic theory foran important class of granular time correlation functions and, 2) the demonstration of itsutility for the evaluation of Helfand and Green-Kubo expressions for Navier-Stokes ordertransport coefficients.The linear kinetic theory is summarized by Eq. (50). Based on corresponding studies forthe elastic limit of this equation, it is expected to have a wide domain of validity with respectto space and time scales, as well as densities. The nature of the approximation, short timefunctional relationship, does not explicitly entail questions of inelasticity so it is expectedto apply as well for a finite range of inelasticity. It encompasses the granular Boltzmannequation, in the low density limit, and the familiar Enskog equation in the elastic limit.The focus here has been on hydrodynamic response, but the theory includes hydrodynamicsbeyond the Navier-Stokes approximation, and even describes very short wavelength non-hydrodynamic behavior, that can be more important at moderate and high densities forgranular fluids. Finally, this kinetic theory applies beyond the set of hydrodynamic fieldsconsidered here. For any observable, z (Γ; r ), that can be written as a sum of single particlefunctions so that e z (Γ , k ) = N X r =1 e i k · q r z ( v r ) , (95)and for the same initial perturbation as considered in Sec. II, the response functions aregiven in the appropriate units by e C γ ( k ; s ) = 1 n h Z d v z ( v ) e − s e Λ( k ) e ψ (1) γ ( v , k ) . (96)This opens the possibility to study a wide range of experimental probes and also fundamentalquestions such as the relationship between fluctuations and response.23he second contribution, evaluation of the formal representations for Navier-Stokes trans-port coefficients, begins the process of exploring the utility of such formal representations,as well as verifying their consistency with earlier Chapman-Enskog based studies. These arelong wavelength properties of the kinetic theory, and therefore a more controlled context forits tests. For example, the Markovian approximation provides a practical context for the in-troduction of velocity correlations associated with the reference homogeneous state [25, 26].Their effect on transport coefficients at strong dissipation is expected to be important buthas not been quantified to date. The kinetic theory scheme developed here also provides thebasis for formulating and assessing more complex theories, such as those describing modecoupling dynamical correlations, which are expected to dominate at very high densities. Fi-nally, the verification of the agreement between the kinetic theory evaluation of the Helfandand Green-Kubo representations here and the earlier Chapman-Enskog method providesfurther support for the implicit assumptions of these complementary formal approaches.An interesting new result, both for normal and granular fluids, is the expression of theGreen-Kubo time-correlation function Ω G ( s ). In the Enskog approximation, the correspond-ing function for the shear viscosity is given by Eq. (91). To interpret this result, the timedependence may be estimated from a leading order cumulant expansion,Ω ηEG ( s ) ≃ Ω ηEG (0) e − s/τ , (97)1 τ = P di,j R d v h ij ( v ) (cid:0) J + ζ (cid:1) γ ij ( v ) P di,j R d v h ij γ ij ( v ) . (98)The corresponding Helfand correlation functions, Ω H ( s ), can be inferred directly from this,Ω ηEH ( s ) ≃ Ω ηEH (0) + Ω ηEG (0) τ (cid:0) − e − s/τ (cid:1) , (99)and the shear viscosity can be identified as η ≃ Ω ηEH (0) + Ω ηEG (0) τ. (100)These results expose the qualitative nature of the time dependence in each case. The re-sulting shear viscosity in these approximations agrees with that obtained by a leading ordersolution to the integral equation (94) as an expansion in Sonine polynomials [11].In conclusion, it is hoped that this work provides a starting point to explore systematicanalytic approximations to the hydrodynamic response of a granular fluid, with the sameattention to detail given in the context of normal fluids. These, together with numericalstudies of exact results, provide a means to understand transport mechanisms in this system.24 X. ACKNOWLEDGEMENTS
The research of A.B. and J.D. was supported in part by the Department of Energy Grant(DE-FG03-98DP00218). A.B. also acknowledges a McGinty Dissertation Fellowship and aIFT Michael J Harris Fellowship from the University of Florida. The research of J.J.B. waspartially supported by the Ministerio de Educaci´on y Ciencia (Spain) through Grant No.BFM2005-01398.
APPENDIX A: LOCAL HOMOGENEOUS STATE
The local
HCS ensemble chosen as the initial perturbation of the HCS, represents a systemdecomposed into spatial cells, each in a HCS with its own local temperature, density, andflow velocity. It is constructed formally as follows. First, the HCS ensemble is determined asthe solution to the homogeneous, stationary Liouville equation (17) in dimensionless form, ρ ∗ h (Γ ∗ ) = ρ ∗ h (cid:18)(cid:26) q rs ℓ , v r − U h v ( T h ) ; r, s = 1 , . . . , N (cid:27)(cid:19) , (A1)where U h , T h , and n h (not shown explicitly) are the flow field, temperature, and particlenumber density characterizing the HCS. Next, a conservative external force is added to Eq.(12) keeping the same L ∗ operator, ( L ∗ − N X r =1 (cid:20) ∂∂ q ∗ r φ ∗ ext ( q r ) (cid:21) · ∂∂ v ∗ r ) ρ ∗′ h = 0 . (A2)Here φ ∗ ext ≡ φ ext / T h , with φ ext ( r ) being the potential associated with the external force.The solution of Eq. (A2) is, therefore, a function of this potential, ρ ∗′ h = ρ ∗′ h (cid:18)(cid:26) q rs ℓ , v r − U h v ( t ) , φ ext ( q r ) T h ; r, s = 1 , . . . , N (cid:27)(cid:19) . (A3)This can be considered as the nonuniform fluid ensemble corresponding to the uniform limit ρ ∗ h , since in general the density will be nonuniform through its functional dependence on φ ext ( r ), n = n [ r | φ ext ] , (A4) φ ext = φ ext [ r | n ] . (A5)The second equality assumes the functional dependence of the density on the external po-tential is invertible so that the potential can be expressed as a functional of the density field.25or normal fluids in the equilibrium Gibbs state, density functional theory assures that thisis the case. In particular, for any chosen density field there is a unique external potentialcreating that field from the uniform state. It will be assumed that these properties holdas well here for the granular fluid, so that Eq. (A3) can be expressed in terms of the localdensity instead of the potential, ρ ∗′ h = ρ ∗′′ h (cid:18)(cid:26) q rs ℓ , v r − U h v ( T h ) , n ( q r ) ℓ d ; r, s = 1 , . . . , N (cid:27)(cid:19) . (A6)With ρ ∗′′ h known from the solution to Eq. (A2 ), the local HCS is constructed by the replace-ments v r − U h v ( T h ) → v r − U h − δ U ( q r ) v [ T h + δT ( q r )] ,n ( q r ) → n h + δn ( q r ) , (A7)to get ρ ∗ ℓh [Γ ∗ | δy ∗ ] ≡ ρ ∗′′ h (cid:18)(cid:26) q rs ℓ , v r − U h − δ U ( q r ) v [ T h + δT ( q r )] , [ n h + δn ( q r )] ℓ d ; r, s = 1 , . . . , N (cid:27)(cid:19) . (A8)Note that the local HCS is no longer a solution to any Liouville equation, but rather issimply a reference ensemble representing an hypothetical HCS with different hydrodynamicparameters in each spatial cell of the fluid. Its construction in the way presented above,supports that interpretation in the sense that ρ ∗ ℓh [Γ ∗ | δy ∗ = constant] = ρ ∗ h (Γ ∗ ; y ∗ h + δy ∗ )so that both { ρ ∗ ℓh [Γ ∗ | δy ∗ ] } δy ∗ =0 = ρ ∗ h (Γ ∗ ; y ∗ h ) and all functional derivatives of ρ ∗ ℓh [Γ ∗ | δy ∗ ]become derivatives of ρ ∗ h (Γ ∗ ; y ∗ h ) at δy ∗ = 0. More explicitly, it is Z d r . . . Z d r p (cid:20) δ p ρ lh [Γ | δy ] δy α ( r ) . . . δy β ( r p ) (cid:21) y = { n h ,T h , } = (cid:20) ∂ p ρ h (Γ; n h , T h , U h ) ∂y α,h . . . ∂y β,h (cid:21) U h = , (A9)where Γ is a point in the phase space associated to the original positions and velocitiesand ρ (Γ) the corresponding density. It is instructive to carry out this construction of thelocal ensemble for the case of a normal fluid. Then Eq. (A1) gives the familiar equilibriumGibbs ensemble, and Eq. (A2) gives the same ensemble with the Hamiltonian modified toinclude the external potential. Finally, the construction in Eq. (A8) gives the familiarlocal equilibrium ensemble used in linear response theory for spatial perturbations of theequilibrium state. In the grand ensemble the dependence on the local density is implicitthrough a local chemical potential µ = µ [ r | n ].26 PPENDIX B: TWO PARTICLE FUNCTIONAL
The reduced distribution functions associated with the solution to the Liouville equation(20) are defined as f ( m ) ( x , . . . , x m , s ) ≡ N !( N − m )! Z dx m +1 . . . Z dx N e − s L ρ ℓh [Γ | δy ] . (B1)As is done in the main text, the asterisk indicating the use of dimensionless variables is leftimplicit. The f ( m ) ’s are clearly not independent functions. For example, f (2) ( x , x , s ) isrelated to f (1) ( x , s ) by ( N − f (1) ( x , s ) = Z dx f (2) ( x , x , s ) . (B2)This implies that f (2) ( x , x , s ) has the representation f (2) ( x , x , s ) = f (1) ( x , s ) f (1) ( x , s ) g (2) ( x , x , s ) , (B3)where the pair correlation function g (2) ( x , x , s ) has the properties g (2) ( x , x , s ) = g (2) ( x , x , s ) , (B4) Z dx f (1) ( x , s ) g (2) ( x , x , s ) = Z dx f (1) ( x , s ) g (2) ( x , x , s ) = N − . (B5)This in turn shows that g (2) ( x , x , s ) is a functional of f (1) ( x , s ). Quite generally then, f (2) ( x , x , s ) can be considered a functional of f (1) ( x , s ), f (2) ( x , x , s ) = f (2) (cid:2) x , x , s | f (1) ( s ) (cid:3) . (B6)However, this functional relationship is not unique. The utility of Eq. (B6) lies in discoveringa choice that leaves the simplest functional dependence. In the low density limit, where itis expected that g (2) ( x , x , s ) tends to unity, this is clearly the case, since the functionalbecomes a constant independent of x , x , and s . More generally, finding an appropriatefunctional form for g (2) ( x , x , s ) requires the detailed analysis of the full many-body prob-lem. 27he corresponding functional for ψ (2) γ ( x , x , s ) defined by Eq. (30) can be computed as ψ (2) γ ( x , x , s ) ≡ (cid:20) δf (2) ( x , x , s ) δy γ ( , (cid:21) δy =0 = Z dx " δf (2) (cid:2) x , x , s | f (1) ( s ) (cid:3) δf (1) ( x, s ) δf (1) ( x, s ) δy γ ( , δy =0 = Z dx (cid:20) δf (2) [ x , x , s | f (1) ( s )] δf (1) ( x, s ) (cid:21) δy =0 ψ (1) γ ( x, s )= Z dx K ( x , x , s ; x ) ψ (1) γ ( x, s ) , (B7)with K ( x , x , s ; x ) given by Eq. (35). This entails the additional requirement thatthe functional f (2) [ x , x , s |· ] is independent of the specific initial fields δy ( r , f (2) (cid:2) x , x , s | f (1) ( s ) (cid:3) depends on these fields only through f (1) ( x, s ). Consequently, K ( x , x , s ; x ) is also independent of such initial data. This is expected, since the collisionoperator M ( s ) constructed from K ( x , x , s ; x ) by means of Eq. (37) should be universal fora wide class of initial conditions. Finally, the functional form for ψ (2) γ ( x , x , s ) is seen to belinear in ψ (1) γ ( s ), while in general f (2) [ x , x , s | f (1) ( s )] is a nonlinear functional of f (1) ( s ). APPENDIX C: INTERPRETATION OF e c λ ( v , k ) The contribution to the action of the collision operator f M ( k ) on e ψ (1) γ ( v , − k , s ) from theterm proportional to e c λ ( v , k ) in Eq. (52) is − X λ e c λ ( v , k ) 1 n h Z d v a λ ( v ) e ψ (1) γ ( v , − k , s ) = − X λ e c λ ( v , k ) e C λγ ( k ; s ) , (C1)where Eq. (28) has been employed. The above contribution depends on only low ordermoments of the dependent variable e ψ (1) γ in the kinetic equation (50), and in fact only thosemoments are of interest for determining the response functions. In this sense, e c λ is a meanfield operator rather than a true collision operator, since its action does not depend directlyon differences in e ψ (1) γ before and after a collision like the first term of (52). Instead, Eq.(49) shows that e c λ reflects an average of collisional effects induced through changes in thecorrelations. To provide some interpretation of this term, consider first the elastic limit.28 . Elastic limit In this case, g (2) ℓh [ x , x | δy ] = g (2) ℓe [ q , q | δn ], independent of the velocities, temperature,and flow field, i.e. it is a function of the spacial coordinates and a functional of the density.Equation (49) becomes c λ ( v , q ) = δ λ n h Z dx T ( x , x ) ϕ ( v ) ϕ ( v ) " δg (2) ℓe [ q , q | δn ] δn ( q ) δn =0 , (C2)where ϕ ( v ) is the Maxwellian. The pair correlation function for a nonuniform fluid, g (2) ℓe [ q , q | δn ], appears in the stationary first BBGKY hierarchy equation ( 32) in the pres-ence of an external potential φ ext associated with the given density (see Appendix A), (cid:20) v · ∂∂ q − (cid:18) ∂∂ q φ ext [ q | n ] (cid:19) · ∂∂ v (cid:21) n ( q ) ϕ ( v )= Z dx T ( x , x ) ϕ ( v ) ϕ ( v ) n ( q ) n ( q ) g (2) ℓe [ q , q | δn ] . (C3)The functional derivative of this equation with respect to δn ( q ) evaluated at at δn = 0gives n h ϕ ( v ) v · ∂∂ q C ( q ) = n h Z dx T ( x , x ) ϕ ( v ) ϕ ( v ) g (2) e ( q ) [ δ ( q ) + δ ( q )]+ n h Z dx T ( x , x ) ϕ ( v ) ϕ ( v ) " δg (2) ℓe [ q , q | δn ] δn ( q ) δn =0 . (C4)where it has been used that g e ( q ) = g ℓe [ q , q |
0] and c ( q ) is the direct correlation functiondefined by [27] n h c ( q − q ) = δ ( q − q ) + 2 n h (cid:20) δφ ext [ q | n ] δn ( q ) (cid:21) φ ext =0 . (C5)The first term on the right hand side of Eq. (C4) can be evaluated using the elastic limit ofthe explicit form for T ( x , x ) given in Eq. (13). The result is Z dx T ( x , x ) ϕ ( v ) ϕ ( v ) g (2) e ( q ) [ δ ( q ) + δ ( q )] = g (2) e ( σ ) ϕ ( v ) δ ( q − σ ) b q · v (C6)Finally, therefore, Eq. (C4 gives n h Z dx T ( x , x ) ϕ ( v ) ϕ ( v ) " δg (2) ℓe [ q , q | δn ] δn ( q ) δn =0 ϕ ( v ) v · ∂∂ q (cid:2) c ( q ) − g (2) e ( σ )Θ( q − σ ) (cid:3) , (C7)where the delta function in Eq. (C5) has been written in terms of the derivative of Θ( q − σ ).Using this into Eq. (C2) gives the desired result, c λ ( v , q ) = δ λ n h ϕ ( v ) v · ∂∂ q (cid:2) c ( q ) − g (2) e ( σ )Θ( q − σ ) (cid:3) . (C8)In the elastic case, it is seen that the contribution shown in Eq. (C1) is the same as thatfor an external force whose potential is c ( q ) − g (2) h ( σ ) Θ ( q − σ ). The direct correlationfunction has a discontinuity at q = σ , with value c ( σ ) = g (2) e ( σ ), so the subtracted thetafunction contribution assures that this potential is continuous.
2. Inelastic collisions
For inelastic collisions, the effects of e c λ ( v , k ) are more complex and more difficult tointerpret. However, a significant difference from the elastic case can be seen already for thesimplest case of k = . From Eq. (48), it follows that e c λ ( v , ) = Z d q c λ ( x , q ) = Z dx T ( x , x ) f (1) h ( v ) f (1) h ( v ) " ∂g (2) h ( x , x ; n h , T h , U h ) ∂y λ,h U h = , (C9)where use has been made of the identity Z d q " δg (2) ℓh [ x , x | δy ] δy λ ( q ) δy =0 = " ∂g (2) h ( x , x ; y ) ∂y λ { n,T, U } = { n h ,T h , } . (C10)This k = limit vanishes for elastic collisions, as it can be seen directly from Eq. (C8), butis nonzero for inelastic collisions. This can be verified in the Enskog approximation, wherevelocity correlations are neglected and, therefore, e c Eλ ( v , ) = δ λ ∂g (2) h ( σ ; n ) ∂n ! n = n h Z dx T ( x , x ) f (1) h ( v ) f (1) h ( v )= δ λ ∂ ln g (2) h ( σ ; n ) ∂n ! n = n h ζ ∂∂ v · h v f (1) h ( v ) i . (C11)The second equality follows from the first hierarchy equation (32) particularized for theHCS, ζ ∂∂ v · h v f (1) h ( v ) i = Z dxT ( x , x ) g (2) h ( x , x ) f (1) h ( v ) f (1) h ( V ) , (C12)30hen velocity correlations are neglected. Thus, it is seen that e c λ ( v , , ) includes changes inthe correlations of collisional effects associated with cooling.More generally, the Enskog approximation for arbitrary k reads e c Eλ ( v , k ) = δ λ Z d q e i k · q c λ ( x , q )= Z dx T ( x , x ) f (1) h ( v ) f (1) h ( v ) Z d q e i k · q " δg (2) ℓh [ q , q | δy ] δn ( q ) δn =0 . (C13) APPENDIX D: HYDRODYNAMIC MODES OF e Λ( ) In this Appendix, some of the details leading to the solution of the eigenvalue problem(60) at k = 0 are given. Consider first the functions e ψ (1) γ ( v , k ), defined in Eq. (54), at k = 0. By translational invariance " δf (1) ℓh [ q , v | δy ] δy γ ( ) δy =0 = " δf (1) ℓh [ q + r , v | δy ] δy γ ( r ) δy =0 , (D1)and so e ψ (1) γ ( v , ) = Z d q " δf (1) ℓh [ x | δy ] δy γ ( ) δy =0 = 1 V Z d q Z d r " δf (1) ℓh [ x | δy ] δy γ ( r ) δy =0 = " ∂f (1) h ( v ; y ) ∂y γ y = y h . (D2)The last equality is a consequence of the construction of the local HCS, assuring that allfunctional derivatives in the homogeneous limit are related with ordinary derivatives of theHCS (see Appendix A). Use of the expression of the operator e Λ, Eq. (51), yields e Λ( ) e ψ (1) γ ( v , ) = f M ( ) " ∂f (1) h ( v ; y ) ∂y γ y = y h + ζ ∂∂ v · v " ∂f (1) h ( v ; y ) ∂y γ y = y h . (D3)31ext, Eq. (52) gives f M ( ) " ∂f (1) h ( v ; y ) ∂y γ y = y h = − Z dx T ( x , x ) g (2) h ( x , x ) (cid:26) ∂∂y γ h f (1) h ( v ; y ) f (1) h ( v ; y ) i(cid:27) y = y h − n h X λ Z dx T ( x , x ) f (1) h ( v ) f (1) h ( v ) " ∂g (2) h ( x , x ; δy ) ∂y λ y = y h × Z d v a λ ( v ) " ∂f (1) h ( v ; y ) ∂y γ y = y h = − (cid:20) ∂∂y γ Z dx T ( x , x ) g (2) h ( x , x ; y ) f (1) h ( v ; y ) f (1) h ( v ; y ) (cid:21) y = y h , (D4)where Eq. (C9) has been employed. The right hand side in the above equation can be furthersimplified by means of Eq. (C12), after changing v r into v r − U , f M ( ) " ∂f (1) h ( v ; y ) ∂y γ y = y h = − (cid:26) ∂∂y γ (cid:20) ζ ( y )2 ∂∂ v · ( v − U ) f (1) h ( v , y ) (cid:21)(cid:27) y = y h , (D5)and Eq. (D3) becomes e Λ( ) e ψ (1) γ ( v , ) = − (cid:20) ∂ζ ( y ) ∂y γ (cid:21) y = y h ∂∂ v · h v f (1) h ( v ) i + ζ (cid:18) ∂ U ∂y γ (cid:19) y = y h · ∂∂ v f (1) h ( v ) . (D6)Realizing that " ∂f (1) h ( v , y ) ∂y y = y h = − ∂∂ v · h v f (1) h ( v ) i (D7)and " ∂f (1) h ( v , y ) ∂ y y = y h = − ∂f (1) h ( v ) ∂ v , (D8)Eq. (D6) is seen to be equivalent to e Λ( ) e ψ (1) γ ( v , ) = (cid:20) ∂ζ ( y ) ∂y γ (cid:21) y = y h e ψ (1)2 ( v , ) − ζ (cid:18) ∂ U ∂y γ (cid:19) y = y h · e ψ (1)3 ( v , ) . (D9)Above, the index 3 and vectorial notation has been used to identify the d componentsassociated to the velocity field.The specific cases following from Eq. (D9) are e Λ( ) e ψ (1)1 ( v , ) = n h (cid:18) ∂ζ ∂n h (cid:19) T h e ψ (1)2 ( v , ) , (D10)32 Λ( ) e ψ (1)2 ( v , ) = ζ e ψ (1)2 ( v , ) , (D11)and e Λ( ) e ψ (1)3 ( v , ) = − ζ e ψ (1)3 ( v , ) . (D12)Equations (D10) and (D11) can be combined to give e Λ( ) " e ψ (1)1 ( v , ) − (cid:18) ∂ ln ζ ∂ ln n h (cid:19) T h e ψ (1)2 ( v , ) = 0 . (D13)The above equations (D11)-(D13) have the form of eigenvalue equations e Λ( ) φ ( β ) ( , v ) = λ ( β ) ( ) φ ( β ) ( , v ) , (D14) β = 1 , . . . , d + 2, with the eigenvalues and eigenfunctions given in Eqs. (61) and (62). APPENDIX E: TRANSPORT COEFFICIENTS
In this appendix, the pressure, the cooling rate, and the Navier-Stokes transport coeffi-cients are evaluated in the kinetic theory approximations developed in the main text. Thesequantities are identified from the linearized phenomenological Navier-Stokes equations foran isolated granular fluid, that in dimensionless units have the form ∂∂s δ e y ∗ β ( k ∗ , s ) + d +2 X γ =1 K ∗ hydβγ δy ∗ γ ( k ∗ , s ) = 0 , (E1)where the transport matrix K ∗ hyd is found to be block diagonal with a “longitudinal” part,corresponding to the fields { δn ∗ , δT ∗ , δU ∗⊥ } , given by K ∗ hyd ( k ∗ ) = − ik ∗ ζ ∗ ∂ ln ζ ∂ ln n h + (cid:0) µ ∗ d − ζ ∗ n (cid:1) k ∗ ζ ∗ + (cid:0) λ ∗ d − ζ ∗ T (cid:1) k ∗ − i (cid:16) p ∗ h d + ζ ∗ U (cid:17) k ∗ − ip ∗ h ∂ ln p h ∂ ln n h k ∗ − ip ∗ h k ∗ − ζ ∗ + h d − d η ∗ + κ ∗ i k ∗ . (E2)The “transverse” components, δU ∗⊥ ,i , decouple from the longitudinal ones, and their trans-port matrix reads K hyd ∗ ( k ) = (cid:18) − ζ ∗ η ∗ k ∗ (cid:19) I, (E3)with I being the unit matrix of dimesnsion d −
1. These expressions include the unspecifiedfunctions p h ( n h , T h ) and ζ h ( n h , T h ), defining the pressure and the cooling rate of the HCS,33espectively, as well as the unknown shear viscosity η ( n h , T h ), the bulk viscosity κ ( n h , T h ),the thermal conductivity λ ( n h , T h ), and the new coefficient, µ ( n h , T h ), associated with thecontribution of the density gradient to the heat transport in a granular fluid. Finally, ζ n ( n h , T h ), ζ T ( n h , T h ), and ζ U ( n h , T h ) are transport coefficients arising from the local coolingrate. The dimensionless form of the cooling rate was defined in Eq. (19), and the definitionsfor the remaining dimensionless quantities are p ∗ h = p h n h T h , η ∗ = ηmn h lv , κ ∗ = κmn h lv ,λ ∗ = λln h v , µ ∗ = µlT h v ,ζ ∗ U = ζ U , ζ ∗ n = n h ζ n lv , ζ ∗ T = T h ζ T lv . (E4)Formal expressions for these parameters have been obtained in ref. [9], by carrying outa linear response analysis. Each of those expressions is evaluated below by means of thekinetic theory developed in the text, for comparison to the results obtained by a Chapman-Enskog solution to the nonlinear Enskog equation [11]. First, some general considerationsthat apply to all of the transport coefficients will be addressed. Dimensionless units will beused in the remaining of this Appendix, and the asterisk will be suppressed for simplicity,as done in the main text.
1. Reduced Form of the Transport Coefficients
The Helfand form for a generic dimensionless transport coefficient χ as obtained in ref.[9] is χ = lim Ω H ( s ) = lim 1 V Z d Γ j (Γ) (1 − P ) e − s ( L− λ ) M (Γ) , (E5)where j (Γ) is a flux associated with the densities of the hydrodynamic variables, and λ denotes one of the eigenvalues in (61). These fluxes have the generic form j (Γ) = N X r =1 j ( x l ) + N X r =1 N X s = r j ( x r , x s ) , (E6) j and j being one and two-particle functions of the phase point, respectively. The adjointfunctions M (Γ) are related to the conjugate densities in the hydrodynamic response definedin Eq. (27), and read M (Γ) = M Z d r b k · r (cid:20) δρ lh [Γ | y ] δy ( r ) (cid:21) y = y h , (E7)34 being some constant. Moreover, P is the projection operator defined by P (Γ) X (Γ) = 1 V d +2 X γ =1 e ψ γ (Γ; ) Z d Γ ′ e a γ (Γ ′ ; ) X (Γ ′ ) , (E8)where the phase functions e a γ (Γ; ) and e ψ γ (Γ; ) are the densities defined by Eqs. (22) andthe Fourier transform of (27), respectively, both evaluated at k = . The generator for thedynamics is L− λ , where λ is one of the hydrodynamic modes at k = identified in Eq. (61).Both the projection operator and the additional time dependence of the term containing λ are necessary to ensure that the long time limit of the correlation function in Eq. (E5) iswell defined.Proceeding as in Sec. VII for the shear viscosity, a reduced expression for these transportcoefficients in the Markov approximation can be obtained. The generic correlation functionΩ H ( s ) of Eq. (E5) becomesΩ H ( s ) ≃ V Z dx J ( x ) (cid:2) − P (1) ( x ) (cid:3) e − s (Λ − λ ) M (1) ( x ) . (E9)The direct flux J ( x ) in this reduced time correlation function is J ( x ) ≡ j ( x ) + Z dx Z dx j ( x , x ) K ( x , x , x ) , (E10)where K ( x , x , x ) is the kernel for the collision operator given by Eq. (44), and M (1) isthe one-particle function in the hierarchy associated with MM (1) ( x ) = N Z dx . . . Z dx N M (Γ) . (E11)The generator for the Markov dynamics Λ is given in (42). Lastly, P (1) is the one particleanalog of P defined in Eq. (E8), P (1) ( x ) X ( x ) ≡ V d +2 X γ =1 e ψ (1) γ ( v , ) 1 n h Z dx a γ ( v ) X ( x ) , (E12)where the e ψ (1) γ ( v , ) are the functions defined in Eq. (D2).Equation (E5) can be transformed into the Green-Kubo form to give χ = Ω H (0) + lim Z s ds ′ ∂∂s ′ Ω H ( s ′ ) ≡ Ω H (0) + lim Z s ds ′ Ω G ( s ′ ) . (E13)The Markov approximation in this representation follows directly from Eq. (E9),Ω H (0) ≃ V Z dx J ( x ) (cid:2) − P (1) ( x ) (cid:3) M (1) ( x ) , (E14)35 G ( s ) ≃ V Z dx J ( x ) (cid:2) − P (1) ( x ) (cid:3) e − s (Λ − λ ) γ ( x ) , (E15)with the conjugate flux γ ( x ) given by γ ( x ) ≡ − (Λ − λ ) M (1) ( x ) . (E16)In this Green-Kubo form, the role of the projection operator can be readily interpreted asfollows. Taking into account that P (1) projects over the subspace spanned by the hydrody-namic eigenfunctions of e Λ( ), the property (cid:2) − P (1) ( x ) (cid:3) e − s (Λ − λ ) = (cid:2) − P (1) ( x ) (cid:3) e − s (Λ − λ ) (cid:2) − P (1) ( x ) (cid:3) , (E17)is obtained. This shows that the presence of the projection operator in Eq. (E15) ensuresthat the generator of dynamics exp[ − s (Λ − λ )] acts on a function that is orthogonal to itsinvariants and, therefore, has a well defined long time limit.The above expressions can be specialized to the granular Enskog approximation by replac-ing the Markovian kernel K in Eq. (E10) with its Enskog approximation, and by replacingthe linear operator Λ by Λ E , defined in Eq. (85). Further, the reduced conjugate functionsin (E11) and Eqs. (E12) in the Enskog approximation become (cid:8) M (1) γ ( x ) (cid:9) ≃ (cid:26) M b k · q f (1) h ( v ) , M b k · q ∂∂ v · h v f (1) h ( v ) i , − M b k · q ∂∂ v f (1) h ( v ) (cid:27) (E18)and (cid:8) ψ (1) γ ( v , ) (cid:9) ≃ (cid:26) f (1) h ( v ) , ∂∂ v · h v f (1) h ( v ) i , − ∂∂ v f (1) h ( v ) (cid:27) , (E19)respectively. Now, f (1) h ( v ) is the one-particle distribution of the HCS obtained from thenonlinear Enskog kinetic equation. In the rest of the appendix, explicit expressions for thehydrodynamic parameters are given in the granular Enskog approximation.
2. Evaluation in the Enskog Approximation
The Helfand and Green-Kubo expressions for the Navier-Stokes transport coefficients ofa hard sphere granular fluid are reported in ref. [9]. From those expressions, the N particlefunctions j (Γ) and M (Γ), together with the eigenvalue λ appearing in Eq. (E5) for each ofthe transport coefficients, can be read off. Then, following the procedure illustrated aboveand, in the context of the shear viscosity in the main text, the results reported below are36btained. Attention is restricted to the transport coefficients that were calculated in ref. [11],which excludes the transport coefficients ζ T and ζ n from the local cooling rate. Moreover,the parameter ℓ defining the length scale is chosen such that n h ℓ = 1 in the following, forthe sake of simplicity. a. The pressure The expression for the pressure has been identified in [9] as p = 1 + (1 + α ) σ V d Z dx Z dx δ ( q − σ )Θ ( − b q · g ) ( b q · g ) f (2) h ( q , v , v ) . (E20)This is the second moment of the normal component of the relative velocity averaged over thetwo-particle distribution at contact. In the Enskog approximation, all velocity correlationsin the two-particle distribution are neglected, and the above expression simplifies to p = 1 + π d/ (1 + α ) σ d d/ d g (2) h ( σ ) . (E21) b. The zeroth order cooling rate ζ The homogeneous dynamics in Eq. (E1) is determined entirely determined by ζ , forwhich the expression ζ = 1 − α V d Z dx Z dx δ ( q − σ )Θ ( b q · g ) | b q · g | f (2) h ( q , v , v ) (E22)was derived in [9]. When the Enskog approximation for the two-particle distribution functionis used, this simplifies to ζ = π ( d − / σ d − g (2) h ( σ )(1 − α )2Γ (cid:0) d +32 (cid:1) d Z d v Z d v g f (1) h ( v ) f (1) h ( v )= Γ ( d/
2) (1 − α ) π / Γ (cid:0) d +32 (cid:1) σ ( p − Z d v Z d v g f (1) h ( v ) f (1) h ( v ) . (E23) c. The Euler transport coefficient ζ U This is a new transport coefficient unique to granular fluids. Its expression, as obtainedin ref. [9], reads ζ U = lim Ω ζ U H ( s ) = Ω ζ U H (0) + lim Z s ds ′ Ω ζ U G ( s ′ ) , (E24)37here Ω ζ U H ( s ) = V − Z d Γ W (Γ) (1 − P ) e − s ( L + ζ ) M ζ U (Γ) (E25)and Ω ζ U G ( s ) = V − Z d Γ W (Γ) (1 − P ) e − s ( L + ζ )Υ ζ U (Γ) . (E26)In the above expressions, W (Γ) is the source term in the microscopic energy balance equationthat characterizes the dissipation due to the inelastic collisions, W (Γ) = 1 − α d N X r =1 N X s = r δ ( q rs − σ )Θ( − b q rs · g rs ) | b q rs · g rs | , (E27)and the conjugate density and flux are M ζ U (Γ) = − N X r q r · ∂∂ v r ρ h (Γ) , (E28)Υ ζ U (Γ) = − (cid:18) L + ζ (cid:19) M ζ U (Γ) . (E29)The Helfand form for this transport coefficient in the Enskog approximation can be obtainedby application of Eq. (E9) with the resultΩ ζ U H ( s ) ≃ V − Z dx J ζ U ( x ) (cid:2) − P (1) ( x ) (cid:3) e − s ( Λ E + ζ ) M (1) ζ U ( x ) , (E30)where J ζ U ( v ) = π d − g (2) h ( σ ) σ d − Γ (cid:0) d +32 (cid:1) d Z d v | v − v | f (1) h ( v ) , (E31) M (1) ζ U ( x ) = − q · ∂∂ v e ψ (1)1 ( v , ) , (E32)and Λ E is the collision operator given in Eq. (85). The Green-Kubo form of this transportcoefficient is determined by Ω ζ U H (0) and Ω ζ U G ( s ). Direct evaluation of the former gives [9]Ω ζ U H (0) = − d (1 − α ) ( p − . (E33)Here p is the pressure in the Enksog approximation, Eq. (E21). The time correlation functionis (see Eq. (E15))Ω ζUG ( s ) = V − Z dx J ζ U ( x ) (cid:2) − P (1) ( x ) (cid:3) e − s ( Λ E + ζ ) γ ζ U ( x ) (E34)with the conjugate flux given by γ ζ U ( x ) = − (cid:18) Λ E + ζ (cid:19) M (1) ζ U ( x ) . (E35)38roceeding as in the case of the shear viscosity, it can be shown that e − s ( Λ E + ζ ) γ ζ U ( v ) = e − s ( J + ζ ) γ ζ U ( v ) , (E36)where J is the operator defined in Eq. (89). The time correlation function in Eq. (E34) canbe expressed as Ω ζ U G ( s ) = Z d v J ζ U ( v ) (cid:2) − P (1) ( x ) (cid:3) e − s ( J + ζ ) γ ζ U ( v ) . (E37)Therefore, this Euler order transport coefficient is given in the Enskog approximation by ζ U = − d (1 − α ) ( p −
1) + Z d v J ζ U ( v ) (cid:2) − P (1) ( x ) (cid:3) D ζ U ( v ) , (E38)where D ζ U ( v ) is the solution to the integral equation (cid:18) J + ζ (cid:19) D ζ U ( v ) = (cid:2) − P (1) ( x ) (cid:3) γ ζ U ( v ) . (E39) d. The bulk viscosity κ The expression for the bulk viscosity reported in ref. [9] is κ = lim Ω κH ( s ) = Ω κH (0) + lim Z s ds ′ Ω κG ( s ′ ) , (E40)where Ω κH ( s ) = − V d Z d Γ tr H (1 − P ) e − s ( L + ζ ) M κ (Γ) , (E41)and Ω κG ( s ) = − V d Z d Γ tr H (1 − P ) e − s ( L + ζ )Υ κ (Γ) . (E42)The direct flux tr H is the trace of the momentum flux given in Eq. (66), and the adjointdensity is the same as that for ζ U , i.e., M κ (Γ) = M ζ U (Γ) = − N X r =1 q r · ∂∂ v r ρ h (Γ) . (E43)Consequently, the conjugate flux is also the same, Υ κ = Υ ζ U , given by (E29). In the Enskogapproximation, the time correlation function Ω κH ( s ) becomesΩ κH ( s ) ≃ V − Z dx J κ ( x ) (cid:2) − P (1) ( x ) (cid:3) e − s ( Λ E + ζ ) M (1) κ ( x ) , (E44)39ith the direct flux given by J κ ( x ) = − v d + 12 d ∂ ln g (2) h ( σ ) ∂n h ( p − − (1 + α ) σg (2) h ( σ )2 d Z dx δ ( q − σ ) Θ ( − b q · g ) | b q · g | f (1) h ( v ) (E45)and M (1) κ ( x ) = M (1) ζ U ( x ), given in Eq. (E32). The instantaneous contribution Ω κH (0) is simplyrelated to Ω ηH (0) in (93) through [9]Ω κH (0) = d + 2 d Ω ηH (0) . (E46)Furthermore, the correlation function Ω κG ( s ) vanishes. This can be seen as follows. Theconjugate flux and the generator of dynamics can be simplified as in the case of ζ U above,so that they become independent of the coordinate q . Then, the direct flux above can besimplified to give J κ ( v ) = − " α ) σ d π d/ g (2) h ( σ )2Γ (cid:0) d (cid:1) d v d − " ∂ ln g (2) h ( σ ) ∂n h + 1 p − d . (E47)Since J κ ( v ) is orthogonal to the subspace spanned by 1 −P (1) ( x ), Ω κG ( s ) vanishes. Therefore,the bulk viscosity in the Enskog approximation is simply κ = d + 2 d Ω ηH (0) = π ( d − / (1 + α ) σ d +1 g (2) h ( σ )2 d Γ (cid:0) d +12 (cid:1) Z d v Z d v | v − v | f (1) h ( v ) f (1) h ( v )= π − / Γ ( d/ σ ( p − (cid:0) d +12 (cid:1) d Z d v Z d v | v − v | f (1) h ( v ) f (1) h ( v ) . (E48) e. The thermal conductivity λ The thermal conductivity is expressed in ref. [9] as λ = lim Ω λH ( s ) = Ω λH (0) + lim Z s ds ′ Ω λG ( s ′ ) , (E49)where Ω λH ( s ) = − ( V d ) − Z d Γ S (Γ) · (1 − P ) e − s ( L− ζ ) M λ (Γ) (E50)and Ω λG ( s ) = − ( V d ) − Z d Γ S (Γ) · (1 − P ) e − s ( L− ζ ) Υ λ (Γ) . (E51)40he direct flux S is the heat flux in the microscopic balance equation for the energy density, S = N X r =1 v r v r + N X r =1 N X s = r s ( x r , x s ) , (E52)with s ( x r , x s ) = (1 + α ) σ δ ( q rs − σ )Θ ( − b q rs · g rs ) ( b q rs · g rs ) ( b q rs · G rs ) b q rs , (E53)where G rs ≡ ( r + v s ) /
2. The conjugate density in Eq. (E50) is M λ (Γ) = − N X r =1 q r ∂∂ v r · [ v r ρ h (Γ)] (E54)and the associated flux is Υ λ (Γ) = − (cid:18) L − ζ (cid:19) M λ (Γ) . (E55)In the Enskog approximation, the time correlation function determining the Helfand formof the thermal conductivity λ isΩ λH ( s ) ≃ V − Z dx J λ ( x ) · (cid:2) − P (1) ( x ) (cid:3) e − s ( Λ E − ζ ) M (1) λ ( x ) , (E56)where J λ ( x ) = v v d − (1 + α ) σg (2) h ( σ ) d Z dx δ ( q − σ ) Θ ( − b q · g ) × ( b q · g ) ( b q · G ) b q f (1) h ( v ) (E57)and M (1) λ ( x ) = − q ∂∂ v · h v f (1) h ( v ) i . (E58)The function J λ defined in Eq. (E57) differs from the generic form given in (E10) by avelocity independent term which does not contribute to Ω λH . The contribution from theinitial correlation function to the Green-Kubo form of this transport coefficient is given byΩ λH (0) = Γ( d )2 π / Γ (cid:0) d (cid:1) σ ( p − Z d v Z d v f (1) h ( v ) f (2) h ( v ) × (cid:20) g G + g ( b g · G ) + 14 g + 32 ( b g · G ) (cid:21) (E59)while the time-dependent correlation function becomesΩ λG ( s ) = − (cid:20) p − d + 2 (cid:21) Z d v v v d · (cid:2) − P (1) ( v ) (cid:3) e − s ( J − ζ ) γ λ ( v ) . (E60)41roceeding as in the case of the shear viscosity in Sec. VII, it is found that γ λ = (cid:18) Λ E − ζ (cid:19) q ∂∂ v · (cid:2) v f h ( v ) (cid:3) . (E61)Thus the thermal conductivity in the Enskog approximation is given by λ = Ω λH (0) − (cid:20) p − d + 2 (cid:21) Z d v v v d · (cid:2) − P (1) ( v ) (cid:3) A ( v ) , (E62)where A is the solution to the integral equation (cid:18) J − ζ (cid:19) A ( v ) = (cid:2) − P (1) ( v ) (cid:3) γ λ ( v ) , (E63)with J the operator defined in Eq. (89). f. The coefficient µ This is a new transport mechanism for granular fluids that arises due to the inelasticityof collisions. As discussed in detail in ref. [9], this transport coefficient consists of twotime correlation functions, one of which can be recognized as the time correlation part ofthe thermal conductivity. Therefore, the quantity in terms of which the time correlationfunction takes the simplest form is the linear combination µ ≡ µ − ∂ ln ζ h ∂ ln n h λ = lim Ω µH ( s ) = Ω µH (0) + lim Z s ds ′ Ω µG ( s ′ ) , (E64)with Ω µH ( s ) = − ( V d ) − Z d Γ S (Γ) · (1 − P ) e − s L M µ (Γ) , (E65)and Ω µG ( s ) = − ( V d ) − Z d Γ S (Γ) · (1 − P ) e − s L Υ µ (Γ) . (E66)The direct flux S above is the heat flux given in Eq. (E52). The conjugate density is M µ = Z d r r (cid:20) δρ lh δn ( r ) − ∂ ln ζ ∂ ln n T δρ lh δT ( r ) (cid:21) y = y h = Z d r r "(cid:18) δρ lh [Γ | y ] δn ( r ) (cid:19) ζ y = y h , (E67)where the functional derivative with respect to the density in the last equality is to be takenat constant cooling rate, as indicated. In the Enskog approximation, the Helfand formbecomes Ω µH ( s ) ≃ V − Z dx J µ ( x ) · (cid:2) − P (1) ( x ) (cid:3) e − s ( Λ E − ζ ) M (1) µ ( x ) . (E68)42ere, J µ ( x ) = J λ ( x ), given in Eq. (E57), and M (1) µ = q (cid:26) f (1) h ( v ) + ∂ ln ζ ∂ ln n h ∂∂ v · h v f (1) h ( v ) i(cid:27) . (E69)The Green-Kubo form is determined from Ω µH (0) and Ω µG ( s ). Direct evaluation shows therelationship Ω µH (0) = − ∂ ln ζ ∂ ln n h Ω λH (0) . (E70)The time correlation function Ω µG ( s ) can be simplified as in the previous cases yieldingΩ µG ( s ) = Z d v J µ ( v ) · (cid:2) − P (1) ( v ) (cid:3) e − s J γ µ ( v ) , (E71)with γ µ,i ( v ) = − v i f (1) h ( v ) − " ∂ ln g (2) h ( σ ) ∂ ln n h Q i h f (1) h ( v ) i − ∂ ln ζ ∂ ln n h γ λ,i . (E72)The Green-Kubo form for the transport coefficient µ is µ − ∂ ln ζ ∂ ln n h λ = − ∂ ln ζ ∂ ln n h Ω λH (0) + Z d vJ λ ( v ) · (cid:2) − P (1) ( x ) (cid:3) C ( v ) , (E73)where C ( v ) is a solution to the integral equation J C ( v ) = (cid:2) − P (1) ( x ) (cid:3) J λ ( v ) . (E74)All of the above results agree in detail (for d = 3) with those obtained in ref. [11] bymeans of the Chapman-Enskog solution of the non-linear Enskog kinetic equation. [1] See, for example, P.C. Martin in Many Body Physics , edited by C. de Witt and R. Balian(Gordon and Breach, New York, 1968).[2] D. Forster,
Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions (Ben-jamin, Reading MA, 1975).[3] J-P. Boon and S. Yip,
Molecular Hydrodynamics (Dover, New York, 1991).[4] J.A. McLennan,
Introduction to Nonequilibrium Statistical Mechanics (Prentice-Hall, NewJersey, 1989).[5] E. Helfand, Phys. Rev. , 1 (1960).
6] J.A. McLennan, in Advances in Chemical Physics, edited by I. Prigogine (Interscience Pub-lishers, Inc., New York, 1963), Vol. 5.[7] P. R´esibois and M. De Leener,
Classical Kinetic Theory of Fluids (Wiley Interscience, NewYork, 1977).[8] J.W. Dufty, A. Baskaran, and J.J. Brey, submitted to Phys. Rev. E.,cond-mat/0612408[9] A. Baskaran, J.W. Dufty, and J.J. Brey, submitted to Phys. Rev. E.,cond-mat/0612409[10] Hydrodynamics derived from inelastic hard sphere models has found widespread use in multi-phase computational fluid dynamics methods, to describe numerous industrial processes in-volving solid particulates in gas fluidized systems (e.g., software packages such as Fluent 2122and MFIX).[11] V. Garz´o and J.W. Dufty, Phys. Rev. E , 5895 (1999).[12] J.W. Dufty, J.J. Brey, and J. Lutsko, Phys. Rev. E , 051303 (2002); J. Lutsko, J.J. Brey,and J.W. Dufty, Phys. Rev. E , 051304 (2002).[13] V. Garz´o and J.W. Dufty, J. Stat. Phys. , 723-744 (2001).[14] N. N. Bogoliubov, Problems in a Dynamical Theory in Statistical Physics , English trans.by E. Gora in
Studies in Statistical Mechanics I , edited by G.E. Uhlenbeck and J. de Boer(North-Holland, Amsterdam, 1962).[15] E.G.D Cohen, Physica , 1025 (1962); ibid. , 1045 (1962); ibid , 63 (1975); J. R. Dorfman in Perspectives inStatistical Physics , edited by H. Raveche (North-Holland, Amsterdam, 1981).[17] B. Alder and T. Wainwright, Phys. Rev. A , 18 (1970).[18] E. Leutheusser, Phys. Rev. A , 2765 (1984).[19] T. P. C. van Noije and M. H. Ernst, in Granular Gases , editde by T. P¨oschel and S. Luding(Springer, New York, 2001).[20] H. van Beijeren and M.H. Ernst, Physica A , 437 (1973); ibid. , 225 (1973); J. Stat. Phys.
125 (1979).[21] W.E. Alley, B.J. Alder, and S. Yip, Phys. Rev. A , 3174 (1984).[22] I.M. de Schepper, E.G.D. Cohen, C. Bruin, J.C. van Rijs, W. Montfrooij, and L.A. de Graaf,Phys. Rev. A , 271 (1988).[23] J.F. Lutsko, Phys. Rev. E , 061211 (2001).[24] J. W. Dufty and J.J. Brey, Phys. Rev. E , 030302(R) (2003); J.J. Brey and J.W. Dufty, hys. Rev. E , 011303 (2005).[25] J.J. Brey, M.I. Garc´ıa de Soria, P. Maynar, and M.J. Ruiz-Montero, Phys. Rev. E , 051301(2004).[26] J.J. Brey, M.J. Ruiz-Montero, P. Maynar, and M.I. Garc´ıa de Soria, J. Phys.: Condens. Matter , S2489 (2005).[27] J-P. Hansen and I. R. McDonald, Theory of Simple Liquids , (Elsevier Academic Press, London,1986)., (Elsevier Academic Press, London,1986).