Hydrodynamic modes, Green-Kubo relations, and velocity correlations in dilute granular gases

INSTITUTE OF PHYSICS PUBLISHING JOURNAL OF PHYSICS: CONDENSED MATTER J. Phys.: Condens. Matter 17 (2005) S2489–S2502 doi:10.1088/0953-8984/17/24/008 Hydrodynamic modes, Green–Kubo relations, and velocity correlations in dilute granular gases JJavierBrey,MJRuiz-Montero, P Maynar and M I Garc´ ıa de Soria F´ ısica Te´ orica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain E-mail: [email protected] and majos[email protected] Received 27 October 2004, in final form 13 December 2004 Published 3 June 2005 Online at stacks.iop.org/JPhysCM/17/S2489 Abstract It is shown that the hydrodynamic modes of a dilute granular gas of inelastic hard spheres can be identified, and calculated in the long wavelength limit. Assuming they dominate at long times, formal expressions for the Navier– Stokes transport coefficients are derived. They can be expressed in a form that generalizes the Green–Kubo relations for molecular systems, and it is shown that theycan also be evaluatedby means of N-particle simulation methods. The formofthehydrodynamicmodestozerothorderinthegradientsisusedtodetect the presence of inherentvelocity correlationsin the homogeneouscoolingstate, even in the low density limit. They manifest themselves in the fluctuations of the total energy of the system. The theoretical predictions are shown to be in agreement with molecular dynamics simulations. Relevant related questions deserving further attention are pointed out. 1. Introduction Rapid granular flows exhibit a phenomenological behaviour that is similar in many ways to that of ordinary fluids. So it is not surprising that the methods of statistical mechanics and kinetic theory are being successfully extended to describe them. The primary feature of grains that must be captured in any idealized model is the inelasticity of collisions. The prototype is asystemofsmooth inelastic hard spheres (or discs in two dimensions) whose energy loss in collisions is characterized by a coefficient of normal restitution. In real granular media, this coefficient depends on the relative velocity of the colliding particles [1], but for the sake of simplicity, and expecting to retain the main qualitative features, it will be considered here as amaterial constant. At themost fundamental level of particle dynamics, the effect of the inelasticity is easily incorporated by changing the binary collision rule. As in the elastic case, the dynamics of the system consists of a sequence of free streaming between collisions and instantaneous velocity changes on binary collisions. Therefore, a statistical mechanics description can be formulated 0953-8984/05/242489+14$30.00 © 2005 IOP Publishing Ltd Printed in the UK S2489 S2490 JJBreyet al by means of the Liouville equation of the system [2, 3]. From a conceptual point of view, this description has the same degree of justification asthat of inherently non-equilibriummolecular systems. Proceeding to a more mesoscopic description, as providedby kinetictheory, is a more difficult and controversial problem. Formally, the same methods as developed for hard sphere gases with elastic collisions can be applied. In particular, in the low density limit, the inelastic Boltzmann equation is derived.Thus, it seems that the basis for this equation is the same as in the elastic case. A different question is the domain in the parameter space of density and coefficient of restitution in which it actually provides an accurate description of the system [4]. The Boltzmann equation is a kinetic equation describing the time evolution of the oneparticledistribution function, thatis an average quantity. It is also possible to study,again in the low density limit, fluctuations and correlations in the positions and velocities of the particles by means of the associated distribution functions.Itisknown that correlations can become quite large in out of equilibrium elastic systems, on the appropriatelong time and length scales. This shouldnotbe interpretedasa signatureofthefailure ofthe Boltzmannequation. The latter only requiresthatthe precollisionalpart of the correlationsof two particles at contact be small [5, 6]. Closed kinetic equations for the fluctuation and correlation functions, valid in the low density limit, can be derived inside the same approximation schemes that lead to the Boltzmann equation. This includes formal expansions in a small parameter [7] and cluster expansions of the reduced distribution functions [8]. Both methods have been considered for granular gases [9, 10]. Nevertheless, while the inelastic Boltzmann equation and its implications have been extensively studied, little is known about fluctuations and correlations in granular gases, even though it is widely recognized that they are crucial for the understandingof the properties of these systems. As mentioned above, granular gases display fluid dynamical behaviour. Therefore, they are expected to obey some kind of hydrodynamic equations, understood as a set of closed equations for the hydrodynamics fields: density, flow velocity, and (granular) temperature. In fact, more or less phenomenological hydrodynamic equations have been successfully applied to many different situations. Anappropriate context in which the issue of the existence of hydrodynamics for dilute granular gases, the form of the hydrodynamic equations, and their range of validity, can be addressed is provided by the inelastic Boltzmann equation. By using an extension of the Chapmann–Enskog procedure, the Navier–Stokes equations have been derived, and formal expressions for the transport coefficients obtained. They are, in principle, valid for arbitrary physical values of the coefficient of restitution [11]. Recently, these expressions have been expressed in the form of time integrals of the correlations of some velocity functions computed in a given reference state, the so-called homogeneous cooling state [12, 13]. Consistency requires that these integrals exist in the long time limit. The above expressionsare the generalizationof the Green–Kubo relations for elastic gases, although with significative differences that could hardly have been guessed apriori. It must be stressed that the Chapmann–Enskog method and the Green–Kubo relations presume the existence of hydrodynamics on the appropriate length and time scales and, therefore, they do not add any fundamental support to the validity of a hydrodynamic description for granular gases. On the other hand, if Navier–Stokes equations exist, the Chapmann–Enskog results provide their explicit form. A deeper physical understanding of the structure of the Green–Kubo relations and their validity follows from the analysis of the linearized Boltzmann equation, and more specifically of the spectrum of the associated operator [13]. It is possible to identify the hydrodynamic part of the spectrum, that is the one determining the form of the Green–Kuborelations. It is given by d+2points that in the elastic limit correspondto the d+2collisional invariants, where dis the dimension of the system. The keyremaining point for the existence of hydrodynamics is whether the corresponding modes Hydrodynamic modes, Green–Kubo relations, and velocity correlations in dilute granular gases S2491 dominatefor long times and long wavelengths [14, 15]. Quite interestingly, the hydrodynamic modes also determine the form of the fluctuations of the hydrodynamic quantities in granular systems, reflecting a generalization of the fluctuation–dissipation theorems. The aim of this paper is to review some recent advances and present some new results along the above ideas. Particular emphasis will be put on the relevance of the hydrodynamic modes for the study of many fundamental and applied issues in granular fluids, as is the case for molecular, elastic systems. In the paper, and especially in the lastsection, some of the many remaining open questions in the context ofthematterdiscussed here will be indicated. 2. The homogeneous cooling state of a dilute granular gas The system considered is composed of Nsmooth inelastic hard spheres (d=3) or discs (d=2) of mass mand diameter σ.Their dynamics consists of free streaming until two particles iand jare in contact, when their velocities change according to the rule vi→v i≡bσvi=vi−1+α 2ˆ σ·vijˆ σ, vj→vj≡bσvj=vj+1+α 2ˆ σ·vijˆ σ, (1) where vij =vi−vjis the relative velocity, ˆσis the unit vector pointing from the centre of particle jto the centre of particle iat contact, and αis the coefficient of normal restitution. It is defined in the range 0 <α⩽1and measures the inelasticity of collisions. Here it will be taken as a constant. From the above dynamics, the Liouville equation governing the time evolution of the Nparticle distribution function of the system can be derived [2, 3]. An equivalent representation is given bytheBBGKY hierarchy for the s-particle distribution functions (s=1,2,...). In the low density limit, a closed equation for the one-particle distribution function, f(r,v,t), can be derived by using the same procedures as for elastic hard particles [2, 3]. This is the inelastic Boltzmann equation, ∂ ∂t+v·∂ ∂rf(r,v,t)=J[f,f],(2) where J[f|f]isthe inelastic Boltzmann collision operator [2, 3, 16–18]. Foranisolatedgranular fluid, no homogenous steady solution is possible. It is trivially seen that the energy of the system decreases monotonicallydue to the inelasticity of collisions. On the other hand, there is a spatially homogeneous state, the so-called homogeneous cooling state (HCS), for which all the time dependence occurs only through the granular temperature. The latter is defined in the usual way from the average kinetic energy. For the one-particle distribution, this implies a scaling of the form [18] fhcs(v,t)=nHv−d 0[Thcs(t)]χhcs(c), v0[Thcs(t)]=2kBThcs(t) m1/2 .(3) HerenHis the number of particles density,kBis the Boltzmann constant,and Thcs(t)is the timedependent temperature of the HCS. Finally, χhcs(c)is an isotropic function of c=v/v0[T(t)]. The use of equation (3) in equation (2) provides an integral equation for χhcs, ζ0 2 ∂ ∂c ·(cχhcs)=Jc[χhcs,χ hcs],(4) where Jcis the scaled Boltzmann collision operator, Jc[χhcs,χ hcs]=dc1T0(c,c1)χhcs(c)χhcs(c1), (5) S2492 JJBreyet al T0(c,c1)=dˆσ[g·ˆσ]g·ˆσα−2b−1 σ(c,c1)−1,(6) with g≡c−c1,the Heaviside step function, and b−1 σan operator replacing the velocities cand c1to its right by the precollisional values obtained by inverting the collision rules given in equation (1). Moreover, the scaling property(3)implies that the time dependence of the temperature is given by the Haff law: Thcs(t)=Thcs(0)1+1 2ζhcs(0)t−2.(7) The cooling rateoftheHCS,ζhcs(t),isgivenby ζhcs(t)=(1−α2)π d−1 2σd−1nHv0(t) 2d+3 2ddcdc1g3χhcs(c)χhcs(c1), (8) and ζ0=ζhcs(t) v0[Thcs(t)],(9) where =1/nHσd−1is proportionalto the mean free path of the particles. Since ζhcs ∝T1/2 hcs , equation (7) becomes independent of theinitial condition in the long time limit, Thcs ∼4(¯ ζt)−2,¯ ζ=ζhcs(t) T1/2 hcs (t) ,(10) indicating that all the homogeneous cooling states of a given system tend to converge in the long time limit [30]. This seems to suggest the existence of some kind of Htheorem for the homogeneous Boltzmann equation for inelastic hard particles, although no proof of such a theorem has been given to date. 3. Linearized Boltzmann equation and hydrodynamic modes Now, the initial value problem for small perturbations around the HCS of an isolated granular gas will be considered. In this case, the Boltzmann equation can be linearized around the HCS, and it will be shown that, in appropriate dimensionless units, its general solution can be related to a linear eigenvalue problem. By considering times over which the deviation of the distribution function from fhcs remains small, we can write f(r,v,t)=fhcs(v,t)+δf(r,v,t), |δf(r,v,t)|fhcs(v,t). (11) In order to eliminate the time dependence of the HCS, it is convenient to introduce the scaled velocity cdefined below equation (3) and the dimensionless quantities l=r ,s=t 0 dtv0(t) ,δχ( l,c,s)=n−1 Hvd 0(t)δ f(r,v,t). (12) The time scale sis proportional to the average of the accumulated number of collisions per particle in the interval (0,t).Then, substitution of equation (11) into the Boltzmann equation, retaining only terms up to first order in the deviations, gives the linear Boltzmann equation ∂ ∂sδχ +c·∂ ∂lδχ =δχ(l,c,s), (13) where the linear operator is given by δχ =Jc[χhcs,δχ]+Jc[δχ,χhcs]−ζ0 2 ∂ ∂c ·(cδχ).(14) Hydrodynamic modes, Green–Kubo relations, and velocity correlations in dilute granular gases S2493 In this representation, the effect of the coolingappears explicitly through the last term on the right-hand side of the above expression. Solutions to equation (13) will be sought in a Hilbert space defined by the scalar product g|h=dcχ−1 hcs(c)g∗(c)h(c), (15) where g∗denotesthecomplexconjugateofg.Itmustbenotedthatthisisastrongerrequirement than the necessary normalization of the distribution function. Since the collision operator  does not affect the space variable, it is enough to consider a single Fourier mode, i.e., δχ(l,c,s)=eik·lδ˜χ(k,c,s), (16) so that the linearized Boltzmann equation becomes ∂ ∂sδ˜χ=(−ik·c)δ˜χ. (17) All linear excitations are determined from the spectrum of the operator −ik·c.The hydrodynamic part of this spectrum is defined as follows. For k=0itisgivenby those eigenvaluesthatcoincidewiththeeigenvaluesoftheexactbalanceequationsforthemomentum density, flow velocity, and temperature obtained from the linearized homogeneousBoltzmann equation, in the same limit k=0, i.e., for homogeneousstates. Assuming that the cooling rate for homogeneous systems can be expressed as a function of the density and the temperature, these eigenvalues are found to be [13] λ1=0,λ 2=ζ0/2λ3=−ζ0/2,(18) the eigenvalue λ2being d-fold degenerated. For finite k,the hydrodynamic modes are defined as the solutions of the equation (−ik·c)ξi(k,c)=λiξi(k,c)(19) that are continuously connected as functions of kto the above special solutions for k=0. Of course, consistency requires that includes the spectrum of the balance equations for k=0. Straightforward calculations show that [13] ξi(c)=λiξi(c), i=1,2,3(20) with the λigiven by equation (18) and ξ1(c)=χhcs(c)+∂ ∂c ·[cχhcs(c)], ξ2(c)=−∂χhcs(c) ∂c,ξ 3(c)=−∂ ∂c ·[cχhcs(c)]. (21) This confirms the existence of the hydrodynamic modes. Their evaluation for finite kis now atechnical problem. To Navier–Stokes order, they have been computed by using perturbation expansion methods [14], and the results agree with those obtained from the Navier–Stokes equations derived by means of the Chapmann–Enskog procedure [11]. Here an alternative approach, designed to derive in a direct way expressions for the Navier–Stokes transport coefficients, will be outlined. As a consequence of being not Hermitian, the eigenfunctions ξiare not orthogonal. It is then useful to introduce a set of functions {¯ ξi}that are biorthogonal to the set {ξi,i=1,2,3}. Aconvenient choice is {ξi}=χhcs(c), cχhcs(c), c2 d+1 2χhcs(c),ξi|ξj=δi,j.(22) S2494 JJBreyet al The general solution to the linearized Boltzmann equation in the Hilbert space can be formally written as δ˜χ(k,c,s)= i˜ai(k,s)χi(k,c)+δ˜χm(k,c,s), (23) where the sum extends over the d+2hydrodynamic modes and the term δ˜χmcontains all the other ‘microscopic’ excitations. Assuming that the hydrodynamic spectrum dominates for long times and small gradients, δ˜χmcan be neglected, and the coefficients ˜aiare given by ˜ai(k,s)=¯ ξi|δ˜χ(k,s),(24) since δ˜χnow lies in the subspace spanned by the eigenfunctions {ξi}.Thesecoefficients are combinationsof the hydrodynamicfields, i.e., they can be expressed in terms of the deviations of the macroscopic density n,temperature Tand fluid velocity u.Itisnowpossible to derive aformalexpression for δ˜χ(k,c,s)to first order in k,andfrom it the Navier–Stokes order for thepressure tensor andthe heat conductivity[13]. Theresults will be given in the next section. Theabove discussionreliesontheassumptionthatthehydrodynamic point spectrumfound here is isolated and smaller than all the other parts of the spectrum, so that it dominatesfor long times and long wavelengths. This is a necessary condition for the existence of hydrodynamics and has been proven in the elastic case. The proof does not extend directly to the inelastic case, although there are some reasons to expect that this is the case even for rather strong inelasticity. A detailed discussion of the basis for a hydrodynamic description of a granular gas in this context is given in [15]. 4. Green–Kubo representation of the transport coefficients The procedure described in the previous section leads to the following expressions for the dissipative part of the pressure tensor ij(r,t)and the heat flux q(r,t)to Navier–Stokes order [13]: ij =−η∂ui ∂rj +∂uj ∂ri−2 dδij∇·u,q=−κ∇T−µ∇n,(25) where ηis the shear viscosity, κthe(thermal) heat conductivity, and µanewtransport coefficient that vanishes in the elastic limit and that will be referred to as the diffusive heat conductivity. They are given by η=nHmv0(T)∞ 0 dse−sζ0/2xy(s)2,xy,(26) κ=nHv0(T)∞ 0 dsesζ0/2x(s)3,x,(27) µ=mv3 0(T)∞ 0 dsesζ0/2−1x(s)3,x+1 2∞ 0 dsx(s)cx.(28) In the above expressions, xy and xare the transversal momentum flux and the heat flux, respectively, xy(c)=cxcy, x(c)=c2−d+2 2cx,(29) while 2,ij and 3,xare ‘modified fluxes’, 2,xy(c)=−cx ∂ln χhcs(c) ∂cy , 3,x(c)=−cx 2d+c·∂ln χhcs(c) ∂c.(30) Hydrodynamic modes, Green–Kubo relations, and velocity correlations in dilute granular gases S2495 The angular brackets denote ‘time-correlation functions’ in the HCS, g(s)h≡dcχhcs(c)g(c,s)h(c), (31) where the time dependence is given by g(c,s)=es˜ (c)g(c), (32) with ˜ (c)=dc1χhcs(c1)dˆσ(g·ˆσ)g·ˆσ[bσ(c,c1)−1](1+P)+ζ0 2c·∂ ∂c.(33) Here, Pis an operator that interchanges the velocities cand c1to its right.Equations (26)– (28) are the generalization for inelastic hard collisions of the dilute limit of the well-known Green–Kubo relations for molecular fluids, reducing to them in the elastic limit [12]. They can be understood as generalized fluctuation–dissipationrelations, in the sense that they relate thetransport coefficients characterizing dissipation in inhomogeneous systems with the time decay of fluctuations in the homogeneous reference state (the HCS). Equivalent expressions for the transport coefficients are obtained by solving the nonlinear Boltzmann equation by theChapman–Enskog method [11, 12] and also by perturbation analysis of the hydrodynamic spectrum of the linear Boltzmann equation [14]. The effect of the inelasticity of collisions manifests itself in several ways. Perhaps the most unexpected among them is that the initial dissipative fluxes xy and xare replaced by the modified fluxes 2,xy and 3,xin the time correlation functions appearing in the expressions of ηand κ.Thisisadirect consequence of the form of the hydrodynamic modes. While in a molecular system they are given by linear combinations of the Maxwellian times 1, c,andc2,ingranular gases they are given by equation (21), and it should be realized that χhcs(c)is not a Maxwellian outside the elastic limit α=1. Since the diffusive heat conductivity µis a peculiarity of granular systems, it seems appropriateto comment aboutits relevance. In general, its contributionto the heat flux is much smaller than the one coming from the usual heat conductivity associated with the temperature gradient, because κµfor all α.Nevertheless, there are some cases where the existence of this contribution manifests itself in a clear way. Consider, for instance, an open granular system in the presence of gravity. The system is fluidized by means of a vibrating bottom plate. In the steady state, the granular temperature profile presents a minimum at a certain height, increasing from there on. This temperature inversion has been observed in computer simulations [19–21] and also in experiments [22–25]. Its theoretical explanation is tied to the existence of the diffusive contribution to the heat flux [19], and the value of the transport coefficient µcan be directly obtainedfrom the measurementof the heat flux at the temperature minimum [26]. For a low density gas, the results are in good agreement with the theoretical predictions discussed in the next section. 5. Simulation study of the Green–Kubo relations The formal expressions of the transport coefficients given by equations (26)–(28) are rather involved and have been evaluated analytically only in some approximations, namely by introducing an expansion in orthogonal functions, keeping only rather arbitrarily the lowest orders contributions [11, 27, 28]. On the other hand, the above Green–Kubo relations can be transformed into a form that is suitable for evaluation by means of N-particle simulation techniques. Two main points must be addressed for that. First, there is the technical difficulty S2496 JJBreyet al 0.2 0.4 0.6 0.8 1.0 1.0 1.2 1.4 1.6 η∗ 0.2 0.4 0.6 0.8 1.0 1.0 1.5 2.0 2.5 κ∗ αα Figure 1. Left: the dimensionless reduced coefficient of shear viscosity η∗as a function of α. Thesymbols are from the DSMC method, while the solid curve is the analytical result obtained in the first Sonine approximation [11]. Right: the same for the reduced coefficient of thermal heat conductivity κ∗. that the operator ˜ ,defining the time dependence of the fluxes, involves the cooling rate ζ0 that, therefore, should be known apriori,anditis not. This can be overcome by making a change in the time scale and realizing the existence of an exact mapping of the HCS onto a steady state [29, 30]. The second issue is to transform equations (26)–(28) into expressions involving a well-defined N-particle dynamics, instead of the effective Boltzmann dynamics defined by ˜ .This can be done by using the same kind of assumptions as are needed to derive the Boltzmann equation and, in particular, assuming that velocity correlations in the HCS are negligible in the low density limit [30]. We have used the direct simulation Monte Carlo (DSMC) method [31] to simulate the N-particle dynamics of a dilute granular gas composed of inelastic hard spheres (d=3). The details of the practical implementation of the above ideas and the simulation technique are discussed in [30] and [32]. Here it will only be stressed that the DSMC method allows us to simulate a homogeneousstate in such a way that it stays homogeneous‘bydefinition’,avoiding thespontaneousdevelopingofspatialinhomogeneities. Thisisimportantinthepresentcontext, since the HCS is known to be unstable against spatial long wave perturbations, even in the low density limit. This is the so-called clustering instability [33, 34]. In the simulations, it is observed that the HCS, or more precisely its equivalent steady representation, is reached after a few collisions per particle. This supports that the hydrodynamic point spectrum identified in section 4 is isolated and smaller than all the other parts ofthe spectrum. Moreover,all the timecorrelationfunctionsappearingin equations(26)– (28) are found to decay exponentially, at least until they decay by two orders of magnitude fromtheir initial values. This has been exploited in order to compute the time integrals. The results obtained for the shear viscosity are shown in figure 1 (left). What is actually plotted is the shear viscosity scaled with its elastic value in the first Sonine approximation, i.e., η∗≡η(T) η0(T),η 0(T)=5(mkBT)1/2 16π1/2σ2.(34) Also plotted in the same figure is the approximate theoretical prediction obtained when computingequation(26) in the first Sonine approximation[11]. It is observed that theviscosity monotonicallyincreaseswiththeinelasticityofcollisions(decreasingα)andthatthetheoretical prediction is in good agreement with the simulation results over all the range of values of α. Hydrodynamic modes, Green–Kubo relations, and velocity correlations in dilute granular gases S2497 0.2 0.4 0.6 0.8 1.0 α 0.0 0.5 1.0 1.5 2.0 µ∗ Figure 2. Thesame as in figure 1, but for the reduced coefficient of diffusive heat conductivity µ∗. The simulation results for the (thermal) heat conductivity are given in figure 1 (right). Similarly to above, the dimensionless plotted quantity is κ∗≡κ(T) κ0(T),κ 0(T)=75kB(kBT)1/2 64(πm)1/2σ2.(35) Again, it is found that the transport coefficient is a monotonic decreasing function of α.However, in this case the first Sonine approximation leads to significant quantitative discrepancies for small values of α.Thismaybe, at least partially, due to the fact that the fluxesappearing in the expression of the heat conductivity involve higher velocity moments than those present in the expression of the shear viscosity. The dimensionless reduced diffusive heat conductivity µ∗has been defined as µ∗≡nµ(T) Tκ0(T),(36) and the results obtained for it are presented in figure 2. As is the case for the thermal heat conductivity, the first Sonine approximation significantly overestimates this transport coefficient for strong dissipation. The analysis of the transport coefficients for a dilute gas of inelastic hard discs (d=2) has been reported in [32]. The results are quite similar, the main difference being that the discrepancies of the first Sonine approximation predictions are stronger for d=2thanfor d=3, as is usually the case. As already pointed out, in order to transform the effective one-particle dynamics, implicit in the linearized evolution operator ˜ ,intoanN-particle dynamics, it has been assumed that velocity correlations are negligible in the HCS of a low density granular gas. Nevertheless, when the DSMC data are carefully analysed,itisseen that for strong inelasticity (α0.6), velocity correlation effects are clearly identified [32]. For instance, their contribution to the initial value of the time-correlation functions appearing in the N-particle representation of the Green–Kubo expressions can be up to 10% of the total value. This raises some fundamental issues. Of course, a first possibility is that the velocity correlations are just an artefact introduced by the simulation method itself, i.e., an inaccuracy of the DSMC method. On the other hand, if they are not, how can these correlations be incorporated in a theoretical description of dilute granular gases? Is it consistent with the description provided by the inelastic Boltzmann equation? A partial answertothesequestions will be given in the next section.