From creep to flow: Granular materials under cyclic shear

When unperturbed, granular materials form stable structures that resemble the ones of other amorphous solids like metallic or colloidal glasses. Whether or not granular materials under shear have an elastic response is not known, and also the influence of particle surface roughness on the yielding transition has so far remained elusive. Here we use X-ray tomography to determine the three-dimensional microscopic dynamics of two granular systems that have different roughness and that are driven by cyclic shear. Both systems, and for all shear amplitudes Γ considered, show a cross-over from creep to diffusive dynamics, indicating that rough granular materials have no elastic response and always yield, in stark contrast to simple glasses. For the system with small roughness, we observe a clear dynamic change at Γ ≈ 0.1, accompanied by a pronounced slowing down and dynamical heterogeneity. For the large roughness system, the dynamics evolves instead continuously as a function of Γ. We rationalize this roughness dependence using the potential energy landscape of the systems: The roughness induces to this landscape a micro-corrugation with a new length scale, whose ratio over the particle size is the relevant parameter. Our results reveal the unexpected richness in relaxation mechanisms for real granular materials.


Introduction
Yielding of amorphous materials is ubiquitous, governing a wide range of phenomena like mechanical failure of metallic glasses [5], complex rheologies of soft glasses [6], or geophysical catastrophes [7].A typical amorphous solid under quasistatic simple shear will deform elastically at small strain and flow plastically beyond the yielding strain.Depending on whether the material is brittle or ductile, the stress-strain curve shows a drop or a smooth crossover to a plateau, respectively.On the particle level, yielding is a cooperative phenomenon of local plastic events [8,9], and computer simulations hint that this phenomenon shares many aspects with a first order phase transition [4,10,11].However, for real systems the precise nature of this transition has not yet been clarified since such experiments are very challenging.
In contrast to standard glasses, for granular systems there is at present no satisfactory understanding of yielding on the level of the particles.It has been suggested that sheared granular solids are marginally stable, i.e., contacts between particles will change irreversibly even under a tiny applied strain, which implies that such systems have no elastic regime [12].This view seems to be at odds with experimental results which found that the stress-strain curve of granular solids under simple shear does in fact resemble the one for a glass [13].One possibility to resolve this discrepancy is to consider a different type of driving, i.e., cyclic shear, since it allows to probe directly the reversibility of the particle trajectories and hence the presence of elasticity in the system.A number of computational studies have in fact used this setup to investigate the particle dynamics in soft-sphere jammed states and model glasses [14][15][16][17][18][19][20][21][22][23][24].These works, as well as related experiments on soft glasses [25][26][27], have indeed revealed a reversible-irreversible transition as the cyclic shear amplitude Γ increases, thus supporting the existence of an elastic behavior if Γ 0.1.However, these findings are in contrast with experimental results for cyclically sheared granular particles in three dimensions (3D) which indicate the absence of an elastic regime, since the particle trajectories are found to be irreversible, at least for the values of Γ considered, i.e., 0.07 ≤ Γ ≤ 0.26 [28,29].
Although a few experimental studies have probed the microscopic dynamics of driven 3D granular materials [30][31][32], none of them has made an explicit connection between this dynamics and the phenomenon of yielding.Establishing this link is important not only for understanding the plasticity of amorphous materials but also for developing reliable granular constitution laws [33].The goal of the present work is therefore to identify this connection.

Cyclic shear experiment
In this work, we experimentally investigate a cyclically sheared granular system using a X-ray tomography technique [28,29,34,35].The system contains ≈14000 50 : 50 bidisperse plastic spherical beads of diameters 5 mm and 6 mm and no crystallization was detected.Particles are placed in a shear cell and cyclically sheared with an amplitude Γ [28,36].The size of the shear cell at shear γ = 0 is 24d × 24d × 24d, where d is the diameter of the small beads, and which in the following will be taken as the unit length.The shear rate is small so that we are in quasistatic conditions with an inertia number less than 10 −3 [37].
Beads are initially placed in the cell, forming a reproducible loose packings, and then compacted by cyclic shear until the steady state is reached (see Extended Data Fig. 1).X-ray tomography scans are taken at γ = 0, with a periodicity between 1 and 10 cycles, from which we extract the microscopic structure of the system and the dynamics of the particles.To improve the statistics of the results, we use 3 ∼ 5 independent realizations for each Γ, and to mitigate finite size effects we exclude the particles located closer than 6d from the boundaries.

Steady-state dynamics
For small Γ the dynamics of the particles is basically isotropic although weak convection caused by gravity emerges for large Γ (see Extended Data Fig. 2).Hence we focus here on the dynamics in the horizontal x-and y-directions.Figure 1(a) presents the mean squared displacement (MSD), δh 2 (∆n) , as a function of the number of cycles ∆n, where δh is the horizontal displacement and .stands for the average over different particles, starting configurations, and realizations.At large ∆n we find the expected diffusive growth while at small ∆n the MSD shows a power-law dependence, with a Γ-independent exponent ≈ 0.65.Such a universal creep dynamics at small ∆n demonstrates that our system has no caging or elastic regime, in agreement with earlier findings [28,29,38].This is also confirmed by the Bulk dynamics of cyclic shear for different shear amplitudes Γ.(a) Mean squared displacement in horizontal directions vs. shear cycle number ∆n.A universal crossover from sub-diffusion MSD ∝ ∆n 0.65 to normal diffusion (dashed lines) is observed, corresponding to the Γ-dependent yielding point.(b) Schematics of how a granular system explores its potential energy landscape.Due to the particle surface roughness, the PEL has not only a structure on the particle scale d, but also a micro-corrugation.Starting in the left metabasin of the PEL, the double arrows indicate the back and forth motions of the system during a shear cycle.Depending on Γ, the system does/(does not) leave the MB during a cycle (blue and orange arrows, respectively).Yielding means that the system has overcome permanently the PEL barrier (on the scale d).(c) Γ-dependence of the diffusion coefficient D (circles), the yielding strain ∆γ c (diamonds), and the strain ∆γ M at which memory is lost (triangles).Note that ∆γ c ≈ 0.4∆γ M for all Γ.Error bars represent the standard deviations from 3 − 5 independent realizations.(d) Memory M , defined in Eq. ( 1), as a function of ∆γ, from which ∆γ M is obtained by M (∆γ M ) = 0.For Γ ≤ 0.175, ∆γ M is estimated from a linear extrapolation of M (∆γ).Color codes are the same as in Fig. 1(a).
absence of a two-step relaxation in the self-intermediate scattering function (see Extended Data Fig. 3).
The absence of caging is related to the fact that granular particles have a rough surface which engenders to the potential energy landscape (PEL) of the system a micro-corrugation on a length scale that is much smaller than d, absent in atomic systems or pure hardspheres, and sketched in Fig. 1(b).This corrugation permits the particles to accumulate an irreversible displacement even at the smallest Γ hence preventing caging.Starting from a local minimum in the PEL, the creep dynamics corresponds thus to the exploration of the local metabasin (MB), i.e., the set of configurations which are on the particle level very similar to the initial one.This motion, affected strongly by memory effects (see below), allows the system to slowly reach the boundary of the MB, i.e., the yielding point, at which many particles will have changed their neighbors.After yielding the memory is lost and the MSD becomes diffusive.
In the following we report the dynamics as a function of the accumulated strain, i.e., ∆γ = 4∆nΓ, which takes into account one part of the expected Γ-dependence of the dynamics.Figure 1(c) shows that the diffusion coefficient D, obtained from the Einstein relation δh 2 = 2D∆γ, is basically constant for Γ 0.1 and starts to grow sharply beyond this threshold.Also the yielding strain ∆γ c , locating the crossover from sub-diffusive to diffusive regime (see Extended Data Fig. 4), displays a maximum at Γ = 0.1.These observations indicate a change in the underlying dynamics as a function of Γ, directly linked to the way how the system (i) explores its MB and (ii) yields at large ∆n.
Figure 1(b) illustrates that for Γ 0.1 the system will stay inside the MB during many cycles before it yields.Since the approach to the MB boundary is slow if Γ 0.1, crossing the boundary will involve many particles at a time, i.e., this yielding is a highly collective process.In contrast to this, if Γ 0.1 the particles can leave the MB during each cycle, notably at maximum strain, but due to friction and particle roughness their trajectories are largely reversible.Since this reversibility is not perfect, some particles will have left the original MB and thus the MB is slowly drained during the cycling.In this case the yielding, occurring once many particles have switched their neighbors, is no longer a cooperative process since it occurs very gradually.In contrast to this, simulations for model glasses [20] or soft-sphere jammed states [15,18], in which the PEL has no micro-corrugation, found that D is only non-zero above a certain threshold in Γ, thus a Γ-dependence that is very different from the one found here.
The mentioned reversibility of the particle motion can be inferred directly from the smallness of the MSD at ∆n = 1 in Fig. 1(a), indicating that the dynamics is not Markovian but has instead a significant memory [39].To quantify this memory we consider the correlation function between the displacements of particle i at two consecutive intervals of length ∆γ: where δx i (t) = x i (t) − x i (t − ∆γ). Figure 1(d) demonstrates that M is significantly positive at small ∆γ, i.e., the displacements are anti-correlated, and M decays with a Γ-dependent rate.For large ∆γ the memory vanishes, i.e., after yielding the dynamics becomes Markovian.We determine the strain ∆γ M at which the memory is lost, M (∆γ M ) = 0, and find that ∆γ M (Γ) tracks ∆γ c (Γ) very well, i.e., ∆γ c ≈ 0.4∆γ M , as shown in Fig. 1(c).Hence the creep dynamics is accompanied by strong memory and once the initial memory has been halved the system yields.

Yielding as a phase transition
The two distinct regimes in the yielding dynamics indicate that the nature of the corresponding dynamic phase transition, studied in simulations with simple shear [4,11], changes with Γ.One standard approach to probe the properties of a phase transition is to use as an order parameter the overlap Q(∆γ) ∈ [0, 1], which measures the similarity of two config- urations separated by ∆γ (see Methods).One expects that Q decreases with increasing ∆γ and the distribution P (Q) allows to identify the nature of the phase transition.Fig- ures 2(a-c) show that P (Q) are close to 1 for small ∆γ/∆γ c , i.e., most particles have not yet moved significantly.With increasing ∆γ/∆γ c , P (Q) shifts to the left before converging to the random distribution.For all ∆γ considered, the width of P (Q) is significantly larger for Γ ≤ 0.1 than for Γ = 0.2.This is quantified in Fig. 2(d) by plotting P (Q) at ∆γ = ∆γ c for different Γ, resulting in two master curves: A wide one for Γ ≤ 0.1 and a narrow one for Γ > 0.1.Such a Γ-dependence is supported by the inset of Fig. 2(d), which demonstrates that the standard deviation of P (Q) strongly drops at Γ ≈ 0.1.
While the P (Q) documented in Ref. [4] showed the double peak structure signalling a first order transition in which the two phases of the system can coexist, here we do not see clear evidence for this.However, it is well known that finite size effects do smear out such a structure, leaving only a single broad peak [4], which might explain the absence of a double peak due to the moderate particle number that we consider.Alternatively one can argue that the free energy barrier between the two phases is small since the micro-corrugation of the PEL allows the existence of many locally stable particle configurations that cannot be clearly assigned to one of the two phases.In other words, the surface tension between the two phases is small, and hence the first order transition is weak.A weak transition hints the vicinity of a critical point at which the transition becomes second order and the dynamics shows a critical slowing down.Upon a further increase of the external parameter (here Γ), the transition ceases to exist.This scenario is indeed compatible with our data in that the total strain ∆γ c at which the yielding occurs shows a maximum, Fig. 1(c), and also the ∆γ-dependence of the order parameter, Q (∆γ), shows a non-monotonic dependence on Γ, Fig. 2(e).(Also here the divergence expected at the critical point is rounded off because the dynamics blurs the two phases.)For Γ ≥ 0.1 the ∆γ-dependence of P (Q) indicates that there is no longer a phase transition but just a smooth crossover between the two phases, Fig. 2(c).Figure 2(f) summarizes the different behaviors as a function of ∆γ and Γ.

Dynamical heterogeneity
To obtain a microscopic understanding of the yielding, we investigate the dynamics on the particle level.Figures 3(a-c) show the distributions of the particle displacements in the x-or y-direction, i.e., the self-part of the Van Hove function G s (δx, ∆γ) [2], for ∆γ = 0.25∆γ c , ∆γ c , and 3∆γ c .Plotting G s for different Γ as a function of the rescaled distance d x = |δx|/ δx 2 roughly results in a master curve.In contrast to thermal systems the shape of this curve is clearly non-Gaussian even for small ∆γ (dotted lines), demonstrating the presence of particle motion that is faster than expected for a Gaussian process.
For d x < 2 the distribution can be fitted well by a Gumbel law [28] (solid lines), where λ = 0.57 characterizes the shape of the Gumbel distribution and B(λ) = 2.72 is a normalization constant.If d x > 2, G s exceeds the Gumbel law and shows an exponential decay (see insets).The deviation of this excess tail from the Gumbel law is quantified in Fig. 3(d), which presents the ratio . This deviation attains a maximum at ∆γ ≈ ∆γ c for all Γ, signaling that at yielding the distribution is widest and hence one has maximal dynamical heterogeneity.The largest excess is found for Γ ≈ 0.1, which is in line with theoretical arguments that non-Markovian processes with a significant memory give rise to a pronounced tail in the Van Hove function [40].
As a direct probe of dynamical heterogeneity, we determine the spatial arrangement of the fastest particles (top 10%) by calculating the number of particles belonging to the largest connected cluster (defined via a nearest neighbor criterion).Figure 4(c) presents this number, normalized by the total number of fast particles, f , and it peaks at ∆γ/∆γ c ≈ 1, i.e., yielding is accompanied by a maximal cooperativity, in agreement with Fig. 3(d).A random choice of 10% of the particles in the sample gives a f ≈ 0.04, well below the values we find here, demonstrating that the observed clustering is indeed significant.(See Extended Data Fig. 5 for the cluster size distribution.)Consistent with the phase diagram in Fig. 2(f), dynamical heterogeneity peaks at yielding and Γ ≈ 0.1 as displayed in Figs. 3 and 4.

Conclusions
Under simple shear conditions, yielding is associated with a drop in the stress-strain curve, signalling the transition from elastic to plastic behavior.For this type of driving, the mechanical response of granular materials will be very similar [3] since the PEL's microcorrugation is irrelevant.In contrast to this, the cyclic shear considered here permits to investigate the highly non-trivial effect of this micro-corrugation and hence to detect mechanical features of the system inaccessible in a simple shear setup.In view of the minimal ingredients needed to generate the micro-corrugation, we expect the surprising creep dynamics and the associated yielding behavior reported here to be generic features of granular materials and hence to be important in a multitude of situations, such as small tremors in geo-sciences or aging of civil engineering structures.
The details of the yielding dynamics, i.e., the nature of the phase transition, will depend on the micro-roughness, shape, as well as the friction coefficient of the particles, since all these parameters influence the micro-corrugation of the PEL and hence the system dynamics.It will thus be important to study how these quantities affect the fractional diffusion behavior of granular materials, or the strength of the memory.How this dependence can be included in the present theoretical approaches is not clear and thus remains a challenge for the future.We also mention that for simple shear it is custom to classify the yielding as either ductile or brittle.Our results show that for granular materials the nature of yielding depends strongly on the driving protocol, i.e., simple shear vs. cyclic shear, and also on the shear amplitude Γ.This dependence, which is absent in more standard disordered materials, suggests that for granular materials such a classification might not be possible.Advancing on these points will lead to a fundamental understanding of granular rheology and thus permit to develop a new holistic view of the failure of complex materials.
Extended Data Fig. 3. Self-intermediate scattering function F s (q, ∆γ).(a) F s (q, ∆γ) for q = 3.5, i.e., the first peak in static structure factor, and different Γ. Solid curves indicate a fit with a stretched exponential F s (q, ∆γ) = exp[−(∆γ/∆γ α ) βα ], where ∆γ α is the relaxation strain scale.Inset: Stretched exponent β α as a function of q.(b) ∆γ α as a function of q for different Γ.The scaling ∆γ α ∝ q −2 for small q (dashed line) indicates the Gaussian dynamics for large length scale, in agreement with the limit β α → 1 for q → 0 shown in the inset of panel (a).The colors of the curves in (a) and (b) are the same as in Fig. 1

FIG. 2 .
FIG. 2. Yielding as a phase transition.(a-c) Overlap function distribution P (Q) as a function of ∆γ/∆γ c for Γ = 0.05, 0.1, and 0.2.P (Q) shifts to smaller values as ∆γ/∆γ c grows (see legends).(d) P (Q) at the yielding point ∆γ = ∆γ c shows the presence of two master curves.Inset: The standard deviation of P (Q) versus Γ has a sharp transition at Γ ≈ 0.1.(e) Q as a function of ∆γ − ∆γ c .The decay is slowest for Γ ≈ 0.1 indicating the presence of a critical slowing down close to a critical point.Inset: Also Q versus ∆γ shows a slowing down at Γ ≈ 0.1.(f) Dynamic phase diagram.For any Γ, the system evolves from a creep (sub-diffusion) to flow (diffusion) regime as ∆γ grows and yields at ∆γ = ∆γ c .This yielding is first-order-like for Γ < 0.1 (purple solid curve) but shows only a continuous crossover for Γ > 0.1 (red dashed curve).Γ = 0.1 corresponds to the critical point (star).In (d) and (e) the color codes are the same as in Fig. 1(a).

FIG. 3 .
FIG. 3. Self-part of the Van Hove function G s for different Γ.The particle displacement is expressed in terms of the normalized distance d x = |δx|/ δx 2 .The color codes are the same as in Fig. 1(a).Panels (a)-(c) correspond to ∆γ = 0.25∆γ c , ∆γ c (Γ), and 3∆γ c (Γ), respectively.Solid and dashed curves are, respectively, the Gumbel [Eq.(2)] and Gaussian laws, showing that G s is non-Gaussian even at small d x .Insets show the ratio between G s and the Gumbel law G g .The vertical dashed lines at d x = 2 mark the distance beyond which G s is no longer described well by the Gumbel law.For d x ≤ 2 the Gumbel law is universal for all Γ.(d) Intensity of the excess exponential tail obtained by averaging G s /G g in the interval d x ∈ [3.5, 4.5].These curves show a maximum at ∆γ ≈ ∆γ c , i.e., at the yielding point G s has the most pronounced tail.

FIG. 4 .
FIG. 4. Dynamical heterogeneity as a function of ∆γ and Γ.(a) Non-Gaussian parameter α 2 as a function of ∆γ for different Γ. Solid curves are exponential fits α 2 = A • exp(−∆γ/∆γ g ) at large ∆γ.(b) The Γ-dependence of ∆γ g tracks the one of ∆γ c (Γ), i.e., ∆γ c ≈ 0.15∆γ g .Inset shows that A(Γ) mildly peaks at Γ ≈ 0.1.(c) Fraction of particles involved in the largest connected cluster (two particles are defined as connected if the center distance is smaller than 1.2 times their average diameter), among the top 10% mobile particles as a function of ∆γ/∆γ c .In (a) and (c) the color codes are the same as in Fig. 1(a).
(a) of the main text.(c) For q = 3.5 one finds that ∆γ α ∝ D −1 for different Γ (dashed line), indicating that F s (q, ∆γ) conveys the same information as the MSD presented in the main text.Error bars represent the standard deviations from different realizations.Extended Data Fig. 4. Scaling of the MSD to determine the yielding strain ∆γ c .(a) Powerlaw exponent ε = d log( δh 2 )/d log(∆γ) as a function of ∆γ/∆γ c .The yielding strain ∆γ c is located at ε = 0.825 (vertical and horizontal dashed lines).The evolution of ε can be described by ε = ε 1 + (ε 0 − ε 1 ) exp(−∆γ/∆γ c ), where ε 0 and ε 1 are fit parameters.(b) Γ-dependence of the parameters ε 0 and ε 1 .This graph confirms that the two exponents in the MSD are independent of Γ. Error bars represent the standard deviations from different realizations.Extended Data Fig. 5. Characterizing the clusters of top 10% fastest particles, i.e., the same threshold used in Fig. 4(c).Color of symbols are the same as in Fig. 1(a) of the main text.(a) Cluster size distribution P (s) for different Γ at the yielding point, i.e., ∆γ = ∆γ c .An approximately universal scaling P (s) ∝ s −2 (dashed line) is found.P (s) from the randomly chosen 10% (crosses with line) deviates strongly from this master curve.(b) The average gyration radius R g (s) as a function of s also shows a universal scaling R g (s) ∝ s 0.5 (dashed line).