Papers in Physics, vol. 3, art. 030004 (2011) Received: 11 May 2011, Accepted: 15 July 2011 Edited by: J-C. Géminard Reviewed by: B. Tighe, Instituut-Lorentz, Universiteit Leiden, Netherlands. Licence: Creative Commons Attribution 3.0 DOI: http://dx.doi.org/10.4279/PIP.030004 www.papersinphysics.org ISSN 1852-4249 Master curves for the stress tensor invariants in stationary states of static granular beds. Implications for the thermodynamic phase space Luis A. Pugnaloni,1∗ José Damas,2† Iker Zuriguel,2‡ Diego Maza2§ We prepare static granular beds under gravity in different stationary states by tapping the system with pulsed excitations of controlled amplitude and duration. The macroscopic state—defined by the ensemble of static configurations explored by the system tap after tap—for a given tap intensity and duration is studied in terms of volume, V , and force moment tensor, Σ. In a previous paper [Pugnaloni et al., Phys. Rev. E 82, 050301(R) (2010)], we reported evidence supporting that such macroscopic states cannot be fully described by using only V or Σ, apart from the number of particles N. In this work, we present an analysis of the fluctuations of these variables that indicates that V and Σ may be sufficient to define the macroscopic states. Moreover, we show that only one of the invariants of Σ is necessary, since each component of Σ falls onto a master curve when plotted as a function of Tr(Σ). This implies that these granular assemblies have a common shape for the stress tensor, even though it does not correspond to the hydrostatic type. Although most results are obtained by molecular dynamics simulations, we present supporting experimental results. I. Introduction The study of static granular systems is of funda- mental importance to the industry to improve the storage of such materials in bulk, as well as to op- timize the packaging design. However, far from yielding such benefits of practical interest, physi- cists have found a fascinating challenge on their way that has stuck them, to a large extent, in the ∗E-mail: luis@iflysib.unlp.edu.ar †E-mail: jdelacruz@alumni.unav.es ‡E-mail: iker@unav.es §E-mail: dmaza@unav.es 1 Instituto de F́ısica de Ĺıquidos y Sistemas Biológicos, (CONICET La Plata, UNLP), Calle 59 No 789, 1900 La Plata, Argentina. 2 Departamento de F́ısica y Matemática Aplicada, Fac- ultad de Ciencias, Universidad de Navarra, Pamplona, Spain. area of static granular matter. Such challenge is finding an appropriate description of the simplest state in which matter can be found: equilibrium [1]. At equilibrium, a sample will explore differ- ent microscopic configurations over time, in such a way that macroscopic averages over large periods will be well defined, with no aging. Moreover, if not only equilibrium but ergodicity is present, av- erages over a large number of replicas at a given time should give equivalent results to time aver- ages [2]. Finally, one expects that all macroscopic properties of such equilibrium states can be put in terms of a few independent macroscopic variables. Then, a thermodynamic description and, hopefully, a statistical mechanics approach can be attempted. Theoretical formalisms based on these assumptions have been used to analyze data from experiments and numerical models. However, some of the foun- dations are still supported by little evidence. 030004-1 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. The use of careful protocols to make a granular sample explore microscopic configurations within a seemingly equilibrium macroscopic state has given us the first standpoint [3–6]. In such protocols, the sample is subjected to external excitations of the form of pulses. Well defined, reproducible time av- erages are found after a transient if the same pulse shape and pulse intensity are applied. However, for low intensity pulses, a previous annealing might be in order since equilibrium is hard to reach —as in supercooled liquids. We will use the expressions equilibrium state, steady state or simply state to refer to the collection of all configurations gener- ated by using a given external excitation after any transient has disappeared. More than twenty years ago, Edwards and Oak- shott [7] put forward the idea that the number of grains N and the volume V are the basic state vari- ables that suffice to characterize a static sample of hard grains in equilibrium. The NV granular en- semble was then introduced as a collection of mi- crostates, where the sample is in mechanical equi- librium, compatible with N and V . However, newer theoretical works [8–13] suggest that the force mo- ment tensor, Σ, (Σ = V σ, where σ is the stress ten- sor) must be added to the set of extensive macro- scopic variables (i.e., an NV Σ ensemble) to ade- quately describe a packing of real grains. In the rest of this paper, we show experimen- tal and simulation evidence that the equilibrium states of static granular packings cannot be only described by V (or equivalently the packing frac- tion, φ, defined as the fraction of the space covered by the grains) nor by Σ. We do this by generat- ing states of equal V but different Σ and states of equal Σ but different φ. We also show that states of equal V may present different volume fluctua- tions. Moreover, we show that states of equal V and Σ display the same fluctuations of these vari- ables, suggesting that no other extensive parameter might be required to characterize the state (apart from N). Finally, but of major significance, we show that the shape of the force moment tensor is universal, in the sense that different states that present the same trace of the tensor actually have the same value in all the components of Σ. II. Simulation We use soft-particle 2D molecular dynamics [14,15]. Particle–particle interactions are controlled by the particle–particle overlap ξ = d − |rij| and the ve- locities ṙij, ωi and ωj. Here, rij represents the center-to-center vector between particles i and j, d is the particle diameter and ω is the particle an- gular velocity. These forces are introduced in the Newton’s translational and rotational equations of motion and then numerically integrated by a ve- locity Verlet algorithm [16]. The interaction of the particles with the flat surfaces of the container is calculated as the interaction with a disk of infinite radius. The contact interactions involve a normal force Fn and a tangential force Ft. Fn = knξ −γnvni,j (1) Ft = −min (µ|Fn|, |Fs|) · sign (ζ) (2) where Fs = −ksζ −γsvti,j (3) ζ (t) = ∫ t t0 vti,j (t ′) dt′ (4) vti,j = ṙij ·s + 1 2 d (ωi + ωj) (5) The first term in Eq. (1) corresponds to a restor- ing force proportional to the superposition ξ of the interacting disks and the stiffness constant kn. The second term accounts for the dissipation of energy during the contact and is proportional to the nor- mal component vni,j of the relative velocity ṙij of the disks. Equation (2) provides the magnitude of the force in the tangential direction. It implements the Coulomb’s criterion with an effective friction fol- lowing a rule that selects between static or dynamic friction. Notice that Eq. (2) implies that the max- imum static friction force |Fs| used corresponds to µ|Fn|, which effectively sets µdynamic = µstatic = µ. The static friction force Fs [see Eq. (3)] has an elastic term proportional to the relative shear dis- placement ζ and a dissipative term proportional to the tangential component vti,j of the relative veloc- ity. In Eq. (5), s is a unit vector normal to rij. The 030004-2 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. elastic and dissipative contributions are character- ized by ks and γs, respectively. The shear displace- ment ζ is calculated through Eq. (4) by integrating vti,j from the beginning of the contact (i.e., t = t0). The tangential interaction behaves like a damped spring which is formed whenever two grains come into contact and is removed when the contact fin- ishes [17]. The particular set of parameters used for the simulation is (unless otherwise stated): µ = 0.5, kn = 10 5(mg/d), γn = 300(m √ g/d), ks = 2 7 kn and γs = 200(m √ g/d). In some cases, we have varied µ and γn in order to control the friction and resti- tution coefficient. The integration time step is set to δ = 10−4 √ d/g. The confining box (13.39d-wide and infinitely high) contains N = 512 monosized disks. Units are reduced with the diameter of the disks, d, the disk mass, m, and the acceleration of gravity, g. Tapping is simulated by moving the confining box in the vertical direction following a half sine wave trajectory [A sin(2πνt)(1−Θ(2πνt−π))]. The excitation can be controlled through the amplitude, A, and the frequency, ν, of the sinusoidal trajectory. We implement a robust criterion based on the sta- bility of particle contacts to decide when the sys- tem has reached mechanical equilibrium [14] before a new tap is applied to the sample. Averages were taken over 100 taps in the steady state and over 20 independent simulations for each value of A and ν. The volume, V , of the system after each tap can be obtained from the packing fraction, φ, as V = Nπ(d/2)2/φ. We measure φ in a rectangu- lar window centered in the center of mass of the packing. The measuring region covers 90% of the height of the granular bed (which is of about 40d) and avoids the area close to the walls by 1.5d. We have observed that φ is sensitive to the chosen win- dow. However, none of the conclusions drawn in this paper are affected by this choice. The stress tensor, σ, is calculated from the particle–particle contact forces as σαβ = 1 V ∑ cij rαijf β ij, (6) where the sum runs over all contacts. The force moment tensor, Σ, is defined as Σ ≡ V σ. During the course of a tap Σ is non-symmetric, however, once mechanical equilibrium is reached in accordance with our criterion, Σ becomes symmet- ric within a very small error compared with the fluctuations of Σ. Although Σ may depend on depth, we have measured the force moment tensor by simply summing over particle–particle contacts in the entire system. The fluctuations of φ (∆φ) and Σ (∆Σ) are cal- culated as the standard deviation in the 100 taps obtained in each steady state. We average φ, Σ, and their fluctuations over 20 independent runs for each steady state and estimate error bars as the standard deviation over these 20 runs. III. Experimental method (a) (b) Figure 1: (a) Schematic diagram of the experimen- tal set-up. A: accelerometer, S: shaker, C: camera, O: oscilloscope, FG: function generator, PC: com- puter. (b) Example of an image of the packing. The region of measurement is indicated by a red square. The experimental set-up is sketched in Fig. 1(a). A quasi 2D Plexiglass cell (28 mm wide and 150 mm high) is filled with 900 alumina oxide beads of di- ameter d = 1±0.005 mm. The separation between the front and rear plates was made 15% larger than the bead diameter. The cell is tapped by an elec- 030004-3 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. tromagnetic shaker (Tiravib 52100) with a trend of half sine wave pulses separated five seconds be- tween them. The tapping amplitude was controlled adjusting the intensity, A, and frequency, ν, of the pulse and was measured by an accelerometer at- tached to the base of the cell. Averages are taken over 500 taps after equilibration. High resolution images (more than 10 MPix) are taken after each tap. To calculate the packing frac- tion, we only consider a rectangular zone at the center of the packing whose limits are at 4 mm from the borders [see Fig. 1(b)]. We determine the centroid of each particle by means of a numerical algorithm with subpixel resolution. Then, we cal- culate the packing fraction by assuming that the 2D projection of each bead corresponds to a disk of diameter d = 1 mm. We estimate the packing fraction with a resolution of ±0.001. As the sep- aration between plates is larger than the particle diameter, small overlaps between the spheres are present in the 2D projections of the images. There- fore, the calculated packing fraction in some dense configurations might result slightly higher than the hexagonal disk packing limit (π √ 3/6). IV. Tapping characterization and asymptotic equilibrium states It is often debated [18, 19] what is the appropri- ate parameter to characterize the external excita- tion used to drive a granular sample. Dijksman et al. [18] proposed a parameter related to the lift- off velocity of the granular bed. Ludewing et al. [19] presented an energy based parameter. Pug- naloni et al. [15, 20] suggested that the factor of expansion induced on the sample would be a suit- able measure [21]. The usual perspective to define a pulse parameter, �, is to achieve a collapse of the φ–� curves as different details of the pulse shape are changed (such as amplitude and duration). Pa- rameters defined in all these previous works fail to make such curves collapse for the data presented in this paper. The main reason for this is that our φ–� curves are non-monotonic (see for example Fig. 6) presenting a minimum whose depth depends on the details of the pulse shape. Therefore, a sim- ple rescale of the horizontal axis does not suffice to collapse the curves. Since we are interested in macroscopic states, the actual pulse used to drive the system is merely a control parameter but not a macroscopic variable that describes the state. Therefore, the external pulse does not need to be described with a simpli- fied quantity. The complete functional form of the pulse can be given instead. In our case, we use a sine pulse and both, the pulse amplitude and fre- quency, are needed to fully describe the excitation. We will employ the usual nondimensional peak ac- celeration Γ ≡ apeak/g = A(2πν)2/g (where g is the acceleration of gravity) and the frequency ν to precisely define the external excitation. A detailed study of the dynamical response of a granular bed to a pulse of controlled intensity and duration can be found in Damas et al. [22]. One major issue in studying equilibrium states is the evidence that one can have, indicating if the system is actually at equilibrium [5]. Since the defi- nition of equilibrium is circular [23], we can simply do our best to check if different properties of the system have well defined means (as well as higher order moments of the distributions) which should not depend on the history of the processes applied to the sample. 1 10 100 1000 Number of taps 0.82 0.84 0.86 0.88 0.9 φ Figure 2: Evolution of φ towards the steady state starting from a disordered configuration initially obtained by strong taps (experimental results). Re- sults for two tapping intensities are shown: Γ = 4.8 (blue), and Γ = 17 (green). In both experiments, the frequency of the pulse is ν = 30 Hz. We have generated configurations corresponding to a particular pulse (of a given shape and inten- sity) by repeatedly applying such pulse to the sys- tem and allowing enough time for any transient to fade. To prove that our samples are at equilib- 030004-4 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. rium, we approach a particular pulse through differ- ent paths —and starting from different initial con- ditions (ordered and disordered configurations)— and confirm that the steady states obtained present equivalent mean values and second moments of the distributions of the variables of interest. In Ref. [26], we showed that the steady states corresponding to excitations of high intensity are reached in a few taps, even if the initial configu- ration corresponds to a very ordered structure. In Fig. 2, we consider the equilibration process from a highly disordered initial configuration. For low tap- ping intensities (blue line in Fig. 2), the packing fraction evolves to the steady state in two stages. Just after switching the tapping amplitude to the reference value, the system rapidly evolves to values of φ close to the final steady state. Beyond this ini- tial convergence, a slower compaction phase takes the system to the final steady state. For high tap- ping intensities, the evolution to the steady state is very rapid. The steady state is reached after about a hundred taps [green line in Fig. 2)]. Therefore, we apply a sequence of at least 1000 taps in all our experiments before taking averages to warrant that the steady state has been reached. In our simula- tions, 400 taps of equilibration were enough. An interesting example of equilibration is pre- sented in Fig. 3. In this case, we present the values of φ and of Tr(Σ) during a very special sequence of pulses obtained in our simulations. The system is initially deposited from a dilute random config- uration with all particles in the air. First, we tap the system 200 times at Γ = 4.9, then 800 times at Γ = 61.5 and finally, 200 times at Γ = 4.9. In the whole run, we keep ν = 0.5 √ g/d. These two values of the pulse intensity have been chosen be- cause they are known to produce packings with the same mean φ in the steady state [26]. However, the mean Σ is clearly different. Notice, however, that the same values of φ, Σ and its fluctuations, ∆φ and ∆Σ, are observed for the steady state of low Γ obtained before and after the 800 pulses of high Γ. Unless otherwise stated, all the results we present in what follows correspond to steady states. We have tested this by obtaining the same states through two different preparation protocols consist- ing of: (i) application of a large number of identi- cal pulses starting from a disordered configuration, and (ii) application of a reduced number of iden- tical pulses after annealing the system from much 0 500 1000 Number of taps 0.82 0.84 0.86 0.88 0.9 0.92 φ ∆φ=0.01185 φ=0.85905 ∆φ=0.01192 φ=0.85959 ∆φ=0.01191 φ=0.85887 (a) 0 500 1000 Number of taps 20 30 40 50 60 70 80 90 100 T r( ) /N [ u n it s o f m g d ] Σ=50.97 ∆Σ=10.7 ∆Σ=6.9 Σ=37.68 ∆Σ=11.1 Σ=51.77 Σ (b) Figure 3: Evolution of φ (a) and Tr(Σ) (b) as the pulse intensity is suddenly changed between two values that produce steady states with the same φ but different Σ in our simulations. The initial 200 taps correspond to Γ = 4.9, the middle 800 pulses to Γ = 61.5 and the final 200 pulses, again, to Γ = 4.9. In all cases, ν = 0.5 √ g/d. The mean values and standard deviations in each section are indicated by arrows. higher pulse intensities. In a few cases, the results (mean values and/or fluctuations) from both pro- tocols did not match. This was an indication that a steady state, if existed, was not reached by one of the protocols (or by both protocols). This hap- pens especially for low intensity taps, which require longer equilibration times. Such cases have been removed from the analysis. V. Volume and volume fluctuations In Fig. 4(a), we plot φ in the steady state as a func- tion of Γ from our experiments for different pulse frequencies [24]. As we can see, there exists a min- 030004-5 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. imum φ at relatively high Γ. A similar experiment on a three-dimensional cell also yielded analogous results [26]. This behavior has also been reported for various models [15, 20, 25]. An explanation for this, based on the formation of arches, has been given in [15]. The position of the minimum in φ shifts to larger Γ if the frequency of the pulse is increased (i.e., if the pulse duration is reduced). 0 10 20 30 40 50 Γ 0.84 0.86 0.88 0.9 0.92 φ 30 45 60 ν (Hz) (a) 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 φ 0 50 100 F re q u e n c y c o u n t (b) Figure 4: (a) Experimental results for the steady state packing fraction, φ, as a function of the re- duced peak acceleration, Γ, for different frequen- cies, ν, of the tap pulse. (b) Histogram of the con- figurations visited by the system for ν = 30 Hz: Γ = 15 (green) and Γ = 28 (red). Although the increment in the packing fraction beyond the minimum is rather small, it is impor- tant to remark that this difference is not an artifact introduced by our experimental resolution. In Fig. 4(b), we show the histogram for the sequence of packings obtained for Γ ' 15 (the minimum pack- ing fraction for ν = 30 Hz) and for Γ = 28 (the largest excitation explored for ν = 30 Hz). Both steady states are statistically comparable; however, it is possible to distinguish different mean values. Since steady states of equal φ obtained at both 0 10 20 30 Γ 0.005 0.006 0.007 0.008 0.009 0.01 ∆φ Figure 5: The steady state volume fraction fluctua- tions, ∆φ, as a function of Γ from our experiments with ν = 30 Hz. The solid line is only a guide to the eyes. sides of the φ minimum are generated via very dif- ferent tap intensities, it is worth assessing if such states are, in fact, equivalent. This can be done by comparing the volume fluctuations of such states. A similar analysis was done in Ref. [27] for states generated with different pulse amplitude and dura- tion in liquid fluidized beds. The fluctuations of φ in the steady state, as mea- sured by the standard deviation ∆φ, are presented in Fig. 5. As we can see, a minimum in the fluctu- ations is just apparent. The position of such min- imum in the fluctuations coincides with the mini- mum in φ. Unfortunately, the resolution of φ in our experiments is around the size of the fluctuations. However, the results from the simulations are not limited in this respect and we address to those for a more reliable assessment of the fluctuations. In Fig. 6(a), we plot φ in the steady state as a function of Γ for our simulations. Although the fluctuations of φ are large [see Fig. 3], its mean value is well defined with a small confidence inter- val (see error bars). For low excitations, φ decreases as Γ is increased. However, beyond a certain value Γmin, the packing fraction grows. The same trend is observed if the tap frequency ν is changed. How- ever, the minimum is deeper for lower ν and its position Γmin shifts to larger values of Γ as ν is increased in coincidence with our experiments (see Fig. 4). Due to the change in the depth of the φ minimum in the φ–Γ curves, a simple rescaling of Γ is unable to collapse the data for different frequen- cies. However, rescaling φ and Γ with φmin and 030004-6 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. Γmin, respectively, yields very good collapse, even between simulated and experimental data [26]. In Fig. 6(b), the volume fraction fluctuations, ∆φ, are plotted as a function of Γ. As we can see, these fluctuations are non-monotonic, as suggested by our experiments (Fig. 5). Non-monotonic vol- ume fluctuations have also been reported in Ref. [6]. For ∆φ, we obtain a minimum and a maximum. We have also observed a maximum in ∆φ for values of Γ below the ones reported here (see for example Refs. [26] and [25]). However, we do not report such low values of Γ in this work and focus on tap- ping intensities that warrant the steady state with a modest number of pulses. 1 10 100 Γ 0.81 0.82 0.83 0.84 0.85 0.86 φ 0.125 0.250 0.500 0.750 1.000 ν (g/d) 1/2 (a) 1 10 100 Γ 0.006 0.008 0.01 0.012 0.014 0.016 ∆φ 0.125 0.250 0.500 0.750 1.000 ν (g/d) 1/2 (b) Figure 6: (a) Simulation results for the volume frac- tion, φ, as a function of the reduced peak accelera- tion, Γ, for different frequencies, ν, of the tap pulse. (b) The corresponding volume fraction fluctuation, ∆φ, as a function of Γ. The value of Γ, at which the fluctuations dis- play a minimum, coincides with Γmin, the value at which the minimum packing fraction, φmin, is ob- tained. The maximum coincides with the inflection point in the φ–Γ curve at higher Γ. Since one ex- pects to find few mechanically stable configurations compatible with a large volume (low φ), it seems reasonable that fluctuations reach a minimum if φ does so. Similarly, there should be few low-volume, mechanical stable configurations which implies that at high φ fluctuations should diminish. Hence, a maximum in ∆φ should be present at intermedi- ate packing fractions. This is more clearly seen in Fig. 7 where we plot the fluctuations as a function of the average value of φ. 0.81 0.82 0.83 0.84 0.85 0.86 φ 0.006 0.008 0.01 0.012 0.014 0.016 ∆φ 0.125 0.250 0.500 0.750 1.000 ν (g/d) 1/2 Figure 7: Density fluctuations as a function of φ for different frequencies of the tap pulse in the sim- ulations. Figure 7 presents two distinct branches: the lower branch corresponds to Γ > Γmin and the up- per branch to Γ < Γmin. For Γ > Γmin, the fluc- tuations corresponding to different tap durations collapse, suggesting that such equal-φ, equal-∆φ states might correspond to unique states (below, we will find that this is not the case). For Γ < Γmin, we obtain states of same φ as some states in the lower branch but presenting larger fluctuations. This is clear evidence that the equal-φ states of the upper and lower branch are indeed distinct and that other macroscopic variables must be used to distinguish one from the other. We have assessed a number of other structural descriptors (coordination number, bond order pa- rameter, radial distribution function). In all cases, equal–φ states from the upper and lower branch of the φ–∆φ curve (Fig. 7) present similar values of the structural descriptors with only subtle dis- crepancies. Although this indicates that the states are not equivalent, it also suggests that such de- scriptors are not good candidates to form a set of macroscopic variables, along with V , to uniquely 030004-7 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. identify a given steady state. In Ref. [26], we have assessed the force moment tensor Σ as a good candidate to complete V and N in describing a stationary state. This has been sug- gested by some theoretical speculations [12]. How- ever, some authors prefer to directly replace V by Σ. Below, we will show that both, V and Σ, are required for the equilibrium states generated in our simulations. Moreover, we will show that only one of the invariants of Σ is necessary (at least in our 2D systems) and that fluctuations of these variables suggest that no other extra macroscopic parameter may be required. VI. Stress tensor Before we focus on the force moment tensor, we will consider the stress tensor, σ, in order to under- stand the phenomenology of the force distribution in our tapped granular beds. We recall that Σ and σ are simply related through Σ ≡ V σ. However, we have to bear in mind that V is not a simple constant since the volume of the system depends, in a nontrivial way, on the shape and intensity of the excitation. In Fig. 8, we show the components of σ as a func- tion of Γ for different ν. As a reference, we show results of our simulations for a frictionless system. In a frictionless system, the shear vanishes and σyy is only determined by the weight of the sample since the Janssen’s effect is not present. As we can see, the frictionless sample presents a constant value of σyy for all Γ. For low Γ, the frictional samples dis- play values of σyy below the frictionless reference. This is a consequence of the Janssen’s effect since part of the weight of the sample is supported by the wall friction. Consequently, in this region, σxy is also positive [see Fig. 8(b)]. However, for each ν, there is a critical value of Γ, Γshear=0. Beyond it, the sample presents an apparent weight above the weight of the packing. In correspondence with this, σxy changes sign and becomes negative. This indicates that, for Γ > Γshear=0, the frictional walls are not supporting any weight. Rather, they pre- vent the packing from expansion by a downward frictional force. As Γ is increased beyond Γshear=0, the packing tends to store most of its stress in the horizontal direction (σxx) while σyy eventually sat- urates. For very intense pulses, the sample expands and lifts off significantly during the tap. When the bed falls back, it creates a very compressed struc- ture with most of the stress transmitted in the lat- eral directions and the wall friction sustaining the system downwards. It is worth mentioning that Γshear=0 is always higher than Γmin (see Fig. 6). 1 10 100 Γ 10 15 20 [u n it s o f m g /d ] 0.125 0.250 0.500 0.750 1.000 0.5 ( =0)µ ν (g/d) 1/2 σ σ σ yy xx (a) 1 10 100 Γ -0.8 -0.6 -0.4 -0.2 0 0.2 [u n it s o f m g /d ] 0.125 0.250 0.500 0.750 1.000 0.5 ( =0) xy σ ν (g/d) 1/2 µ x y (b) Figure 8: (a) Diagonal components of the stress tensor, σ, as a function of Γ for different frequen- cies ν of the tap pulse (the upper set of curves cor- responds to σyy and the lower set to σxx). (b) Off diagonal component of the stress tensor, σxy. The horizontal line corresponds to σyy [in panel (a)] and to σxy [in panel (b)] from simulations of a packing of frictionless disks. VII. The force moment tensor mas- ter curve Since the stress is not an extensive parameter, the force moment tensor is generally used to character- ize the macroscopic state [12,13]. Therefore, we will use Σ in the rest of the paper. Let us simply remark that since V presents a non-monotonic response to Γ, the curves in Fig. 8 present a somewhat differ- 030004-8 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. ent shape if Σ is plotted instead of σ. Particularly, Σyy does not display a minimum at low Γ, as the one observed for σyy in frictional disks but a mono- tonic increase. In Fig. 9, we show the trace of Σ as a function of Γ. There is a clear monotonic increase of Tr(Σ) as Γ is increased. Moreover, for a given Γ, if the frequency of the excitation pulse is increased, a significant reduction in the force moment tensor is observed. 1 10 100 Γ 28 30 32 34 36 38 40 T r( ) /N [ u n it s o f m g d ] 0.125 0.250 0.500 0.750 1.000 Σ ν (g/d) 1/2 Figure 9: Trace, Tr(Σ), of the force moment tensor as a function of Γ for different frequencies ν of the tap pulse. In Fig. 10, we plot the components of Σ as a func- tion of its trace for all the steady states generated. We can see that all data for different Γ and ν col- lapse into three master curves. This indicates that if two equilibrium states present the same Tr(Σ), all the components of Σ are also equal. Here, we point out a relevant piece of information that will be discussed in the next section. Two states may present equal force moment tensor but differ in vol- ume. This means that many points collapsing in Fig. 10 correspond to states of different φ. There- fore, at equilibrium, irrespective of the structure of the sample, two states with the same trace in Σ will present equal Σ. In a liquid at equilibrium, the stress tensor is diagonal and all elements along the diagonal are equal. This hydrostatic property allows us to know the full stress tensor if we only know the hydro- static pressure (i.e., if we only know the trace of the tensor). In our granular samples, the force mo- ment tensor can also be known if the trace is known. However, the shape of the tensor in static packings under gravity is defined by the three master curves of Fig. 10. 28 30 32 34 36 38 40 Tr( )/N [units of mgd] 10 12 14 16 18 20 /N [ u n it s o f m g d ] 0.125 0.250 0.500 0.750 1.000 Σ Σ Σ Σ ν(g/d) 1/2 yy xx (a) 28 30 32 34 36 38 40 Tr( ) [units of mgd] -0.8 -0.6 -0.4 -0.2 0 0.2 /N [ u n it s o f m g d ] 0.125 0.250 0.500 0.750 1.000 Σ Σ x y ν (g/d) 1/2 (b) Figure 10: (a) Diagonal components of the force moment tensor, Σ, as a function of Tr(Σ) for dif- ferent frequencies of the tap pulse (the upper set of curves corresponds to Σyy and the lower to Σxx). (b) Off diagonal component of the force moment tensor, Σxy. To our knowledge, there is no previous specula- tion that this property must hold for static granular packings. A more detailed study on the extent of this commonality of the shape of the force moment tensor will be pursued in a future paper [28]. How- ever, we show some suggestive preliminary results below. In Fig. 11, we show the components of Σ as a function of Tr(Σ) for a range of samples of differ- ent materials, for different tapping intensities, and for different tapping frequencies. As we can see, there is reasonable collapse of the data onto the same three master curves shown in Fig. 10. This is an indication that these master curves may be uni- versal and enclose a rather fundamental underlying property (inaccessible to us at this point) of static granular beds. 030004-9 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. 25 30 35 40 Tr( )/N [units of mgd] 8 10 12 14 16 18 20 / N [ u n it s o f m g d ] 0.125 0.5 300 0.250 0.5 300 0.500 0.5 300 0.750 0.5 300 1.000 0.5 300 0.250 0.0 300 0.250 0.2 300 0.250 1.0 300 0.250 5.0 300 0.250 50.0 300 0.250 0.5 30 0.250 0.5 3000 ν µ γ Σ Σ Σ Σ yy xx (a) 20 25 30 35 40 Tr( )/N [units of mgd] -0.6 -0.4 -0.2 0 0.2 0.4 / N [ u n it s o f m g d ] 0.125 0.5 300 0.250 0.5 300 0.500 0.5 300 0.750 0.5 300 1.000 0.5 300 0.250 0.0 300 0.250 0.2 300 0.250 1.0 300 0.250 5.0 300 0.250 50.0 300 0.250 0.5 30 0.250 0.5 3000 ν µ γ Σ Σ x y (b) Figure 11: (a) Diagonal components of force mo- ment tensor, Σ, as a function of Tr(Σ) for different frequencies of the tap pulse, different friction co- efficient and different restitution (the upper set of curves corresponds to Σyy and the lower to Σxx). (b) Off diagonal component of the force moment tensor, Σxy. The dashed lines are only a guide to the eyes. VIII. Force moment tensor fluctua- tions Since we have shown in the previous section that only the trace of Σ suffices (along with the master curves) to describe the full force moment tensor, we will now focus on this invariant and its fluctu- ations. In Fig. 12, the fluctuations of Tr(Σ) are plotted as a function of Γ. We obtain a single min- imum, in contrast with the minimum and maxi- mum observed in ∆φ. Interestingly, the states with minimum ∆Tr(Σ) correspond to the states where the minimums φ and ∆φ are reached for each ν. However, unlike ∆φ, the depth of the minimum in ∆Tr(Σ) is fairly independent of ν. It is unclear why the force moment fluctuations should present a minimum. Provided that the minimum of ∆Tr(Σ) coincides with the minimum of ∆φ, it can be spec- ulated that a reduced number of geometric con- figurations can accommodate a limited number of force configurations. We have seen that all individ- ual components of Σ present the same minimum in their fluctuations, however, the actual values in ∆Tr(Σ) are dominated by ∆Σxx which takes values five times larger than ∆Σyy. 1 10 100 1000 Γ 2 4 6 8 10 T r( ) /N [ u n it s o f m g d ] 0.125 0.250 0.500 0.750 1.000 ∆ Σ ν(g/d) 1/2 (a) 28 30 32 34 36 38 40 Tr( )/N [units of mgd] 2 4 6 8 10 T r( ) /N [ u n it s o f m g d ] 0.125 0.250 0.500 0.750 1.000 Σ Σ ∆ ν(g/d) 1/2 (b) Figure 12: (a) Fluctuations of the trace of the force moment tensor as a function of Γ for different fre- quencies of the tap pulse. (b) Fluctuations of the trace of the force moment tensor as a function of Tr(Σ). If we plot ∆Tr(Σ) in terms of the average value Tr(Σ) [see Fig. 12(b)], we can see that the curves collapse on top of each other over a wide range of Tr(Σ). However, some deviations are apparent at very low and very high forces. Although the fair collapse of the curves suggests that states of equal- Σ may correspond to the same equilibrium states, we will see in the next section that many of these equilibrium states are distinguishable through the 030004-10 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. volume. The inflection observed at high Tr(Σ) cor- responds to the change of regime observed in Fig. 8(a) where the vertical stress saturates and most of the contact forces are directed in the x–direction. IX. The thermodynamic phase space As we have suggested, the mean volume of static granular samples is not sufficient to describe the equilibrium state since states of equal-φ may present distinct fluctuations. On the other hand, the force moment tensor seems to be able to serve as a standalone descriptor since states of equal Σ do generally present the same Σ fluctuations. How- ever, states of equal-Σ may present different vol- umes. In Fig. 13, we plot the loci of the equilibrium states generated in our simulations in a hypothet- ical φ–Σ thermodynamic phase space. As we can see, states of equal V but different Σ are obtained as well as states of equal Σ and different V . 28 30 32 34 36 38 Tr( )/N [units of mgd] 0.81 0.82 0.83 0.84 0.85 0.86 φ 0.125 0.250 0.500 0.750 1.000 Σ ν (g/d) 1/2 Figure 13: Phase space φ–Tr(Σ). Loci visited in the simulations for different frequencies of the tap pulse. We can ask now if these two state variables suffice to fully describe the equilibrium states. A hint that this may be the case is given by the fact that states generated with different Γ and ν but that corre- spond to the same state in the φ–Σ plot display the same fluctuations of these variables. In Fig. 14, we have highlighted some pairs of neighboring states. We can see that such states do also present similar fluctuations [Fig. 14(b)]. In contrast, states which are distant in the φ–Σ plot present distinct fluctua- tions even if they correspond to an equivalent mean volume or an equivalent mean force moment tensor (see states joined by solid lines in Fig. 14). 28 30 32 34 36 38 Tr( )/N [units of mgd] 0.81 0.82 0.83 0.84 0.85 0.86 φ Σ (a) 3 4 5 6 7 Tr( )/N [units of mgd] 0.008 0.01 0.012 0.014 ∆ φ Σ∆ (b) Figure 14: (a) Phase space φ–Tr(Σ). (b) Fluctu- ations of the state variables. Selected neighboring states are colored in pairs for comparison. Distant states of equal V or Tr(Σ) are joined by thick lines. Let us point out here, that the sole coincidence of fluctuations in the macroscopic variables is not a rigorous proof that the set of chosen variables is a full complete set of thermodynamic param- eters. Future explorations of these systems may confirm or disprove that N, V and Tr(Σ) are a full set of macroscopic, extensible variables able to describe all equilibrium states. Meanwhile, it is clear that for moderate tapping intensities, around which the minimum in φ is observed, the approxi- mation of a simple NV or NΣ ensemble is not war- ranted in view of the large discrepancies between the curves generated with different pulse frequen- cies in Fig. 13. 030004-11 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. X. Concluding remarks We have studied steady states of mechanically sta- ble granular samples driven by tap like excitations. We have varied the external excitation by changing both, the pulse amplitude and the pulse duration. We have considered the macroscopic extensive vari- ables V (volume) and Σ (force moment tensor), and their fluctuations. From the results, we can draw the following conclusions: • There seems to be a rather robust set of master curves for Σαβ which implies that the knowl- edge of Tr(Σ) suffices to infer the other com- ponents of the force moment tensor. • The equilibrium states cannot be only de- scribed by V or Σ, apart from the number of particles N. • The equilibrium states seem to be well de- scribed by the set NV Tr(Σ). There exists a number of points to be considered in view of these findings. Here, we mention a few that may serve as starting points for future direc- tions of research: • What is the extent to which the Σ master curves are applicable? Is this dependent on the dimensionality of the system, the excita- tion procedure, the chosen contact force law, etc? • What is the dynamics during a single pulse that leads to the appearance of the φ mini- mum? Is this minimum present in states gen- erated with other types of pulses like fluidiza- tion or shear? The ubiquity of this minimum in simulation models [15, 20] suggests that it might be found in numerous conditions. • Are the fluctuations shown in Figs. 7 and 12 the definitive phenomenological equation of states? Other authors have so far found mono- tonic density fluctuations [27] or concave up density fluctuations [6]. • How much of the φ–Tr(Σ) plane can be ex- plored by changing material properties? • Are there other excitation protocols (such as shearing) that may give rise to steady states that are thermodynamically equivalent to the ones obtained by tapping? • Is it possible to construct a phenomenological entropy function from the equations of states (Figs. 7 and 12) by simple integration of a Gibbs–Duhem-like equation? Let us bear in mind that Fig. 7 is multiply valued. It is worth stressing that if two ensembles gen- erated by arbitrary excitation protocols —such as tapping or shearing— happen to present the same mean values (and fluctuation) for all macroscopic variables, then such macroscopic states should be considered thermodynamically identical. However, it may be the case that a given protocol produces a narrow range of macroscopic states that can be eventually described with a reduced set of macro- scopic variables. Acknowledgements - LAP acknowledges dis- cussions with Massimo Pica Ciamarra. JD ac- knowledges a scholarship of the FPI program from Ministerio de Ciencia e Innovacin (Spain). This work has been financially supported by CONICET (Argentina), ANPCyT (Argentina), Project No. FIS2008-06034-C02-01 (Spain) and PIUNA (Univ. Navarra). [1] H B Callen, Thermodynamics and an Intro- duction to Thermostatistics, 2nd ed., Wiley- VCH, New York (1985). [2] R K Pathria, Statistical Mechanics, 2nd ed., Butterworth-Heinemann, Oxford (1996). [3] E R Nowak, J B Knight, M L Povinelli, H M Jeager, S R Nagel, Reversibility and irre- versibility in the packing of vibrated granular material, Powder Tech. 94, 79 (1997). [4] P Richard, M Nicodemi, R Delannay, P Ribire, D Bideau, Slow relaxation and compaction of granular systems, Nature Mater. 4, 121 (2005). [5] Ph Ribière, P Richard, P Philippe, D Bideau, R Delannay, On the existence of stationary states during granular compaction, Eur. Phys. J. E 22, 249 (2007). 030004-12 Papers in Physics, vol. 3, art. 030004 (2011) / L. A. Pugnaloni et al. [6] M Schröter, D I Goldman, H L Swinney, Sta- tionary state volume fluctuations in a granular medium, Phys. Rev. E 71, 030301(R) (2005). [7] S F Edwards, R B S Oakeshott, Theory of pow- ders, Physica A 157, 1080 (1989). [8] J H Snoeijer, T J H Vlugt, W G Ellenbroek, M van Hecke, J M J van Leeuwen, Ensemble the- ory for force networks in hyperstatic granular matter, Phys. Rev. E. 70, 061306 (2004). [9] S F Edwards, The full canonical ensemble of a granular system, Physica A 353, 114 (2005). [10] S Henkes, C S OHern, B Chakraborty, Entropy and temperature of a static granular assembly: An ab initio approach, Phys. Rev. Lett. 99, 038002 (2007). [11] S Henkes, B Chakraborty, Stress correlations in granular materials: An entropic formula- tion, Phys. Rev. E. 79, 061301 (2009). [12] R Blumenfeld, S F Edwards, On granular stress statistics: Compactivity, angoricity, and some open issues, J. Phys. Chem. B 113, 3981 (2009). [13] B P Tighe, A R T van Eerd, T J H Vlugt, En- tropy maximization in the force network en- semble for granular solids, Phys. Rev. Lett. 100, 238001 (2008). [14] R Arévalo, D Maza, L A Pugnaloni, Identi- cation of arches in two-dimensional granular packings, Phys. Rev. E 74, 021303 (2006). [15] L A Pugnaloni, M Mizrahi, C M Carlevaro, F Vericat, Nonmonotonic reversible branch in four model granular beds subjected to vertical vibration, Phys. Rev. E 78, 051305 (2008). [16] J Schäfer, S Dippel, D E Wolf, Force schemes in simulations of granular materials, J. Phys. I (France) 6, 5 (1996). [17] H Hinrichsen, D Wolf, The physics of granular media, Willey-VCH, Weinheim (2004). [18] J A Dijksman, M van Hecke, The role of tap duration for the steady-state density of vibrated granular media, Eur. Phys. Lett. 88, 44001 (2009). [19] F Ludewig, S Dorbolo, T Gilet, N Vandewalle, Energetic approach for the characterization of taps in granular compaction, Eur. Phys. Lett. 84, 44001 (2008). [20] P A Gago, N E Bueno, L A Pugnaloni, High intensity tapping regime in a frustrated lattice gas model of granular compaction, Granular Matter 11, 365 (2009). [21] Notice however, that this expansion parameter is not based on the properties of the mechani- cal excitation itself but on the sample response to it. [22] J Damas et al., (unpublished) [23] According to Callen [1], a system is at equilib- rium if Thermodynamics holds for such state. [24] Notice, that in Fig. 1(c) of our previous work [26], the maximum value of Γ reported was half of the value shown here. This discrepancy is due to a deficient filtering of the accelerom- eters used in Ref. [26]. The accelerations re- ported here have been corrected and checked against measurements done with a high speed camera. [25] C M Carlevaro, L A Pugnaloni, Steady state of tapped granular polygons, J. Stat. Mech. P01007 (2011). [26] L A Pugnaloni, I Sánchez, P A Gago, J Damas, I Zuriguel, D Maza, Towards a relevant set of state variables to describe static granular pack- ings, Phys. Rev. E 82, 050301(R) (2010). [27] M Pica Ciamarra, A Coniglio, M Nicodemi, Thermodynamics and statistical mechanics of dense granular media, Phys. Rev. Lett. 97, 158001 (2006). [28] L. A. Pugnaloni, et al., (unpublished). 030004-13