Deformation and flow of amorphous solids: a review of mesoscale elastoplastic models

Deformation and flow of amorphous solids: a review of mesoscale elastoplastic models


The deformation and flow of disordered solids, such as metallic glasses and concentrated emulsions, involves swift localized rearrangements of particles that induce a long-range deformation field. To describe these heterogeneous processes, elastoplastic models handle the material as a collection of ‘mesoscopic’ blocks alternating between an elastic behavior and plastic relaxation, when they are too loaded. Plastic relaxation events redistribute stresses in the system in a very anisotropic way. We review not only the physical insight provided by these models into practical issues such as strain localization, creep and steady-state rheology, but also the fundamental questions that they address with respect to criticality at the yielding point and the statistics of avalanches of plastic events. Furthermore, we discuss connections with concurrent mean-field approaches and with related problems such as the plasticity of crystals and the depinning of an elastic line.


Frequently used notations

Macroscopic shear stress
Macroscopic yield stress
Local shear stress
Local yield stress
Shear strain
Shear rate
EPM Elastoplastic model
MD Molecular dynamics
rhs (lhs) right-hand side (left-hand side)


Figure 1: Overview of amorphous solids. From left to right, top row: cellular phone case made of metallic glass (1); toothpaste (2); mayonnaise (3); coffee foam (4); soya beans (5). Second row: a transmission electron microscopy (TEM) image of a fractured bulk metallic glass () by X. Tong et. al (Shanghai University, China); TEM image of blend (PLLA/PS) nanoparticles obtained by miniemulsion polymerization, from L. Becker Peres et al. (UFSC, Brazil); emulsion of water droplets in silicon oil observed with an optical microscope by N. Bremond (ESPCI Paris); a soap foam filmed in the lab by M. van Hecke (Leiden University, Netherlands); thin nylon cylinders of different diameters pictured with a camera, from T. Miller et al. (University of Sydney, Australia). The white scale bars are approximate. Just below, a chart of different amorphous materials, classified by the size and the damping regime of their elementary particles. At the bottom: some popular modeling approaches, arranged according to the length scales of the materials for which they were originally developed. STZ stands for the shear transformation zone theory of Langer (2008), and SGR for the soft glassy rheology theory of Sollich et al. (1997).

19th-century French Chef Marie-Antoine Carême (1842) claims that ‘mayonnaise’ comes from the French verb ’manier’ (‘to handle’), because of the continuous whipping that is required to make the mixture of egg yolk, oil, and vinegar thicken. This etymology may be erroneous, but what is certain is that the vigorous whipping of these liquid ingredients can produce a viscous substance, an emulsion consisting of oil droplets dispersed in a water-based phase. At high volume fraction of oil, mayonnaise even acquires some resistance to changes of shape, like a solid; it no longer yields to small forces, such as its own weight. Similar materials, sharing solid and liquid properties, pervade our kitchens and fridges: Chantilly cream, heaps of soya grains or rice are but a couple of examples. They also abound on our bathroom shelves (shaving foam, tooth paste, hair gel), and in the outside world (sand heaps, clay, wet concrete), see Fig. 1 for further examples. All these materials will deform, and may flow, if they are pushed hard enough, but will preserve their shape otherwise. Generically known as amorphous (or disordered) solids, they have no more in common than what the etymology implies: their structure is disordered, that is to say, deprived of regular pattern at “any” scale, as liquids, but they are nonetheless solid. So heterogeneous a categorization may make one frown, but has proven useful in framing a unified theoretical description (Barrat and de Pablo, 2007). In fact, the absence of long range order or of a perceptible microstructure makes the steady-state flow of amorphous solids simpler, and much less dependent on the preparation and previous deformation history, than that of their crystalline counterparts. A flowing amorphous material is therefore a relatively simple realization of a state of matter driven far from equilibrium by an external action, a topic of current interest in statistical physics.

A matter of clear industrial interest, the prediction of the mechanical response of such materials under loading is a challenge for Mechanical Engineering, too. This problem naturally brings in its wake many questions of fundamental physics. Obviously, it is not exactly solvable, since it involves the coupled mechanics equations of the elementary constituents of the macroscopic material; this is a many-body problem with intrinsic disorder and very few symmetries. Two paths can be considered as alternatives: (i) searching for empirical laws in the laboratory, and/or (ii) proposing approximate, coarse-grained mathematical models for the materials. The present review is a pedagogical journey along the second path.

Figure 2: Scientific position of elastoplastic modeling.

Along this route, substantial assumptions are made to simplify the problem. The prediction capability of models hinges on the accuracy of these assumptions. Following their distinct interests and objectives, different scientific communities have adopted different modeling approaches. Material scientists tend to include a large number of parameters, equations and rules, in order to reproduce different aspects of the material behavior simultaneously. Statistical physicists aspire for generality and favor minimal models, or even toy models, in which the parameter space is narrowed down to a few variables. At the interface between these approaches, “elastoplastic” models (EPM) consider an assembly of mesoscopic material volumes that alternate between an elastic regime and plastic relaxation, and interact among themselves. As simple models, they aim to describe a general phenomenology for all amorphous materials, but they may also include enough physical parameters to address material particularities, in view of potential applications. They rely on simple assumptions to connect the microsopic phenomenology to the macroscopic behavior and therefore have a central position in the endeavor to bridge scales in the field (Rodney et al., 2011). To some extent, EPM can be compared to classical lattice models of magnetic systems, which permit the exploration of a number of fundamental and practical issues, by retaining a few key features such as local exchange and long range dipolar interactions, spin dynamics, local symmetries, etc., without explicit incorporation of the more microscopic ingredients about the electronic structure.

This review aims to articulate a coherent overview of the state of the art of these EPM, starting in Sec. I with the microscopic observations that guided the coarse-graining efforts. We will discuss several possible practical implementations of coarse-grained systems of interacting elastoplastic elements, considering the possible attributes of the building blocks (Sec. II) and the more technical description of their mutual interactions (Sec. III). Section IV is then concerned with the widespread approximations of the effect of the stress fluctuations resulting from these interactions. In Sec. V and VI, we describe the understanding of the macroscopic response of amorphous solids to a shear deformation that can be gained from the study of EPM. Section VII focuses on more microscopic and statistical features, notably the temporal and spatial organization of stress fluctuations in ‘avalanches’, while Sec. VIII gives a short perspective on the much less studied phenomena of creep and aging. The review ends on a discussion of the relation between EPM and several other descriptions of mechanical response in disordered systems, in Sec. IX, and some final outlooks.

I General phenomenology

i.1 What are amorphous solids?

From a mechanical perspective, amorphous solids are neither perfect solids nor simple liquids.

Albeit solid, some of these materials are made of liquid to a large extent and appear soft. Nevertheless, at rest they preserve a solid structure, for example, challenging gravity, or offering elastic resistance to deformations, and will flow meekly only if a sufficient load is applied to them. Accordingly, in the rheology of complex fluids (Bonn et al., 2017), they are often referred to as “yield stress materials”. Foams and emulsions, that is, densely packed bubbles or droplets dispersed in a continuous liquid phase, owe their solidity to the action of surface tension, which strives to restore the equilibrium shape of their constituent bubbles or droplets upon deformation. Their elastic moduli are then approximately given by the surface tension divided by the bubble or droplet size, which may range from tens of microns to several millimeters; a few hundred Pascal would be a good order of magnitude. Colloidal glasses, on the other hand, are dense suspension of solid particles of less than a micron in size, which makes them light enough for Brownian agitation to impede sedimentation. They rely on entropic forces to maintain their reference structure; values between 10 and 100 Pa are often encountered for their shear moduli.

Poles apart from these soft solids, “hard” amorphous solids comprise oxyde or metallic glasses, as well as glassy polymers. They are typically made of much smaller particles than their soft counterparts. Indeed, very roughly speaking, the elastic moduli are inversely proportional to the size of the constituents. (Granular media, in which the elastic moduli depend on the material composing the grains and the applied pressure, are obviously an exception to this rule of thumb.) For instance, the atoms that compose the metallic or silica glasses live in the Angström scale, and these materials have very large Young moduli, of order (somewhat below for silicate glasses, sometimes above for metallic glasses). These atomic glasses are obtained from liquids when temperature is lowered below the glass transition temperature while crystallization is avoided. To do so, high cooling rates of typically are required for metallic glasses (Greer, 1995; Greer and Ma, 2007), whereas values below may be used for oxide glasses. After a certain amount of deformation, brittle materials will break without incurring significant plastic (irretrievable) deformation (the typical example would be silica glasses, although this has been questioned (Lacroix et al., 2012)), whereas ductile materials will deform plastically before breaking.

i.2 What controls the dynamics of amorphous solids?

Figure 3: Schematic macroscopic response of amorphous solids to deformation. (Left) evolution of the shear stress with the imposed shear strain , with a stress overshoot . In the event of fracture, the stress dramatically drops down. (Right) Steady-state flow curve, i.e., dependence of the steady-state shear stress on the shear rate , represented with semi-logarithmic axes. If the material shear bands, a stress plateau is generally observed.

Another distinction regards the nature of the excitations that can alter the structural configuration of the system.

Athermal systems

When the elementary constituent sizes are large enough ) to neglect Brownian effects (thermal fluctuations), the materials are said to be athermal. Dry granular packings, dense granular suspensions, foams, and emulsions (see Fig. 1) belong in this category. An external force is required to activate their dynamics and generate configurational changes. Typical protocols for externally driving the system include: shearing it by rotating the wall of a rheometer (Barnes et al., 1989), deforming it by applying pressure in a given direction, or simply making use of gravity if the material lies on a tilted plane (Coussot and Boyer, 1995). Rheometers control either the applied torque or the angular velocity of the rotating part. In the former case, the applied macroscopic shear stress is kept fixed, at a value on a rotating cylinder of radius and height (Fardin et al., 2014), while one monitors the resulting shear strain or shear rate if the material flows steadily. Conversely, strain-controlled experiments consist in imposing or and monitoring the stress response .

Now, how do amorphous solids respond to such external forces? For small applied stresses , the deformation is elastic, i.e., mostly reversible (see Fig.3[left]). Submitted to larger stresses, the material shows signs of plastic (irreversible) deformation; but the latter ceases rapidly, unless overcomes a critical threshold known as yield stress. For , the material yields. This process can culminate in macroscopic fracture; for brittle materials like silica glass, it always does so. Contrariwise, most soft amorphous solids will finally undergo a stationary plastic flow. The ensuing flow curve in the steady state is often fitted by a Herschel-Bulkley law


with (see Fig. 3[right]).

The transition between the solid-like elastic response and the irreversible plastic deformation is known as yielding transition. Statistical physicists often regard it as an example of a dynamical phase transition, an out-of-equilibrium phenomenon with characteristics similar to equilibrium phase transitions (Lin et al., 2014b, 2015; Jaiswal et al., 2016).

Thermal systems

On the other hand, thermal fluctuations may play a role in materials with small enough () elementary constituents, such as colloidal and polymeric glasses, colloidal gels, silicate and metallic glasses. Still, these materials are out of thermodynamic equilibrium and they do not sample the whole configuration space under the influence of thermal fluctuations. It follows that different preparation routes (and in particular different cooling rates) tend to produce systems with different mechanical properties. Even the waiting time between the preparation and the experiment matters, because the system’s configuration evolves meanwhile, through activated events: this is the aging phenomenon. In particular, the high cooling rates used for quenching generate a highly heterogeneous internal stress field in the material (Ballauff et al., 2013). In some regions, particles manage to rearrange geometrically, minimizing in part the interaction forces among them, but many other regions are frozen in a highly strained configuration. Slow rearrangements will take place at finite temperature and tend to relax locally strained configurations (“particles break out of the cages made by their neighbors”), along with the stress accumulated in them.

That being said, the elastic moduli are usually only weakly affected by the preparation route, i.e., the cooling rate (Ashwin et al., 2013) and the waiting time (Divoux et al., 2011b), while other key features of the transient response to the applied shear are often found to depend on it. This sensitivity to preparation particularly affects the overshoot in the stress vs. strain curve, depicted in Fig. 3 and used to define the static yield stress , and is observed in experiments (Divoux et al., 2011b) as well as numerical simulations (Rottler and Robbins, 2005; Patinet et al., 2016). In soft materials that can undergo stationary flow, this issue may be deemed secondary; the flow creates a nonequilibrium stationary state, and the memory of the initial preparation state is erased after a finite deformation. On the other hand, in systems that break at finite deformation, the amount of deformation before failure is of paramount importance, and so is its possible sensitivity to the preparation scheme, due to different abilities of the system to localize deformations (see Sec. V).

Potential Energy Landscape

The Potential Energy Landscape (PEL) picture offers an illuminating perspective to understand the changes associated with aging in thermal systems (Goldstein, 1969; Doliwa and Heuer, 2003b, a). The whole configuration of the system (particle coordinates and, possibly, velocities) is considered as a “state point” that evolves on top of a hypersurface representing the total potential energy. Despite the high dimension of such a surface (proportional to the number of particles), it can be viewed as a rugged landscape, with hills and nested valleys; the number of local minima generally grows exponentially with (Wales and Bogdan, 2006). Contrary to crystals, glassy (disordered) states do not minimize the free energy of the system; aging thus consists in an evolution towards lower-energy states (on average) through random, thermally activated jumps over energy barriers, or more precisely saddle points of the PEL. As the state point reaches deeper valleys, the jumps become rarer and rarer; the structure stabilizes, even though some plasticity is still observed locally (Ruta et al., 2012).

Under external loading, the system responds on much shorter time scales than for aging. Accordingly, some thermal systems may be treated as athermal, for all practical purposes. Nonetheless, interesting physical behavior emerges when these two time scales compete (either because the temperature is high enough, or because the system is driven very slowly) (Rottler and Robbins, 2005; Shi and Falk, 2005; Johnson and Samwer, 2005; Vandembroucq and Roux, 2011; Chattoraj et al., 2010).

i.3 Jagged stress-strain curves and localized rearrangements

The contrasting inelastic material responses to shear, ranging from failure to flow, may give the impression that there is a chasm between “hard” and “soft” materials. They are indeed often seen as different fields, plasticity for hard solids versus rheology for soft materials. Nevertheless, the gap is not so wide as it looks. Indeed, some hard solids may flow plastically to some extent without breaking, while soft solids retain prominent solid-like features under flow at low enough shear rates, unlike simple liquids.

To start with, consider the macroscopic response to a constant stress (or shear rate ) of a foam (Lauridsen et al., 2002) or a metallic glass (Wang et al., 2009): Instead of a smooth deformation, the strain evolution with time (or stress evolution , respectively) is often found to be jagged. The deeper the material lies in its solid phase, the more “serrated” the curves (Dalla Torre et al., 2010; Sun et al., 2012). Serrated curves are not specific features of the deformation of amorphous solids; they are observed in all “stick-slip” phenomena. Indeed, the system is repetitively loaded until a breaking point, where an abrupt discharge (energy release) occurs. Interestingly, this forms the basis of the elastic rebound theory proposed by Reid (1910) after the 1906 Californian earthquake. Other elementary examples include pulling a particle with a spring of finite stiffness in a periodic potential, a picture often used in crystalline solids to describe the motion of dislocations - the elementary mechanism of plasticity. In the plastic flow of amorphous solids, potential energy is accumulated in the material in the form of elastic strain, until some rupture threshold is passed. At this point, a plastic event occurs, with a release of the stored energy and a corresponding stress drop. The precise nature of the plastic event that gives rise to the stress drop, however, will strongly depend on the scale of observation.

Once again, the PEL perspective is enlightening: in this perspective, the external driving imposes a (usually time-dependent) constraint on the regions that the state point can visit. Mathematically, this is enforced by means of a Lagrange multiplier, which effectively tilts the potential into


where is the volume of the system and the macroscopic stress. As the imposed macroscopic strain increases, some of the barriers surrounding the state point subside, until one of them flattens so much that the system can slide into another valley without energy cost. This marks the onset of a plastic event; for smooth potentials, close to the topological change at , Gagnon et al. (2001); Maloney and Lacks (2006) demonstrated that the barrier height scales as


Note that this instability can be anticipated if thermal fluctuations are present. In summary, in the PEL, deformation is a succession of barrier-climbing phases (elastic loading) and descents. (For a discussion on the properties of the PEL of a model glass, see, e.g., (Doliwa and Heuer, 2003a)) . The first step in building a microscopic understanding of the flow process is to identify the nature of these plastic events.

But what can be said about the microscopic deformation of atomic or molecular glasses when the motion of atoms or molecules remains virtually invisible to direct experimental techniques? In the 1970s, inspiration was brought by the better known realm of crystals. As early as 1934, with the works of Orowan, Polanyi and G. I. Taylor, it was known that the motion of crystalline defects (dislocations) is the main lever of their (jerky) deformation. Could similar static structural defects be identified in the absence of a regular structure? The question has been vivid to the present day, so that it is fair to say that, at least, they are are much more elusive than in crystals. But the main inspiration drawn from research on crystals was in fact of more pragmatic nature: Bragg and Nye (1947) showed that “bubble rafts”, i.e., monolayers of bubbles, could serve as upscaled models of crystalline metals and provide insight into the structure of the latter. The lesson was simple: If particles in crystals are too small to be seen, let us make them larger. Some thirty years later, the idea was translated to disordered systems by Argon and Kuo (1979), who used bidisperse bubble rafts as models for metallic glasses. Most importantly, they observed the prominence of singularities in the deformation, more precisely, rapid rearrangements involving a few bubbles. Princen and Kiss (1986) suggested that the elementary rearrangement in such systems was a local topological change of the foam structure known as T1 event and involving four bubbles (see Fig. 4a).

Figure 4: Localized rearrangements. (a) T1 event in a strained bubble cluster. From (Biance et al., 2009). (b) Sketch of a rearrangement. From (Bocquet et al., 2009). (c) Instantaneous changes of neighbors in a slowly sheared colloidal glass. Adapted from (Schall et al., 2007). Particles are magnified and colored according to the number of nearest neighbors that they lose.


Since then, evidence for such swift localized rearrangements has been amassed in very diverse systems, both experimentally and numerically. Here, we shall simply list some of the works that followed Argon and Kuo (1979)’s and early investigations on foams and emulsions (Princen, 1983, 1985; Princen and Kiss, 1986, 1989):

- simple numerical glass models like Lennard-Jones glasses (Falk and Langer, 1998; Maloney and Lemaître, 2004, 2006; Tanguy et al., 2006) and other systems (Gartner and Lerner, 2016),

- numerical models of metallic glasses (Srolovitz et al., 1983; Rodney and Schuh, 2009),

- numerical models of silicon glasses (amorphous silicon) (Fusco et al., 2014; Albaret et al., 2016),

- numerical models of polymer glasses (Smessaert and Rottler, 2013; Papakonstantopoulos et al., 2008)

- dense colloidal suspensions (Schall et al., 2007; Chikkadi and Schall, 2012; Jensen et al., 2014),

- dense emulsions (Desmond and Weeks, 2015),

- dry and wet foams (Debregeas et al., 2001; Kabla and Debrégeas, 2003; Biance et al., 2009, 2011),

- granular matter (Amon et al., 2012, 2013; Le Bouil et al., 2014; Denisov et al., 2016).

Admittedly, the details of these rearrangements do vary between the systems (see below). But, in all cases, they are the essential events whereby the macroscopic deformation is transcribed into the material structure, beyond the elastic response; their essential characteristic - as compared to the crystalline case - is their strong spatial localization. In the following, we shall refer to these events, which are the building bricks of EPM, as “plastic events”1. Since these rearrangements must contribute to the deformation, they will retain part of the symmetry of the externally imposed shear and can thus be modeled as localized shear deformations (or “shear transformations”), if variations are overlooked, whether correctly or not.

Quantitative description

Although these rearrangements can sometimes be spotted visually, a more objective and quantitative criterion for their detection is desirable. Making use of the inelastic nature of these transformations, Falk and Langer (1998) pioneered the use of , a quantity that measures how non-affine the local displacements around a particle are. More precisely, the relative displacement of neighboring particles between successive configurations is computed, and compared to the one that would result from a locally affine deformation; is the deviation obtained by optimizing the parameters of the local affine deformation to minimize the deviation. This quantity has been used heavily since then (Schall et al., 2007; Chikkadi et al., 2011; Chikkadi and Schall, 2012; Jensen et al., 2014), Generally speaking, a very strong localization of events is observed, with spatial maps of that consist of a few active regions of limited spatial extension, separated by regions of (locally) affine and elastic deformation.

Other indicators of nonaffine transformations have also been used. For instance, different observables, including the strain component along a neutral direction (say, if the applied strain is along in a 3D system) (Schall et al., 2007), the field of deviations from the uniform transformation of particle positions, the count of nearest-neighbor losses (Chikkadi and Schall, 2012) or the identification of regions with large (marginal) particle velocities, are also good options to detect rearrangements. Up to differences in their intensities, these methods were shown to provide similar information about shear transformations in slow flows of colloidal suspensions (Chikkadi and Schall, 2012). Alternative methods take advantage of the irreversibility of plastic rearrangements, by reverting every strain increment imposed on the system () in a quasistatic shear protocol and comparing the reverted configuration with the original one (Albaret et al., 2016). Differences will be seen in the rearrangement cores (which underwent plasticity) and their surroundings (which were elastically deformed by the former). To specifically target the anharmonic forces active in the core, shear can be reverted partially, to harmonic order, by following the Hessian upstream instead of performing a full nonlinear shear reversal (Lemaître, 2015).

Some reservations should now be made with respect to the picture of clearly separated localized transformations. First, the validity of the binary picture distinguishing elastic and plastic regions has been challenged for hard particles, such as grains (Bouzid et al., 2015a). More generally, near the jamming transition, the complexity of particle motion and the spatial extent of low-energy vibration modes may jeopardize the accuracy of this vision (Andreotti et al., 2012). It is also clear that as the temperature or the shear rate are increased and the material departs from solidity, thermal or mechanical noise may wash out the picture of well isolated, localized events. Nevertheless, it has recently been argued that localized rearrangements can still be identified at relatively high temperatures. For instance, these rearrangements leave an elastic imprint in supercooled liquids via the elastic field that each of them induce; this imprint is revealed when one studies suitable stress or strain correlation functions (Chattoraj and Lemaître, 2013; Lemaître, 2014; Illing et al., 2016).


These quantitative indicators of microscale plasticity have brought to light substantial variations and differences between actual rearrangements and idealized shear transformations. Even though EPM will generally turn a blind eye to this variability, let us shortly mention some of its salient features:

First, the size of plastically rearranging regions varies from a handful of particles in foams, emulsions and colloidal suspensions (for instance, about 4 particles in a sheared colloidal glass, according to Schall et al. (2007)) to a couple of dozens or hundreds in metallic glasses (10 to 30 in the numerical simulations of Fan et al. (2015), 25 for the as-cast glass and 34 for its annealed counterpart in the indentation experiment of Choi et al. (2012), 200 to 700 in the shearing experiments of Pan et al. (2008)). Note that, for metallic glasses, the indicated sizes are not backed out by direct experimental evidence, but are based on activation energy calculations.

Albaret et al. (2016) proposed a detailed numerical characterization of plastic rearrangements in their amorphous silicon model by fitting the particle displacements during plastic events with the expected (Eshelby) elastic fields around ellipsoidal transformations zones with spontaneous deformation (where the tensor was fitted). Although rearrangements seem to have a typical linear size, around , they found that the most robust quantity is actually the product of with the inclusion volume . Furthermore, the diagonal components of (dilation or compression) only represent about 5% of the deviatoric components (shear), which confirms the prevalent shear nature of the transformation. It should also be mentioned that the diagonal components were either of positive or negative, i.e., either of dilational or of compressional nature depending on the specificities of the implemented potential: Plastic rearrangements are not always dilational. Finally, the authors of the study were able to reproduce the stress vs. strain curve on the basis of the (strain-dependent) shear modulus and the fitted local elastic strain releases . This confirms that localized plastic rearrangements emitting an Eshelby field are the unique elementary blocks of the plastic response.

Secondly, the shape of the rearrangements is also subject to variations. In quiescent systems rearrangements through string-like motion of particles seem to be more accessible (Keys et al., 2011), even though shear transformations have also been claimed to be at the core of structural relaxation in deeply supercooled liquids (Lemaître, 2014). The application of a macroscopic shear clearly favors the latter type of rearrangements. Albeit facilitated by the driving, these transformations may nonetheless be predominantly activated by thermal fluctuations in thermal systems (Schall et al., 2007). There is some (limited) indication that the characteristics of the rearranging regions change as one transits from mechanically triggered events to thermally activated ones, for instance with a visible increase in the size of the region in metallic glass models (Cao et al., 2013).

Thirdly, owing to the granularity of the rearranging region (which is not a continuum!), the displacements of the individual particles in the region do not strictly coincide with a shear transformation, i.e., (where generically refers to a particle position); incidentally, this is the major reason why the observable detects plastic rearrangements.

Structural origins of rearrangements

Figure 5: Local structural determinants for the onset of a rearrangement. (a) Contour maps of the particle-based participation ratio in the 1% softest vibrational modes for two numerical samples of metallic glass (Ding et al., 2014). (b) Correlation between the locations of future plastic rearrangements and diverse local properties in an instantaneously quenched binary glass model. The following properties have been considered: local yield stress (), participation in the soft vibrational modes (PF), lowest shear modulus (), local potential energy (PE), short-range order (SRO), local density (). From (Patinet et al., 2016). (c) Snapshot of the configuration of a compressed granular pillar, with particles colored from gray to red according to their value. Particles identified as soft by the SVM have thick black contours. From (Cubuk et al., 2015).

What determines a region’s propensity to undergo a rearrangement? Microstructural properties underpinning the weakness of a region (i.e., how prone to rearranging it is) have long been searched. In the first half of the 20th century efforts were made to connect viscosity with the available free volume per particle, notably by using (contested) experimental evidence from polymeric materials (Batschinski, 1913; Fox Jr and Flory, 1950; Doolittle, 1951; Williams et al., 1955). The idea that local variations of free volume control the local weakness have then been applied widely to systems of hard particles (metallic glasses, colloidal suspensions, granular materials) (Spaepen, 1977). Falk and Langer (1998)’s Shear Transformation Zone theory originally proposed to distinguish weak zones prone to shear transformations on the basis of this criterion. Hassani et al. (2016) have invalidated criteria based on the strictly local free volume but showed that a nonlocal definition distinctly correlates with the deformation field, as do potential-energy based criteria (Shi et al., 2007). Paying closer attention to the microstructure, Ding et al. (2014) proved the existence of correlations between rearrangements and geometrically unfavored local configurations (whose Voronoi cell strongly differs from an icosahedra) in model binary alloys.

Linear properties have also been considered, with the hope that regions that are soft in terms of their linear response would also be weak in their nonlinear response. Regions with low elastic shear moduli were indeed shown to concentrate most of the plastic activity (Tsamados et al., 2009), even though no yielding criterion based on the local stress or strain is valid uniformly throughout the material (Tsamados et al., 2008). Focusing on vibrational properties, Widmer-Cooper et al. (2008) provided evidence that in supercooled liquids the particles that vibrate most in the lowest energy modes (i.e., those with a high participation fraction in the softest modes, where is arbitrarily fine-tuned), are more likely to rearrange. This holds true at zero temperature (Manning and Liu, 2011) and also for metallic and polymer glasses (Smessaert and Rottler, 2015) (see Fig. 5). Note that this enhanced likelihood should be understood as a statistical correlation, beating random guesses by a factor of 2 or 3 or up to 7 in some cases, rather than as a systematic criterion. However, in the cases where the rearrangement spot is correctly predicted, the soft-mode-based prediction for the direction of motion during the rearrangement is fairly reliable (Rottler et al., 2014).

If one is allowed to probe nonlinear local properties, then Patinet et al. (2016) showed that predictions based on the local yield stress, numerically measured by deforming the outer medium affinely, outperform criteria relying on the microstructure and the linear properties, as indicated in Fig. 5c.

Leaving behind traditional approaches, a couple of recent papers showed that it is possible to train an algorithm to recognize the atomic-scale patterns characteristic of a glassy state and spot its ”soft” regions  (Cubuk et al., 2015; Schoenholz et al., 2016; Cubuk et al., 2016). In this Machine Learning approach, rather than focusing on typical structure indicators, a large number of ‘features’ quantified for each particle is used, concretely M=166 ‘structure functions’, indicating e.g. the radial and angular correlations between an atom and its neighbors (Behler and Parrinello, 2007). Adopting both an experimental frictional granular packing and a bidisperse glass model, the authors focused on the identification of local softness and its relation with the dynamics of the glass transition. First, with computationally costly shear simulations and measurements of nonaffine displacements via , the particles that ‘move’ (i.e., break out of the cages formed by their neighbors) are identified as participating in a plastic rearrangement and used to train a Support Vector Machine (SVM) algorithm. Each particle’s environment is handled as a point in the high-dimensional vector space parameterized by the structure functions and the algorithm identifies the hyperplane that best separates environments associated with ‘moving’ particles and those associated with ‘stuck ones’ in the training set. Once trained, the algorithm is able to predict with high accuracy if a particle will ‘move’ or not when the material is strained, depending on its environment in the quiescent configuration, prior to shear.

i.4 Nonlocal effects

Once a rearrangement is triggered, it will deform the medium over long distances, in the same way as an earthquake is felt a large distance away from its epicenter. This may trigger other rearrangements at a distance, which rationalizes the presence of nonlocal effects in the flow of disordered solids. Importantly, this mechanism relies on the solidity of the medium, which is key to the transmission of elastic waves.

These long-range interactions and the avalanches that they may generate justify the somewhat hasty connection sketched above between the serrated macroscopic stress curves and the abrupt localized events at the microscale. The problem is that in the thermodynamic limit any one of these micro-events should go unnoticed macroscopically. For sure, the thermodynamic limit is not reached in some materials, notably those with large constituents, such as foams and grains, but also in nanoscale experiments on metallic glasses and numerical simulations. On the other hand, in the bulk, without collective effects and avalanches involving a large number of plastic events, the impact of microscopic events on the macroscopic response could not be explained. Since mesoscale plasticity models intend to capture these collective effects, a description of the interactions at play is required.

Idealized elastic propagator

Let us start by focusing on the consequence of a single shear transformation. Its rotational part can be overlooked because its effect is negligible in the far field, as compared to deformation, represented by the linear strain tensor , where stands for the displacement. Recall that a shear deformation, say in two dimensions, consists of a stretch along the direction , in polar coordinates, and a contraction along the perpendicular direction. The induced displacement field simply mirrors this symmetry, with displacements that point outwards along and inwards along . This leads to a dipolar azimuthal dependence for and a four-fold (‘quadrupolar’) one for its symmetrized gradient . More precisely, by imposing mechanical equilibrium on the stress , viz.,

in an incompressible medium () with a linear elastic law, (where the superscript denotes the deviatoric part), Picard et al. (2004) derived the induced strain field in two dimensions,


Here, only one of the strain components is expressed, but the derivation is straightforwardly extended to a tensorial form (Nicolas and Barrat, 2013a; Budrikis et al., 2017). Experiments on colloidal suspensions (Schall et al., 2007; Jensen et al., 2014) and emulsions (Desmond and Weeks, 2015) as well as numerical works (Kabla and Debrégeas, 2003; Maloney and Lemaître, 2006; Tanguy et al., 2006) have confirmed the relevance of Eq. 4, as illustrated in Fig. 6.

Exact induced field and variations

The strain field of Eq. 4 is valid in the far field, or for a strictly pointwise shear transformation. Yet, the response can be calculated in the near field following Eshelby (1957)’s works, by modeling the shear transformation as an elastic inclusion bearing an eigenstrain , i.e., spontaneously evolving towards the deformed configuration . This handling adds near-field corrections to Eq. 4.

Describing a plastic rearrangement with an elastic eigenstrain is imperfect in principle, but the difference mostly affects the dynamics of stress relaxation (Nicolas and Barrat, 2013a). In fact, Eshelby’s expression perfectly reproduces the average displacement field induced by an ideal circular shear transformation in a 2D binary Lennard-Jones glass (Puosi et al., 2014), although significant fluctuations around this mean response arise because of elastic heterogeneities. The numerical study was extended to the deformation of a spherical inclusion in 3D, and to the nonlinear regime, by Priezjev (2015).

Besides elastic heterogeneity, further deviations from the Eshelby response result from the difference between an actual plastic rearrangement and the idealized shear transformation considered here. For instance, Cao et al. (2013) report differences between the medium or far-field response to rearrangements in the shear-driven regime as opposed to the thermal regime and that only the former quantitatively obey Eshelby’s formula. It might be that the dilationnal component of the rearrangement, discarded in the ideal shear transformation, is important in the thermal regime.

Figure 6: Average stress redistribution around a shear transformation. (a) Experimental measurement in very dense emulsions. Adapted from (Desmond and Weeks, 2015). (b) Average response to an imposed shear transformation obtained in atomistic simulations with the binary Lennard-Jones glass used by Puosi et al. (2014). (c) Simplified theoretical form, given by Eq. 4. From (Martens et al., 2012). Note that the absolute values are not directly comparable between the graphs and that, in subfigures b and c, the central block is artificially colored.

The salient points discussed above in the rheology of amorphous solids seem to build a coherent scenario, consisting of periods of elastic loading interspersed with swift localized rearrangements of particles. These plastic events may interact via the long-range anisotropic elastic deformations that they induce. These elements are the phenomenological cornerstones of the EPM described in the following section.

Ii The building blocks of elastoplastic models (EPM)

ii.1 General philosophy of the models

The simplicity and genericity of the basic flow scenario described above has led to the emergence of multiple, largely phenomenological, coarse-grained models. These models are generally described as “elasto-plastic” or sometimes “discrete automaton” or “mesoscopic” models for amorphous plasticity. They incorporate the basic flow scenario by decomposing the system into “mesoscopic” blocks (presumably of the typical size of a rearrangement) in which the elastic behavior is interrupted by plastic events. With a few exceptions, they are implemented on a regular lattice, so they are effectively a subclass of discrete automata evolving according to predefined rules. Schematically, a model can be specified by the following set of rules (Rodney et al., 2011):


a (default) linear elastic response of each mesoscopic block,


a local yield criterion that determines the onset of a plastic event (),


a redistribution of the stress during plasticity that gives rise to long range interactions among blocks,


a recovery criterion that fixes the end of a plastic event (),

where the activity is defined as if the block is purely elastic, and otherwise.

Rules R2 and R4 determine the dynamical rules controlling the manner in which a region switches from elastic to plastic and conversely, for instance by specifying the rates associated with the transitions

whereas R1 and R3 define the mechanical response of the material for a given set of plastically active regions (ongoing plastic events). This is specified by an equation of evolution of the stress carried by block i (where is a d-dimensional vector denoting the lattice coordinates of the block); to fix ideas, if the sample is subjected to simple shear at a rate , this equation might read, in scalar version,


Here, the stress increment per unit time is the sum of an elastic contribution from the external driving, obeying Hooke’s law (viz., with the shear modulus), a local plastic relaxation if the site is currently plastic (i.e., if ), and a redistribution of the stress released by nonlocal plastic events (at positions j), mediated by the propagator . Locally, the plastic deformation of an active block () is constrained by the elastic deformation of the embedding medium, hence a relaxation efficiency . Note that, in translationally invariant geometry, the propagator becomes independent of space, viz., and .

In essence, EPM aspire to follow in the footsteps of the successes of simplified lattice models in describing complex collective phenomena in condensed matter and statistical physics. The assumption is that most microscopic details are irrelevant with respect to the main rheological properties and that the physics can be condensed into a few relevant parameters or processes. Several reasons could be put forward to favor their use over more realistic modeling approaches, e.g.,

  1. to assess the validity of a proposed theoretical scenario and, in particular, to identify the key physical processes in the rheology,

  2. to provide an efficient simulation tool giving access to (otherwise inacessible) large statistics or long-time runs,

  3. to provide a simple route towards the derivation of macroscopic equations and to bridge the gap between rheological models (constitutive laws) and statistical physics models (sandpile models, depinning models, Ising-like models).

The substantial variations in the physical ingredients incorporated in distinct EPM, notably with respect to R2 (yield criterion) and R4 (plastic event duration), seem to be a strong blow to the first objective. But one should bear in mind that these materials are so diverse that a given process (e.g., thermal activation) may be negligible in some of them, and paramount in others. Perhaps less intuitive is the role played by the experimental conditions and the observables under consideration in determining the physical ingredients that need to be retained in an EPM. Let us bolster this statement with a couple of examples. Keeping track of previous configurations (e.g., past yield stresses) might be vital in order to describe oscillatory shear experiments in which the system performs a cycle between a few configurations (Fiocco et al., 2014), whereas it is irrelevant for steady shear. Also, potentially universal aspects of the yielding transition are expected to be relatively insensitive to the precise rules, while the details of the flow pattern will obviously be more affected. Thus, as noted in (Bonn et al., 2017), one should not only select the relevant ingredients in a model only in light of the intrinsic importance of these effects (as quantified for instance by dimensionless numbers), but also depending on their bearing on the investigated properties.

In the following, we list the physical processes that are put in the limelight in the diverse EPM and indicate for what type of materials and in what conditions they are of primary importance.

ii.2 Thermal fluctuations

How relevant are thermal fluctuations and their effect on the motion of particles in the description? This question boils down to the distinction between thermal materials and athermal ones exposed in Section I.2.

It is widely believed that thermal fluctuations largely contribute to the activation of plastic events in metallic and molecular glasses, as well as in colloidal systems made of small enough colloids. In the latter systems, Schall et al. (2007)’s estimation of the activation energy indicates that transformations are mostly thermally activated, but with a stress-induced bias towards one direction. This will impact the choice of the yield criterion (R2 above). EPM focusing on thermal materials (Bulatov and Argon, 1994a; Ferrero et al., 2014) will set a yield rate based on a stress-biased Arrhenius law for thermal activation, viz.,


where is an attempt frequency, represents the height of the (smallest) potential barrier impeding the rearrangement, tilted by the application of an external stress . Recalling Eq. 2, , one immediately recovers the expression of used by Eyring (1935) to calculate the viscosity of liquids if and are treated as independent parameters, whereas locally imposing a Hookean relation between stress and elastic strain, viz., , leads to the -scaling of the tilt used, e.g., in Sollich et al. (1997)’s Soft Glassy Rheology model.

On the other hand, thermal activation plays virtually no role in foams (Ikeda et al., 2013) and granular materials. Consequently, EPM designed for athermal materials (Hébraud and Lequeux, 1998; Chen et al., 1991) will favor a pure threshold-based yield criterion, viz.,

where is the local stress threshold for yielding; a deterministic yield criterion is recovered in the limit . Incidentally, note that, in this one-dimensional tilt vision, directional considerations in the PEL are handled somewhat light-heartedly; in theory, there is no reason why the loading should push the system towards the saddle point.

As far as rheology is concerned, the athermal approximation is conditioned by the possibility to neglect the activation rate with respect to the driving rate. If we consider a sample sheared at rate , schematically this condition reads


which in essence is similar to the limit of large Péclet number , where is the particle size and is the single-particle diffusivity in the dilute limit (Ikeda et al., 2013). The latter condition is however more (and, possibly, excessively) stringent because in the dense system the diffusivity is much smaller than in the dilute limit.

Now, some subtleties ought to be mentioned. An athermal system may very well be sensitive to temperature variations, through changes in their material properties (e.g., dilation): For example, the fact that Divoux et al. (2008) reported creep motion for a granular heap submitted to cyclic temperature variations does not mean that thermal activation is important, but rather points to dilational effects. Secondly, as already stressed, the relevance of thermal fluctuations may depend on the considered level of detail: It has been argued that they may precipitate the emergence of avalanches by breaking nano-contacts between grains in very slowly sheared systems (Zaitsev et al., 2014), but it is very dubious that this may impact steady-shear granular rheology.

ii.3 Driving

Suppose that the material deforms under the action of some external driving; how important are the specificities of the driving conditions?

Stress or strain driven

As reported in Section I.2, the driving may be stress (fixed ) or strain-controlled (fixed ). Numerical simulations have mostly considered strain-driven situations. In EPM, this affects R1(elastic response) and R3(stress redistribution), i.e., the mechanical response for a fixed set of plastic blocks. In a strain-driven protocol, the elastic response (R1) is generally obtained by converting the macroscopic driving into local stress increments , where is the time step and is the current macroscopic strain rate, and the stress redistribution operated by the elastic propagator keeps the macroscopic strain fixed. Stress-controlled setups have been somewhat less studied in the framework of EPM but examples can be found in (Picard et al., 2004; Lin et al., 2014b, 2015; Homer and Schuh, 2009). In this case, the propagator should keep the macroscopic stress constant. In practice, this often comes down to changing its 0-Fourier mode, which controls the mean value.


EPM often focus on steady shear situations, in which case . But time-dependent driving protocols (or ) are also encountered, in particular step shear , oscillatory shear , which gives access to linear rheology for small . In creeping flows subjected to , the shear rate eventually decays to zero, often as a power law (Leocmach et al., 2014). For further discussions on creep, see Sec. VIII.

Symmetry of the driving

Plastic events are biased towards the direction of the external loading. If the latter acts uniformly on the material, it is convenient to focus on only one stress component, thus reducing the stress and strain tensors to scalars. In particular, for simple shear conditions, with a displacement gradient (in the linear approximation in two dimensions), one may settle with the component. It differs from pure shear, in that the latter is rotationless, whereas the former involves a rotational part , but has the same principal strains (eigenvalues) . These deformations are encountered locally whenever volume changes can be neglected; the cone-and-plate, plate-plate, and Taylor-Couette rheometers (Larson, 1999) used to probe the flow of yield-stress fluids fall in this category. For metallic glasses and other hard materials, uniaxial compression tests (i.e., in the bulk, with ) and tension () are often performed.

Even though in several of these situations the macroscopic is more or less uniform and acts mostly on one component of the (suitably defined) stress tensor, the other components reach finite values because of stress redistribution. Full tensorial approaches may then be justified (Bulatov and Argon, 1994a; Homer and Schuh, 2009; Sandfeld and Zaiser, 2014). Recently, the influence of a tensorial, rather than scalar, description on the flow and avalanche properties in these cases was evaluated; it was found to be insignificant overall (Nicolas et al., 2014b; Budrikis et al., 2017), and the effect of dimensionality to be weak (Liu et al., 2016). The reader is referred to Sec. VII for more details. However, there exist a wide range of experimental setups in which the loading is intrinsically heterogeneous, in particular the bending, torsion, and indentation tests on hard glasses (see (Budrikis et al., 2017) for an implementation of these tests in a finite-element-based EPM) or the microchannel flows of dense emulsions (Nicolas and Barrat, 2013a). The exploitation of EPM in heterogeneous driving conditions appears to be a promising new avenue.

Driving rate

To resolve the flow temporally, the simplest approach is a Eulerian method, which computes the strain increments on all blocks at each time step from Eq. 5. Kinetic Monte-Carlo methods have also been employed and are particularly efficient in stress-controlled slow flows, insofar as the long elastic loading phases without plastic events are bypassed: Activation rates are calculated for all the blocks using a refined version of Eq. 6 and the time lapse before the next plastic event is deduced from the cumulative rate (Homer and Schuh, 2009). If the flow is even slower, rate effects may deliberately be overlooked. Indeed, a number of models consider the limit of vanishing shear rate , where the material is essentially a quiescent solid undergoing intermittent plastic events due of the loading (Baret et al., 2002; Talamali et al., 2012; Lin et al., 2014b). In these extremal models, the algorithm identifies the least stable site at each step and increases the applied stress enough to destabilize it. From this single destabilization an avalanche of plastic events may ensue. Connection to real time is lost. Extremal models are the lattice-based counterpart of quasistatic atomistic simulations, in which a small strain increment is applied at each step and the system then relaxes athermally to the local energy minimum (Maloney and Lemaître, 2004).

ii.4 Rearrangement duration and material time scales

In experiments as well as atomistic simulations, rearrangements are seen to occur very rapidly, so much so that they are considered instantaneous in several rheological models, i.e., in Eq. 5. On the contrary, in various models the finite duration of plastic events plays a major role in the -dependence of the rheology (Picard et al., 2005; Martens et al., 2012; Nicolas et al., 2014a; Liu et al., 2016) or in the intrinsic relaxation of the system (Ferrero et al., 2014).

Suppose that, under slow driving, a rearrangement takes a typical time . For overdamped dynamics, one expects this time scale to be the ratio between an effective microscopic viscosity and the elastic shear modulus (Nicolas and Barrat, 2013a), viz.,

while for underdamped systems is associated with the persistence time of localized vibrations. If the time scale competes with the driving time scale , where is the local yield strain, then plastic events will be disrupted by the driving. The rate dependence of the macroscopic stress may then stem from this disruption (Nicolas et al., 2014a). On the contrary, if the driving is too slow to allow such competition, viz., , then individual rearrangements can be considered instantaneous as far as the rheology is concerned, but the avalanches of rearrangements (i.e., the series of plastic events that would still be triggered by an initial event, were the driving turned off) might still take a finite time, controlled by the signal propagation time within the avalanche. Since the size of avalanches is expected to diverge as in the athermal limit, they may be affected by the driving when . While the delays due to shear wave propagation are generally left aside in EPM, some works have bestowed them a central role in the finite shear-rate rheology (Lemaître and Caroli, 2009; Lin et al., 2014b) and endeavored to represent this propagation in a more realistic way (Nicolas et al., 2015; Karimi and Barrat, 2016). Sec. VI and VII will provide more details on these aspects. Note that the true quasistatic limit is reached when

and the athermal criterion of Eq. 7 is satisfied, i.e., the Péclet number is very large. In that case, the material remains in mechanical equilibrium at all times and its trajectory in the PEL is rate-independent.

ii.5 Spatial disorder in the mechanical properties

Figure 7: Spatial variations of the mechanical and configurational properties of glasses. (a) Maps of the weaker local shear modulus in a 2D Lennard-Jones glass. Black (white) represents values larger (smaller) than the mean value. Distances are in particle size (Lennard-Jones) units. From (Tsamados et al., 2009). (b) and (c) Maps of the local contact-resonance frequency, which is related to the indentation modulus, measured by atomic force acoustic microscopy in (b) a bulk metallic glass and in (c) a crystalline sample of . The latter clearly appears more homogeneous in terms of mechanical properties. The radius of contact is of order . From (Wagner et al., 2011).

Glasses, and more generally amorphous solids, are known to display heterogeneous mechanical properties. Indeed, there have been both experimental and numerical reports on the heterogeneity of the local elastic moduli (see Fig. 7) and the energy barriers (Tsamados et al., 2008; Zargar et al., 2013) on the mesoscale. Yet, the extent to which this disorder impacts the rheology remains unclear. This uncertainty is reflected in EPM. Some models feature no such heterogeneity (Hébraud and Lequeux, 1998; Picard et al., 2005), while it is central in others (Langer, 2008). In the latter case, heterogeneity is generally implemented in the form of a disorder on the yield stresses or energy barriers. Let us mention a couple of examples. In Sollich et al. (1997)’s Soft Glassy Rheology model, energy barriers are exponentially distributed and the exponential dependence of activation rates on the energy barrier (Eq. 6) leads to a transition from Newtonian to non-Newtonian rheology for broader energy distributions (see Sec. IV.4.1 for more details on the model). In their EPM centered on metallic glasses, Li et al. (2013) modify the free energy required for the activation of an event depending on the free volume created during previous rearrangements. Amorphous composite materials, i.e., materials featuring meso/macro-inclusions of another material, can be described as a patchwork of regions of high and low yield stresses (Tyukodi et al., 2016a) or high and low elastic moduli (Chen and Schuh, 2015). In the latter case, macroscopic effective shear and bulk moduli can be derived.

More generally, for single phase materials, the survey of the above results gives the impression that disorder has bearing on the rheology when thermal activation plays an important role. On the other hand, the impact of a yield stress disorder may be less important in athermal systems. In fact Agoritsas et al. (2015) showed that disorder is irrelevant in the mean-field description of athermal plasticity originally proposed by Hébraud and Lequeux (1998), in the low shear rate limit; more precisely, it only affects the coefficients in the rheological law.

Yielding Reference Features \arraybackslashRemarks \arraybackslashProposed for
Activated Bulatov and Argon (1994a) et seq. \arraybackslashPropagator computed on hexagonal lattice \arraybackslashamorphous solids, in particular glasses and glass-forming liquids
Homer and Schuh (2009) et seq. \arraybackslashStress redistribution computed with Finite Elements \arraybackslashmetallic glasses
Ferrero et al. (2014) \arraybackslashPl. events of finite duration \arraybackslashamorphous solids
Sollich et al. (1997) \arraybackslashEffective activation temperature accounts for mechanical noise \arraybackslashsoft materials (foams, emulsions, etc.)
Threshold Chen et al. (1991) \arraybackslashPropagator computed on square spring network \arraybackslashearthquakes
Baret et al. (2002) et seq., Talamali et al. (2011), Budrikis and Zapperi (2013) et seq. \arraybackslashUniform distribution of barriers; coupled to a moving spring (,); or stress controlled with extremal dynamics () or adiabatic driving () \arraybackslashamorphous solids, notably metallic glasses
Dahmen et al. (2011) \arraybackslash‘Narrow’ distribution of thresholds; static and dynamic thresholds differ; mean-field approach \arraybackslashgranular matter and akin
Hébraud and Lequeux (1998) \arraybackslashFinite yield rate above threshold; stress redistributed as white noise \arraybackslashsoft materials (dense suspensions)
Picard et al. (2005), Martens et al. (2012) \arraybackslashFinite yield rate above threshold; pl. events of finite duration \arraybackslashamorphous solids
Nicolas et al. (2014a) et seq \arraybackslashPl. events end after finite strain \arraybackslashsoft athermal amorphous solids
Lin et al. (2014b) \arraybackslashStress- and strain- control protocols \arraybackslashsoft amorphous solids
‘Continuous’ approaches Onuki (2003b) \arraybackslashDynamical evolution on a periodic potential; dipolar propagator due to opposite dislocations \arraybackslash2D crystalline and glassy solids
Jagla (2007) \arraybackslashDynamical evolution on random potential; propagator computed via compatibility condition \arraybackslashamorphous solids
Marmottant and Graner (2013) \arraybackslashOverdamped evolution in a periodic potential; pl. events of finite duration; no stress redistribution \arraybackslashfoams
Legend – Barrier distribution: Single value. Distributed (exponentially, unless otherwise specified).
             Instantaneous plastic events. Quadrupolar elastic propagator.
Table 1: Classification of some of the main EPM in the literature.

ii.6 Spatial resolution of the model

On a related note, one may wonder how important it is to resolve spatially an EPM, or equivalently, in what cases one may settle with a mean-field approach blind to spatial information. Clearly, there are situations in which mean field makes a bad candidate, in particular when the driving or flow is macroscopically heterogeneous, when the focus is on spatial correlations (Nicolas et al., 2014c) or on critical properties (Lin et al., 2014b; Liu et al., 2016). But a mean field analysis could suffice in many other situations. Indeed, Martens et al. (2012) showed that the flow curve obtained with their spatially resolved EPM can be predicted on the basis of mean-field reasoning. Thus, the details of the spatial correlations only had limited effect on the macroscopic stress. Similarly, Ferrero et al. (2014)’s EPM-based simulations confirmed mean-field predictions by Bouchaud and Pitard (2001) regarding thermal relaxation of amorphous solids in some regimes; but not without finding discrepancies in others. In the latter regimes, spatial correlations thus seemed to play a significant role.

The discussion about whether spatial resolution is required to describe global quantities is not settled yet. It has been argued that, owing to the long range of the elastic propagator (which decays radially in d dimensions), mean-field arguments should generically hold in amorphous solids (Dahmen et al., 1998, 2009). However, it has been realized that the non-convex nature of the propagator (alternatively positively and negatively along the azimuthal direction) undermines this argument (Budrikis and Zapperi, 2013) and results in much larger fluctuations than the ones produced by a uniform stress redistribution (Talamali et al., 2011; Nicolas et al., 2014b; Lin et al., 2014a). Mean-field predictions have been tested against the results of lattice-based models simulations of a sheared amorphous solid close to (or in) the limit of vanishing driving, with a focus on the statistics of stress-drops or avalanches, and non-mean-field exponents were found for the power-law distribution of avalanche sizes (Talamali et al., 2011; Budrikis and Zapperi, 2013; Lin et al., 2014b; Liu et al., 2016). This question is addressed in greater depth in Sec. VII.

In this review, we will put the spotlight on spatially resolved models, which are not exactly solvable in general and require a numerical treatment. When relevant, we will discuss how a mean-field treatment can be performed to obtain analytical results.

ii.7 Bird’s eye view of the various models

To conclude this section, some of the main EPM are classified in Table 1.

Iii Elastic couplings and the interaction kernel

A key feature of EPM is to allow plastic events to interact via an elastic deformation field, which can generate avalanches. In this respect, the choice of the elastic interaction kernel may significantly impact the results of the simulations (Martens et al., 2012; Budrikis and Zapperi, 2013). This fairly technical section presents the various idealizations of the interaction kernel that have been used in the literature on amorphous solids, by increasing order of sophistication. We endeavor to relate this level of sophistication with the nature of the developments that were sought.

iii.1 Sandpile models and first-neighbor stress redistribution

Figure 8: (a) Sketch of the discrete 1D Burridge and Knopoff model. From (Carlson et al., 1994). (b) and (c) Chen et al. (1991)’s spring network model. (b) Sketch of the effect of a bond rupture in the model. (c) Distribution of avalanche sizes, in terms of number of broken bonds. Adapted from (Chen et al., 1991).

The watershed between the models for earthquakes and general avalanches and the more recent EPM is often fuzzy. In fact, the latter literally burgeoned on avalanche and earthquake-ridden scientific grounds.

As a paradigmatic earthquake model, consider the celebrated model by Burridge and Knopoff (1967), whose main features are concisely reviewed in (Carlson et al., 1994). It focuses on the fault separating two slowly moving tectonic plates. This region is structurally weak because of the gouge (crushed rock powder) it is made of; thus, failure tends to localize along this fault. In the model, the contact points across the fault are represented by massive blocks and the compressive and shear forces acting along it are modeled as springs, as sketched in Fig. 8a. Due to these forces, the initially pinned (stuck) blocks may slide during avalanches. More precisely, in the continuous, nondimensional 1D form, the displacement at time of the material at position reads


Here, the left-hand side (lhs) is related to inertia, the second derivative on the right-hand side (rhs) is of compressive origin, and the loading term due to the motion of the plate as well as the displacement contribute to a shear term. Finally, is a velocity-dependent frictional term. Had Coulomb’s law of friction been used, it would have been constant for , but the original model assumed velocity weakening, i.e., a decrease of with . At , the function is degenerate, which allows static friction to exactly cancel the sum of forces on the rhs of Eq. (8), so the blocks remain pinned at a fixed position until the destabilizing forces exceed a certain threshold. Phenomenologically, simulations of the model show frequent small events (with a power-law distribution of cumulative slip) and rare events of large magnitude, in which the destabilization of a number of sites close to instability results in a perturbation of large amplitude (Otsuka, 1972; Carlson et al., 1994).

Important in the above model is the effect of the pinning force at . It entails that the destabilizing action caused by the depinning of a site (via the diffusive term in Eq. 8) is fully screened by its neighbors, unless they yield too. Such first-neighbor redistribution of strain is readily simulated using cellular automata, which can be interpreted as sandpile models: Whenever a column of sand, labelled , gets too high with respect to its neighbors (say, for convenience, whenever ), some grains at its top are transferred to the neighboring columns, with the following discharge rules in two dimensions:


where is the height difference between columns. The sandpile is loaded by randomly strewing grains over it in a quasistatic manner. The study of these systems soared in the late 1980s and early 1990s, whence the overbearing concept of self-organized criticality emerged (Bak et al., 1987). According to the latter, the avalanches naturally drive the sandpiles toward marginally stable states, with no characteristic lengthscale for the regions on the verge of instability, hence the observation of scale-free frequency distributions of avalanche sizes. As an aside, let us mention that this approach has not been used only for earthquakes (Carlson and Langer, 1989; Sornette and Sornette, 1989; Bak and Tang, 1989; Ito and Matsuzaki, 1990) and avalanches in sandpiles, it has also been transposed to the study of integrate-and-fire cells (Corral et al., 1995) and forest fires (Chen et al., 1990), inter alia.

In seismology, these models have been fairly successful in reproducing the Gutenberg and Richter (1944) statistics of earthquake. This empirical law states that the frequency of earthquakes of (energy) magnitude


where is the energy release, in a given region obeys the power law relation, , where , or equivalently

For the sake of accuracy, we ought to say that there exist several earthquake magnitude scales besides that of Eq. 10. They roughly coincide at not too large values; in fact, is not the initial Richter scale. More importantly, the value of the exponent depends on the considered earthquake catalog, notably on the considered region. For sandpile-like models, various exponents have been reported: in 2D and in 3D, with no effect of a disorder on the yield stresses (Bak and Tang, 1989), whereas the exponent for the mean-field democratic fiber bundle close to global failure is (see Section IX.3). More extensive numerical simulations have led to the values (Lübeck and Usadel, 1997), or (Chessa et al., 1999), for the 2D Bak et al. (1987) sandpile model.

Olami et al. (1992) modified the model to make the redistribution rule of Eq. 9 non-conservative. In this sandpile picture, this would correspond to a net loss of grains, which seems unphysical; but in Burridge and Knopoff (1967)’s block-and-spring model the non-conservative parameter simply refers to the fraction of strain which is absorbed by the driving plate during an event, instead of being transferred to the neighbors. Interestingly, as non-conservativeness increases, criticality is maintained, insofar as the avalanche distribution remains scale-free, even though the critical exponent gradually gets larger. Only when less than 20% of the strain is transferred to the neighbors does a transition to an exponential distribution occur. The dynamics then become more and more local with increasing dissipation, until the blocks completely stop interacting, when the redistribution is purely dissipative.

However, unlike the redistribution of grains in the sandpile model, elastic interactions are actually long-ranged, as we wrote in Section I.4. In particular, in the deformation of amorphous solids, no pinning of the region surrounding an event can be invoked to justify the restriction of the interaction to the first neighbors.

iii.2 Networks of springs

Accordingly, a more realistic account of the long-ranged elastic propagation is desirable. Unfortunately, the complexity of the bona fide Eshelby propagator obtained from Continuum Mechanics hampers its numerical implementation and use, so most studies have relied on simplified propagators, which share some similarities with Eshelby’s.

First, in the spirit of the classical description of a solid as an assembly of particles confined to their positions by the interactions with their neighbors, the material was modeled as a system of blocks connected by “springs” of stiffness and potential energy

where is the displacement of block i. Note that this expression for the potential energy entails noncentral forces, so that the “springs” can bear shear forces; some details about the difference with respect to networks of conventional springs are presented in Section IX.3. The pioneering steps towards EPM followed from the application of such spring network models to the study of rupture. For this purpose, each bond is endowed with a random threshold, above which it yields and redistributes the force that it used to bear. In their study of a 2D triangular lattice with central forces, Hansen et al. (1989) measured the evolution of the applied force with the displacement ; this evolution starts with a phase of linear increase, followed by a peak and a smooth decline until global failure. The curves for different linear lattice sizes roughly collapsed onto a master curve if and were rescaled by . In addition, just before failure, the distribution of forces in the system was “multifractal”, with no characteristic value.

Chen et al. (1991) considered a square lattice of blocks and “springs”, sketched in Fig. 8b. The rupture of a spring triggers the release of a dipole of opposite point forces (generating vorticity) on neighboring blocks. This differs somewhat from the force quadrupole corresponding to an (irrotational) local shear (see Section III.3), but also leads to an anisotropic shape. Contrary to Hansen et al. (1989), they allowed broken springs to instantly regenerate to an unloaded state, after the redistribution of their load. Physically, this discrepancy parallels a change of focus, from brittle materials to earthquakes, for which the external loading due to tectonic movements is assumed to be by far slower than the healing of bonds. For a quasistatic increase of the load, the model displays intermittent dynamics and scale-free avalanches, and a power-law exponent was reported in 2D, in semi-quantitative agreement with the Gutenberg-Richter earthquake statistics.

iii.3 Pointwise idealization of the Eshelby propagator


An alternative to block-and-spring models, rooted in Continuum Mechanics and popularized by Picard et al. (2004), consists in simplifying Eshelby (1957)’s calculations of the elastic propagator by considering the pointwise circular limit of a 2D shear transformation. The latter then contributes to the equation of mechanical equilibrium as a source term , viz.,


The surrounding medium is supposed to be linear elastic and uniform, so that its elastic stress reads


where is the pressure and is the linear strain tensor. For simplicity, incompressibity is assumed, . The solution of Eq. 11 is then well known in hydrodynamics and involves the Oseen-Burgers tensor , with the identity matrix, viz.,


Finally, to get the source term , one assumes that a shear transformation located at the origin locally shifts the unstrained configuration by an amount : the configuration that cancels the shear stress in is no longer , but , in the limit of a rearrangement of vanishing volume . In this sense, can be regarded as an eigenstrain that generates a source term in Eq. 11. In the unbounded 2D plane, setting coordinates such that , the response to in terms of xy-component of the stress reads


where are polar coordinates. This field is shown in Fig. 6c. Reassuringly, in the far field this coincides with the response to a cylindrical Eshelby inclusion.

As a short aside, let us mention a variant to these calculations, which puts in the limelight the connection with deformation processes in a crystal. This variant is reminiscent of Eshelby’s a cut-and-glue method, in which an ellipsoid is cut out of the material, deformed, and then reinserted. Following earlier endeavors by Ben-Zion and Rice (1993), Tüzes et al. (2017) carved out a square around the rearrangement, instead of an ellipsoid, displaced its edges to mimick shear, and then glued it back. This is tantamount to inserting four edge dislocations in the region and also yields an Eshelby-like quadrupolar field.

Rather than focusing on unbounded media, it is convenient to work in a bounded system with periodic boundary conditions and with a general plastic strain field . Switching to Fourier space (), the counterpart of Eq. 14 is then


Note that the frame is sometimes defined such that ; in this case, . In practice, the system will generally be discretized into a (square) lattice, which allows one to use a Fast Fourier Transform routine and restrict the considered wavenumbers to .

Besides, because of dissipative forces, quantified by an effective viscosity , the plastic strain rate in the shear transformation cannot be infinite and a rearrangement will last for a finite time (see Section II.4). Therefore, in each numerical time step , the plastic strain implemented in Eq. 15 will be the strain increment during that step. This amounts to saying that, locally, the rearrangement occurs gradually, even though the redistribution of stress to the rest of the medium is instantaneous (because mechanical equilibrium was assumed, so that there is no time dependence in the elastic propagator in Eq. 15).

Issues with this approximation and possible remediations

The idealized elastic propagator in Eq. 15 brings on some technical issues. First, its slow () radial decay in space raises convergence problems in periodic space. Indeed, the fields created by the periodic images of each plastic event have to be summed, but the sum converges only conditionally in real space, i.e., depends on the order of summation. This is reflected by the singularity of near . In polar crystals, such a difficulty also arises, when computing the Madelung energy, but may be solved with the Ewald (1921) method. Here, we make use of the conserved quantities to state that in a stress-controlled system and in a strain-controlled system. Another possibility is to sum the images in an arbitrary order that is compatible with convergence. These distinct implementations match in the far field, but differ in the near field, which leads to different organizations for the flow (Budrikis and Zapperi, 2013).

Second, on a periodic lattice, one should in principle compute the periodic sum

if, at the lattice nodes, one wishes the backward discrete Fourier transform of to coincide with the solution for an unbounded medium. However, the high-frequency components in , due to the spurious singularity of at (Eq. 15), make the periodic sum diverge. In practice, wavenumbers outside the first Brillouin zone are plainly discarded, which comes down to solving Eqs. 11-12 on the periodic lattice, rather than in the continuum. Nevertheless, spurious fluctuations in the response field are sometimes observed; the problem is mitigated by using a finer grid and smoothing the obtained field (Nicolas et al., 2014b).

Variations: Soft modes and lattice symmetries; tensoriality; convection

All in all, many technical details of the implementation of the elastic propagator appear to affect the spatial organization of the flow (Talamali et al., 2011), but leave the qualitative picture and (apparently) the scaling laws unaltered. However, an aspect that seems to be crucial is the need to preserve the eigenmodes of the propagator associated with eigenvalue (energy) zero. These so called soft modes are the fields such that

and, since they cost no elastic energy, their deployment is favored by the dynamics (Tyukodi et al., 2016b). Their importance is further explained in Section IV.3. It turns out that the eigenmodes of in Eq. 15 are simply the Fourier modes (plane waves); among these, the soft modes are those with wavevectors that make an angle with respect to the principal direction of the plastic strain tensor .

In particular, if the material is loaded under simple shear with velocity direction and velocity gradient along , the emergence of a uniform shear band along should produce no elastic stress in the medium, at least if such a band emerges uniformly. However, misaligned lattice axes, not directed along or , would not be compatible with such a shear band (which would then have sawtooth-like edges) and artificially suppress the soft modes (Tyukodi, 2016). More generally, the use of a regular lattice in EPM may be questioned, insofar as the localization of plastic events is sensitive to variations of stress redistribution in the near field (Budrikis and Zapperi, 2013), but the scalings of, e.g., avalanche sizes seem to be mostly insensitive to these details. Indeed, these details do not affect the long-range elastic propagator.

On another note, the foregoing calculations focused on the xy-shear stress component, because of the macroscopic stress symmetry, thus promoting a scalar description. It is straightforward to generalize the reasoning to a fully tensorial form; but it turns out that, for setups with uniform loading, the tensorial extension has virtually no effect (Nicolas et al., 2014b). At most, tenuous differences are reported in the statistics of avalanches of plastic events when moving from a scalar model to a tensorial one: the values of the critical exponents at the yielding transition reported for scalar models (Sandfeld et al., 2015) are close to those obtained in the corresponding tensorial models (Budrikis et al., 2017). Similarly, moving from 2D to 3D does not introduce qualitative changes and scaling relations are preserved (Budrikis et al., 2017; Liu et al., 2016). Also, the periodic boundary conditions can be substituted by no slip boundary conditions at a wall, via the image method (Picard et al., 2004), for instance to model flow in a microchannel (Nicolas and Barrat, 2013a). Translational invariance of the propagator is then broken and plastic events relax stress faster, for a given eigenstrain, if they occur close to the walls. Changes in the boundary conditions used to compute the redistributed stresses affect the spatial organization of the flow, but not the critical properties at the yielding transition (Sandfeld et al., 2015). Finally, despite the convenience of using a fixed lattice grid with static elasto-plastic blocks, physically these blocks should be advected by the flow. In a bounded medium, a coarse version of advection can be implemented by incrementally shifting the blocks along the streamlines without altering the global shape of the lattice (Nicolas and Barrat, 2013a). On the other hand, with periodic boundary conditions, the deformation of the frame results in the shift of the periodic images with respect to the simulation cell; advection thus requires to compute the elastic propagator afresh, in the deformed frame (Nicolas et al., 2014b).

iii.4 Finite-Element-based approaches

Figure 9: Average displacement field induced by a shear transformation in an underdamped elastic medium, computed with a basic Finite Element routine. The plotted snapshots correspond to different delays after the transformation was (artificially) triggered at the origin: (a) , (b) , (c) . Red hues indicate larger displacements. Adapted from (Nicolas et al., 2015)

Albeit computationally more costly, Finite-Element (FE)-based computations of stress redistributions overcome some limitations of the foregoing approaches and offer more flexibility. The FE method solves the continuum mechanics equation within each element of a meshgrid by interpolating the local strain and stress from the values of the displacements and point forces at the nodes of the element.

If mechanical equilibrium is maintained at all times, the stress redistributed by a shear transformation can be computed by equilibrating the elastic stress generated by the eigenstrain borne by a given element (where is the stiffness tensor). Using a triangular mesh refined around this eigenstrain-bearing element, Sandfeld et al. (2015) demonstrated that the computed stress field coincides with the elastic propagator of Eq. 14 in the pointwise rearrangement limit. But these researchers also found that a coarser mesh made of uniform square elements gives results that are almost as good, except in a near-field region of a handful of sites’ radius. The flexibility of the method was then exploited to study the quasistatic deformation of the system beyond the periodic boundary conditions, e.g., in a bounded medium and with free surfaces, and with inhomogeneous loading conditions (indentation, bending, etc.). Universal, but non-mean-field, statistics of avalanches of plastic events were reported in these diverse conditions (Budrikis et al., 2017) (also see Section VII).

In an earlier endeavor (Homer and Schuh, 2009; Homer et al., 2010), each shear transformation zone consisted of several elements of a triangular mesh which all bore an eigenstrain. As the size of this zone increases, the redistributed stress field accurately converged to the theoretical Eshelby field. Zones made of 13 elements were deemed quite satisfactory in this respect. Homer and Schuh (2010) later extended the approach to 3D. Dynamics were brought into play via the implementation of an event-driven (Kinetic Monte Carlo) scheme determining the thermal activation of shear transformations, in the wake of the pioneering works of Bulatov and Argon (1994a). The cooling of the system, its thermal relaxation and its rheology under applied stress were then studied. Macroscopically homogeneous flows were observed at low stresses and/or high temperatures, whereas the strain localized at low temperature for initially unequilibrated (zero residual stress) systems, which was not necessarily supported by experimental data. More systematic strain localization at low temperature was found by Li et al. (2013), who incorporated the processes of free volume creation during plastic rearrangements and subsequent free volume annihilation (see Section V.3.2).

The capabilities of FE methods were further exploited by Nicolas et al. (2015) to go beyond the assumption of elastic homogeneity of the material. For this purpose, the mesoscale elastic constants in regions of 5 particles in diameter were measured in an atomistic glass model; the local shear moduli were found to be broadly distributed, with relative fluctuations of around 30% and marked anisotropy (i.e., one direction of shear being much weaker than the other one). Introducing this heterogeneity of shear moduli in the EPM sufficed to capture the sample-to-sample fluctuations of the elastic response to an artificially triggered shear transformation observed in Molecular Dynamics simulations (Puosi et al., 2014); accounting for anisotropy was less critical. Besides, inertial and viscous terms were not omitted in the FE description, so that mechanical equilibration was not instantaneous and shear waves were seen to propagate in the transient, as in the atomistic simulations. The natural inclusion of inertia in FE was also exploited by Karimi et al. (2017) to analyze the effect of inertia on the universal avalanche statistics and compare with atomistic simulations directly. (Note that the effect of a delay in signal propagation had already been contemplated in an effective way by Lin et al. (2014a), while, for the same purpose, Papanikolaou (2016) introduced a pinning delay in his EPM based on the depinning framework.) It was then possible to investigate the influence of the damping strength on the rheology of the elastoplastic system, which was indeed done by Karimi and Barrat (2016). Using a Maxwellian fluid description for blocks in the plastic regime and an unstructured mesh, these researchers found trends qualitatively very similar to what is observed in Molecular Dynamics when the friction coefficient is varied.

iii.5 Continuous approaches based on periodic potentials

Notwithstanding their variable sophistication, all above methods rest on a clearcut distinction between plastic rearrangements and elastic deformations. This binary distinction is relaxed in continuous approaches based on a free energy functional that depends on the mechanical strain variables, generically denoted by here. Since the energy increase due to elastic loading is cyclically interrupted by plastic rearrangements (associated with discrete jumps between valleys of the energy landscape), involves periodic functions of the plastic deformation or the total deformation . For instance, Marmottant and Graner (2013) made use of the following effective potential,

where is an elastic modulus, is a yield strain and is the period of the pinning potential. If this prescription is coupled with a dynamical equation of the form

with the characteristic relaxation timescale (leading to the Prandtl–Tomlinson model for stick–slip), a serrated stress vs. strain curve is obtained under constant driving. The finite time needed by the plastic deformation to jump between energy valleys implies that, at high driving rates, will not be able to instantaneously jump between, say, and . Therefore, the elastic strain will continue increasing in the valley around for some time, although the criterion for the onset of plasticity has already been met, which is similar to having a finite latency time prior to relaxation once the threshold is exceeded in Picard et al. (2005)’s model. Similar equations of motion in a random potential have been proposed for solid friction; the occurrence of stick-slip dynamics owes to the “pinning” of the system in one potential valley, up to some threshold, while there exists another stable position (Tyukodi et al., 2016b).

To go beyond the mean-field level, this type of continuous approach can be resolved spatially. In an inspiratonial study, Onuki (2003a) introduced an elastic free energy of the form


where is the bulk modulus and the volumetric strain as well as the shear strains and are explicit functions of the displacement field . Here, is an arbitrarily chosen function that is invariant under rotations of the reference frame (because a 2D triangular lattice is assumed) and periodic in its arguments. Introducing in the equation of motion


where is the density, is the viscosity and is a random stress tensor due to thermal fluctuations, suffices to obtain qualitatively realistic stress vs. strain curves. The framework was then extended to study the effect of an interplay between the volumetric strain and the density , and to capture the elastic effects of edge dislocations, if the material is crystalline (Onuki, 2003b).

If the strain components , and are handled as independent primary variables, instead of being functions of as in Eq. 17, then their compatibility as components of a strain tensor should be ensured by the Saint-Venant condition

This constraint is implemented by means of a Lagrange multiplier in the total free energy , viz., . A minimization of yields equations of motion such as , for a generic strain component , that complete the definition of the EPM (Jagla, 2007). If the function entering the free energy in Eq. 16 is assumed to depend quadratically on (as in linear elasticity), while its nonlinearity with respect to preserves the possibility of plastic relaxation, then it can be shown analytically that, owing to the Saint-Venant condition, a plastic strain along gradually gives rise to an elastic strain with a quadrupolar structure (Kartha et al., 1995). Once mechanical equilibrium is reached, the final strain field produced by this local eigenstrain coincides in the incompressible limit with the elastic propagator of Eq. 15 (Jagla, 2017). The approach was extended recently by Jagla (2017) to study the influence of different choices for the effective plastic disorder potentials in the flow curve and critical exponents.

Iv Mechanical noise and its approximations

The previous section has shed light on the modeling of the elastic propagator, i.e., the effect of a single rearrangement on the surrounding elastic medium. In practice, however, several rearrangements may occur simultaneously, and the rate of stress increment received by a given block (say, site ) at time is then a sum of contributions from many sites, i.e., using Eq. 5, , where denotes the plastic activity of site . Due to its fluctuating nature, this quantity is often referred to as mechanical noise. By rewriting Eq. 5 as


one can readily see that, in combination with the external loading and the dynamical rules governing , the mechanical noise signal fully determines the local stress evolution. All “one-point” properties (such as the flow curve, the density of plastic sites, the distribution of local stresses, etc.) can be obtained by averaging the local properties at over time. This shows the central role of in determining these properties. Unfortunately, this signal is complex, as it stems from interacting plastic events throughout the system; nevertheless, mean-field approaches suggest to substitute it with a simpler “mean” field.

iv.1 Uniform redistribution of stress

The mechanical noise can be split into :

  • a constant background , which contributes to a drift term in Eq. 18, and

  • zero-average fluctuations .

Owing to the infinite range and slow decay of the elastic propagator ( in d-dimensional space, see Sec. I.4), site is significantly coupled to a large number of other sites. This large connectivity has led some researchers to overlook the fluctuations in favor of the average drift term. Along these lines, in the framework of Picard’s EPM, which features a constant rate of yield above a uniform threshold and a constant rate of elastic recovery, viz.,


Martens et al. (2012) averaged Eq. 18 over time and found a mean-field analytical expression for the flow curve, which reproduces the simulation results to a large extent. It also correctly predicts the destabilization of the homogeneous flow leading to shear-banding for a range of model parameters, in particular at large .

In fact, the neglect of fluctuations would be rigorously justified if the system were infinite and the non-convexity of the propagator were left aside. The latter criterion is for instance fulfilled in a simple quasistatic model in which sites yield past a threshold and redistribute the released stress () uniformly to the other sites (Dahmen et al., 1998), viz.,

The simplicity of the model allows analytical progress. A first approach consists in treating the distances to the threshold as independent variables in the system and sorting them in ascending order (). An avalanche will persist as long as the stress increment due to the yielding of the most unstable site suffices to make the second most unstable fail, viz., . Using an argument along these lines in a model featuring disorder in the yield thresholds () and post-failure weakening (i.e., when site i yields, the threshold is restored to a lower value ), Dahmen et al. (1998) were able to rationalize the existence of a regime of power-law distributed avalanches and a regime of runaway, system-spanning avalanches.

Alternatively, owing to the similarity of the simplified problem with force-driven depinning, one can make use of the machinery developed in the latter field. Transversal scaling arguments and renormalization group expansions (Fisher et al., 1997) then allow one to derive scalings for different properties of the system in the quasistatic limit, such as the size of avalanches. Note that this method was initially applied to the depinning problem and to earthquakes. Only later on was it claimed to be much more general and to have bearing on very diverse systems exhibiting intermittent dynamics or “crackling noise” (Sethna et al., 2001), in particular the yielding transition of amorphous solids (see Sec. VII). Recently, these mean-field scaling predictions about avalanche sizes, shapes, and dynamics have been used to fit experimental data, in metallic glasses subjected to extremely slow uniaxial compression (Antonaglia et al., 2014; Dahmen et al., 2009) as well as in compacted granular matter (Denisov et al., 2016).

iv.2 Random stress redistribution

Deviations from uniform mean field

Notwithstanding this success, the underpinning of the foregoing mean-field approach has been called into question. Theoretically, the argument based on the long range of the interactions is undermined by the fact that these interactions are sometimes positive and sometimes negative (Budrikis and Zapperi, 2013). A crude estimate of the ratio of fluctations over mean value of the stress increments points to a divergence of this ratio at low shear rates and therefore to the failure of the mean-field theory, according to Ginzburg and Landau’s criterion (Nicolas et al., 2014b). Numerically, some lattice-based simulations do indeed reveal departures from mean-field predictions for the critical exponents (Budrikis and Zapperi, 2013; Lin et al., 2014a; Liu et al., 2016). For instance, in these simulations, near , the distribution of avalanche sizes follows a power law with an exponent that deviates from the value predicted by mean field (see Sec. VII for details).

The Hébraud-Lequeux model

To improve on the hypothesis of a constant mean field , fluctuations of the mechanical noise need to be accounted for. In the crudest approximation, they can be substituted by random white noise , with . This turns Eq. 18 into a biased Brownian walk for the local stresses, in the elastic regime . Hébraud and Lequeux (1998)’s model was developed along these lines. The ensuing stochastic equation (Eq. 18 with and ) can be recast into a probabilistic Fokker-Planck-like equation operating on the distribution of local stresses , viz.,


where the diffusive term on the rhs arises from the fluctuations acting on , with a coefficient assumed to be proportional to the number of plastic sites , viz., . The first term on the rhs of Eq. 20 is a drift term, which amalgamates with ; the last two terms correspond to the failure of overloaded sites (above ) on a timescale and their rebirth at due to stress relaxation. The resulting mean-field equations can be solved in the limit of vanishing shear rates (Olivier, 2011; Agoritsas et al., 2015). For a coupling constant , diffusion vanishes at low shear rates, with , a yield stress is obtained and the average stress obeys , with , in the limit of slow shear rates. For , the system behaves like a Newtonian liquid.

Fraction of sites close to yielding

The diffusive term introduced in Eq. 20 impacts the distribution of sites close to yield, i.e., at distances from the yield threshold , where . On these short distances, or, equivalently, in the limit of short time scales , the back-and-forth diffusive motion over typical distances prevails over the drift in the random walk. Therefore, for , determining the distribution is tantamount to finding the concentration of Brownian particles near an absorbing boundary at (yield): the well-known solution is a linear vanishment of the concentration near , viz., for (Lin et al., 2014a; Lin and Wyart, 2016). This result ought to be compared with for drift-dominated problems, such as depinning. Lin et al. (2014a) further claim that this discrepancy is at the origin of the differences in scaling behavior between the depinning transition ( with ) and the flow of disordered solids ( with generally).

Still, the singularity of the propagator casts doubt on the Gaussian nature of the random stress increments and would rather suggest a broad density function

in the limit of sparse plastic events, with an upper cut-off proportional to the volume of a rearranging region. For such a heavy-tailed distribution, the biased random Brownian walk of is replaced by a Lévy flight of index for . On the basis of a simple extremal model, Lemaître and Caroli (2007) demonstrated that this change altered the avalanche statistics as well as the distribution of distances to yield . To be explicit, their model was based on plastic yield above a uniform yield strain , which resets the local stress to zero and increments the stresses at other sites by random values drawn from ). Lin and Wyart (2016) explored a related, but more general (Hébraud-Lequeux-like) model analytically and confirmed the impact of the non-Gaussianity of the noise. These authors derived an asymptotic expression for , which scales as for . The exponent is a non-universal exponent that depends on the loading and the amplitude of the noise and supplements the other two exponents characterizing the depinning transition.

iv.3 Validity of the above “mean-field” approximations

The foregoing paragraphs have presented distinct levels of “mean-field” approximations. Now we enquire into their range of validity and record the results in Table 2.

Uniform mean field

Clearly, the neglect of fluctuations in the constant mean-field approach is sensible only in the drift-dominated regime, i.e., when on the considered time scales , with the notations of Sec. IV.1. With interactions that change signs, this excludes vanishing shear rates or too small time windows . But at high shear rates, this approach appears to correctly predict the avalanche scaling exponents in the EPM studied by Liu et al. (2016).

White-noise fluctuations

Complemented with Gaussian fluctuations, the approximation is valid beyond the drift-dominated regime. In fact, any mechanical noise signal with (i) finite average and finite variance, in particular , and (ii) no significant time correlations, can be replaced by Gaussian white noise in Eq. 18 (Lin and Wyart, 2016). Accordingly, the universality class of the Hébraud-Lequeux model encompasses all models based on similar rules for plasticity and where the mechanical noise fulfils the above criteria (i) and (ii). In particular, for coupling parameters such that the diffusivity goes to zero at , their flow curves will follow a Herschel-Bulkley behavior with in the low shear rate limit. This holds true in the presence of disorder on the local yield thresholds (Agoritsas et al., 2015) and for plastic events that do not relax the local stress strictly to zero, but to a low random value (Agoritsas and Martens, 2017). On the other hand, should the shear modulus of elastic blocks or the relaxation time display a shear-rate dependence, the exponent will deviate from (Agoritsas and Martens, 2017).

Besides, even in an elastic system with sparse plastic events, where the formula for the elastic propagator suggests , the finite-variance constraint (i) on could be fulfilled. Indeed, the elastic propagator is but a far-field approximation and the large values predicted in the near field have not been observed numerically or experimentally. Consequently, there will always be an upper cut-off in the predicted -distribution, which may justify the Gaussian noise approximation on long time scales. In this regard, it is interesting to note that the extremal model studied by Lemaître and Caroli (2007) seemed to capture the avalanche scaling exponent measured in atomistic simulations better with Gaussian fluctuations than with the heavy-tailed fluctuations suggested by the elastic propagator.

Heavy-tailed fluctuations

Nevertheless, if it so happens that the fluctuation distribution is heavy-tailed, the breadth of the mechanical noise distribution should not be overlooked, as it modifies the expected distribution of local distances to yield . Yet, in dimensions and , taking into account the heavy tail does not suffice to capture the exponent (defined by ) measured in lattice-based simulations relying on the genuine elastic propagator (Lin and Wyart, 2016). This points to the failure of criterion (ii) above, i.e., the importance of temporal correlations in the mechanical noise signal ; these correlations are notably due to avalanches and are lost in the mean-field reasoning. Only for a higher dimension does the mean-field prediction for come closer to the value measured in the lattice-based model, which suggests an upper critical dimension . This claim may be contrasted with Chen et al. (1991)’s early speculation of an upper critical dimension of 3 for the applicability of constant mean field in their model. Leaving aside mean-field concerns for a minute, we find quite noteworthy that the exponents measured in the lattice-based EPM are quite compatible with their (indirectly) measured value in atomistic simulations in the quasistatic regime, in 2D and 3D (Lin et al., 2014a).

Structure of the elastic propagator and soft modes

Coming back to the comparison between EPM and their approximations, we note that the mean-field predictions are recovered (even for ) if the elastic propagator is shuffled, that is, if the coupling between sites i and j is given by , where is a random permutation of indices which changes at each time step, instead of being given by . This shows that the temporal correlations in the mechanical noise signal arise because of the spatial structure of . Of particular importance in this structure, claim Talamali et al. (2012) and Tyukodi et al. (2016b), are the soft deformation modes of the propagator (the uniform shear bands described in Sec. III.3), that create no elastic stress in the material. To clarify this importance, the authors focused on the evolution of the cumulative plastic strain in extremal dynamics and recast the EPM equation of motion (Eq. 5) into a depinning-like equation (also see Sec. IX.2), viz.,

where is a viscosity, is the local stress threshold, and if and 0 otherwise. The deformation of a disordered solid in d dimensions is then regarded as the advance of a -dimensional elastic hypersurface in a d+1-dimensional space, where the additional dimension is . Under the influence of the driving, the elastic hypersurface moves forward along , and, in so doing, gets deformed owing to the disorder in the yield thresholds seen by different points on the hypersurface.

Despite the similar framework, Tyukodi et al. (2016b) showed that EPM will considerably deviate from elastic depinning problems because of the existence of soft modes in the EPM kernel , while nontrivial soft modes are prohibited by the definite positiveness of the depinning propagator. As time goes on, the width of the elastic hypersurface (where the overbar denotes a spatial average and the brackets indicate an ensemble average over the disorder) saturates in the depinning problem. This saturation is due to the higher elastic stresses released by regions of higher , which destabilize regions of lower and therefore act as restoring forces to homogenize over the hypersurface. On the contrary, in EPM, (the variance of grows endlessly by populating the soft modes of plastic deformation, which generate no elastic restoring force, and diverges in a diffusive fashion at long times.

iv.4 A mechanical noise activation temperature?

The Soft Glassy Rheology model (SGR)

The Soft Glassy Rheology model of Sollich et al. (1997) proposed an alternative way to handle mechanical noise fluctuations . In the SGR spirit, these random stress “kicks” operate as an effective temperature that can activate plastic events, in the same way as thermal fluctuations do. Accordingly, the fluctuation-induced diffusive term in Eq. 20 is replaced by an Arrhenius law to describe activated effects. More precisely, in SGR, the material is divided into mesoscopic regions that evolve in a landscape of traps whose depths are randomly drawn from an exponential distribution (Bouchaud, 1992)

Here, is a material parameter that will be set to unity. To hop from trap to