Deformation and flow of amorphous solids: An updated 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.

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 loaded above a threshold. 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. 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 seem to 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 N 1 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. 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 we describe the current understanding of strain localization based on the study of EPM. Section VI focuses on the statistical marks of criticality encountered when the system is driven extremely slowly, especially in terms of the temporal and spatial organization of stress fluctuations in 'avalanches', while Sec. VII describes the bulk rheology of amorphous solids in response to a shear deformation. Section 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.
These sections are largely self-contained and can thus be read separately. Sections I and II are both particularly well suited as entry points for newcomers in the field, while Sections III and IV might be more technical and of greater relevance for the experts interested in the implementation of EPM. Finally, Sections V to VIII focus on applications of the models to specific physical phenomena and are largely independent from each other.

I. GENERAL PHENOMENOLOGY
A. 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, and will flow only if a sufficient load is applied to them. Accordingly, in the rheology of complex fluids In the event of material failure, which is generally preceded by strain localization, the stress dramatically drops down. (b) Steady-state flow curve, i.e., dependence of the steady-state shear stress Σss on the shear rateγ, represented with semi-logarithmic axes. If the flow is split into macroscopic shear bands, a stress plateau is generally observed. (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, are solid because surface tension strives to restore the equilibrium shape of their constituent bubbles or droplets upon deformation. Their elastic moduli, i.e., the coefficients of proportionality between the elastic strain and stress, are then approximately given by the surface tension divided by the bubble or droplet radius (between a few microns and several millimeters); a few hundred Pascal would be a good order of magnitude. Colloidal glasses, on the other hand, are dense suspensions of solid particles of less than a micron in diameter, which makes them light enough for Brownian agitation to impede sedimentation. They rely on entropic forces to maintain their reference structure and typically have shear moduli of the order of 10 − 100 Pa. 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 linear 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 vague 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 100 GPa (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 10 5 − 10 6 K · s −1 are required for metallic glasses (Greer, 1995;Greer and Ma, 2007), whereas values below 1K · s −1 may be used for oxide glasses. After a certain amount of deformation, brittle materials will break without incurring significant plastic (irretrievable) deformation, whereas ductile materials will deform plastically before breaking. We will discuss connections between these forms of deformation in Sec. V.

B. What controls the dynamics of amorphous solids?
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 ( 1µm) 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 T or the angular velocity Ω of the rotating part. In the former case, the applied macroscopic shear stress is kept fixed, at a value Σ = T 2πhR 2 on a rotating cylinder of radius R and height h (Fardin et al., 2014), while one monitors the resulting shear strain γ or shear rateγ if the material flows steadily. Conversely, strain-controlled experiments impose γ(t) oṙ γ and monitor the stress response Σ(t). How do amorphous solids respond to such external forces? For small applied stresses Σ, the deformation is elastic, i.e., mostly reversible. Submitted to larger stresses, the material shows signs of plastic (irreversible) deformation; but the latter ceases rapidly, unless Σ overcomes a critical threshold Σ y known as yield stress (see Fig.3(a)). For Σ > Σ y , 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 stationary plastic flow. The ensuing flow curve Σ = f (γ) in the steady state is often fitted by a Herschel-Bulkley law with n > 0 (see Fig. 3(b)).
The transition between the solid-like elastic response and the irreversible plastic deformation is known as the yielding transition. Statistical physicists are inclined to regard it as a dynamical phase transition, an out-of-equilibrium phenomenon with characteristics similar to equilibrium phase transitions (Jaiswal et al., 2016;Lin et al., 2015Lin et al., , 2014b.

Thermal systems
On the other hand, thermal fluctuations may play a role in materials with small enough ( 1 µm) 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 mechanically distinct systems. Even the waiting time between the preparation and the experiment matters, because the system's configuration evolves meanwhile, through activated events. The evolution of the mechanical properties with the time since preparation, usually making the system more solid, is called aging. 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 Σ max . It is observed in experiments (Divoux et al., 2011b) as well as numerical simulations (Patinet et al., 2016;Rottler and Robbins, 2005). In soft materials amenable to 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 (Doliwa and Heuer, 2003a,b;Goldstein, 1969). The whole configuration of the system (particle coordinates and, possibly, velocities) is considered as a 'state point' Γ that evolves on top of a hypersurface V (Γ) representing the total potential energy. Despite the high dimension of such a surface (proportional to the number N of particles), it can be viewed as a rugged landscape, with hills and nested valleys; the number of local minima generally grows exponentially with N ( 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).
External driving restricts the regions of the PEL that can be visited by the state point to, say, those with a (usually time-dependent) macroscopic strain γ. Mathematically, this constraint is enforced by means of a Lagrange multiplier, which effectively tilts V (Γ) into where Ω 0 is the volume of the system and Σ the macroscopic stress. The system's dynamics are then controlled by ∂V σ /∂Γ, instead of ∂V /∂Γ, which results in major changes, as we shall see next. Typically, driven systems respond on much shorter times than (quiescent) aging ones. Accordingly, some thermal systems may be treated as athermal, for all practical purposes. Nonetheless, interesting physical behavior emerges when the aging and driving time scales compete, either because temperature is high or because the driving is slow (Chattoraj et al., 2010;Johnson and Samwer, 2005;Rottler and Robbins, 2005;Shi and Falk, 2005;Vandembroucq and Roux, 2011).

C. 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 evolution of strain γ(t) with time (or stress Σ(t)) 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). This type of curves is not specific to the deformation of amorphous solids. It observed in all 'stick-slip' phenomena, in which the system is repeatedly 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 irregularities (defects) in the structure called dislocations -the elementary mechanism of plasticity. In the plastic flow of amorphous solids, potential energy V 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.
From the PEL perspective, the energy accumulation phase coincides with the state point smoothly tracking the evolution of the local minimum in the effective potential V σ (Γ, γ) [Eq.
(2)], as γ increases. Meanwhile, some effective barriers subside, until one flattens so much that the system can slide into another valley without energy cost. This topological change in V σ at a critical strain γ = γ c is a saddle-node bifurcation and marks the onset of a plastic event. For smooth potentials, close to γ c , the effective barrier height scales as (Gagnon et al., 2001;Maloney and Lacks, 2006) Note that the instability can be triggered prematurely if thermal fluctuations are present. In summary, in the PEL, deformation is a succession of barrier-climbing phases (elastic loading) and descents. 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 and molecules remains virtually invisible to direct experimental techniques? In the 1970s, inspiration came from the better known realm of crystals. As early as 1934, with the works of Orowan, Polanyi and Taylor, the motion of dislocations was known to be 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 at least fair to say that, should they exist, such defects would be more elusive than in crystals (we will come back to this question later in this chapter). In fact, the main inspiration drawn from research on crystals was of more pragmatic nature: Bragg and Nye (1947) showed that bubble rafts, i.e., monolayers of bubbles, could serve as upscaled models for 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 transposed to disordered systems by Argon and Kuo (1979), who used bidisperse bubble rafts as models for metallic glasses. Most importantly, they observed prominent singularities in the deformation: rapid rearrangements involving a few bubbles. Princen and Kiss (1986) suggested that the elementary rearrangement in these 2D systems was a topological change involving four bubbles, termed T1 event (see Fig. 4a).  (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.
These are the essential events whereby the irreversible macroscopic deformation is transcribed into the material structure. As such, they are the building bricks of EPM and we will often refer to them as 'plastic events' 1 . Since they must contribute to the externally imposed shear deformation, they will retain part of its symmetry and can thus be idealized as localized shear deformations, or 'shear transformations' (ST). Admittedly, in reality, their details do somewhat vary across systems (see below ), but, compared to the crystalline case, their strong spatial localization is a generic and remarkable feature.
At this stage, no doubt should be left as to the solid nature of the materials considered here. Our review is concerned with materials that are clearly solid at rest, and excludes systems at the fringe of rigidity, such as barely jammed packings of particles or emulsions. The latter are most probably governed by different physics, in which spatially extended rearrangements may occur under vanishing loading (Müller and Wyart, 2015).

Quantitative description
Although 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 D 2 min , a quantity that measures how nonaffine the local displacements around a particle are. More precisely, the relative displacements of neighboring particles between successive configurations are computed, and compared to the ones that would result from a locally affine deformation F ; D 2 min is the minimal square difference obtained by optimizing the tensor F . This quantity has been used heavily since then (Chikkadi and Schall, 2012;Chikkadi et al., 2011;Jensen et al., 2014;Schall et al., 2007), Generally speaking, very strong localization of events is observed at low enough shear rates, with spatial maps of D 2 min 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, yz if the applied strain is along xy in a 3D system) (Schall et al., 2007), the field of deviations of particle displacements with respect to a strictly uniform deformation, the count of nearest-neighbor losses (Chikkadi and Schall, 2012) or the identification of regions with large (marginal) particle velocities (Nicolas and Rottler, 2018), are also good options to detect rearrangements. Up to differences in their intensities, these methods were shown to provide similar information about STs 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 about the picture of clearly separated localized transformations. First, the validity of the binary vision distinguishing elastic and plastic regions has been challenged for hard particles, such as grains (Bouzid et al., 2015a). de Coulomb et al. (2017 thus recently contended that, as particles become stiffer, this binary picture fades into complex reorganizations of the network of contacts via cooperative motion of the particles. It is very tempting to relate these changes to the surge of delocalized low-energy excitations in emulsions  and packings of frictionless spheres (Andreotti et al., 2012;Wyart et al., 2005) as pressure is reduced and they approach jamming from above: These spatially extended and structurally complex soft modes are swept away upon compression and leave the floor to more localized excitations at higher pressure and energies.
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 induces; this imprint is revealed when one studies suitable stress or strain correlation functions (Chattoraj and Lemaître, 2013;Illing et al., 2016;Lemaître, 2014).

Variations
The foregoing quantitative indicators of microscale plasticity have brought to light substantial variations and differences between actual rearrangements and idealized STs. Even though EPM will generally turn a blind eye to this variability, let us shortly mention some of its salient features.
First, the sizes of rearranging regions vary 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 and therefore strongly tied to T g (Johnson and Samwer, 2005;Zhang et al., 2017b). Albaret et al. (2016) proposed a detailed numerical characterization of plastic rearrangements in atomistic models for amorphous silicon by fitting the actual particle displacements during plastic events with the elastic displacement halos expected around spherical STs, that is, by finding the size and spontaneous deformation of the inclusion which, upon embedding in an elastic medium, generates displacements that best match the observed ones (see Sec. I.D).
Although rearrangements seem to have a typical linear size, around 3 • A, they found that the most robust quantity is actually the product of with the inclusion volume V in . Furthermore, the mean strain Tr( )/3 is either positive or negative depending on the specificities of the implemented potential (thus evidencing either a local dilation or a local compression) and represents only about 5% of the shear component, which confirms that shear prevails in the transformation. The orientations of the STs, i.e., the directions of maximal shear, were studied in greater detail in a more recent work, in 2D, where a fairly broad distribution of these orientations around the macroscopic shear direction was reported (Nicolas and Rottler, 2018). Finally, Albaret et al. (2016) 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 proves that localized plastic rearrangements surrounded by an elastic halo are the unique elementary carriers 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 STs 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, in thermal Comparison between different indicators based on the local structure to predict future plastic rearrangements. (a) Contour maps of the participation ratio in the 1% softest vibrational modes for two numerical samples of a binary metallic glass (Ding et al., 2014). (b) Correlation between the locations of future rearrangements and diverse local properties in an instantaneously quenched binary glass model. The following properties have been considered: local yield stress (τy), participation in the soft vibrational modes (PF), lowest shear modulus (2µ I ), local potential energy (PE), short-range order (SRO), local density (ρ). From (Patinet et al., 2016). (c) Identification of soft particles (thick black contours) by a machinelearning algorithm (SVM), in a compressed granular pillar. Particles are colored from gray to red according to their actual nonaffine motion (D 2 min value). From (Cubuk et al., 2015).
systems these STs may nonetheless be predominantly activated by thermal fluctuations (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 an ST, i.e., r → r+ ·(r − r c ) (where r generically refers to a particle position); incidentally, this is the major reason why the observable D 2 min detects plastic rearrangements.

Structural origins of rearrangements
What determines a region's propensity to rearrange? Local 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 V f per particle, notably by using (contested) experimental evidence from polymeric materials (Batschinski, 1913;Doolittle, 1951;Fox Jr and Flory, 1950;Williams et al., 1955). The idea that local variations of V f 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 STs 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. Beyond this particular example, the question of how to decipher the weakness of a region from its local microstructure remains largely open.
Looking beyond locally available information, the linear response of the whole system has also been considered, with the hope that linearly soft regions will 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). One should mention that, albeit a local property, the local shear modulus is best evaluated if the response of the global system is computed. Indeed, approximations singling out a local region and enforcing an affine deformation of the outer medium are sensitive to the size of the region and overestimate the shear modulus, because they hinder nonaffine relaxation (Mizuno et al., 2013).
Focusing on vibrational properties, Brito and Wyart (2007) and Widmer-Cooper et al. (2008) provided evidence that in hard-sphere glasses as well as in supercooled liquids the particles that vibrate most in the M lowest energy modes (i.e., those with a high participation fraction in the M softest modes, where M is arbitrarily fine-tuned), are more likely to rearrange (note that for hard spheres the vibrational modes were computed using an effective interaction potential). 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 .
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., 2015Schoenholz 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 D 2 min , 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 highdimensional 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.

D. 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 shear 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, if the sample is large compared to the ST size, the impact of microscopic events on the macroscopic response could not be explained without collective effects and avalanches involving a large number of plastic events. 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 ST. 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 = ∇u+∇u 2 , where u stands for the displacement. Recall that a shear deformation, say (r ≈ 0) = 0 1 1 0 in two dimensions (2D), consists of a stretch along the direction θ = π 4 [π], in polar coordinates, and a contraction along the perpendicular direction. The induced displacement field u simply mirrors this symmetry, with displacements that point outwards along θ = π 4 [π] and inwards along θ = 3π 4 [π]. This leads to a dipolar azimuthal dependence for u and a four-fold ('quadrupolar') one for its symmetrized gradient . More precisely, by imposing mechanical equilibrium on the stress Σ, viz., 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 ST 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 (b-c) the central blocks are artificially colored.
in an incompressible medium (∇ · u = 0) with a linear elastic law, Σ dev ∝ dev (where the superscript denotes the deviatoric part), Picard et al. (2004) derived the induced strain field in 2D, Here, only one of the strain components is expressed, but the derivation is straightforwardly extended to a tensorial form (Budrikis et al., 2017;Nicolas and Barrat, 2013a). Experiments on colloidal suspensions (Jensen et al., 2014;Schall et al., 2007) 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 ST. Yet, the response can be calculated in the near field following Eshelby (1957)'s works, by modeling the ST as an elastic inclusion bearing an eigenstrain , i.e., spontaneously evolving towards the deformed configuration . This handling adds near-field corrections to Eq. (4), whose analytical expression is derived by Jin et al. (2016) and Weinberger et al. (2005) for an ellipsoidal inclusion in 3D on the basis of a method based on Green's function, which is probably more accessible than Eshelby (1957)'s original paper [see Jin et al. (2017) for the 2D version of the problem].
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 ST 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 then extended to the deformation of a spherical inclusion in 3D, and to the nonlinear regime, by Priezjev (2015). (Also see  for the nonlinear consequences of artificially triggered STs in a 2D glass).
Besides elastic heterogeneity, further deviations from the Eshelby response result from the difference between an actual plastic rearrangement and the idealized ST 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; only the former would quantitatively obey Eshelby's formula. It might be that the dilational component of the rearrangement, discarded in the ideal ST, is important in the thermal regime.
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)
A. 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 'mesoscopic' models for amorphous plasticity, or sometimes 'discrete automata'. To mimick the basic flow scenario, the material is split into 'mesoscopic' blocks, presumably of the typical size of a rearrangement. These blocks are loaded elastically (R1) until a yield condition is met (R2), at which point their stress relaxes plastically and is redistributed to other blocks via the emission of an elastic stress field (R3); finally, they become elastic again (R4). Accordingly, EPM hinge on the following set of predefined rules (Rodney et al., 2011): R1. a (default) elastic response of each mesoscopic block, R2. a local yield criterion that determines the onset of a plastic event (n : 0 → 1), R3. a redistribution of the stress during plasticity that gives rise to long-range interactions among blocks, R4. a recovery criterion that fixes the end of a plastic event (n : 1 → 0), where the activity n is defined as n = 0 if the block is purely elastic, and n = 1 otherwise.
To make it more concrete, consider the class of EPM pioneered by Picard et al. (2005). To fix rules R2 and R4, which define when a region switches from elastic to plastic and conversely, the model specifies the rates governing the transitions n : 0 ↔ 1, whereas R1 and R3 are implemented via the following equation of evolution for the stress σ i carried by block i (where i is a d -dimensional vector denoting the lattice coordinates of the block), Here, the stress incrementσ i per unit time is the sum of three contributions. First, the externally applied simple shear contributes a uniform amountσ = µγ following Hooke's law, with µ the shear modulus andγ the shear rate. Second, nonlocal plastic events (at positions j with n j = 1) release a stress of order σ j over a time scale τ , which is transmitted by the elastic propagator G ij . Third, if the site is currently plastic (i.e., if n i = 1), its stress σ i relaxes. Spontaneously, it would do so at a rate τ −1 , but, owing to the constraint of the surrounding elastic medium, the efficiency of the process is not optimal and the rate drops to |G 0 |/τ , with 0 < |G 0 | < 1. (During this relaxation, elastic stress increments are still received, as are they in Maxwell's visco-elastic fluid model.) In many regards, the stress evolution described by Eq. (5) is overly simplified: It focuses on one stress component whereas stress is a tensor, it assumes instantaneous linear transmission of stress releases in a uniform medium, the local stress relaxation rate is constant, etc. In Sec. III, we will detail how these approximations have been relaxed in part in some EPM and we will see how the phenomenological set of rules R.1-4 can emerge in the framework of Continuum Mechanics (Sec. IV.F). But, for the time being, we find it important to delve into the philosophy of these models and reflect on their goals.
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. Central is the assumption that most microscopic details are irrelevant for the main rheological properties and that only a few relevant parameters or processes really matter. Several reasons could be put forward to favor their use over more realistic modeling approaches, e.g., 1. to assess the validity of a theoretical scenario and identify the key physical processes in the rheology, 2. to provide an efficient simulation tool giving access to (otherwise inaccessible) large statistics or long-time runs, 3. to facilitate 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).
That distinct EPM highlight distinct physical ingredients seems 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. Less intuitive is perhaps the role of the experimental conditions and the observables of interest in determining the important physical ingredients, but here come a couple of enlightening examples: There is no need to keep track of previous configurations (e.g., past yield stresses) to study steady shear, whereas this might be crucial for oscillatory shear experiments in which the system cycles between a few configurations (Fiocco et al., 2014). Also, potentially universal aspects of the yielding transition show little to no sensitivity to the precise EPM rules, while the latter affect the details of the flow pattern. 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 properties of interest.
In the following, we list the physical processes that are put in the limelight in diverse EPM and indicate for what type of materials and in what conditions they are of primary importance.

B. Thermal fluctuations
How relevant are thermal fluctuations and the associated Brownian motion of particles? This question brings us back to the distinction between thermal materials and athermal ones exposed in Section I.B.
It is widely believed that thermal fluctuations largely contribute to the activation of plastic events in metallic and molecular glasses, as well as in suspensions of small enough colloids. For a suspension of 1.5 µm-large colloids, Schall et al. (2007) thus argue on the basis of an estimate of the activation energy that transformations are mostly thermally activated, 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) set a yield rate based on a stress-biased Arrhenius law for thermal activation, viz., where ν 0 is an attempt frequency, V σ is the height of the (smallest) potential barrier hindering the rearrangement. Recalling from Eq. (2) that the potential is tilted by the stress σ, i.e., V σ = V − Ω 0 σγ, one immediately recovers the expression of ν(σ) used by Eyring (1935) to calculate the viscosity of liquids if σ and γ are treated as independent parameters. To the contrary, if the local stress and elastic strain are related by Hooke's law, viz., σ ∝ γ, one finds the γ 2 -scaling of the tilt used, e.g., in Sollich et al. (1997)'s Soft Glassy Rheology model (Sec. IV.D.1).
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 (Chen et al., 1991;Hébraud and Lequeux, 1998) favor a a binary yield criterion, viz., where σ y is the local stress threshold for yielding; a deterministic yield criterion is recovered in the limit ν 0 → ∞. Incidentally, note that, in this 1D tilt picture, the existence of favorable directions in the PEL is handled somewhat light-heartedly. Indeed, in a high-dimensional PEL, the direction in which the loading pushes the system may differ from that of the lowest saddle point; this fact may be particularly relevant for small systems, whose granularity is more apparent.
As far as rheology is concerned, thermal activation can be neglected if it does not trigger rearrangements much below the local yield stress σ y . By requiring that the activation rate ν(σ) in Eq. (6) match the driving rateγ/γ y (rescaled by the yield strain γ y ) at a stress close to σ y , we find that the athermal approximation is conditioned on This criterion bears some resemblance with the limit of large Péclet number Pe ≡γ a 2 /D, where a is the particle size and D is the single-particle diffusivity in the dilute limit (Ikeda et al., 2013), but duly takes into account the cage constraints which restrict diffusion in a dense system. 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 observation of creep by Divoux et al. (2008) in a granular heap submitted to cyclic temperature variations does not underscore a possible importance of thermal activation, 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 expedite 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.

C. Driving
Suppose that the material deforms under external driving; how important are the specific driving conditions?

Driving protocol
Numerical simulations have mostly considered strain-controlled (fixedγ), rather stress-controlled (fixed Σ), situations (see Sec. I.B). For strain-driven protocols, the stress redistribution (R3) operated by the elastic propagator G in EPM keeps the macroscopic strain fixed. Meanwhile, the elastic response (R1) is generally obtained by converting the macroscopic driving into local stress increments µγ(t)dt, where dt is the time step andγ(t) is the current macroscopic strain rate. EPM often focus on steady shear situations, in which caseγ(t) = cst. But time-dependent driving protocolsγ = f (t) [or Σ = f (t)] are also encountered, in particular step shearγ(t) = γ 0 δ(t) and oscillatory sheaṙ γ(t) = γ 0 cos (ωt), which gives access to linear rheology for small γ 0 .
Stress-controlled setups have received less attention in the frame of EPM but examples can be found in (Homer and Schuh, 2009;Jagla, 2017b;Lin et al., 2015Lin et al., , 2014bPicard et al., 2004). In this case, the 0-Fourier mode of the elastic propagator G is adjusted so that G keeps the macroscopic stress constant (see Sec. III.C.2). In creeping flows subjected to σ(t) = cst < σ y ,γ(t) eventually decays to zero, often as a power law (Leocmach et al., 2014). Creep is further discussed in Sec. VIII.

Symmetry of the driving
Plastic events are biased towards the direction of the external loading (Nicolas and Rottler, 2018). 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 ∇u = 0 γ(t) 0 0 (in the linear approximation in 2D), one may settle with the xy component of the linear strain tensor = ∇u+∇u 2 .
It has the same principal strains (eigenvalues) ± γ(t) /2 as pure shear, ∇u = 0 γ(t) /2 γ(t) /2 0 , but involves a rotational part ω = 0 γ(t) /2 − γ(t) /2 0 , whereas the latter is rotationless. 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., σ(t) = σ(t) 1 0 0 0 in the bulk, with σ(t) < 0) and tension (σ(t) > 0) are often performed (Priezjev and Makeev, 2017). Even though in several of these situations the macroscopic loading 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 (Budrikis et al., 2017;Nicolas et al., 2014b), and the effect of dimensionality to be weak . The reader is referred to Sec. VI 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). For these heterogeneous driving conditions, EPM has emerged as a promising alternative to atomistic simulations. Nevertheless, further upscaling is needed whenever the length scale associated with the heterogeneous driving is very large compared to the particle size.

D. Driving rate and material time scales
To resolve the flow temporally, the simplest approach is an 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. They are particularly efficient in stress-controlled slow flows, insofar as the long elastic loading phases without plastic events are bypassed: The activation rate ν i is calculated for every block i using a refined version of Eq. (6) and the time lapse before the next plastic event is deduced from the cumulative rate ν = i ν i (Homer and Schuh, 2009).
In various models the finite duration of plastic events plays a major role in theγ-dependence of the rheology Martens et al., 2012;Nicolas et al., 2014a;Picard et al., 2005) or in the intrinsic relaxation of the system (Ferrero et al., 2014). Suppose that, under slow driving, a rearrangement takes a typical time 2 τ pl . For overdamped dynamics, one expects this time scale to be the ratio between an effective microscopic viscosity η eff and the elastic shear modulus µ (Nicolas and Barrat, 2013a), viz., while for underdamped systems τ pl is associated with the persistence time of localized vibrations. If τ pl competes with the driving time scale τγ ≡ γ y /γ, where γ y 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).
At even lower driving rates, one reaches a regime where τ pl τγ and individual rearrangements become insensitive to the driving. But the latter may still affect avalanches of rearrangements (i.e., the series of plastic events that would still be triggered by an initial event, were the driving turned off). Indeed, since the size of avalanches is expected to diverge in an athermal system in the limits of vanishing shear rate and infinite system size L, their duration τ av (L), bounded below by the signal propagation time between rearranging regions, may become arbitarily large asγ → 0. While most EPM turn a blind eye to the delays due to shear wave propagation, some works have bestowed them a central role in the finite shear-rate rheology (Lemaître and Caroli, 2009;Lin et al., 2014b) and there have been endeavors to represent this propagation in a more realistic way in EPM (Karimi and Barrat, 2016;Nicolas et al., 2015) (see Sec. III.D). Sections VI.D.3 and VII.B will provide more details on the influence of vanishingly small shear rates on the flow curve.
The quasistatic limit is reached when τ pl τγ → 0 and τ av (L) τγ → 0 (8) and the athermal criterion of Eq. (7) is satisfied. In that case, the material remains in mechanical equilibrium at all times and its trajectory in the PEL is rate-independent. Atomistic simulations can then be simplified by applying a small strain increment at each step and letting the system relax athermally to the local energy minimum (Maloney and Lemaître, 2004). The EPM counterparts of these quasistatic simulations are called extremal or quasistatic models and have been studied intensively (Baret et al., 2002;Jagla, 2017b;Lin et al., 2014b;Talamali et al., 2012). In these 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. The material time scales are then naturally brushed aside, while connection to real time is lost.

E. Spatial disorder in the mechanical properties
Glasses, and more generally amorphous solids, are mechanically heterogeneous. Indeed, there have been both experimental and numerical reports on the heterogeneity of the local elastic moduli (see Fig. 7) and the energy barriers on the mesoscale (Tsamados et al., 2008;Zargar et al., 2013). 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 plays a central role in others (Langer, 2008;Sollich et al., 1997). 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. Sollich et al. (1997)'s Soft Glassy Rheology model introduces exponentially distributed energy barriers; this translates into an even broader distribution of activation rates via Eq. (6) and leads to a transition from Newtonian to non-Newtonian rheology as temperature is reduced (see (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 units. From (Tsamados et al., 2009). (b-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 (PdCuSi) and (c) a crystal (SrTiO3). The latter clearly appears to be mechanically more homogeneous. The radius of contact is of order 10 nm. From (Wagner et al., 2011).
Sec. IV.D.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. Finally, amorphous composite materials, i.e., materials featuring meso/macro-inclusions of another material, can be modeled 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: It only affects the coefficients of the rheological law, and not the functional shape.
F. Spatial resolution of the model On a related note, how important is it to spatially resolve an EPM? In what cases can one 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 even on some critical properties (Lin et al., 2014b;. But a mean field analysis could suffice in many other situations. For example, Martens et al. (2012) showed that the flow curve obtained with their spatially resolved EPM at finite shear rates can be predicted on the basis of mean-field reasoning, whereas spatial correlations and avalanches are thought to impact the macroscopic stress at vanishing shear rates . 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 r −d in d dimensions), mean-field arguments should generically hold in amorphous solids (Dahmen et al., 1998(Dahmen et al., , 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 (Lin et al., 2014a;Nicolas et al., 2014b;Talamali et al., 2011). Meanfield 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 nonmean-field exponents were found for the power-law distribution of avalanche sizes (Budrikis and Zapperi, 2013;Lin et al., 2014b;Talamali et al., 2011). This question is addressed in greater depth in Sec. VI.
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.

G. Bird's eye view of the various models
To conclude this section, some of the main EPM are classified in Table I.

Marmottant and Graner (2013)
Overdamped evolution in a periodic potential; pl. events of finite duration; no stress redistribution foams Symbols -Barrier distribution: single value / distributed (exponentially, unless otherwise specified).

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 (Budrikis and Zapperi, 2013;Martens et al., 2012). 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 relate this level of sophistication with the nature of the developments that were sought.
A. Sandpile models and first-neighbor stress redistribution EPM owe much to the quake-ridden scientific grounds on which they burgeoned at the beginning of the 1990s, marked by the advent of seminal models for earthquakes and avalanches. As a paradigmatic example for earthquakes, 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 its length. 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 U (x, t) at time t of the material at position x 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 vt due to the motion of the plate as well as the displacement −U contribute to a shear term. Finally, φ(U ) is a velocity-dependent frictional term. Had Coulomb's law of friction been used, it would have been constant for |U | = 0, but the original model assumed velocity weakening, i.e., a decrease of |φ(U )| with |U |. AtU = 0, the function φ is degenerate, which allows static friction to exactly cancel the sum of forces on the rhs of Eq. (9), so the blocks remain pinned at a fixed position U until the destabilizing forces ξ 2 ∂ 2 U ∂x 2 + vt 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 (Carlson et al., 1994;Otsuka, 1972). Important in the above model is the effect of the pinning force φ atU = 0. It entails that the destabilizing action caused by the depinning of a site (via the diffusive term in Eq. (9)) 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 (i, j), gets too high with respect to its neighbors (say, for convenience, whenever σ i,j 4), some grains at its top are transferred to the neighboring columns, with the following discharge rules in 2D: 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 concept of self-organized criticality emerged (Bak et al., 1987). According to the latter, 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 (Bak and Tang, 1989;Carlson and Langer, 1989;Ito and Matsuzaki, 1990;Sornette and Sornette, 1989) 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 E is the energy release, in a given region obeys the power law relation, log P (m m 0 ) −bm 0 + cst, where b 0.88, or equivalently For the sake of accuracy, we ought to say that there exist several earthquake magnitude scales besides that of Eq. (11). They roughly coincide at not too large values; in fact, M e is not the initial Richter scale. More importantly, the value of the exponent b ∈ [0.8, 1.5] depends on the considered earthquake catalog, and notably on the considered region. For sandpile-like models, various exponents have been reported: τ ≈ 1 in 2D and τ ≈ 1.35 in 3D, with no effect of disorder of the yield stresses (Bak and Tang, 1989), whereas the exponent for the mean-field democratic fiber bundle close to global failure is τ = 3 /2 (see Sec. IX.C). More extensive numerical simulations led to the values τ 1.30 (Lübeck and Usadel, 1997), or τ 1.27 (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. (10) 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 p(E) 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 transfer 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.D. 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.

B. 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 variants thereof.
First, in the spirit of the classical description of a solid as an assembly of particles confined to their positions by interactions with their neighbors, the material was modeled as a system of blocks connected by "springs" of stiffness κ and potential energy where u i 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.C. 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 F with the displacement u; this evolution starts with a phase of linear increase, followed by a peak and a smooth decline until global failure. The F (u) curves for different linear lattice sizes L roughly collapsed onto a master curve if F and u were rescaled by L − 3 /4 . 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. In passing, note the subtle difference with respect to the double force-dipole (often called quadrupole) describing irrotational local shear (see Sec. III.C), which has a distinct anisotropy. Contrary to Hansen et al. (1989), they allowed broken springs to instantly regenerate to an unloaded state, after 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 τ = 1.4 was reported in 2D, in semi-quantitative agreement with the Gutenberg-Richter earthquake statistics.

C. Elastic propagators
After the discrete vision promoted by block-and-spring models, let us momentarily turn to a continuous description of the amorphous solid. The free energy of the material can be expressed in terms of the displacement field u(r) as where the free energy density ψ depends on the strain tensor (r) = ∇u(r)+∇u(r) 2 (in the linear approximation), because rigid transformations cost no energy. In the quiescent system, ψ reaches its minimum in the reference strain state (r) = 0, viz., where σ = dψ/d is the stress. But a plastic rearrangement taking place at, say, r = 0 will shift the reference state to = pl in a small region around 0. Because of its embedding, this region cannot deform freely. Therefore, the plastic strain pl (more generally known as eigenstrain) will induce a nonlocal elastic response in the surrounding medium, which was worked out by Eshelby (1957) for ellipsoidal inclusions in a linear elastic solid. Picard et al. (2004) simplified the calculation of the elastic response by supposing that the shear transformation (ST) has vanishing linear size a → 0. Prior to the ST, the material is linear elastic, viz.,

Pointwise transformation in a uniform medium
where p is the pressure, and incompressible, ∇ · u = 0. After the ST, the reference state is shifted and Eq. (13) turns into where the source term generated by the plastic strain at the origin reads f (r) = −2µ∇ · pl a d δ (r) . The solution of Eq. (15) is well known in hydrodynamics and involves the Oseen-Burgers tensor O(r) = 1 8πµr I + r⊗r r 2 in 2D, with I the identity matrix, viz., In the unbounded 2D plane, setting coordinates such that pl = 0 0 0 0 , the response to f (r) in terms of xycomponent of the stress reads where (r, θ) are polar coordinates. This field is shown in Fig. 6c. Reassuringly, in the far field it coincides with the response to a cylindrical Eshelby inclusion. As a short aside, let us mention a variant to these calculations, which underscores the connection with deformation processes in a crystal. This variant is reminiscent of Eshelby's a cut-and-glue method, whereby 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 pl (r). Switching to Fourier space (r ↔ q ≡ (q x , q y )), the counterpart of Eq. (17) is then Note that the frame is sometimes defined such that pl = . 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 q x , q y = 2πn L , n ∈ −L 2 , . . . , L 2 . Besides, because of dissipative forces, quantified by an effective viscosity η eff , the strain rate˙ in the ST cannot be infinite and a rearrangement will last for a finite time τ pl ∼ η eff/µ (see Sec. II.D). Therefore, in each numerical time step, the plastic strain pl implemented in Eq. (18) will be the strain increment δ pl during that step. This amounts to saying that, locally, dissipative forces make the rearrangement gradual, while stress is redistributed instantaneously to the rest of the medium (because of the assumption of mechanical equilibrium), so that there is no time dependence in the elastic propagator in Eq. (18).

Technical issues with pointwise transformations and possible remediations
The idealized elastic propagator in Eq. (18) brings on some technical issues. Firstly, its slow (∝ r −d ) radial decay 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 G (q) near q = 0. 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 G (q = 0) = 0 in a stress-controlled system and G (q = 0) = −1 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).
Secondly, 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 G sum (q) to coincide with the solution G ∞ (r) for an unbounded medium. However, the high-frequency components in G (q), due to the spurious singularity of G ∞ (r) at r = 0 [Eq. (18)], make the periodic sum diverge. In practice, wavenumbers outside the first Brillouin zone ] − π, π] d are plainly discarded, which comes down to solving Eqs. (14)-(15) 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 G(q) associated with zero energy. These so called soft modes (or null modes) pl , satisfying ∀q, G (q) pl (q) = 0, cost no elastic energy; their deployment is thus favored by the dynamics (Tyukodi et al., 2016b). Their importance is further explained in Sec. IV.C. It turns out that the eigenmodes of G(q) in Eq. (18) are the Fourier modes (plane waves); among these, the soft modes are those with wavevectors q making an angle ± π /4 with respect to the principal direction of the plastic strain tensor pl .
In particular, under simple shear with velocity direction x and velocity gradient along y, the emergence of a uniform shear band along x should produce no elastic stress in the medium, at least if such a band emerges uniformly. However, misaligned lattice axes (not directed along x or y) are incompatible 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). The scalings of avalanche sizes, however, seem to be mostly insensitive to these details. Indeed, these details do not affect the long-range interactions between blocks.
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 FIG. 9 Average displacement field induced by an ST 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) ∆t = 2, (b) ∆t = 10, (c) ∆t = 1000. Red hues indicate larger displacements. Adapted from (Nicolas et al., 2015) tensorial form; but it turns out that, for setups with uniform loading, the tensorial extension has virtually no effect (Nicolas et al., 2014b). The statistics of avalanches of plastic events found in tensorial models and in scalar models are similar, up to tenuous differences: 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;. Another refinement consists in taking into account the possible anisotropy of the solid. Cao et al. (2018) argued that anisotropic shear moduli (µ 2 = µ 3 ) may result from the shear softening of well annealed solids just before their macroscopic failure. The displacement field induced by an ST, where θ = 0 denotes the principal direction of positive stretch, then tends to concentrate along 'easy' axes as the anisotropy parameter δ = (µ 3 − µ 2 )/(µ 3 + µ 2 ) increases. Translational invariance may also be broken, if the system is confined between walls, as in a microchannel, instead of being periodic. If there is no slip at the wall, a method of images allows the derivation of the elastic propagator for the bounded medium (Nicolas and Barrat, 2013a;Picard et al., 2004). Plastic events are found to relax stress faster, for a given eigenstrain, if they occur close to the walls. Such changes in boundary conditions 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).

D. Approaches resorting to Finite-Element methods
Albeit computationally more costly, Finite-Element (FE)-based computations of stress redistribution overcome some limitations of the foregoing approaches and offer more flexibility. The FE method solves the continuum mechanics equation associated with the free energy of Eq. (12) by interpolating the strain and stress σ within each element of a meshgrid from the values of the displacements and point forces at the nodes of the element. As far as EPM are concerned, the default elastic response of each block is generally assumed linear, so that the free energy density in Eq. (12) reads ψ( ) = 1 2 · C · and σ = C · . If mechanical equilibrium is maintained at all times, as in Eq. (13), the response to an ST is obtained by equilibrating the elastic stress C pl that it releases. Using a triangular mesh refined around the ST-bearing element, Sandfeld et al. (2015) demonstrated that the computed stress field coincides with the elastic propagator of Eq. (17) in the limit of a pointwise ST. 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 elements' 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 Sec. VI).
In an earlier endeavor Homer and Schuh, 2009), each ST 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.  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 STs, 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 Sec. V.C.2).
The capabilities of FE methods were further exploited by Nicolas et al. (2015) to go beyond the assumption of elastic homogeneity and capture the fluctuations caused by elastic disorder, notably those evidenced in the response to an ideal ST in MD simulations of a binary glass (Puosi et al., 2014). For this purpose, stiffness matrices C were measured locally in that model glass, in mesoscale regions of 5 particles in diameter. The authors found a broad distribution p(µ) of local shear moduli, with a relative dispersion of around 30%, and marked anisotropy (i.e., one direction of shear being much weaker than the other one). The dispersion of p(µ) explains the fluctuations observed in MD. Indeed, an FE-based model in which the shear modulus of each element was randomly drawn from p(µ) displayed comparable fluctuations in the response to an ideal ST, triggered by suitably moving the nodes of a chosen element. On the other hand, accounting for anisotropy was less critical. Furthermore, in the FE description inertial and viscous terms were restored in Eq. (15). This gives access to the transient elastic reponse, involving the propagation of shear waves. Exploiting this opportunity, Karimi et al. (2017) analyzed the effect of inertia on the avalanche statistics and compared it with results from atomistic simulations. (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 MD when the friction coefficient is varied.

IV. MEAN-FIELD TREATMENTS OF MECHANICAL NOISE
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. The rate ζ i (t) of stress increment experienced by a given block (say, site i) at time t is then a sum of contributions from many sites, i.e., using Eq. (5), ζ i (t) = j =i n j G ij σj τ , where n j denotes the plastic activity of site j. 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 n i , the mechanical noise signal {ζ i (t)} 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 i over time (and over i if ergodicity is broken). This shows the central role of {ζ i (t)} 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.

A. Uniform redistribution of stress
The mechanical noise can be split into : • a constant background ζ i , which contributes to a drift term µγ eff i ≡ µγ + ζ i in Eq. (20), and • zero-average fluctuations δζ i (t).
Owing to the infinite range and slow decay of the elastic propagator (∝ r −d in d-dimensional space, see Sec. I.D), site j is significantly coupled to a large number of other sites. This large connectivity has led some researchers to brush aside fluctuations in favor of the average drift term. Along these lines, in the framework of Picard's EPM, which features a constant rate τ −1 of yield above a uniform threshold and a constant rate τ −1 res of elastic recovery, viz., Martens et al. (2012) averaged Eq. (20) over time and found an analytical expression for the flow curve, which reproduces the simulation results to a large extent, at least at reasonably large shear rates. It also correctly predicts the destabilization of the homogeneous flow leading to shear-banding for a range of model parameters, in particular at large τ res . In fact, the neglect of fluctuations would be rigorously justified if the system were infinite and the propagator G were positive. The latter criterion is for instance fulfilled in a simple quasistatic model in which sites yield past a threshold σ y and redistribute the released stress (δσ i ) uniformly to the other N − 1 ≈ N sites (Dahmen et al., 1998), viz., The simplicity of the model allows analytical progress. A first approach consists in treating the distances x i = σ y − σ i to the threshold σ y as independent variables in the system and sorting them in ascending order (i → 1, 2, . . .). An avalanche will persist as long as the stress increment δσ1 N due to the yielding of the most unstable site suffices to make the second most unstable fail, viz., δσ1 N > x 2 . Using an argument along these lines in a model featuring disorder in the yield thresholds (σ y → σ y,i ) and post-failure weakening (i.e., when site i yields, the threshold is restored to a lower value σ y,i (t + 1) < σ y,i (t)), 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. 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). We will come back to this point in Sec. VI.

Deviations from uniform mean field
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). The ratio of fluctations over mean value of the stress increments, estimated by Nicolas et al. (2014b) in an EPM, diverges at low shear ratesγ, which points to the failure of the meanfield theory, according to Ginzburg and Landau's criterion. 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;. For instance, in these simulations, nearγ → 0, the distribution of avalanche sizes S follows a power law P (S) ∼ S −τ with an exponent τ that deviates from the τ = 3 /2 value predicted by mean field (see Sec. VI for details).

The Hébraud-Lequeux model
To improve on the hypothesis of a constant mean fieldγ eff , fluctuations of the mechanical noise need to be accounted for. In the crudest approximation, they can be substituted by random white noise ζ (w.n.) (t), with ζ (w.n.) = 0. This turns Eq. (20) into a biased Brownian walk for the local stresses, in the elastic regime n i = 0. Hébraud and Lequeux (1998)'s model was developed along these lines. The ensuing stochastic equation (Eq. (20) with ζ i (t) →γ eff +ζ (w.n.) (t) and τ → 0) can be recast into a probabilistic Fokker-Planck-like equation operating on the distribution P (σ, t) of local stresses σ, viz., where the diffusive term D ∂ 2 P ∂σ 2 on the rhs arises from the fluctuations acting on σ i , with a coefficient D(t) assumed to be proportional to the number of plastic sites Γ(t) ≡ τ −1 c |σ |>σy P (σ , t)dσ , viz., D(t) = αΓ(t). The first term on the rhs of Eq. (22) is a drift term, which amalgamatesγ eff withγ; the last two terms correspond to the failure of overloaded sites (above σ y ) on a timescale τ c and their rebirth at σ = 0 due to stress relaxation. The resulting mean-field equations can be solved in the limit of vanishing shear ratesγ (Agoritsas et al., 2015;Olivier and Renardy, 2011). For a coupling constant α < 1 /2, diffusion vanishes at low shear rates, with D ∝γ, a yield stress Σ y > 0 is obtained and the average stress obeys Σ Σ y + kγ 1 /2 , with k > 0, in the low shear rate limit. For α > 1 /2, the system behaves like a Newtonian liquid.

Fraction of sites close to yielding
The diffusive term introduced in Eq. (22) impacts the distribution of sites close to yield, i.e., at distances x ≡ |σ| − σ y 1 from the yield threshold σ y . On these short distances, or, equivalently, in the limit of short time scales ∆t, the back-and-forth diffusive motion over typical distances ∝ √ ∆t prevails over the drift in the random walk. Therefore, forγ → 0, determining the distribution P (x) is tantamount to finding the concentration of Brownian particles near an absorbing boundary at x = 0 (yielding): The well-known solution is a linear vanishing of the concentration near x = 0, viz., P (x) ∼ x for x ≈ 0 (Lin et al., 2014a;. This result ought to be compared with P (x) ∼ x 0 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 [v ∝ F − F c ) β with β < 1] and the flow of disordered solids [γ ∝ Σ − Σ Y ) β with β > 1 generally].

C. 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 II.

Uniform mean field
Neglecting fluctuations in the constant mean-field approach makes sense in the drift-dominated regime, i.e., when |γ eff |∆t | ∆t 0 δζ(t )dt | on the considered time window ∆t, with the notations of Sec. IV.A. With interactions that change signs, this excludes vanishing shear rates or too small ∆t. But at higher shear rates, this approach appears to correctly predict the avalanche scaling exponents in the EPM studied by .

White-noise fluctuations
Complemented with Gaussian fluctuations, the approximation is valid beyond the drift-dominated regime. In fact, if the distribution of global mechanical noise δζ (i) is Gaussian-distributed and (ii) has no significant time correlations, the noise fluctuations δζ can be replaced by Gaussian white noise in Eq. (20) . Note that, if a single plastic event releases stress fluctuations with a distribution w 1 such as w 1 (δζ) ∝ |δζ| −1−k with k > 2, condition (i) will be fulfilled as soon as condition (ii) is (i.e., events are uncorrelated).
Provided that the mechanical noise fulfils the above criteria (i) and (ii), all models based on similar rules for plasticity thus fall in the universality class of the Hébraud-Lequeux model in the limit of large systems. In particular, for coupling parameters α such that the diffusivity D(t) goes to zero atγ → 0, their flow curves will follow a Herschel-Bulkley behavior Σ = Σ y + kγ n with n = 1 /2 in the low shear rate limit. This holds true in the presence of disorder on the local yield thresholds σ y (Agoritsas et al., 2015) and for plastic events that do not relax the local stress strictly to zero, but to a low random value . On the other hand, should the shear modulus of elastic blocks or the relaxation time depend onγ, the exponent n will deviate from 1 /2 .

Heavy-tailed fluctuations
However, the decay of the elastic propagator as G ∼ cos 4θ r d casts doubt on the Gaussian nature of the random stress increments δζ and would rather suggest a broad density function for the mechanical noise, in the limit of sparse plastic events, with an upper cut-off δζ M proportional to the volume of a rearranging region. For such a heavy-tailed distribution, the biased random Brownian walk of σ j is replaced by a Lévy flight of index k = 1 for σ j . 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 yielding P (x). To be explicit, their model was based on plastic yielding above a uniform yield strain γ y , which resets the local stress to zero and increments the stresses at other sites by random values drawn from w 1 (δζ). Somewhat surprinsingly, the use of a Gaussian distribution w 1 gave better agreement with quasistatic atomistic simulations than a heavy-tailed distribution, in terms of the scaling of avalanches with system size. Further insight is gained by understanding that the distribution w 1 in Eq. (23) describes the instantaneous stress released by a single event, whereas a material region will yield under the cumulated effect of a sum of such contributions. Assuming that this sum is random,  arrived at the following probabilistic equation , v > 0 drives sites towards their yielding point under the action of stress, and the distribution w δγ accounts for the sum of stress releases (drawn from w 1 (δζ)) taking place over δγ: w ∼ w 1 if k < 2; otherwise, it is a Gaussian. This model is a variant of the Hébraud-Lequeux model which differs from it in that (i) broad distributions w 1 (δζ) are allowed, (ii) time is measured in terms of plastic strain γ, and (iii) sites that cross the yield threshold σ y = 1 abruptly relax to 0 as γ is incremented. In this framework, the authors derived an asymptotic analytical expression for the distribution P (x) (see Sec. IV.B.3), which scales as x θ for x 1. Interestingly, if the noise distribution is heavy-tailed, with k < 2, its breadth k affects the value of θ > 0 in this mean-field model; θ turns into a non-universal exponent that depends on the loading and the amplitude of the noise, in line with atomistic simulations, and it supplements the other two exponents characterizing the depinning transition. While these predictions accurately match the value measured in lattice-based EPM relying on the genuine elastic propagator in dimension d = 4, some discrepancy was noticed in 2D and to a lesser extent 3D. This points to the progressive failure of the criterion of no-correlation as dimensionality is lowered. These results can be confronted 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).
More recently, Aguirre and Jagla (2018) underlined the need to improve on the foregoing approximation of mechanical noise as a random sum of single events. In reality, the noise results from avalanches displaying a linear spatial structure. A schematic argument taking into account the statistics of avalanche sizes S [namely, P (S) ∼ S −τ , see Sec. VI] then suggests to use the cumulative noise distribution W (δσ) ∼ |δσ| −1−k with k = 3 − τ (due to multiple avalanches), instead of the instantaneous single-event distribution w 1 (δζ). If the density of stability P (x) ∼ x θ coincides with the probability of presence of random walker subjected to the noise δσ near the aborbing boundary x = 0, they further speculate that θ = k − 1, which is consistent with their data.

Structure of the elastic propagator and soft modes
Coming back to the validity of mean-field estimates, we note that the latter become accurate (even for d < 4) if the elastic propagator G ij is shuffled, that is, replaced by G τ (i)τ (j) , where τ is a random permutation of indices which changes at each time step . This shows that the temporal correlations in the mechanical noise signal are due to the spatial structure of G. Of particular importance in this structure, claim Talamali et al. (2012) and Tyukodi et al. (2016b), are the soft deformation modes of the propagator (recall that these are the uniform shear bands described in Sec. III.C which create no elastic stress in the material). To clarify this importance, the authors focused on the evolution of the cumulative plastic strain pl in extremal dynamics and recast the EPM equation of motion [Eq. (5)] into a depinning-like equation (also see Sec. IX.B), viz., where η = µτ is a viscosity, σ y is the local stress threshold, and P(x) = x if x > 0 and 0 otherwise. The deformation of a disordered solid in d dimensions is then regarded as the advance of an elastic hypersurface in a (d+1)-dimensional space, where the additional dimension is pl . Under the influence of the driving, the hypersurface moves forward along pl . In so doing, it gets deformed as a result of the disorder in the thresholds σ y seen by different sites. Although EPM and depinning models share a formally comparable framework, Tyukodi et al. (2016b) showed that their results will differ because of the existence of soft modes in the EPM kernel G, while nontrivial soft modes are prohibited by the definite positiveness of the depinning propagator. As time goes on, the width W ≡ pl − pl 2 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 pl , which destabilize regions of lower pl and therefore act as restoring forces to homogenize pl over the hypersurface. On the contrary, in EPM, W (the variance of pl ) 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.
D. 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 {δζ(t)}. In the SGR spirit, these random stress 'kicks' operate as an effective temperature x that can activate plastic events, in the same way as thermal fluctuations do. Accordingly, the diffusive term in Eq. (22) is replaced by an Arrhenius rate 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) ρ(E) ∝ exp(−E/E g ).
Here, E g is a material parameter that will be set to unity. The external driving facilitates hops from trap to trap (over the local energy barrier E) by elastically deforming each region at a ratel =γ, where l is the local strain. This lowers the barrier: E → E − 1 2 kl 2 . (The stiffness parameter k is such that kl is the local stress.) Finally, SGR assumes that the random 'kicks' due to mechanical noise activate hops at a rate ω(E, l) given by an Arrhenius law, viz., where ω 0 is the attempt frequency and x quantifies the intensity of the mechanical noise. After a hop, l is set back to zero and a new trap depth E is randomly picked from ρ. The low-γ rheology that emerges from this simple model is quite interesting. As the effective temperature x decreases, the system transits from a Newtonian regime Σ ∝γ, for x > 2, to a power-law regime Σ ∝γ x−1 for 1 < x < 2. A yield stress emerges for x < 1 and the stress follows the Herschel-Bulkley law Σ − Σ y ∝γ 1−x . Indeed, for x < 1, the ensemble average of the time spent in a trap, viz., 10 Sketches illustrating the difference between thermal fluctuations ξT and mechanical noise ξ pl ; the variable l represents the local strain configuration. ξT are thermal kicks within a fixed potential energy landscape (PEL), whereas σ(t) and its fluctuations ξ pl tilt the local PEL up and down From (Agoritsas et al., 2015).
diverges atγ = 0. The system ages and falls into deeper and deeper traps on average; it follows that there is no typical material time for the relaxation of the cumulated stress. Moreover, the wealth of timescales afforded by an Arrhenius law also leads to interesting linear viscoelastic properties, in accordance with experimental data on colloidal pastes and emulsions.

Mechanical noise v. thermal fluctuations
However, the grounds for using an effective activation temperature to describe the effect of mechanical noise have been contested in recent years (Agoritsas et al., 2015;Nicolas et al., 2014a). The bone of contention is that, contrary to thermal fluctuations, mechanical noise fluctuations persistently modify the energy landscape of the region, insofar as the plastic events that trigger them are mostly irreversible.
More precisely, the argument runs as follows. The motion of particles in a mesoscopic region, of unit volume, is controlled by their interaction energy V 0 , subject to some constraint enforced at the boundary of the region. This constraint tilts the potential is the time-dependent stress imposed by the outer medium and l(t) is the strain associated with the internal configuration (see Fig. 10). While the region responds elastically, l(t) is a sensible proxy for its configuration. Assuming overdamped dynamics (with viscosity η), one can then write Thermal fluctuations, denoted by ξ T here, differ from mechanical noise in that they have a short-lived effect: in the case of white noise. Exceptional sequences of fluctuations ξ T are then required to climb up and overcome energy barriers. In contrast, changes to σ(t), due to either the external driving or distant plastic events, are persistent (hence cumulative). Even if one subtracts a constant drift term from σ(t), as we did in Sec. IV.A, the effect of mechanical noise fluctuations ξ pl (t) ≡ t 0 δζ (t ) dt is cumulative, viz., where the autocorrelation function C(∆t) ≡ δζ(t)δζ(t + ∆t) was assumed to decay quickly to zero. It follows that, under the sole influence of ξ pl , the energy barrier V σ flattens out after a diffusive time T ∼ (max dV /dγ) 2 ≡ σ 2 y , hence much faster than the Arrhenius law [Eq. (6)] encountered in activated processes.
The diffusive growth of local stress fluctuations with time has been confirmed by molecular dynamics simulations of model glasses at least at very low shear rates, where Puosi et al. (2015) have reported that (ξ pl (t + ∆t) − ξ pl (t)) 2 ∝ ∆t. One should however mention that in this context the observation of stochastic resonance induced by mechanical noise in lattice-Boltmann simulations of emulsions is puzzling (Benzi et al., 2015). The question of the mechanical noise has also sparked intense debate on the experimental side. In granular matter, it is now clear that locally shearing a region of the sample can affect distant, unsheared regions: The applied shear facilitates the penetration of an intruder (Nichol et al., 2010) or the motion of a rodlike probe (Reddy et al., 2011), presumably by agitating the grains in the distant region, as if they were thermally agitated. An Eyring-like activation picture using the magnitude of force fluctuations as temperature may indeed account for the observed fluidization. But Bouzid et al. (2015b) have argued that other nonlocal models can replicate this observation as well. Studying a related effect, Pons et al. (2015) have shown that applying small oscillatory stress modulations to a granular packing subjected to a small loading can dramatically fluidize it: Steady flow is then observed even though the loading is below the yield stress. This effect presumably stems from the cumulative impact of the stress modulations; the secular enhancement of the fluidity in the proposed rationalization is at odds with the expectations for any activated process.

E. Connection with the diffusion of tracers
Rather than the local stresses, many experimental works have access to observables related to particle displacements, in particular dynamic light scattering or particle tracking techniques. It is thus interesting to be able to connect the TABLE II Synthetic view of the distinct types of fluctuations at play and the methods with which they can be handled.  Applicable 'mean field' treatment None known so far Reasoning on P (σ) taking into account w(δζ)

Fluctuation-dominated regime
Eq. on P (σ) with dif--fusive term due to δζ Uniform mean-field approx. may be valid

References and examples
Lin and Wyart (2016)  local stress dynamics to the diffusion of tracer particles. Single events as well as plastic avalanches are expected to contribute to the tracers' motion even far away from the plastic zone due to their long-range effects, as sketched in Fig. 11(a) (Lemaître and Caroli, 2009;Nichol et al., 2010). The displacement field induced by a single plastic event can be calculated using Eq. (16) and is displayed in Fig. 11(b). To mimic diffusion, Martens et al. (2011) introduced imaginary tracers that follow the displacement field generated by the ongoing plastic events and were able to rationalize the relation between the nonaffine part of the self-diffusion coefficient D and dynamical heterogeneities (characterized by the four-point stress susceptibility χ 4 ), as shown in Fig. 11(c). In particular, D was found to decrease asγ − 1 /2 at highγ, while at low shear rates saturation at a value depending on the linear system size L as L 3 /2 was reported, contrasting with the D ∼ L scaling measured in atomistic simulations (Maloney and Robbins, 2008). Nicolas et al. (2014b) later argued that including convection in the EPM favored linear structures in the flow along the velocity direction and altered the dependence of χ 4 on the system size, but no proper scaling of the data could be achieved. The dispute was very recently settled by Tyukodi et al. (2018): It turns out that the roughly linear scaling of D with L observed in atomistic simulations is also present in EPM, but when one looks at short strain intervals ∆γ. More precisely, Tyukodi et al. (2018)'s extremal model yields a D ∼ L 1.05 scaling, which is robust to variations of implementation, and can be explained by considering that in each window ∆γ a roughly linear cascade of plastic events is either present or not. At longer times, the diffusivity enters the regime previously identified (D ∼ L 1.6 ), which is very sensitive to the presence of soft modes in the elastic propagator (see Sec. III.C.3).
Quantities comparable to the self-intermediate scattering function in purely relaxing systems are also accessible, as discussed in Sec. VIII.A.

F. Continuous approaches based on plastic disorder 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, which intend to stay closer to the schematic dynamics in the Potential Energy Landscape (PEL) outlined in Sec. I.C. The PEL is reduced to a free energy functional akin to that of Eq. (12), except that the strain tensor is conveniently traded off for new strain variables: one volumetric strain e 1 = xx+ yy 2 and two 'shear' strains e 2 = xx− yy 2 and e 3 = xy , in 2D (5 in 3D, e 2 , . . . , e 6 ). The multiple valleys in the PEL between which the system jumps during plastic rearrangements are reflected by multiple equilibrium values for the shear strains in the free energy functional.
For instance, Marmottant and Graner (2013) focused on the shear strain e 3 , split it into a cumulated plastic strain p and a complementary elastic strain, and proposed a minimalistic mean-field model based on an effective energy U eff that depends periodically on p , viz.
where E is an elastic modulus, y is a yield strain and 0 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 p to jump between energy valleys implies that, at high driving rates, p will not be able to instantaneously jump between, say, p , has appeared. This is similar to having a finite latency time before 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 a free energy of the form where B is the bulk modulus and e 1 , e 2 , and e 3 are explicit functions of the displacement field u(r). Here, F is an arbitrarily chosen function that is invariant under rotations of the reference frame θ → θ + π /3 (because a 2D triangular lattice is assumed) and periodic in its arguments. Introducing F in the equation of motion where ρ is the density, η 0 is the viscosity and σ R 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 e 1 and the density ρ, and to capture the elastic effects of edge dislocations, if the material is crystalline (Onuki, 2003b).
To avoid keeping track of the displacement field u, one may handle the strain components e 1 , e 2 and e 3 as independent primary variables, writing for instance F [e 1 , e 2 , e 3 ] = dr B e 1 (r) 2 + µ e 2 (r) 2 + V (e 3 (r)), if the only deviation from linear elasticity is borne by e 3 and encoded in a 'plastic disorder potential' V . Close to the reference state, V will not deviate much from linear elasticity, viz., V (e 3 ) ≈ µe 2 3 , but more globally it should have a corrugated shape that allows the system to reach new equilibria after each ST. The compatibility of (e 1 , e 2 , e 3 ) as differentials of a displacement field should then be ensured by the Saint-Venant condition This constraint is implemented by means of a Lagrange multiplier in the total free energy F , viz., F → F + λS (Jagla, 2007). It couples the different strain components. In particular, in an incompressible linear elastic solid (B → ∞ and V (e 3 ) = µe 2 3 ), a plastic strain arising at r = 0 (for example, if the potential V is chosen different at this position) will eventually unfold into the 'quadrupolar' elastic field given by Eq. (18) (Cao et al., 2018;Kartha et al., 1995). But, contrary to binary EPM, this unfolding is not instantaneous. Instead, it is generally governed by overdamped dynamics,ė i (r) ∝ −δF/δe i (r) (i = 1, 2, 3). An additional difference with respect to binary EPM is that the potential V affects the tangential shear modulus µ 3 = 1 /2V (e 3 ) as e 3 varies and may therefore alter the destabilization dynamics. In this case, the system becomes elastically heterogeneous, which precludes the use of Green's functions to calculate stress redistribution. Jagla (2017a) examined the influence of different functional choices for V on the flow curve and critical exponents, and reported differences between smooth and cuspy potentials.
In this chapter, different levels of detail in the description of the elastic interactions have been considered. We will see that the specific form of these interactions may impact the low-shear-rate rheology (see Sec. VII) and the local stress fluctuations (discussed in the next chapter), while many flow properties at high shear rates do not require as exquisite a description.

V. STRAIN LOCALIZATION: FROM TRANSIENT HETEROGENEITIES TO PERMANENT SHEAR BANDS
In Sec. I, we brought to light latent similarities in the deformation of amorphous solids. These, however, should not mask the widely different macroscopic consequences of applying shear to these materials. The elastoplastic viewpoint helps to understand these differences in a common framework.

A. Two opposite standpoints
In the common sense, there is a chasm between (i) foams and other soft solids, that flow, and (ii) metallic or silicate glasses that break/fracture after a certain amount of deformation (see Fig. 12b(right)).
To start with the far end of the latter category, perfectly brittle materials will deform elastically and then break, without going through a stage of plastic deformation. In daily life, this situation is exemplified by the soda-lime glasses routinely used to make windowpanes, bottles, etc., and more generally silicate glasses. Nevertheless, at small scales plastic deformations, resulting in a denser material, were revealed in indentation experiments with a diamond tip (Yoshida et al., 2007) as well as experiments of uni-axial compression of micropillars of amorphous silica (Lacroix et al., 2012) (which overall behaves comparably to soda-lime glass (Perriot et al., 2011)) and simulations of extended shear (Rountree et al., 2009). However, in many situations, plasticity plays virtually no role, in particular when failure is initiated by a crack: No evidence of plasticity-related cavities was seen by Guin and Wiederhorn (2004) (also see references therein) and, with the help of simulations, Fett et al. (2008) claimed that the surface displacements experimentally observed at crack tips are compatible with theoretical predictions discarding plasticity. (It should however be mentioned that a minority of works support the existence of plasticity near the crack tip).
In metallic glasses, global failure is preceded by substantial plastic deformation. The latter is generally localized in thin shear bands, that appear as clear bands in post-mortem scanning electron micrographs. These bands are typically 10 to 50nm or even 100 nm-thin (Bokeloh et al., 2011;Schuh et al., 2007), i.e., much thinner than the adiabatic shear bands encountered in crystalline metals and alloys, which are about 10 − 100 µm-thick. Despite these plastic deformations, brittleness remains a major industrial issue for metallic glasses. Added to their cost and the difficulty of obtaining large samples, this drawback may outshine their advantageous mechanical properties, such as their high elastic limit . As a consequence, much effort has been devoted to improving their ductility.
By contrast, foams, emulsions and various other soft solids can undergo permanent shear flow without enduring irretrievable damage. This conspicuous discrepancy with hard molecular glasses can however be lessened by noticing that, even among soft solids, the flow sometimes localizes in shear bands (Bécu et al., 2006;Lauridsen et al., 2004). Still, the distinction between hard solids that deform and break and soft solids that deform and flow is overly caricatural. The case of gels, which consist of long entangled (and often cross-linked) chains, demonstrates that soft solids, too, may break upon deformation. But, then, what distinguishes a material that flows from one that fails? What determines whether the deformation will be macroscopically localized in shear bands or homogeneous (on the macroscopic scale)?

The shear-banding instability from the standpoint of rheology
To start with, let us consider the rheological perspective. Shear-banding in complex fluids is interpreted as the consequence of the presence of an instability in the constitutive curve, i.e., the flow curve Σ 0 = f (γ) that would be obtained if the flow were macroscopically homogeneous. Indeed, it is easy to show that homogeneous flow in decreasing portions of the constitutive curve is unstable to perturbations and gives in to co-existing bands. The actual flow curve displays a stress plateau Σ (γ) = cst forγ between two valuesγ l andγ h . Shear localization corresponds to the particular caseγ l 0, i.e., that of a non-flowing band. In other words, it will occur if the constitutive curve already starts decreasing atγ = 0.
Accordingly, the shear-banding criterion based on the slope of the constitutive curve Σ 0 (γ) can be studied at the mean-field level (see, for instance, Coussot and Ovarlez (2010)'s analysis of a simple model). Incidentally, note that this is somehow counterintuitive, given the manifest spatial heterogeneity associated with the phenomenon. Nevertheless, mean-field calculations obviously leave aside the spatial organization of the flow (its banded structure), which hinges on the shape of the elastic propagator in simulations: In EPM, with similar dynamical rules, the banded flow structure obtained with the long-ranged elastic propagator of Eq. (18) is not preserved if the propagator is substituted by a stress redistribution to the first-neighbors, even if the latter is anisotropic (Martens et al., 2012).
The simple criterion based on the steady-state constitutive curve needs to be somewhat adjusted for amorphous solids, which often exhibit aging effects. Then, the yield stress of the quiescent material may vary with the waiting time since preparation (Varnik et al., 2003). Consequently, even if the flow curve obtained by ramping downγ from a high value is strictly monotonic, shear-banding may arise in non presheared samples. This will happen if an initially undeformed band gradually solidifies and thus further resists deformation, while the rest of the material is sheared. The solid band is 'trapped' in its solid state because of the aging at play (Martin and Hu, 2012;Moorcroft et al., 2011).

The mechanics of bands in a solid
Turning to the viewpoint of solid-state mechanics, as emphasized in Sec. III.C, uniform strain bands inclined by ± π /4 with respect to the principal directions of the strain tensor are soft modes of the elastic propagator [Eq. (18)], which means that they do not generate elastic stresses in the system. Should there be a weak stripe in the material (in the sense of low elastic moduli or low yield thresholds), it will then be energetically beneficial to accommodate part of the macroscopic strain in it in the form of a slip line. Such an energy-based argument is especially relevant in a quasistatic protocol in which the system always reaches the local energy minimum between strain increments. If the stripe in which the strain localizes displays ideal plasticity, the macroscopic stress-strain curve Σ = f (γ) stops increasing due to the banding instability.
But this continuum-based approach ignores the granularity of the material at the scale of plastic rearrangements by postulating the spontaneous and synchronous creation of a strain band all at once. Contrasting with this postulate, some experimental evidence in colloidal glasses (Chikkadi et al., 2011) and granular matter (Amon et al., 2012;Le Bouil et al., 2014) indicates that shear bands actually consist of disconnected, non-simultaneous localized plastic rearrangements, as implemented in EPM. Therefore, only on average is a strain band uniform; its granularity (as a patchwork of localized plastic rearrangements) as well as the time fluctuations in its plastic activity have no reason to be overlooked. The sequential emergence of the band may explain its sensitivity to details in the implementation of the elastic propagator (Talamali et al., 2012).
Taking the granularity of the band into account, Dasgupta et al. (2013Dasgupta et al. ( , 2012 proposed to explain the existence and the direction of shear bands by an argument based on the minimization of the elastic energy of a collection of Eshelby inclusions in a uniform elastic medium over their possible configurations in space. The neglect of the elastic heterogeneity of glasses in the reasoning was justified by the authors by the specific consideration of carefully quenched (hence, more homogeneous) glasses. An additional concern could be raised as regards the use of a global one-step minimization, whereas plastic events occur sequentially and the elastic deformation field in the material evolves during the process. Nevertheless, in a similar endeavor, Karimi and Barrat (2018) rationalized the observed deviation between the direction of the macroscopic shear band in a deformed granular medium and that of the microscopic correlations between rearrangements. The authors contended that the former direction is the direction of maximal instability with respect to the Mohr-Coulomb failure criterion, rather than that of maximal increase of the shear stress.
More generally, the strain bands described in the context of solids probably differ from the long-lived or permanent shear bands observed experimentally in steadily sheared materials. The former might be more accurately referred to as transient 'slip lines' and some reports of 'shear bands' in atomistic simulations should probably rather be interpreted as slip lines, as already noted by Maloney and Lemaître (2006). However, it has been suggested that the transient banding instability can act as a precursor to the formation of a shear band (Fielding, 2014).
In fact, transient banding is a matter of interest per se, as it can be long-lived (Divoux et al., 2010). Moorcroft and Fielding (2013) proposed a way to rationalize its occurrence on the basis of a generic banding criterion involving the transient constitutive curves Σ 0 = f (γ, γ), where γ is the cumulative strain since shear startup, in a fictitious system constrained to deform homogeneously. The rheological criterion dΣ0 dγ < 0 is recovered at infinite times γ → ∞, while a purely elastic banding instability is predicted if A ∂Σ0 ∂γ +γ ∂ 2 Σ0 ∂γ 2 < 0, with a model-dependent prefactor A > 0, provided that the material is sheared much faster than it can relax (γ → ∞). In the light of this, the authors claim that there is a generic tendency to transient banding in materials that exhibit a stress overshoot in shear startup.
Coincidence between a stress overshoot and an emerging shear band has effectively been noticed in EPM and atomistic simulations (Lin et al., 2015;Ozawa et al., 2018;Popović et al., 2018). Moreover, EPM unveiled the major impact of the system's preparation on these features. Within a depinning-like model, Ozawa et al. (2018) proved that a broad initial distribution P 0 (x) of site stabilities (i.e., distances to yield, as defined in see Sec. IV.B.3), reflecting poor annealing, suppresses the overshoot. Contrariwise, the latter and the associated stress drop grow as P 0 (x) becomes more sharply peaked, as expected for a well aged glass. At some point, the stress-strain curve even becomes discontinuous, as the system undergoes a spinodal instability. A critical point separates the continuous and discontinuous regimes. The authors observed very good agreement between these mean-field results and atomistic simulations of ultrastable glasses (Ozawa et al., 2018). In a parallel paper, Popović et al. (2018) showed that these features survive in a spatially-resolved EPM. In addition, they related the presence of run-away avalanches to the curvature of the distribution of site stabilities P (x) at x = 0. The stress drops caused by these avalanches grow as failure is approached and can thus signal the imminence of failure. But there may exist an alternative mechanism for failure (without diverging avalanches), namely the nucleation of a shear band in a fortuitously weak region. We will discuss an earlier work by Vandembroucq and Roux (2011) in Sec. V.C, where we review the mechanisms promoting shear banding.
For the time being, let us enquire about the fate of the material after the overshoot. Ozawa et al. (2018) and Popović et al. (2018) argue that the transition from a continuous stress-strain curve to a discontinuous one, as the initial preparation P 0 (x) is varied, marks a transition from a ductile response to a brittle one. Accordingly, the discontinuity in the stress curve is interpreted as irreversible material failure, in the very spirit of fiber bundle models, where fibers break irreversibly (see Sec. IX.C.2). In this case, the finite stress signal displayed by their systems (EPM and atomistic simulations) after failure must be considered spurious: The model is no longer valid after failure.
On the other hand, if the material does preserve some cohesion after the stress drop accompanying the overshoot, one may wonder whether the transient band will convert into a steady-state band, under homogeneous loading. What is required for this purpose is a mechanism that explains how the transient 'slip lines', instead of being dispersed, concentrate in the same region of the shear-banded material as time goes on. The distinction between the situation at finite strains and in the steady state should perhaps be emphasized. The first-order yield transition in the statistics of low-energy barriers observed by Karmakar et al. (2010) at a finite strain γ c is not necessarily associated with a firstorder (banding) transition in the steady-state flow curve Σ (γ). Similarly, Jaiswal et al. (2016) numerically observed that, in a batch of finite-size samples subjected to a strain γ, about half of the samples will have irreversibly yielded when γ = γ c , while the other half come back to their initial configuration upon unloading; but it is not straightforward to conclude from this interesting observation that, if one stitched a 'yielding' sample together with a 'recovering' one, a shear band would localize in the 'yielding' part at longer times.

B. Spatial correlations in driven amorphous solids
EPM help bridge the time and length scale gap between transient slip lines and permanent shear-banding. At short to intermediate time scales and under slow enough driving, the organization of the flow is complex and exhibits strong intermittency and marked spatial correlations between rearrangements even in driven amorphous solids that are not susceptible to macroscopic shear localization.

Spatial correlations
The spatial extent of correlations in the flow can be quantified by cooperativity or correlation lengths ξ in bulk flows, brought within reach by the computational efficiency of EPM. The Kinetic Elastoplastic (KEP) Theory of Bocquet et al. (2009), an extension of the Hébraud-Lequeux model (see Sec. IV.B.2) that includes heterogeneities, predicts a decrease of ξ with the shear rate as in contrast with Lemaître and Caroli (2009)'s theoretical prediction ξ ∼γ −1 /2 in 2D, beyond which independent avalanches are supposedly triggered. Simulations of homogeneous shear flow in spatially resolved EPM have generally shown results departing from the ξ ∼γ −1 /4 scaling. Picard et al. (2005) reported a correlation length that scales withγ −1 /2 in 2D (see Fig.13a), on the basis of a study of the variations of the average stress drop δσ withγ for different linear system sizes L; indeed, the data can be collapsed onto a master curve by rescalingγ into L 2γ . Nicolas et al. (2014b) related this scaling to the average spacing between simultaneous plastic events, which scales asγ −1 /d in d dimensions, and several definitions of correlation lengths were shown to follow this dependence in EPM. The variable sign of the elastic propagator enters the reasoning, insofar as plastic events are able to screen each other, because the sign of their contributions may differ.
Nevertheless, theγ −1 /d scaling is not generic. In particular, the correlation length derived from the four-point stress correlator G 4 (r), exploited by Martens et al. (2011) (see Fig.13b), is more sensitive to the avalanche shape and was shown to depend on the chosen EPM dynamical rules. Along similar lines, Roy (2015) disks in 2D point to a sensitivity of the correlation length ξ D derived from finite-time particle diffusion to the damping scheme; more precisely, they measured ξ D ∼γ −1 /3 for mean-field drag and ξ D ∼γ −1 /2 if the drag force depends on the relative particle velocities.
Below the yield stress, Lin et al. (2015) claim that the system is critical, with system-spanning avalanches in the transient, which is supported by a study of the cutoffs in the avalanche size distributions in EPM simulations. This implies a diverging correlation length ξ = ∞ in the whole Σ < Σ y phase -not unlike what is seen in 2D dislocation systems, at all applied stresses (Ispánovity et al., 2014). However, the divergence ξ → ∞ observed in athermal EPM, e.g., in the quasistatic limitγ → 0, will be strongly cut off in systems at a finite temperature, where thermal noise stifles the correlations (Hentschel et al., 2010). More generally, one should say that EPM tend to overestimate the absolute magnitude and and extent of the correlations between plastic events, e.g., compared to particle-based simulations (Nicolas et al., 2014c). We surmise that the overestimation is due not only to the neglect of elastic heterogeneities, but also to the regular lattice generally used in EPM, which standardizes the interactions between blocks.

Cooperative effects under inhomogeneous driving
Correlations in the flow dynamics are found in macroscopically homogeneous flows, but their impact is most conspicuous when the loading or the flow is inhomogeneous over the correlation length scale. In these situations, marked cooperative effects are generated. This is the case in pressure-driven flows through a narrow channel, of transverse width w (w ≈ 100 µm for microchannels). In this geometry, the streamline-averaged shear stress Σ varies linearly across the channel, from zero in the center (in 2D, but also in 3D) to w 2 ∇p at the wall. Therefore, a 'plug' of advected, but unsheared material is expected in the central region where |Σ| < Σ y , for yield-stress fluids.
These expectations were reshaped following seminal experiments on concentrated oil-in-water emulsions by Goyon et al. (2008). Indeed, compared to the predictions from the bulk rheology, the observed 2D velocity profiles are more rounded and, overall, the flow is enhanced. Thus, there is no unique relation between the local strain rate and the local stress (Goyon et al., 2010): the rheology is nonlocal. Using a similar system, Jop et al. (2012) demonstrated the existence of finite strain rate fluctuations δγ(r) ≡ γ(r) 2 − γ(r) 2 in the plug, which reach their minimum at the channel center. Numerical simulations of athermal soft disks confirmed the impact of confinement on the rheology: In 2D periodic Poiseuille flows, which put side by side Poiseuille flows of alternate directions, the wall stress below which flow stops substantially increases with decreasing channel width w . All these phenomena clearly arise because of the interactions between streamlines subjected to different stresses, via the elastic fields generated by plastic events.
The remarkable effect of spatial correlations on inhomogeneous flows is often rationalized by means of a nonlocal term in the equation controlling the fluidity f . This variable, defined here as the inverse viscosityγ/σ, is thought of as a proxy for the rate of plastic events. Owing to the symmetry of the propagator, the leading-order correction to the local fluidity involves the Laplacian ∇ 2 f . The steady-state fluidity diffusion equation thus reads where ξ is a cooperative length and f b is the fluidity in a bulk flow subject to stress Σ. The KEP model of Bocquet et al. (2009) provides a formal justification for Eq. (28) to linear order in f , with ξ ∼ (Σ − Σ y ) −1 /2 , by accounting for the mechanical noise generated in the immediate vicinity of plastic events. In fact, using a constant value of ξ (for each material) in Eq. (28) already provides very good fits of the experimental curves, not only for concentrated emulsions (Goyon et al., 2008) and lattice-Boltzmann simulations thereof (Benzi et al., 2014), but also for polymer microgels (Carbopol) (Geraud et al., 2013). In the case of emulsions, ξ vanishes below the jamming point, and reaches up to 3 to 5 droplet diameters (20 − 30 µm) in the very dense limit (Goyon et al., 2008). Similarly, for Carbopol samples, ξ measures 2 to 5 structural sizes (sizes of optical heterogeneities) (Géraud et al., 2017).
However, the fitting is highly sensitive to the adjustment of the fluidity f w at the wall (Geraud et al., 2013), which limits the accuracy of the experimental measurement of ξ. This difficulty highlights the value of EPM for testing the validity of theoretical predictions. In EPM descriptions of channel flow, the driving term µγ is set to zero in Eq. (5); flow arises on account of the initially imposed transverse stress profile Σ (y). Channel walls are accounted for by a no-slip boundary condition. This adds a correction to Eq. (18) for the elastic propagator, which can be calculated via a method of images (Nicolas and Barrat, 2013a), and leads to a faster local relaxation for plastic events near walls [also see (Hassani et al., 2018)]. Combined with appropriate dynamical rules, the model semi-quantitatively reproduces the shear rate fluctuations δγ(r) observed by Jop et al. (2012) in the plug as well as the moderate deviations of the velocity profiles from the bulk predictions witnessed with smooth walls, provided that the EPM block size corresponds to around 2 droplet diameters (see Fig. 13c). The fluidity diffusion equation, Eq. (28) either with ξ = cst or ξ ∼γ −1 /4 , captures the EPM fluidity profiles reasonably well, albeit imperfectly. Gueudré et al. (2017) took a closer look at the decay ofγ(y) in a region (y > 0) subject to Σ < Σ y contiguous to a sheared band (y < 0). They found that EPM results obey a scaling relation involving a length scale ξ(Σ) ∼ (Σ − Σ y ) −ν but the scaling exponents differ from mean-field predictions and are also inconsistent with KEP-based Eq. (28). In particular, for Σ ≈ Σ y ,γ(y) is argued to decay algebraically with y > 0 instead of exponentially. Gueudré et al. (2017) further claim that pressure-driven flows display larger finite-size effects than simple shear flows. Indeed, the finite size L of a system shifts the critical stress for flow initiation by ∆Σ ∝ L −1 /ν in a homogeneous setup (so that ξ(Σ) = L at initiation). On the other hand, in a pressure-driven flow, the length scale entering the critical stress for flow cessation should not be the system size L, but the (much smaller) width of the sheared band near the wall.
The description of nonlocal effects by Eq. (28) has also been applied to granular matter, which generically display heterogeneous flow and shear bands (Kamrin and Koval, 2012). To do so, the fluidity was redefined asγ/µ, because the rheology of dry frictional grains is best expressed in terms of the inertial number I (a rescaled shear rate) and its dependence on the friction µ ≡ Σ/P (with P the pressure). The resulting model successfully captures cooperative effects and accounts for the global velocity profile observed in discrete element simulations of a simple shear flow with gravity, a gravity-driven flow in a channel (Kamrin and Koval, 2012) as well as the flow of a granular layer on an inclined plate, which is sensitive to the thickness of the layer (Kamrin and Koval, 2012). Nevertheless, the validity of the definition of a 'granular fluidity', which is not an intrinsic state variable (because of the denominator Σ or µ), has been questioned; employing another variable would also lead to an exponential decay of the flow away from an actively sheared zone (Bouzid et al., 2015a). Other suggestions for the fluidity variable f that should enter a diffusive equation include the ratio between the 'static' and the 'fluid' part of the stress tensor (Aranson and Tsimring, 2006) and the inertial number I, which Bouzid et al. (2015b) claim to best match their discrete-element simulations, in particular regarding the necessary continuity of f = I at the interface between differently-loaded regions.

Cooperative effects due to boundaries
Coming back to emulsions, Goyon et al. (2008)'s observations indicate that the flow deviates much more from the bulk predictions, with an enhanced fluidization, when smooth walls are replaced by rough walls. Further experimental studies on regularly patterned surfaces show that the wall fluidization enhancement varies nonmonotonically with the height of the (steplike) asperities, for asperities smaller than the droplet diameter, as does the wall slip velocity (Mansard et al., 2014). These strong deviations in the presence of rough walls exceed by far what is found in EPM. This points to another physical origin than the coupling to regions subject to higher shear stresses. Since wall slip was experimentally observed, it has been suggested that the 'collisions' of droplets against surface asperities, as they slide along the wall, are the missing source of plastic activity; adding sources of mechanical noise along the walls in EPM can indeed capture the experimental features (Nicolas and Barrat, 2013b). Derzsi et al. (2017) experimentally confirmed the presence of roughness-induced scrambles at the wall and, with the help of lattice-Boltzmann simulations, the ensuing increase in the rate of plastic rearrangements near rough walls.
C. Alleged causes of permanent shear localization or fracture Several EPM have been able to reproduce permanent strain localization (Bulatov and Argon, 1994b;Coussot and Ovarlez, 2010;Jagla, 2007;Li et al., 2013;Martens et al., 2012;Nicolas et al., 2014b;Vandembroucq and Roux, 2011). In these cases, after a transient, the plastic activity will typically concentrate in a narrow region of space, generally a band, that may slowly diffuse over time. Hints at the ingredients suspected of causing this phenomenon come from its observation in certain (but not all) EPM and for a certain range of parameters only. Suspicions particularly target the rules for yielding or for elastic recovery, as we will see. Of course, relating these somewhat abstract rules to microscopic physical properties may not be straightfoward. Therefore, the interpretation remains mostly qualitative, with very few detailed comparisons so far between microscopic calculations and EPM rules.
To start with, one notices that large applied stresses Σ Σ y are incompatible with localization. Indeed, the applied stress then exceeds the local yield stresses: Plastic events pervade the system, which globally flows in a viscous manner. In other words, large loadings fluidize the material, consistently with experimental observations (Divoux et al., 2012).
On the other hand, at lower stresses (hence, lower shear rates), plastic events are sparser and may hit the same regions over and over again, provided that the latter are strongly or durably weakened by these events. Meanwhile, in the rest of the material, the driving term in Eq. (5) is compensated by the nonlocal contributions due to a band of plastic events, i.e., a uniform relaxation (Martens et al., 2012). The general cause for localization thus evidenced is the insufficient healing of regions following rearrangements (Nicolas et al., 2014b). In the following, we look into the distinct possible origins of this weakening.

Long breakdowns (rearrangements), slow recovery
Coussot and Ovarlez (2010) rationalized shear-banding in jammed systems by considering the formation and breakage of particle clusters. Locally, these events delimit periods of solid and liquid behavior, with elastic stresses σ el = µγt and σ el = 0, respectively, that comes on top of a constant viscous stress ηγ. Here, µ is the shear modulus andγt is the local strain. On the basis of a mean-field argument, they showed that if the liquid-like phase lasts longer than η /µ, then the flow curve becomes nonmonotonic, which is the hallmark of shear-banding. The idea was elaborated by Martens et al. (2012), who used a spatially resolved EPM of the Picard type with a variable rearrangement ('healing') time τ res as a parameter, with the notation of Eq. (21). Their findings confirmed the formation of shear bands in space for large τ res , associated with the emergence of nonmonotonicity in the macroscopic flow curve. The banded flow shares many properties with systems at a first-order transition in which different phases coexist; the shear rate is well defined (independent of the driving) inside the band and there is an interface with the nonflowing phase. This spatial organization in the form of a band is intrinsically related to the band being a soft mode of the propagator (see Sec. III.C); this would not be possible without its long range and its anisotropy.
Attractive interactions in adhesive colloidal systems (Irani et al., 2014) and directional bonds in molecular systems are tentative candidates for possible microscopic origins of long rearrangements, i.e., long time delays before the destabilized region reaches another stable configuration. Similarly, deactivating potential forces for a finite 'pinning delay' after yielding enhances strain localization (Papanikolaou, 2016).

Influence of initial stability (aging) and shear rejuvenation (softening):
In (thermal) amorphous solids, with age comes strength and above all stability (see Sec. I.B.2). Yielding will then be more abrupt. Indeed, letting a system age in the absence of strain favors strain localization, or even fracture, as indicated by experimental (Rogers et al., 2008) and numerical (Shi and Falk, 2005) data. The EPM proposed by Vandembroucq and Roux (2011)  a plastic event was shifted by an amount δ with respect to the initial P (σ y ), to mimic the lower structural temperature of the pristine material. For large enough negative δ (strain weakening), the first regions to yield are rejuvenated to a state with lower threshold, so that the system gets trapped in a banded structure. The bands thus created are localized and pinned in space if the elementary slip distance is small; otherwise, larger slip events are created, enhancing nonlocal effects and making bands less stable and more diffusive. Nicolas et al. (2014b) introduced a healing process in this picture, by allowing the blocks that have just become elastic again to age and gradually recover higher energy barriers, viz., where k is the rate of recovery at which the energy barrier rises from its post-yielding value E min y to the asymptotic value E ∞ y . For low enough recovery rates k, shear localization was observed. However, the localized behavior tends to fade away whenγ reaches very small values. This may be paralleled with the recovery of a homogeneous flow in the dense colloidal suspension studied by Chikkadi et al. (2011) for shear rates below a certain value, which allow the strained system to structurally relax before further deformation.
Along similar lines, Li et al. (2013) implemented a process of free volume creation and annihilation in a finiteelement-based EPM designed to describe the deformation of the metallic glass Vitreloy 1. In their model, free volume is created by the dilation accompanying a shear transformation (ST) and is annihilated gradually in strictly local diffusional events. The activation of STs, in turn, is facilitated by a local excess of free volume. Simulations relying on a kinetic Monte Carlo scheme for the dynamics showed that at low temperatures the deformation localizes in bands and that the variations of free volume are critical for this localization. A parallel can obviously be drawn between the creation of free volume during STs and the lowering of yield stresses in other EPM. There is perhaps an even stronger connection with the plasticity-induced enhancement of the local effective temperatures x and χ in variants of the Soft Glassy Rheology model (Fielding et al., 2009) and the Shear Transformation Zone theory (Manning et al., 2007), respectively. In these models, the effective temperature evolves locally in response to three processes: (i) rises due to plastic rearrangements, (ii) relaxation to a steady-state value, and (iii) diffusion in the shear gradient direction. For a range of parameters, the homogeneous temperature profile is unstable and a high-temperature shear band emerges in the midst of a low-temperature unsheared background.
Another approach to account for the competition between local relaxation and driving-induced plastic events was proposed by Jagla (2007). In his continuous model (see Sec. IV.F), the system relaxes via a slow drift of the local energy landscape seen by a given site towards lower energies. Sites whose evolution towards potential minima is not interrupted by plastic deformations benefit from this local 'structural relaxation'. Their elastic energy decreases and the local yield stress increases; while their plastically active counterparts have no time to undergo structural relaxation, and their yield stress remains consequently low. Again, this leads to a nonmonotonic flow curve in a mean-field analysis, and to strain localization at lowγ.
To what extent precisely these strain localization mechanisms are connected with the weakening-induced runaway (system-spanning) events observed in Fisher et al. (1997)

Shear bands like it hot
A temperature rise ∆T has been experimentally evidenced during the operation of shear bands in metallic glasses (Lewandowski and Greer, 2006;Zhang et al., 2007). The dominant view is that it is however not the initial cause of the shear banding observed at low strain rates, as ∆T is small in this case. Still, local heating may result in the recrystallisation of the material, with associated changes in its mechanical properties (presumably more brittle behavior). Such effects are obviously not included in EPM, and probably better described at the level of macroscopic equations as a thermomechanical instability. The discussion above is therefore only relevant for the initiation of the instability and for systems in which thermal effects are weak.
A related mechanism leading to a nonmonotonic flow curve, first identified in MD simulations (Nicolas et al., 2016) and then also seen in finite-element-based EPM (Karimi and Barrat, 2016), is at play when one enters the underdamped regime. At a given strain rateγ, long-lived inertial vibrations can then be sustained, because of which the yield threshold may be exceeded earlier than if mechanical equilibrium had been maintained the whole time. In other words, the energy dissipated in the underdamped flow remains long enough in the relevant degrees of freedom to activate plastic relaxation. In MD simulations this facilitation was shown to be quantitatively described by a heating effect, whereby a more strongly damped system is heated to a strain-rate-dependent temperature T (γ). For strongly underdamped systems, this leads to a nonmonotonic constitutive curve and to the formation of shear bands, if the system is large enough 3 .

VI. CRITICAL BEHAVIOR AND AVALANCHES AT THE YIELDING TRANSITION
Amorphous solids retain complex solid-like properties under continuous flow, but the onset of flow is of particular interest from a physical viewpoint owing to the critical behavior that may come along with this transition. Far from being a weakness, the simplified description provided by EPM (which were originally phenomenological models) represents an asset for the study of these critical properties. In this section we review the thriving literature about the statistics of avalanches close to the yielding transition.

A. Short introduction to out-of-equilibrium transitions
Statistical physics is largely concerned with phase transitions, whereby some properties of a system abruptly change upon the small variation of a control parameter. The paradigmatic example of an equilibrium phase transition is the Ising model, which consists of spins positioned on a lattice and interacting with their first neighbors. This model describes the ferromagnetic to paramagnetic transition of a magnet as the temperature T rises above a critical temperature T c ; the transition is marked by the presence of correlated domains of all sizes and the vanishing of the magnetization m (the "order parameter") as Quite interestingly, the critical exponents, β and its kin, are shared by many other, a priori unrelated systems: The latter are said to belong to the same universality class as the Ising model. These ideas extend beyond equilibrium, but fewer methods are available to deal with the dynamical phase transitions encountered out of equilibrium. In this respect, it is worth noting that the Herschel-Bulkley constitutive law can be recast into an expression analogous to Eq. (29), viz., This yielding transition is receiving more and more attention as an example of transition in a driven system. Still, the existence of a critical behavior asγ → 0 is not unanimously accepted: While some authors reasonably argue that thermal fluctuations will wash out criticality at any T > 0 (Hentschel et al., 2010), others claim that the material's state should become independent ofγ once the driving gets slower than any internal relaxation rate (Langer, 2015), perhaps overlooking the possibility that the latter (for instance, mediated by the propagation of shear waves) may become unboundedly long for increasingly large systems under slow driving. Besides, among the defenders of criticality, there have been lively discussions as to whether it belongs to the same universality class as the depinning transition of driven elastic lines (Sec. IX.B).

Avalanches in sandpile models
As models featuring threshold dynamics and a toppling rule, EPM are also connected to the somewhat simpler sandpile models, introduced in Sec. III.A. Let us clarify some concepts using the latter class of systems.
Simulations of 2D sandpile models display avalanches of grains of duration T (number of iterations to reach stability) and size S (total number of transferred grains). Because of the rules governing grain transfer, these avalanches are compact structures, unlike those observed in EPM. In a nutshell, at vanishing deposition rate, the cumulative distributions of S and T exhibit power-law scalings, with a cut-off at large scales due to the finite size of the system, viz., where τ > 0 and τ > 0 are critical exponents, f and g are fast decaying functions, and the positive exponents d f and z are called the fractal dimension of the avalanches and the dynamical exponent, respectively. This means that small avalanches are more frequent than larger ones, but in such a fashion that no typical or characteristic size can be established, which has been called self-organized criticality. Let us note that the extremal dynamics used to trigger avalanches can be substituted by a very slow (quasistatic) uniform loading of the columns of sand if some randomness is introduced in the stability thresholds. In this sense, self-organized criticality in the sandpile model simply exposes the criticality associated with the dynamical phase transition undergone by the loaded system. Different regimes of avalanches can be seen when the deposition rate is varied or, somewhat equivalently, when one inspects them at different frequencies ω. At large frequencies ω, independent, non-overlapping avalanches are observed. As ω decreases, the avalanches start interacting. In this regime, their overlaps cut off the correlation lengths of single avalanches, but due to mass conservation during grain transfer, the scale-free behavior is preserved. On long time scales, i.e., for low ω, the observed features are typical of discharge events, whereby the whole sandpile becomes unstable after having been loaded (Hwa and Kardar, 1992).

Stress drops and avalanches in EPM
Similarly to the instabilities in sandpile models, the plastic events occurring in EPM can trigger avalanches of successive ruptures. To facilitate the comparison with experiments or atomistic simulations, these avalanches are usually quantified by looking at the time series of the macroscopic stress σ(t) and, more specifically, at the stress drops ∆σ associated with plastic relaxation (in this chapter, we use a lowercase symbol σ for the stress to underscore that it is an intensive variable). Close to criticality, the duration T of these drops and their extensive size S ≡ ∆σL d in a system of volume L d in d dimensions most often display statistics formally similar to Eq. (31), viz.
where the upper cut-offs S cut and T cut entering the scaling functions f and g will typically depend on system size, e.g. S cut ∝ L d f . In the following, we will pay particular attention to the possible impact of the peculiarities of the quadrupolar stress redistribution in EPM, notably its fluctuating sign, on the avalanche statistics.

B. Avalanches in mean-field models
Shortly after the emergence of the first EPM, mean-field approximations were exploited to determine the statistics of avalanches. Most of these approaches assume a uniform redistribution of the stress released by plastic events, as exposed in Sec. IV.A. An exponent τ = 3 /2 is then consistently found in the avalanche size scaling of Eq. (32). For instance, Sornette (1992) proposed to map the Burridge-Knopoff model for earthquakes, introduced in Sec. III.A, onto a fiber bundle which carries a load equally shared among unbroken fibers (see Sec. IX.C). At criticality the extremal load needed to make the weakest surviving fiber break fluctuates; more precisely, it performs an unbiased random walk. An avalanche of failures lasts as long as this extremal load remains below the initial load, so its size is given by the walker's survival time close to an absorbing boundary, whence an exponent τ = 3 /2. If deformation starts farther from the critical point, a larger exponent is found, τ = 5 /2. A posterior, but widely celebrated model for heterogeneous faults in earthquakes was proposed by Fisher et al. (1997), and later applied to the deformation of crystals by Dahmen et al. (2009) and more recently to the deformation of granular matter (Dahmen et al., 2011) and amorphous solids (Antonaglia et al., 2014). Here, the problem is directly mapped onto an elastic line depinning problem (see Sec. IX.B). Once again, above an upper critical dimension that decreases with the interaction range, the model yields the mean-field exponent τ = 3 /2. But if a post-yield weakening mechanism is introduced or if stress pulses due to inertial effects are present, the power-law regime only holds for small avalanches, while larger ones trigger runaway events that span the whole system and result in a bump at a characteristic size in the avalanche statistics.
Much more recently, there have been endeavors to extend mean-field approaches in order to account for the nonpositiveness of the redistributed stress, which undermines the mean-field reasoning. For instance, the Hébraud-Lequeux model introduced in Sec. IV.B.2 features an additional diffusive term acting on local stresses. Jagla (2015) studied avalanches in a discrete variant of this model and reported on subtleties that are absent from depinning problems. Indeed, if avalanches are artificially triggered by picking a random block and destabilizing it, the problem can yet again be mapped onto a survival problem for an unbiased random walk, similarly to the fiber bundle, and the mean-field exponent τ 3 /2 is obtained. (In passing, with a random-kick protocol of the sort, Lin et al. (2014a) had arrived at a similar result, in two EPM variants.) But now consider the physically more relevant protocol of quasistatic loading, where stresses are uniformly increased until a block is destabilized. The foregoing result still holds in the depinning case, because the distribution of local stresses σ is fairly homogeneous close to the yield point σ y , viz. p(σ ≈ σ y ) cst, and thus insensitive to the stress shift induced by the uniform loading. By contrast, stress fluctuations in disordered solids deplete local stresses close to σ y , so much so that p(σ − y ) = 0 and p(σ) varies substantially close to σ y . Accordingly, significantly smaller exponents τ 1.1 − 1.2 are both predicted and observed numerically in that case (Jagla, 2015). Furthermore, the power law is cut off at a value S cut that depends on the distance to criticality and on the system size. An extension of these results to heavy-tailed distributions of stress fluctuations  is still pending.

Experiments
Various experimental settings have been designed to characterize avalanche statistics in deformed amorphous solids in the last decade, even though experiments are still trailing behind the theoretical predictions and numerical computations in this area. Let us mention examples of such works. Lauridsen et al. (2002) sheared a foam in a Couette cell and investigated its plastic behavior. The distributions P (S) of normalized stress drops S (plotted in Fig. 15a) were shown to follow a power law at three different shear-rates, with an apparent exponent τ 0.8 in Eq. (32). This value was reported to be consistent with the bubble model of Durian (1997), but contrasts with other theoretical predictions, as we will see. It should however be noted that the power law was fitted over barely a decade in S.
At the other end of the softness spectrum, the compression of millimetric metallic glass rods was studied by Sun et al. (2010) and the stress drops were analyzed. Again, P (S) follows a power law regime over one decade of experimental measurements, but this time with exponents in the range τ ∈ [1.37, 1.49], as can be seen in Fig. 15b. Among several works that came in the wake of this seminal paper, Antonaglia et al. (2014)'s compression experiments of microsamples were argued to be compatible with the mean-field prediction P (S) ∼ S −3 /2 . Following the same approach, Tong et al. (2016) reported exponents in the range τ ∈ [1.26, 1.6] for four different samples of a Cu 50 Zr 45 Ti 5 alloy.
A granular packing subject to the simultaneous application of pressure and shear was also shown to display stress drops with power-law statistics by Denisov et al. (2016). The power-law exponents, which seem to lie in a relatively broad range in Fig. 15c, were not fitted, but, upon rescaling, were reported to be in good agreement with the meanfield value τ = 3 /2. It remains uncertain to what extent the value reported in this work and in the other ones may have been influenced by the large body of literature claiming that the deformation of (a large variety of) amorphous materials belongs to a unique universality class, the one describing the depinning of an elastic line (Dahmen et al., 2009;Dahmen, 2017). Also note that the two decades of raw (non-cumulative) stress drops over which Denisov et al. (2016) could collect data, at least for the smaller strain rates, make granular packings particularly promising experimental test systems. Barés et al. (2017)'s recent study of a sheared bidisperse mixture of photoelastic particles confirms these promises. From the gradient of the image intensity, they quantified the local pressure acting on each grain, hence the energy stored in it, and tracked the fluctuations of the global energy. This allowed them to define avalanches as spontaneous energy drops, with a dissipated power E related to granular rearrangements. The researchers reported power-law distributions of avalanches, P (E) ∼ E −τ with τ = 1.24 ± 0.11, in a range dependent on the threshold used to filter the signal and spanning over three decades in E.

Atomistic simulations
In parallel to experiments, stress drops have been analyzed in atomistic simulations of the deformation of glassy materials. In a 2D packing of soft spheres, Maloney and Lemaître (2004) measured power-law distributed energy drops with an exponent τ = 0.5 − 0.7 comparable to that obtained in Durian (1997)'s foam experiments. On the contrary, exponential distributions of stress drops and energy drops were then reported in athermal systems of particles interacting with three distinct potentials in 2D (Maloney and Lemaître, 2006), but also with a more realistic potential for a metallic glass in 3D (Bailey et al., 2007). All these studies were however limited to fairly small system sizes. Using larger systems in 2D and 3D, Salerno and Robbins (2013) found power-law distributed energy drops and stress drops, with distinct values for the exponent τ in the overdamped regime and the underdamped one, and in 2D and 3D. In the overdamped case, the value is identical in 2D and 3D, τ = 1.3 ± 0.1. We also mention that, opposing the rather consensual view of scale-free avalanches and non-trivial spatiotemporal correlations, Dubey et al. (2016) suggested that the characteristics of the stick-slip behavior stemmed from trivial finite-size effects.

D. Avalanche statistics in EPM
The large amount of statistics afforded by EPM can enlighten the debate about the criticality of the yielding transition and the existence (or not) of a unique class of universality by overcoming the uncertainty and limitations of some experimental measurements. In the last years, EPM have tended to challenge the strict amalgamation of the yielding transition with the depinning one.  (Talamali et al., 2011). (b) In strain-controlled simulations with the 'image sum' implementation of the elastic propagator kernel (see Sec. III.C.2), with fits to Eq. (33). Notice that the driving springs are much stiffer than in panel (a). From (Budrikis and Zapperi, 2013). (c) Under extremal dynamics, for systems of different sizes L. The exponent reported for the unscaled curves (inset) is τ 1.2, while the rescaled curve shown in the main plot was fitted with τ 1.36. From (Lin et al., 2014b).

Avalanche sizes in the quasistatic limit
Avalanches are most easily defined in the limit of quasistatic driving, in which the external load is kept fixed during avalanches (Sec. II.C). Applying extremal dynamics to a 2D EPM, Talamali et al. (2011) defined an avalanche size S as the number of algorithmic steps ∆t during which the external stress Σ remains lower than Σ start − k∆t, as though the system were driven by a slowly moving spring of stiffness k. Quite interestingly, driving the system by pulling on it with a moving spring is equivalent to a strain-controlled protocol in the limit k → ∞, while a stress-controlled protocol is recovered in the opposite limit k → 0 (Popović et al., 2018). Talamali et al. (2011)'s numerical simulations displayed a scale-free distribution P (S) ∝ S −τ with τ = 1.25 ± 0.05 cut by a Gaussian tail (Fig.16a). It was made explicit that this result is at odds with the mean-field exponent τ = 3 /2. On the other hand, the measured value is similar to that measured by Durin and Zapperi (2000) (τ 1.27) for one class of Barkhausen avalanches, due to the motion of ferromagnetic domain walls under an applied magnetic field, and to that predicted for this effect using a model of elastic line depinning with anisotropic (dipolar, but positive) interactions (Zapperi et al., 1998). Of at least equal relevance is the similarity with the avalanche size exponent τ 1.25 found when simulating differential equations (Bonamy et al., 2008) or cellular automata (Laurson et al., 2010) to describe the interfacial growth of a crack in a heterogeneous medium. Indeed, the alignment of plastic events along the Eshelby 'easy' axes was seen as an effective dimensional reduction, leading to avalanches belonging to a quasi 1D problem with positive interactions decaying as r −2 , similarly to the interfacial crack growth model of Bonamy et al. (2008).
A couple of years later, Budrikis and Zapperi (2013) exploited a closely related EPM, with randomly distributed stress thresholds, to investigate the effect of two distinct implementations of periodicity for the long-range elastic propagator G defined in Sec. III.C.2. A first series of simulations focuses on the nonstationary plastic activity below the macroscopic yield stress Σ y , by adiabatically increasing the applied stress Σ. Overloaded blocks yield simultaneously; their strain is increased by dγ = 0.1 and a new local yield stress is drawn. For Σ Σ y , avalanche distributions P (S) are found to decay as exponentials (or compressed exponentials); but for stresses closer to Σ y , a short powerlaw regime appears. The distributions can be fitted by Le Doussal and Wiese (2012)'s first-order correction to the mean-field prediction for depinning (see Eq. 35), but with τ 1.35 instead. The tails of P (S) collapse upon rescaling with a cutoff depending on the distance to the critical point S cut ∝ (Σ y − Σ) −1 /σ with 1 /σ 2.3. In a second series of simulations, apparently inspired by Talamali et al. (2011), the system was pulled by a spring of stiffness k ∈ [0.1, 1], moving adiabatically, hence an external stress Σ = k(γ tot − γ), where γ tot is the position of the spring and γ is the plastic strain. In this case, the statistics can be improved because the system eventually reaches a critical steady state. The avalanches size distributions show a larger power-law regime and a little 'bump', and the authors fit them with the empirical shape where all coefficients are free. Two very close, but not strictly identical, exponents τ 1.35 (again) were measured for different implementations of the propagator; the precise values of τ can be found in Fig. 16b and Table III. These values somewhat differ from Talamali et al. (2011)'s measurement, presumably because larger spring constants k were used. Still, they definitely deviate from the mean-field value, too. The authors also considered the avalanche durations T , measured in algorithmic steps and fitted them to the power law P (T ) ∝ T −τ , with τ 1.5. Joining these researchers, Sandfeld et al. (2015) tested the robustness of these avalanche statistics to (i) variations of the boundaries, (ii) different computations of stress redistribution and (iii) finite-size effects. To do so, they used an eigenstrain-based finite element method with different types of meshgrids, and found that these variations have no influence on the critical exponents. Lin et al. (2014a) implemented two slightly different automata based on the Hébraud-Lequeux model but embedded in finite dimensions. In stress-controlled simulations at Σ ∼ Σ y , in which sites are randomly 'kicked' to trigger an avalanche, they found an exponent τ larger than that of the quasistatic simulations described above: They measured τ ≈ 1.42 in both model variants. This may be a consequence of the random-kick protocol [see Sec. VI.B and (Jagla, 2015)]. Yet, later on, Lin et al. (2014b) reported τ 1.36 in 2D and τ 1.43 in 3D, for the same protocol. Besides, power-law distributions were reported for the avalanche durations, with exponents τ 1.6 in 2D and 1.9 in 3D. In parallel, extremal dynamics were implemented and yielded smaller exponents for the same models, τ 1.2 in 2D and τ 1.3 in 3D, closer to previous quasistatic approaches, even though not devoid of finite-size effects.  Mean field (Fisher et al., 1997) [depinning] 3 /2 2 -0 2 (Jagla, 2015) [Hébraud-Lequeux like] 1.1 − 1.2 --1 -Legend -: † Estimated from the avalanches shape. * Obtained using the τ exponents from the random-kick protocol.

Connection with other critical exponents
A discussion on the density of zones close to yielding and its connection with the critical exponents was opened up by Lin et al. (2014a). Denoting x ≡ σ y − σ the distance to threshold of local stresses, a stark contrast was emphasized between depinning-like models, with only positive stress increments and p(x) ∼ x 0 for small x, and EPM, where a pseudo-gap emerges at small x, viz., p(x) ∼ x θ with θ > 0. In Lin et al. (2014a)'s stress-controlled simulations with randomly 'kicked' sites, identical values of θ were obtained in two variants of the model embedded in 2D (θ 0.6) and 3D (θ 0.4), whereas the stress-strain curves differed (see Table III for the slightly smaller values of θ measured using extremal dynamics).
Shortly afterwards, Lin et al. (2014b) proposed to link p(x) with P (S), in a scaling description of the yielding transition. Their scaling argument can be summarized as follows. Starting from Eq. (32), one obtains ∆σ ∝ L d f (2−τ )−d . Now, in a stationary situation, on average this stress drop must balance the stress increase that is applied  2D and (b) in 3D systems, whose size L is varied. From (Lin et al., 2014b). (c) In a 3D system, as the shear rate is varied. The inset shows the rescaled distribution of avalanche sizes. From .
to trigger an avalanche. Among the L d sites, the one with the smallest x, x min , will start the avalanche, so ∆σ ∝ x min .
Identifying the two expressions of ∆σ leads to which is supported by their EPM simulations (notably with the random-kick protocol). In this regard, the discrepancy was once again underscored between the depinning transition (with fractal dimensions d f ≥ d typically, due to the compactness of the avalanches, and a velocity-force exponent β 1) and the yielding transition (with typically d f < d and a rheological exponent β 1 in Eq. (30)). Generalized scaling relations encompassing both transitions were put forward (see Lin et al. (2014b)-Supporting Information).

At finite strain rates
Seeking to narrow the gap between experiments and EPM,  analyzed the EPM stress signal with methods mimicking the experimental ones and studied the effect of varying the applied shear rateγ. At very lowγ, avalanches are power-law distributed with an exponent τ 1.28 in 2D and τ 1.25 in 3D, cut off by finite size effects with d f = 0.90 and 1.3, respectively. These results coincide very well with MD simulations in the quasistatic limit and support the nascent convergence towards an avalanche size exponent τ 1.25 in 2D or 3D EPM, deviating from the (depinning) mean-field value 3 /2. Much more tentatively, there might be a downward trend of τ with increasing dimensions, which would be compatible with Jagla (2015)'s suggestion τ 1.1−1.2 above the upper critical dimension.
Interestingly,  observe a systematic crossover towards higher values of τ when the shear rate is increased, so that τ reaches τ 1.5 at intermediateγ, before entering the high-γ regime of pure viscous flow. At the same time, the external driving starts to dominate over the signed stress fluctuations originating from mechanical noise; this nudges the system into a depinning-like scenario, with an exponent θ in p(x) ∼ x θ decreasing towards zero asγ reaches finite values both in 2D and 3D. Note that the same effect occurs in a stress-controlled system, as soon as the imposed stress gets perceptibly lower than the yield stress (Budrikis et al., 2017). Similarly to pulling the system with a stiff spring (large k), increasing the shear rate generates simultaneous uncorrelated plastic activity in the system, which leads to larger τ , closer to 1.5. Overall, applying a finite shear rate does not destroy the criticality of avalanche statistics; but it affects the critical exponents and eventually produces more trivial effective statistics.

Insensitivity to EPM simplifications and settings
At present, technical difficulties still hamper a clear discrimination between theoretical predictions through experiments. The simplifications used in the models thus need to be carefully examined. Budrikis et al. (2017) investigated the effect of the scalar approximation of the stress (see Sec. II.C.2) by comparing the results of a scalar model to those of a finite-element-based fully tensorial model, under different deformation protocols (uniaxial tension, biaxial deformation, pure shear, simple shear) and in both 2D and 3D. Irrespective of the dimension, and (most of the) loading and boundary conditions, a universal scaling function is observed for the avalanche distribution, shown in Fig. 18 and coinciding with Le Doussal and Wiese (2012)'s proposal Heterogeneous deformations, such as bending and indentation, were also considered and yielded similar values for τ . Nevertheless, the cutoff value is different from the homogeneous cases. This is not unexpected: An independent length scale enters the problem and the yield stress Σ y used to measure exponent σ is not universal. Also, while the observed τ value was nearly identical in the different (homogeneous) loading cases when treated separately, some range of variation was observed for exponent σ ∈ [1.53, 2.05]. Finally, the average avalanche size was related to its duration T via S ∼ T γ with γ = 1.8 ± 0.1.
A possible explanation for the insensitivity of avalanche statistics to the aforementioned aspects may lie with the quasi-1D geometry of the avalanches, resulting from the quadrupolar propagator. Most cooperative phenomena thus appear to be controlled by the stress component along one direction, and a scalar description may be sufficient in this respect. (Scalar models do indeed reproduce the same power-law exponent and evidence a fractal dimension d f ≈ 1 in 2D and 3D, as shown in Fig.18c).

Effects of inertia
Without the assumption of instantaneous stress redistribution, stress waves are expected to propagate throughout the system (see Sec. III.D and Fig. 9), in a ballistic way or a diffusive one depending on the damping. This is not described by the traditional elastic propagator G of Eq. (17), but finite-element based EPM have recently made it possible to account for inertial effects (Karimi and Barrat, 2016). Karimi et al. (2017) exploited this type of model to study Salerno and Robbins (2013)'s claim, based on extensive atomistic simulations in the quasistatic regime, that inertial effects drive the system into a new (underdamped) class of universality. At odds with this claim, but consistently with results of sandpile models (Khfifi and Loulidi, 2008;Prado and Olami, 1992) and seismic fault models (Carlson and Langer, 1989), they found that inertial effects destroy the universal, scale-free avalanche statistics. A characteristic hump (or secondary peak) of large events emerges in the avalanche size distribution P (S), similarly to Fisher et al. (1997)'s findings. In Karimi et al. (2017)'s work, both the relative weight and the scaling with the system size of this peak are controlled by the damping coefficient Γ, a dimensionless parameter that quantifies the relative impact of dissipation. The effective fractal dimension d f (Γ) of avalanches, albeit dependent on the damping, satisfies a scaling relation with the exponent θ (Γ) defined by p(x) ∼ x θ (Γ) .

FIG. 19
Avalanche shapes in experiments and EPM. (a) Experimental avalanche shapes, for avalanches of fixed duration (top) and fixed sizes (bottom) in a bulk metallic glasses (BMG) and a granular system. From (Denisov et al., 2017). (b) Avalanche shapes at different fixed durations in a strain-controlled EPM simulation at fixedγ. From . (c) Avalanche shapes at fixed sizes. From (Budrikis et al., 2017).
These results are compatible with Papanikolaou (2016)'s phenomenological description of inertial effects, which are accounted for by a temporary vanishing of elasticity after local plastic events (plastic delay): Simulations of the model showed the appearance of a hump of large events in P (S), an increase of the exponent τ , as well as the emergence of dynamical oscillations, accompanied with strain localization.

Avalanche shapes
In addition to their duration and size, further insight has been gained into the avalanche dynamics by considering their average temporal signal, i.e., the 'shape' of the bursts. This observable can be determined experimentally with higher quality (Antonaglia et al., 2014;Denisov et al., 2017). Avalanche shapes have thus been estimated for various systems displaying crackling noise; examples include earthquakes (Mehta et al., 2006), plastically deforming crystals (Laurson et al., 2013), and the Barkhausen noise (Mehta et al., 2002;Papanikolaou et al., 2011).
In the latter example, the magnetization of a film changes mostly changes via the motion of domain walls 4 ; its rate of change is recorded as a time series V (t). When the film thickness, which controls the long-range dipolar interactions, is such that mean field is valid, the average shape V (t|T ) of avalanches of duration T is well described in the scaling limit by an inverted parabola (Papanikolaou et al., 2011), viz., Since oftentimes mean field does not hold, a generalized functional form was proposed by Laurson et al. (2013): Here, the shape factor γ is also the exponent that controls the scaling between size and duration (S ∼ T γ ), since S(T ) is nothing but the integral T 0 V (t|T )dt. γ and the parameter a control the asymmetry (a > 0 refers to positive skewness); the mean-field formula is recovered for γ = 2 and a = 0. As the interaction range increases from local to infinite, the university-class parameters evolve from γ 1.56, a 0.081 to γ 2.0, a 0.01. Dobrinevski et al. (2015) provided an analytic formalization for this expression as a one-loop correction around the upper critical dimension; these authors also computed the shape of avalanches of fixed size S. The need for this generalization beyond mean field was confirmed by Durin et al. (2016). In the deformation of amorphous solids, the inverted-parabola shape predicted by mean field was shown to provide a satisfactory description of some experiments, e.g. in metallic glasses (Antonaglia et al., 2014) and granular matter (Denisov et al., 2017), even though the agreement is not perfect (see Fig.19a). On the contrary, Barés et al. (2017)'s granular experiments point to a clear asymmetry in the shape of long avalanches.
On the EPM side,  studied the effect of finite shear ratesγ on the avalanche shape. By sorting the avalanches according to their duration T , at fixedγ, they found that short avalanches are noticeably more asymmetric and display faster velocities at earlier times (positive skewness, see Fig.19b). For larger T , it is argued that avalanches most likely result from the superposition of uncorrelated activity, which leads to more mean-field like results. This would explain the gradually more symmetric shapes observed for increasing T (see the evolution of the asymmetry parameter in the inset of Fig. 19b). In the quasistatic limit, asymmetric stress-drop shapes are then expected. Indeed, at lowγ, fits with Eq. (37) give a non-mean-field value γ 1.58 in both 2D and 3D. This feature gradually disappears at largerγ. Budrikis et al. (2017) extended the study to different loading conditions, in 2D and 3D, and measured values for γ in the range [1.74, 1.87]. Contrary to 's findings with a scalar EPM, they saw clearly asymmetric avalanches with positive skewness only in the bending and indentation protocols, and not (visibly, at least) in the tension and shear simulations. In addition, the shapes obtained by sorting the avalanches according to their sizes (see Fig.19c) collapsed well with the scaling form proposed by Dobrinevski et al. (2015), with a shape exponent γ 1.8 (note the difference with respect to the mean-field value γ = 2).

VII. STEADY-STATE BULK RHEOLOGY
In this Chapter we redirect the focus to materials that flow rather than fail. This is the relevant framework for foams, dense emulsions, colloidal suspensions, and various other soft glassy materials exhibiting a yield stress. Experimentally, in the absence of shear-banding, their flow curve (which characterizes the steady-state macroscopic rheology) is generally well described by the Herschel-Bulkley law [Eq. (1)] Σ = Σ y + Aγ n ; the exponent n typically lies in the range 0.2-0.8, often around 0.5, perhaps closer to 0.3 for foams (Bécu et al., 2006;Bonn et al., 2017;Jop et al., 2012). Such a nontrivial dependence on the shear rateγ proves that the rheology of these materials cannot be understood as a mere sequence ofγ-independent elastic loading phases interspersed withγ-independent plastic events (Puglisi and Truskinovsky, 2005). Instead, the violation of the quasistatic conditions of Eq. (8) at finite shear rates implies that the specific dynamical rules implemented in EPM will affect the rheology, at odds with the situation observed in the quasistatic limit. In particular, we will see that the local yielding and healing dynamics (notably via rules R2 and R4 in Sec. II.A) play a crucial role in determining the flow properties at finiteγ. More generally, different flow regimes will be delineated, depending on the material time scales at play.
A. Activation-based (glassy) rheology v. dissipation-based (jammed) rheology The rheology of glasses was long thought to be tightly connected with that of jammed systems such as foams (Liu and Nagel, 1998). But the contrast between the role played by thermal activation in the former and the importance of dissipative processes for the latter has pointed to prominent differences. The particle-based simulations of Ikeda et al. (2012Ikeda et al. ( , 2013, in particular, contributed to disentangling glassy (thermal) and jammed (athermal) rheologies. They identified the time scales of the relevant processes at play, namely Brownian motion and dissipation, in addition to the driving, and were able to separate the thermal rheology associated with the former from the dissipation-based one. As the shear rateγ is increased, the driving starts to interrupt the thermally triggered plastic relaxation and dissipation starts to dominate the rheology. The idea of a transition from a thermal to an athermal regime was bolstered by experiments on microgel colloidal suspensions, which are impacted by thermal fluctuations close to the transition to rigidity (glass transition), but obey jamming-like scalings further from the transition (Basu et al., 2014). Already in the first EPM, the competition between the driving and the realization of plastic events was emphasized. Since then, it has been implemented in different ways.
For flows dominated by thermal relaxation, it makes sense to consider EPM of the type of Bulatov and Argon (1994a)'s and Homer and Schuh (2009)'s, as well as Sollich et al. (1997)'s SGR model (presented in Sec. IV.D.1), in which plastic events are activated at a rate given by an Arrhenius law [Eq. (6)]. For instance, in SGR, where the Arrhenius law is controlled by a fixed effective temperature x, asγ increases, blocks can accumulate more elastic strain before a plastic event is activated. The macroscopic stress thus increases withγ. It does so in a non-universal way: The flow curve at lowγ follows a Herschel-Bulkley law [Eq. (1)] with exponent 1 − x for x < 1.
Most other EPM dedicated to the study of steady-state rheology consider systems close to the athermal regime, in particular foams and dense emulsions of large droplets, which in practice undergo negligible thermal fluctuations. This will be the focus of the rest of this chapter. These athermal systems will depart from the quasistatic conditions of Eq. (8) because the driving time scale γ y /γ competes with either the time scale τ pl of individual rearrangements, or, at lower shear rates, the duration τ av of avalanches. A rough upper bound for the latter is given by the propagation time (if applicable) across the system of linear size L, viz., τ av ∼ L z , where the exponent z depends on the damping regime. Concretely, at vanishing shear rate, the flow consists of well separated avalanches, some of which span the whole system. These avalanches are gradually perturbed and cut off asγ is increased, while higher shear rates add further more local corrections. The external driving thus hampers the relaxation of the system at increasingly local scales as it gets faster. Therefore, at least two scaling regimes could be seen asγ is varied (Bonn et al., 2017), and, indeed, two regimes were observed by  in their EPM simulations, as shown in Fig. 21(a). We will discuss these regimes separately.
B. Athermal rheology in the limit of low shear rates At vanishingly low shear ratesγ the nonlinear response of athermal materials is anchored in the critical dynamics discussed in Sec. VI. Accordingly, there have been many recent endeavors to deduce the Herschel-Bulkley exponent n (which controls the variations of Σ whenγ → 0) from other critical properties.
To start with, the phenomenological mean-field model proposed by Hébraud and Lequeux (1998) and discussed in Sec. IV.B.2 leads to n = 1 /2, a typical experimental value. Recall that, in this model, the local stresses in the system drift and diffuse due to endogenous Gaussian white noise, and yield at a finite rate τ −1 c above a local threshold σ y . The growth of the macroscopic stress withγ mirrors the associated decrease of the relative yielding rate (γτ c ) −1 , which makes the boundary σ = σ y more permeable. The result n = 1 /2 is robust to several variations of the yielding rules (Olivier, 2010;Olivier and Renardy, 2011), notably the inclusion of a distribution of yield stresses (Agoritsas et al., 2015), and tensorial generalizations of the model for multidimensional flows (Olivier and Renardy, 2013). On the other hand, it varies if a shear-rate dependence is introduced in the elastic modulus or the local restructuring time (Fig. 21(b)  ).
Although this model can fit several aspects of athermal rheology, the assumption of Gaussian mechanical noise fluctuations has been debated. Indeed, the distribution of stress releases due to a single plastic event is heavy-tailed, w 1 (δζ) ∼ |δζ| −1−k with k = 1 for the elastic propagator (see Sec. IV.C.3). Accordingly, more cautious approximations have been propounded. Assuming that the stress received by a block is a sum of random stress increments drawn from w 1 in a variant of the Hébraud-Lequeux model, Lin and Wyart (2018) showed that the Herschel-Bulkley exponent n, while equal to 1 /2 for any k > 2, rises to 1 /k for mechanical noises characterized by k ∈ [1, 2]. For the physically relevant value k = 1, an exponent n = 1 is thus found, up to logarithmic corrections. To derive these results, the authors perturbed the density of stabilities P (x) around its critical state at the yield stress (x = |σ| − σ y is the local distance to yield). However, the mean-field values thus obtained are larger than those measured in the corresponding spatially-resolved EPM in 2D (n 0.66) and 3D (n 0.72), although the discrepancy seems to shrink as the dimensionality is increased.
Besides, Lin and Wyart (2018) argued that the flow exponent β = 1 /n should obey the hyperscaling relation where the exponent z relating the duration and the size of avalanches (T ∼ L z ) is claimed to be close to 1 and the fractal dimension d f of avalanches can be expressed with Eq. (34) [also see (Lin et al., 2014b)]. Equation (38) yields flow exponents β somewhat larger than typical experimental values and than those actually measured in EPM with instantaneous elastic propagation [Eq. (17)]; it was claimed that the difference originates from a better account of the finite speed of shear waves. So far, we have seen that, within mean-field models, the dynamics of shear wave propagation and the (heavy-tailed) statistics of mechanical noise fluctuations may affect the low-shear-rate rheology, and that the finite dimension of space introduces deviations from mean-field predictions due to correlations in the noise. Another ingredient of the models is worth studying: the way blocks soften when they are destabilized -or, in other words, the plastic disorder potential V (e 3 ) of Eq. (27), where e 3 ∼ γ is the shear strain. We recall that binary EPM rely on linear elasticity [V (e 3 ) ∝ e 2 3 ] within the elastic regime, which is equivalent to V being a concatenation of parabolas joined by cusps. For the problem of elastic line depinning (Sec. IX.B), it is known that such a cuspy potential will not give the same results as a smooth potential in mean field, because destabilization is very slow atop a smooth hill (which has a flat crest), whereas it is instantaneous at a cusp. The discrepancy vanishes in finite dimensions, because sites are destabilized by abrupt 'kicks' anyway. Aguirre and Jagla (2018) suggest that the situation differs widely for the yielding transition. Within a 2D continuous approach based on a plastic disorder potential [see Sec. IV.F], they separate the mean background of the mechanical noise from its zero-average fluctuations δσ(t) and arrive at an equation of motion which schematically reads where η is a viscosity and k > 0 can be interpreted as the constant of a spring connecting the current strain e 3 to a driven 'wall' at w(t). This equation can be compared with Eq. (20), but it should be noted that the fluctuations δσ are here cumulative, i.e., integrated since t = 0. As mentioned in Sec. IV.C.3, Aguirre and Jagla (2018) argue that spatial correlations (avalanches) affect the distribution of these fluctuations, W (δσ) ∼ δσ −1−k ; the value k = 1 expected for uncorrelated (and unbounded) stress releases should thus be substituted by k = 1.5 if one considers objects as extended as avalanches. Once again, we should note that this result is heavily impacted by the non-positiveness of the elastic propagator, which undermines purely mean-field arguments. Equation (39) describes the motion of a particle pulled by a spring on a corrugated potential; it is a stochastic Prandtl-Tomlinson equation. This model was worked out by Jagla (2018), who obtained the following scaling relation where α = 1 for the cuspy parabolic potentials V and α = 2 if V is smooth. We also recall the speculation k = θ + 1 from Sec. IV.C.3. In the case of parabolic potentials, the proposed scaling relation is nicely obeyed by the values of k, θ, and β 1.51 measured in their simulations, as well as those found by  in 2D [β 1.54(2)] and 3D. Contrary to θ or k, the flow exponent β is thus found to explicitly depend on α, i.e., the shape of the potential (Jagla, 2017a). The scaling relation (40) seems to involve fewer parameters of the problem than Eq. (38); one should nonetheless bear in mind that in depinning problems the relevant effective potential V entering mean-field reasoning nontrivially depends on different properties of the system. Besides the choice of a specific potential shape, an alternative way to model the different destabilization speeds is to introduce stress-dependent transition rates (Jagla, 2017b). This, too, yields diverse exponents β. In this regard, note that a dependence of the flow exponent on the specific form of the viscous dissipation was reported in particle-based simulations .
Despite the very promising recent works in this direction, no firm theoretical consensus has been reached yet regarding the flow exponent β and how universal it is. This exponent clearly has a value distinct from that encountered in elastic depinning problem and, as we already discussed at length in Sec. VI, the mechanical noise fluctuations induced by the alternate sign of the elastic propagator most probably play a prominent role in these deviations close to criticality. As one departs from the low-shear-rate limit, among other corrections, the mechanical noise properties are expected to vary. Uncorrelated events will start to occur simultaneously, presumably leading to a more Gaussian noise distribution (see Table II), and scalings different from those obtained atγ → 0 may be expected.

C. Athermal rheology at finite shear rates
Rheology is concerned not only with the onset of flow atγ → 0, but with the whole range ofγ > 0. The regime of finite driving rates was already targeted by early EPM, including that of Picard et al. (2005) (see Sec. IV.A). In this athermal model, elastoplastic blocks can only yield when their stress exceeds a uniform local threshold σ y . Yielding is then a stochastic local process and so is the subsequent elastic recovery; these processes have fixed rates τ −1 and τ −1 res , respectively [see Eq. (21)]. Figure 20(a) shows the resulting monotonic flow curve for τ = τ res = 1. Its shape broadly matches that of many experimental flow curves, but more quantitatively the simulated rheology does not follow a Herschel-Bulkley law: It crosses over to a Newtonian regime at stresses Σ only slightly above the yield stress Σ y . This is due to the postulated elastic stress accumulation above the threshold σ y for a fixed duration τ on average.
These seemingly oversimplified yielding and healing rules have been refined since then. To make the picture more realistic, Nicolas et al. (2014a) opted for an instantaneous triggering of plastic events at σ y ; they also introduced a yield stress distribution. In their model, the event lasts for a fixed local strain 'duration' γ c . Therefore, the local dissipation process can be disrupted by the external driving, which contributes to the local strain. This mimicks the fact that, contrary to any problem of depinning on a fixed substrate, a deforming region in a solid will not wade through the same potential landscape at different driving rates. The ensuing flow curves are more compatible with experimental ones and are well described by a Herschel-Bulkley law, as shown in Fig. 20(c).  noticed that this model actually exhibits a transition from a low-shear-rate regime, characterized by a Herschel-Bulkley exponent n 0.65, to a regime with an exponent n 0.5 over the range Σ ∈ [1.06, 1.6] in units of Σ y , beyond which further corrections to scaling set in.
It is worth noting that many experimental soft systems exhibit a qualitative change of their flow behavior when the adhesion properties of their constituents are modified, with e. g. higher propensity to shear-banding when an emulsion is loaded with bentonite, creating attracting links between droplets (Ragouilliaux et al., 2007) [also see (Bécu et al., 2006)]. This discovery prompted the idea that there exist different classes of jammed systems depending on microscopic interactions. Coussot and Ovarlez (2010) suggested that adhesion results in longer local restructuring events. Martens et al. (2012)'s EPM-based studies confirmed that long plastic events lead to a nonmonotonic constitutive curve and the formation of permanent shear bands in the unstable parts of the flow curve, as discussed in detail in Sec. V.C and shown in Fig. 20(b).
Thus, EPM and experiments highlight the sensitivity of the finite-shear-rate rheology to the specific microscopic interactions between particles or dynamical rules at play. On the other hand, quite interestingly, this finite-shear-rate regime appears to be amenable to mean-field approaches oblivious to correlations in the flow. While the latter are pivotal for avalanches, stronger driving decorrelates plastic events. The mechanical noise felt in a given region, then, results from the superposition of a large number of events and its distribution acquires a Gaussian shape . Concomitantly, as we noted in Sec. IV.C.1, mechanical noise fluctuations play a less important role. As a result, one should not be particularly surprised that Martens et al. (2012) succeeded in reproducing the finite-shear-rate rheology of Picard et al. (2005)'s model with a mean-field approach discarding fluctuations. In the same vein, we remark that overall flow curves are only moderately altered by finite-size effects, whether it be in particle-based  or EPM computations [see Fig. 20(a)].

D. Strain-driven vs. stress driven protocols
Most EPM works consider strain-controlled protocols (defined in Sec. II.C). Some of the counterexamples are given by Lin et al. (2014b) [see Fig. 21(c)] and the recent work by Jagla (2017a). Another example of stress-imposed modeling is the numerical work in Liu (2016)'s PhD thesis. In a section dedicated to the transient dynamics prior to fluidization, a stress-controlled EPM is introduced. To this end, the internal stress resulting from plastic events is separated from the externally applied stress field, which can be chosen arbitrarily.
In this type of protocols, depending on the initial condition, two types of stationary solutions are obtained, namely, steady flow and a dynamically frozen state. Under athermal conditions the system may always reach a configuration with large local yield stresses, in which the dynamics gets stuck, even if the applied stress Σ is larger than the dynamical yield stress Σ y . The smaller Σ and the smaller the system size, the more likely becomes the visiting of such an absorbing state. But if a flowing stationary state is reached for a given time and granted that the mechanical  . The data suggest two different scaling regimes: Close to criticality the Herschel-Bulkley exponent is n = 0.65, whereas at highγ n tends towards 1/2. (b) Flow curves for the same EPM at relatively highγ with a shear-rate-dependent local shear modulus G0(γ) ∼γ ψ 1 and plastic event 'duration' γc(γ) ∼γ −ψ 2 . These dependences introduce corrections to the Herschel-Bulkley exponent, n = (1 + ψ1 − ψ2)/2 (pay attention to the horizontal axis on the plot). From . (c) Flow-curves obtained from stress-imposed EPM simulations in 2D. The dashed line is a fit to a Herschel-Bulkley law with n ≈ 0.66. From (Lin et al., 2014b).
properties do not show history dependence (Narayanan et al., 2017), strain-controlled and stress-controlled protocols yield identical flow curves .

VIII. RELAXATION, AGING AND CREEP PHENOMENA
So far EPM have mostly been exploited to investigate the macroscopic flow behavior and flow profiles (Sec. V), characterize stationary flow (Sec. VII), or study fluctuations and correlations in the steady flow close to criticality, where one finds scale-free avalanches (Sec. VI). Still, some works, however few, are concerned with relaxation, aging, and creep phenomena. This section is dedicated to both the dynamics in the temperature assisted relaxation (aging) of disordered systems and to the transient dynamics under loading (creep), prior to yielding or complete arrest. The latter phenomenon can be either an athermal process, provided that the stress load is above, but close to, the yielding point, or thermally assisted creep, in response to a load below the dynamical yield stress.

A. Relaxation and aging
A striking feature in the theory of viscous (glassforming) liquids is their response to an external perturbation, close to the glass transition: They do not exhibit an exponential structural relaxation, with a simple time scale, but a stretched exponential relaxation. More specifically this means that the temporal behavior of the response function R(t) (e.g., the response in stress Σ(t) to the application of a strain step at time t = 0) can often be described by the so called Kohlrausch-Williams-Watt (KWW) function In this expression, b typically takes a value between 0 and 1, which stretches the exponential relaxation. This was ascribed to the formation of dynamical heterogeneities close to the glass transition, thus producing separately relaxing domains and leading to a broad distribution of relaxation times (Macedo and Napolitano, 1967), hence a stretched exponential relaxation (Bouchaud, 2008;Campbell et al., 1988). With this picture in mind, it came as a surprise that a series of dynamical light-scattering measurements on colloidal gels showed the opposite behavior, namely, a compressed exponential structural relaxation, characterized by an exponent b > 1 (Cipelletti et al., 2000(Cipelletti et al., , 2003Ramos and Cipelletti, 2001). More recent experiments using X-ray photon correlation spectroscopy have found that this feature is not specific to gels (Orsi et al., 2012), but also arises  (Ferrero et al., 2014). (c) Relaxation of thin Zr67Ni33 metallic glass ribbons with time, measured by the the decay of the X-ray photon intensity autocorrelation g2 at T = 373K, for different waiting times tw and wavectors Q (shown in the inset). The characteristic relaxation time τ (Q, tw) was determined by fitting R = g2 − 1 to the KWW form of Eq. (41), which yielded a shape parameter b 1.8 ± 0.08. From (Ruta et al., 2013). in supercooled liquids (Caronna et al., 2008), colloidal suspensions (Angelini et al., 2013) and even in hard amorphous materials like metallic glasses (Ruta et al., 2013(Ruta et al., , 2012. Although this anomalous relaxation was observed ubiquitously in experimental systems, it took more than a decade to reproduce dynamics with compressed exponential decay in molecular-scale simulations, until Bouzid et al. (2017) and Chaudhuri and Berthier (2017) eventually reported such dynamics in microscopic models for gels. The main obstacle had been to probe the right parameter range, notably with respect to temperature and also length scales.
From the outset, Cipelletti et al. (2000) suggested that the faster than exponential relaxation stems from the elastic deformation fields generated by local relaxation events. Shortly afterwards Bouchaud and Pitard (2001) put forward a mean-field model based on the assumption of elasticity to explain this anomalous relaxation. In case this explanation is correct, EPM should be the ideal tool to test it (Ferrero et al., 2014). In a quiescent system, the driving term vanishes in Eq. (5), which turns intoσ where pl j (t) denotes the local plastic deformation at site j and the other notations were defined below Eq. (5). As before, this equation describes the response of the surrounding medium to local relaxation events. Here, only thermally activated processes are relevant, and their modeling is inspired by the the trap model of Denny et al. (2003) and Sollich et al. (1997)'s Soft Glassy Rheology model [SGR, see Eq. (24)], with an Arrhenius-like yielding rate for sites below the threshold, viz., while sites with |σ| > σ y yield instantaneously. In Eq. (42), the signs correspond to the direction of the yielding event, σ y is a local yield stress, κ is a dimensional prefactor, and T the ambient temperature. Such models confirm the dependence of the shape parameter b of structural relaxation on the dimensionality of the system, which Bouchaud and Pitard (2001)'s mean-field arguments predict to be b = 3 2 in 3D and b = 2 in 2D. Moreover, in EPM insight into the microscopic dynamics can be gained by following the motion of tracers advected by the elastic displacement field, as explained in Sec. IV.E. This led Ferrero et al. (2014) to distinguish three dynamical regimes in 2D, namely (I) ballistic, (II) subdiffusive and (III) diffusive. In the ballistic regime (see Fig. 22), compressed relaxation was found, with a shape parameter b ≈ 2. The subdiffusive regime was ascribed to correlations in the relaxation dynamics, a feature that has not been reported in experiments. This disagreement can either be due to oversimplifications of the model or to the fact that experiments are usually performed in 3D, and not 2D. Preliminary EPM studies in 3D observed ballistic motion at short times, with a compressed exponent b = 3 2 , followed by a diffusive regime 5 .
There remain many other open questions that could be addressed by EPM. For instance the q-dependence of the experimental intermediate scattering functions S(q, t) (Cipelletti et al., 2003) cannot be captured in EPM at present, but could be included by implementing hybrid models that consider smaller-scale dynamics as well. Besides, the self and the intermediate part of S(q, t) cannot be distinguished in EPM yet, because the tracers do not interact, but the two may differ in reality. Other questions include the 3D dynamics and the possibility of intermittency in time as well as spatial correlations of the localized relaxation events.

B. Creep
Another field that has stimulated much experimental work in the last years (Bonn et al., 2017) but few rationalization attempts at the mesoscale is creep. The definition of creep is somewhat ambiguous. In some contexts it may refer to stationary motion at a vanishingly small velocity, in particular the creep dynamics of a driven elastic manifold over a disordered landscape at finite temperature , but also the flow of a granular medium subjected to a constant stress Σ Σ y supplemented with an additional small cyclic stress modulation (Pons et al., 2016). But here we will restrict our attention to the traditional definition in material science, namely, the slowdown of deformation prior to failure, fluidization or complete arrest, under load Σ. This load is usually comparable to, or smaller, than the material yield stress Σ y and creep can in principle be both of thermal and athermal nature.
For Σ > Σ y the usual response of most dense soft glassy materials can be separated into three regimes (in polymeric systems, five regimes are listed by Medvedev and Caruthers (2015)). Primary creep corresponds to a first slowdown of the dynamics, with a gradual decrease of the (initially high) strain rateγ. The deformation rate is roughly constant in the secondary creep regime but abruptly shoots up in the tertiary regime, which ultimately culminates in macroscopic failure or fluidization. The measured macroscopic quantities are usually the time-dependentγ(t) and the fluidization or failure time τ f (Divoux et al., 2011a;Skrzeszewska et al., 2010).
Creep is observed in many experimental systems, from crystalline and amorphous solids to soft materials. In the former materials, a power-law slowing down of the deformation rate with an exponent close to or slightly less than 2 /3 is often reported (Miguel et al., 2002), viz., This law is commonly called Andrade creep and hints at a possible universality of the dynamics. However, experiments and simulations on creep in amorphous systems have found a variety of power-law exponents for the decay ofγ(t) in primary creep, ranging between −1 /3 (Bauer et al., 2006) and −1.0 (the latter value corresponding to logarithmic creep γ(t) ∼ ln(t)), with a multitude of values in-between (Ballesta and Petekidis, 2016;Chaudhuri and Horbach, 2013;Divoux et al., 2011b;Landrum et al., 2016;Leocmach et al., 2014;Sentjabrskaja et al., 2015). Bonn et al. (2017) extensively reviewed the literature on the topic. Scaling results for the fluidization (or failure) time τ f also vary and basically fall in two classes. Among other works, Divoux et al. (2011b) found a power-law scaling of τ f , defined as the time to reach a homogeneous stationary flow, viz., τ f ∼ (Σ y − Σ) −β , where β varies between 4 and 6. On the other hand, other works defined τ f as the duration of the rapid increase ofγ(t) at the end of secondary creep and reported an inverse exponential dependence τ f ∼ exp Σ0 Σ , where a characteristic stress scale Σ 0 has been introduced (Gibaud et al., 2010;Gopalakrishnan and Zukoski, 2007;Lindström et al., 2012).
Thus, rather than a universal behavior, experiments suggest a multitude of dependencies, notably on the preparation protocol prior to the application of the step stress (quench or pre-shear), on temperature, age and also on the dominant physical process at play during creep. In some systems the initial creep regime appears to be completely reversible and one expects the creep to be a result of visco-elasticity. Accordingly, Jaishankar and McKinley (2013) were able to reproduce the experimental power-law creep in Acacia gum solutions using a modified Maxwell model featuring fractional time derivatives. On the other hand, on the basis of molecular dynamics simulations, Shrivastav et al. (2016) claim that the power-law creep in a variety of glassy systems can be related to a percolation dynamics of mobile regions, thus plasticity, which would render EPM particularly suitable to tackle the open questions in the field. Among the 'hot topic' highlighted by Bonn et al. (2017), the detection of precursors that may point to incipient failure stands as the Atlantis in many disciplines from material science to engineering and geology. Using a lattice-based EPM, Bouttes and Vandembroucq (2013) made a first endeavor to address thermal creep and showed its strong dependence on initial conditions and the impact of aging on the creep behavior. In the model, each site is assigned an energy barrier E 0 (renewed after every plastic event) in the stress-free configuration, with a uniform distribution of E 0 . The elastic stress redistributed by plastic events via the usual elastic propagator [Eq. (18)] biases this potential. The plastic activation probabilities are analogous to Eq. (42), with an Arrhenius-like law, and are resolved with a kinetic Monte-Carlo algorithm. The resulting creep dynamicsγ(t), studied in pure shear, depend on the applied stress Σ and temperature T and all display an apparent exponent suggestive of logarithmic creep (see Fig. 23a). Besides, the fluidization time τ f is found to decrease with increasing Σ and T . Merabia and Detcheverry (2016) explored the transient thermal creep that occurs upon application of a stress step, prior to steady flow, at relatively high temperatures. Within an EPM, they also resorted to a kinetic Monte-Carlo scheme and Arrhenius-type plastic activation rates, but they used a non-uniform distribution of intrinsic trap depths ρ(E 0 ). With an exponential distribution ρ(E 0 ) ∼ exp[−αE 0 ] (leaving aside a lower cutoff), the model is formally similar to the SGR model (see Sec. IV.D.1), but here the temperature parameter is interpreted as the room temperature, instead of an effective noise temperature, and samples are assumed to be thermally equilibrated before stress is applied (αk B T > 1). Contrary to Bouttes and Vandembroucq (2013), the simulated creep does not always slow down logarithmically. Instead, a power-law decay γ(t) ∼ t α−1 is observed, for 1 < α < 2, in agreement with a mean-field analysis; it tends to logarithmic creep as α → 1. Merabia and Detcheverry (2016) also considered a Gaussian distribution ρ(E 0 ). In that case, the steady-state flow curve grows logarithmically, Σ ∼ ln(γ). Regarding the creep regime before steady state, the cumulative strain contains a term that grows linearly in time and the fluidization time τ f follows the inverse exponential dependence on Σ [i.e., τ f ∼ exp Σ0 Σ , see Fig. 23c)] found in experiments on carbopol black gels by Gibaud et al. (2010). The latter result is robust to variations of the Gaussian half-peak width.
The authors also tried different stress propagators of short range character, besides the quadrupolar (Eshelby-like) one. It turns out that their mean-field predictions agree best with the simulations with a short-range propagator and an exponential distribution of energy barriers, whereas there is a systematic offset in the creep exponent with respect to the more realistic quadrupolar propagator. This is somewhat counter-intuitive because increasing the interaction range usually leads to a more mean-field-like behavior.
An alternative mean-field approach is based on the Hébraud-Lequeux model (see Sec. IV.B.2). Its initial purpose was not to describe aging, but Sollich et al. (2017) have shown that the modeled systems that age under zero stress rapidly freeze into a preparation-dependent state; the initial stress does not fully relax. Within the same framework, Liu et al. (2018) studied athermal creep under a load Σ ≡ σ > Σ y and they, too, reported a strong dependence on the preparation. The initial distribution of stresses P(σ, t = 0) was taken as a proxy for the sample age insofar as, in real systems, aging results in stress relaxation and thus a narrower distribution P(σ, t = 0). For Σ slightly above the yield stress and long aging, there is first a power-law decayγ(t) ∼ t −µ (µ > 0) to a minimal value and then an acceleration up to the steady-state value. This evolution is consistent with several experimental measurements in bentonite suspensions and colloidal hard-sphere systems. But, contrary to expectations, the model exhibits a parameter-dependent (thus, non-universal) power-law exponent µ. Within the model, the first creep regime is dominated by the plastic activation of sites that have not yielded yet, which become rarer and rarer, until the memory of the initial configuration is lost and steady-state fluidization is achieved. This occurs at a fluidization time τ f that decreases as Σ increases, but in a non-universal way.
In conclusion, these few seminal papers proposing a mesoscopic approach to creep leave room for further exploration with EPM, for instance about the universality (or not) of the long-time response in thermal and athermal systems. It would also be interesting to determine if precursors can be defined to predict failure and, once the validity of EPM is established, to upscale the mesoscopic approach into a valid macroscopic description of the creep response.

IX. RELATED TOPICS
Amorphous solids seem to form a specific class of materials. However, the phenomenology exposed in the previous chapters suggests underlying theoretical connections with other problems. And, indeed, EPM are related to a spectrum of other models, notwithstanding physical differences, in particular in the interaction kernels. This section reviews, and attempts to compare to EPM, some of these related approaches, from mesoscale models for crystalline plasticity and elastic line depinning to fiber bundles, fuse networks and random spring models. The ample connections with seismology, hinted at in Sec. III, and tribology (Jagla, 2018;Lastakowski et al., 2015;Persson, 1999) -the latter being plausibly mediated by fracture mechanics (Svetlizky and Fineberg, 2014;Svetlizky et al., 2017) -will not be discussed here.
A. Mesoscale models of crystalline plasticity

Crystal plasticity
Like amorphous solids, driven crystalline materials respond elastically to infinitesimal deformations, via an affine deformation of their structure, but undergo plastic deformation under higher loading. To be energetically favorable, plastic deformation increments must somehow preserve the regular stacking of atoms. The question is whether it saves energy to jump to the closest regular structure ('switch neighbors'), rather than to keep on with the affine deformation of the current structure. For a perfect crystal, such a criterion would predict an elastic limit of around 5%.
Real crystals actually have a much lower elastic limit because they harbor structural defects, which were created at the stage of their preparation and which play a key role in the deformation. These defects in the regular ordering take the form of dislocations and grain boundaries separating incompatible crystalline domains. Dislocations are line defects obtained by making a half-plane cut in a perfect crystal and mismatching the cut surfaces before stitching them back together. Dislocations are similar to creases on a carpet in that they can glide across the crystal (and occasionally "climb" when they encounter a defect), thereby generating slip planes, in the same way as creases can be pushed across the rug to move it gradually without having to lift it as a whole. Grain boundaries also promote deformation; in these regions, gliding is facilitated by the mismatch-induced weakness of the local bonds. On the other hand, the presence of impurities, e.g., solute atoms in the crystal, may pin a dislocation at some location in space until it is eventually freed by a moving dislocation, which results in a dent in the stress vs. strain curve; this is the so called Portevin-Le Chatelier effect.
The stress field around a dislocation is well known (it decays inversely proportionally to the distance to the line) and the attractive or repulsive interactions between dislocations can also be rigorously computed. As a matter of fact, the elastic propagator used in EPM can be regarded as the stress field induced by four edge dislocations whose Burgers vectors sum to zero (Ben-Zion and Rice, 1993;Ispánovity et al., 2014;Tüzes et al., 2017). However, owing to the vast lengthscales separating the individual dislocation from the macroscopic material, it is beneficial to coarse-grain the description to the mesoscale, by considering the dislocation density field.

Models and results
Mesoscale dislocation models, which exist in several variants (Field Dislocation Model, Continuum Dislocation Dynamics), bear formal similarities with EPM.
Noticing that the plastic deformation induced by crystallographic slip generates an elastic stress field τ int (r) (via the very same elastic propagator as in EPM), Zaiser and Moretti (2005) separated this internal stress τ int (r) from the aspects more specific to dislocations and crystals and arrived at the following equation in 2D: where B, D, and G are material constants, τ ext is the externally applied stress, and ρ is the dislocation density. The last two terms on the rhs have no strict counterparts in EPM; they account for the mechanisms generated by interactions between dislocations that alter the stress required to set a dislocation in motion. The third term is a homogenizing term while the fourth one is a (ρ-dependent) fluctuating term; its dependence on the plastic strain γ may be used to effectively describe strain hardening effects due to the multiplication of dislocations. In EPM, such effects would belong to the rules that govern the onset of a plastic event. Armed with this model, the authors then studied the slip avalanches in order to explain the experimentally observed deformation patterns consisting of slip lines and bands, echoing the endeavors in this direction on the EPM side. They found scaling exponents for such avalanches that are comparable, but not strictly equal, to the mean-field exponents for the depinning problem; this difference is not unexpected, owing to the fluctuating sign of their elastic propagator, which is identical to the EPM one (see Sec. VI). Also, large avalanches are cut off due to strain hardening, which is one possible explanation for the macroscopic smoothness of the deformation. Contrasting with this macroscopically smooth situation, the deformation dynamics may feature strong intermittency, which points to collective effects. Power-law-distributed fluctuations have recently been evidenced in the acoustic emissions as well as in the stress vs. strain curves of loaded crystals (Weiss et al., 2015;Zhang et al., 2017a). These fluctuations may be "mild", with bursts superimposed on a relatively constant, seemingly uncorrelated fluctuation background, which is the case for many bulk samples, especially those with an fcc (face-centered cubic) structure. On the other hand, intermittency becomes dominant in hcp (hexagonal close-packed) crystals and in smaller samples, where large bursts dominate the statistics. Samples with fewer defects also tend to have "wilder" fluctuations. A meanfield rationalization of these phenomena considers the density ρ m of mobile dislocations and expresses its evolution with the strain γ as where A is a nucleation rate, C is the rate of annihilation of dislocation pairs, and D controls the intensity of the white noise ξ (Weiss et al., 2015). Notice that the latter is multiplied by ρ, owing to the long-ranged interactions between dislocations; the presence of multiplicative mechanical noise makes collective cascade effects possible. Such a model allowed the authors to capture the distinct types of fluctuations in the dynamics, from mild to wild, depending on the noise intensity D. More recently, Valdenaire et al. (2016) rigorously coarse-grained a fully discrete 2D dislocation picture into a continuum model centered on a kinetic equation for the dislocation density, with superficial similarities with the EPM equation of motion, Eq. (5).

Relation to EPM
Although the microscopic defects and the microscopic deformation mechanisms differ between crystals and disordered solids, the macroscopic phenomenology and, to some extent, the mesoscopic one share many similarities: Microscopic defects interact via long-range interactions and their activity is, in some conditions, controlled by temperature. Globally, the dynamics are highly intermittent at low shear rates and involve scale-invariant avalanches, as indicated, inter alia, by acoustic emission measurements on stressed ice crystals (Miguel et al., 2001). This intermittency is generically known as crackling noise (Sethna et al., 2001) and does not connect EPM only to crystal plasticity, but also to the fields of seismology and tribology.
The phenomenological similarity is paralleled by a proximity in the models. In some EPM, the stress redistributed by a shear transformation is actually described as the effect of a combination of dislocations (Ben-Zion and Rice, 1993;Ispánovity et al., 2014;Tüzes et al., 2017). Conversely, quadrupolar interactions may be directly implemented in mesoscale models of crystal plasticity, for instance in Eq. (2) of (Papanikolaou et al., 2012). More generally, the basic equations of evolution in the two fields look very much alike, and models sometimes seem to have bearing on both classes of materials (Shiba and Onuki, 2010). Rottler et al. (2014) numerically investigated the transition between the dislocation-mediated plasticity of crystals and the shear-transformation-based deformation of amorphous solids. They found that the directions of the nonlinear displacements under strain could be well predicted from the low-frequency vibrational modes and that polycrystals already behave comparably to glasses, despite their regular structure at the grain scale.
Nevertheless, the connection between crystals and disordered solids should not be overstated. Even though flow defects ("soft spots") in the latter might to some extent persist over rearrangements , on no account can they be assimilated to well identified structural defects moving through a crystal. Following from this discrepancy are the facts that, contrary to plastic rearrangements, dislocations are strongly dependent on the preparation of the material (which determines the dislocation density), and may be pinned by defects, annihilate through the merger of partials ("opposite" defects) or multiply.

The classical depinning problem
In several systems, an interface is driven through a disordered medium by a uniform external force. This interface can be a magnetic or ferroelectric domain wall, the water front (contact line) in a wetting problem, the fracture front, or even charge density waves and arrays of vortices in superconductors. In all these cases, the interplay between the quenched disorder (e.g., due to impurities) and the elastic interactions along the interface is at the root of a common phenomenology and a universal dynamical response.
If the external force is weak, the interface will advance and soon get pinned and unable to advance any further. If the force is strong enough, instead, the interface will overcome even the largest pinning centers, reaching a steady state of constant velocity. This is the well documented dynamical phase transition known as depinning. Beyond the transition itself, the literature now also describes the equilibrium configuration of the elastic line, several variations of the problem (short/long-range elasticity, different disorder types, etc.), thermally activated dynamical regimes and, in general, tackles the transport problem and its relation with the geometry of the interface. The interested reader is referred to one of the following self-contained works or reviews: (Agoritsas et al., 2012;Chauve et al., 2000;Ferrero et al., 2013;Fisher, 1998;Kolton et al., 2009).

Models
The most celebrated model to describe the depinning problem is the quenched Edwards-Wilkinson (QEW) equation. A d=1-dimensional interface without overhangs is driven by an external pulling force f . In the overdamped limit, its local displacement at time t, described by a single-valued function, h(x, t), obeys where c∇ 2 h(x, t) represents the elastic force due to the surface tension, the (quenched) disorder induced by impurities is encoded in the pinning force F p (x, h) and thermal fluctuations are included as a Langevin thermal noise ξ(x, t).
In general, two different kinds of disorder are considered: random bond disorder, in which the pinning potential is short-range correlated in the direction of motion (< V (h, i)V (h , j) >= δ ij δ hh ), and random-field disorder, where the pinning force is short-range correlated (thus generating correlations of the potential in the direction of motion, Of course, the QEW model just mentioned is minimal. Some of its variants take into account additional ingredients. For example, charge density waves and vortices involve a periodic elastic structure, in fracture and wetting the elastic interactions are long-ranged, and anharmonic corrections to elasticity or anisotropies could also be relevant. These features would call for a rewriting of Eq. (44) into a more general form involving an elastic interaction energy H el Remarkably, all these different problems, grouped in a few distinct universality classes, share the same basic physics, discussed in the following.

Phenomenology
The velocity-force characteristics <ḣ >= v(f ) is well known for the depinning problem (see Fig. 24a); the information conveyed by this "equation of state" is enriched by a vast analytical and numerical knowledge of universal properties at three special points: (i) equilibrium, i.e., f = 0; (ii) depinning, i.e., f = f c at T = 0; and (iii) fast-flow f f c . Around these points, at vanishing temperature, the steady-state interface h(x) displays a self-affine geometry (in the sense that it is invariant under dimensional rescaling, viz., h(ax) ∼ a ζ h(x)) above a microscopic length scale, with characteristic roughness exponents: (i) ζ eq , (ii) ζ dep , and (iii) ζ f f .

FIG. 24
The depinning picture. (a) Connection between transport and geometry in depinning. From (Ferrero et al., 2013). (i) Snapshot of a domain wall in a 2D ferromagnet. (ii) Typical velocity-force characteristics. (iii) Crossover lengths opt and av representing the optimal excitation and the deterministic avalanches, respectively. (iv) Geometric crossover diagram. (b) Steady-state structure factor S(q) of the line in the limit of vanishing temperature for different forces (curves are shifted for clarity). Adapted from (Kolton et al., 2006). (c) Comparison of the depinning and yielding critical transitions in a correspondence (v ↔γ), (f ↔ Σ).
Turning to transport properties, at equilibrium, the mean velocity is zero and the dynamics is glassy. When the applied force approaches zero, macroscopic movement can be observed only at finite temperatures and at very long times. Collective rearrangements on a scale of size opt ( opt → ∞ as f → 0) are needed in order to overcome barriers E b ( opt ) growing as E b ∼ θ opt , with θ > 0 a universal exponent related to the roughness by θ = d−2+2ζ eq . This is the creep regime. At the zero temperature depinning transition the velocity vanishes as v(f, Approaching f c from above the motion is very jerky and involves collective rearrangements of a typical longitudinal size av that diverges at f c . The avalanche size S, defined as the area covered by the moving interface, has power-law statistics, viz., At finite temperature, the sharp depinning transition is rounded, the velocity behaves as v(f c , T ) ∼ T ψ and the size av is finite at the transition. In the fast-flow regime f f c , the response is linear, viz., v ∼ f . Here impurities generate an effective thermal noise on the interface. Therefore, the fast-flow roughness corresponds to the Edwards-Wilkinson roughness ζ f f = (2 − d)/2.
One of the remarkable lessons learned from this simple model is the possibility to relate transport and geometry. If the applied force f lies in between two of the above mentioned reference points, the interface geometry [in particular the roughness exponent, see Fig. 24(b)] depends on the observation scale and its relative position compared to the characteristic lengths ( opt , av ,. . . ). Granted that one knows the functional dependencies of these characteristic lengths with f and the velocity-force characteristics for a given system, transport properties (which intrinsically pertain to the dynamics) can be deduced from the static interface geometry, and vice-versa.

Similarities and differences with EPM
The manifest qualitative similarity between the yielding transition and the depinning one has enticed many researchers to look for a unification of these theories. The analogy has promoted the vision of yielding as a critical phenomenon and has given rise to interesting advances, but, in our opinion, the (misguided) belief in a strict equivalence of the problems has been deceptive in some regards.
To stay on firm ground, a formal approach consists in finding an EPM analog to the depinning equation, Eq. (45). In the stress-controlled situation (with applied stress Σ ext ), Weiss et al. (2014) (Eq. S3 of the Supplemental Information) and Tyukodi et al. (2016b) thus proposed to substitute the EPM equation of motion [Eq. (5)] with where U el [ pl ] ≡ − 1 2 drdr pl (r)G(r − r ) pl (r ), with G the elastic propagator, and P(x) denotes the positive part of x (x if x > 0, 0 otherwise). In so doing, the deformation of an amorphous solid is mapped onto a problem of motion through an (abstract) disordered space for the pl -manifold pulled by the 'force' Σ ext . The positive part P in Eq. (47) creates genuine threshold dynamics; it has no direct counterpart in the depinning equation but was argued by Tyukodi et al. (2016b) not to be a core dissimilarity between yielding and depinning.
This formal similarity between the two classes of phenomena seems to buttress the application of results from the depinning problem (hence mean field, owing to the long range of the elastic propagator) to the question of, e.g., avalanche statistics in disordered solids (see Sec. VI). However, the following differences must be borne in mind.
First, and perhaps foremost, as often mentioned along the present review, the interaction kernel in depinning problems is positive, whereas the quadrupolar elastic propagator G used in EPM has positive and negative bits. This has profound consequences on the critical behavior at the yielding transition observed in EPM, in particular with respect to the possibility of strain localization and the avalanche statistics. Furthermore, while in depinning v vanishes at f c as v ∼ (f − f c ) β with typically β < 1, the strain rateγ does so at the yielding transition asγ ∼ (Σ − Σ c ) β with β > 1, as schematically shown in Fig. 24(c). Note that, if the systems were at equilibrium, this difference in the value of β would imply a change in the order of the continuous phase transition. Other consequences can be deduced from the general scaling relations proposed by Lin et al. (2014b) [Supplementary Information], which are claimed to encompass the depinning and the yielding cases (not all these relations are strictly obeyed in finite-dimensional EPM): Here, d f is the fractal dimension of the avalanches, z is the dynamical exponent, ν is the exponent controlling the divergence of the correlating length at the transition and α k is the dimension of the elastic interaction kernel. In EPM α k = 0 and d f < d so that β > 1. In depinning, α k = 2 for short-ranged elasticity and α k = 1 for long-ranged elasticity, θ = 0 and d f ≥ d.
Secondly comes the question of the nature of the disorder in the pinning force F p . In elastic depinning models, regardless of how realistic the chosen correlations of F p are, the origin of the disorder is generally extrinsic. More precisely, it reflects the disorder of the substrate on which the elastic manifold advances, hence F p = F p (h). On the other hand, in the yielding phenomenon, as stated by Papanikolaou (2016), 'the pinning disorder for every particle originates in the actual interface that attempts to depin (other nearby particles); a disordered solid pins itself during deformation'. Therefore, it is inaccurate to consider that F p only depends on the local value of pl . In particular, a given system will not encounter the same pinning forces F p along its deformation between, say, pl = 0 and pl = 1 if it is sheared slowly and if it is sheared fast. Typically, at high shear rates, the potential energies of the inherent structures of the material are higher (as evidenced by the variations of potential energies of the inherent structures with the shear rate, in atomistic simulations). This dependence should impact theγ = f (Σ) curve at finite shear rates.
Lastly, the EPM equation of motion [Eq. (5)] cannot always be reduced to an expression akin to Eq. (47), because of the memory effects contained in the plastic activity variable n.
Let us now mention a subclass of problems that may be more closely related to EPM: the so called "plastic depinning". This phenomenon is observed for example in particle assemblies driven over random substrates whenever irreversible plastic deformations actually occur, or in charge density wave problems. Unfortunately, this comparison has been much less exploited by the amorphous solids community, even though the connection was very recently pointed out in Reichhardt and Reichhardt (2016)'s review.
To conclude on the topic, there undoubtedly remains much to be learned from the 30 + years of studies on depinning phenomena. Some intriguing open questions left from this comparison are the following: Are the transport properties of driven amorphous solids related to geometrical properties, as they are in elastic manifolds? Is it possible, for example, to infer from a picture at which strain-rate a dense emulsion is being sheared? C. Fiber bundle, fuse networks and continuum models for the study of cracks and fracture

Brief introduction to cracks and fracture
In partial overlap with the scope of EPM, the question of the failure of hard solids under loading, e.g. in tension, has attracted much attention over the last centuries. Pioneering in this respect, as recalled by Alava et al. (2006), is Leonardo da Vinci's observation that, if one loads a metal wire in tension with a weight, it will fail more readily if it is longer, for the same cross section; this runs counter to basic continuum mechanics predictions for a uniform medium. In fact, the failure of brittle solids, in particular rocks, is ascribed to the growth and propagation of pre-existing cracks (at the scale of the crystalline grains constituting the material) or, more generally, defects.
If one considers an individual crack in a homogeneous medium, according to Griffith (1921)'s criterion, its growth hinges on a competition between a surface energy term averse to the opening of solid-air interfaces and an elastic energy term favoring its growth and thereby reducing the elastic energy stored in the bulk. For example, for a single elongated elliptic crack of length a in a 2D medium, the sum of these competing terms reads where E is the Young modulus of the material, γ is the interfacial energy, and Σ is the applied stress. Thus, the evolution of the crack depends on the sign of the derivative dE T da (Alava et al., 2006). However, cracks very seldom have so simple a geometric shape. Roughly speaking, owing to the presence of heterogeneities, the crack will zigzag around hard spots. This will result in undulations and protrusions in the post-mortem fracture surface, which exhibits a selfsimilar (fractal) pattern: If the surface height at a point (x, z) is denoted by h(x, z), the root mean square fluctuation w(l) of the height in a region of size ∆x ≈ ∆z ≈ l obeys where ζ ⊥ is the (out-of-plane) Hurst exponent, or roughness exponent. Interestingly, this exponent seems to be weakly sensitive to the material or the loading, with values centered around ζ ⊥ 0.8 and early claims of universality (Bouchaud et al., 1990). The fractal dimension d f of the surface is then related to ζ ⊥ via d f = 3 − ζ ⊥ for 3D fracture. While the material is being fractured, the crack propagates along a rough, scale-invariant frontline (see Fig. 25a), characterized by the in-plane roughness exponent ζ . Roughness bears practical importance, since it modifies the scaling of the surface energy term. Let us mention two subtleties. First, the exponents ζ and ζ ⊥ are not independent (Ertaş and Kardar, 1994). Second, ζ ⊥ might in fact mix two distinct exponents, insofar as Ponson et al. (2006)'s fracture experiments on silica and aluminium alloys hint at anisotropic height variations in the fracture plane, with distinct behaviours along the front line and along the crack propagation direction.
In addition to being spatially nontrivial, the propagation of the crack front also displays marked variations in time. The associated dynamics is highly intermittent and involves avalanches of events which span a broad range of energies. Indeed, the crackling noise emitted during these events has a power-law power spectrum, for instance in composite materials (Garcimartin et al., 1997). For instance, the crack produced when tearing apart two sandblasted Plexiglas sheets stuck together through annealing undergoes a stick-slip motion at small scales that is reminiscent of dry solid friction (Måløy and Schmittbuhl, 2001), which in turn may tell us about earthquake dynamics (Svetlizky and Fineberg, 2014).
At this stage, a discrepancy with respect to soft solids ought to be mentioned: In (rock) fracture, the microruptures very generally do not have time to heal on the time scale of the deformation; without recovery process, the material is thus permanently damaged. However, the crack velocity may still have an influence on the dynamics of the process owing to the finite duration of the avalanches.

Fiber bundles
Arguably, the simplest way to model fracture is to consider two blocks bound by N aligned fibers. These fibers share the global load and break irreversibly when their elongation x exceeds a randomly distributed threshold; this is the basis of fiber-bundle models (Herrmann and Roux, 2014). In democratic fiber bundles, the load of broken fibers is redistributed equally to all survivors. Analytical progress is possible in this intrinsically mean-field model.
In particular, it is easy to show that, on average, when the bundle is stretched by x (with x = 0 the reference configuration), a fraction C(x) of fibers have broken, where C(x) is the cumulative distribution of thresholds, and the total load (normalized by the initial number of fibers, of stiffness κ each) readsf (x) = κx [1 − C(x)]. It follows that the maximum strength per fiber of the bundle is, on average, If one pulls on a given bundle, however, the load f will not evolve along the smooth average profilef (x), but along a rugged profile {f (x k ) ≡ κx k [1 − C(x k )] , k = 1 . . . N } due to the randomness of the thresholds x 1 x 2 . . . x N , sorted according to the order of failure. The f (x k ) thus perform a random walk in "time" k with a time-dependent bias f (x k+1 ) − f (x k ) (Sornette, 1992). If, starting from a stable situation, the rupture of the k-th bond leads to S additional failures, viz., f (x k+i ) < f (x k ) for i between 1 and S (but not for i = S + 1), an avalanche of size S will occur under fixed load. Noting that (i) this is a problem of first return for the walker f (x k ) [or, equivalently, of survival close to the absorbing boundary f = f (x k ) − ], and that (ii) close to global failure f ≈ f c the random walk is unbiased, i.e., f (x k+1 ) − f (x k ) = 0, Sornette (1992) showed that the distribution of avalanche sizes s obeys More precisely, for a uniform distribution of thresholds between x l and 1, the distribution reads where the cutoff size S cut ≡ 1 2(1−2x l ) 2 diverges at the critical point x l = 1 /2 (Pradhan et al., 2005). For x l < 1 /2, the fiber fails gradually as the loading is increased, whereas for x l > 1 /2 failure occurs all at once. A parallel may here be drawn with the discussion about the brittle-to-ductile transition in amorphous solids in Sec. V. For x l 1 /2, the power law with exponent τ = 3 /2 of Eq. (52) is recovered for S S cut , whereas for S S cut the random walk of the f (x k ) is biased upward and a steeper power law is obtained, with an exponent 5 /2. The scaling p(S) ∼ S −5 /2 is found generically if all avalanches since the start of the deformation (x = 0) are taken into account (Hemmer and Hansen, 1992). The gradual shift to an exponent τ = 3 /2 then signals imminent failure. Interestingly, the power-law behavior fades out in favor of a much faster decay of p(S) if the load released by broken fibers is redistributed locally to the first neighbors only, instead of being shared by all intact fibers (Kloster et al., 1997).

Fuse networks
Unfortunately, the picture promoted by mean-field or 1D fiber bundles is incapable of describing the heterogeneous and anisotropic propagation of cracks. Extending the approach to higher dimensions, fuse networks connect lattice nodes (say, nodes i and j) by fuses of conductance K ij that break past a threshold x ∈ [0, 1], thereby burning the fuse (K ij → 0). To take an example, the distribution p of the thresholds can be set as a power law, p(x) ∼ x θ with θ > 0. The voltages V i are imposed at two opposite edges of the system, as depicted in Fig. 25c. The Hamiltonian of the system reads where the sum runs over all adjacent nodes (i, j). Note that, if the K ij are constant, then the model can be viewed as a discretization of Poisson's equation in the vacuum, ∇ 2 V =0. Fuse networks are thus closer to EPM than fiber bundles, insofar as the stress redistribution when one fuse burns (in the pristine network) is strongly anisotropic, with a shielding of the current fore and aft and an enhancement sideways (Barthelemy et al., 2002;Rathore, 2016). It can then be understood that failure occurs along a line of burnt fuses, the "crack" line, provided that there is finite disorder (θ > 0) and the network is large (Shekhawat et al., 2013). Besides, in a 2D fuse network, Hansen et al. (1991) computed a roughness exponent ζ approximately equal to 0.7 for weak disorder, not far from experimental values for fractured surfaces ζ ⊥ ≈ 0.8 (note that ζ = ζ ⊥ in 2D). The intact region appears in black and the image-processed front line is shown in (b). The roughness along the propagation (z ) direction has a power-law spectrum characterized by the roughness exponent ζ . From (Schmittbuhl and Måløy, 1997). (c) Sketch of the random fuse network and (d-e) failure process for distinct probability density functions for the thresholds, p(x) ∼ x θ . The crack has been colored in red. From (Shekhawat et al., 2013).
The expression of the Hamiltonian in Eq. (54) evokes a random bond Ising model; the equivalence is formally exact if the voltages are restricted to the values ±1, and the thresholds are infinite, thus making bonds unbreakable (perfect ductility). These differences are not negligible in any way. Indeed, the interactions between nodes are thereby much reduced, in spatial extent and magnitude; by contrast, in random fiber or fuse models, the impact of breaking a bond is magnified close to failure, owing to the small number of intact bonds which will share the load. Nevertheless, the process of fracture can be mimicked in the random Ising models by imposing spin +1 (-1) on the left (right) edges of the sample and monitoring the interface line between the +1 and -1 domains. Rosti et al. (2001) studied the probability that this interface passes through an artificial "notch", i.e., a segment in which the bond strengths K ij have been set to zero, and observed a transition from low to high probabilities as the notch length was increased above a disorder-dependent threshold value. Similar results were obtained in experiments in which sheets of papers with pre-cut notches were torn.

Spring models
From a mechanical perspective, should one replace the voltage V i in Eq. (53) with the displacement u i at node i, viz., the interpretation of the Hamiltonian as the energy of a network of random springs of stiffness K ij will become apparent. The x, y, and z components of the dispacements in H nc decouple, so that model is actually scalar (De Gennes, 1976). However, it features noncentral forces: the force exerted by j on i is not aligned with e ij . A more consistent description of a network of nodes connected by harmonic springs relies (to leading order) on the Hamiltonian On a triangular lattice, with bonds of uniform strength K ij = 1, the continuum limit of this Hamiltonian represents an isotropic elastic medium with a Poisson ratio of 1 /3 in 2D and 1 /4 in 3D (Monette and Anderson, 1994). As bonds are gradually removed in a random fashion, the initially rigid system transitions to a non-solid state with vanishing elastic moduli at a critical bond fraction p c . Such a transition is also observed with the models based on the scalar Hamiltonian H nc or the noncentral Hamiltonian H nc , although at a distinct fraction p c . Somewhat surprisingly, the scalings of the shear and bulk moduli with the fraction of bonds p around p c differ between the H c and H nc -based models; the discrepancy stems from the distinct symmetries, in the same way as the Heisenberg model differs from the Ising model (Feng and Sen, 1984). The distinction subtly differs from the dichotomy between scalar and tensorial EPM, in that the EPM propagator is always derived from the same constitutive model (tensorial continuum elasticity); the scalar description simply discards some tensorial components at the end of the day. Regarding the avalanches of ruptures close to the point of global failure, i.e., under loading f ≈ f c , Zapperi et al. (1999) claimed that both the random fuse network of Eq. (53) and the central-force spring model of Eq. (55) (supplemented with bond-bending forces) fall in the universality class of spinodal nucleation, in that the avalanche sizes S are distributed according to p(S) ∼ S −τ Φ [S (f c − f )] , where τ = 3 /2 and Φ is a scaling function. Zaiser et al. (2015) also found that fuse networks yielded results similar to spring models regarding the initiation of failure, with localized correlations in the damage patterns.
We conclude this section on spring models with a historical note referring to the fact that such models had in fact been pioneered by De Gennes (1976) to tackle the "converse problem", namely, gel formation (e.g., through cross-linking): Instead of gradually destroying bonds, he cranked up the fraction p of bonds by randomly connecting pairs of neighbours until bonds percolated throughout the system; this occurred at a critical fraction p c , supposedly corresponding to gel formation. In any event, the nature of the transition associated with the random depletion (or creation) of bonds, which pertains to percolation, is distinct from what is observed in random fuse or spring networks. In the latter models, the disorder in the yield thresholds bestows critical importance to the spatial redistribution of stresses following ruptures. This distinction is at the origin of different scaling relations, e.g., between the failure force and the system size (Hansen et al., 1989).

Beyond random spring models
Refinements have been suggested to bring random fuse (or spring) networks closer to models of material deformation and fracture. First, the irreversible breakage of the fuses past a threshold mirrors perfectly brittle fracture. At the opposite end, perfect plasticity is mirrored by the saturation of the fuse intensity past a threshold. But a continuum of possibilities can be explored between these extreme cases, whereby the conductivity of the fuse is decreased to mimic partial weakening, similarly to what can be done in EPM.
Another limitation of the models stems directly from the description of the bonds on a regular lattice: let alone the presence of soft modes in several cases, the (H c -based) central-force model, discretized on a triangular lattice, displays an anisotropic tensile failure surface (despite an isotropic linear response), with an anisotropy ratio of 50% (Monette and Anderson, 1994). These deficiencies can be remedied in part by complementing the spring-stretching energies in H c with bond-bending energies. This refinement leads to an isotropic elastic medium with adjustable Poisson coefficient and a more isotropic failure surface.
As with EPM, the following step in the endeavor to refine the description led to the introduction of a finiteelement approach, which relies on a continuum description down to the scale of one mesh element. The equations of inhomogeneous elasticity are solved and a damage (of magnitude D) is introduced by reducing the local elastic constant E → (1 − D)E whenever the local stress exceeds a threshold value. The process can evolve into avalanches, and eventually to a vanishing of the elastic resistance through the propagation of a fracture through the system. Incidentally, this mechanism had first been implemented by Zapperi et al. (1997) using a fuse network with damage operating on the fuse resistances; the model displayed scale-free behavior with power-law distribution of event sizes, P (S) ∼ S −τ with τ 1.2. Amitrano et al. (1999) refined the modeling approach by using a pressure-modified (Mohr-Coulomb) criterion for the onset of plasticity, viz., C + σ n tan φ − σ < 0, where C represents the cohesion of the material, σ n and σ are the normal and shear stresses, respectively. A transition from brittle failure with very localized damage (at low internal friction angle φ, i.e., little sensitivity of the yield criterion to pressure) to ductile with diffuse damage (at large φ) was observed. At low φ the damage around a single event is similar to the stress redistribution considered in EPM, while for large φ it becomes much more directional. The transition from ductile to brittle shares qualitatively similarities with the strain localization transition, but the control parameter is different from those discussed in Sec. V.C.
In the case of large φ and brittle failure, a description of compressive failure under uniaxial stress as a critical phenomenon analogous to depinning was proposed by Girard et al. (2010) and elaborated by Weiss et al. (2014). The interpretation in terms of a criticality notably affords a detailed description of size effects on the critical stress (Girard et al., 2012).

X. OUTLOOK
In the last ten years, EPM have become an essential theoretical tool to understand the flow of solids. Starting from elementary models intended to reproduce earthquake dynamics, they have blossomed into more refined approaches that have helped rationalize many experimentally observed features, at least at a qualitative level, and unveil new facets of the rheology of these materials. Future developments in the field can be expected in a number of directions, following current experimental and theoretical interests.
In rheology, considerable attention has recently been devoted to the study of transient regimes. For instance, one can mention the study of the load curves at fixed shear rateγ, that can exhibit stress overshoots depending on the initial preparation and non-trivial scalings of the time to reach the stationary state withγ. Other examples include creep under imposed stress, the dynamics of relaxation and the residual stresses after sudden cessation of the driving, and oscillatory regimes. In the latter category, the Large Amplitude Oscillatory Strain (LAOS) protocol probes the nonlinear behavior and the frequency dependent one at the same time, and therefore involves a complex interplay between plastic deformation and internal relaxation. Reproducing the complex response of particular systems under such protocols is particularly challenging for simple models. Several issues could be investigated within the framework of EPM, such as the onset of tracer diffusion as the amplitude of the oscillatory strain is increased, or the fatigue behavior leading to failure. Recently, it was suggested that the LAOS protocol could induce strain localization in systems with a monotonic flow curve, based on a study of a spatially resolved version of the soft glassy rheology model, presented in Sec. IV.D.1 Fielding, 2016, 2018). Creep (see Sec. VIII) is an equally challenging phenomenon; a recent mean-field EPM illustrated its very strong sensitivity to the initial conditions (Liu et al., 2018).
A more unexpected emerging avenue is the study of systems with internal activity, such as living tissues or dense cell assemblies. The general ideas exploited for the description of amorphous systems can indeed be expanded to incorporate new types of events, such as cell division (assimilated to a local anisotropic dilation) and cell death (local isotropic contraction). At the mean-field level, Matoz-Fernandez et al. (2017) and others conducted a first analysis along these lines. For further exploration of the collective behavior resulting from the interplay between cell division, apoptosis, locomotion, and contractility, as well as the mechanosensitivity of these processes, an EPM describing all these ingredients at the same time would be invaluable. At present, new experimental tools are providing information on the statistical fluctuations in such systems, which will allow to calibrate these models.
From the viewpoint of statistical physics, the yielding transition described by EPM stands as a new type of dynamical phase transition, with specificities that are still to be understood in extent. Considerable efforts have been devoted to the theoretical study of the related problem of the depinning transition (Sec. IX.B). In the latter case, (mostly) exact exponents, scaling functions, and avalanche shapes were derived using scaling analysis and renormalization techniques. For the yielding transition, the (slow) process of consensus building has not converged yet, but there are reasons to believe that the results on avalanche statistics obtained in the depinning problem cannot be directly transposed to this field, because the propagator controlling stress redistribution is partly negative, which affects the density of sites close to yielding. Whether this features only induces an effective dimensional reduction, leaving us in a well known universality class but for d < D, or whether it exhibits a completely distinct set of critical exponents, still needs to be clarified. Scaling relations between critical exponents have been proposed (Aguirre and Jagla, 2018;Lin et al., 2014b) and tested in diverse EPM, but analytical calculations beyond mean field are scant. Recent efforts to relate EPM to (better known) problems of motion through a disordered landscape open new vistas for the understanding of yielding and transport properties under slow driving (Jagla, 2017b), but there is still no consensual theory explaining the flow exponents (the low-shear-rate rheology). The situation is somewhat similar on the experimental side: The depinning phenomenon has benefited from a very detailed experimental characterization in various systems (magnetic domain walls, contact lines, vortices), including avalanche statistics and shapes, which has permitted comparison to the theory. Amorphous plasticity is not on quite so good a footing, with only a few attempts to characterize the distribution of stress drops in deformed systems. The situation is however improving, thanks to several recent efforts, e.g. those combining mechanical deformation and confocal microscopy in colloidal glasses.
The foregoing discussion is related to the critical aspects of the yielding phenomenon, discussed in Sec. VI and VII. In a number of real systems (Bonn et al., 2017), the onset of flow is in fact discontinuous and implies a coexistence between flowing and immobile states. EPM and other theoretical studies have proposed possible mechanisms that may influence the continuous or discontinuous character of the transition (see Sec. V). Nevertheless, it turns out to be experimentally difficult to control the transition in a systematic way by changing some experimental parameter. Wortel et al. (2016)'s work on weakly vibrated granular media represents a notable exception, insofar as the intensity of external shaking could be used to continuously tweak the flow curve towards nonmonotonicity. The ensuing emergence of a critical point at a finite driving rate has scarcely been addressed in the literature and could be analyzed with the EPM approach. Similar systems of vibrated grains have also permitted the experimental realization of a Gardner transition (Seguin and Dauchot, 2016), a transition which may be important for the theory of glasses and which has been associated with shear yielding (Urbani and Zamponi, 2017). On a related note, cutting-edge atomistic simulations suggest that the ductility or brittleness of the yielding phenomenon hinges on the initial preparation of the glass, rather than the microscopic interactions between particles or the dynamics (Ozawa et al., 2018); this puts EPM in the forefront for the study of these questions, but requires them to establish adequate proxys for the initial stability of the glass.
These prospective lines of research have hardly been explored using EPM. So, for all our efforts to articulate a comprehensive view of the state of the art here, we can only wish that this review will soon need to be updated with insightful results in these new avenues.