Understanding an instability in vibrated granular monolayers

PHYSICAL REVIEW E 99, 032903 (2019) Understanding an instability in vibrated granular monolayers P. Maynar,*M. I. García de Soria, and J. Javier Brey Física Teórica, Universidad de Sevilla, Apartado de Correos 1065, E-41080, Sevilla, Spain (Received 15 June 2018; revised manuscript received 21 January 2019; published 21 March 2019) We investigate the dynamics of an ensemble of smooth inelastic hard spheres confined between two horizontal plates separated by a distance smaller than twice the diameter of the particles, in such a way that the system is quasi-two-dimensional. The bottom wall is vibrating and, therefore, it injects energy into the system in the vertical direction and a stationary state is reached. It is found that if the size of the plates is small enough, the stationary state is homogeneous. Otherwise, a cluster of particles is developed. The instability is understood by using some effective hydrodynamic equations in the horizontal plane. Moreover, the theoretical prediction for the size of the system above which it becomes unstable agrees very well with molecular dynamics simulation results without any fitting parameter. DOI: 10.1103/PhysRevE.99.032903 I. INTRODUCTION A granular system is an ensemble of macroscopic particles, named grains, that interact via inelastic collisions. When two grains collide, because of their macroscopic character, part of the kinetic energy associated to the center of mass of the particles is dissipated, exciting other degrees of freedom. Hence, besides the fact that granular matter is ubiquitous in nature, these systems are very interesting from a theoretical point of view because they are intrinsically out of equilibrium [1]. Experimentally, it is relatively easy to generate stationary states by just vibrating the walls of the container of the system. In these states, the energy injected in the system is compensated by the energy dissipated at collisions. The price to be paid is that the steady state is highly inhomogeneous, as can be understood from the generalized Navier-Stokes equations describing the dynamics of the system [2,3]. An exception to this is a granular system placed in a vertically vibrated quasi-two-dimensional (Q2D) shallow box. In this case, the stationary state can be considered as homogeneous when projected in the horizontal plane and, if observed from above, the system behaves as a two-dimensional fluid. Actually, in recent years, many experiments have been performed under these conditions [4–12]. Experiments can be carried out with or without a top lid, with gravity in the last case being responsible for the Q2D confinement. Typically, the system is vibrated sinusoidally in the direction perpendicular to the plates with frequency wand amplitude Athat is always much smaller than the height of the system, H, i.e., AH.Mostof the experiments are performed with smooth plates, although in some cases a rough one is used [5]. On average, particles gain energy in the vertical direction through collisions with the vibrating walls and it is transferred to the horizontal degrees of freedom via collisions between particles, where energy is also dissipated as a result of inelasticity. Experiments show that, for a wide range of the parameters, the system reaches a homogeneous stationary state. On the other hand, it is also *[email protected] observed that the system becomes unstable by increasing the density and/or varying the parameters controlling the vibrating walls. In Refs. [4,13], in which there is no top lid and the Q2D confinement is due to gravity, it is found that the instability depends on the dimensionless parameter ≡Aw2 g,with g being the gravitational acceleration. On the other hand, when the confinement is due to the presence of a top lid and the effect of gravity can be neglected, as in Refs. [6,7], the relevant parameter is the dimensionless amplitude A/σ . In both cases, there is a regime in which an aggregate surrounded by a more dilute hotter phase is formed. Let us also mention that, depending on the averaged density, the coexistence is between a solidlike phase and a liquidlike phase [4,12] or between a liquid and a gas [6,7]. It is worth pointing out that both situations probably correspond to fully different phase transitions. In recent years, there have been many efforts to understand this phenomenology. For example, it is known that, in the context of hydrodynamics, the instabilities are triggered by a negative compressibility in the associated horizontal equation of state [6,12,13] and some of the phases have been characterized [9,10]. Nevertheless, there are still many points that are not clear. In particular, what are the essential ingredients leading to this kind of instability? Why is the compressibility negative? Can the equation of state in the horizontal plane be derived from a microscopic point of view? In any case, what is the mechanism that stabilizes and destabilizes the system? In order to tackle these questions, some two-dimensional (2D) models have been considered. As the system is Q2D, it is expected that the actual three-dimensional (3D) dynamics could be substituted by an effective 2D dynamics in which energy is injected by some kind of mechanism as, for example, in the the stochastic thermostat [14]orintheso-called model [15]. Nevertheless, both models have been extensively studied, and it has been found that the homogeneous stationary state that is reached in the long time limit is always stable [16–18]. Very recently, a new 2D model has been formulated in which the homogeneous stationary state is unstable [19], but the proposed microscopic dynamic is defined in terms of some parameters that, in principle, must be fitted. 2470-0045/2019/99(3)/032903(6) 032903-1 ©2019 American Physical Society MAYNAR, GARCÍA DE SORIA, AND BREY PHYSICAL REVIEW E 99, 032903 (2019) Z σ vp H FIG. 1. Schematic representation of the model. The objective in this paper is to introduce a simple model that captures the phenomenology of the experiments and that lets us address the aforementioned questions in the simplest situations. We will consider that the system is confined between two walls and gravity can be neglected. We will also restrict ourselves to sufficiently dilute systems. In this sense, our theory is closer to describe the experimental setup considered in Refs. [6,7], but probably requires some modifications for the description of the one of Refs. [4,13]. The remainder of the paper is organized as follows. In the next section, the model is introduced. It is shown that, for a wide range of the parameters, the homogeneous stationary state that is reached is stable and that it turns unstable when increasing the density. In Sec. III, the homogeneous phase is studied and characterized by the horizontal and vertical temperatures. The stability of this homogeneous phase is studied in Sec. IV and the origin of the instability is clarified. Finally, Sec. Vcontains some final comments and a summary of the conclusions. II. THE MODEL To properly describe the energy transfer from the vertical to the horizontal direction, we consider a simple 3D model: an ensemble of Ninelastic hard spheres of mass mand diameter σ, confined between two flat planes located at z=0 and z=H. It will be assumed that H<2σso that particles cannot jump on to each other and the system can be considered to be Q2D (see Fig. 1). The plates are square shaped of area L×Lwith LHand periodic boundary conditions are used in the horizontal directions. The collisions between particles are inelastic, characterized by a constant (independent of the relative velocity) coefficient of normal restitution, α(0 <α⩽ 1, being the collisions elastic for α=1). The collisions with the two walls are elastic and the top wall is at rest, so that the horizontal components of the velocity of a particle colliding with it does not change and vz→−vz. The bottom wall will be modeled as a sawtooth wall of velocity vp, i.e., all particles colliding with the wall find it with a constant velocity vp in the direction of positive z. The position of the wall can be taken as fixed since the amplitude of the wall motion is assumed to be very small as compared with H. Hence, when a particle collides (elastically) with the bottom wall, its horizontal velocity remains unchanged, while vz→2vp−vz. This kind of collision always injects energy into the system. Note that friction between particles and the wall is neglected and the total horizontal momentum of the system is a constant of the motion. In this way, experimental situations in which friction of the particles with the walls plays a dominant FIG. 2. Snapshot of a configuration of a MD simulation where the cluster shows up (the system is observed from above). The values of the parameters of the system are given in the main text. role are beyond the range of application of the theory being developed here. We have performed molecular dynamics (MD) simulations of the model, finding that it exhibits the main features observed in the experiments: For a wide class of values of the parameters, the system reaches a homogeneous stationary state and, when increasing the density at fixed size, a dense cluster surrounded by a gas is developed. In Fig. 2, a snapshot of a configuration of a MD simulation where the cluster is formed is shown (the system is seen from above). The parameters of the simulations are N=2000, α=0.9,H= 1.5σ, L=115σ(so that the effective 2D density is N L2≈ 0.15σ−2), and vp=0.01[2T(0) m]1/2, where T(0) is the initial horizontal granular temperature. The granular temperature is defined as usual by setting the Boltzmann constant equal to unity [1]. Note that the results are similar to the ones of Refs. [6,7] in which a liquidlike phase surrounded by a hotter gas is observed, although in these references the height of the system is of the order of 10σ. III. THE HOMOGENEOUS PHASE To study the homogeneous dilute phase, we will admit, in the same spirit as in Ref. [13], the existence of a closed description in terms of the 2D density, n, the horizontal temperature, T, and the vertical temperature, Tz. The idea is that since energy is injected in the vertical direction and transferred to the horizontal directions through collisions, this is the minimal number of variables to understand the dynamics of the homogeneous system. Total momentum in the plane does not play any role in this context, as it is a constant of the motion. In the low-density limit, assuming that the system is very thin, i.e., ≡H−σ σ1, and that collisions between particles are nearly elastic, 1 −α1, the evolution equations are dT dt =2π mnσT1 2−(1 −α)T+2 3(Tz−T),(1) dT z dt =−4 3π m2nσT1 2(Tz−T)+2vp σ Tz.(2) 032903-2 UNDERSTANDING AN INSTABILITY IN VIBRATED … PHYSICAL REVIEW E 99, 032903 (2019) Let us discuss briefly each term in the above equations. The first term on the right-hand side of Eq. (1) describes the energy loss due to inelastic collisions between particles in a 2D system [1]. The second term in the same equation is the kinetic energy transfer from the vertical degree of freedom to the horizontal ones that is proportional to the difference of temperatures times the thermal horizontal velocity. We have taken the elastic limit reported in Ref. [20] from kinetic theory tools. This term also appears in Eq. (2) but with a different sign and multiplied by 2 due to the difference in the definition of vertical and horizontal temperatures in terms of the kinetic energy. Finally, the second term in the right-hand side of Eq. (2) describes the energy injection due to the vibrating wall. In this case, the term can be evaluated exactly as it is proportional to vptimes the pressure of the granular gas just above the vibrating wall in the direction perpendicular to it [21]. In fact, these equations can be derived starting from a Boltzmann kinetic equation for the monolayer and carrying out an expansion in approximating the velocity distribution function by a Gaussian with two temperatures [22]. The only difference is that the energy transfer terms have a nontrivial dependence on αneglegible in the quasielastic limit being considered here. Equations (1) and (2) admit only one stationary solution, Tsand Tz,s, that can be easily calculated as a function of the inelasticity, density, and dimensionless heigh, .FromEq.(1), the quotient of the stationary temperatures, γ≡Tz,s Ts,is γ=1+3(1 −α) 2(3) that is density independent. Moreover, it is always equal or greater than one, with equipartition holding only in the elastic limit. The stationary horizontal temperature is now obtained from Eq. (2), Ts=m 4πγvp (1 −α)nσ22 .(4) Remarkably, it is proportional to n−2, so that the pressure in the stationary state, defined as p(n)≡nTs(n), goes as n−1, and the compressibility of the nonequilibrium steady state is negative. Note that this dependence is a direct consequence of the particular way in which energy is injected and transferred in the monolayer. The divergence of the stationary temperature for n→0 comes from the fact that energy is always injected in a particle-wall collision. To test the accuracy of the predictions given by Eqs. (3) and (4), we have performed MD simulations, measuring the horizontal and vertical temperatures. The simulations were run for systems with N=500,nσ2=0.03, vp=0.01[2T(0) m]1/2, two different values of the height, =0.5 and =0.2, and different values of the coefficient of normal restitution in the range 0.7⩽α⩽0.95. The initial condition was taken to be a Gaussian with two temperatures. For these values of the parameters, it is found that, independently of the initial condition, the system always reaches a homogeneous stationary state. In Fig. 3, the quotient of the stationary temperatures, γ, measured in the simulations is plotted as a function of the inelasticity for =0.5 (circles) and =0.2 (squares). The solid and dashed lines are the corresponding theoretical 0.7 0.8 0.9 1 α 0 10 20 30 γ FIG. 3. Ratio of the stationary temperatures, γ≡Tz,s Ts, measured in the simulations as a function of the inelasticity for =0.5 (circles) and =0.2 (squares). The solid and dashed lines are the corresponding theoretical predictions given by Eq. (3). predictions given by Eq. (3). As can be seen, the agreement between the MD simulation results and the theoretical prediction is excellent and only some discrepancies are obtained for the stronger inelasticities. This was expected as the theoretical predictions are obtained for mild inelasticities. In Fig. 4,the stationary horizontal temperature scaled with mv2 pis plotted as a function of the inelasticity for =0.5. The circles are the simulation results and the solid line is the theoretical prediction given by (4). The same is represented in the inset for =0.2. Again, the agreement between the simulation results and the theoretical prediction is very good. This is remarkable since there is not any adjustable parameter. Notice that the scale of the dimensionless quantity Ts mv2 pis of order 105 and 107for =0.5 and =0.2 respectively. This is because there are two small dimensionless parameter, and nσ2, that contribute to the increase of the stationary temperature, as can be seen in Eq. (4). A stability analysis of Eqs. (1) and (2) shows that the stationary solution is linearly stable [22]. In fact, the matrix 0.7 0.8 0.9 1 α 0 1×105 2×105 3×105 4×105 5×105 Ts/(mvp 2) 0.7 0.8 0.9 1 α 1.2×107 1.6×107 2.0×107 Ts/(mvp 2) FIG. 4. Stationary horizontal temperature scaled with mv2 pas a function of the inelasticity for =0.5. The circles are the simulation results and the solid line is the theoretical prediction given by (4). The same is represented in the inset for =0.2. 032903-3 MAYNAR, GARCÍA DE SORIA, AND BREY PHYSICAL REVIEW E 99, 032903 (2019) associated to the linear dynamics of δT≡T−Tsand δTz≡ Tz−Tz,shas two eigenvalues, one of which is always the slowest one and, in addition, goes to zero in the elastic limit. Then, there is a timescale in which the dynamics is dominated by this slowest mode and a description in terms of only one temperature, T, is valid. In the following, we will call it “homogeneous hydrodynamic” timescale. In this timescale, the evolution equation for the temperature is d dsδT≃−βδT,(5) where we have introduced the dimensionless time ds = nσ2Ts mdt and βis the absolute value of the slowest eigenvalue that is a known function of αand . We do not give its explicit expression here since it is not very illuminating. IV. THE NONHOMOGENEOUS CASE Now, let us suppose that there are gradients in the horizontal plane. The previous results make plausible a description in terms of the 2D density, n(r,t), 2D flow velocity, u(r,t), and horizontal temperature, T(r,t). The evolution equations for these fields are assumed to be of the form ∂n ∂t=−∇·(nu),(6) ∂u ∂t=−u·∇u−1 nm ∇·P,(7) ∂T ∂t=−G(n,T)−u·∇T−1 n(∇u:P+∇·q),(8) where Pis the pressure tensor and qis the heat flux. In the low-density limit and to Navier-Stokes order, we assume Pij =nT δij −η(∂uj ∂xi+∂ui ∂xj−δij∇·u), and q= −κ∇T−μ∇n, where ηis the shear viscosity, κis the (thermal) heat conductivity, and μis an additional transport coefficient called diffusive heat conductivity that couples the heat flux with the density gradients and that is peculiar to inelastic collisions [2]. The unknown term G(n,T) describes the homogeneous evolution of the temperature. All the gradient contributions coming from the cooling and heating rate are neglected. The system of Eqs. (6)–(8) admits a homogeneous stationary solution, n(r,t)=ns,u(r,t)=0, and T(r,t)=Ts.The two constants nsand Tsmust be related by Eq. (4), which has to be equivalent to G[ns,Ts(ns)] =0. The objective now is to investigate whether this state is linearly stable. To do that, it is convenient to introduce the dimensionless deviations of the fields around the homogeneous stationary solution ρ≡n−ns ns, w≡m 2Tsu, and θ≡T−Ts Ts, as functions of the dimensionless timescale, s, and dimensionless space variable l≡nsσr.Let us also introduce the Fourier components of these functions through yk≡dle−ik·ly(l) and let us decompose wkinto its parallel, wk,|| ≡wk·k k, and transversal wk,⊥≡wk·ˆ k⊥ components (ˆ k⊥is a unit vector perpendicular to k). The evolution equation for the transversal component is ∂ ∂swk,⊥= −ηk2wk,⊥, that is decoupled from the rest of Fourier components, that verify ∂ ∂s⎛ ⎝ ρk wk,|| θk ⎞ ⎠=L⎛ ⎝ ρk wk,|| θk ⎞ ⎠,(9) with L=− ⎛ ⎝ 0ik 0 i 2kηk2i 2k 2β+μk2ik β+κk2 ⎞ ⎠,(10) where the dimensionless transport coefficients, η≡σ √2mTsηs, κ≡m 2Tsσκs, and μ≡m 2 nsσ T3/2 s μs, have been introduced (the subindex sin the bare transport coefficients indicates that they are evaluated in the homogeneous stationary state). To obtain Eq. (9), the needed quantities ∂G(ns,Ts) ∂Tsand ∂G(ns,Ts) ∂ns have been identified by taking into account the analysis of the homogeneous phase made previously. In effect, ∂G(ns,Ts) ∂Ts= nsσ2Ts mβdue to Eq. (5), while ∂G(ns,Ts) ∂nshave been calculated using that G(ns,Ts)=0, so that ∂G(ns,Ts) ∂ns=2Ts ns ∂G(ns,Ts) ∂Ts. The stability of the system depends on the properties of the matrix Lgiven by Eq. (10). The eigenvalues of Lare the three roots, {λi}3 i=1, of the following algebraic equation: λ3+[β+(η+κ)k2]λ2+(1 +βη+ηκk2)k2λ −[β+(μ−κ)k2]k2 2=0.(11) Here, it is seen that a mode can vanish for a finite wave number, k=kc≡β κ−μ. In fact, it can be seen by general arguments that this is the only unstable mode for small k [15]. This occurs because the compressibility in the stationary state is negative, i.e., dp(ns) dns<0, and then the mode goes as λ≃−1 2Ts dp(ns) dnskto first order in k. More specifically, in our case the unstable mode goes as λ≃1 √2k. This mode becomes stable for k>kcdue to heat conduction (note that kcdepends also on μ), so that the condition of negative compressibility for the instability is necessary but not sufficient. In Fig. 5, the real part of the eigenvalues of Lare plotted for =0.5 and α=0.9 as functions of k. Similar results are obtained for other values of and α. The unstable (positive) mode is clearly identified. For the transport coefficients, we have taken their 2D elastic values [23], consistent with the previous analysis. Hence, a finite system of area L×Lcan be stable or unstable depending of whether L<Lc≡2π nsσkcor L>Lc. Note the strong analogy with the shear instability of the free cooling case [24], although, as seen in Fig. 5, the structure of the modes is very different. The theory developed here can be generalized for moderated densities. At the Enskog level, position correlations can be taken into account in an effective way by multiplying all the terms excepting the wall contribution in Eqs. (1) and (2) by the pair correlation function at contact, χ(ns). In this case, performing a similar analysis, the expression for the critical wave number is kc=βχ(ns) κ−μ(1 +2nsdlog χ(ns) dns). We have performed MD simulations to check the validity of this expression. Concretely, we have taken α=0.9,=0.5 and 032903-4 UNDERSTANDING AN INSTABILITY IN VIBRATED … PHYSICAL REVIEW E 99, 032903 (2019) 0 0.1 0.2 0.3 0.4 0.5 k -0.2 -0.15 -0.1 -0.05 0 0.05 λ kc FIG. 5. Real part of the eigenvalues of Las functions of kfor =0.5andα=0.9. kand λare measured in the dimensionless units defined in the main text. two different densities, nsσ2=0.15 and nsσ2=0.3. For the first density and taking the approximate expression for the equilibrium χgiven in Ref. [23], the theoretical prediction is such that the system is supposed to be unstable above a number of particles Nc≡4π2 nsσ2k2 c∼1650. We have performed MD simulations finding that for N=1300 and 1400 the system is stable, for N=1700 and 1800 it is unstable, while it is hard to say anything reliable for 1500 ⩽N⩽1600 as the system fluctuates between the homogeneous and the inhomogeneous phases. For the second density (nsσ2=0.3), it is Nc∼620, and for N=500 and 600 the steady state was observed to be stable, while for N=800 and 900 the system developed the instability. For N=700 no conclusion could be reached from the simulations. Hence, we can say that the theory predicts well the stability and instability of the system, as compared with MD results. V. CONCLUSIONS To sum up, we have introduced a simple model that allows us to explain from a microscopic point of view the origin of an instability appearing in some Q2D granular systems. Essentially, the idea is that there is one mechanism that destabilizes the system in the steady state (negative compressibility) and another one that stabilizes it (dissipation due to heat conduction). For small gradients (large systems), the first mechanism dominates the second one, making the system unstable, while for large gradients (small systems), the perturbation decays faster and the system is linearly stable. The fact that the compressibility is negative follows from the particular way in which energy is injected in the vertical direction and transferred to the horizontal direction through collisions. Hence, the instability in the model is more a question of “size” that of density and it is expected that a critical length can also be identified in the experiments. Of course, the instability can be tuned by many other aspects such as friction with the walls, inelasticity of particle-wall collisions, or gravity, to mention a few, but in our opinion the essential ingredients have been identified. On the other hand, the coexistence between the dense and dilute phases has not been treated as the formalism is no longer valid for very high densities. Nevertheless, it seems that the simplicity of the model would allow us to make progress in this direction. Also, the model developed here cannot be applied, in principle, to explain the high-density transition (solid-dense fluid coexistence) observed in some experiments [10,12]. ACKNOWLEDGMENTS This research was supported by the Ministerio de Economía, Industria y Competitividad (Spain) through Grant No. FIS2017-87117-P (partially financed by FEDER funds). [1] I. Goldhirsch, Rapid granular flows, Annu. Rev. Fluid Mech. 35,267 (2003). [2] J. J. Brey, J. W. Dufty, C. S. Kim, and A. Santos, Hydrodynamics for a granular flow at low density, Phys. Rev. E 58,4638 (1998). [3] J. J. Brey and D. Cubero, in Granular Gases,editedbyT. Pöschel and S. Luding (Springer-Verlag, Berlin, 2001). [4] J. S. Olafsen and J. S. Urbach, Clustering, Order, and Collapse in a Driven Granular Monolayer, Phys. Rev. Lett. 81,4369 (1998). [5] A. Prevost, D. A. Egolf, and J. S. Urbach, Forcing and Velocity Correlations in a Vibrated Granular Monolayer, Phys.Rev.Lett. 89,084301 (2002). [6] K. Roeller, J. P. D. Clewett, R. M. Bowley, S. Herminghaus, and M. R. Swift, Liquid-Gas Phase Separation in Confined Vibrated Dry Granular Matter, Phys.Rev.Lett.107,048002 (2011). [7] J. P. D. Clewett, K. Roeller, R. M. Bowley, S. Herminghaus, and M. R. Swift, Emergent Surface Tension in Vibrated, Noncohesive Granular Media, Phys. Rev. Lett. 109,228002 (2012). [8] P. Melby, F. Vega Reyes, A. Prevost, R. Robertson, P. Kumar, D. A. Egolf, and J. S. Urbach, The dynamics of thin vibrated granular layers, J. Phys.: Condens. Matter 17,S2689 (2005). [9] J. S. Olafsen and J. S. Urbach, Two-Dimensional Melting Far from Equilibrium in a Granular Monolayer, Phys.Rev.Lett.95, 098002 (2005). [10] G. Castillo, N. Mújica, and R. Soto, Fluctuations and Criticality of a Granular Solid-Liquid-Like Phase Transition, Phys. Rev. Lett. 109,095701 (2012). [11] M. Guzmán and R. Soto, Critical phenomena in quasi-two- dimensional vibrated granular systems, Phys. Rev. E 97,012907 (2018). [12] M. G. Clerc, P. Cordero, J. Dunstan, K. Huff, N. Mújica, D. Risso, and G. Varas, Liquid-solid-like transition in quasione-dimensional driven granular media, Nat. Phys. 4,249 (2008). [13] E. Khain and I. S. Aranson, Hydrodynamics of a vibrated granular monolayer, Phys.Rev.E84,031308 (2011). [14] D. R. M. Williams and F. C. MacKintosh, Driven granular media in one dimension: Correlations and equation of state, Phys. Rev. E 54,R9 (1996). 032903-5 MAYNAR, GARCÍA DE SORIA, AND BREY PHYSICAL REVIEW E 99, 032903 (2019) [15] R. Brito, D. Risso, and R. Soto, Hydrodynamic modes in a confined granular fluid, Phys. Rev. E 87,022209 (2013). [16] M. I. García de Soria, P. Maynar, and E. Trizac, Linear hydrodynamics for driven granular gases, Phys.Rev.E87,022201 (2013). [17] V. Garzó, M. G. Chamorro, and F. Vega Reyes, Transport properties for driven granular fluids in situations close to homogeneous steady states, Phys.Rev.E87,032201 (2013). [18] J. J. Brey, V. Buzón, M. I. García de Soria, and P. Maynar, Stability analysis of the homogeneous hydrodynamics of a model for a confined granular gas, Phys.Rev.E93,062907 (2016). [19] D. Risso, R. Soto, and M. Guzmán, Effective two-dimensional model for granular matter with phase separation, Phys.Rev.E 98,022901 (2018). [20] J. J. Brey, M. I. García de Soria, and P. Maynar, Boltzmann kinetic equation for a strongly confined gas of hard spheres, Phys.Rev.E96,042117 (2017). [21] J. J. Brey and M. J. Ruiz-Montero, Vibrated granular gas confined by a piston, Phys.Rev.E79,031305 (2009). [22] J. J. Brey, M. I. García de Soria, and P. Maynar (unpublished). [23] D. Gass, Enskog theory for a rigid disk fluid, J. Chem. Phys. 54, 1898 (1971). [24] I. Goldhirsch and G. Zanetti, Clustering Instability in Dissipative Gases, Phys.Rev.Lett.70,1619 (1993). 032903-6