Methods.

We employ a numerical implementation inspired by the experimental setup used for the Chicago experiment. Our pile consists of 60,000 bidisperse spheres 1 to 1.02 mm in diameter contained in a vertical cylinder of 2.4-cm diameter. The resulting column height is about 11.5 cm. A minute bidispersity, in equal proportion, was introduced in order to prevent crystallization and the grain diameter ratio was chosen to avoid segregation.

We use MD, in particular the implementation of the soft spheres (25) model provided by the LIGGGHTS (LAMMPS improved for general granular and granular heat transfer simulations) (26) open source software. Within LIGGGHTS, we apply a Hertz model for the grain–grain and grain–wall contact forces (with Young’s modulus $$ Y=6\times 10^6{\mathrm{N}\mathrm{m}}^{-2}$$, Poisson ratio $$ n=0.3$$, restitution coefficient $$ ϵ=0.5$$, and friction coefficient $$ \mu =0.5$$). Particle densities equal 2,500 kg/$$ {\mathrm{m}}^3$$.

Taps consist of moving the whole container along a semisinusoidal wave pulse of amplitude $$ A$$ and angular frequency $$ \omega $$. A new tap will be applied to the system only after the system has attained a mechanically stable state. The system will be considered mechanically static when the kinetic energy of the pile falls below $$ 1.6^{-8}$$ J. This cutoff value was deemed sufficient after examination of the system dynamics and its kinetic energy evolution. The amplitude $$ A$$ of the tap becomes the control parameter while $$ \omega $$ remains constant, equal to $$ \omega =134$$ Hz. Although there is an open discussion (272829–30) about the proper parameter to characterize tap energy, following the most standard usage results will be reported as function of the dimensionless acceleration $$ \mathrm{\Gamma }={\omega }^2\times A/g$$, where $$ g$$ is the acceleration due to gravity.

Because of gravity, density is not uniform along the height of the system (14, 30, 31). For this reason, the analysis is performed for three narrow horizontal slices at different heights of the silo, at z = 3 to 4 cm, z = 4.5 to 5.5 cm, and z = 6 to 7 cm. To avoid wall effects, only the 2-cm-diameter inner region is considered. Then, each slice contains of the order of 3,500 particles. Apart from the quantitative differences imposed by gravity, all regions behave consistently, and, for simplicity, results corresponding only to the central slice (z = 4.5 to 5.5 cm) are shown here. Analysis corresponding to the remaining slices can be found in SI Appendix.

Local “grain densities” are obtained by dividing grain volume by their corresponding Voronoi volume, obtained using Voro++ (32) open source software. The reported density (or packing fraction) $$ \varphi $$ for the system is then obtained by averaging local grain densities of the grains over the observed region. Ten independent realizations are performed for each experiment.

Annealing Simulations.

Using a stepwise protocol inspired by the Chicago experiments, we obtain a “reversible branch” (14) that will serve as a reference for the annealing simulations. It follows a stepwise ramp-down of the tap acceleration $$ \mathrm{\Gamma }$$, initiated at a high value such as to avoid the “nonreversible” density branch (28). Starting with a loosely poured configuration at $$ \mathrm{\Gamma }=12$$, a set of 150 taps are applied at each $$ \mathrm{\Gamma }$$, followed by a decrease with a wide step, $$ -\mathrm{\Delta }\mathrm{\Gamma }$$, for another set of 150 taps, and so on. Only the density for the last 100 taps at each $$ \mathrm{\Gamma }$$ is averaged over, meant to avoid transients. Fig. 1 shows the packing fraction $$ \varphi $$ for the middle part of the system obtained in this manner as a function of $$ \mathrm{\Gamma }$$, averaged over 10 independent realizations.

Our annealing simulations, also shown in Fig. 1, employ instead a protocol of gradual changes after each tap, reducing $$ \mathrm{\Gamma }$$ minutely at a fixed $$ -\dot{\mathrm{\Gamma }}$$. In this annealing protocol, we start from the stationary state created by the previous protocol at a high value of $$ {\mathrm{\Gamma }}_{\mathrm{\text{high}}}=12$$. Then, $$ \mathrm{\Gamma }$$ is reduced by $$ -\dot{\mathrm{\Gamma }}$$ each tap until reaching well below $$ \mathrm{\Gamma }<1$$, where the density plateaus into a terminal density $$ ⟨{\varphi }_{\mathrm{\Gamma }\to 0}⟩$$ that depends on $$ -\dot{\mathrm{\Gamma }}$$. We repeat these simulations for various $$ -\dot{\mathrm{\Gamma }}$$. In Fig. 1, Inset we plot $$ ⟨{\varphi }_{\mathrm{\Gamma }\to 0}⟩$$ as a function of $$ -\dot{\mathrm{\Gamma }}$$. It varies rapidly, albeit systematically, with a sublinear exponent and an apparently much lower “ideal” equilibrium density $$ ⟨{\varphi }_{\mathrm{\Gamma }\to 0}⟩$$ for $$ -\dot{\mathrm{\Gamma }}\to 0$$, although both are difficult to determine from fitting these data.

The densities during gradual annealing as well as during the stepped Chicago protocol track each other closely for high accelerations $$ \mathrm{\Gamma }$$ and, as a function of is ramp-down with $$ -\dot{\mathrm{\Gamma }}$$, they start deviating from each other at lower $$ \mathrm{\Gamma }$$. This suggests a transition into a glassy, nonequilibrium regime. This behavior closely resembles glass transitions found in corresponding studies of polymers and many other glasses (16, 18, 2223–24), and it is consistent with previous studies finding “nonergodicity” in these systems at low values of $$ \mathrm{\Gamma }$$ (21, 33). In particular, at very low $$ -\dot{\mathrm{\Gamma }}\to 0$$, densities are reached that are above the ostensibly reversible branch, projecting an even more dense ideal protocol in the limit than is found in the Chicago pile (14, 15). (Note that we plot increasing density downward to conform with most of the previous literature on annealing of glasses, which concerns the loss of free energy, corresponding here to a loss of free volume and gravitational potential energy that amounts to the complement, $$ 1-\varphi $$, of the density.) The entire annealing behavior for a granular pile, as exhibited in Fig. 1, clearly warrants more detailed investigation in the future. Here, we merely note the existence of an actual glass-like transition.

Aging Simulations.

For the remainder of this paper we focus on the aging protocol indicated by dashed arrows in Fig. 1. The aging dynamics is induced by a hard quench ($$ -\dot{\mathrm{\Gamma }}=\infty $$) from an equilibrated state at $$ {\mathrm{\Gamma }}_{\mathrm{\text{high}}}=12$$ into the glassy regime at $$ {\mathrm{\Gamma }}_q=1.3$$, just above $$ \mathrm{\Gamma }=1$$, below which the acceleration in a tap is insufficient to excite further changes. At $$ {\mathrm{\Gamma }}_q$$, the acceleration of the pile is high enough that relaxation can proceed but sufficiently low inside the glassy regime to avoid equilibration on any experimentally reasonable timeline. We explore the aging dynamics for a sequence of $$ 2^{12}$$ taps.

In Fig. 2, we plot the evolution of the average density $$ \varphi $$ within the observation region as a function of taps. Consistent with ref. 14, $$ \varphi $$ increases logarithmically with the number of taps during aging. This is also analogous to state variables in many other glassy (thermal) systems (5, 1617–18, 34). As minute amounts of free volume and gravitational potential energy get dissipated via friction and collisions, we also find a commensurate increase in the average number of supporting contacts $$ ⟨z_c⟩$$ holding each grain in position after a tap (red squares in Fig. 2).

Fig. 2.

Average density $$ \varphi $$ of the granular pile (left axis) and average number of contacts, $$ ⟨z_c⟩$$ (right axis), as a function of number of taps, after a quench from $$ {\mathrm{\Gamma }}_{\mathrm{\text{high}}}=12$$ to $$ {\mathrm{\Gamma }}_q=1.3$$, averaged over 10 runs. Solid symbols show the log-binned average for each variable.

Identifying and counting events that facilitate such irreversible changes is far from straightforward. We define as such an event the moment at which a given grain for the first time increases its coordination number from $$ z_c\to z_c+1$$. Although its coordination number may still decrease occasionally, achieving a new record in contacts for the first time signifies an irreversible structural change in the contact network of that grain and the forces holding it in place. We argue that such an event permanently modifies the energetic landscape. Fig. 3 shows the decelerating rate at which these record-breaking events occur, essentially hyperbolically, $$ \sim 1/t$$, over four decades as a function of time $$ t$$ since the quench. Assuming that each “event” dissipates approximately an equal amount of potential energy, then the accumulated decline in potential energy—and correspondingly, the increase in density—as a logarithm in time follows quite naturally as the integral of this rate, as we will argue below.

Fig. 3.

Number of new record events per unit time (tap), where an event is accounted when a particle increases its coordination number from $$ z_c\to z_c+1$$ for the first time, averaged over 10 runs. The full symbols show the log-binned averaged.

It proves prohibitive to relate such a topological change in the contact network to an immediate, adjacent increase in density. Due to gravity, discrete topological changes in one place may often lead to a net increase in local density spread widely and much farther upstream. However, the accumulative effect of these events reveals itself microscopically in the net downward drift of all grains individually, as Fig. 4 shows. Unlike for colloidal systems in the absence of gravity, where diffusion drives relaxation and dissipation of free energy (5, 13), no discernible mean-square displacement within the horizontal plane is present in the granular pile. Instead, particles predominantly drift downward in intermittent steps, ever more slowly, as those tracks in Fig. 4 indicate. Such downward displacements signify the loss of free energy, the main mechanism by which this system relaxes and density increases.

Fig. 4.

Tracking of the displacement for a few randomly selected particles in the observed region (dashed box) inside the pile. The color coding is logarithmic in time (taps) and indicates a decelerating downward drift along increasingly rugged (intermittent) trajectories.

To examine the history of particle trajectories in more detail, we study the downward drift as a function of two times, namely the displacement achieved at $$ t$$ taps relative to the height attained after a number of $$ t_w$$ taps beyond the quench. In equilibrium, such a measure would be translationally invariant and merely depend on the difference $$ t-t_w$$. An aging system is characterized by its history dependence via the “waiting time” $$ t_w$$. The overall decelerating dynamics is reflected in the fact that the later $$ t_w$$, the fewer particles manage to advance during the subsequent “lag time” $$ \tau =t-t_w$$. However, as Fig. 5 shows, such waiting-time dependence assumes a very specific form for the granular pile, similar to other glassy systems (13, 35, 36). There, we plot the average downward displacement for all grains in the observed region as a function of lag time $$ \tau $$ for a geometric progression of waiting times. Notably, the data collapse reasonably well when time is simply rescaled by the age of the quench process itself and plotted as function of $$ \tau /t_w=\frac{t}{t_w}-1$$, as Fig. 5, Inset shows.*

Fig. 5.

Plot of the downward displacement averaged over all grains within the observation region as a function of lag time $$ \tau =t-t_w$$, relative to each grain’s position at time $$ t_w$$ after the quench. (Inset) Collapse of the data in the main panel when plotted as function of $$ \tau /t_w$$.