Master curves for the stress tensor invariants in stationary states of static granular beds. Implications for the thermodynamic phase space

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, \Sigma. 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 \Sigma, 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 \Sigma may be sufficient to define the macroscopic states. Moreover, we show that only one of the invariants of \Sigma is necessary, since each component of \Sigma falls onto a master curve when plotted as a function of Tr(\Sigma). 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 fundamental importance to the industry to improve the storage of such materials in bulk, as well as to optimize the packaging design. However, far from yielding such benefits of practical interest, physicists have found a fascinating challenge on their way that has stuck them, to a large extent, in the 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 different 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, averages over a large number of replicas at a given time should give equivalent results to time averages [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 foundations are still supported by little evidence. 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][4][5][6]. In such protocols, the sample is subjected to external excitations of the form of pulses. Well defined, reproducible time averages 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 generated by using a given external excitation after any transient has disappeared.
More than twenty years ago, Edwards and Oakshott [7] put forward the idea that the number of grains N and the volume V are the basic state variables that suffice to characterize a static sample of hard grains in equilibrium. The N V granular ensemble was then introduced as a collection of microstates, where the sample is in mechanical equilibrium, compatible with N and V . However, newer theoretical works [8][9][10][11][12][13] suggest that the force moment tensor, Σ, (Σ = V σ, where σ is the stress tensor) must be added to the set of extensive macroscopic variables (i.e., an N V Σ ensemble) to adequately describe a packing of real grains.
In the rest of this paper, we show experimental and simulation evidence that the equilibrium states of static granular packings cannot be only described by V (or equivalently the packing fraction, φ, defined as the fraction of the space covered by the grains) nor by Σ. We do this by generating states of equal V but different Σ and states of equal Σ but different φ. We also show that states of equal V may present different volume fluctuations. Moreover, we show that states of equal V and Σ display the same fluctuations of these variables, 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 − |r ij | and the velocitiesṙ ij , ω i and ω j . Here, r ij represents the center-to-center vector between particles i and j, d is the particle diameter and ω is the particle angular velocity. These forces are introduced in the Newton's translational and rotational equations of motion and then numerically integrated by a velocity 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 F n and a tangential force F t . where The first term in Eq. (1) corresponds to a restoring force proportional to the superposition ξ of the interacting disks and the stiffness constant k n . The second term accounts for the dissipation of energy during the contact and is proportional to the normal component v n i,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 following a rule that selects between static or dynamic friction. Notice that Eq. (2) implies that the maximum static friction force |F s | used corresponds to µ|F n |, which effectively sets µ dynamic = µ static = µ. The static friction force F s [see Eq. (3)] has an elastic term proportional to the relative shear displacement ζ and a dissipative term proportional to the tangential component v t i,j of the relative velocity. In Eq. (5), s is a unit vector normal to r ij . The elastic and dissipative contributions are characterized by k s and γ s , respectively. The shear displacement ζ is calculated through Eq. (4) by integrating v t i,j from the beginning of the contact (i.e., t = t 0 ). The tangential interaction behaves like a damped spring which is formed whenever two grains come into contact and is removed when the contact finishes [17].
The particular set of parameters used for the simulation is (unless otherwise stated): µ = 0.5, k n = 10 5 (mg/d), γ n = 300(m g/d), k s = 2 7 k n and γ s = 200(m g/d). In some cases, we have varied µ and γ n in order to control the friction and restitution 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 stability of particle contacts to decide when the system 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 rectangular 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 window. 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 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 symmetric 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 calculated 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.  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 diameter 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-tromagnetic shaker (Tiravib 52100) with a trend of half sine wave pulses separated five seconds between them. The tapping amplitude was controlled adjusting the intensity, A, and frequency, ν, of the pulse and was measured by an accelerometer attached to the base of the cell. Averages are taken over 500 taps after equilibration.

III. Experimental method
High resolution images (more than 10 MPix) are taken after each tap. To calculate the packing fraction, 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 calculate 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 separation between plates is larger than the particle diameter, small overlaps between the spheres are present in the 2D projections of the images. Therefore, 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 appropriate parameter to characterize the external excitation used to drive a granular sample. Dijksman et al. [18] proposed a parameter related to the liftoff velocity of the granular bed. Ludewing et al. [19] presented an energy based parameter. Pugnaloni et al. [15,20] suggested that the factor of expansion induced on the sample would be a suitable 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). Parameters 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 simple 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 simplified 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 frequency, are needed to fully describe the excitation. We will employ the usual nondimensional peak acceleration Γ ≡ a peak /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 definition 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. We have generated configurations corresponding to a particular pulse (of a given shape and intensity) by repeatedly applying such pulse to the system and allowing enough time for any transient to fade. To prove that our samples are at equilib-rium, we approach a particular pulse through different paths -and starting from different initial conditions (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 configuration corresponds to a very ordered structure. In Fig. 2, we consider the equilibration process from a highly disordered initial configuration. For low tapping 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 initial convergence, a slower compaction phase takes the system to the final steady state. For high tapping 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 simulations, 400 taps of equilibration were enough.
An interesting example of equilibration is presented 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 configuration 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 because 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 consisting of: (i) application of a large number of identical pulses starting from a disordered configuration, and (ii) application of a reduced number of identical pulses after annealing the system from much higher pulse intensities. In a few cases, the results (mean values and/or fluctuations) from both protocols 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 happens 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 function of Γ from our experiments for different pulse frequencies [24]. As we can see, there exists a min-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). Although the increment in the packing fraction beyond the minimum is rather small, it is important 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 packing 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 A similar analysis was done in Ref. [27] for states generated with different pulse amplitude and duration in liquid fluidized beds. The fluctuations of φ in the steady state, as measured by the standard deviation ∆φ, are presented in Fig. 5. As we can see, a minimum in the fluctuations is just apparent. The position of such minimum in the fluctuations coincides with the minimum 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 interval (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. However, 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 frequencies. However, rescaling φ and Γ with φ min and Γ 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 volume 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 tapping intensities that warrant the steady state with a modest number of pulses. The value of Γ, at which the fluctuations display a minimum, coincides with Γ min , the value at which the minimum packing fraction, φ min , is obtained. The maximum coincides with the inflection point in the φ-Γ curve at higher Γ. Since one expects 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 intermediate packing fractions. This is more clearly seen in Fig. 7 where we plot the fluctuations as a function of the average value of φ.  Figure 7 presents two distinct branches: the lower branch corresponds to Γ > Γ min and the upper branch to Γ < Γ min . For Γ > Γ min , the fluctuations 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 parameter, 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 discrepancies. Although this indicates that the states are not equivalent, it also suggests that such descriptors are not good candidates to form a set of macroscopic variables, along with V , to uniquely 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 suggested by some theoretical speculations [12]. However, 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 understand 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 function 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 display 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 prevent 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 saturates. For very intense pulses, the sample expands and lifts off significantly during the tap. When the bed falls back, it creates a very compressed structure with most of the stress transmitted in the lateral directions and the wall friction sustaining the system downwards. It is worth mentioning that Γ shear=0 is always higher than Γ min (see Fig. 6).

VII. The force moment tensor master curve
Since the stress is not an extensive parameter, the force moment tensor is generally used to characterize 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-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 monotonic 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. In Fig. 10, we plot the components of Σ as a function of its trace for all the steady states generated. We can see that all data for different Γ and ν collapse 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 volume. This means that many points collapsing in Fig. 10 correspond to states of different φ. Therefore, 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 hydrostatic pressure (i.e., if we only know the trace of the tensor). In our granular samples, the force moment 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. To our knowledge, there is no previous speculation 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]. However, 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 different 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 universal and enclose a rather fundamental underlying property (inaccessible to us at this point) of static granular beds.

VIII. Force moment tensor fluctuations
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 fluctuations. In Fig. 12, the fluctuations of Tr(Σ) are plotted as a function of Γ. We obtain a single minimum, in contrast with the minimum and maximum 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 speculated that a reduced number of geometric configurations can accommodate a limited number of force configurations. We have seen that all individual 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 . 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 volume. The inflection observed at high Tr(Σ) corresponds 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. However, states of equal-Σ may present different volumes. In Fig. 13, we plot the loci of the equilibrium states generated in our simulations in a hypothetical φ-Σ thermodynamic phase space. As we can see, states of equal V but different Σ are obtained as well as states of equal Σ and different V . 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 correspond 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). 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 parameters. 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 approximation of a simple N V or N Σ ensemble is not warranted in view of the large discrepancies between the curves generated with different pulse frequencies in Fig. 13.

X. Concluding remarks
We have studied steady states of mechanically stable 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 variables 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 knowledge of Tr(Σ) suffices to infer the other components of the force moment tensor.
• The equilibrium states cannot be only described by V or Σ, apart from the number of particles N .
• The equilibrium states seem to be well described by the set N V 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 directions of research: • What is the extent to which the Σ master curves are applicable? Is this dependent on the dimensionality of the system, the excitation procedure, the chosen contact force law, etc?
• What is the dynamics during a single pulse that leads to the appearance of the φ minimum? Is this minimum present in states generated with other types of pulses like fluidization 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 monotonic density fluctuations [27] or concave up density fluctuations [6].
• How much of the φ-Tr(Σ) plane can be explored 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 generated 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 macroscopic variables.