Transversal inhomogeneities in dilute vibrofluidized granular fluids

Transversal inhomogeneities in dilute vibrofluidized granular fluids J. Javier Brey, M. J. Ruiz-Montero, F. Moreno, and R. Garcı ´a-Rojo Fı ´sica Teo ´rica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain 共Received 23 January 2002; published 21 June 2002兲 The spontaneous symmetry breaking taking place in the direction perpendicular to the energy flux in a dilute vibrofluidized granular system is investigated, using both a hydrodynamic description and simulation methods. The latter include molecular dynamics and direct Monte Carlo simulation of the Boltzmann equation. A marginal stability analysis of the hydrodynamic equations, carried out in the WKB approximation, is shown to be in good agreement with the simulation results. The shape of the hydrodynamic profiles beyond the bifurcation is discussed. DOI: 10.1103/PhysRevE.65.061302 PACS number共s兲: 45.70.Mg, 81.05.Rm, 47.20.⫺k I. INTRODUCTION Granular materials are assemblies of macroscopic particles dissipating their energy through inelastic collisions 关1兴. They exhibit a very rich phenomenology that is only partially understood. One of the most peculiar behaviors of granular systems, which has attracted a lot of attention in recent years, is their tendency to spontaneously develop strong spatial inhomogeneities. In many different situations, the density shows a sharp profile that is not induced by the boundary conditions. This phenomenon is often referred to as a clustering effect 关2兴, since high density regions coexist in the system with regions where the density is very low. In vibrated granular systems, clustering effects show up in many cases as a spontaneous symmetry breaking in the direction parallel to the vibrating wall. Consider a gas enclosed in a box that is being supplied energy through a vibrating wall located at x⫽0. There are no other external forces acting on the system. The box is divided into two equal compartments by a wall along the xaxis starting at a certain distance from the vibrating wall. At sufficiently low average density, the hydrodynamic fields are symmetric on both sides of the partition, but above a critical average density, which depends on the value of the restitution coefficient, an asymmetry in the number of particles at each side of the container occurs 关3兴. This asymmetry has been shown to be associated with a bifurcation of the solution of the hydrodynamic equations describing the state of the system. A similar symmetry breaking has been observed in a system in the presence of a gravitational force acting in the x direction. In this case, the system is unbounded for x⬎0, and the partition has a hole at a certain height. Again, an asymmetry in the number of particles in the two compartments develops if a control parameter, dependent on the amplitude of the vibration and the degree of inelasticity, is larger than a critical value 关4,5兴. Symmetry breaking in the direction parallel to the vibrating wall has also been observed in systems without any partition of the container. Sunthar and Kumaran 关6兴have reported molecular dynamics simulation results showing the presence of convection rolls and phase separation into coexisting dense and dilute regions in a granular system in the presence of gravity. The phase separation takes place on the surface of the vibrating wall and, as already indicated, in the direction perpendicular to the energy flux. No theoretical explanation for this phenomenon is provided in Ref. 关6兴, although the effect of the different parameters controlling the behavior of the system is discussed in detail, on the basis of the simulation results. For a two-dimensional closed system in the absence of gravity, a transversal continuous spontaneous symmetry breaking has also been predicted 关7兴. Not at all surprisingly, the gradients are now sharper next to the elastic wall, opposite the energy source. This is consistent with the positions that the holes must have in systems with a separating wall in order to observe symmetry breaking with and without a gravitational field acting on the system. The work by Livne et al. 关7兴is restricted to the nearly elastic limit and it is based on a numerical marginal stability analysis of the hydrodynamic equations. The predictions of this analysis are compared with numerical solutions of the hydrodynamic equations with the appropriate boundary conditions. In this paper, the bifurcation predicted in Ref. 关7兴will be considered again. There are several reasons for that. First, hydrodynamic equations derived from the Boltzmann equation for smooth inelastic hard disks and valid, in principle, for arbitrary inelasticity will be used, thus somewhat extending the previous results. However, it must be pointed out that, in steady states of granular systems such as the one considered here, there is a coupling between gradients and inelasticity. As a consequence, for these states small gradients imply in practice also small inelasticity. Secondly, instead of a numerical analysis of the stability of the solutions of the hydrodynamic equations, an analytical study, based on the WKB approximation, will be presented here. One of the main advantages of this approach is that the dimensionless control parameter governing the bifurcation phenomenon is clearly identified. A third motivation for the present work is to report simulation results, both from molecular dynamics and also from Monte Carlo simulation of the Boltzmann equation, showing the existence of the predicted transition. Since these simulation techniques do not contain any externally introduced hydrodynamic concept, they provide a direct proof of the existence of continuous symmetry breaking, and a test of the theoretical predictions. Attention will also be paid to the form of the hydrodynamic profiles beyond the bifurcation. This leads to a deeper understanding of the development of the instability. In any case, it is clear that the work in Ref. 关7兴opened the way to more systematic investiPHYSICAL REVIEW E, VOLUME 65, 061302 1063-651X/2002/65共6兲/061302共10兲/$20.00 ©2002 The American Physical Society65 061302-1 gations of this instability, such as the one in this paper. The plan of the paper is as follows. In Sec. II, the hydrodynamic description of the one-dimensional state of a low density vibrofluidized granular gas will be briefly summarized. The analytical expressions for the hydrodynamic profiles are given. This state is the starting point for the marginal stability analysis developed in Sec. III. Linearization of the hydrodynamic equations around the one-dimensional state leads to a second order linear differential equation. The WKB solution of the closed problem posed by this equation and the corresponding boundary conditions is built up. This requires consideration of three different cases, depending on the values of the parameters characterizing the system. From the WKB solution, the marginal stability curve follows easily. Simulation results are presented and compared with the theoretical predictions in Sec. IV, where an order parameter characterizing the transition is defined. A good agreement is found. The last section contains some final remarks, as well as a comparison of the results derived here with those in Ref. 关7兴in the common range of applicability, namely, nearly elastic collisions and very low density. II. BASIC EQUATIONS AND THE REFERENCE STATE For a steady state without macroscopic flows, the balance equations for a two-dimensional gas of smooth inelastic hard disks of mass mand diameter ␴ have the form “•P⫽0, 共1兲 “•q⫹nT ␨ ⫽0, 共2兲 where nand Tare the number density and granular temperature 共with Boltzmann’s constant set equal to unity兲, respectively. In the low density limit, for a gas described by the inelastic Boltzmann equation, and to lowest order in the gradients 共Navier-Stokes order兲, the pressure tensor Pand heat flux q, for the steady state under consideration, are given by 关8,9兴 P⫽pI,共3兲 q⫽⫺ ␬ “T⫺ ␮ “n,共4兲 Ibeing the unit tensor, p⫽nT the hydrodynamic pressure, ␬ the heat conductivity, and ␮ a transport coefficient that has no analog in the elastic case. These transport coefficients are proportional to the elastic heat conductivity ␬ 0(T), ␬ ⫽ ␬ *共 ␣ 兲 ␬ 0共T兲, ␮ ⫽ ␮ *共 ␣ 兲T ␬ 0共T兲 n,共5兲 ␬ 0共T兲⫽2 ␴ 冉 T ␲ m 冊 1/2 .共6兲 Finally, the cooling rate in the same approximation is related to the elastic shear viscosity ␩ 0by ␨ ⫽ ␨ *共 ␣ 兲p ␩ 0,共7兲 ␩ 0共T兲⫽1 2 ␴ 冉 mT ␲ 冊 1/2 .共8兲 The functions ␬ *, ␮ *, and ␨ *depend only on the constant coefficient of normal restitution ␣ characterizing the inelasticity of collisions. Their explicit expressions are given in Refs. 关9兴and 关10兴, and will not be reproduced here. By using the above expressions we get from Eqs. 共1兲and 共2兲 ⳵ p ⳵ x⫽ ⳵ p ⳵ y⫽0, 共9兲 共 ␬ *⫺ ␮ *兲 冋 ⳵ ⳵ x 冉 ␬ 0 ⳵ T ⳵ x 冊 ⫹ ⳵ ⳵ y 冉 ␬ 0 ⳵ T ⳵ y 冊 册 ⫺ ␨ *p2 ␩ 0 ⫽0. 共10兲 Next, we specify the boundary conditions. We consider that the system has Nparticles enclosed in a rectangular box with dimensions Lxand Ly. The wall located at x⫽0 is vibrating and, therefore, supplies energy to the system. This energy is needed in order to keep and sustain a fluidized steady state. The other three walls, located at x⫽Lx,y⫽0, and y⫽Ly, respectively, are at rest. For the sake of simplicity, collisions of particles with all four walls are assumed to be elastic. Then the mathematical boundary conditions to be imposed are 冉 ⳵ T ⳵ x 冊 x⫽Lx ⫽ 冉 ⳵ T ⳵ y 冊 y⫽0 ⫽ 冉 ⳵ T ⳵ y 冊 y⫽Ly ⫽0, 共11兲 ⫺关 ␬ *共 ␣ 兲⫺ ␮ *共 ␣ 兲兴 冋 ␬ 0共T兲 ⳵ T ⳵ x 册 x⫽0 ⫽Q.共12兲 Equations 共11兲express that the heat flux must vanish at the immobile walls, while Eq. 共12兲is the energy balance at the vibrating wall. The quantity Qis the rate of energy input through this wall per unit of length. Its calculation in terms of the parameters defining the motion of the wall has been addressed in several works. Here we will consider the simplest possibility, namely, that the wall moves in a sawtooth manner with velocity vb关10–12兴. This is a good approximation to more realistic motions, as long as the characteristic frequency of vibration of the wall is much larger than the collision rate of the gas molecules in its vicinity. In addition, it will be assumed that the amplitude of the vibration is much smaller than the mean free path of the particles of the gas next to it, so collective motions in the system are not being generated. Under the above conditions, it is 关11兴 Q⫽pvb.共13兲 The closed mathematical problem defined by Eqs. 共9兲–共13兲 admits a y-independent solution T⫽TR(x) that has been discussed in detail in Ref. 关10兴. The existence of this onedimensional solution was previously noticed by Grossman et al. 关13兴. In the following, we will refer to this solution as the reference state and its properties will be characterized with a subindex R. The reference temperature profile is BREY, RUIZ-MONTERO, MORENO, AND GARCI ´A-ROJO PHYSICAL REVIEW E 65 061302 061302-2 TR共 ␰ x兲⫽T0,R 冋 cosh共 ␰ x *⫺ ␰ x兲 cosh ␰ x * 册 2 ,共14兲 where ␰ xis a dimensionless scaled length defined by ␰ x⫽ 冑 a共 ␣ 兲 冕 0 xdx ␭R共x兲,共15兲 with ␭R(x) being the local mean free path, ␭R共x兲⫽1 2 冑 2 ␴ nR共x兲,共16兲 and a共 ␣ 兲⫽ ␲␨ * 16共 ␬ *⫺ ␮ *兲.共17兲 Moreover, ␰ x *is the value of ␰ xfor x⫽Lx, i.e., ␰ x *⫽2 冑 2a共 ␣ 兲 ␴ N Ly.共18兲 Therefore, ␰ x *is proportional to the number of monolayers of particles perpendicular to the xaxis at rest, ␴ N/Ly. Finally, the uniform pressure of the reference state is pR⫽T0,R共2 ␰ x *⫹sinh2 ␰ x *兲 8 冑 2a共 ␣ 兲 ␴ Lxcosh2 ␰ x *,共19兲 and the temperature T0,Rof the gas next to the vibrating wall is given by T0,R⫽ 冉 e tanh ␰ x * 冊 2 ,e共 ␣ 兲⫽ 冉 2am ␲ 冊 1/2 vb ␨ *.共20兲 The simplicity of the above results is a consequence of the limiting kind of motion of the vibrating wall considered, whose only effect is to transfer energy to the grains, without inducing any periodic motion in the system. In some previous studies 关13,14兴, a ‘‘thermal’’ wall, instead of a vibrating one, was considered at x⫽0. By definition, particles which collide with a thermal wall leave it with a velocity distribution corresponding to the temperature of the wall. Although this kind of wall is far from reality for granular systems, its consideration might be useful for comparison purposes. The only change to be made in the above discussion in order to apply it to a system with a thermal wall is to replace the boundary condition given in Eq. 共12兲by the requirement that the temperature T0,Rhas the value determined by the wall. As a consequence, the expressions for the hydrodynamic profiles remain the same, while Eq. 共20兲does not apply in this case. III. MARGINAL STABILITY ANALYSIS Our aim now is to investigate whether the system described in the preceding section exhibits another steady state, in addition to the reference one. Thus, we introduce a perturbation ␦ T(x,y)by T共x,y兲⫽TR共x兲⫹ ␦ T共x,y兲,共21兲 with ␦ T(x,y)ⰆTR(x), and search for a solution of this form to Eqs. 共9兲and 共10兲with the boundary conditions 共11兲and 共12兲. Substitution of Eq. 共21兲into Eq. 共10兲and linearization in ␦ Tyields 冋 ⳵ 2 ␦ ␰ x 2⫹2⫺1 2 冉 ⳵ lnTR ⳵␰ x 冊 2 ⫹ 冉 ␭R ␭ ¯ 冊 2 ⳵ 2 ⳵␰ y 2 册 ␦ T ⫽4TR pR ␦ p,共22兲 where ␦ phas been defined by p⫽pR⫹ ␦ p, ␰ xis the dimensionless scale introduced in Eq. 共15兲, and ␰ y⫽ 冑 a共 ␣ 兲y ␭ ¯ .共23兲 Here ␭ ¯ is the average mean free path, ␭ ¯ ⫽1 2 冑 2 ␴ n ¯ ,共24兲 n ¯ ⫽N/LxLy. Next, we consider factorized solutions of Eq. 共22兲, ␦ T共 ␰ x, ␰ y兲⫽ ␾ 共 ␰ x兲 ␸ 共 ␰ y兲.共25兲 The conditions at the boundaries to be satisfied by the functions ␾ and ␸ are 冉 ⳵ ␾ ⳵␰ x 冊 ␰ x⫽ ␰ x * ⫽ 冉 ⳵␸ ⳵␰ y 冊 ␰ y⫽0 ⫽ 冉 ⳵␸ ⳵␰ y 冊 ␰ y⫽ ␰ y * ⫽0, 共26兲 1 2 冉 ⳵ lnTR ⳵␰ x 冊 ␰ x⫽0 ␾ 共0兲⫹ 冉 ⳵ ␾ ⳵␰ x 冊 ␰ x⫽0 ⫽0, 共27兲 where ␰ y *⫽2 冑 2a共 ␣ 兲 ␴ N Lx.共28兲 Moreover, Eqs. 共11兲and 共12兲also imply that it must be ␦ p ⫽0, i.e., the pressure is not changed by the small perturbation. When Eq. 共25兲is substituted into Eq. 共22兲, use of separation of variables leads to the equations 1 ␸ 共 ␰ y兲 d2 ␸ 共 ␰ y兲 d ␰ y 2⫽⫺k2,共29兲 冋 d2 d ␰ x 2⫹2⫺1 2 冉 dlnTR d ␰ x 冊 2 ⫺ 冉 kn ¯ TR pR 冊 2 册 ␾ 共 ␰ x兲⫽0, 共30兲 TRANSVERSAL INHOMOGENEITIES IN DILUTE... PHYSICAL REVIEW E 65 061302 061302-3 where kis the constant of separation. The solution of Eq. 共29兲satisfying the corresponding boundary conditions in Eq. 共26兲is ␸ ⫽Acosk ␰ y,共31兲 with Aan arbitrary constant. Moreover, the values of kare restricted to k⫽ ␲ q/ ␰ y *,qbeing an integer. It is important to realize that the low density limit does not imply either ␰ x * Ⰶ1or ␰ y *Ⰶ1. In fact, both quantities can be large in a very dilute system. An additional restriction to be required to ␦ T(x,y) is that the total number of particles in the system, N, is conserved. It is easily seen that this condition is equivalent to 冕 0 ␰ x *d ␰ x 冕 0 ␰ y *d ␰ y ␦ T共 ␰ x, ␰ y兲 TR共 ␰ x兲⫽0. 共32兲 Equation 共25兲with ␸ ( ␰ y) given by Eq. 共31兲guarantees that this equality is satisfied. An equation having a structure similar to Eq. 共30兲was obtained in Ref. 关7兴for a dense system in the limit of nearly elastic collisions, by employing approximate constitutive relations introduced by Grossman et al. 关13兴. While the equation in 关7兴was solved using numerical techniques, here we will use a Wentzel-Kramers-Brillouin 共WKB兲approximation 关15兴to investigate the possible solutions 关16兴. First we use Eqs. 共14兲and 共19兲to rewrite Eq. 共30兲as d2 ␾ 共 ␰ x兲 d ␰ x 2⫹f共 ␰ x兲 ␾ 共 ␰ x兲⫽0, 共33兲 with f共 ␰ x兲⫽2 cosh2共 ␰ x *⫺ ␰ x兲 ⫺16k2 ␰ x *2cosh4共 ␰ x *⫺ ␰ x兲 共2 ␰ x *⫹sinh2 ␰ x *兲2, 共34兲 which is a monotonically increasing function of ␰ xin the whole interval 0⭐ ␰ x⭐ ␰ x *. To construct the WKB exponential approximation of Eq. 共33兲it is necessary to consider three different ranges of parameters, which will be discussed separately in the following. (a) The function f( ␰ x)is positive everywhere in the system. This is equivalent to f(0)⬎0or k⭐2 ␰ x *⫹sinh2 ␰ x * 2 冑 2 ␰ x *cosh3 ␰ x *.共35兲 Then, the WKB solution is oscillatory, ␾ 共 ␰ x兲⫽c1 冑 g共 ␰ x兲exp 冋 i 冕 0 ␰ xd ␰ x ⬘g共 ␰ x ⬘兲 册 ⫹c2 冑 g共 ␰ x兲exp 冋 ⫺i 冕 0 ␰ xd ␰ x ⬘g共 ␰ x ⬘兲 册 ,共36兲 where c1and c2are constants, and g共 ␰ x兲⫽ 冑 f共 ␰ x兲.共37兲 Imposing the boundary condition at ␰ x⫽ ␰ x *leads to ␾ 共 ␰ x兲⫽c 冑 g共 ␰ x兲cos 冕 ␰ x ␰ x *d ␰ x ⬘g共 ␰ x ⬘兲,共38兲 with canother constant. When the boundary condition at ␰ x ⫽0, Eq. 共27兲, is also required, the consistency condition tan 冕 0 ␰ x *d ␰ xg共 ␰ x兲⫽3tanh ␰ x * g3共0兲cosh2 ␰ x *共39兲 follows. The above equation determines the possible values of the parameters for which a WKB solution of the differential equation 共33兲exists in the region under study. If a thermal wall is considered at x⫽0, the boundary condition 共27兲must be replaced by ␾ (0)⫽0, and instead of Eq. 共39兲we find cos 冕 0 ␰ x *d ␰ xg共 ␰ x兲⫽0, 共40兲 i.e., 冕 0 ␰ x *d ␰ xg共 ␰ x兲⫽共2q⫹1兲 ␲ 2,共41兲 where qis an arbitrary integer. (b) The function f( ␰ x)is negative everywhere in the system. This is the case if f( ␰ x *)⭐0. Therefore, the region of parameters being considered is defined by k⭓k1⫽2 ␰ x *⫹sinh2 ␰ x * 2 冑 2 ␰ x *.共42兲 The exponential WKB approximation in this case reads ␾ 共 ␰ x兲⫽b1 冑 h共 ␰ x兲exp 冋 冕 0 ␰ xd ␰ x ⬘h共 ␰ x ⬘兲 册 ⫹b2 冑 h共 ␰ x兲exp 冋 ⫺ 冕 0 ␰ xd ␰ x ⬘h共 ␰ x ⬘兲 册 ,共43兲 with h共 ␰ x兲⫽ 冑 ⫺f共 ␰ x兲,共44兲 and b1and b2arbitrary constants. Imposing the boundary conditions 共26兲and 共27兲leads to the relationship h共0兲3tanh 冕 0 ␰ x *d ␰ xh共 ␰ x兲⫽3tanh ␰ x * cosh2 ␰ x *.共45兲 For a thermal wall at x⫽0, the above equation is substituted by BREY, RUIZ-MONTERO, MORENO, AND GARCI ´A-ROJO PHYSICAL REVIEW E 65 061302 061302-4 cosh 冕 0 ␰ x *d ␰ xh共 ␰ x兲⫽0, 共46兲 and, therefore, there is no WKB solution in this range of parameters. (c) The function f( ␰ x)changes sign in the interval 0 ⬍ ␰ x⬍ ␰ x *. It must be f(0)⬍0 and f( ␰ x *)⬎0. This corresponds to the kinterval 2 ␰ x *⫹sinh2 ␰ x * 2 冑 2 ␰ x *cosh3 ␰ x *⭐k⭐2 ␰ x *⫹sinh2 ␰ x * 2 冑 2 ␰ x *.共47兲 Since f( ␰ x) exhibits in this case a zero in the integration range, we have to consider separately the regions ␰ x⬍aand ␰ x⬎a, where ais the turning point, i.e., f(a)⫽0. In the former region, the WKB solution is given by an expression of the form 共43兲, while in the latter it has the form 共36兲.A global solution is constructed by matching both WKB approximations through the connection formulas expressing the connection between the oscillatory and exponential behaviors. Since it is easily seen that the turning point is a simple 共first-order兲zero, a standard application of the theory suffices 关15兴. Then, imposing the boundary conditions gives the following equation to be verified by the solutions in this range: tan 冉 ␪ ⫺ ␲ 4 冊 3tanh ␰ x *⫺h3共0兲cosh2 ␰ x * 3tanh ␰ x *⫹h3共0兲cosh2 ␰ x * ⫽1 2exp 冋 ⫺2 冕 0 ad ␰ xh共 ␰ x兲 册 ,共48兲 where ␪ ⫽ 冕 a ␰ x *d ␰ xg共 ␰ x兲.共49兲 The result for a thermal wall at x⫽0 is given in this case by tan 冉 ␪ ⫺ ␲ 4 冊 ⫽1 2exp 冋 ⫺2 冕 0 ad ␰ xh共 ␰ x兲 册 .共50兲 Once we have derived the equations determining the solutions to the boundary value problem, the strategy to build up the marginal stability curve is as follows. Given a value of ␰ x *, the first question is to see whether there is a bifurcation from the one-dimensional solution, i.e., whether Eq. 共30兲has a solution. In case the answer to this question is affirmative, the smallest value of ␰ y *for which there is a solution gives the point of the marginal stability curve corresponding to that value of ␰ x *. For larger values of ␰ y *, the one-dimensional solution is not stable. Suppose a solution of Eq. 共33兲exists for a given wave number k. This value is compatible with many 共infinite兲values of ␰ y *, the smallest one corresponding to the choice q⫽1 in the relationship between the possible values of kand ␰ y *关see below Eq. 共31兲兴. Moreover, the larger the value of k, the smaller the value of ␰ y *associated. The conclusion is that the stability curve is determined by the largest value of kfor which the problem has a solution. From a practical point of view, we have to start by looking for solutions belonging to case bdiscussed above. The values ␰ y *⫽ ␲ /kobtained from the solutions of Eq. 共45兲are the critical values, i.e., they define points on the marginal stability curve. Nevertheless, there is no solution for every value of ␰ x *belonging to the kinterval defining case b. More precisely, there is no solution at all for a thermal wall, while it is a simple matter to show that the existence of a solution in this region for a vibrating wall requires that 关h(1)共0兲兴3tanh 冕 0 ␰ x *d ␰ xh(1)共 ␰ x兲⭐3tanh ␰ x * cosh2 ␰ x *,共51兲 where h(1)共 ␰ x兲⫽lim k→k1 h共 ␰ x兲.共52兲 Equation 共51兲gives ␰ x *⯝0.555. Therefore, solutions for larger values of ␰ x *must be found, if they exist, in a different region of values of k, namely, in that considered in case c. With this procedure, it is an easy task to find the largest value of k共smallest value of ␰ y *) for which the eigenvalue problem defined by Eq. 共30兲has a solution for each value of ␰ x *,in the WKB approximation. The marginal stability curve for a vibrating wall obtained in this way is given in Fig. 1, where we have plotted the critical values ⌬cof the asymmetry parameter ⌬⬅Ly/Lx⫽ ␰ y */ ␰ x *, as a function of ␰ x *. The solid line is the WKB solution corresponding to the so-called case b, while the dashed line is from case c. Note that the solutions from both regions of parameters match rather smoothly at ␰ x *⯝0.555. Moreover, the critical asymmetry ⌬cis a FIG. 1. Marginal stability curve ⌬c( ␰ x *). The solid and dashed lines are the WKB predictions for a system driven by a vibrating wall, while the dotted line is for a system with an isothermal wall. Also shown are the results from DSMC 共filled symbols兲and MD 共open symbols兲simulations. For a given value of ␰ x *, the system exhibits transversal inhomogeneities in the steady state above the marginal curve. TRANSVERSAL INHOMOGENEITIES IN DILUTE... PHYSICAL REVIEW E 65 061302 061302-5 monotonically decreasing function of ␰ x *, growing very fast when ␰ x *tends to zero, as expected since the transversal inhomogeneities must disappear in the elastic limit. For a thermal wall, we have obtained the result that there is no solution in region b. Therefore, we must search for solutions belonging to the range cof parameters. Analysis of Eq. 共50兲leads to the result that the equation has no solution for ␰ x *ⱗ1.181. Moreover, Eq. 共41兲has no solution in that region either, implying that, in the WKB approximation, the transition to the transversally inhomogeneous state requires, in the case of systems driven by a thermal wall, that the inelasticity and the number of monolayers at rest be not small. The dotted line in Fig. 1 shows the marginal stability curve for a thermal wall. In the limit of large ␰ x *, the curve overlaps with the one for a vibrated system. IV. SIMULATION RESULTS To test the above theoretical results and also to investigate the nature of the predicted symmetry breaking, two more fundamental descriptions of the system, via the direct simulation Monte Carlo 共DSMC兲method and molecular dynamics 共MD兲simulation, respectively, have been considered. Although both are based on a dynamical simulation of the system, their nature is rather different. The DSMC method provides an algorithm to obtain numerical solutions of the Boltzmann equation for given initial and boundary conditions 关17兴. Therefore, it relies on the validity of a kinetic theory description, namely, the one given by the Boltzmann equation. On the other hand, no hydrodynamic approximations are introduced into the description, so that the validity of a hydrodynamic level of description, as provided for instance by the Navier-Stokes equations, is not taken for granted in this approach. The MD simulations follow the motion of the particles of the system as a sequence of free motions and binary collisions 关18兴, i.e., by direct application of Newton’s equations of motion. Therefore, MD provides the more basic description of the evolution of the system. In this context, it is important to stress that although the DSMC method also uses ‘‘particles’’ at a formal level, they are not real particles, but fictitious ones, which are employed to mimic the ideal dynamics described by the Boltzmann equation. In fact, the ideal nature of the particles in the DSMC method allows a very high numerical accuracy, since the number of particles used does not affect at all the physical state being simulated. In particular, the above number is not related to the actual density of the system. In the following, we will report simulation results obtained for the system studied in the previous sections, restricting ourselves to the vibrating wall. We have used the two methods DSMC and MD because they complement one another. For instance, comparison of the results obtained by both methods provides a test of the validity of the inelastic Boltzmann equation to describe a system under the physical conditions used in the MD simulations. Since the onedimensional state has been discussed in detail elsewhere 关10,13兴, attention will be focused here on the behavior of the system in the vicinity of the bifurcation. In the MD simulations, the density of course plays a relevant role. As the interest here is in the low density limit, small values of the surface fraction ␯ ⫽N ␲␴ 2/4LxLy, typically of the order of 10⫺2, have been used. For fixed given values of ␯ , ␣ , and Lx, a set of simulations have been run corresponding to different widths Lyof the system, starting from a small enough value. This means that each set of simulations corresponds to the same value of ␰ x *, as defined in Eq. 共28兲, while they differ in the asymmetry parameter ⌬. The simulations of the Boltzmann equation have been carried out with a similar systematic, the main difference being that the density plays no role in them. The other parameter needed to specify the simulations is the velocity of the vibrating wall vb. We have verified that, in agreement with the theory developed in the previous sections, the formation of transversal inhomogeneities is not altered by modifying this velocity, as long as it is large enough to fluidize the complete system. All the simulations started from a spatially homogeneous configuration, with a Gaussian velocity distribution, corresponding to an arbitrary temperature. The simulation is then followed until the system reaches a steady state, in which the several monitored properties of the system 共mean kinetic energy, density fluctuations, and hydrodynamic profiles兲become time independent. Once the system is in the steady state, all statistical averages of interest are accumulated. For the purpose here, the density and temperature profiles provide the relevant information. Let us describe what is observed in a set of simulations, namely, we are going to present DSMC data from systems with ␣ ⫽0.95 and Lx ⫽10␭ ¯ . This corresponds to ␰ x *⫽1.015. In Fig. 2 the steady two-dimensional density profile for Ly⫽20␭ ¯ (⌬⫽2) is shown in a three-dimensional plot. It is seen that no gradients in the ydirection are present, i.e., the system is in the onedimensional reference state. The transversal homogeneity apFIG. 2. Three-dimensional plot of the stationary density profile obtained by the DSMC method for a system with ⌬⫽2 and ␰ x * ⫽1.015. The density is normalized by its average value n ¯ , and the lengths with the average mean free path ␭ ¯ . The system does not exhibit gradients in the ydirection, i.e., it is in the reference state 共below the marginal stability curve兲. BREY, RUIZ-MONTERO, MORENO, AND GARCI ´A-ROJO PHYSICAL REVIEW E 65 061302 061302-6 pears even clearer in Fig. 3, where the density is plotted as a function of yfor several fixed values of x. When the width Lyis increased keeping all the other parameters fixed, a critical value shows up such that gradients in the y-direction spontaneously develop in the system for larger widths. An example of this is given in Fig. 4, where the density surface for the same parameters as in Fig. 2, except that now Ly⫽26.5␭ ¯ , is plotted. A gradient in the direction parallel to the vibrating wall is clearly identified, becoming more pronounced next to the opposite wall, as illustrated in Fig. 5. The density gradients in this case are relatively small, the maximum variation of the density in the ydirection being of the order of 5%. In fact, in all the simulations, both by DSMC and MD, it has been found that on increasing Lya continuous transition from the onedimensional state to a state with weak inhomogeneities in the transversal direction occurs. Moreover, the transversal density profile near the transition exhibits, as in the case of Fig. 5, a wavelength equal to twice the width of the system. In the language used in Sec. III, what is observed is a perturbation with q⫽1(k⫽ ␲ / ␰ y *), consistently with our theoretical analysis. If the width of the system is increased further, the transversal inhomogeneities grow very fast and, of course, the results from the linear marginal stability analysis do not apply. Simulations show that the density becomes very large in one of the corners of the system, away from the vibrating wall. The rapid increase of the gradients and the large value of the density in a localized region of the system lead us to conclude that in this regime a cluster of particles is formed 关7兴. Of course, for such a region of parameters the hydrodynamic profiles obtained from DSMC and MD simulations are quite different, as the former considers the particles as points while the latter assigns them a finite diameter, implying that the density is bounded by the close-packing value. Once the appearance of transversal inhomogeneities has been observed by visual inspection, it is desirable to have a ‘‘quantitative’’ criterion to establish whether the system is or is not transversally homogeneous. This is equivalent to identifying an order parameter to characterize the transition. Since there may be gradients in both directions, the identification of such a parameter is not at all a trivial task. The one we have chosen is defined as follows. First, we introduce the dimensionless quantity ␳ x(y)by ␳ x共y兲⫽n共x,y兲 Ly ⫺1 冕 0 Lydyn共x,y兲 ⫺1. 共53兲 If the system is transversally homogeneous, ␳ x(y) is independent of both xand y, and equal to zero. When transversal inhomogeneities are present, it depends of course on ybut, as can be guessed from Figs. 4 and 5, it also depends on x. FIG. 3. Density profiles as a function of yfor several fixed values of x, for the same system as in Fig. 1. The curves correspond, from bottom to top, to x⫽Lx/4, Lx/2, 3Lx/4, and Lx. FIG. 4. Three-dimensional plot of the density profile obtained by the DSMC method for the same system as in Fig. 2, with the only difference that now Ly⫽26.5␭ ¯ . The system is now above the marginal stability curve and transversal gradients are clearly observed. FIG. 5. Density profiles as a function of yfor several fixed values of xfor the same system as in Fig. 4. The curves correspond, from bottom to top, to x⫽Lx/4, Lx/2, 3Lx/4, and Lx. TRANSVERSAL INHOMOGENEITIES IN DILUTE... PHYSICAL REVIEW E 65 061302 061302-7 Nevertheless, the simulation data indicate that this dependence is rather weak, at least near the transition. As an example, in Fig. 6 the function ␳ x(y) has been plotted for the same system as in Figs. 4 and 5. The different lines correspond to four different values of x, equally separated, in the interval 关3Lx/4,Lx兴. This is the region where the transversal gradients are larger. From the figure it follows that the x dependence is essentially scaled out in the definition of ␳ x(y). It must be mentioned, however, that if the whole range of variation of xis considered, some dependence on x shows up. In any case, the departure from zero of the average value of ␳ x(y), ␳ ¯ (y), over a certain xinterval, next to the vibrating wall and not too large to avoid xdependence, provides a good criterion to distinguish transversally inhomogeneous systems from homogeneous ones. The results to be discussed in the following have been obtained using the interval 关3Lx/4,Lx兴. The above discussion, the theoretical analysis, and the numerical results, like those illustrated in Fig. 6, suggest that a good order parameter may be given by the absolute value of the first Fourier component 兩 f1 兩 of ␳ ¯ (y). In fact, we have verified that the absolute value of the Fourier transform of this quantity exhibits an abrupt maximum for the first component when transversal gradients begin to build up in the system. The behavior of 兩 f1 兩 as a function of the asymmetry ⌬in the vicinity of the transition point, is shown in Fig. 7 for the same parameters considered in Figs. 2–6. First of all, it must be noted that 兩 f1 兩 varies in a continuous way through the transition, although it grows very fast when one goes into the inhomogeneous region. This is the typical behavior of the order parameter of a nonequilibrium second order phase transition 关19兴. The continuous character of 兩 f1 兩 introduces some arbitrariness in the determination of the transition point ⌬c. We have made the choice, somewhat arbitrarily but consistently, that the transition takes place when 兩 f1 兩 becomes an order of magnitude larger than its typical value in the reference state, which is determined by the noise level. For the systems used in the DSMC method, this latter value is of the order of 10⫺3, so that a system has been considered as transversally inhomogeneous when 兩 f1 兩 ⬃10⫺2, which implies deviations from homogeneity of the order of 1%. This leads in the case of Fig. 6 to an estimation ⌬c⫽2.5⫾0.1 for the critical asymmetry. By changing the initial parameters of the system and repeating the above procedure, ⌬chas been computed for different values of ␰ x *. The results from the DSMC simulation are represented by the filled symbols in Fig. 1. The agreement between the theoretical predictions and the simulation data is rather good, although a systematic deviation appears, larger for smaller ␰ x *. When evaluating the comparison, it must be taken into account that the simulations only provide an upper bound for ⌬c. When the system is very close to the transition point, the time required to go from the transversally homogeneous state to the inhomogeneous one may be too large to observe the transition during the simulation time. In any case, it is fair to say that the hydrodynamic equations and the WKB approximation provide an accurate description of what is observed in the simulations. The results discussed up to this point were obtained from DSMC simulations. Just to illustrate how MD simulations lead to an equivalent scenario, we present next some results for a set of MD simulations with ␣ ⫽0.925, ␯ ⫽10⫺2, and Lx⫽100 ␴ . For these values, it is ␰ x *⫽0.281. In Fig. 8 a three-dimensional plot of the density profile is shown for ⌬ ⫽6.4. The system exhibits gradients in the transversal direction, increasing again as we move away from the vibrating wall. However, the density gradients are very small, the maximum deviation from homogeneity being of the order of 10%. That means that the system is close to the transition, probably in the region where a linear approximation around the steady state still provides an accurate description. FIG. 6. DSMC results for the dimensionless function ␳ x(y), defined by Eq. 共53兲, for the same system as in Figs. 4 and 5. The different curves correspond to equally spaced values of xin the interval 关3Lx/4,Lx兴. FIG. 7. DSMC results for dimensionless order parameter 兩 f1 兩 as a function of the asymmetry ⌬for a system with ␰ x *⫽1.015 共the dashed line has been included as a guide for the eye兲. A second order nonequilibrium bifurcation is clearly identified. BREY, RUIZ-MONTERO, MORENO, AND GARCI ´A-ROJO PHYSICAL REVIEW E 65 061302 061302-8 When the asymmetry is increased further, gradients in the perpendicular direction become sharper, and a state with a sharply peaked density is observed. This is illustrated in Fig. 9 for ⌬⫽7. In Fig. 10, the quantity ␳ ¯ is shown for different simulations belonging to the set we are considering, i.e., they differ only in the value of ⌬. The continuous transition from the reference state to the transversally asymmetric one is clearly observed, as well as the dramatic increase of the transversal gradients when the system gets well inside the unstable region. In conclusion, MD results are in full qualitative agreement with the DSMC ones. Even more, the critical values of the asymmetry parameters obtained from MD are in very good quantitative agreement with those following from DSMC calculations, as seen in Fig. 1, where they are represented by the open symbols. This provides a proof of the validity of the kinetic theory description as given by the Boltzmann equation for dilute inelastic gases. It is worth emphasizing that points in Fig. 1 were obtained from DSMC and MD simulations by changing the values of ␣ and ␴ N/Lyto sample different values of ␰ x *. The fact that ⌬cobtained in this way varies smoothly with ␰ x *supports the theoretical prediction that the dependence on the different parameters of ⌬coccurs through it. This has also been ratified by considering two sets of simulations having different values of both ␣ and ␴ N/Ly, but leading to the same value of ␰ x *. The same critical asymmetry parameter was found, supporting the scaling predicted by the theory. V. DISCUSSION AND CONCLUDING REMARKS In this paper, the spontaneous transversal symmetry breaking in a vibrated granular fluid predicted by Livne et al. 关7兴was further investigated for low density systems. It was shown that the transition takes places both in MD simulations and in systems described by the Boltzmann equation. Moreover, there is a reasonably good agreement between the theoretical predictions, following from a marginal stability analysis of the hydrodynamic description of the system, and the results from MD and also from the numerical solution of the Boltzmann equation by the DSMC method. This refers to the values of the critical asymmetry as a function of the control parameter, and also to the shape of the hydrodynamic profiles in the vicinity of the symmetry breaking. To characterize the transition, an order parameter quantifying the initial setup of transversal inhomogeneities has been introduced. In terms of this parameter, the transition presents the features of a second order nonequilibrium phase transition. In this context, it is worth mentioning that no subcritical bifurcations have been observed in the simulations. Moreover, for states well inside the instability curve, the density profile exhibits a characteristic nonlinear ␭/2 shape. On the other hand, in Ref. 关7兴, nonlinear twodimensional states inside the linear stability region were found at high densities from the numerical solutions of the hydrodynamic equation. Of course, there is no contradiction in this, since our analysis was restricted to low density gases. In Ref. 关7兴an analytical expression for the marginal stability curve in a certain limit is given. In our notation, the limit considered is ␰ x *Ⰶ1, and the expression reads ⌬⯝1.6/ ␰ x *2, where we have neglected subleading terms in the density. An FIG. 8. Three-dimensional plot of the stationary density profile obtained by MD simulation. The particles are disks of diameter ␴ , the area fraction is ␯ ⫽10⫺2, the coefficient of restitution ␣ ⫽0.925, Lx⫽100 ␴ , and ⌬⫽6.4. Small transversal gradients are present, and the system is in the vicinity of the bifurcation. FIG. 9. The same as in Fig. 8 but with ⌬⫽7. Quite large transversal gradients growing in the xdirection are identified. The system is above the marginal stability curve. FIG. 10. MD results for the dimensionless function ␳ ¯ (y), for several values of the asymmetry ⌬, as indicated in the figure. All the other parameters of the system are the same as in Figs. 8 and 9. The transition to a state inhomogeneous in the ydirection is clearly identified. TRANSVERSAL INHOMOGENEITIES IN DILUTE... PHYSICAL REVIEW E 65 061302 061302-9