Universality of slip avalanches in flowing granular matter

The search for scale-bridging relations in the deformation of amorphous materials presents a current challenge with tremendous applications in material science, engineering and geology. While generic features in the flow and microscopic dynamics support the idea of a universal scaling theory of deformation, direct microscopic evidence remains poor. Here, we provide the first measurement of internal scaling relations in the deformation of granular matter. By combining macroscopic force fluctuation measurements with internal strain imaging, we demonstrate the existence of robust scaling relations from particle-scale to macroscopic flow. We identify consistent power-law relations truncated by systematic pressure-dependent cutoff, in agreement with recent mean-field theory of slip avalanches in elasto-plastic materials, revealing the existence of a mechanical critical point. These results experimentally establish scale-bridging relations in the flow of matter, paving the way to a new universal theory of deformation.

U nifying scaling relations in the deformation of solids have been a long-standing challenge of material science and engineering. Such universal scaling relations are very attractive, as they provide a naturally scale-bridging framework connecting macroscopic stress-strain response to microscopic single-particle fluctuations in a single theory of deformation. Recent work on amorphous materials highlights generic features of the flow, mechanics and microscopic dynamics of a broad range of materials from metallic glasses to soft glasses [1][2][3][4][5] to granular [6][7][8] and crystalline [9][10][11] matter, lending credence to the idea of a universal theory of deformation 12 . While scaling approaches 13 that describe the scale-free flow of materials have been developed recently, direct experimental proof of an underlying internal scaling of strain over the full-length-scale range remains elusive. This would require linking the macroscopic applied force and its fluctuations to the microscopic internal fluctuations of plasticity on length scales from system size down to particle scale, a highly challenging task prohibitively difficult in conventional atomic solids.
Granular and soft materials offer the advantage that the internal flow field can be imaged conveniently with optical techniques, providing access to important microscopic quantities such as internal displacements and strain fields, which are hardly accessible in molecular systems. In particular, granular matter with particle sizes of the order of millimetres allows the motion of individual particles to be tracked accurately using bulk imaging techniques [14][15][16] . At the same time, macroscopic stress-strain measurements allow the applied force and its fluctuations to be monitored with exquisite time resolution. The combination of both offers a unique experimental opportunity to bridge length scales from single particle to macroscopic scale, making granular materials prime candidates to elucidate generic, scale-bridging aspects of flow.
A recent mean-field model 13 provides a universal description of material flow in terms of slip avalanches coupled by their induced internal elastic field: the material slips locally at weak spots, which are elastically coupled to other weak spots throughout the material resulting in highly correlated slip avalanches. The resulting scalebridging model offers a promising generic framework to bridge particle-scale dynamics to macroscopic stress-strain response. Recent highly sensitive force measurements on crystalline nanopillars 17 and amorphous metals 7 have indeed detected signatures of these scaling relations in the applied force fluctuations; however, the crucial internal scaling relations bridging down to microscopic displacements remain experimentally elusive, and the validity of this generic concept remains unclear. Furthermore, the universality of this approach, that is, its generality and applicability to a wide range of materials including soft and granular matter remains an open issue.
Here, we provide the first direct experimental measurement of just these scale-bridging relations in the flow of granular matter. By combining macroscopic force measurements with direct imaging of the internal strain distribution, we unveil surprisingly consistent scaling relations from macroscopic force response down to microscopic fluctuations. We demonstrate that at all scales, fluctuations exhibit consistent power-law distributions and correlations, truncated by systematic cutoff that grows with the applied confining pressure. These measurements give strong experimental evidence of the proximity of a mechanical critical point. We show that these scaling relations agree with predictions of the universal mean-field theory of slip avalanches in elasto-plastic materials 13 . These results for the first time demonstrate experimentally the existence of robust internal scaling relations in granular flow, bridging from microscopic strain to macroscopic stress and accounted for by a universal mean-field theory. Hence, these measurements lay the ground for a generic description of material flow within a universal, scale-bridging framework.

Results
Force fluctuations. To combine both macroscopic and microscopic measurements, we use a shear-cell set-up that links force measurements by built-in pressure sensors with simultaneous particle-scale internal imaging by laser sheets (Fig. 1a). This allows us to track fluctuations of the applied force, while imaging the internal strain distribution over the full-granulate volume (Methods section). We do this as a function of applied load that exerts a well-defined confining pressure on the top layer of the granulate. By tilting the side walls, we shear the granulate uniformly at constant (low) rate up to a maximum strain of 0.2, starting from a well-defined initial state. To increase statistics, we average over 10 shear cycles.
Using this set-up, we can detect pronounced force fluctuations as shown in Fig. 1b: sharp force drops follow continuous force increases, revealing internal relaxations of the granulate that release some of the applied force. We define a single relaxation event from a monotonic force drop with amplitude 410 À 2 N and short duration o10 À 1 s. Using this definition, we detect B10 4 relaxation events in a total sequence of 10 shear experiments. We define their size s from the magnitude DF of a sharp force drop (see Fig. 1b  frequency, D(s), as a function of size s reveals power-law distributions D(s)Bs À t truncated by systematic strain-rate dependent cutoff as shown in Fig. 2a. This cutoff moves to smaller size as the shear rate increases, suggesting that fluctuations may overlap and a relaxation event cannot complete before a new one starts, leading to truncation of large relaxations according to s max $ _ g À m (ref. 18). Indeed, we can collapse all data by rescaling avalanche sizes by _ g m , with m ¼ 0.5, as shown in Fig. 2a, inset. Furthermore, all curves exhibit a robust scaling exponent of tB1.5, as is demonstrated by the collapse in the inset. The scaling collapse applied here and below allows extrapolation of critical behaviour away from the critical point; this is particularly useful in finite-size systems, where correlation lengths close to the critical point would exceed the system size 17,19 .
Our measurements probe the initial loading of the granulate, where the applied force increases steadily with strain. This allows us to investigate how the fluctuations depend on the applied force magnitude. We sort fluctuations by the increasing force magnitude and plot separate distributions in Fig. 2b. Clearly, all distributions exhibit consistent power-law decay with a cutoff that grows with the applied force magnitude. To evaluate this growth of fluctuations, we compute average avalanche sizes os4 that are less affected by the poor statistics at large s. These grow monotonically with the applied force magnitude as shown in Fig. 2b, inset. We find a power-law dependence with exponent À 1 on approaching the critical force F c B56 ± 3 N, as shown by the double-logarithmic representation of os4 as a function of F À F c . This critical force is consistent with the apparent saturation value of the force at large strain as shown in Fig. 1b. Using this value of F c we can indeed collapse the force-dependent data shown in Fig. 2b by rescaling the avalanche sizes by (1 À oF4/F c ) 1/s , with sE0.6 as shown in Fig. 2c. This scaling collapse, which extrapolates the critical behaviour close to the critical point, indicates that the granular flow develops truly critical behaviour on approaching F c .
To explore this critical behaviour in more detail, we vary the rigidity of the granulate by changing the applied load. By releasing some of the confining pressure, we lower the rigidity of the granulate 20-23 , and we expect smaller stress build-up and hence smaller force fluctuations to occur 22,24 . This is indeed what we find as shown in Fig. 2d: as the confining pressure decreases, the cutoff of the power-law shifts to smaller sizes, indicating smaller fluctuations, similar to the applied shear force dependence in Fig. 2b. We can collapse all curves using the rescaling sP À m P with m P ¼ 0.4 (Fig. 2d, inset). A similar exponent is observed when we collapse the avalanche sizes ( Fig. 2b) with respect to the growing applied force, where we obtain m F B0.5 suggesting that the dependence on both shear force magnitude and confining pressure can be accounted for by similar scaling relation. This consistent critical scaling suggests a simple scaling model to account for the observed force fluctuations.
We apply recent mean-field theory of a stick-slip model of deformation 8 , where critical behaviour arises from the interplay of local slip and long-range elastic interactions. The model assumes that the material has weak spots where it slips when the stress exceeds a local failure stress. Because all weak spots are elastically coupled, a slipping weak spot can trigger other weak spots to also fail, resulting in a slip avalanche 13 . In the quasistatic limit, where the material is sheared slowly enough so that every  Figure 2 | Scaling of force fluctuations. (a) Probability distribution, D(s), as a function of force fluctuation size, s, for different shear rates (see legend) at constant confining pressure P ¼ 9.6 kPa. These data suggests truncated power laws DðsÞ $ s À t expð À s=x _ g Þ with x _ g $ _ g m . Inset shows these data collapse for t ¼ 1.5 and m ¼ 0.5. (b) Force dependence of the size distribution (see legend) at constant confining pressure P ¼ 9.6 kPa and shear rate _ g ¼ 1:82Â10 À 3 s À 1 . Inset: mean avalanche size as a function of force difference F c À F. The double-logarithmic plot indicates divergence hsip(F c À F) a with a ¼ À 1.02 on approaching F c ¼ 56 N. (c) Data collapse of the force-dependent avalanche size distribution shown in b, DðsÞs t $ s 1 À oF4 Fc 1=s , with t ¼ 1.5, F c ¼ 56 and sB0.6. (d) Size distribution D(s) for different confining pressures (see legend) at constant shear rate _ g ¼ 9:08Â10 À 4 s À 1 . Inset shows data collapse DðsÞs t $ sP À m P , with t ¼ 1.5 and m P ¼ 0.4. avalanche has time to complete, the model predicts that the probability density distribution D(s, F) of slip sizes s occurring at an applied force F follows a power-law with a force-dependent maximum size cutoff s max p(F c À F) À 1/s , hence, where, F c is a critical force, above which the material cannot sustain any load. The exponents t ¼ 3/2 and s ¼ 1/2 are detailindependent 'universal' scaling exponents. At finite shear rate, avalanche overlap leads to truncation of large avalanches according to s max / _ g À m with m ¼ 2 for steady-state flow. The predicted exponents t ¼ 3/2 and s ¼ 0.5 are indeed in good agreement with our observations (Fig. 2), while the value m ¼ 2 predicted for steady-state deviates substantially from the measured value m ¼ 0.5; such deviations are, however, expected when the experiments are not in steady state, as is the case here and in ref. 7. Yet, the increasing applied force magnitude allows us to test the predicted force dependence of slip sizes. To do so, we use D(s, F) from equation (1) to compute average slip sizes os4 ¼ R s D(s, F)ds, obtaining os(F)4p(F c À F) (t À 2)/s (ref. 7). The model thus predicts a power-law dependence with exponent (t À 2)/(s) ¼ À 1. The slope of the data in Fig. 2b (inset) is aB À 1.02, in excellent agreement with this prediction. The critical force F c ¼ 56 ± 3N is also very reasonable, given that the applied force in our measurements Fo55N. Thus, overall the agreement with this simple model is remarkable.
Internal strain imaging. The advantage of using granular matter to elucidate scale-bridging relations is that we can directly image the internal strain distribution. Using laser sheets, we visualize individual particles in the bulk of the granulate and track their motion over the entire strain cycle. We then determine, for each particle, the local strain from the displacement of the particle relative to its nearest neighbours (Methods section). A two-dimensional rendering of the most relevant strain component E xy is shown in Fig. 3a, where colour indicates the magnitude and sign of E xy , see colour bar. High shear strain concentrates in connected clusters that span the imaged volume, indicating correlated slip events accumulated over time. However, the deformation remains overall homogeneous and no shear banding occurs, as dictated by the rigid tilting walls, see Fig. 1a. We focus on the top 20% highest strain particles and follow them as a function of strain using fixed strain intervals of 4% to capture always the incremental strain. Plotting the average cluster size os cl 4 of high-strain regions as a function of applied strain (Fig. 3b) clearly reveals that the clusters grow with the applied strain. When we plot os cl 4 as a function of the applied force, we find indeed a power-law divergence on approaching the critical force F c B56 N, similar to the average avalanche size in Fig. 2b inset, as shown by the double-logarithmic representation of os cl 4 versus F À F c in Fig. 3c. This behaviour is robust within reasonable limits of the threshold strain used to define the high-strain particles. This direct correspondence to the avalanche size with the same value of the critical force gives independent microscopic evidence of the increasingly critical state on approaching F c . We further investigate the structure of these highly active clusters. We compute, for each cluster, the radius of gyration R 2 g ¼ 1 2s 2 cl P ij r i À r j À Á 2 from relative particle positions r i À r j in the cluster, and determine the scaling of R g as a function of cluster size s cl . We find that the gyration radius scales with cluster size as R g $ s 1=d f cl , with the fractal dimension d f ¼ 2.5 ± 0.1, which again remains robust on varying the threshold strain and considering up to the top 2% high-strain particles. This value lies very close to the value 2.53 of threedimensional (3D) percolation, suggesting that the accumulated avalanches span the entire system.
To get full insight into the microscopic avalanche evolution, we measure directly the internal scaling relations of strain using the spatial correlation function C Exy ðDrÞ ¼ E xy ðr þ DrÞ À E xy À Á E xy ðrÞ À E xy À Á =w 2 that correlates strains at locations separated by Dr. Here, angular brackets denote averaging over all particles and w 2 ¼ E xy À E xy À Á 2 D E , the s.d. The resulting correlation functions are anisotropic, exhibiting the longest correlation length along the shear direction, as shown in the inset of Fig. 3a. This correlation function visualizes directly the correlations between slip events that underlie the collective avalanches. Along the shear direction, the correlations decay closely to a power-law with cutoff that grows with the applied strain, see Fig. 4a. This behaviour is analogous to that of the force fluctuations, where the avalanche size grows with force magnitude (Fig. 2b).
The increase of the correlation length is even more pronounced with growing confined pressure-at the largest pressure, the correlation length grows up to the full-system size, see Fig. 4b.
Such increase of power-law cutoff is in qualitative agreement with the growth of force fluctuations in Fig. 2d. To test this quantitatively, we determine the full-correlation volume. We observe that in the other two spatial directions, correlation lengths hardly change; consequently, the growth shown in Fig. 4b reflects directly the bulk growth of the avalanches. Hence, if it is indeed the growth of these internal correlations that underlies the growth of macroscopic force fluctuations, we expect the same scaling collapse to apply to both. This is precisely what we find when we rescale Dr/d in the shear direction by P À m S (inset of Fig. 4b): excellent collapse is obtained for m S ¼ 0.4, the same exponent as for the macroscopic force fluctuations, demonstrating their common origin. We hence identify experimentally the underlying mechanism behind the scaling of force fluctuations: in the highly constrained granular material, microscopic strain fluctuations are strongly correlated, leading to power-law correlations from particle scale to system size for the largest confining pressure. The slope of this power-law correlation, lB0.8, is indeed reasonably close to the mean field value l ¼ 1 along the direction of greatest slip 25 . It is also close to the numerical value obtained in simulations on Lenard-Jones systems in quasistatic shear 5 that report l ¼ 1.18 along the direction of greatest displacement. These results demonstrate the robust internal scaling of microscopic strain underlying the macroscopic fluctuations, consistent with the mean-field model.
Surprisingly, we observe also strong time correlation of the microscopic deformation, in addition to its spatial correlation. To investigate this in detail, we determine the typical persistence time of the strain activity of a particle. We correlate, for each particle, its instantaneous strain with that at a later time to compute C E xy ðDgÞ ¼ hðE xy ðg þ DgÞ À hE xy i g Þ Á ðE xy ðgÞ À hE xy i g Þi g =w 2 with w 2 ¼ hðE xy À hE xy i g Þ 2 i g . The resulting time correlation function (Fig. 4c) approaches a power-law decay, again truncated by a pressure-dependent cutoff. Remarkably, the activity of a particle is correlated over half of the straining cycle, indicating strong memory. Such memory can arise from local shear-induced dilation 26 that facilitates successive shear events. We investigate this shear-dilation coupling by using the full-strain tensor to compute the local dilation from the normal strains according to The resulting normalized correlation coefficient between shear and dilation shows indeed significant coupling (Fig. 4c, inset), indicating that shear-dilation coupling plays an important role in the observed internal time coherence.

Discussion
In conclusion, our scale-bridging measurements provide the first direct experimental evidence of critical internal scaling relations behind the critical scaling of force fluctuations in flowing particulate matter. These scaling relations have the form of those known from the physics of critical phenomena: the correlation length grows with increasing applied force to diverge at the critical force F c , where the material yields. Hence, macroscopic critical fluctuations originate from a hierarchical internal strain distribution with diverging correlation length. While the granular system allows us to conveniently image and measure these scaling relations, they should be generic to the deformation of solid matter including both amorphous 7 and crystalline materials 27 , as the model is generic, based only on elasticity and local failure. The excellent agreement of the model predictions with our scalebridging measurements lays the ground for a new understanding of yielding and flow central to fields from engineering to material science to geology, delineating a new universal theory of plastic deformation.

Methods
Experiments. As granular system, we use polymethyl methacrylate spheres with diameter of d ¼ 1.5 mm and polydispersity of B5%. We fill the particles into a shear device with transparent, tiltable side walls with built-in pressure sensors (Fig. 1a). The shear cell has dimensions of 10 Â 10 Â 10 cm 3 , containing about 3 Â 10 5 particles. A top plate charged with additional weights is used to confine the granulate vertically, exerting a constant normal force between 10 and 100 N on the top layer of the granulate. After a fixed pre-shear protocol generating a reproducible initial packing, the granulate is sheared at a constant rate between _ g ¼ 3:6Â10 À 4 and _ g ¼ 1:8Â10 À 2 s À 1 to a total strain of g ¼ 20%. We measure the applied force at a frequency of 500 Hz with an accuracy of ±10 À 2 N, and a maximum force on each sensor of 55 N. To visualize the motion of the individual granular particles we use particles with a diameter of 4 mm, and add an index- correlations. These data suggests C Exy $ ðDrÞ À l expð À Dr=x P Þ with x P $ P m S . Inset shows data collapse obtained for l ¼ 0.8 and m S ¼ 0.4. (c) Autocorrelation of particle strain activity as a function of applied strain. Colour and symbols correspond to the confining pressures in b. These data suggests C Exy $ ðDgÞ À n with nE0.9. Inset shows correlation between the shear strain component E xy and dilatation, averaged over all particles. matching solution of water, NaI and fluorescein that makes the particles appear as dark spots on a bright background 14 . Individual particles are imaged in a 100 mm by 70 mm by 85 mm volume using laser sheets (Fig. 1a). 3D image stacks consisting of 180 sections with a separation of 0.389 mm are acquired every 70 s during the entire shear cycle.
Analysis. To reconstruct particle positions and trajectories from 3D image stacks we use algorithms based on 28 to track particle centres with an accuracy of ± 0.01 mm in the x-and y-, and ± 0.03 mm in the z direction 14 . To determine the local strain, we follow individual particle trajectories and identify the nearest neighbours of each particle as those that are in a range corresponding to the first minimum of the pair correlation function. This range corresponds to 1.5 particle diameters. We then compare the nearest neighbor vectors d n ¼ r n À r 0 , at applied strains of g and g þ Dg. Here, the vector r 0 indicates the position of the central particle and the index n refers to the nearest-neighbour particles. We then find the best affine deformation tensorẑ that transforms the nearest-neighbour vectors over the strain interval Dg by minimizing the mean square difference D 2 min g; Dg ð Þ¼ P n ðd n ðgÞ Àẑd n g À Dg ð Þ Þ 2 (ref. 29), where D 2 min represents the local non-affine deformation. The symmetric part ofẑ corresponds to the strain tensorÊ of the central particle. We concentrate our analysis on the shear direction-shear gradient (x À y) plane, taking into account the tensor component E xy , on which we base our correlation function analysis, see Fig. 4c. To correlate shear and dilation components, we take E xy and the sum of the three normal strains, DV ¼ E xx þ E yy þ E zz as measure of local dilation. We then compute correlations according to C Exy $DV ¼ E xy À E xy À Á Á DV À DV h i ð Þ =y Exy y DV , where y Exy and y DV denote s.d.'s of E xy and DV, and angular brackets denote averaging over all particles.
To determine force fluctuations from the measured force data, we identify changes with respect to the monotonically increasing force, see Fig. 1. Stress relaxation events correspond to sharp force drops; we thus take into account only monotonic events with negative force derivative and short duration o10 À 1 s. We take the magnitude of these abrupt force changes to be the size s of a relaxation event. The size probability distribution then follows a power-law D(s)Bs À t , as shown in Fig. 2a. To demonstrate the universality of the distributions, we collapse them onto master curves by scaling the abscissa with s Á _ g m and the ordinate axis with D(s) Á s t , see insets of Fig. 2a.