# Scaling description of the yielding transition in soft amorphous solids at zero temperature

See allHide authors and affiliations

Edited by David A. Weitz, Harvard University, Cambridge, MA, and approved August 21, 2014 (received for review April 8, 2014)

## Significance

Yield stress solids flow if a sufficiently large shear stress is applied. Although such materials are ubiquitous and relevant for industry, there is no accepted microscopic description of how they yield. Here we propose a scaling description of the yielding transition that relates the flow curve, the statistics of the avalanches of plasticity observed at threshold, and the density of local zones that are about to yield. Our description shares some similarity with the depinning transition that occurs when an elastic manifold is driven through a random potential, but also presents some striking differences. Numerical simulations on a simple elasto-plastic model find good agreement with our predictions.

## Abstract

Yield stress materials flow if a sufficiently large shear stress is applied. Although such materials are ubiquitous and relevant for industry, there is no accepted microscopic description of how they yield, even in the simplest situations in which temperature is negligible and in which flow inhomogeneities such as shear bands or fractures are absent. Here we propose a scaling description of the yielding transition in amorphous solids made of soft particles at zero temperature. Our description makes a connection between the Herschel–Bulkley exponent characterizing the singularity of the flow curve near the yield stress Σ_{c}, the extension and duration of the avalanches of plasticity observed at threshold, and the density *P*(*x*) of soft spots, or shear transformation zones, as a function of the stress increment *x* beyond which they yield. We argue that the critical exponents of the yielding transition may be expressed in terms of three independent exponents, θ, *d*_{f}, and *z*, characterizing, respectively, the density of soft spots, the fractal dimension of the avalanches, and their duration. Our description shares some similarity with the depinning transition that occurs when an elastic manifold is driven through a random potential, but also presents some striking differences. We test our arguments in an elasto-plastic model, an automaton model similar to those used in depinning, but with a different interaction kernel, and find satisfying agreement with our predictions in both two and three dimensions.

Many solids will flow and behave as fluids if a sufficiently large shear stress is applied. In crystals, plasticity is governed by the motion of dislocations (1, 2). In amorphous solids, there is no order, and conserved defects cannot be defined. However, as noticed by Argon (3), plasticity consists of elementary events localized in space, called shear transformations, in which a few particles rearrange. This observation supports that there are special locations in the sample, called shear transformation zones (STZs) (4), in which the system lies close to an elastic instability. Several theoretical approaches of plasticity, such as STZ theory (4) or soft glassy rheology (5), assume that such zones relax independently or are coupled to each other via an effective temperature. However, at zero temperature and small applied strain rate *r*), argued to decay as a power law of distance and to display a fourfold symmetry (21), as supported by observations (22⇓–24). This perturbation may trigger novel plastic events and lead to avalanches. However, even within this picture, the relationship between the avalanche dynamics and the singularity of the flow curves remains debated (7, 25).

It is tempting to seek progress by building a comparison between the yielding transition and the much better understood depinning transition that occurs when an elastic interface of dimension *d* is driven in a *d* + 1 random environment (2, 26). The role that transverse displacements play in depinning corresponds to the local accumulated plastic strain *F*_{c}, the velocity *v* vanishes nonanalytically *v* ∼ (*F* − *F*_{c})^{β} and the interplay between disorder and elasticity at threshold leads to broadly distributed avalanches corresponding to jerky motions of the interface. Much more is known about the depinning transition: it is a dynamical critical point characterized by two independent exponents related to avalanche extension and duration (26, 27). These exponents have been computed perturbatively with the functional renormalization group (28⇓–30) and evaluated numerically with high precision (31, 32). The comparison between these two phenomena has led to the proposition that the yielding transition is in the universality class of mean-field depinning (33, 34). However, experiments find a reological exponent β > 1 against the β ≤ 1 predicted for elastic depinning, and numerical simulations display intriguing finite size effects that differ from depinning (9, 10, 35⇓–37).

Formally, elasto-plastic models are very similar to automaton models known to capture the depinning transition well (17); the key difference lies in the interaction kernel *x* from instability, *P*(*x*), is singular with *P*(*x*) ∼ *x*^{θ}, unlike in depinning for which θ = 0. As we shall recall, this singularity naturally explains the finite size effects observed in simulations. In this article, we argue that once this key difference with depinning is taken into account, the analogy between these two phenomena is fruitful and leads to a complete scaling description of the yielding transition. In particular, we find that the Herschel–Bulkley exponent is related to avalanche extension and duration via Eq. **9** and that the avalanche statistics may be expressed in terms of three independent exponents: θ, *d*_{f}, and *z*, characterizing, respectively, the density of shear transformations, the fractal dimension of the avalanches, and their duration.

## Definition of Exponents

Several studies (see Table 2) characterized the yielding transition with several exponents, which we now recall.

### Flow Curves.

Rheological properties are singular near the yielding transition. Herschel and Bulkley (38) noticed that for many yield stress materials, *n*, such that

In contrast to depinning, one finds β > 1 in the yielding transition, as we explain below. Our analysis below focuses on the regime (Σ − Σ_{c})/Σ_{c} ≪ 1; effects not discussed here are expected to affect the flow curves at larger stresses (39, 40).

### Length Scales.

Near the yielding transition, the dynamics become more and more cooperative and are correlated on a length scale ξ:

### Avalanche Statistics.

At threshold Σ = Σ_{c}, the dynamics occur by avalanches whose size we define as *S* ≡ Δγ*L*^{d}, where Δγ is the plastic strain increment due to the avalanche and *L*^{d} is the volume of the system. The normalized avalanche distribution ρ(*S*) follows a power law:*L*, this distribution is cut off at some value *S*_{c}, where the linear extension of the avalanche is of order *L*, enabling definition of the fractal dimension *d*_{f}:

A key exponent relates length and time scales; *z* characterizes the duration *T* of an avalanche whose linear extension is *l*:

### Density of Shear Transformations.

If an amorphous solid is cut into small blocks containing several particles, one can define how much stress *x*_{i} needs to be applied to the block *i* before an instability occurs. The probability distribution *P*(*x*) is a measure of how many putative shear transformations are present in the sample (20). Near the depinning transition, a similar quantity can be defined, and in that case it is well known that *P*(*x*) ∼ *x*^{0} (26). We have argued (20) that it must be so when the interaction kernel *x*_{i} always decreases with time until *x*_{i} < 0, when the block *i* rearranges. Thus, nothing in the dynamics allows the block *i* to forecast an approaching instability, and no depletion or accumulation is expected to occur near *x*_{i} = 0. By contrast, for the yielding transition, the sign of *x*_{i} jumps both forward and backward, performing some kind of random walk. Because *x* = 0 acts as an absorbing boundary condition (as the site is stabilized by a finite amount once it yields), one expects that depletion may occur near *x* = 0 (20, 37, 41). In ref. 20, we argued that *P*(*x*) indeed must vanish at *x* = 0 if the interaction is sufficiently long range (particularly if |*r*^{d}, as is the case for the yielding transition); otherwise, the system would be unstable: a small perturbation at the origin would cause extensive rearrangements in the system. Thus, the yielding transition is affected by an additional exponent θ that does not enter the phenomenology of the depinning problem:*d* = 3 and θ ∼ 0.6 for *d* = 2 (20), as we confirm here with improved statistics.

The definitions of the relevant exponents are summarized in Table 1 and their values, as reported in the literature, in Table 2.

## Scaling Relations

We now propose several scaling relations, which essentially mirror arguments made in the context of the depinning transition (*SI Text*), with the additional feature that *P*(*x*) is singular.

### Stationarity.

Consider the application of a quasistatic strain in a system of linear size *L*, as represented in Fig. 2. The stress Σ fluctuates as the result of avalanches: an avalanche of size *S* leads to a stress drop proportional to the plastic strain Δγ ∼ *S*/*L*^{d}, of average 〈*S*〉/*L*^{d}. Using Eqs. **3** and **4** and assuming 2 > τ >1, one gets

Between plastic events, the system loads elastic energy, and stress rises by some typical amount ΔΣ. ΔΣ is limited by the next plastic event and thus is inversely proportional to the rate at which avalanches are triggered. Although one might think that if the system were twice as large a plastic event would occur twice as soon (implying ΔΣ ∼ 1/*L*^{d}), this is not the case, and ΔΣ depends on system size with a nontrivial exponent (9, 10, 35, 36). As argued in refs. 35 and 37, one expects ΔΣ to be of order *x*_{min}, the weakest site in the system. If the *x*_{i} are independent, this implies that

Supposing that in a stationary state, the average drop and jump of stress must be equal leads to *SI Text*, a similar but not identical relation also holds for the depinning transition (26, 42).

### Dynamics.

A powerful idea in the context of the depinning transition is that avalanches below threshold and flow above threshold are intimately related (26). Above threshold, the motion of the interface may be thought as consisting of several individual avalanches of spatial extension ξ, acting in parallel. We propose the same image for the yielding transition. If so, the strain rate

### Statistical Tilt Symmetry.

If flow above Σ_{c} consists of independent avalanches of size ξ, then the avalanche-induced fluctuations of stress on that length scale, δΣ, must be of order_{c}. Eq. **10** then leads to**11** may apply at the yielding transition. A similar relation holds for depinning of an interface if the elasticity is assumed to be linear, a nontrivial assumption underlying Eq. **10**. In that case, it can be derived using the so-called statistical tilt symmetry. In *SI Text*, we discuss evidence that linearity applies at the yielding transition, enabling us to use this symmetry to derive Eq. **11**.

Overall, the scaling relations Eqs. **7**, **9**, and **11**) allow to express the six exponents we have introduced in terms of three, which we choose to be θ, *d*_{f}, *z*. The corresponding relations are indicated in Table 1.

## Elasto-Plastic Model

The phenomenological description proposed above may apply to real materials with inertial or overdamped dynamics, as well as to elasto-plastic models, although the yielding transition in these situations may not lie in the same universality class (36, 43). In what follows, we test our predictions in elasto-plastic models, implemented as in ref. 20, whose details are recalled here.

We consider square (*d* = 2) and cubic (*d* = 3) lattices of unit lattice size with periodic boundary conditions, where each lattice point *i* may be viewed as the coarse-grained description of a group of particles. It is characterized by a scalar stress σ_{i}, a local yield stress *i* becomes plastic, which occurs at a rate 1/τ_{c} if the site is unstable, defined here as *σ*^{th} = 1. τ_{c} is the only time scale in the problem, and it defines our unit of time.

When plasticity occurs, the plastic strain increases locally and the stress is reduced by the same amount *L*^{d}. When desired, we model this effect by modifying the interaction kernel as follows:

In our model, the average plastic strain is defined as *x*) is the Heaviside function. Above Σ_{c}, the system will reach a steady state with a finite _{c}, however, the system may stop spontaneously. When this happens, to generate a new avalanche, we trigger the dynamics by giving very small random kicks to the system (chosen to conserve stress on every line and column) until one site becomes unstable.

This elasto-plastic model essentially is identical to the automaton models introduced in ref. 44 in the context of the depinning transition, in which the role of the plastic strain *u*_{i}. The only qualitative difference is the form of

## Numerical Estimation of Critical Exponents

### Flow Curves and Length Scales.

We first implement the extremal dynamics protocol: the average stress decreases by 1/*L*^{d} after each plastic event during avalanches, and increases again to generate a new active site at the beginning of a new avalanche. The corresponding stress–plastic strain curves shown in Fig. 2 allow us to estimate the critical stress Σ_{c} and the correlation length exponent ν from the fluctuations δΣ of Σ at different sizes:_{c}(*L*)〉 is the mean stress and δΣ(*L*) the SD at a given size *L*, and *k*_{1}, *k*_{2} are nonuniversal constants. From our data (*SI Text*), we obtain Σ_{c} = 0.5221 ± 0.0001, ν = 1.16 ± 0.04 for *d* = 2 and Σ_{c} = 0.5058 ± 0.0002, ν = 0.72 ± 0.04 for *d* = 3. These quantities also may be extracted reliably from finite strain rate measurements, as shown in *SI Text*.

We then compute the flow curve at a fixed strain rate. The stress is adjusted to keep the fraction of unstable sites fixed. The determination of the exponent β is very sensitive to the value of Σ_{c}. Using the values obtained from the previous analysis, we find β = 1.52 ± 0.05 for *d* = 2 and β = 1.38 ± 0.03 for *d* = 3, as shown in Fig. 3.

### Avalanche Statistics.

Avalanche statistics may be investigated using extremal dynamics and Fig. 4. As documented in *SI Text*, this method leads for our largest system size to τ ∼ 1.2 for *d* = 2 and τ ∼ 1.3 for *d* = 3. This measure appears to have large finite size effects, however. We find that such effects are diminished if we work instead at constant stress Σ and consider ρ(*S*, Σ) as *SI Text*, we find using this method that τ = 1.36 ± 0.03 for *d* = 2 and τ = 1.45 ± 0.05 for *d* = 3.

Next, we evaluate the fractal dimension *d*_{f} and the dynamical exponent *z* using extremal dynamics. Here, the avalanche cutoff *S*_{c} corresponds to an avalanche of linear extension ∼*L* so that for large systems, one expects *h* is some function and *H*(*x*) = *x*^{−τ}*h*(*x*). This collapse is checked in Fig. 4 and leads to *d*_{f} = 1.10 ± 0.04 for *d* = 2 and *d*_{f} = 1.50 ± 0.05 for *d* = 3 (for the collapse, we used the values of τ measured with the constant stress protocol). Error bars are estimated by considering the range of exponents for which the collapse is satisfactory. To measure *z*, we record the duration *T* of each avalanche and compute the duration distribution ρ(*T*) for different system sizes. These distributions are cut off at some *T*_{c}, corresponding to the duration of avalanches of spatial extension *L*, so that *T*_{c} ∼ *L*^{z}. As shown in the right panels of Fig. 4, we indeed find a good collapse *z* = 0.57 ± 0.03, τ′ ∼ 1.6 for *d* = 2 and *z* = 0.65 ± 0.05, τ′ ∼ 1.9 for *d* = 3.

### Density of Shear Transformations.

In elasto-plastic models, it is straightforward to access the local distance to thresholds *x*_{i} and to compute its distribution *P*(*x*) (20). Here we recall these results with improved statistics. We fix the stress at Σ_{c} and let the system evolve for a long enough time such that γ ≫ 1. The dynamics occasionally stop; at that point, we measure *P*(*x*) and average over many realizations. As shown in Fig. 5, we find θ = 0.57 ± 0.01 for *d* = 2, θ = 0.35 ± 0.01 for *d* = 3, where the error bar is from the error estimation of linear fit. Although in experiments *P*(*x*) is hard to access, the system size dependence of the average increment of stress where no plasticity occurs should be accessible, and follows *Insets*, ΔΣ is computed via extremal dynamics, leading to slightly smaller exponents θ ≃ 0.50 for *d* = 2 and θ ≃ 0.28 for *d* = 3, a difference presumably resulting from corrections to scaling.

## Comparison with Molecular Dynamics and Experiments

Although elasto-plastic models are well-suited to test theories, they make many simplifications, and thus may not fall in the universality class of real materials. One encouraging item is our estimate of θ, which is very similar to the value extracted from finite size effects in molecular dynamics (MD) simulations using overdamped dynamics, as reported in Table 2. This is consistent with our finding (20) that θ (and τ) is independent of the choice of dynamical rules in our model, which however, may dramatically affect the dynamics. Concerning the latter, our choice that the interaction is instantaneous in time, although still long range, likely will affect the exponents *z* and β. We expect that if a more realistic time-dependent interaction kernel *z* will satisfy *z* ≥ 1. According to Eq. **9**, this will lead to larger values of β, in agreement with experiments.

The scaling relations for τ and ν in Table 1 appear to be supported by MD simulations. In ref. 36 for overdamped dynamics, *d*_{f} = 0.9 and θ ∼ 0.54 for *d* = 2, whereas *d*_{f} = 1.1 and θ ∼ 0.43 for *d* = 3, leading to τ ∼ 1.2 in both 2*d* and 3*d*, which compares well with their measured value τ = 1.3 ± 0.1. In *d* = 2, all numerics (36, 51) report *d*_{f} ∼ 1, leading to ν ∼ 1, as observed in ref. 7. In *d* = 3, there is some disagreement on the value of *d*_{f}: refs. 51 and 52 report *d*_{f} ∼ 1.5 as we do in our elasto-plastic model, in disagreement with ref. 36, for which *d*_{f} < *d*/2, implying that ν < 2/*d*. It would be useful to resolve this discrepancy, because in the depinning problem, when ν < 2/*d*, another length scale enters the scaling description, which affects particularly finite size effects (44, 53). In this situation, however, we expect our scaling description to be unchanged if ν is meant to characterize the correlations of the dynamics for Σ > Σ_{c}.

## Conclusion

We have proposed a scaling description of stationary flow in soft amorphous solids, and it is interesting to reflect whether this approach can be applied to other systems. Plasticity in crystals shares many similarities with that of amorphous solids, and the far-field effect of a moving dislocation essentially is identical to the effect of a shear transformation (2). Thus, we expect that the stability argument of ref. 20 regarding the density of regions about to yield also applies in crystals, leading to a nontrivial exponent θ in that case too. Our scaling relations thus may hold in crystals, although the formation of structures such as domain walls might strongly affect the yielding transition.

Avalanches of plasticity are seen in granular materials in which particles are hard (12, 13). However, we believe that at least for overdamped systems, this behavior is only transient and that the elasto-plastic description does not apply to such materials under continuous shear. Some of us have argued that in this case, a picture based on geometry applies (54, 55), which also leads to a diverging length scale but of a different nature (56).

The scaling relations proposed here do not fix the values of the exponents, particularly that of θ. To make progress, it is tempting to seek a mean-field description of this problem that would apply beyond some critical dimension. Current mean-field models in which the interaction is random and does not decay with distance lead to θ = 1 (20, 41). However, the anisotropy is lost in this view, and the fact that θ diminishes as *d* increases when anisotropy is considered suggests that a mean-field model that includes anisotropy is needed. Such a model would be valuable in building a hydrodynamic description of flow that would apply, for example, to slow flow near walls (57, 58), a problem for which current descriptions do not include the role of anisotropy (25, 59).

## Acknowledgments

We thank Eric DeGiuli, Alaa Saade, Le Yan, Gustavo Düring, Bruno Andreotti, Mehdi Omidvar, Stephan Bless, and Magued Iskander for discussions. M.W. acknowledges support from New York University Poly Seed Fund Grant M8769, National Science Foundation (NSF) Chemical, Bioengineering, Environmental, and Transport Systems Grant 1236378, NSF DMR Grant 1105387, and Materials Research Science and Engineering Centers Program of the NSF Division of Materials Research (DMR)-0820341 for partial funding.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: mw135{at}nyu.edu.

Author contributions: J.L., E.L., A.R., and M.W. designed research, performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1406391111/-/DCSupplemental.

## References

- ↵
- ↵
- ↵
- ↵
- ↵.
- Sollich P

- ↵
- ↵.
- Lemaitre A,
- Caroli C

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Schall P,
- Weitz DA,
- Spaepen F

- ↵
- ↵
- ↵
- ↵.
- Nattermann T,
- Stepanow S,
- Tang L,
- Leschhorn H

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Lemaître A,
- Caroli C

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Lerner E,
- Düring G,
- Wyart M

- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics