Global stability and H-theorem in lattice models with non-conservative interactions
GGlobal stability and H -theorem in lattice models with non-conservative interactions C. A. Plata ∗ and A. Prados † F´ısica Te´orica, Universidad de Sevilla, Apdo. de Correos 1065, Sevilla 41080, Spain (Dated: April 17, 2018)In kinetic theory, a system is usually described by its one-particle distribution function f ( r , v , t ),such that f ( r , v , t ) d r d v is the fraction of particles with positions and velocities in the intervals( r , r + d r ) and ( v , v + d v ), respectively. Therein, global stability and the possible existence of anassociated Lyapunov function or H -theorem are open problems when non-conservative interactionsare present, as in granular fluids. Here, we address this issue in the framework of a lattice model forgranular-like velocity fields. For a quite general driving mechanism, including both boundary andbulk driving, we show that the steady state reached by the system in the long time limit is globallystable. This is done by proving analytically that a certain H -functional is non-increasing in the longtime limit. Moreover, for a quite general energy injection mechanism, we are able to demonstratethat the proposed H -functional is non-increasing for all times. Also, we put forward a proof thatclearly illustrates why the “classical” Boltzmann functional H B [ f ] = (cid:82) d r d v f ( r , v , t ) ln f ( r , v , t ) isinadequate for systems with non-conservative interactions. Not only is this done for the simplifiedkinetic description that holds in the lattice models analysed here but also for a general kineticequation, like Boltzmann’s or Enskog’s. I. INTRODUCTION
In thermodynamics and statistical mechanics, globalstability of the equilibrium state is usually proven byintroducing a Lyapunov functional [1]. This Lyapunovfunctional of the probability distribution function (PDF)has the following three properties: (i) it is bounded frombelow, (ii) it monotonically decreases with time and (iii)its time derivative equals zero only when the PDF is theequilibrium one. Therefore, in the long time limit, theLyapunov functional must tend to a finite value and thusits time derivative vanishes. As a consequence, any PDF,corresponding to an arbitrary initial preparation, tendsto the equilibrium PDF: the equilibrium state is irre-versibly approached and said to be globally stable.The first example of such a Lyapunov functional isthe renowned Boltzmann H -functional. In the Boltz-mann description, the nonequilibrium behaviour of adilute gas is completely encoded in the one-particlevelocity distribution function f ( r , v , t ). By introduc-ing the Stosszahlansatz or Molecular Chaos hypothesis,Boltzmann derived a closed non-linear integro-differentialequation for f ( r , v , t ) governing its time evolution [2].Also, for a spatially homogeneous state, he showed thatthe functional H B [ f ] = (cid:82) d v f ( v , t ) ln f ( v , t ) has the threeproperties of a Lyapunov functional. This H -theoremshows that all solutions of the Boltzmann equation tendin the long time limit to the Maxwell velocity distribu-tion and irreversibility naturally stems from a molecu-lar picture [3, 4]. Interestingly, a key point for derivingthe H -theorem is the reversibility of the underlying mi-croscopic dynamics. In an inhomogeneous situation, onehas to consider the spatial dependence of the one-particledistribution function f ( r , v , t ), and the above functional ∗ [email protected] † [email protected] must be generalised to H B [ f ] = (cid:90) d r d v f ( r , v , t ) ln f ( r , v , t ) . (1)Provided that the walls of the gas container are smooth,in the sense that there is no energy transport throughthem, it can be also shown that this is a non-increasingLyapunov functional [5].Another example of a Lyapunov functional canbe found in the realm of Markovian stochastic pro-cesses. Therein, the stochastic process X ( t ) is com-pletely determined by the conditional probability den-sity P | ( X, t | X , t ) of finding the system in state X at time t , given it was in state X at time t , andthe probability density P ( X, t ) of finding the systemin state X at time t [6]. Both probability densitiessatisfy the same evolution equation, named the masterequation, but with different initial conditions: one al-ways has that P | ( X, t | X , t ) = δ ( X − X ), whereas P ( X, t ) = P ini ( X ), with P ini ( X ) corresponding to the(arbitrary) initial preparation. When the stochastic pro-cess is irreducible or ergodic, that is, every state can bereached from any other state by a chain of transitionswith non-zero probability, there is only one stationarysolution of the master equation. In physical systems,this steady solution must correspond to the equilibrium-statistical-mechanics distribution P eq ( X ). What is more,a Lyapunov functional can be constructed in the follow-ing way, H [ P ] = (cid:90) dXP eq ( X ) g (cid:20) P ( X, t ) P eq ( X ) (cid:21) , (2)where g ( x ) is any positive-definite convex function( g (cid:48)(cid:48) ( x ) ≥ H -theorem for master equations rely only on the ergodicityof the underlying microscopic dynamics: it is not neces-sary to assume that detailed balance, which is connectedwith the microscopic reversibility, holds [6]. a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r The most usual choice for g is g ( x ) = x ln x − x + 1,which leads to H [ P ] = (cid:90) dXP ( X, t ) ln (cid:20) P ( X, t ) P eq ( X ) (cid:21) . (3)The physical reason behind this choice is the “extensive-ness” of H [ P ]: if the system at hand comprises two inde-pendent subsystems A and B , so that dX ≡ dX A dX B and P ( X ) = P A ( X A ) P B ( X B ), one has that H [ P ] = H A [ P A ] + H B [ P B ]. It is to consider −H as a nonequi-librium entropy S that this extensiveness is desirable: inthis way, the non-increasing behaviour of H leads to anon-decreasing time evolution of S . Moreover, H [ P ] re-mains invariant upon a change of variables Y = f ( X ), asemphasised in Refs. [7, 8].Although the Boltzmann equation is not a masterequation, we may wonder why the expressions for H B in Eq. (1) and H [ P ] in Eq. (3) are different. Specifically,we may wonder why not writing H [ f ] = (cid:90) d r d v f ( r , v , t ) ln (cid:20) f ( r , v , t ) f eq ( v ) (cid:21) (4)instead of H B [ f ]. Up to now, we have been implic-itly considering the “classic” problem with elastic col-lisions between particles, in which the system eventuallyreaches thermodynamic equilibrium. Therein, the answeris trivial: since ln f eq ( v ) is a sum of constants of motion, H [ f ] − H B [ f ] is constant and both are utterly equivalent.Whether there exists an extensive H -functional or notis an important question in nonequilibrium statisticalphysics. If the answer were positive, it would make itpossible to define a non-equilibrium entropy − H thatmonotonically grows for all times, extending the Clausiusinequality. In general, the system at hand does not reachequilibrium but a nonequilibrium steady state. Thus, theequilibrium distribution f eq in H has to be substitutedwith the stationary one f st . In this context, the fieldof granular fluids is a benchmark for intrinsically out-of-equilibrium, dissipative, systems: the microscopic dy-namics is not time-reversible because collisions betweenparticles are inelastic, but a nonequilibrium steady statecan be attained if some driving mechanism injects energyinto the system.In granular fluids, the functionals H [ f ] and H B [ f ] areno longer equivalent, since ln f st is not a sum of constantsof motion. Indeed, for granular gases described by the in-elastic Boltzmann equation [9, 10], there are some resultsthat hint at H B not being a Lyapunov functional. Withinthe first Sonine approximation, it has been proven thatthe time derivative of H B does not have a definite signin the linear approximation around the steady state [11].Moreover, Marconi et al. have numerically shown that H B is non-monotonic and even steadily increases fromcertain initial conditions [7]. They have also put forwardsome numerical evidence (further reinforced by Garc´ıade Soria et al. [8]) in favour of H being a “good” Lya-punov functional. Notwithstanding, only spatially homo-geneous situations, in which the r -dependence of f and thus the integration over r may be dropped, have beenanalysed in Refs. [7, 8].Some years ago, a simplified model for a granular-like velocity field was introduced to study correlations ingranular gases [12]. Very recently, a variant of this modelon a one-dimensional lattice has been proposed to mimicthe velocity component along the shear direction [13],and both its hydrodynamic limit and finite size effectshave been analysed [13–15]. This model has been shownto retain a relevant part of the granular phenomenology:the shear instability of the homogeneous cooling state,the existence of boundary driven steady states such as theCouette and Uniform Shear Flow (USF) states, the renor-malisation of the cooling rate due to fluctuations close tothe shear instability, etc. Other properties thereof, whenit is driven by a mechanism resembling collisions witha randomly moving inelastic wall, have been studied inRef. [16]. At the N -particle level, the dynamics of thesystem is governed by a master equation, which is anal-ogous to the Kac equation [17], that leads to a “kinetic”equation at the one-particle level, which is analogous tothe Boltzmann equation. In the latter, the collision term,although being simpler than that in the Boltzmann equa-tion, remains a non-linear integro-differential one [14].It must be recalled that an analytical proof of eitherglobal stability or the H -theorem is currently unavailableat the level of the kinetic description for granular gases.This is true even for simple collision terms, such as thosecorresponding to hard-spheres or the cruder Maxwell par-ticle model (where the collision rate is considered to bevelocity-independent), which are considered in the pio-neering work in Refs. [7, 8]. Therefore, it seems worthinvestigating this subject in simplified models, for whichanalytical calculations are more feasible.Our main goal here is to investigate the global stabilityand the possibly associated H -theorem in the above classof lattice models. Unlike the approach in Refs. [7, 8], wedo not restrict ourselves to spatially homogeneous situ-ations but consider the whole space and velocity depen-dence of the one-particle PDF f ( r , v , t ). Specifically, weintroduce a general energy injection mechanism, in whichthe system may be driven both through the boundariesand in the bulk. We show that, under quite general con-ditions, the steady state is globally stable: independentlyof the initial preparation, the system always ends up inthe steady state. Interestingly, it is not necessary to havean H -theorem to prove this: it suffices to show that H is decreasing in the long time limit, not for all times.In this sense, the situation is analogous to the proof ofthe tendency towards the equilibrium curve in systemswhose dynamics is governed by master equations withtime-dependent transition rates [18–24].Our proof of global stability also enables us to show theinadequacy of Boltzmann’s H B as a candidate for Lya-punov functional in inelastic systems. Not only is thisdone for the simplified models considered in the paper,but for a general collision term that does not conserveenergy in collisions. Therefore, this result also appliesto the inelastic Boltzmann or Enskog equations used ingranular fluids. The main idea is that the sign of dH B /dt can be reversed by a suitable choice of the initial PDF,and thus cannot have a definite sign. In this respect,our result generalises that in Ref. [11], which was de-rived within the first Sonine approximation of the inelas-tic Boltzmann equation, to an arbitrary collision kernelwith non-conservative interactions.Having proved global stability by showing that H is anon-increasing functional for long times, a natural ques-tion remains. Is it H a Lyapunov function, that is, anon-increasing functional for all times? There does notseem to be a unique proof, valid for any driving mecha-nism, even within our simplified model. Nevertheless, wehave been able to derive a specific proof for a quite gen-eral driving mechanism, which includes as limiting casesboth the sheared system, in which the steady state is theUSF state, and the uniformly heated system by meansof the so-called stochastic thermostat [25–33]. The proofis based on a suitable expansion of the one-particle PDFin Hermite polynomials, which is a generalisation of theusual Sonine expansion of kinetic theory.The paper is organised as follows. In Sec. II, we brieflyintroduce the model, its dynamics and the continuumlimit. Section III is devoted to the proof of the global sta-bility of the nonequilibrium steady states, for a generalenergy injection mechanism. The inadequacy of Boltz-mann’s H B as a Lyapunov functional for inelastic systemsis discussed in Sec. IV. Later, in Sec. V, we consider someconcrete physical situations in our model, which includethe sheared and the uniformly heated systems. Therein,we show that H [ f ] is a monotonically decreasing Lya-punov functional. Finally, Sec. VI gives the main con-clusions of the paper. Some technical details, which areomitted in the main text, are given in the Appendices. II. THE MODEL: DYNAMICS ANDCONTINUUM LIMIT
Here, we present the general class of models that wasintroduced in Ref. [13], focusing on the continuum de-scription obtained in the large system size limit [14].Specifically, our system is defined on a 1d lattice: at eachlattice site l , there is a particle with velocity v l . Thus, ata given time τ , the configuration of the system is com-pletely determined by v ≡ { v , ..., v N } . The dynamicsproceeds through inelastic nearest-neighbour binary col-lisions: each pair ( l, l + 1) collides inelastically with acharacteristic rate ω − , independently of their relativevelocity (the so-called Maxwell-molecule model [34]) andthe state of the other pairs. We introduce the operatorˆ b l that transforms the pre-collisional velocities into thepost-collisional ones,ˆ b l v l = v l − α ( v l − v l +1 ) , (5a)ˆ b l v l +1 = v l +1 + α ( v l − v l +1 ) , (5b) where α is the normal restitution coefficient, with 0 <α ≤ stochastic thermostat [25–33]. Specifically, fora short time interval, the change of the velocity due tothe heating is given by∆ v i ( τ ) | noise ≡ v i ( τ + ∆ τ ) − v i ( τ ) | noise = ξ i ( τ ) − N N (cid:88) j =1 ξ j ( τ ) ∆ τ, (6)where ξ i ( t ) are Gaussian white noises, verifying (cid:104) ξ i ( τ ) (cid:105) noise = 0 , (cid:104) ξ i ( τ ) ξ j ( τ (cid:48) ) (cid:105) noise = χδ ij δ ( τ − τ (cid:48) ) , (7)for i, j = 1 , . . . , N . Above, χ is the amplitude of thenoise, and (cid:104)· · ·(cid:105) noise denotes the average over the differ-ent realisations of the noise. Note that this version ofthe stochastic thermostat conserves total momentum, anecessary condition to have a steady state [30, 33].We define P N ( v , τ ) as the probability density of findingthe system in state v at time τ . The stochastic process v ( τ ) is Markovian and the equation governing the timeevolution of P N ( v , τ ) has two contributions. First, wehave a master equation contribution stemming from col-lisions [14, 15] ∂ τ P N ( v , τ ) | coll = ω N (cid:88) l =1 (cid:34) P N (ˆ b − l v , τ ) α − P N ( v , τ ) (cid:35) , (8)in which the operator ˆ b − l is the inverse of ˆ b l , that is,it changes the post-collisional velocities into the pre-collisional ones when the colliding pair is ( l, l + 1). Sec-ond, there is a Fokker-Planck contribution stemmingfrom the stochastic forcing [7, 8] [35] ∂ τ P N ( v , τ ) | noise = χ N (cid:88) i,j =1 (cid:18) δ ij − N (cid:19) ∂ ∂v i ∂v j P N ( v , τ ) . (9)The time evolution of P N ( v , τ ) is obtained by combiningEqs. (8) and (9), that is, ∂ τ P N ( v , τ ) = ∂ τ P N ( v , τ ) | coll + ∂ τ P N ( v , τ ) | noise . (10)In this work, we focus on the evolution of quantitiesthat can be written in terms of the one-particle distribu-tion function, namely P ( v ; l, τ ) = (cid:90) d v P N ( v , τ ) δ ( v l − v ) . (11)All the one-site velocity moments can be calculated from P , (cid:104) v nl ( τ ) (cid:105) ≡ (cid:90) + ∞−∞ dv v n P ( v ; l, τ ) . (12)The first two moments give the hydrodynamic fields: theaverage velocity u l ( τ ) and granular temperature T l ( τ )[36], which are defined by the relations u l ( τ ) ≡ (cid:104) v l (cid:105) , T l ( τ ) ≡ (cid:104) v l ( τ ) (cid:105) − u l ( τ ) . (13)Here, we do not write the evolution equations on the lat-tice for either P or the hydrodynamic fields ( u and T ),since they are not necessary for our present purposes.The unforced case ( χ = 0) can be found in Ref. [14].However, we would like to stress that the evolution equa-tion for P is not closed, since the collision term involvesthe two-particle distribution function P ( v, v (cid:48) ; l, l + 1 , τ ).As usual in kinetic theory, one can write a closed equationfor P after introducing the Molecular Chaos assumption,that is, P ( v, v (cid:48) ; l, l + 1 , τ ) = P ( v ; l, τ ) P ( v (cid:48) ; l + 1 , τ ) + O ( N − ). In other words, one assumes that the correla-tions at different sites are of the order of N − and thusnegligible in the large system size limit.The continuum limit of the model is introduced forlarge system size N (cid:29)
1, in which we expect the aver-age velocity u l and temperature T l to be smooth func-tions of space and time. This is expressed mathemati-cally by defining “hydrodynamic” continuous space andtime variables by x = l/N and t = ωτ /N , respectively[14]. Note that 0 ≤ x ≤ t ≥
0. In the continuumlimit, the one-particle distribution function also becomesa smooth function of x and τ , P ( v ; x, t ) ≡ P ( v ; l = N x, τ = N t/ω ).From now on, we use the usual notation in kinetictheory f ( x, v, t ) ≡ P ( v ; x, t ). The physical picture isstraightforward: f ( x, v, t ) dxdv gives the fraction of thetotal number of particles with positions between x and x + dx and velocities between v and v + dv . We have that (cid:82) + ∞−∞ dv f ( x, v, t ) = 1 for all x and t , since there is nomass transport in the system. The time evolution of f isgoverned by the non-linear integro-differential (pseudo-Boltzmann) equation [14] ∂ t f = ∂ x f + ν ∂ v { [ v − u ( x, t )] f } + ξ ∂ v f, (14)where u ( x, t ) is the local average velocity, ν is the macro-scopic dissipation coefficient and ξ is the macroscopicnoise strength, which are respectively given by ν = (1 − α ) N , ξ = χN ω . (15)This shows that the microscopic noise strength χ mustscale as N − in order to have a finite contribution inthe continuum limit. Of course, for ξ = 0, we recoverthe kinetic equation for the case in which there is nostochastic forcing, see Ref. [14]. The N -scaling of themacroscopic dissipation coefficient ν is similar to thatfound in the dissipative version of the Kipnis-Marchioro-Presutti model [37–41].The average velocity u ( x, t ) and granular temperature T ( x, t ) are the continuum limit of u l and T l defined in Eq. (13), u ( x, t ) ≡ (cid:104) v (cid:105) ( x, t ) , T ( x, t ) ≡ (cid:104) v (cid:105) ( x, t ) − u ( x, t ) , (16)where the velocity moments are given by (cid:104) v n (cid:105) ( x, t ) = (cid:82) dv v n f ( x, v, t ). From the kinetic equation for f ( x, v, t ),one can derive the evolution equations of u and T , ∂ t u = ∂ xx u, (17a) ∂ t T = − νT + ∂ x T + 2 ( ∂ x u ) + ξ. (17b)On the one hand, Eq. (17a) is a diffusion equation forthe average velocity, which expresses the conservationof total momentum. On the other hand, the tempera-ture equation (17b) contains a purely dissipative term − νT that stems from the inelasticity of collisions andalways contributes to “cooling” the system, a diffusiveterm ∂ x T , a viscous heating term 2 ( ∂ x u ) , and finallythe term corresponding to the uniform heating ξ . Ofcourse, either the kinetic equation for f or the averageequations for ( u, T ) must be complemented with suitableboundary conditions in each physical situation. A. Non-equilibrium steady states and boundaryconditions
We are interested in driven cases, in which there isan input of energy that balances (in average) the en-ergy loss in collisions, so that the system eventuallyreaches a steady state. These non-equilibrium steadystates (NESS) are described by the corresponding sta-tionary solutions f st ( x, v ) of the kinetic equation, whichverify 0 = ∂ x f st + ν ∂ v { [ v − u st ( x )] f st } + ξ ∂ v f st . (18)where u st ( x ) = (cid:82) dv vf st ( x, v ) is the stationary averagevelocity profile. To be concrete, we consider two cases: asystem that is (a) sheared and (b) uniformly heated.First, let us consider a sheared system: there is nostochastic forcing, ξ = 0, and the driving is introducedby imposing a velocity difference (“shear”) between theleft and right edges of the system. At the level of thehydrodynamic description, the corresponding boundaryconditions are u (1 , t ) = u (0 , t ) + a, u (cid:48) (1 , t ) = u (cid:48) (0 , t ) , (19a) T (1 , t ) = T (0 , t ) , T (cid:48) (1 , t ) = T (cid:48) (0 , t ) , (19b)which are said to be of Lees-Edwards type [42]. We haveused (cid:48) to denote spatial derivative. The imposed shearallows the viscous heating term, which is proportional to( ∂ x u ) , to compensate for the energy dissipation term, − νT . The boundary conditions for the one-particle dis-tribution function read f (1 , v, t ) = f (0 , v − a, t ) , f (cid:48) (1 , v, t ) = f (cid:48) (0 , v − a, t ) , (20a)from which Eq. (19) directly follow. Equation (20) hasa simple physical interpretation: particles that leave thesystem through its right edge with velocity v are rein-serted through its left edge with velocity v − a .The steady state for the sheared system is known asthe USF state, which has a linear velocity profile and ahomogeneous temperature, u st ( x ) = a (cid:18) x − (cid:19) , T st = 2 a ν . (21)For our simplified model, the stationary PDF is Gaus-sian, f st ( x, v ) = (2 πT st ) − / exp (cid:20) − ( v − u st ( x )) T st (cid:21) . (22)An extensive investigation of the sheared system, at thelevel of the average hydrodynamic equations, can befound in Ref. [14].Second, we address the uniformly heated system, inwhich there is no shear, a = 0, but there is stochasticforcing, ξ (cid:54) = 0. In this case, we have the usual periodicboundary conditions. In particular, for the PDF we have f (1 , v, t ) = f (0 , v, t ) , f (cid:48) (1 , v, t ) = f (cid:48) (0 , v, t ) . (23)In the steady state, the system is homogeneous: there isno average velocity and the temperature is uniform, u st ( x ) = 0 , T st = ξν . (24)The corresponding stationary PDF is also Gaussian, f st ( (cid:1) x, v ) = (2 πT st ) − / exp (cid:20) − v T st (cid:21) . (25)With this “stochastic thermostat” forcing, the systemremains homogeneous for all times if it is initially so, as isalso the case of a inelastic gas of hard particles describedby the inelastic Boltzmann equation [25]. III. GLOBAL STABILITY
In this section, we analyse the global stability of thenonequilibrium stationary solutions of the kinetic equa-tion (14) submitted to quite a general class of boundaryconditions. Following the discussion in the introduction,we define the H -functional as H [ f ] = (cid:90) dx dvf ( x, v, t ) ln (cid:20) f ( x, v, t ) f st ( x, v ) (cid:21) . (26)Let us consider the time evolution of H [ f ]. It is directlyobtained that dHdt = (cid:90) dx dv ∂ t f ln (cid:18) ff st (cid:19) = (cid:90) dx dv L f ln (cid:18) ff st (cid:19) , (27) where L stands for the nonlinear evolution operator onthe rhs of the kinetic equation (14), that is, ∂ t f = L f .Now we note the following property: if we define ∆ f = f − f st to be the deviation of the PDF from the steadystate, the linear terms in the deviation vanishes, sinceboth factors in the integrand of (27) are equal to zero for f = f st . This is a desirable property: were it not true,the sign of dH/dt could be reversed for initial conditionsclose enough to the steady state by simply reversing theinitial value of ∆ f . Thus, the existence of an H -theoremwould be utterly impossible, see also next Section.Then, we can write dHdt = (cid:90) dx dv L f ln (cid:18) ff st (cid:19) − (cid:90) dx dv L f st f − f st f st . (28)Now, the idea is to split the operator L into the threecontributions on the rhs of Eq. (14): first, the diffusiveone; second, the one proportional to ν , which is intrinsi-cally dissipative; and third, the one proportional to thenoise strength ξ : L diff , L inel and L noise , respectively. Ac-cordingly, we have that the time derivative of H has threecontributions, dHdt = dHdt (cid:12)(cid:12)(cid:12)(cid:12) diff + dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel + dHdt (cid:12)(cid:12)(cid:12)(cid:12) noise , (29)obtained by inserting into Eq. (28) the relevant part ofthe evolution operator L . Note that, although L f st = 0,in general L diff f st (cid:54) = 0, L inel f st (cid:54) = 0 and L noise f st (cid:54) = 0.After some tedious but easy algebra, a summary ofwhich is given in Appendix A, the following expressionsare derived. Firstly, for the diffusive term, dHdt (cid:12)(cid:12)(cid:12)(cid:12) diff = − (cid:90) dx dv f ( ∂ x ln f − ∂ x ln f st ) ≤ . (30)Secondly, for the inelastic term, proportional to ν , dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel = − ν (cid:90) dx ( u − u st ) (cid:90) dvf ∂ v ln f st . (31)Finally, the noise term, proportional to ξ , reads dHdt (cid:12)(cid:12)(cid:12)(cid:12) noise = − ξ (cid:90) dx dv f ( ∂ v ln f − ∂ v ln f st ) ≤ . (32)These results, and the following throughout this section,are valid for a quite general set of boundary conditions,leading to the cancellation of all the boundary terms aris-ing after integrating by parts, as detailed in Appendix A.This set includes but is not limited to the Lees-Edwardsand periodic boundary conditions corresponding to thesheared and uniformly heated situations, respectively.For instance, they also apply to the Couette state, inwhich the system is driven by keeping its two edges attwo (in general, different) fixed temperatures T L and T R .The inelastic term dH/dt | inel in Eq. (31) does not havea definite sign in general. Therefore, it is the inelas-tic term that prevents us from proving H to be a non-increasing function of time. It must be stressed that thediffusive, inelastic and noise contributions to dH/dt inEqs. (30)-(31) come exclusively from the diffusive, noiseand inelastic contributions in the kinetic equation, re-spectively, only once the linear terms has been subtractedas is done in Eq. (28): see Appendix A for details.Despite the above discussion, global stability of thesteady state can be established without proving an H -theorem. The key point is the following: the long timelimit of dH/dt is non-positive and thus H has a fi-nite limit, since it is bounded from below. Therefore, dH/dt tends to zero in the long time limit and it canbe shown that this is only the case if f ( x, v, ∞ ) ≡ lim t →∞ f ( x, v, t ) = f st ( x, v ).The average velocity u ( x, t ) satisfies a diffusive equa-tion (17a), and thus it irreversibly tends to the steadyprofile corresponding to the given boundary condi-tions in the long time limit. Therefore, u ( x, ∞ ) ≡ lim t →∞ u ( x, t ) = u st ( x ) and taking into account Eq. (31),lim t →∞ dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel = 0 ⇒ lim t →∞ dHdt ≤ . (33)Since H [ f ] is bounded from below, the only possibility islim t →∞ dHdt = 0 , (34)and all the contributions to dH/dt in Eqs. (30)-(32) van-ish in the long time limit. The vanishing of Eq. (30)imposes that f ( x, v, ∞ ) = f st ( x, v ) φ ( v ), where φ ( v ) is anarbitrary function of v . For ξ (cid:54) = 0, Eq. (32) implies that φ ( v ) must be a constant, independent of v , and normal-isation yields φ ( v ) = 1. For ξ = 0, Eq. (32) identicallyvanishes but it can be also shown that φ ( v ) = 1 by usingthe kinetic equation in the limit as t → ∞ . Therefore,for arbitrary ξ , including ξ = 0, we have that f ( x, v, ∞ ) = f st ( x, v ) . (35)This completes the proof. The steady distribution f st ( x, v ) is globally stable: each time evolution f ( x, v, t )(corresponding to a given initial condition) tends to it inthe long time limit. IV. INADEQUACY OF H B AS A LYAPUNOVFUNCTIONAL FOR NON-CONSERVATIVESYSTEMS
Here we show that Boltzmann’s H B [ f ] cannot be usedto build a Lyapunov functional for intrinsically dissipa-tive systems, in agreement with the numerical results byMarconi et al. [7]. Not only do we prove it for the sim-plified models considered here, but for a general kineticequation in which energy is not conserved in collisions,such as the inelastic Boltzmann or Enskog equations. Tokeep the notation simple, we still write ∂ t f = L f , butnow L stands for the evolution operator in the consid-ered kinetic description, which is nonlinear in general. First, we restrict ourselves to homogeneous situationsand thus drop the integral over x , H B [ f ] = (cid:90) dvf ln f (36a) dH B dt = (cid:90) dv ∂ t f ln f = (cid:90) dv L f ln f, (36b)Also, we consider a system that is initially close to thesteady state, such that we can expand everything in pow-ers of ∆ f = f − f st . Then, L f ≡ L ( f st + ∆ f ) = (cid:8)(cid:8)(cid:8)(cid:42) L f st + L lin ∆ f + O (∆ f ) , (37)in which L lin is the linearised evolution operator. Ne-glecting O (∆ f ) terms, the linear approximation arises, dH B dt (cid:12)(cid:12)(cid:12)(cid:12) lin = (cid:90) dv ( L lin ∆ f ) ln f st = ddt (cid:104) ln f st (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) lin . (38)On the one hand, the linear contribution vanishes in theelastic case: ln f st is a sum of constants of motion, whichare unchanged by the linearised kinetic operator. Then, H B can be a candidate for a Lyapunov functional. Onthe other hand, only mass and linear momentum are con-served for non-conservative interactions. Thus, no longeris ln f st a sum of conserved quantities, and dH B dt (cid:12)(cid:12)(cid:12)(cid:12) lin (cid:54) = 0 . (39)Therefore, by changing the initial sign of ∆ f = f − f st ,which can always be done, the initial sign of dH B /dt isreversed and H B cannot be a Lyapunov functional.In Fig. 1, we show the evolution of H B in our ki-netic model. We consider a uniformly heated system,so that the system remain homogeneous for all times,as described in Sec. II A. Two different initial conditionsare considered, corresponding to Gaussian distributionswith zero average velocity but non-steady values of thetemperature, specifically 1 . T st and 0 . T st . We can seehow, in agreement with our discussion, not only is oneof the functionals increasing, but also it can be obtainedas the mirror image of the decreasing one through thestationary value. Technical details about the simulationare provided in Appendix B.Taking into account the specific (Gaussian) shape ofthe steady PDF for the uniformly heated system, as givenby Eq. (25), the time derivative of H B in Eq. (38) reducesto dH B dt = − T st d (cid:104) v (cid:105) dt . (40)Since the plots in Fig. 1 correspond to evolutions of thesystem for which u ( x, t ) ≡ (cid:104) v (cid:105) = T and, consistently, the H B -curve corresponding to aninitial value of the temperature that is higher (lower)than the steady one monotonically increases (decreases). FIG. 1. Evolution of the functional H B for two differ-ent initial conditions in a uniformly heated system. Bothsimulations start from the stationary Gaussian shape butwith a homogeneous temperature slightly shifted from thestationary one: (i) T ( t = 0) = 1 . T st (circles) and (ii) T ( t = 0) = 0 . T st (triangles). As predicted by the linearapproximation, both functionals are symmetric with respectto the stationary value, which is schematically illustratedby plotting their mean value (dashed line). A system with N = 330 sites has been considered, with ν = 20 and ξ = 50.The plots correspond to averages over 3000 trajectories. This is consistent with the situation found in Ref. [11],in which the uniformly heated granular gas described bythe inelastic Boltzmann equation was investigated withinthe first Sonine approximation. Therein, the entropy pro-duction was shown to have linear terms in the deviationsof the temperature and the excess kurtosis. Also, ourresult is consistent with the numerical results in Ref. [7]for several collision models. The above argument alsoproves why H B is not non-increasing for an elastic sys-tem immersed in a heat bath at a temperature differentfrom the initial temperature of the gas, as also observedin Ref. [7]. Although ln f st is conserved in collisions, theevolution operator includes a term coming from the in-teraction with the bath that does not conserve the kineticenergy, and again dH/dt | lin (cid:54) = 0, making it impossible for H B to be a Lyapunov functional.In spatially non-homogeneous situations, the main dif-ference is that an additional integral over x is present,both in H B and, consequently, dH B /dt . There is no rea-son to expect this integral over space to make dH/dt | lin vanish, since one still has that dH B dt (cid:12)(cid:12)(cid:12)(cid:12) lin = ddt (cid:104) ln f st (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) lin , (41)and, in general, ln f st is not a sum of constants of motion.In fact, again the sign of dH B /dt | lin is reversed when∆ f → − ∆ f , similarly to the homogeneous case. Wehave numerically checked this prediction for the shearedsystem, with the resulting evolution of H B being com-pletely similar to that for the uniformly heated case in Fig. 1, and thus is not shown here. V. H -THEOREM FOR SOME SPECIFIC NESS Here we prove that the functional H [ f ] is monoton-ically decreasing for all times in some specific physicalsituations. Our proof applies both to the sheared anduniformly heated systems described in Sec. II A. To beas general as possible, we consider a system that is bothheated and sheared: a (cid:54) = 0 and ξ (cid:54) = 0. In this situa-tion, the boundary conditions for the PDF are given byEq. (20), which lead to Eqs. (19a) and (19b) for the av-erages u ( x, t ) and T ( x, t ).The steady solution of the hydrodynamic equations is u st ( x ) = a (cid:18) x − (cid:19) , T st = 2 a + ξν . (42)On the one hand, the average velocity has a linear pro-file, similarly to the situation in the USF state. On theother hand, the temperature remains homogeneous butits steady value has two contributions, one coming fromthe shear and the other from the stochastic thermostat.The viscous heating 2( ∂ x u ) and uniform heating ξ termscancel the cooling term − νT for all x . The stationary so-lution of the kinetic equation is readily found: f st ( x, v ) = (2 πT st ) − / exp (cid:34) − ( v − u st ( x )) T st (cid:35) , (43)that is, the Gaussian distribution corresponding to thehydrodynamic fields in Eq. (42). Of course, the USF stateand NESS of the uniformly heated system in Sec. II Acan be easily recovered as particular cases of Eq. (43):for ( a (cid:54) = 0 , ξ = 0) and ( a = 0 , ξ > H -theorem, that is, the existence of a nonequi-librium entropy ensuring the monotonic approach of theone-particle PDF to the steady state. Our starting pointis the following expansion of the one-particle PDF in Her-mite polynomials, f ( x, v, t ) = 1 (cid:112) πT ( x, t ) exp (cid:20) − [ v − u ( x, t )] T ( x, t ) (cid:21) × (cid:34) ∞ (cid:88) n =3 γ n ( x, t ) H n (cid:32) v − u ( x, t ) (cid:112) T ( x, t ) (cid:33)(cid:35) , (44)which is known as the Gram-Charlier series [43–46].Therein, u ( x, t ) and T ( x, t ) are the (exact) average veloc-ity and temperature stemming from the hydrodynamicequations for the considered distribution. The above ex-pansion is suggested by the Gaussian shape of the sta-tionary PDF in Eq. (43). Now we define c = v − u ( x, t ) (cid:112) T ( x, t ) , ˜ f ( x, c, t ) = (cid:112) T ( x, t ) f ( x, v, t ) , (45)From the orthogonality relation of the Hermite polyno-mials [47], it is readily obtained that γ n ( x, t ) = 1 n ! (cid:90) dc H n ( c ) ˜ f ( x, c, t ) . (46)Also, we could write γ n as a combination of moments ofthe distribution.Some comments on the Gram-Charlier expansion arepertinent. First, note that n ≥ γ = γ = 0because the zero-th order Gaussian contribution exactlygives the first two moments u ( x, t ) and (cid:104) v (cid:105) ( x, t ) = u ( x, t ) + T ( x, t ). Second, if f ( x, v, t ) were symmetricwith respect to v = u , that is, (cid:104) ( v − u ) n +1 (cid:105) = 0 for all n ∈ N , only even values of n would be present in thesum and one would end up with the usual expansion inSonine-Laguerre polynomials of kinetic theory. Finally, it is worth stressing that the series (44) converges forfunctions such that the tails of ˜ f ( x, c, t ) approach zerofaster than e − c / for c → ±∞ [46, 48, 49].After a lengthy but straightforward calculation, whichis summarised in Appendix C, it is shown that dHdt = A ( t ) + B ( t ) , with both A ( t ) , B ( t ) ≤ . (47)The expressions for A ( t ) and B ( t ) are A ( t ) = − (cid:90) dx T (cid:34)(cid:18) u (cid:48) T − u (cid:48) st T st (cid:19) + ξ (cid:18) T − T st (cid:19) (cid:35) (48)and B ( t ) = − √ π (cid:90) dx dc e − c / (cid:80) ∞ n =3 γ n H n ( c ) (cid:40) T (cid:48) T H ( c ) + ∞ (cid:88) n =3 γ (cid:48) n H n ( c ) − ∞ (cid:88) n =3 γ n u (cid:48) √ T nH n − ( c )+ ∞ (cid:88) n =3 γ n T (cid:48) T [ H n +2 ( c ) + nH n ( c )] (cid:41) − ξ √ π (cid:90) dx dc e − c / (cid:80) ∞ n =3 γ n H n ( c ) 1 T (cid:34) ∞ (cid:88) n =3 γ n nH n − ( c ) (cid:35) . (49)We recall that the prime denotes spatial derivative.Therefore, dH/dt ≤ H -theorem holds for the sheared and heatedsystem. Rigorously, our proof holds for those PDFs suchthat the above Hermite expansion converges. Note thatthe proof remains valid for the approach to any NESS,whose PDF is a Gaussian with a homogeneous tempera-ture, independently of the corresponding boundary con-ditions. In Sec. III, we have already demonstrated that dH/dt only vanishes for f ( x, v, ∞ ) = f st ( x, v ), but thesame result can be rederived here in a different way. Byimposing that both A ( t ) and B ( t ) vanish in the long timelimit and making use of the hydrodynamic equations forthe averages, it can be shown that u ( x, ∞ ) = u st ( x ), T ( x, ∞ ) = T st and γ n ( x, ∞ ) = 0, ∀ n ≥ A. Numerical results for the USF state
Here we consider the sheared system, and we numeri-cally check our theoretical predictions. Throughout thissection, we use the values of the parameters ν = 20, a = 5and ξ = 0 (there is no stochastic forcing).Firstly, in Figure 2, we show the evolution of the distri-bution and the H -functional from a Gaussian initial con-dition with the steady velocity profile u ( x,
0) = u st ( x )but a higher temperature, T ( t = 0) = 7 T st . In panel(a), we depict the velocity distribution at x = 1 / κ = (cid:104) [ v − u ( x )] (cid:105) / (cid:104) [ v − u ( x )] (cid:105) − x = 1 /
4, consistently with theGaussian shape. In panel (b), it is clearly observed thatthe H -functional is monotonically decreasing with time.Secondly, we study the relaxation to the USF statefrom another initial preparation, for which the veloc-ity profile u ( x,
0) is different from the stationary but T ( x,
0) = T st . The numerical results are shown in Fig. 3,and for the sake of simplicity we use again an initialGaussian distribution. Specifically, we use u ( x,
0) = u st ( x )+4 . πx ). Here, the departure from the Gaus-sian shape is evident, and thus we have not plotted thekurtosis. Consistently with our theoretical prediction, weget again a monotonous relaxation of H towards its nullstationary value.Finally, we consider situations for which the above pre-sented proof is not rigorously applicable. As stated be-fore, the Gram-Charlier series does not converge whenthe tails of the distribution decay to zero slower than thesquare root of the Gaussian. Nevertheless, when all thecoefficients γ n defined in Eq. (46) exist and are finite, westill expect the H -theorem to hold. We illustrate thissituation with an initial exponential distribution; specif- FIG. 2. (Color online) Relaxation towards the USF state.The initial condition is Gaussian, centered in u st ( x ) and withvariance T = 7 T st . (a) Velocity distribution function at x = 1 / H functional, in clear agreement with the proven H -theorem. Data is averaged over 6000 trajectories in a systemwith N = 660 sites, ν = 20, and a = 5. Solid lines cor-respond to the (theoretical) Gaussian distributions for theplotted times, except for the last time in which it representsthe theoretical steady distribution. ically, we consider f ( x, v,
0) = 1 (cid:112) T ( t = 0) exp (cid:34) √ | v − u ( x, t = 0) | (cid:112) T ( t = 0) (cid:35) , (50)with u ( x, t = 0) = u st ( x ) + 4 . πx ) and T ( t = 0) =0 . T st . Consistently with our expectation, we can see inFig. 4 that indeed the H -functional also monotonicallydecreases. FIG. 3. (Color online) The same plots as in Fig. 2, but start-ing from a different initial condition. Now, the initial PDFis a Gaussian centered in u ( x,
0) = u st ( x ) + 4 . πx ) andwith variance T ( t = 0) = T st . In (a), solid lines correspond tothe theoretical PDFs for the initial time and the steady state.In (b), H decreases again monotonically towards its steadyvalue, consistently with our theoretical prediction. B. Numerical results in the uniformly heatedsystem
To conclude, we put forward the results of simulationsfor the uniformly heated system. Specifically, our sim-ulations have been done for ν = 20, a = 0 (no shear)and ξ = 50. In order not to overload the reader with toomany examples, we only present the more complex casein Fig. 5: the relaxation towards the steady state froman initial exponential distribution, as given by Eq. (50).In particular, we consider that u ( x, t = 0) = 4 . πx )and T ( t = 0) = 0 . T st . Note that the perturbation fromthe steady values is the same as in Fig. 4 for the shearedcase. Again, we observe the monotonic relaxation of H towards the stationary value, consistently with our theo-retical result, even for a initial distribution for which theGram-Charlier series does not converge.0 FIG. 4. (Color online) The same plots as in Fig. 2, butstarting from an initial PDF with a divergent Gram-Charlierseries. Concretely, the plots correspond to an exponential ini-tial distribution centered in u ( x, t = 0) = u st ( x )+4 . πx )and with T ( t = 0) = 0 . T st . VI. CONCLUSIONS
Within a simplified model for a granular-like velocityfield, we have analytically shown that the nonequilib-rium steady state that the system eventually reaches inthe long time limit is globally stable. This has been donefor quite a general situation, in which energy may be in-jected into the system both through the boundaries andby a heating mechanism that acts in the bulk. The proofis valid both for spatially homogeneous situations (suchas the uniformly heated system) and inhomogeneous sit-uations (such as the USF or Couette states).The proof of global stability is based on showing thatthe H -functional H = (cid:82) dx dv f ln( f /f st ) (the Kullback-Leibler divergence between the time evolving one-particlePDF f ( x, v, t ) and its value in the stationary state f st ( x, v ) [50]) is non-increasing in the infinite time limit.Thus, we do not need H to be a “good” Lyapunov func- FIG. 5. (Color online) Numerical results for the uniformlyheated system. The plots are analogous to those in Fig. 2,(a) time evolution of the PDF and (b) time evolution of the H -functional. The system is initially prepared with an expo-nential PDF centered in u ( x, t = 0) = 4 . πx ) and with T ( t = 0) = 0 . T st . Data is averaged over 3000 trajectories ina system with N = 330 sites. tional for all times in order to prove global stability.In conclusion, global stability and the validity of an H -theorem do not seem to be unavoidably tied.Moreover, we have analytically shown that the Boltz-mann functional H B = (cid:82) dx dv f ln f cannot be, ingeneral, a Lyapunov functional for systems with non-conservative interactions. Close to the steady state, wehave proven that dH B /dt contains non-vanishing termsthat are linear in the deviations ∆ f = f − f st . Therefore,a reversal of the sign of ∆ f entails a reversal of the signof dH B /dt . This general analytical proof on the inade-quacy of H B as a Lyapunov functional is in agreementwith previous results in some specific cases [7, 11].We have also succeeded in demonstrating that the H -functional is non-increasing and thus a “good” Lyapunovfunctional for some specific driving mechanisms. Our1proof is not restricted to spatially homogeneous situa-tions and is applicable to two relevant physical cases:the approach to (i) the USF state and (ii) the NESScorresponding to the uniformly heated case. The proofinvolves a suitable expansion of the one-particle PDFin Hermite polynomials, which is a generalisation of thewell-known Sonine-Laguerre expansion in kinetic theory.Although the proof is only rigorous for PDFs having aconvergent series expansion, we expect it to remain validfor more general PDFs. In fact, we have numericallyvalidated this expectation in some specific situations.The analytical results presented here are thus in agree-ment with the numerical evidence in Refs. [7, 8] and ad-vance the understanding of this field in a twofold way.First, an analytical proof, which was lacking, is providedfor a simplified model. Second, spatially inhomogeneoussituations are considered, both in the time evolution andin its steady state.Some limitations of our results have to be under-lined, though. First, the simplifications introduced inthe model make it impossible to address the problemof the stability of the homogeneous cooling state in theundriven system at the level of the kinetic equation, asalready discussed in Ref. [14]. Second, for the drivingmechanisms for which we can analytically prove that H is a “good” Lyapunov functional, the steady distributionis exactly Gaussian. Nevertheless, we think that thisis not a fundamental point and expect that the kind ofexpansion-based proof presented here may be extendedto other situations. A particularly appealing case is theapproach to the Couette NESS, for which the stationaryPDF is non-Gaussian in the model [14].On a different note, our kinetic equation (14) showssome resemblance to evolution equations for the one-particle PDF found in other physical contexts, such as theVlasov equation in plasma physics or astrophysics [51–53]or the non-linear (in the distribution function) Fokker-Planck equation for systems of infinitely many couplednon-linear oscillators exhibiting phase transitions [54]. Itis a drift-like term depending on a certain average of thePDF that all these different problems share. Both forthe Vlasov and the non-linear Fokker-Planck equations,the existence of a Lyapunov functional has been provedby considering a variant of the functional H [ f ] definedin Eq. (4) [55–57]. Thus, an interesting prospect is toinvestigate if this kind of approach may be extended toour class of models with non-conservative interactions.Our work also opens the door to applying the ideas de-veloped in this paper to more complex models, closer toreal non-conservative systems, like granular gases. Thepioneering numerical work in Refs. [7, 8] strongly suggeststhat the H -functional is a “good” Lyapunov functionalfor granular fluids. It seems worth trying to analyticallyprove that this is indeed the case for the inelastic Boltz-mann equation, at least for some specific situations. Ifnothing else, one would like to be able to show that thelong time solutions are globally stable by showing that H is asymptotically non-increasing, similarly to what has been done here. ACKNOWLEDGMENTS
We acknowledge the support of the Spanish Ministeriode Econom´ıa y Competitividad through Grant FIS2014-53808-P. Carlos A. Plata also acknowledges the supportfrom the FPU Fellowship Programme of the Spanish Min-isterio de Educaci´on, Cultura y Deporte through GrantFPU14/00241.
Appendix A: Derivation of the expression for dH/dt in a general driven state
Let us consider the three contributions to dH/dt inEq. (29). We start with the diffusive one, dHdt (cid:12)(cid:12)(cid:12)(cid:12) diff = (cid:90) dx dv L diff f ln (cid:18) ff st (cid:19) − (cid:90) dx dv ff st L diff f st , (A1)where L diff f = ∂ x f and we have used that (cid:82) dx dv ∂ x f st vanishes identically. Integrating by parts the first termon the rhs of Eq. (A1), the result is (cid:90) dv ∂ x f ln (cid:18) ff st (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) − (cid:90) dx dv f ∂ x ln f ( ∂ x ln f − ∂ x ln f st ) . (A2)Also integrating by parts the second term, one obtains − (cid:90) dv ff st ∂ x f st (cid:12)(cid:12)(cid:12)(cid:12) + (cid:90) dx dv f ∂ x ln f st ( ∂ x ln f − ∂ x ln f st ) . (A3)We assume that the boundary terms are equal to zero,that is, (cid:90) dv (cid:20) ∂ x f ln (cid:18) ff st (cid:19) − ff st ∂ x f st (cid:21) = 0 . (A4)This is obviously true for Lees-Edwards and periodicboundary conditions [58]. Summing the two contribu-tions to the diffusive term above, we have dHdt (cid:12)(cid:12)(cid:12)(cid:12) diff = − (cid:90) dx dv f ( ∂ x ln f − ∂ x ln f st ) , (A5)which is Eq. (30) of the main text.The noise term is treated along the same lines as above,but integrating by parts in v instead of x , since L noise f = ξ ∂ v f . There in, the boundary terms vanish if f and f st tend to zero fast enough for v → ±∞ , and dHdt (cid:12)(cid:12)(cid:12)(cid:12) noise = − ξ (cid:90) dx dv f ( ∂ v ln f − ∂ v ln f st ) , (A6)which is Eq. (32).2Now we focus on the inelastic contribution, dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel = (cid:90) dx dv L inel f ln (cid:18) ff st (cid:19) − (cid:90) dx dv ff st L inel f st , (A7)in which L inel f = ν ∂ v [( v − u ) f ]. Then, dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel = ν (cid:90) dx dv ∂ v [( v − u ) f ] ln (cid:18) ff st (cid:19) − ν (cid:90) dx dv ∂ v [( v − u st ) f st ] f − f st f st . (A8)Again, integrating by parts in v (here we do not writethe boundary terms at v → ±∞ ), the first term on therhs of Eq. (A8) is − ν (cid:90) dx dv ( v − u ) f ( ∂ v ln f − ∂ v ln f st ) , (A9)whereas the second term gives ν (cid:90) dx dv ( v − u st ) f ( ∂ v ln f − ∂ v ln f st ) . (A10)Summing up these two contributions, and taking intoaccount that both u and u st do not depend on v , dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel = ν (cid:90) dx ( u − u st ) (cid:90) dv f ( ∂ v ln f − ∂ v ln f st )(A11)Since (cid:82) dv f ∂ v ln f ≡ (cid:82) dv ∂ v f = 0, this leads to Eq. (31). Appendix B: Simulation strategy
In the simulations, so as to generate a trajectory ofthe stochastic process, we proceed as follows. (i) A pair( l, l + 1) is chosen at random and undergoes the inelasticcollision described by Eq. (5), (ii) all the particles are sub-mitted to the stochastic thermostat according to (6) and(7), and (iii) time is incremented by δτ = − ( N ω ) − ln x ,with x being a homogeneously distributed random num-ber in (0 ,
1) [59–62]. This cycle (random choice of a pairand noise interaction followed by a time increment) isrepeated until time exceeds some maximum time t max . Regarding the measurements of f ( x, v, t ), we sampleboth position and velocity spaces by defining N x bins ofwidth ∆ x and N v bins of with ∆ v . Of course, the product N x ∆ x = 1, the whole lattice, whereas N v ∆ v gives therange of velocities bounded by the cutoffs v min and v max .In our simulations, we control that the contribution tothe PDFs coming from velocities outside the consideredinterval [ v min , v max ] is negligible. With such a binning,we build up an histogram and therefrom the distributionfunction f ( x, v, t ), which is represented by a N x × N v matrix for each time t . Both H and H B are computedby numerically replacing the integral over x and v withsums over the prescribed bins. Appendix C: Derivation of Eq. (47)
In order to derive Eq. (47), we have to substitute theGaussian stationary solution (43) and the Gram-Charlierseries (44) into the three contributions to dH/dt , givenby Eqs. (30), (31) and (32).For the inelastic term, it is readily obtained that dHdt (cid:12)(cid:12)(cid:12)(cid:12) inel = ν T st (cid:90) dx ( u − u st ) . (C1)For the diffusive and noise terms, the key ideas are achanging the integration over velocities from v to c =( v − u ) / √ T and the use of the recursion relations andthe orthogonality property of the Hermite polynomials[47]. Working along these guidelines, we arrive at dHdt (cid:12)(cid:12)(cid:12)(cid:12) diff = − (cid:90) dx T (cid:18) u (cid:48) T − u (cid:48) st T st (cid:19) − u (cid:48) T (cid:90) dx ( u − u st ) + B ( t ) , (C2) dHdt (cid:12)(cid:12)(cid:12)(cid:12) noise = − ξ (cid:90) dx T (cid:18) T − T st (cid:19) − ξ T (cid:90) dx ( u − u st ) + B ( t ) , (C3)where B and B are, respectively, the first and the sec-ond term in Eq. (49). The sum of the factors multiplying (cid:82) dx ( u − u st ) vanishes by taking into account the equa-tion for the (spatially homogeneous) stationary temper-ature. Therefore, the sum of the remaining terms leadsright to Eq. (47). [1] A. M. Lyapunov, Int. J. Cont. , 531 (1992).[2] L. Boltzmann, Lectures on Gas Theory (Dover, 1995).[3] J. L. Lebowitz, Phys. Today , 32 (1993).[4] J. L. Lebowitz, Physica A , 1 (1993).[5] S. Chapman and T. G. Cowling, The Mathematical The-ory of Non-Uniform Gases , 3rd ed. (Cambridge Univer-sity Press, Cambridge, UK, 1970).[6] N. G. Van Kampen,
Stochastic processes in physics and chemistry , Vol. 1 (Access Online via Elsevier, 1992).[7] U. M. B. Marconi, A. Puglisi, and A. Vulpiani, J. Stat.Mech. , P08003 (2013).[8] M. I. Garc´ıa de Soria, P. Maynar, S. Mischler, C. Mouhot,T. Rey, and E. Trizac, J. Stat. Mech. , P11009 (2015).[9] T. P¨oschel and S. Luding, eds.,
Granular Gases , LectureNotes in Physics, Vol. 564 (Springer, Berlin, 2001).[10] C. Villani, J. Stat. Phys. , 781–822 (2006). [11] I. Bena, F. Coppex, M. Droz, P. Visco, E. Trizac, andF. van Wijland, Physica A , 179 (2006).[12] A. Baldassarri, U. Marini Bettolo Marconi, andA. Puglisi, EPL , 14 (2002).[13] A. Lasanta, A. Manacorda, A. Prados, and A. Puglisi,New J. Phys. , 083039 (2015).[14] A. Manacorda, C. A. Plata, A. Lasanta, A. Puglisi, andA. Prados, J. Stat. Phys. , 810 (2016).[15] C. A. Plata, A. Manacorda, A. Lasanta, A. Puglisi, andA. Prados, J. Stat. Mech. , 093203 (2016).[16] V. V. Prasad, S. Sabhapandit, A. Dhar, and O. Narayan,Phys. Rev. E , 022115 (2017).[17] M. Kac, in Proc. 3rd Berkeley Symp. Mathematics Statis-tics Probability , Vol. 3 (Berkeley, CA: University of Cal-ifornia Press, 1956) p. 171.[18] J. J. Brey and A. Prados, Phys. Rev. E , 1541 (1993).[19] J. J. Brey and A. Prados, Phys. Rev. B , 984 (1994).[20] J. J. Brey, A. Prados, and M. J. Ruiz-Montero, J. Non-Cryst. Solids , 371 (1994).[21] M. O. Vlad and J. Ross, J. Phys. Chem. B , 8756(1997).[22] M. O. Vlad, F. Moran, and J. Ross, J. Phys. Chem. B , 4598 (1998).[23] A. Prados, J. J. Brey, and B. S´anchez-Rey, Physica A , 277 (2000).[24] B. Earnshaw and J. Keener, SIAM J. Appl. Dyn. Syst. , 220 (2010).[25] T. P. C. Van Noije and M. H. Ernst, Granul. Matter ,57 (1998).[26] J. M. Montanero and A. Santos, Granul. Matter , 53(2000).[27] M. I. Garc´ıa de Soria, P. Maynar, and E. Trizac, Mol.Phys. , 383 (2009).[28] M. I. Garc´ıa de Soria, P. Maynar, and E. Trizac, Phys.Rev. E , 051301 (2012).[29] I. Pagonabarraga, E. Trizac, T. P. C. van Noije, andM. H. Ernst, Phys. Rev. E , 011303 (2002).[30] P. Maynar, M. I. Garc´ıa de Soria, and E. Trizac, Eur.Phys. J. Spec. Top. , 123 (2009).[31] A. Prados and E. Trizac, Phys. Rev. Lett. , 198001(2014).[32] E. Trizac and A. Prados, Phys. Rev. E , 012204 (2014).[33] V. V. Prasad, S. Sabhapandit, and A. Dhar, EPL ,54003 (2013).[34] E. Ben-Naim and P. L. Krapivsky, in Granular Gas Dy-namics , Lecture Notes in Physics, Vol. 624, edited byT. P¨oschel and N. Brilliantov (Springer, Berlin, 2003)pp. 65–94.[35] This kind of Fokker-Planck kernel is similar to the oneappearing in the Fokker-Planck equation for biomoleculesin which their total length L = (cid:80) Nl =1 η l , where η l is thelength of each of the modules of the molecule, is con-trolled and kept constant. This is quite logical, sincein both situations there is an additive conservation law: here it is the total momentum (cid:80) Nl =1 v l that is conserved,whereas in biomolecules the constant of motion is thetotal length [63].[36] Note that the density of the model is fixed, there is nomass transport in the system.[37] C. Kipnis, C. Marchioro, and E. Presutti, J. Stat. Phys. , 65 (1982).[38] A. Prados, A. Lasanta, and P. I. Hurtado, Phys. Rev.Lett. , 140601 (2011).[39] A. Prados, A. Lasanta, and P. I. Hurtado, Phys. Rev. E , 031134 (2012).[40] P. I. Hurtado, A. Lasanta, and A. Prados, Phys. Rev. E , 022110 (2013).[41] A. Lasanta, P. I. Hurtado, and A. Prados, Eur. Phys. J.E , 35 (2016).[42] A. W. Lees and S. F. Edwards, J. Phys. C , 1921 (1972).[43] P. L. Chebyshev, Acta Math. , 305 (1890).[44] C. V. L. Charlier, Ark. Math. Astr. och Phys. , 1 (1905-06).[45] F. Y. Edgeworth, Cambridge Philos. Soc. , 36 (1905).[46] D. L. Wallace, Ann. Math. Stat. , 635 (1958).[47] M. Abramowitz and I. A. Stegun, eds., Handbook ofMathematical Functions (Dover, New York, 1972).[48] H. Cram´er, in
Proceedings of the Sixth ScandinavianCongress of Mathematicians, Copenhagen (1925) pp.399–425.[49] G. Szeg¨o, Amer. Math. Soc. (1939).[50] T. M. Cover and J. M. Thomas, Elements of InformationTheory (New York: Wiley, 2006).[51] P.-H. Chavanis, J. Sommeria, and R. Robert, Astrophys.J. , 385 (1996).[52] P.-H. Chavanis, Ann. N. Y. Acad. Sci. , 120 (1998).[53] Y. Y. Yamaguchi, J. Barr´e, F. Bouchet, T. Dauxois, andS. Ruffo, Physica A , 36 (2004).[54] R. C. Desai and R. Zwanzig, J. Stat. Phys. , 1 (1978).[55] M. Shiino, Physics Letters A , 302 (1985).[56] M. Shiino, Phys. Rev. A , 2393 (1987).[57] L. L. Bonilla, J. A. Carrillo, and J. Soler, Phys. Lett. A , 55 (1996).[58] For the Couette state, in which the PDF at the bound-aries is Gaussian with zero average velocity and a giventemperature T B for all times, the first term is identicallyzero and the second vanishes because (cid:82) dv f st ( x, v ) = 1for all x .[59] A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comput.Phys. , 10 (1975).[60] D. T. Gillespie, J. Comput. Phys. , 403 (1976).[61] A. Prados, J. J. Brey, and B. S´anchez-Rey, J. Stat. Phys. , 709 (1997).[62] S. A. Serebrinsky, Phys. Rev. E , 037701 (2011).[63] L. L. Bonilla, A. Carpio, and A. Prados, Phys. Rev. E91