Dynamics of Annihilation I : Linearized Boltzmann Equation and Hydrodynamics Mar´ıa Isabel Garc´ıa de Soria,1 Pablo Maynar,2,3 Gr´egory Schehr,2 Alain Barrat,2 and Emmanuel Trizac1 1Universit´e Paris-Sud, LPTMS, UMR 8626, Orsay Cedex, F-91405 and CNRS, Orsay, F-91405 2Laboratoire de Physique Th´eorique (CNRS UMR 8627), Bˆatiment 210, Universit´e Paris-Sud, 91405 Orsay cedex, France 3F´ısica Te´orica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain (Dated: February 2, 2008) 8 We study the non-equilibrium statistical mechanics of a system of freely moving particles, in 0 whichbinaryencounterslead eithertoanelastic collision ortothedisappearance ofthepair. Such 0 asystemofballisticannihilationthereforeconstantlyloosesparticles. Thedynamicsofperturbations 2 aroundthefreedecayregimeisinvestigatedfromthespectralpropertiesofthelinearizedBoltzmann n operator,thatcharacterize linearexcitationsonalltimescales. Thelinearized Boltzmannequation a issolvedinthehydrodynamiclimit byaprojection technique,whichyieldstheevolutionequations J fortherelevantcoarse-grainedfieldsandexpressionsforthetransportcoefficients. Wefinallypresent 5 theresults of Molecular Dynamics simulations that validate thetheoretical predictions. 1 PACSnumbers: 51.10.+y,05.20.Dd,82.20.Nk ] h c I. INTRODUCTION e m Understanding the differences and similarities between a flow of macroscopic grains and that of an ordinary liquid - t is anactive field ofresearch[1, 2]. From a fundamental perspective,it is tempting to drawa correspondencebetween a t the grainsofthe formerandthe atomsofthe latterinordertomakeuse ofthe powerfultoolsofstatisticalmechanics s toderivealargescaledescriptionforthe variousfields ofinterest,suchasthe localdensity ofgrains. Akeydifference . t a betweenagranularsystemandanordinaryliquidisthatcollisionsbetweenmacroscopicgrainsdissipateenergy,dueto m the redistribution of translationalkinetic energy into internal modes. This simple fact has far reaching consequences [2,3],butalsoposesana prioriseriousproblemconcerningthevalidityoftheprocedureleadingtothehydrodynamic - d description. Indeed, the standardapproachretains in the coarse-graineddescription only those fields associatedwith n quantities that are conserved in collisions (such as density and momentum). There is however good evidence –both o numerical and theoretical– that in the granular case, a relevant description should include the kinetic temperature c field, defined as the kinetic energy density ([2, 4] and references therein), which is therefore not associated to a [ conserved quantity. 1 Ourpurposehereistotestahydrodynamicdescriptionwithsuitablecoarse-grainedfields,foramodelsystemwhere v notonly the kinetic energyis notconservedduringbinary encounters,but alsothe number of particlesandthe linear 9 momentum. The ballistic annihilationmodel [5–10]providesa valuable candidate: inthis model, eachparticle moves 9 2 freely (ballistically) until it meets another particle; such binary encounters lead to the annihilation of the colliding 2 pairofparticles. Inaddition,weintroduceaparameter0 p 1thatmaybethoughtofasameasureofthedistance ≤ ≤ . to equilibrium, so that anensemble of sphericalparticlesin dimension d undergoingballistic motion either annihilate 1 upon contact (with probability p) or scatter elastically (with probability 1 p). For the corresponding probabilistic 0 − 8 ballisticannihilationmodel,theChapman-Enskog[11]schemewasappliedrecently[12]. Thehydrodynamicequations 0 were derived and explicit formulas for the transport coefficients obtained. Our goal here is two-fold. First, we would : liketoshedlightonthe contextandlimitations ofthe derivation,by obtainingthe hydrodynamicdescriptiondirectly v from the linearized Boltzmann equation. Second, we aim at putting to the test the theoretical framework thereby i X obtained by careful comparison with numerical simulations of the annihilation process. For granular gas dynamics, r the same program is quite complete, although challenges remain [1, 2]. The objective here is to initiate a similar a formulation for the ballistic annihilation model in view of a more stringent test of the hydrodynamic machinery. The paper is organized as follows. We start in section II with a reminder of results derived in Refs [8, 9]. The kinetic description adopted is that of the Boltzmann equation, since it has been shown that for p = 1 (all collision events leading to annihilation), the underlying molecular chaos closure provides an exact description at long times, provided space dimension d is strictly larger than 1 [9]. We may assume that the same holds for an arbitrary but non vanishing value of p, since the density is then still a decreasing function of time. The focus is here on an unforced system, which is characterized by an algebraic decay with time of the total density and kinetic energy density (homogeneous decay state) [8, 9]. More precisely, we are interested in small perturbations around this state, so that the Boltzmann equation will be subsequently linearized. After having identified the operator that generates the dynamics of fluctuations, attention will be paid in section III to its spectral properties. This will provide the basis for finding in section IV the evolution equations for the hydrodynamic fields (i.e. those chosen for the coarse- 2 grained description) and for obtaining explicit formulas for the transport coefficients. Finally, our predictions will be confronted in section V against extensive Molecular Dynamics simulations. Such a comparison is an essential step in testing the foundations of the hydrodynamic treatment. II. THE BOLTZMANN EQUATION APPROACH TO THE HOMOGENEOUS DECAY STATE A. Non-linear description The Boltzmann equation describes the time evolution of the one particle distribution function f(r,v ,t). For a 1 system of smooth hard disks or spheres of mass m and diameter σ, which annihilate with probability p or collide elastically with probability 1 p, it has the form − ∂ +v f(r,v ,t)=pJ [f f]+(1 p)J [f f], (1) 1 1 a c ∂t ·∇ | − | (cid:18) (cid:19) where the annihilation operator J is defined by [9] a J [f g]= σd−1 dv dσˆΘ(v σˆ)(v σˆ)f(r,v ,t)g(r,v ,t). (2) a 2 12 12 1 2 | − · · Z Z The elastic collision operator J reads [13, 14] c J [f g]=σd−1 dv dσˆΘ(v σˆ)(v σˆ)(b−1 1)f(r,v ,t)g(r,v ,t), (3) c | 2 12· 12· σ − 1 2 Z Z with v =v v , Θ the Heaviside step function, σˆ a unit vector joining the centers of the two particles at contact 12 1 2 and b−1 an ope−rator replacing all the velocities v and v appearing in its argument by their precollisionalvalues v∗ σ 1 2 1 and v∗, given by 2 b−1v =v∗ = v (v σˆ)σˆ, (4) σ 1 1 1− 12· b−1v =v∗ = v +(v σˆ)σˆ. (5) σ 2 2 2 12· We assume that the system can be characterized macroscopically by coarse grained (hydrodynamic-like) fields, that we define as in standard Kinetic Theory in terms of the local velocity distribution function f(r,v,t) n(r,t) = dvf(r,v,t), (6) Z n(r,t)u(r,t) = dvvf(r,v,t), (7) Z d m n(r,t)T(r,t) = dv V2f(r,v,t), (8) 2 2 Z where n(r,t), u(r,t), and T(r,t) are the local number density, velocity, and temperature, respectively. We have introduced here V = v u, the velocity of the particle relative to the local velocity flow. We stress that the − temperature defined has a kinetic meaning only, but lacks a thermodynamic interpretation. It seems natural to considerthesefields,astheyaretheusualhydrodynamicalfieldsoftheequilibriumsystem(withp=0). Itishowever not obvious at this point that restricting our coarse-graineddescription to the above three fields provides a relevant and consistent framework. A major goal of this paper is to provide strong hints that this is indeed the case. We will show in particular that closed equations can be obtained for these fields in the appropriate time and length scales, under reasonable assumptions. The Boltzmann equation (1) admits a homogeneous scaling solution f in which all the time dependence is em- H bedded in the hydrodynamic fields, with the further simplification that those fields are position independent. The existence of this regime could notbe shownrigorously,but, numerically,such a scaling solutionquickly emerges from an arbitrary initial condition [8, 9]. It has the form [9] n (t) v f (v,t)= H χ (c), with c= , (9) H v (t)d H v (t) H H 3 where 1/2 2T (t) H v (t)= (10) H m (cid:20) (cid:21) isthe“thermal”(root-mean-square)velocityandχ (c)isanisotropicfunctiondependingonlyonthemodulusc= c H | | of the rescaledvelocity. By taking moments in the Boltzmann equation and using the scaling (9), it can be seen that the homogeneous density and temperature obey the equations [12] ∂n (t) H = pν (t)ζ n (t), (11) H n H ∂t − ∂T (t) H = pν (t)ζ T (t), (12) H T H ∂t − where we have introduced the collision frequency of the corresponding hard sphere fluid in equilibrium (with same temperature and density) ν (t)= nH(t)TH1/2(t)σd−1 8πd−21 (13) H m1/2 (d+2)Γ(d/2) and the dimensionless decay rates ζ and ζ , that are functionals of the distribution function n T γ pζ = dc dc T(c ,c )χ (c )χ (c ), (14) n 1 2 1 2 H 1 H 2 −2 Z Z γ 2c2 pζ = dc dc 1 1 T(c ,c )χ (c )χ (c ). (15) T 1 2 1 2 H 1 H 2 −2 d − Z Z (cid:18) (cid:19) In these expressions, γ is a quantity that does not depend on time, which reads 2n (t)v (t)σd−1 (d+2)√2Γ(d/2) H H γ = = , (16) ν (t) 4π(d−1)/2 H and the binary collision operator T(c ,c ), that should not be confused with the temperature, takes the form 1 2 T(c ,c )= dσˆΘ(c σˆ)(c σˆ)[(1 p)b−1 1]. (17) 1 2 12· 12· − σ − Z Finally, we can write anequationfor the scaleddistribution function χ (c) in terms of the coefficients andoperators H defined above ∂ p (dζ 2ζ )+ζ c χ (c )=γ dc T(c ,c )χ (c )χ (c ). (18) T − n T 1· ∂c H 1 2 1 2 H 1 H 2 (cid:20) 1(cid:21) Z The operator b−1 in the last equation is defined again by equation (4), but substituting (v ,v ) by (c ,c ). σ 1 2 1 2 Althoughanexactandexplicitsolutionofequation(18)isnotknown,its behavioratlargeandsmallvelocitieshas been determined [8, 9]. In this work we will use the approximate form of the distribution function in the so-called first Sonine approximation (an expansion around a Gaussian functional form, see Appendix A), which is valid for velocitiesinthethermalregion,andallthefunctionalsofχ (c), thatisthedecayratesandthetransportcoefficients, H will be evaluated in this approximation [8, 15]. B. Linearized Boltzmann Equation In the remainder, we consider a situation where the system is very close to the homogeneous decay state, so that we can write f(r,v ,t)=f (v ,t)+δf(r,v ,t), δf(r,v ,t) f (v ,t). (19) 1 H 1 1 1 H 1 | |≪ Substitution of equation (19) into the Boltzmann equation (1), keeping only linear terms in δf, yields ∂ +v δf(r,v ,t) 1 1 ∂t ·∇ (cid:18) (cid:19) =p J [δf f ]+J [f δf] +(1 p) J [δf f ]+J [f δf] , (20) a H a H c H c H { | | } − { | | } 4 Given the scaling form of f (Eq. (9)), it is convenient to introduce as well the scaled deviation of the distribution function, δχ, as follows n (t) δf(r,v ,t)= H δχ(r,c ,τ) . (21) 1 v (t)d 1 H Moreover,Eqs. (11) and (12) suggest to use the dimensionless time scale τ defined by 1 t τ = dt′ν (t′), (22) H 2 Z0 which counts the number of collision per particle in the time interval [0,t]. Combining Eq. (13) together with Eq. (11, 12) yields immediately ν (t)=(1/ν (0)+p(ζ +ζ /2)t)−1 and thus H H n T 1 τ = log[1+ν (0)p(ζ +ζ /2)t]. (23) H n T p(2ζ +ζ ) n T In this time scale τ (22), these equations (11, 12) are easily integrated, yielding n (τ)=n (0)exp( 2pζ τ), T (τ)=T (0)exp( 2pζ τ), (24) H H n H H T − − andpower law behaviorsin time t, nH(t) t−2ζn/(2ζn+ζT) andTH(t) t−2ζT/(2ζn+ζT) atlargetime t 1 . It proves also convenient to introduce Fourier com∝ponents [with the notation∝hk = drexp−ik·rh(r)] so that≫the evolution equation for a general k component of δχ is, in the τ timescale, R ∂ δχk(c1,τ)=[Λ(c1) ilH(τ)k c1]δχk(c1,τ). (25) ∂τ − · In this equation, the time dependent length scale l =2v (τ)/ν (τ) is proportional to the instantaneous mean free H H H path (l (τ) n−1(τ), see Eq. (13)) and the homogeneous scaled Boltzmann linear operator reads H ∝ H Λ(c )h(c )=γ dc T(c ,c )(1+ )χ (c )h(c ) 1 1 2 1 2 12 H 1 2 P Z ∂ +p(2ζ dζ )h(c ) pζ c h(c ). (26) n− T 1 − T 1· ∂c 1 1 In this expression,the permutation operator interchanges the labels of particles 1 and 2 and subsequently allows 12 P for more compact notations. In the present representation, all the time dependence due to the reference state is absorbed in the mean free path, obtained from l (τ) n−1(τ) as H ∝ H l (τ)=l (0)exp(2pζ τ), (27) H H n which, as expected, is an increasing function of time. C. Linearized Hydrodynamic Equations around the homogeneous decay state Let us define the relative deviations of the hydrodynamic fields from their homogeneous values by δn(r,τ) ρ(r,τ) = dcδχ(r,c,τ), (28) ≡ n (τ) H Z δu(r,τ) w(r,τ) = dccδχ(r,c,τ), (29) ≡ v (τ) H Z δT(r,τ) 2c2 θ(r,τ) = dc 1 δχ(r,c,τ), (30) ≡ T (τ) d − H Z (cid:18) (cid:19) where δy(r,τ) y(r,τ) y (t) denotes the deviation of a local macroscopic variable, y(r,τ), from its homogeneous H ≡ − decay state value, y (t). Taking velocity moments in the Boltzmann equation (25), we obtain the linearized balance H 5 equation for the k components of the hydrodynamic fields ∂ 2pζn ρk+ilH(τ)k wk p δζn[δχk]=0, (31) ∂τ − · − (cid:18) (cid:19) ∂ p(2ζn+ζT) wk ∂τ − (cid:20) (cid:21) i + lH(τ)k(ρk+θk)+ilH(τ)k Π[δχk] p δζu[δχk]=0, (32) 2 · − ∂ 2p(ζn+ζT) θk ∂τ − (cid:20) (cid:21) 2 2pζTρk+i lH(τ)k (wk+φ[δχk]) p δζT[δχk]=0. (33) − d · − Here, we have introduced the traceless pressure tensor and the heat flux as Π[δχk] = dc∆(c)δχk(c,τ), (34) Z φ[δχk] = dcΣ(c)δχk(c,τ), (35) Z where ∆ and Σ are defined as c2 ∆ (c) = c c δ , (36) ij i j ij − d d+2 Σ(c) = c2 c, (37) − 2 (cid:18) (cid:19) and the functionals p δζn[δχ]=γ dc1 dc2T(c1,c2)(1+ 12)χH(c1)δχk(c2,τ), (38) P Z Z p δζu[δχ]=γ dc1 dc2c1T(c1,c2)(1+ 12)χH(c1)δχk(c2,τ), (39) P Z Z 2c2 p δζT[δχ]=γ dc1 dc2 1 1 T(c1,c2)(1+ 12)χH(c1)δχk(c2,τ). (40) d − P Z Z (cid:18) (cid:19) The previous analysis therefore amounts to obtaining a set of complicated equations expressing the evolution of the hydrodynamic fields as a function of the rescaledhomogeneous distribution function χ and the perturbationδχ. In H orderto obtaina closedset ofequationsfor the hydrodynamic fields (31),(32), (33),we needtherefore to expressthe functionals Π, φ, δζ , δζ and δζ , in terms of the hydrodynamic fields themselves. We will see in the next section n u T that, as long as we can treat l (τ)k as a small parameter and if the linear Boltzmann operator has some specific H properties, it is possible to carry out this program and to close the linear hydrodynamic equations. However, since the meanfree path l (τ) increaseswith time (Eq. (27)), the requirementof a smalll (τ)k is necessarilylimited to a H H time window depending onbothkandthe probabilityofannihilationp. Anupper bound forthis window isprovided by the time when the mean free path becomes of the order of the system size. III. SOLUTION OF THE LINEARIZED BOLTZMANN EQUATION In this Section we explore the solutions to the linearized Boltzmann equation (25) and establish some properties of the homogeneous linear Boltzmann operator that will be essential for the coarse-grained description. From the expressionof the linearizedBoltzmannequation, we canidentify the operatorΛ ik cl (τ) asthe “generatorof the H − · dynamic” of δχk. As we are interested in the solutions of this equation in the hydrodynamic regime (large enough scales), it is convenient to study first the eigenvalue problem of the homogeneous linear Boltzmann operator. The inhomogeneous term will be treated perturbatively later on. 6 A. Hydrodynamic Eigenfunctions of Λ Let us consider the eigenvalue problem of the homogeneous linear Boltzmann operator Λ Λ(c)ξ (c)=λ ξ (c). (41) β β β Finding all the solutions of this equation is an insurmountable task. Nevertheless, it is possible to obtain some particular solutions, which will turn out to be the relevant ones in the hydrodynamic regime. The problem will be posed in a Hilbert space of functions of c with scalar product given by g h = dcχ−1(c)g∗(c)h(c), (42) h | i H Z where g∗ denotes the complex conjugate of g. Of particular interest here are the eigenfunctions and eigenvalues associatedwith linear hydrodynamics. Following [16, 17], we use the fact that the homogeneous decay state is parameterized by the hydrodynamic fields n , T and H H u . Writing the Boltzmann equation satisfied by χ and differentiating it with respect to these fields allows then to H H obtainthreeexactrelationsfromwhichonecanextracteigenfunctionsofthe linearizedBoltzmanncollisionoperator. In Appendix B, we show that the functions ∂ ξ (c) = χ (c)+ [cχ (c)], (43) 1 H ∂c · H ∂ ξ (c) = zχ (c) [cχ (c)], (44) 2 H − ∂c · H ∂ ξ (c) = χ (c), (45) 3 −∂c H with z =2ζ /ζ , are solutions of Eq. (41), with eigenvalues n T λ =0, λ = p(ζ +2ζ ), λ =pζ , (46) 1 2 T n 3 T − respectively, λ being d-fold degenerate. Although we cannot prove in general that these eigenvalues are indeed the 3 hydrodynamic ones (i.e the upper part of the spectrum), we will assume that this is the case ; the self-consistency of theapproachandcomparisonwithnumericalsimulationswillvalidatethisassumption. Interestingly,intheparticular case of Maxwell molecules where the full spectrum of Λ may be computed exactly (see Appendix C), it appears that theabove“hydrodynamic”modesdominateatlongtimes,providedthatp<1/4. Forlargervaluesofp,the“kinetic” mode with largest eigenvalue decays slower than one of the three “hydrodynamic” modes. As a consequence of the non-hermitian character of the operator Λ, the functions ξ are not orthogonal β β=1,...,3 { } with respect to the scalar product defined in (42). They are nevertheless independent and, in order to define the projection onto the subspace spanned by these functions, it is necessary to introduce a set of functions ξ¯ β β=1,...,3 { } verifying the biorthonormality condition ξ¯β ξβ′ =δβ,β′. (47) h | i Although the set ξ¯ is not unique, a convenient choice is given by β β=1,...,3 { } 2+z z c2 ξ¯(c) = χ (c), (48) 1 H 2(1+z) − 1+z d (cid:20) (cid:21) 1 1 c2 ξ¯(c) = + χ (c), (49) 2 H 2(1+z) 1+z d (cid:20) (cid:21) ξ¯ (c) = cχ (c). (50) 3 H Indeed, the functions ξ¯ have to be linear combinations of χ (c), cχ (c) and c2χ (c) to ensure that β β=1,...,3 H H H projectionofδχk ontot{he}ξ¯β yieldsthecoarse-grainedfieldsρk,θk andwk,orcombinationsthereof. Thefunctions ξ¯ span a dual s{ubs}pace of that spanned by the eigenfunctions and for any linear combination of the β β=1,...,3 { } hydrodynamic modes 3 g(c)= a ξ (c), (51) β β β=1 X 7 the coefficients a are given by β a = ξ¯ g = dcχ−1(c)ξ¯ (c)g(c). (52) β h β| i H β Z In particular,the projectionof the distribution function δχk onthe subspace spannedby the functions ξβ is givenby the coefficients 1 z 1 1 ξ¯β δχk = ρk θk, ρk+ θk,wk . (53) {h | i} 1+z − 2(1+z) 1+z 2(1+z) (cid:26) (cid:27) Notably, these coefficients are simply linear combinations of the hydrodynamic fields linearized around the homoge- neous decay state. B. Projection of the Linearized Boltzmann Equation on the Hydrodynamic Subspace In this section, we study the linearized Boltzmann equation on the hydrodynamic subspace. Let us define the projectors 3 Ph(c)= ξ¯ h ξ (c), (54) β β h | i β=1 X and P =1 P. (55) ⊥ − so that any function can be decomposed as h(c)=Ph(c)+P h(c). (56) ⊥ In the definition (54) we are considering the functions (43)-(45) and (48)-(50) defined above. Letusnowconsiderthefunctionδχk. IfweapplytheprojectorsP andP⊥ toequation(25),weobtainthefollowing relations ∂ P(Λ ilHk c)P Pδχk = PilHk cP⊥δχk+PΛP⊥δχk, (57) ∂τ − − · − · (cid:20) (cid:21) ∂ P⊥(Λ ilHk c)P⊥ P⊥δχk = P⊥ilHk cPδχk, (58) ∂τ − − · − · (cid:20) (cid:21) where we have used that P ΛP =0 , (59) ⊥ which is obtained straightforwardlysince ξ are right-eigenfunctions of Λ. We note however that the ξ¯ are β β β=1,...,3 { } not left-eigenfunctions of Λ, so that PΛP =0. This also means that P and Λ do not commute. ⊥ 6 Equations (57) and (58) for the functions Pδχk and P⊥δχk are coupled. Nevertheless, we shall see that, under certain conditions, we can decouple the equation for Pδχk in the long time limit. If we solve formally equation (58), we obtain τ P⊥δχk(c,τ)=G0(τ)P⊥δχk(c,0) dτ′Gτ′(τ τ′)P⊥ilH(τ′)k cPδχk(c,τ′), (60) − − · Z0 where we have introduced the operator Gτ′(τ τ′) defined from − d Gτ′(τ τ′)=P⊥[Λ(c) ilH(τ)k c]P⊥Gτ′(τ τ′), Gτ′(0)=1. (61) dτ − − · − If the hydrodynamiceigenvalues of the operatorΛ are separatedenough fromthe rest of the spectrum, the firstterm on the right hand side of (60) decays with the “non hydrodynamic” modes, faster than the second one. We can then write τ P⊥δχk(c,τ) dτ′Gτ′(τ τ′)P⊥ilH(τ τ′)k cPδχk(c,τ τ′). (62) ≈− − − · − Z0 8 and we see, by substituting equation (62) in (57), that we obtain an involved but closed equation for Pδχk. It is worth pointing out that we have not proved scale separation, but assumed it in order to derive (62). For an explicit discussion of the scale separation assumption in a similar but somewhat simplified context, we refer to Appendix C, already alluded to above. Thesetofhydrodynamicequations(31)-(33)havebeenobtainedthroughtheprojectionoftheBoltzmannequation ontothe hydrodynamicsubspace. Itnowappearsthat, inthe hydrodynamictime scale,the use ofEq. (62)willallow us to close these equations by substituting the distribution function by its decomposition in terms of the projectors, δχk =Pδχk+P⊥δχk. This is the aim of the next section. IV. LINEAR HYDRODYNAMIC EQUATIONS IN NAVIER-STOKES ORDER In this section we shall use the decomposition of δχk into its hydrodynamic part, Pδχk, and non-hydrodynamic part, P⊥δχk, to close the linear hydrodynamic equations (31)-(33). We shall do so in Navier-Stokes order, that is, in the long time limit and in second order in the gradients (order k2). Let us first introduce Pδχk in the linear pressure tensor and in the heat flux vector. Here the calculation is straightforwardand we obtain Π[Pδχk]=0, φ[Pδχk]=0, (63) because the functions χ (c)∆(c) and χ (c)Σ(c) are orthogonal to the subspace spanned by the hydrodynamic H H eigenfunctions ξ (c) . Turning our attention to the other functionals δζ , δζ and δζ , the calculations β β=1,...3 n u T { } become somewhat lengthy, and we show the details in the Appendix D. We obtain δζn[Pδχk] = 4ζnρk ζnθk, (64) − − δζu[Pδχk] = 2ζnwk, (65) − δζT[Pδχk] = 4ζTρk (3ζT +2ζn)θk. (66) − − The negative signs occurring on the right hand side of these relations account for the fact that a fluctuation with a local enhanced density will induce an increased collision rate, hence a faster density decay. The same remark holds for temperature or local velocity flow fluctuations. WenowhavetocalculatethecontributionofP⊥δχk tothesamefunctionals,tosecondorderink. Thisrequiresthe knowledge of P⊥δχk to first order in k since the heat flux and pressure tensor enter the balance equations (31)-(33) throughtheir gradients and are already weightedby a factork. However,it shouldbe noted that for consistency, the decay rates should be computed to second order in the gradients (see Eqs. (31)-(33)). We shall nevertheless restrict to first order, henceforth neglecting the various terms of order two that symmetry allows (such as 2n and 2T for δζ and δζ , or as 2u for δζ ). We will further comment this approximationbelow. To leading or∇der,we ha∇ve that n T u Gτ−τ′(τ′) eP⊥ΛP⊥∇τ′, so that we obtain from equation (62) ≈ P⊥δχk(c,τ) ≈ τ dτ′eP⊥ΛP⊥τ′P⊥ilH(τ τ′)k cPδχk(c,τ τ′) ≈− − · − Z0 τ lH(τ) dτ′eP⊥ΛP⊥τ′e−2pζnτ′P⊥ik cPδχk(c,τ τ′), (67) ≈− · − Z0 where we have used that lH(τ) e2pζnτ (Eq. (27)). We now have to relate Pδχk(τ τ′) to Pδχk(τ), and to be ∝ − consistentwith the approximationmade above,we also haveto restrictto leading orderin k. In doing so, Markovian equations for the fields will be derived. From Eq. (57), we get 3 Pδχk(c,τ τ′) e−PΛPτ′Pδχk(c,τ)= e−λβτ′ ξ¯β δχk(τ) ξβ(c). (68) − ≈ h | i β=1 X Substituting (68) in (67), we obtain an equation for P⊥δχk to first order in k 3 τ P⊥δχ(k1)(c,τ)= lH(τ) ξ¯β δχk(τ) dτ′eP⊥(Λ−2pζn−λβ)P⊥τ′P⊥ik cξβ(c), (69) − h | i · β=1 Z0 X 9 where P⊥δχk = P⊥δχ(k1) + (k2). The pressure tensor and the heat flux up to first order in the gradients of the O fieldsarenowobtainedbysubstitutingequation(69)intoequations(34)and(35). Takingintoaccountthesymmetry properties of the system, the resulting expressions can be written in the form 2 Πij[P⊥δχ(k1)] = ilH(τ)η˜(τ) kjwi,k+kiwj,k+ k wkδij , (70) − d · (cid:20) (cid:21) φ[P⊥δχ(k1)] = ilH(τ)k[κ˜(τ)θk+µ˜(τ)ρk]. (71) − Equation(70)istheexpectedNavier-Stokesexpressionforthepressuretensor,involvingtheshearviscositycoefficient η˜, but equation (71) contains, besides the usual Fourier law characterized by the heat conductivity κ˜, an additional contribution proportional to the density gradient and with an associated transport coefficient µ˜. This latter term is analogous to the one appearing in granular gases [3, 18]. The expression of the (time dependent) transport coefficients are d 1 η˜(τ)= dc∆ (c)F (c,τ)= dc∆ (c)F (c,τ), (72) xy 3,xy d2+d 2 ij 3,ij Z − i,j Z X 1 µ˜(τ)= dcΣ(c)[F (c,τ)+F (c,τ)], (73) 1 2 d(1+z) Z 1 κ˜(τ)= dcΣ(c)[ zF (c,τ)+F (c,τ)], (74) 1 2 2d(1+z) − Z where we have introduced the functions τ F (c,τ) = dτ′eP⊥(Λ−2pζn−pζT)τ′P c ξ (c), (75) 3,ij ⊥ i 3,j Z0 τ F (c,τ) = dτ′eP⊥(Λ−2pζn)τ′P cξ (c), (76) 1 ⊥ 1 Z0 τ F (c,τ) = dτ′eP⊥(Λ+pζT)τ′P cξ (c), (77) 2 ⊥ 2 Z0 and in the second equality of equation (72), we have summed over all the i, j, taking into account the symmetry of the linearized Boltzmann operator. Similarly, we calculate the deviations of the decay rates to first order in the gradients of the fields by substituting equation (69) into equations (38), (39) and (40). Taking into account the symmetry properties, we arrive at δζ [P δχ(1)] = 0, (78) n ⊥ k δζu[P⊥δχ(k1)] = ilH(τ)k[ζu,ρ(τ)ρk+ζu,θ(τ)θk], (79) δζ [P δχ(1)] = 0. (80) T ⊥ k The expression for the coefficients are γβ ζ (τ)= dc dc χ (c )c (c +c ) [F (c ,τ)+F (c ,τ)], (81) u,ρ 1 2 H 1 12 1 2 1 2 2 2 d(1+z) · Z Z γβ ζ (τ)= dc dc χ (c )c (c +c ) [ zF (c ,τ)+F (c ,τ)], (82) u,θ 1 2 H 1 12 1 2 1 2 2 2 2d(1+z) · − Z Z with β =π(d−1)/2/Γ[(d+1)/2] the d-dimensional solid angle. At this point, it is important to note that the transport coefficients defined in equations (72)-(74) and (81)-(82) are time-dependent, and this dependence is governed by the F functions. The (exponential) integrands appearing i in the definitions of the F decay with the non-hydrodynamic (kinetic) modes, as a consequence of the action of i the projector P . From our assumption of hydrodynamic versus kinetic scale separation, all kinetic eigenvalues are ⊥ smallerthan the smallesthydrodynamic eigenvalueofΛ, whichis λ = pζ 2pζ . This ensures the convergenceof 2 T n − − theintegrals(75)-(77)forτ . Inorderforthe transportcoefficientstoreachtheirτ limitfasterthananyof →∞ →∞ 10 thehydrodynamictime scales,weneedmoreoverthe morestringentconditionthatthe fastestkineticmodeisatleast separated by a pζ gap from λ = pζ 2pζ : under this condition, the time dependence of the exponential term T 2 T n in the integral giving F is fast eno−ugh s−o that the transport coefficients, that depend on the F functions through 2 i (72)-(74), can be considered as constants on the hydrodynamic time scale. With this proviso in mind, it is possible to set τ in the integrals (75)-(77) and the time-independent transport coefficients obtained in this section are → ∞ then equivalent to those calculated in reference [12] by the Chapman-Enskogmethod. We recallin Appendix A their expressions in the first order Sonine approximation. Finally, if we substitute the expressions derived above for the fluxes, equations (70)-(71), and the decay rates, equations(78)-(80),andwetakeintoaccountthatinthehydrodynamictimescalewecansubstituteallthecoefficients by their τ limit, we obtain the following closed equations for the linear deviation of the hydrodynamic fields →∞ ∂ ∂τ +2pζn ρk+ilH(τ)kwk||+pζnθk =0, (83) (cid:18) (cid:19) ∂ ∂τ −pζT +lH2 (τ)η˜k2 wk⊥ =0, (84) (cid:20) (cid:21) ∂ 2(d 1) ∂τ −pζT + d− lH2 (τ)η˜k2 wk|| (cid:20) (cid:21) i + lH(τ)k[(1 2pζu,ρ)ρk+(1 2pζu,θ)θk]=0, (85) 2 − − ∂ 2 ∂τ +pζT + dlH2 (τ)κ˜k2 θk (cid:20) (cid:21) 2 2 + 2pζT + dlH2 (τ)µ˜k2 ρk+idlH(τ)kwk|| =0, (86) (cid:20) (cid:21) where wk|| and wk⊥ are the longitudinal and transversal parts of the velocity vector defined by wk|| =wk kˆ, wk⊥ =wk wk||kˆ, (87) · − and kˆ is the unit vector along the direction given by k. Equation (84) for the shear mode is decoupled from the other equations and can be readily integrated. If we introduce a non-dimensionalwavenumber k˜ =l (0)k, scaledby the meanfree path atthe time origin,we obtainthe H explicit solution η˜k˜2 wk⊥(τ)=exp pζTτ e4pζnτ 1 wk⊥(0). (88) " − 4pζn − # (cid:0) (cid:1) Interestingly, depending on k˜, the perturbation may initially increase if pζ η˜k˜2 > 0. For long times however, the T − exponential e4pζnτ always dominates the linear term pζTτ and the perturbation decays. The other three fields, namely, the density ρk, temperature θk, and the longitudinal velocity wk||, obey the system of coupled linear equations ∂ ρk ρk wk|| =M(τ) wk|| , (89) ∂τ   ·  θk θk     where the time-dependent matrix is 2pζ il (τ)k pζ n H n − − − M(τ)= il (τ)k(1 2pζ ) pζ 2(d−1)l2 (τ)η˜k2 il (τ)k(1 2pζ ) . −2 H − u,ρ T − d H −2 H − u,θ  2pζ 2l2 (τ)µ˜k2 i2l (τ)k pζ 2l2 (τ)κ˜k2 − T − d H − d H − T − d H   (90) We note here that this matrix differs from Eq. (59) of Ref. [12], where the analysis amounts to overlooking the time dependence of the mean free path, so that all entries of the hydrodynamic matrix exhibit the same time dependence. ThedifferenttimedependencespresentinEq. (90)renderthestabilityanalysismoredifficult. Forlongtimeshowever, the eigenvaluesofthe matrix M(τ)arealwaysnegativeandthe perturbationsa priori decay. Acaveat isnevertheless