Spatiotemporal chaotic unjamming and jamming in granular avalanches

We have investigated the spatiotemporal chaotic dynamics of unjamming and jamming of particles in a model experiment – a rotating drum partially filled with bidisperse disks to create avalanches. The magnitudes of the first Lyapunov vector δu(t) and velocity v(t) of particles are directly measured for the first time to yield insights into their spatial correlation Cδu,v, which is on statistical average slightly larger near the unjamming than the value near the jamming transition. These results are consistent with the recent work of Banigan et al (Nature Phys. 2013), and it is for the first time to validate their theoretical models in a real scenario. v(t) shows rich dynamics: it grows exponentially for unstable particles and keeps increasing despite stochastic interactions; after the maximum, it decays with large fluctuations. Hence the spatiotemporal chaotic dynamics of avalanche particles are entangled, causing temporal correlations of macroscopic quantities of the system. We propose a simple model for these observations.

the spatiotemporal chaotic fluctuations of velocities of individual particles. As a result, the global Lyapunov vector, its linear growth rate, and the global velocity of the system are strongly correlated temporally.

Results
We first analyzed the correlation C du,v between du(t) and v(t) (The details about the measurement of du(t) can be found in the Methods). The results are shown in Fig. 1(b-f), where it shows the spatial distributions of du(t) and v(t) at t 5 1.3 s (in c and d) and at t 5 2.4 s (in e and f) respectively. Here the avalanche starts around 0.8 s (the onset of the unjamming transition) and finishes around 2.7 s (the onset of the jamming transition). The stable particles are painted in blue, which are excluded in the calculation of C du,v . The correlation is defined as , where the summation is www.nature.com/scientificreports SCIENTIFIC REPORTS | 5 : 8128 | DOI: 10.1038/srep08128 over different particle i and du (or v) stands for the average of the quantity. The results of C du,v in the unjamming (and respectively the jamming) regime are plotted in the upper panel of (b) (and respectively, the lower panel of (b)). Note that the definitions of these two regimes will be discussed in detail later. Also note that du in these two regimes are computed using different ideal trajectories (See Methods for detail). In the upper panel of (b), on average C du,v gradually decreases as a function of time; whereas in the lower panel of (b), it first remains flat and then increases rapidly and finally remains flat again. Both curves show fluctuations around 0.1, which are slightly larger in the lower panel than in the upper panel in (b). These curves allow us to extrapolate the values of C du,v at the unjamming and jamming transitions. We note that there is a plateau on the C du,v curve near the unjamming or the jamming transition as seen in Fig. 1(b). This is a common feature of all runs. We first fit the data points in Fig. 1(b) using an error function of the form a * er f(b(x 2 c)) 1 d in the neighbouring regimes of the plateau as drawn using red solid lines in the figure. Here er f() is the error function. We measure the values of the correlation function C du,v around the unjamming and jamming transitions using the plateau values of the fitting. The results are summarized in Table 1 Besides the spatial correlation between du and v, we also find strong temporal correlations between the global Lyapunov vector of the system dU(t), its linear growth rate g(t), and the global velocity of the system V(t), as shown in Fig. 2. Here dU t ð Þ~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi X i du i t ð Þ 2 q describes the deviation from the ideal trajectory of the whole system in phase space at time t, g t ð Þ~d Note that du i (t) are measured differently in the unjamming and jamming regimes by referring to different ideal trajectories (See Methods for detail) such that there are two different curves of dU(t) as shown in Fig. 2(a). Here the blue solid line represents dU(t) in the unjamming regime and the blue dashed line extends the calculation to t 5 2.7 s; similarly, the black solid line represents dU(t) in the jamming regime and the black dashed line extends the calculation to t 5 0.8 s. In the above definitions, the summation is over the set of unstable/mobile particles before time t. The number of unstable/mobile particles N(t) itself changes dramatically with time, as displayed in the inset of Fig. 2(b). After t < 0.8 s, N(t) grows exponentially to its maximum value around t 5 1.5 s and then starts to decrease exponentially, with the presence of some random fluctuations. Before t < 0.8 s, there are only a few unstable/mobile particles, essentially rattlers, from time to time. We define the unjamming and jamming regimes according to the dynamics of V(t): before (respectively, after) V(t) reaches the peak it is the unjamming (respectively, jamming) regime. By definition, there are two separate curves of g(t) in the unjamming and jamming regimes, as shown in the insets of Fig. 2(c, e) respectively, where g(t) and V(t) are strongly correlated in the unjamming and jamming regimes. In addition, dU(t) and V(t) are strongly correlated in both regimes as well, as shown in the main panels of Fig. 2(d, f), where the insets plot dU versus V for dU obtained from the extended dashed curves in Fig. 2(a) just for comparison. Results of the temporal correlations of all ten independent runs have been summarized in Table 2. To quantify the linear correlations of these macroscopic variables, we fit the scattered data points using linear fit. The degree of the linear correlations are characterized by the correlation coefficient -the last number of each cell of the table. A coefficient of 1.0 or 21.0 means a perfect linear correlation between two variables. On average the correlation coefficients are around 0.9, indicating strong linear correlations between g 2 V and V 2 dU in both unjamming and jamming regimes.

Discussion
In order to understand the above results, a crucial piece of information is the velocity v i (t) of a single particle i, as shown in Fig. 3(a, b). Besides the two trivial stable regimes where v i (t) is zero, we distinguish three regimes on this curve as time evolves: (1) a rapid exponential growth regime where the particle loses its stability; (2) an inertial regime where the particle continues accelerating till reaching the maximum, accompanied with fluctuations due to random interactions between the particle and others; (3) a frictional regime where the particle loses its kinetic energy and this regime is very bumpy with a lot of fluctuations. In Fig. 3(c), we plot the number ratio n/N for particles in the above three regimes. From 0.9 s to 1.5 s the particles in the exponential and inertial regimes dominate and n/N is non-zero in the exponential regime while the system is rapidly diverging from the unstable fixed point of the unjamming transition with a continuous replenishment of fresh unstable particles. Between 1.5 s and 2 s the inertial and frictional regimes coexist. From 2 s to 2.7 s all the moving particles are in the frictional regime while the whole system is converging towards the stable fixed point of the jamming transition. In Fig. 3(d), we also plot the spatial correlation C du,v of particles in their individual and joint set of the exponential and inertial regimes.
To understand the results presented in the above Figs. (1-3), we have proposed a mean-field model, as discussed in detail in the Supplementary Material. First, this model allows us to gain some physical insights to understand the spatial correlation between du and v in the unjamming and the jamming regimes. When a particle becomes unstable, its velocity v i (t) grows rapidly. As a result, the particle deviates quickly from its original position in real space. In a comparison, at every time instant t the contribution of the ideal trajectory to the measurement of du i (t) is negligible since the particle would follow the ideal circular motion in a speed much slower than its actual speed. Therefore, ignoring the fluctuations of the moving direction under the mean-field approximation, du i (t) is approximately an integration of v i (t), which correlates strongly with v i (t) in the spatial domain in the unjamming regime. As shown in Fig. 3(c), in the jamming regime, v i (t)9s of moving particles are incrementally evolving into the frictional regime, which is much more erratic with large fluctuations, hard to be fully captured by a simple mean field model. This might be the reason that the correlation becomes slightly smaller compared with that in the unjamming Table 1 | The values of the correlation C du,v at the unjamming and jamming transitions of ten independent experimental runs as specified by the avalanche number n 0 5 1,2, …, 10. The last two columns are the mean and standard deviation of all the ten runs avalanche number n 0 with the stable particles of the neighbouring regions along the pathway of the moving particle-mainly at the downstream of the inclination-to gain stability, creating a cascade of local jamming transitions as the entire system converges rapidly to the stable fixed point. As a result, N(t) decays rapidly in an exponential form. We cannot fully explain the linear correlations between dU, g and V in  the jamming regime since the expression of v i in the jamming regime is obscured by large fluctuations. However, We postulate that the exponential decay of N(t) might be the dominant factor leading to the linear correlations between dU, g and V. As a quantitative comparison, we find that the theoretical values of dU and V before time t 1 agree with the experimental measurement reasonably well though g is a little off compared with the real data, which is more sensitive to the parameter values used in the model.    In conclusion, we have designed a novel experiment which allows us to successfully measure the first Lyapunov vectors du in the unjamming and jamming regimes of granular avalanche processes. This allows us to study the unjamming and jamming transitions from the dynamical systems theory perspective for the first time in the experiment. At the unjamming transition, when particles become unstable, the velocity v of each individual particle grows exponentially fast such that the contributions from the ideal trajectories to the Lyapunov vectors are negligible compared with the real trajectories. Hence du and v are strongly correlated in the spatial domain. When the system rapidly escapes from the unstable fixed point at the unjamming transition to converge to the stable fixed point at the jamming transition, the strong interactions of particles at the surface layers have produced large spatiotemporal chaotic fluctuations, causing a rapid exponential increase and decrease of the number of avalanche particles. As a result, there are strong fluctuations in the velocity dynamics of individual particles in the inertial and frictional regimes compared with the exponential regime, causing a weaker spatial correlation between du and v at the jamming transition. We also have observed strong temporal correlations between the global Lyapunov vector dU, its linear growth rate g, and the global velocity V, which can be explained reasonably well using a mean-field model. Compared to the recent work of Banigan et al, our results are consistent with their numerical findings, providing supporting experimental evidence for their modeling. The unique setting of our experiment provides new physical insights to explain various correlations in a simple and intuitive way. The drastic difference between two systems in terms of driving suggests the universality of the results. One important question for further study is how to connect the dynamical instability of the system, such as the Lyapunov vector, to the geometrical packing or the force structure of the system. Despite that the grain-scale instabilities and stresses are not correlated 1 , success has been achieved to predict, at least statistically, the local rearrangements of particles from the analysis of soft spots 5,9 . How to integrate the nonlinear response and the linear response of the system will be something important to investigate in the future.

Methods
The experimental techniques. Our system is essentially a rotating drum as sketched in Fig. 1(a). It is consisted of a thin cylinder of a diameter of 80 cm, with a rotation axis perpendicular to the direction of the gravity. The rotation speed is typically slow, e.g. 15 minutes per revolution. The cylinder is hollow with a gap of 8 mm between two flat circular plates made of transparent Plexiglas with surfaces coated to reduce the accumulation of electrostatic charges. Inside the cylinder, it is filled with a monolayer of photoelastic disks up to a height slightly over 2 3 of the radius of the cylinder. These disks, 6.35 mm thick and with a total number of 736, are bidisperse with a large size of 1.4 cm and a small size of 1.2 cm in diameters. The disks are randomly distributed in space with a number ratio of 151 to avoid crystallization. During the experimental run, we have observed no particle segregation. The disks are made of the PSM4 materials manufactured by Vishay. The experiment has been repeated for ten times following the identical protocols and the results of independent runs are similar. So here we will present results in two different groups -results from one randomly selected experimental run to show the details of the dynamics and the statistics of quantitative measures from all ten independent runs to emphasize the common characteristics.
The measurement of the first Lyapunov vector. In our experiment, a key physical quantity is the first Lyapunov vector. In general, the computation is rather complex involving the time reversal of the dynamical trajectory of a system in phase space. First, one prepares an initial state and then sets the system to run at time t 5 0 with zero perturbation so that the system evolves freely in phase space. This allows one to identify the ''ideal'' trajectory of the system in phase space, i.e. the trajectory of free evolution of the system under zero perturbation. Second, once the ''ideal'' trajectory is found, one can reverse the time to return the system to the initial state in phase space. Afterwards, one applies a perturbation and records the new evolution trajectory of the system to compute the difference between the new trajectory and the ''ideal'' one. The difference depends on the perturbation. One has to repetedly apply the time reversal to set the system back to the initial state, and then to apply a different perturbation and let the system evolve in order to find the first Lyapunov vector,i.e. the trajectory where the perturbation grows the fastest. This protocol is much easier to be implemented in computer simulations. However, it is a great challenge in the real granular experiment: (1) one could not have the initial state of a granular system prepared exactly the same; (2) once an unjamming or jamming transition takes place, it is impossible to time-reverse the system to the initial state since the time reversible symmetry is broken [32][33][34][35][36] ; besides, the ''ideal'' trajectory is often untraceable. Fortunately, in our system the ''ideal'' trajectory is predictable even after the system under-goes an unjamming or a jamming transition due to external perturbations. We define the unjamming transition of particles as the onset of the avalanche where some particles at the top surface layers lose rigidity and start moving. Similarly, the jamming transition is defined as the flowing particles start coming to rest so that the avalanche comes to the end. One nice property of our system is that during the unjamming and the jamming transition, the ''ideal'' trajectory of each individual particle is exactly predictable! If there were no perturbation to trigger the avalanche to take place in advance, each unjammed particle would follow a circular trajectory. Note that there is a subtle difference between the trajectory of a particle in phase space and in real space where the trajectory in phase space includes the velocity contribution. In Fig. 4(a-b), we only draw a schematic diagram of a particle's trajectories in real space; the velocity of either the ''ideal'' or the actual trajectory can be measured in a straightforward way.
In order to make progress, we have made two assumptions. The first assumption is with regard to the perturbations. External perturbations are unavoidable in a real experiment. Let's suppose that after time t 0 , the existence of some perturbations makes the real trajectory diverge from the ''ideal'' one, where the system would reach the maximum angle of repose under zero perturbations. We assume that the existence of perturbations makes the avalanche happen in advance, which is consistent with the experimental observation that the critical angle of repose has a wide range of about 10 degrees 22 . Since the duration of the avalanche is short, typically 1 to 2 seconds, and the rotation speed is slow, i.e. at 15 min per revolution, we believe that the avalanche will finish before the system could have reached the maximum angle of the repose under zero perturbation when it had followed the ''ideal'' trajectory.
The second assumption is with regard to the first Lyapunov vector. In our experiment, it is observed that the trajectories of the avalanche particles have been extremely deviated from the ideal ones. Hence we assume that under perturbations, the difference between the actual trajectory and the ''ideal'' one corresponds to exactly the first Lyapunov vector by natural selection. www.nature.com/scientificreports In contrast to the unjamming transition, where small perturbations grow exponentially, in the jamming transition, small perturbations decrease exponentially to drive the system to a stable fixed point. Therefore, if we apply a time reversal after the avalanche, each particle in the system would essentially follow a circular trajectory in the counter clockwise direction. This can be verified by reversing the rotation of the system counter-clockwise after the avalanche, where the time reversible trajectory indeed follows nicely the circular trajectories. Based on these observations, we believe that the irregular trajectories of moving particles near the end of the avalanche can be treated as the deviation from the ''ideal'' circular trajectories due to perturbations. Here we define the first Lyapunov vector near the jamming transition as the deviation of the real trajectory from the circular ideal one. The above two assumptions are also applicable following a similar argument.