Papers in Physics, vol. 8, art. 080001 (2016) Received: 2 November 2015, Accepted: 22 December 2015 Edited by: C. S. O’Hern Reviewed by: A. Baule, Queen Mary University of London, UK. Licence: Creative Commons Attribution 3.0 DOI: http://dx.doi.org/10.4279/PIP.080001 www.papersinphysics.org ISSN 1852-4249 Ergodic–nonergodic transition in tapped granular systems: The role of persistent contacts Paula A. Gago,1–3 Diego Maza,4 Luis A. Pugnaloni1, 2∗ Static granular packs have been studied in the last three decades in the frame of a modified equilibrium statistical mechanics that assumes ergodicity as a basic postulate. The canon- ical example on which this framework is tested consists in the series of static configurations visited by a granular column subjected to taps. By analyzing the response of a realistic model of grains, we demonstrate that volume and stress variables visit different regions of the phase space at low tap intensities in different realizations of the experiment. We show that the tap intensity beyond which sampling by tapping becomes ergodic coincides with the forcing necessary to break all particle–particle contacts during each tap. These results imply that the well-known “reversible” branch of tapped granular columns is only valid at relatively high tap intensities. I. Introduction Granular matter is ubiquitous in nature. How- ever, due to the complexity of the real particle– particle interactions, the standard approaches of continuum mechanics and thermodynamics are still limited in providing meaningful descriptions of the states in which these systems can be. Edwards and Oakeshot introduced a tentative approach inspired by the ideas of equilibrium statistical mechanics to formally describe the global properties of a static ∗E-mail: luis.pugnaloni@frlp.utn.edu.ar 1 Dpto. Ingenieŕıa Mecánica, Facultad Regional La Plata, Universidad Tecnológica Nacional, Av. 60 Esq. 124, 1900 La Plata, Argentina. 2 Consejo Nacional de Investigaciones Cient́ıficas y Técnicas, Argentina. 3 Current address: Department of Earth Science and En- gineering, Imperial College London, South Kensington Campus, London SW7 2AZ, UK. 4 Departamento de F́ısica y Matemática Aplicada, Facul- tad de Ciencias, Universidad de Navarra, Navarra, Spain. granular pack. Since the introduction of this the- ory —where the entropy of the systems is governed by the spatial disorder of the grains [1]—, a num- ber of studies have used it to frame the interpre- tation of the results of specific experiments. The most relevant case is the so-called “Chicago exper- iment”, where a column of grains was repeatedly tapped following an annealing-type protocol [2, 3]. The main outcome of this experiment is that a sta- tionary state can be reached, where the mean vol- ume fraction, φ, is a well defined function of the tap amplitude, Γ. Others have also obtained seemingly reproducible states without the need of annealing [4]. However, it has been shown recently, by sim- ulation of frictionless grains, that these stationary states are not necessarily ergodic [5]. At low Γ, different members of an ensemble of steady-states prepared with a well defined protocol may sample a different region of the phase space, as the fluctu- ations of φ indicate. In this paper, we demonstrate that not only the volume but also the force moment tensor, Σ, are sampled in a non-ergodic fashion and that ergod- 080001-1 Papers in Physics, vol. 8, art. 080001 (2016) / P. A. Gago et al. icity seems to be recovered if all particle–particle contacts are lost during each tap. This sets a clear limit to the range of driving forces able to generate a sequence of configurations for which the Edwards framework can be applied. II. Numerical protocol We simulated using the LAMMPS package [6] a quasi-two-dimensional cell containing N = 1000 spherical particles of diameter d. The cell is 1.1 d thick and 27.8 d wide (the granular column is about 35 layers deep) to have a one to one representa- tion of a previously introduced experimental device [7, 8]. We use a model for soft frictional spheres de- scribed in Refs. [9, 10]. The normal component, Fn, of the contact interaction is given by an elas- tic repulsive force proportional to the overlap of the interacting spheres and a dissipative term pro- portional to the normal component of the relative velocity. The tangential term, Ft, implements an elastic shear force and a damping force. The shear force takes into account the cumulated tangential displacement between the particles while they re- main in contact. Whenever Ft > µFn (µ is the fric- tion coefficient), this lower dynamic friction force is used. In this work we use the same interaction parameters as in Ref. [11, 12]. The wall–particle interaction is defined with the same parameters as the particle–particle force. Tapping is simulated by imposing an external vertical motion to the cell. This pulse is a single sinusoidal cycle A sin(ωt). We fix ω = 2π×33 Hz and use the tap amplitude A as control parameter. The tap intensity is character- ized by Γ = Aω2/g. The mechanical equilibrium after each tap is deemed achieved if the kinetic en- ergy of each particle has fallen (in average) below 10−6mgd. Where m is the mass of one particle and g the acceleration of gravity. We study 20 independent realizations of a de- creasing ramp of the tap amplitude. We initially fill the cell by placing the spheres at random posi- tions before letting them deposit under the action of gravity. In each realization, we decrease Γ in small steps, from Γ = 20.0 down to Γ = 0.8, and apply 200 taps for each Γ. Note that for Γ < 1.0, the column of grains does not detach from the base during a tap. The 200 taps at each value of Γ are enough to reach a steady-state. We do not observe any drift of the mean values of φ or Σ after the initial 100 taps, which we will discard later in our analysis. Finally, we also study a cyclic anneal- ing protocol: starting from the final configuration at Γ = 0.8 for each of the former realizations; the tap amplitude is cyclically increased and decreased every 200 taps in the range 0.8 < Γ < 5.5 in or- der to compare the steady-states reached with an alternative method. III. Data analysis To measure the packing fraction we use the 2D Voronoi tessellation (implemented in [13]) of the x–z plane projection of the particle positions, dis- regarding the third coordinate on the thin direc- tion of the cell. Then, we associate to each particle a “local volume” fraction by dividing the particle area by the corresponding Voronoi area. In order to avoid boundary effects, we disregard particles closer than 2d to the lateral walls. Following the rec- ommendations in Ref. [14], we analyze horizontal slices of the granular column 15 d thick measured at approximately the same depth with respect to the free surface in order to retrieve unbiased results for the force moment tensor due to the uneven free sur- face. Averaging over the N particles contained in the slice of interest, we obtain the volume fraction of each static configuration. To obtain the steady- state φ corresponding to a given Γ, we averaged this quantity over the last 100 configurations ob- tained for each tap intensity. We also obtain the force moment tensor σ αβ i = ∑ c rαc f β c of each parti- cle in the slice of interest. Here, the sum runs over all the contacts c of the particle, −→r c is the vector from the center of the grain to the contact c and −→ f c is the corresponding contact force. We apply the same averaging protocol used for φ to obtain the force moment tensor for the configuration and Σ for the steady-state of a given realization and Γ. IV. Results In Fig. 1(a) the ensemble average of the steady- state φ (i.e., averaged over the 20 independent re- alizations) is displayed as a function of Γ. The error bars correspond to the standard deviation over the 20 realizations. As observed by a number of au- thors, the curve seems to be very well defined with 080001-2 Papers in Physics, vol. 8, art. 080001 (2016) / P. A. Gago et al. 0.895 0.9 1 2 3 φ Γ 0.84 0.88 0.92 φ (a) 16 20 24 28 1 2 3 T r (Σ ) N -1 [m d g ] Γ 16 20 24 28 5 10 15 20 T r (Σ ) N -1 [m d g ] Γ (b) Figure 1: Ensemble average of the steady-state packing fraction φ (a) and trace of the force mo- ment tensor Tr(Σ) (b) as a function of Γ. The error bars correspond to the standard deviation over the 20 averaged realizations. The insets show the same results and two of the 20 independent realizations (dashed lines) in the low-tap-intensity region. The error bars on the single realization data correspond to the estimated standard error of the mean. independent realizations falling within a very nar- row range of φ values for any given Γ. In the past, this led to the conclusion that this was a truly re- versible process, where lowering or raising Γ would lead to the same steady-state φ. In the inset of Fig. 1(a) two of the independent realizations are shown for the low-tap-intensity region. From this picture, it is clear that steady-states correspond- ing to a given Γ can differ from one realization to another. Notice that in the inset the error bars for the two isolated realizations correspond to the standard error of the mean (SEM), which gives an estimate of the uncertainty of the mean value re- ported rather than the size of the φ fluctuations. For these two realizations, although the mean φ seems to agree within the estimated error for inten- sities Γ > 1.5, it is clear that they are different for 1.0 1.3 0.900 0.901 0.902 0.903 0.904 0.905 φ (a) 0 500 1000 1500 2000 Number of taps 16 18 20 22 24 26 28 T r ( Σ ) N −1 [m gd ] (b) Figure 2: φ (a) and Tr(Σ) (b) as a function of the tap number for two of the 20 independent realiza- tions at Γ = 0.8. The notched boxes and “vio- lin” diagrams shown suggest that both realizations can be hardly considered as representing the same steady-state. low Γ. This is consistent with the findings of Pail- lusson and Frenkel [5] for frictionless spheres under event-driven simulations. However, in our simula- tions we are able to extract the stress state of the system as well as the history of the contacts. These reveal valuable information, as we discuss below. In Fig. 1(b), we show the trace Tr(Σ) of Σ aver- aged over all 20 realizations as a function of Γ. As before, the error bars indicate the standard devi- ation over the 20 realizations. As we can see, the variability of the mean stress is significantly large at low Γ. This is not due to the large fluctuations during a given series of taps but to the variations observed from one realization of the protocol to the other. Indeed, the inset in Fig. 1(b) shows that, for low Γ, the mean values of Tr(Σ) for two of the 20 realizations have a relatively small SEM (i.e., the fluctuations in each realization are small). How- ever, realizations differ from each other. The dif- ference between results corresponding to different realizations become much more evident here than in the case of the φ–Γ plot. We suggest that the stress tensor may be more sensitive and then more suitable to sense if ergodicity is fulfilled in exper- imental data. Overall, as Γ is decreased, differ- ent realizations explore non-overlapping ranges of volume/stress. Therefore, temporal averages (on 080001-3 Papers in Physics, vol. 8, art. 080001 (2016) / P. A. Gago et al. 0.82 0.84 0.86 0.88 0.9 φ (a) 0.901 0.905 1 1.5 1.9 φ Γ 18 21 24 27 1 2 3 4 5 T r (Σ ) N -1 [m d g ] Γ (b) Figure 3: Steady-state φ (a) and Tr(Σ) (b) as a function of Γ corresponding to a full annealing pro- tocol. Starting from the filled circles (red) increas- ing ramp, and following by filled squares (blue), open circles (green) and open squares (magenta). The error bars correspond to the estimated SEM. a single time series) do not match with ensemble averages (over the realizations). Since the number of taps we have explored for each Γ may be small to assume that the steady- state has been properly sampled, we carried out 2000 additional taps at Γ = 0.8 for each realiza- tion. Since some of the signals are not normally distributed, we confirmed the stationarity of these states by using a non-parametric test at a level of significance of 5% [15, 16]. Two of the normally distributed realizations are displayed in Fig. 2 as a function of the tap number. Comparing the cor- responding notched boxes and “violin” diagrams of both signals, it is clear that the states do not match. Hence, the reversible branch found in Ref. [3] is not so at low Γ in our case since truly station- ary states with distinguishable φ and Tr(Σ) may be obtained on different realizations of the annealing protocol. Of course, this is harder to detect in φ since the dispersion between realizations is much smaller compared with the range of φ values ob- tained at different Γ. The former results also confirm the hypothesis that non-ergodicity is present in a typical tapping protocol beyond the special case reported in Ref. [5], where the steady-states obtained did not follow any annealing-type protocol. Hence, we observe this non-ergodic behavior even after annealing the system from high tapping strengths. To stress this point and assess if the speed of the annealing may prevent the system from reaching a unique steady- state on each realization, we apply a slower cyclic annealing protocol (similar to the one introduced in [2]) to each of the 20 final states at low Γ in order to reproduce the “reversible branch”. In Fig. 3 we display a sequence of two successive up and down ramps applied to one of the 20 initial realizations using Γ-steps about one half of those used in Fig. 1. Although in the scale used for φ the steady-state packing fraction seems reversible, a close inspection shows that the states have distinguishable φ at low Γ [see Fig. 3(a) and the corresponding inset]. This is much more evident when the stress is analyzed [see Fig. 3(b)]. In order to set a criterion to decide if the steady- states are not ergodic for a given Γ, we show in Fig. 4 the p-values for the Kruskal–Wallis [17] one-way analysis of variance performed on the 20 realiza- tions at each Γ. This simple non-parametric test allows for the rejection of the null hypothesis that all 20 data series are drawn from a unique distri- bution (which does not need to be normal), hence that they correspond to a unique steady-state. If the p-value is significant (in our case p > 0.01), then we cannot rule out the possibility that the 20 series come from the same steady-state. As we can see from Fig. 4(a), the test run on the data for φ indicates that for Γ < 5.0, the null hypothesis must be rejected and therefore there exist at least two out of the 20 steady-states that are not the same. However, for higher Γ, the test is significant and then the 20 realizations may correspond to the same steady-state. Interestingly, when the test is run on Σ [see Fig. 4(b)], the steady-state seems to be unique for all 20 realizations if Γ > 3.75. Al- though differences between realizations are simpler to detect on visual inspection of Tr(Σ), it is actu- ally φ that sets a higher threshold for the Γ val- ues needed to ensure an ergodic steady-state (i.e., Γ > 5.0). The previous results indicate that the steady- states sampled at low tap intensities do not only 080001-4 Papers in Physics, vol. 8, art. 080001 (2016) / P. A. Gago et al. 0.83 0.85 0.87 0.89 â��10 -11 â��10 -6 â��10 -2 â��10 0 φ p -v a lu e (a) p-value (φ) φ 17 20 23 5 10 15 20 25 â��10 -11 â��10 -6 â��10 -2 â��10 0 T r (Σ ) N -1 [m d g ] p -v a lu e Γ (b) p-value (Tr(Σ)) Tr(Σ) Figure 4: P-value (up-triangles, right axis) as a function of Γ for the Kruskal–Wallis [17] one-way analysis of variance for φ (a) and Tr(Σ) (b). The horizontal dotted line corresponds to the signifi- cance level used (1%). The black circles correspond to the φ and Tr(Σ) data from Fig. 1. 0 1 2 3 3 6 9 12 15 18 21 C /C 0 [ % ] Γ (a) p-value (φ ) p-value (Tr(Σ)) C/C0 3 6 9 12 15 18 21 â��10 -11 â��10 -6 â��10 -2 â��10 0 p -v a lu e Γ (b) Figure 5: (a) Percentage C/C0 × 100 of persistent contacts (black circles) as a function of Γ, averaged over 5 taps on 6 independent realizations. The er- ror bars correspond to the standard deviation over the 6 realizations. Up-triangles (green) correspond to the p-values in Fig. 4 for φ and down-triangles (magenta) to Tr(Σ). (b) Same as (a) for frictionless grains. depend on Γ but on the particular history of each realization. Notice that this goes beyond the his- tory dependent out-of-equilibrium trajectories al- ready reported in tapped systems [18] since here we are focusing on the steady-states. One may hypoth- esize that the constraints imposed by the contacts is one of the reasons for the non-ergodic behavior at low tap. If a contact persists from one tap to the next, the contact force after coming back to rest will depend on the history of all contacts of that particular grain. In order to test this idea, we analyze the evolution of all contacts during each tap and identify those that persist (i.e., contacts that did not break at any time during the pulse of energy). Figure 5 shows the average ratio C/C0 of persistent contacts, C, to the total number of contacts, C0, as a function of Γ. For this calcu- lation, each contact was tracked during the final 5 taps for each Γ on 6 of the independent realiza- tions and only grains that fall within the layer of interest, as discussed above, where included in the analysis. The percentage of persistent contacts is very small but non-zero up to Γ ≈ 5.0. As it is expected, when Γ is increased sufficiently all the contacts are broken and new ones are made during each tap (resulting in C/C0 = 0). This transition coincides with the value of Γ where the realizations seem to sample the same steady-state (see the p- values included in Fig. 5). Therefore, when small taps are applied, the aging of some of the contacts seem to lead the system to sample different regions of the phase space during independent realizations. However, if all contacts are made anew at each tap, the sampling becomes compatible with the idea of ergodicity introduced in Fig. 4. In order to general- ize this result, we also simulated frictionless grains. Interestingly, the same conclusion drawn for fric- tional grains is true for frictionless ones: different realizations seem to sample the same steady-state only if all contacts are made anew upon each tap [see Fig. 5 (b)]. V. Conclusions Our analysis of the steady-states of tapped granu- lar systems indicate that these states are history- dependent for tap intensities below a certain threshold. This is in contradiction with the general assumption that macroscopic time averages —such as the volume fraction— can be recovered when the amplitude of the perturbation applied to the 080001-5 Papers in Physics, vol. 8, art. 080001 (2016) / P. A. Gago et al. system is tuned back and forth. The differences between independent realizations become particu- larly noticeable in the stress distribution. These findings show that the postulates of the equilib- rium statistical thermodynamics may not be al- ways fulfilled to describe the steady state of static granular systems (see also Ref. [19] for a discus- sion on the Boltzmann distribution failure for an analytically solvable model). Focusing on tap in- tensities that warrant that all contacts are made anew after each tap may allow exploring the avail- able phase space in agreement with the ergodic hy- pothesis. However, gentle perturbations deserve an approach that includes memory effects to suitably describe the states. In that sense, non-equilibrium thermodynamic approaches may be a suitable al- ternative [20]. Further research on such alterna- tive formalisms, the effect of other types of forcing mechanisms (e.g., shear), and possible extensions to other complex systems (e.g., active matter) be- come necessary. Acknowledgements - This work has been partially supported by projects PICT-2012-2155 ANPCyT (Argentina), FIS2011-26675 MINECO (Spain) and PIUNA (Universidad de Navarra). [1] S F Edwards, R B Oakeshot, Theory of pow- ders, Physica 157D, 1091 (1989). [2] E R Nowak, J B Knight, E Ben-Naim, H M Jaeger, S R Nagel, Density fluctuations in vi- brated granular materials, Phys. Rev. E 57, 1971 (1998). [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 Technol. 94, 79 (1997). [4] 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). [5] F Paillusson, D Frenkel, Probing ergodicity in granular matter, Phys. Rev. Lett. 109, 208001 (2012). [6] S Plimpton, P Crozier, A Thompson, LAMMPS-large-scale atomic/molecular mas- sively parallel simulator, Sandia National Lab- oratories (2007). [7] L A Pugnaloni, I Snchez, 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 (2010). [8] S Ardanza-Trevijano, I Zuriguel, R Arévalo, D Maza, Topological analysis of tapped granular media using persistent homology, Phys. Rev. E 89, 052212 (2014). [9] N V Brilliantov, F Spahn, J-M Hertzsch, T Pöschel, Model for collisions in granular gases, Phys. Rev. E. 53, 5382 (1996). [10] L E Silbert, D Ertaş, G S Grest, T C Halsey, D Levine, S Plimpton, Granular flow down an inclined plane: Bagnold scaling and rheology, Phys. Rev. E 64, 051302 (2001). [11] L A Pugnaloni, J Damas, I Zuriguel, D Maza, Master curves for the stress tensor invariants in stationary states of static granular beds. Im- plications for the thermodynamic phase space, Papers in Physics 3, 030004 (2011). [12] 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). [13] C Rycroft, VORO++: A three-dimensional Voronoi cell library in C++, Chaos 19, 041111 (2009). [14] P A Gago, L A Pugnaloni, D Maza, Relevance of system size to the steady-state properties of tapped granular systems, Phys. Rev. E 91, 032207 (2015). [15] W Constantine, D Percival, Fractal: Frac- tal Time Series Modeling and Analysis. R package version 2.0-0, http://CRAN.R- project.org/package=fractal (2014). [16] R Core Team, R: A language and environ- ment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria (2013). 080001-6 Papers in Physics, vol. 8, art. 080001 (2016) / P. A. Gago et al. [17] W H Kruskal, W A Wallis, Use of ranks in one-criterion variance analysis, J. Am. Stat. Assoc. 47, 583 (1952). [18] Ch Josserand, A V Tkachenko, D M Mueth, H M Jaeger, Memory effects in granular ma- terial, Phys. Rev. Lett. 85, 3632 (2000). [19] R M Irastorza, C M Carlevaro, L A Pug- naloni, Exact predictions from the Edwards en- semble versus realistic simulations of tapped narrow two-dimensional granular columns, J. Stat. Mech. P12012 (2013). [20] G Lebon, D Jou, JC Casas-Vazquez, Un- derstanding Non-equilibrium Thermodinam- ics, Springer-Verlag, Berlin Heidelberg (2008). 080001-7