Hydrodynamics of granular particles on a line

PHYSICAL REVIEW E 97, 062905 (2018) Hydrodynamics of granular particles on a line Andrea Baldassarri and Andrea Puglisi Istituto dei Sistemi Complessi - CNR and Dipartimento di Fisica, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Rome, Italy Antonio Prados* Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain (Received 30 January 2018; published 18 June 2018) We investigate a lattice model representing a granular gas in a thin channel. We deduce the hydrodynamic description for the model from the microscopic dynamics in the large-system limit, including the lowest finite-size corrections. The main prediction from hydrodynamics, when finite-size corrections are neglected, is the existence of a steady “uniform longitudinal flow” (ULF), with the granular temperature and the velocity gradient both uniform and directly related. Extensive numerical simulations of the system show that such a state can be observed in the bulk of a finite-size system by attaching two thermostats with the same temperature at its boundaries. The relation between the ULF state and the shocks appearing in the late stage of a cooling gas of inelastic hard rods is discussed. DOI: 10.1103/PhysRevE.97.062905 I. INTRODUCTION Statistical physics offers a systematic computational scheme for averages and fluctuations, suitable for systems of particles at thermal equilibrium [1]. A similar scheme for systems far from equilibrium is lacking and represents an open challenge [2]. An example of a nonequilibrium statistical system pervading many human activities is given by granular fluids, an assembly of inelastic hard particles agitated by some external driving mechanism [3–5]. The fluidized state of granular matter is an excellent testing ground for kinetic theory [6–8]. One of the most widely used models in granular kinetic theory is a system composed of inelastic smooth hard particles [9–11]. In this context, the inelastic Boltzmann (or Enskog) equation has been shown to be a powerful tool. However, the structure of the collision term for inelastic hard particles makes the corresponding kinetic equation a tough mathematical problem, for which exact solutions and rigorous results are difficult to obtain. In view of the above issue, models that simplify the collision term have been proposed to make their analytical investigation easier, while preserving similar physics. This is the spirit of the Maxwell-like collision model, either elastic [12]or inelastic [13], in which the collision rate is assumed to be independent of the relative velocity. Some experts claim that “What harmonic oscillators are for quantum mechanics, and dumbbells for polymer physics, is what elastic and inelastic Maxwell models are for kinetic theory” [14]. Indeed, many rigorous results in the inelastic case have been derived, both in the freely cooling [13–21] and the uniformly heated [13,22,23] cases. As pointed out in a review article by Villani [24], the most relevant questions in inelastic Maxwell models have been solved in homogeneous situations and, rather than going for *[email protected] refinements, the current priority is to deal with inhomogeneous states. The nonconservation of energy has several, physically relevant, implications. First, the distribution function is in general non-Gaussian [25,26], even for the stationary states reached when some external mechanism injects energy into the system. Second, and most importantly, kinetic theory establishes a link between the microscopic and the macroscopic, hydrodynamic, descriptions, by making it possible to derive the latter from the former [6,10,11]. Along this route, the nonconservation of energy introduces another timescale that makes more delicate the critical requirement on separation of scales needed to accomplish it [27–29]. However, there are many situations, typically in dimensions larger than one and for dilute quasielastic systems, in which hydrodynamic equations fairly reproduce real granular experiments qualitatively [30]or even quantitatively [31,32]. Here we aim at building and verifying the hydrodynamic equations for an idealized one-dimensional (1d) lattice model [15,33] with a Maxwell-like collision rule. The model simplifies the dynamics of a granular fluid in a 1d channel with negligible density fluctuations, a condition that reasonably holds in not too dilute systems. A strong connection between this idealized model and an inelastic 1d gas of hard rods [34] has been established in previous studies [15]. Specifically, there appears a nontrivial correspondence between the velocity profiles of the lattice model and the shocklike structures of the granular gas of hard rods in the asymptotic cooling regime. The 1d lattice model has been further studied in [35,36] and also generalized to two dimensions [37]. In spite of the above described strong connection between 1d granular gases and the 1d Maxwell model, the hydrodynamic equations of the latter have been neither written nor analyzed. Then, in the present study, first we derive the hydrodynamic equations from the microscopic dynamics. Afterwards, we focus on a solution of them: the steady uniform 2470-0045/2018/97(6)/062905(11) 062905-1 ©2018 American Physical Society BALDASSARRI, PUGLISI, AND PRADOS PHYSICAL REVIEW E 97, 062905 (2018) longitudinal flow (ULF). We argue that such a state resembles the shock profiles observed in the late cooling regime of many 1d granular systems. Also, we incorporate finite-size corrections and compare them with simulations of the microscopic dynamics. Section II is devoted to a brief overview of previous studies on granular models in 1dand their comparison with hydrodynamic theories. Our model is introduced in Sec. III;inthesame section its hydrodynamic equations are derived and the ULF solution is discussed. Section IV presents the numerical results and their comparison with theory. Finite-size corrections and boundary layers are discussed in Sec. V. Conclusions are drawn in Sec. VI. Appendices illustrate the technical aspects of the derivation of the hydrodynamics (Appendix A) and of the boundary layer calculations (Appendix B), and give also some numerical results for the local velocity distributions and spatial correlations (Appendix C). II. 1D GRANULAR HYDRODYNAMICS: A BRIEF OVERVIEW Despite the progress made during the last decades in the realm of granular kinetic theory and hydrodynamics, the validity of a hydrodynamic description in 1d systems is still under debate. It should be stressed that 1d systems essentially differ from higher-dimensional setups in which the gradients have a well-defined direction: particles cannot hop over their nearest neighbors, which hinders the necessary “mixing” to obtain a continuum description. We briefly review some of the main results on 1d granular hydrodynamics below. The 1d version of the granular Navier- Stokes equations reads Dtρ=−ρ∂xu, (1a) ρDtu=−∂xP, (1b) ρDtT=−2P∂ xu−2∂xq−ρζT, (1c) where the density ρ, the mean velocity u, the granular temperature T, the pressure P, the heat flux q, and the cooling rate ζ depend on (x,t). See, for instance, [38] for their microscopic definitions. As usual, we have used the notation Dt≡∂t+u∂x for the material derivative. One of the first tests of 1d granular hydrodynamics was carried out in [39]. Therein, a gas of inelastic hard rods was confined by one or two thermostatted walls. Numerical simulations revealed an incompatibility between Eqs. (1) and the observed steady hydrodynamic profiles. Later, the same system has been demonstrated to lack a pure thermodynamic limit, i.e., infinite size at constant restitution coefficient [40,41]. Subsequent studies of the 1d granular gas focused on the cooling state, that is, with periodic boundary conditions and no external driving. The 1d gas of hard rods develops strong inhomogeneities, which include clusters and shocks, for long enough times [34]. Extensive numerical simulations [42]have revealed that this 1d gas becomes then indistinguishable from a perfectly “sticky” gas [34,43,44], regardless of its actual inelasticity, a conjecture analyzed also in higher dimensions [45,46] and for wet particles [47]. Interestingly, the (inviscid) Burgers equation Dtu=μ∂2 xu[48,49], with μ→0, describes the sticky gas [50–52]. Notwithstanding, to the best of our knowledge, an analytical justification of the Burgers equation in 1d granular gases has not been obtained yet. A first step in this direction would be a neat derivation of 1d granular hydrodynamics. 1d continuum equations can appear also in higher dimensions, when only one direction develops a gradient, as in [53] for a steady case under gravity or in cooling systems with high aspect ratio without gravity [54–56]. In the latter case, granular hydrodynamics predicts that the so-called flow by inertia Dtu=0 sets in during the highly inhomogeneous stage of cooling. Therefore, there also appear shocks in u(x,t), which imply a singularity in ρ(x,t). In addition, molecular dynamics simulations of this channel system have made clear that such singularities are approached following the hydrodynamic predictions, that is, the flow by inertia scenario [55]. Eventually, when close-packing is almost reached, the flow by inertia picture breaks down [56,57]. III. THE MODEL AND ITS HYDRODYNAMICS Here, we study the 1d lattice model first introduced in Refs. [15,33] to mimic the evolution of a 1d granular gas, but neglecting density fluctuations. The model is a 1d space-time discretized cellular automaton, where we have one unit-mass particle at each site l,l=1,...,N. Particles are characterized by their (scalar) “velocity” vl, which evolves according to a collisional dynamics. Note that particles do not move but remain on the same site; in this way density fluctuations are neglected. At each discrete time step, a random pair of adjacent sites is chosen and their corresponding velocities, say vand ˜v,are updated according to v=v+1+α 2(˜v−v),˜v=˜v−1+α 2(˜v−v),(2) if vand ˜vsatisfy the so-called kinematic constraint, that is, if the corresponding particles are “approaching” and not ‘moving apart’, see Fig. 1.InEq.(2), vand ˜vare the postcollisional velocities, and α∈(0,1) is the restitution coefficient. Momentum is conserved in collisions, v+˜v=v+˜v, but energy is not: K−K=1 2(α2−1)(˜v−v)2<0, with K=v2+˜v2. The choice of the (candidate) colliding pair is uniformly distributed among all nearest-neighbor pairs, as for “pseudo- Maxwell” molecules [12,13,15,24]. The dropping of the kinematic constraint leads to a completely different class of models without kinematic constraint [58–66]. This choice is justified by a different physical picture: the on-site velocities on the lattice are thought as transverse, not longitudinal, velocities. This interpretation is confirmed by FIG. 1. Allowed (green, left) and forbidden (red, right) collision, according to the kinematic constraint. Arrows show (precollisional) velocities of a pair of (candidate) colliding particles. Green pairs actually collide; red pairs do not collide. 062905-2 HYDRODYNAMICS OF GRANULAR PARTICLES ON A LINE PHYSICAL REVIEW E 97, 062905 (2018) the equivalence of the corresponding hydrodynamic equations with those for the shear mode of inelastic gases in higher dimension [60,61]. Also, variants of the model without kinematic constraint have been employed to analyze the dynamical evolution of social and economic systems [20,67–72]. A physical comment on our adoption of nonmoving particles on a lattice is in order. This may be understood as if we were adopting a Lagrangian coordinate in the actual 1d granular gas, by characterizing the particle positions by their index iinstead of their real positions xi. The corresponding Lagrangian density profile, defined by the mass per unit of space measured by i, would be uniform. Within this physical picture, the partial time derivative ∂tin our model is expected to play the role of the material derivative Dtin the 1d granular gas. A. Microscopic evolution of velocity and energy At discrete time p, the pair {yp,yp+1}that may undergo a collision is chosen at random and the velocity on site lat time p+1 is given by vl,p+1−vl,p =−jl,p +jl−1,p,(3) where we have defined the microscopic momentum flux from site lto site l+1 at time p: jl,p =δyp,l (vl,p −vl+1,p)1+α 2(vl,p −vl+1,p).(4) Therein, δij is Kronecker’s delta, which identifies the colliding pair, and (x) is Heaviside’s step function, which imposes the kinematic constraint. The random integer ypis uniformly distributed in [1,L], where Lis basically equal to Nbut depends on the boundary conditions. For a thermostatted system L=N+1. Obviously, jl,p only differs from zero when the pair (l,l +1) actually collides. Note that no external volume forces (like gravity, for instance) are being considered, but they could be incorporated by adding a term fl,p to the right-hand side of Eq. (3). The evolution of the kinetic energy is obtained by squaring Eq. (3), which after some algebra yields v2 l,p+1−v2 l,p =α2−1 4[δyp,l(vl,p −vl+1,p) ×(vl,p −vl+1,p)2+δyp,l−1(vl−1,p −vl,p) ×(vl−1,p −vl,p)2]−Jl,p +Jl−1,p.(5) The microscopic energy current is Jl,p =(vl,p +vl+1,p)jl,p. The sink terms in squared brackets on the right-hand side of Eq. (5) stem from the inelasticity of collisions and lead to a monotonic decrease of the total energy. The last two terms are a discrete spatial derivative and correspond to the energy flux between neighboring sites, which is already present in the conservative case α=1. B. Equations for the velocity and temperature fields in the continuum limit The average (over realizations) velocity field ul,p =vl,p evolves according to ul,p+1−ul,p =−jl,p −jl−1,p.(6) To compute the last average, we introduce a local equilibrium assumption P2(vl,vl+1)≃1 2πTl,pTl+1,p e−(vl−ul,p)2 2Tl,p −(vl+1−ul+1,p )2 2Tl+1,p .(7) We also assume that ul,p and Tl,p are smooth functions of l in the large system size limit L1, more specifically that the discrete derivatives ul+1,p −ul,p =O(L−1) and Tl+1,p − Tl,p =O(L−1), which gives (vl,p −vl+1,p)(vl,p −vl+1,p)∼Tl,p/π. (8) This is a “modified” pressure; the usual ideal gas equation of state ρT is replaced by √Tbecause (i) collisions do not occur at a rate proportional to √T, but at constant rate and (ii) density ρis uniform, ρ=1. Therefore, ul,p+1−ul,p ∼−1+α 2L(Tl,p/π −Tl−1,p/π).(9) Now we introduce a continuum limit by defining spatial and temporal variables as x=l, t =2p,  =L−11,(10) to obtain ∂tu(x,t)=−1+α 2∂xT(x,t)/π +O().(11) We now discuss the average of Eq. (5) for the kinetic energy, which contains both dissipative and transport terms. Within the local equilibrium approximation, the average of the dissipative or cooling term is α2−1 4L(Tl,p +Tl−1,p)=(α2−1) 2Tl,p +O((α2−1)2), (12) which in the continuum limit reads (α2−1) 2T(x,t)+O((α2−1)2).(13) This linear cooling, instead of the typical T3/2behavior, also stems from the constant collision rate. The average of the energy current, again under local equilibrium, reads Jl,p=(1 +α)u l,pTl,p/π +O(2), so that the average of the last two terms in Eq. (5) goes in the continuum limit to −(1 +α)2∂x[u(x,t)T(x,t)/π]+O(3).(14) Comparing Eqs. (13) and (14), we see that they are of the same order only when α2−1=O(): this choice makes them of order 2, consistently with the scaling of the finite-time difference on the left-hand side of Eq. (5), t =2. Therefore, we introduce a “macroscopic inelasticity” ν=L(1 −α2)/2⩾0.(15) Taking into account that v2=u2+T, we get ∂tT(x,t)=−νT(x,t)−2T(x,t)/π ∂xu(x,t)+O(), (16) where we have also made use of 1 +α∼2forL1 with constant ν.Moreover,Eq.(11) becomes ∂tu(x,t)=−∂xT(x,t)/π +O().(17) 062905-3 BALDASSARRI, PUGLISI, AND PRADOS PHYSICAL REVIEW E 97, 062905 (2018) Equations (17) and (16) are the hydrodynamic equations of the model to the lowest order. They coincide with the hydrodynamic Eqs. (1) after identifying Dt→∂t(Lagrangian coordinate), ρ→1, P→√T/π[Eq. (8)], ζ→ν, and ∂xq= 0. The absence of dissipative transport (viscous stress and heat flow) makes these equations close relatives of the so-called ideal granular hydrodynamics [56,57,73]. C. Steady uniform longitudinal flow In the large-system size limit, the stationary solution of Eqs. (17) and (16) can be obtained. It suffices to ask that the velocity flow in the middle of the system is 0 [74] to get the following steady ULF profiles: Ts(x)=T0,u s(x)=ν 2πT01 2−x,(18) in which the constant T0remains undetermined at this level of description [75]. The steady ULF in Eq. (18) does not require the presence of thermostats at the boundaries; it is selfsustained. Nevertheless, as demonstrated numerically below, one can recover the above steady ULF profile at the system bulk in a finite-size system by attaching two identical thermostats at the boundaries. Also, note that this ULF, at difference with the one discussed in Ref. [38], is incompressible because no mass flow is allowed in our lattice system by definition. Quite strikingly, regions with a linear velocity profile with negative slope are also observed in the Lagrangian coordinate during the formation of shocks, in numerical simulations of the inelastic 1d gas of hard rods [15,34], in the inelastic hard disk gas in a channel [55], and in the lattice model [15], all in the cooling regime at long times. This suggests an interesting connection between the ULF and the shocks characterizing granular cooling in 1d (or quasi-1d) systems. IV. NUMERICAL RESULTS Numerical simulations of the model have been carried out and compared with the previous theory. A system with N particles is initialized with uncorrelated, normally distributed random velocities. We introduce thermostats at its boundaries, at sites 0 and N+1, to make the system reach a steady state. Then, at each time step p, the uniformly distributed random integer ypchoosing the candidate colliding pair {yp,yp+1}is drawn between 0 and N, that is, L=N+1. The velocities for the peripheral sites of indexes 0,L are randomly and independently drawn from a normal distribution with zero average and unit variance [76]. Due to the kinematic constraint, the collision described by Eq. (2) only takes place if vyp− vyp+1>0. Otherwise, velocities remain unchanged. A large number of independent long runs has been performed to get average stationary profiles of velocity and energy. Typically, our data correspond to 105runs of more than 102L2steps each, starting from independent normally distributed velocities. Our numerical results for the average temperature profiles are shown in Fig. 2for different values of the macroscopic inelasticity ν. From them, we measure the temperature of the largest system at mid-position, denoting it by T1/2≡T(x= 1/2). Later, we compare T1/2with the bulk temperature T0 FIG. 2. Average profile for the steady temperature, for different values of the macroscopic inelasticity ν. The solid (red) lines stand for the estimated value of the temperature T1/2in the system bulk. of Eq. (18), which predicts a constant-temperature profile, in quite good agreement with our observed numerical profiles. Numerical results for the average velocity profiles are reported in Fig. 3. Again, the comparison with the theory is satisfactory: in the bulk, the profiles are linear with a slope almost perfectly matching the one given by the value of T1/2 in Fig. 2. Both for the temperature and the average velocity, deviations from the theoretical ULF profiles are apparent near the boundaries, especially for ν=2. See below for a closer look at this issue. V. FINITE-SIZE CORRECTIONS AND BOUNDARY LAYERS In Appendix A, we show how to incorporate finite-size corrections. Still, we do so assuming local equilibrium, in order FIG. 3. Average profile for the stationary velocity, for the same values of νconsidered in Fig. 2. The solid (red) lines stand for the theoretical expression for us(x)inEq.(18), with T0givenbythe numerical estimate T1/2for the bulk temperature in Fig. 2. 062905-4 HYDRODYNAMICS OF GRANULAR PARTICLES ON A LINE PHYSICAL REVIEW E 97, 062905 (2018) FIG. 4. Average temperature profiles T(x) in the stationary state for different values of ν(blowup of data in Fig. 2). Red lines stand for the numerical estimate T1/2for the bulk temperature, see Fig. 2, whereas black lines give our analytical prediction T0,asgivenby Eq. (20). to estimate in the simplest way the effect of these higherorder terms. Specifically, these O() corrections introduce second-order spatial derivatives in the evolution equations, and these viscous terms make it possible to accommodate all the boundary conditions. The considered thermostats impose that u(0,t)=u(1,t)=0,T(0,t)=T(1,t)=1.(19) The smallness of the viscous terms brings to bear two boundary layers close to the thermostats at x=0,1; see Appendix Bfor details. Therein, we determine the unknown bulk temperature T0by asymptotic matching [77], which yields ϕ(T0)=√π 4ν, with ϕ(T0)≡1 3+2 3T−3/2 0−T−1 0.(20) The function ϕ(T0) decreases monotonically from infinity to zero as T0is varied from zero to unity. Therefore, Eq. (20) tells us that T0is a monotonically decreasing function of the macroscopic inelasticity ν, with lim ν→0T0=1,lim ν→∞T0=0.(21) These two limit results are expected on a physical basis: in the elastic limit ν→0, the system should be in equilibrium at the temperature of the heat baths whereas in the strongly dissipative limit ν→∞all the energy is dissipated in the boundary layers before reaching the bulk, which is then at zero temperature [78]. A more precise inspection of the temperature profiles is shown in Fig. 4, aimed to emphasize the deviation from the uniform bulk profile in Eq. (18). Interestingly, the estimate for the bulk temperature T0stemming from our finite-size analysis, as given by Eq. (20), is quite close to our numerical estimation T1/2. Most importantly, the success of our theory is good for very different values of ν, both small and large. The absolute error |T0−T1/2|increases very slowly with ν, ranging from 0.01 for ν=0.2to0.02forν=20. The relative error is also reasonable, remaining under ten percent for ν2. In Appendix C, we look into possible sources of discrepancy between simulations and theory. We focus on violations of local equilibrium, Eq. (7): both non-Gaussian local velocity distributions and nearest-neighbor correlations. Our preliminary conclusion is that the former are more important than the latter. VI. CONCLUDING REMARKS We have analyzed a 1d lattice model, which by construction does not contain density fluctuations. Our derivation of hydrodynamics can be thought as if carried out in an off-lattice system but using a Lagrangian coordinate at the microscopic level. The resulting Eqs. (17) and (16), in which finite-size corrections have been neglected, tell us that a self-sustained ULF steady state appears in the bulk. Therein, the granular temperature is uniform and the average velocity has a linear profile, whose slope is directly related to the temperature. In order to observe such a steady state in a finite-size system, thermostats at the boundaries must be introduced, as shown in our numerical simulations which fairly compare with our theory. The dependence of the bulk temperature on the temperature of the thermostats is only obtained after incorporating finite-size corrections to the theory. Our work demonstrates that 1d granular hydrodynamic equations, whose validity is still under debate, can be derived by means of a proper continuum limit that keeps a finite macroscopic inelasticity, but in which the microscopic dynamics is quasielastic. Our main assumption is the local equilibrium approximation, which has been shown to be valid in other simple models for a wide range of the system parameters [60,61,79]. How to improve upon the current results by going beyond local equilibrium is an open perspective for future work. An emerging interesting conjecture is whether the system in the cooling regime is still well described by our hydrodynamic equations. In fact, the ULF discussed here seems to be similar to the linear velocity profiles in shock regions of 1d cooling granular gases, specifically for either hard rods on a line [34] or hard disks in a channel geometry [55]. A verification of this conjecture, together with a route to connect our theory with the Burgers or flow-by-inertia equations, is currently under investigation. ACKNOWLEDGMENTS We acknowledge useful discussions with Dario Villamaina. Antonio Prados acknowledges the support of Universidad de Sevilla’s VI Plan Propio de Investigación through Grant No. PP2018/494. APPENDIX A: HYDRODYNAMIC EQUATIONS WITH FINITE-SIZE CORRECTIONS Here we derive the hydrodynamic equations, incorporating into them finite-size, O(L−1), corrections. First, we analyze the evolution equation for the velocity: our starting point is Eq. (6), in which the average momentum current at site lis exactly 062905-5 BALDASSARRI, PUGLISI, AND PRADOS PHYSICAL REVIEW E 97, 062905 (2018) given as jl,p=(1 +α) 2(vl,p −vl+1,p)(vl,p −vl+1,p),=L−1,(A1) since ypis an independent, uniformly distributed, stochastic integer choosing the specific pair that collides. We recall that Lis the total number of possible colliding pairs, whose relation to the number of sites Ndepends on the boundary conditions. The average on the right-hand side of the above equation is done over all the velocities {vl,p,vl+1,p}, by assuming the local equilibrium approximation written in Eq. (7) of the main text. The result is (vl,p −vl+1,p)(vl,p −vl+1,p)=ul,p −ul+1,p 2erful,p −ul+1,p 2(Tl+1,p +Tl,p)+1+exp −(ul,p −ul+1,p)2 2(Tl+1,p +Tl,p)Tl+1,p +Tl,p 2π, (A2) where erf(z) is the error function defined by erf(z) =2 √πz 0 dt e−t2.(A3) Now we go to a continuum limit in space by assuming that both ul,p and Tl,p vary smoothly with the site index l. Then, we define a continuous spatial variable x=l, with =L−11, and make the mapping ul,p →u(x;p), Tl,p →T(x;p) and, consistently, ul±1,p →u(x±;p), Tl±1,p →T(x±;p). Of course, this is also done for the current, jl,p →j(x;p) and jl±1;p→j(x±;p). In addition, we expand all the functions evaluated at x±in powers of the small parameter , which transforms the exact momentum balance Eq. (6)into u(x;p+1) −u(x;p)=−(1 +α)2 2∂xT(x;p) π− 2∂2 xu(x;p)+O2.(A4) It is clearly seen that the evolution of u(x;p), being proportional to 2, is very slow in discrete time. This suggests the introduction of a continuous timescale t=pt, with t =2[80]. Over this scale, u(x;p+1) −u(x;p)=t ∂tu(x,t)+O(4) and we thus finally have ∂tu(x,t)=1+α 2−∂xT(x,t) π+ 2∂2 xu(x,t)+O2.(A5) Next, we repeat the above procedure for the balance of energy equation. The corresponding expressions are much lengthier than those for the average velocity above and thus we only give the final expressions for the averages of the dissipative and flux terms on the right-hand side of Eq. (5). For the dissipative term, we have (α2−1) 4(vl,p −vl+1,p)(vl,p −vl+1,p)2+(vl−1,p −vl,p)(vl−1,p −vl,p)2 =(α2−1) 2T(x;p)−2T(x;p) π∂xu(x;p)+O(α2−1)3,(A6) again within the local equilibrium approximation. For the flux terms, Jl,p −Jl−1,p=(1 +α)2∂xu(x;p)T(x;p) π−1 4∂2 xT(x;p)+2∂x[u(x;p)∂xu(x;p)]+O2.(A7) The evolution equation for the energy eis then e(x;p+1) −e(x;p)=(α2−1) 2T(x;p)−2T(x;p) π∂xu(x;p) −(1 +α)2∂xu(x;p)T(x;p) π−1 43∂2 xT(x;p)+2∂x[u(x;p)∂xu(x;p)] +O((α2−1)3)+O(4).(A8) Going again to the continuous time variable, ∂te(x,t)+O(2)=α2−1 2T(x,t)−2T(x,t) π∂xu(x,t)−(1 +α)∂xu(x,t)T(x,t) π +1+α 4∂2 xT(x;p)+2∂x[u(x;p)∂xu(x;p)]+O(2)+O(α2−1).(A9) 062905-6 HYDRODYNAMICS OF GRANULAR PARTICLES ON A LINE PHYSICAL REVIEW E 97, 062905 (2018) In order to have a consistent limit over the continuous time t,α2−1 must be of the order of . Hence, we introduce the macroscopic inelasticity ν, as defined in Eq. (15), which is assumed to be of order unity. This implicitly assumes that the underlying microscopic dynamics is quasielastic, since 1 −α∼ν =ν/L 1[81]. With this definition, ∂te(x,t)=−νT(x,t)−(2 −ν)∂xu(x,t)T(x,t) π+2νT(x,t) π∂xu(x,t)+1 2∂2 xT(x,t)+∂x[u(x,t)∂xu(x,t)]+O(2). (A10) Finally, the hydrodynamic equations for the average velocity and the temperature up to order =L−1are obtained by taking into account the definition of the macroscopic inelasticity (15) and the identity e(x,t)=u2(x,t)+T(x,t), which yields ∂tu(x,t)=−∂xT(x,t) π+ 2ν∂ xT(x,t) π+∂2 xu(x,t)+O2,(A11a) ∂tT(x,t)=−νT(x,t)−2T(x,t) π∂xu(x,t)+3νT(x,t) π∂xu(x,t)+1 2∂2 xT(x,t)+[∂xu(x,t)]2+O(2).(A11b) Note that the parameter can be understood as the ratio of a microscopic length (1 lattice site) to a macroscopic length (the total length L). In this way, the condition 1 can be conceived as the usual small Knudsen number condition for the validity of a continuum description. A more detailed look at this issue can be found in the following appendix. APPENDIX B: BOUNDARY LAYER CALCULATIONS Now we solve the hydrodynamic equations to the lowest order, but incorporating the boundary layers close to the system edges. Thus, we consider Eqs. (A11) with the boundary conditions in the half interval [0,1/2]: u(0,t)=0,T(0,t)=1,(B1a) u(1/2,t)=0,∂ xT(x,t)|x=1/2=0.(B1b) The boundary conditions at x=1/2stemfromu(x,t) and T(x,t) being an odd and an even function with respect to x= 1/2, respectively. We focus on stationary solutions of the hydrodynamic equations. Then, we can simplify the notation by introducing ≡∂x. Using the terminology in Ref. [77], the “outer” solutions {uout(x),Tout(x)}satisfy the equations Tout π =0,−νTout −2Tout πu out =0.(B2) These equations are first order in space and thus we can only impose the boundary conditions at one of the end points of the half interval [0,1/2], specifically those at the system center x=1/2. Then, Tout =T0,u out(x)=−ν 2πT 0x−1 2,(B3) where T0remains undetermined. The explicit form of the outer (bulk) solutions allows us to have a more precise look at the necessary condition to have a meaningful continuum description. Making use of Eq. (B3), we can readily identify a macroscopic length ξ=(ν√π2)−1 [the inverse of the factor in front of the expression for uout(x), apart from √T0that carries with it the velocity dimensions]. The microscopic length scale in the continuum is =L−1, corresponding to one lattice site, and then the sought condition would be  ξ=ν√π 2L=(1 −α2)√π 41,(B4) where we have again used the definition of ν,Eq.(15). This is nothing but a condition of quasielasticity for the microscopic dynamics, a restriction that is already known for other lattice models [60,61,79]. Several comments are in order. First, note that this condition has already been imposed when we have derived the hydrodynamic equations, specifically when we have assumed that the macroscopic inelasticity νis of the order of unity. Second, this is physically reasonable: the magnitude of the gradients in granular systems is often typically controlled by the inelasticity, not by the imposed boundary conditions. Therefore, the microscopic inelasticity 1 −αmust be small in lattice models to ensure that the gradients are not so large that the lattice is unavoidable “seen” and a continuum description is not possible. Third, despite the quasielasticity of the microscopic dynamics, the observed macroscopic behavior does not correspond to a quasielastic granular fluid, because the macroscopic inelasticity ν(also called the cooling rate) is not small. A different distinguished limit is obtained within the boundary layer, with the scaling x=X, d dx =−1d dX.(B5) It is Xthat the “inner” solution within the boundary layer, which we denote by {uinn(X),Tinn(X)}, depends on. The dominant terms in the balance equations become −d dXTinn π+1 2 d2 dX2uinn =0,(B6a) −2Tinn π d dXuinn +1 2 d2 dX2Tinn +d dXuinn2 =0.(B6b) 062905-7 BALDASSARRI, PUGLISI, AND PRADOS PHYSICAL REVIEW E 97, 062905 (2018) We have to solve the above equations with the boundary conditions (B1a)atX=0, that is, uinn(X=0) =0,T inn(X=0) =1,(B7) and the matching conditions [77] lim X→∞uinn(X)=lim x→0uout(x)=ν 4πT0,(B8a) lim X→∞Tinn(X)=lim x→0Tout(x)=T0.(B8b) The latter conditions assure that the solution in the boundary layer smoothly match the outer solution in the bulk. These matching conditions make it possible to determine the unknown bulk temperature T0, as shown below. Equation (B6a) is integrated straightforwardly to give d dXuinn −2Tinn π=−2T0 π,(B9) where the constant on the right-hand side has been obtained by matching the outer and inner solutions to the lowest order; note that duinn/dX →0forX→∞. Substitution of (B9)into (B6b) yields, after some simple algebra, d2 dX2θ+8 π(1 −√θ)=0,θ=Tinn T0 .(B10) A first integral can be directly derived from this equation, 1 2d dXθ2 +8 πθ−2 3θ3/2=D, (B11) where Dis a constant. Again, Dis obtained by matching arguments in the limit as X→∞, for which we have that θ→1 and dθ/dX →0. Therefore, D=8 3πand d dXθ=− 4 √π1 3+2 3θ3/2−θ. (B12) We have chosen the minus sign on the right-hand side because dθ/dX must be negative in the boundary layer, since the bulk temperature T0<1 as a consequence of the dissipative character of the dynamics. We do not need to solve Eq. (B12) to obtain T0, which is our main goal. Going back to Eq. (B9), we can rewrite it as d dXuinn =2T0 π(√θ−1) (B13) and combining it with Eq. (B12), duinn =−√T0 2 √θ−1 1 3+2 3θ3/2−θ dθ. (B14) This equation allows us to calculate uinn as a function of θ, taking into account that uinn(X=0) =0 and θ(X=0) = T−1 0, uinn(θ)=T01 3+2 3T−3/2 0−T−1 0−1 3+2 3θ3/2−θ. (B15) Now we impose the matching conditions in Eq. (B8), that is, {uinn →ν√πT0/4,θ →1}in the limit as X→∞. Hence, one can write that √π 4ν=ϕ(T0),ϕ(T0)=1 3+2 3T−3/2 0−T−1 0.(B16) The function ϕ(T0) on the right-hand side is a monotonically decreasing function of T0(recall that 0 ⩽T0⩽1), lim T0→0ϕ(T0)=∞,lim T0→1ϕ(T0)=0,dϕ(T0) dT0 <0.(B17) Thus, Eq. (B16) is the desired expression for the bulk temperature as a function of the macroscopic inelasticity ν, since it univocally gives T0for each value of ν. In Fig. 4, a good agreement is shown between the prediction for T0obtained here and numerical simulations. The discrepancies remain quite small for ν2, becoming only larger for the highly dissipative case ν=20. It must be taken into account that the bulk temperature T0has been assumed to be of order unity in our theory, whereas Eq. (B16) implies that it becomes very small for high ν, specifically T0∼32 3π2/3 ν−4/3,ν1.(B18) This means that a more elaborate theory, corresponding to a different dominant balance in the hydrodynamic equations, might be necessary in the highly dissipative limit ν1. See also the next appendix for other possible sources of discrepancy between our theory and the numerical results as νincreases. One can also examine the role of the boundary layers in Fig. 4. The boundary layers are barely noticeable for the smaller values of ν,ν=0.2 and ν=0.02, for which the bulk temperature is close to unity. In addition, these small values of νnecessarily bring about a stronger noise in the averages, since the bulk temperature is larger. As νincreases, the bulk temperature decreases and so do the fluctuations while the boundary layers become more visible. In Figs. 2and 4, it seems that the boundary layer is wider for ν=2 than for ν=20. This is reasonable and consistent with the behavior observed in other models in the strongly dissipative limit ν1[79,82], in which the width of the boundary layer algebraically decreases with νand vanishes in the limit as ν→∞. Besides the above calculations show that, as a function of the system size, the width of the boundary layer in the xvariable is expected to be of order L−1. This scaling is not clearly confirmed in our simulations, which seem to indicate that this width does not go to zero in the limit as L→∞, especially for the highest macroscopic inelasticity ν=20. This discrepancy might also be mended by a more elaborate theory in the highly dissipative case; see also next section. Finally, some kind of “boundary resistance” should also be noted, that is, a difference between the actual value of T at the boundary and the value imposed by the thermostat in the simulation (T=1). This phenomenon is known to appear in nonlinear transport problems [83,84] and has been already observed in other models [79]. APPENDIX C: NON-GAUSSIAN LOCAL DISTRIBUTIONS AND SPATIAL CORRELATIONS Here we investigate the validity of the local equilibrium approximation. A first check is obtained by measuring the local 062905-8 HYDRODYNAMICS OF GRANULAR PARTICLES ON A LINE PHYSICAL REVIEW E 97, 062905 (2018) FIG. 5. Local velocity distributions P(v,x) in the stationary state for different values of νand different positions. All the plots corresponds to the largest considered system size N=1600. Top row: Position close to the left boundary, x→0andx=0.1. Bottom row: Bulk position x=0.5. velocity pdf. This is done in Fig. 5. In the top row, we consider positions close to the boundaries, and deviations from Gaussianity are already apparent for ν=0.2, with an anomalous but localized peak at v=0. These discrepancies become even more patent as νincreases, with the emergence of asymmetric tails in the pdf. The pdf in the bulk of the system is presented in the bottom row of the same figure. Again, discrepancies arise and are evident, with an analogous localized peak at v=0for ν=0.2 that splits into two symmetric peaks for larger ν1 [85,86]. Note that the tails remain symmetric with respect to v=0 in the bulk, however. Our conclusion is that some of the observed discrepancies as νincreases stem from the non-Gaussianities described above. Certainly, the non-Gaussianities are more important at the boundary layers than in the bulk, where the deviations seem to be milder, especially for the not-so-large inelasticity ν=2. Anyhow, non-Gaussianities may be responsible for the deviations of the actual bulk temperature T1/2from the theoretical description T0for ν1. It has to be taken into account that we have assumed that the local equilibrium approximation holds both in the bulk and at the boundary layers for deriving Eq. (20). It is remarkable that the relative error between T0and T1/2remains under ten percent for ν=2, despite the large discrepancies at the boundaries that include asymmetry with respect to v=0. Also, we have looked into the nearest-neighbor correlations dl=vl−1vl−vl−1vl, which have been assumed to vanish upon writing the local equilibrium approximation (7). This assumption is consistent with the molecular chaos hypothesis, in which these correlations are assumed to be of order L−1. Taking into account that we have incorporated O(L−1) corrections into our theory, these correlations are another possible source for discrepancies and should be investigated. In Fig. 6, the numerical evaluation of the nearest-neighbor correlations is displayed. It is clearly seen that dlis always different from zero FIG. 6. Nearest-neighbor velocity correlation dlin the stationary state as a function of x=l/L, for different system sizes while keeping constant the macroscopic inelasticity ν=19. 062905-9