Yielding, shear banding and brittle failure of amorphous materials

Widespread processes in nature and technology are governed by the dynamical transition whereby a material in an initially solid-like state then yields plastically. Major unresolved questions concern whether any material will yield smoothly and gradually (ductile behaviour) or fail abruptly and catastrophically (brittle behaviour); the roles of sample annealing, disorder and shear band formation in the onset of yielding and failure; and, most importantly from a practical viewpoint, whether any impending catastrophic failure can be anticipated before it happens. We address these questions by studying the yielding of slowly sheared athermal amorphous materials, within a minimal mesoscopic lattice elastoplastic model. Our contributions are fourfold. First, we elucidate whether yielding will be ductile or brittle, for any given level of sample annealing. Second, we show that yielding comprises two distinct stages: a pre-failure stage, in which small levels of strain heterogeneity slowly accumulate, followed by a catastrophic brittle failure event, in which a crack quickly propagates across the sample via a cooperating line of plastic events. Third, we provide an expression for the slowly growing level of strain heterogeneity in the pre-failure stage, expressed in terms of the macroscopic stress-strain curve and the sample size, and in excellent agreement with our simulation results. Fourth, we elucidate the basic mechanism via which a crack then nucleates and provide an approximate expression for the probability distribution of shear strains at which failure occurs, as determined by the disorder inherent in the sample, expressed in terms of a single annealing parameter, and the system size.


I. INTRODUCTION
Amorphous materials [1,2] include yield stress fluids [3,4] and soft glassy materials [5] such as dense colloids, emulsions, foams and microgels, as well as hard materials such as molecular and metallic glasses [6,7]. When experiencing low loads or small deformations, such materials typically behave in a solid-like way. At higher loads or larger deformations, they then yield plastically. Indeed, numerous processes in nature and technology are governed by the dynamical transition whereby a material in an initially solid-like state yields, whether suddenly or gradually. Examples include the restarting of a pipeline of waxy crude oil [8]; the rising of bubbles in radioactive sludge [9]; the spreading of fresh cement [10,11]; the yielding of food during chewing [12]; the catastrophic material failure of metallic glasses [7]; geological processes such as landslides and lava flows [13,14]; and the reshaping of biological tissue under internal stresses caused by active processes such as cell division or cell death [15][16][17].
Alongside the complex fluids and soft solids just described, some notably similar phenomenology is seen in harder materials such as polymeric and metallic glasses [6,[44][45][46][47]. These likewise show initially solidlike behaviour at low strains, followed by yielding with a strong associated tendency for strain localisation and shear banding [7]. But whereas yield stress fluids can then flow indefinitely post-yield without losing their ability to return to a solid state when later unloaded, in metallic glasses the formation of shear bands typically leads to brittle yielding and catastrophic sample failure, thereby limiting the material's strength. Different scenarios for its onset have been discussed [7,48,49]. These include homogeneous nucleation [50], in which local plastic events triggered within the material under shear spatially cooperate to form a system-spanning shear band, once a percolation threshold for their link-up has been achieved as a function of increasing stress. The stress needed for this can be rather large, however, so in prac-arXiv:2103.06782v2 [cond-mat.stat-mech] 22 Nov 2022 tice shear bands often nucleate heterogeneously at stress concentrations, for example around any indentations on the sample's surface [51][52][53]. The important role of sample size in modifying shear banding formation and brittle yielding has also been studied [54,55]. In particular, the yield strength of metallic glasses has been found to increase with decreasing sample size: a phenomenon suggested to arise from the reduced number of sites for shear band nucleation as system size decreases.
Theoretical studies of yielding in amorphous materials have ranged from microscopic through mesoscopic and up to the continuum level. In the context of continuum modelling of the rheology of complex fluids, it can be predicted within a minimal set of generalised constitutive assumptions that a state of initially homogeneous shear must generically become unstable to the formation of shear bands once the maximum of stress as a function of strain is attained, and the material first starts to flow [56][57][58][59][60][61][62]. For samples subject instead to a constant applied shear stress, shear bands are likewise generically predicted to form as the initial creep regime terminates and the sample yields. These predictions have been confirmed experimentally [18][19][20][21][22][23][24][25], by molecular simulations [63][64][65][66][67][68][69] and in several specific constitutive models of complex fluids rheology [61,62,[70][71][72][73][74][75][76][77].
At the level of microscopic and mesoscopic modelling, much attention has been devoted to understanding the dynamical properties of the finally flowing shear state, in which the stress has attained a (statistically) steady state as a function of applied strain. (Any such state may or may not be attainable in practice for a given sample, being sometimes pre-empted by catastrophic sample failure.) A unifying picture has emerged in which local plastic events triggered within the material by the applied shear cooperate to form avalanches that flicker intermittently across the system, with power-law statistics. Their dependence on strain rate, system size and temperature has been carefully characterised [78][79][80][81][82][83][84][85][86][87][88][89].
Beyond this steadily flowing state, attention is increasingly turning to understanding the initial onset of yielding after a shear is first switched on. Several approaches have been put forward, studying dynamical yielding within a replica field theory [90][91][92]; as a critical point in an elastoplastic model [74,93]; as a directed percolation transition [64,65]; within a random first order transition theory for the glass transition [94]; as a Gardner transition [95]; as a spinodal point [96]; and within particle simulations that seed initial weak spots [97]. Yielding in oscillatory shear has been studied using particle simulations [98][99][100], in elastoplastic models [101] and in energy landscape models [102]. Yielding following creep under an applied load was studied using particle simulations in Ref. [103]. Microscopic precursors to yielding have recently been observed experimentally in soft materials [104][105][106].
With this backdrop, it is clear that yielding represents a problem of central importance to several areas of physics, with the observation of common phenomenolo-gies across multiple classes of material stimulating a search for universal explanations. In the twin contexts of soft matter physics and fluid dynamics, yielding is crucial to the question of how complex fluids flow. In statistical physics, it represents a non-equilibrium phase transition that is only just starting to be understood from first principles. In materials physics, it is core to understanding a material's ultimate strength in principle, and to improving material performance in practice. In active matter, its role in the reshaping of biological tissues remains to be elucidated.
A question of key significance to the practical performance of a material is whether it will yield smoothly and gradually ("ductile" behaviour) or instead fail abruptly and catastrophically ("brittle" behaviour) [107][108][109]. The influence of annealing and disorder in the sample prior to shear, and its role in determining the formation of shear bands during yielding and failure, is increasingly being appreciated [69,[107][108][109][110][111]. And from a practical viewpoint, perhaps the most important question is whether any impending sudden catastrophic failure can be anticipated, before it actually occurs, in terms of an identifiable material property. In this work, we address these questions within a lattice elastoplastic model [2] that contains only minimal assumptions, and should therefore capture the mechanics of quasistatically sheared athermal amorphous materials in a universally generic way.
Our first contribution will be to carefully elucidate whether yielding is ductile or brittle, for any given level of sample annealing prior to shear. For highly annealed samples we find brittle yielding for all system sizes. In contrast, for poorly annealed samples we demonstrate an important dependence of the nature of yielding on the size of the sample being sheared, with ductile yielding for small samples and brittle yielding only for sizes larger than have previously been studied theoretically. We thereby resolve an apparent contradiction between the recent work of Refs. [107,108], on the one hand, and that of Ref. [109] on the other.
Indeed, recent studies of mean field elastoplastic models of quasistatically sheared athermal amorphous materials [107,108] have suggested ductile and brittle yielding to be separated by a "random critical point". In this scenario, the underlying stress-strain relation, Σ(γ), is posited to display a qualitative change in form as the degree to which the sample has been annealed prior to shear increases: for poorly annealed samples, Σ(γ) has an overshoot but not an overhang, and yielding proceeds in a smoothly ductile way, whereas for highly annealed samples, it additionally develops an overhang, leading to catastrophic brittle yielding, with the stress dropping precipitously once this overhang is reached. Ductile and brittle yielding are thereby separated by a critical point, at which the stress-strain curve switches between these two qualitatively different forms with increasing annealing. Particle simulations and studies of two-dimensional (2D) lattice elastoplastic models beyond mean field have been argued to be consistent with this scenario [107].
In contrast, calculations performed beyond mean field within an athermal 1D elastoplastic model [109] reported a discontinuous stress drop for any level of sample annealing, however small, provided a stress drop indeed exists. They thereby predicted that yielding will always be brittle in slowly sheared athermal amorphous materials, contradicting the critical point scenario of Refs. [107,108].
Our work demonstrates that the scenario of Ref. [109], explored in that work in a simplified 1D elastoplastic model, also obtains in the limit of infinite system size in 2D lattice elastoplastic models comprising N × N elastoplastic elements. We achieve this by simulating a range of sample sizes that has not (to this author's knowledge) been attained in any previous study, from N = 64 to N = 4096. We thereby suggest that evidence for the random critical point scenario of Refs. [107,108], beyond mean field, may be a finite size effect, potentially amenable to further finite size analysis along the lines of that in Ref. [112].
Our next contribution will be to show that yielding comprises two distinct stages: a pre-failure stage, in which small levels of strain heterogeneity slowly accumulate across the sample, followed (in large or highly annealed samples) by a catastrophic brittle failure event, in which a macroscopic shear band suddenly propagates across the sample.
In the pre-failure regime, we provide an exact analytical expression for the slowly growing level of strain heterogeneity across the sample, expressed in terms of the macroscopically measurable stress-strain curve and the sample size. We show this to be in excellent agreement with our simulation results. We also perform the counterpart analysis for a mean field elastoplastic model, again obtaining excellent agreement between our analytical predictions and simulation results.
Finally, we shall elucidate the mechanism of catastrophic material failure, in which a shear band nucleates and spreads quickly across the sample. Our basic observation is that a single element first yields plastically and elastically propagates its stress to other elements via the Eshelby propagator, creating further yielding events in turn, which cooperatively percolate along a line [113]. In simulating a 2D lattice with periodic boundary conditions, our findings necessarily pertain to homogeneous nucleation. Although in practice failure often occurs via heterogeneous nucleation of shear bands at stress concentrations around surface imperfections, elucidating the basic physical mechanisms of homogeneous nucleation is an essential pre-requisite to progress in also understanding heterogeneous nucleation, and to determining the ultimate theoretical limit of material strength in the absence of surface imperfections.
We provide an expression for the probability distribution of shear strains at which a shear band appears and material failure occurs, in terms of the disorder inherent in the sample, as determined by the degree of annealing prior to shear. This expression shows good agreement with our simulation results for highly annealed samples, with progressively less good agreement for less well annealed samples. Importantly, this analytical result for the strain at which catastrophic failure occurs depends only on a single parameter characterising the level of initial sample annealing, together with the system size. Our findings thereby provide a possible route to predicting the brittle failure of highly annealed amorphous materials, before it actually occurs. This is an important prediction that we hope will stimulate further work to validate or disprove it, both experimentally and in particle simulations of soft and hard amorphous materials.
For highly annealed samples, we demonstrate a weak progression towards increasing sample strength with decreasing sample size, consistent with the trend observed in metallic glasses [55]. The finding of shear band nucleation in our 2D elastoplastic model is consistent with the insight of Wyart and co-workers in Ref. [108], who outlined arguments for band nucleation in an elastoplastic model based on classical ideas of fracture mechanics [114].
The paper is structured as follows. In Sec. II we introduce two elastoplastic models to be studied: a 2D lattice model with force balance, and its mean field counterpart. These models have been widely studied in the existing literature. In Sec. III we present our novel results for the predictions of these models for yielding in athermal quasistatic shear. We start in Sec. III A by showing that, beyond mean field, yielding is always brittle in the limit of large system size. In contrast, for small systems (or in mean field), yielding is brittle only for highly annealed samples, and instead appears ductile for poor annealing. In Sec. III B we elucidate the physics of the "pre-failure" regime, in which small amounts of strain heterogeneity slowly accumulate as a result of local plastic events that are initially largely uncorrelated. Sec. III C concerns the catastrophic brittle failure event that follows. Here a shear band of highly correlated plastic events nucleates and spreads quickly across the sample, leading to a discontinuous drop in the macroscopic stress. In Sec. IV, we discuss the implications of our results for experiments and molecular simulations. Sec. V gives conclusions and perspectives for future work. Before describing the model to be simulated, we make the following remarks to ensure that our terminology is unambiguous.
First, we use the term "brittle" to characterize abrupt yielding in which the stress drops suddenly from its maximum value that pertains just before yielding, to a final post-yield value, during a single event in which the strain becomes localized along a system-spanning shear band. We call yielding "ductile" if the fall in stress happens more gradually, or as a series of several smaller sudden drops, as a function of the applied strain. This terminology is consistent with that as used in Ref. [115]. We note, however, that the term brittle is sometimes reserved to describe yielding via crack propagation in the absence of significant plasticity, whereas significant plasticity does occur inside a forming shear band. In places where we use the term "brittle", therefore, readers concerned with this distinction may prefer to read it as "quasi-brittle".
Second, we use the term "shear banding" both in the sense typically adopted in the rheology literature (to mean a sustained coexistence of macroscopic bands flowing with different shear rates), and as used in the metallic glass literature (to mean the formation of a highly localised band of deformation), noting that these phenomena indeed also share many notable similarities. In any place where the distinction is not clear from the context, we shall clarify this explicitly. We note further that the formation of a shear band in studies with periodic boundary conditions (as the vast majority of theoretical works use) is likely in practice to be associated with the formation of a crack with new free surfaces, in an experimental system where the boundaries can move apart.
Third, we use the term "nucleation" to denote the onset of shear banding in metallic glasses for consistency with the use of the term in that literature, rather than necessarily to mean the thermally activated crossing of a barrier leading to the formation of a new phase or pattern. (These phenomena clearly nonetheless share notable similarities in their stochastic nature.)

II. ELASTOPLASTIC MODEL
Because amorphous materials lack any long range crystalline order, their behaviour cannot be understood in terms of (for example) the motion of defects in a background lattice structure. Instead, it has been attributed to the generic presence of structural disorder (think of a disordered arrangement of foam bubbles, for example), metastability (for foam bubbles to rearrange, the soap films must first stretch, which typically incurs an energy cost much greater than k B T ), and broken ergodicity (because of these energy barriers, the bubbles cannot rearrange in the absence of an externally applied shear).
A popular modelling approach is to distill these ingredients into a mesoscopic description that bridges the gap in lengthscales between a material's constituent microscopic substructures (foam bubbles, emulsion droplets, etc), and its emergent macroscopic flow behaviour. The basic idea is to consider any macroscopic sample of material to comprise many mesoscopic elastoplastic elements, each representing a patch of several (e.g.) foam bubbles. Each element is assumed to load elastically in shear, in between intermittent local plastic events. In any plastic event, a local energy barrier to particle rearrangements is surmounted, and stress is released.
In what follows we shall simulate the behaviour in shear of a 2D lattice elastoplastic model [2,116], which has a single elastoplastic element on each of N ×N lattice sites. Each element is assigned a local elastic shear strain l, a corresponding local shear stress kl, and a local yield energy E. In between local plastic yielding events, the strain of each element affinely follows the globally applied strain γ, corresponding to elastic loading withl =γ. For any element, this elastic loading continues until its strain energy reaches the yield threshold, 1 2 kl 2 = E. Beyond this threshold, the element yields stochastically at a rate τ 0 .
When any element yields, it adopts a new local strain chosen at random from a Gaussian distribution l w . The small value of l w is not important to the physics that follows and we set l w = 0.02 throughout. (As discussed below, we work in units in which the local yield strain l c ≡ 2E/k = 1 for all elements.) This yielding of an element at any lattice site models a local plastic yielding event within the material: physically, this might correspond to the sudden rearrangement of a few particles, say. Just after any such local yielding event, the stress across the 2D lattice will not be divergence free: i.e., force balance will be violated. Force balance is then immediately re-imposed via the Eshelby stress propagator in 2D [117,118]. In this way, the change in plastic stress experienced locally at any lattice site is propagated elastically to the surrounding medium. Specifically, the propagated elastic displacement field u resulting from a plastic event that results in a local stress change S is given at wavevector q in Fourier space by [117] The associated elastic stress follows as k 2 (∇u + ∇u T ). As is common practice in the literature on elastoplastic models, we consider only the shear component of this.
We consider throughout the limit of quasistatic shear, γτ 0 = 0, such that local yielding in fact occurs instantaneously once threshold is reached. Furthermore, in having assumed that the yielding rate of any element is zero until the yield threshold is reached, we are considering the athermal limit of zero temperature. Accordingly, the protocol studied here is that of "athermal quasistatic shear" (AQS), as considered widely in the literature on sheared amorphous materials. We assume inertia to be negligible throughout and, as is standard in studies of elastoplastic models, consider the stress propagation that follows any local yielding event to be well approximated by that in a homogeneous elastic medium.
To initialise the sample prior to shear, we assign each element a local shear strain l, chosen at random from a distribution that has zero mean and standard deviation l 0 . We then immediately impose force balance across the lattice. This imposition of force balance reduces the standard deviation of this initial distribution of local strains by a factor equal to the square root of the sum over all lattice sizes of the square of the Eshelby propagator. Throughout most of the paper, our numerical results will pertain to a Gaussian distribution of initial local strains. To demonstrate that our conclusions are robust with regards the shape of this distribution, however, we shall also perform calculations with a beta distribution of initial local strain values. To make a connection with some recent particle based simulations [119,120], we shall furthermore consider a variant of our model in which all elements have zero initial local strain and instead have different local yield stress values, distributed with Weibull statistics.
Small (resp. large) values of the parameter l 0 are taken to correspond to highly (resp. poorly) annealed samples. In a physical preparation protocol, a high degree of annealing could correspond (say) to equilibration at a low temperature (only just above the glass transition temperature), or to slow cooling from an initially high to zero temperature, or to ageing for a long time at a low temperature (below the glass transition temperature), before quenching to zero temperature. (Conversely, poor annealing would correspond to equilibration at a higher temperature before quench, faster cooling, or a shorter ageing time.) We do not model any of these annealing processes explicitly, but rather take the value of the single parameter l 0 as a proxy for the level of annealing, as is standard practice in studies of elastoplastic models.
So far, we have described the 'full', spatially aware elastoplastic model, with explicit force balance via an Eshelby propagator, in two spatial dimensions (2D). To simulate a mean field version of this model, we randomly shuffle the location of all elements on the lattice after each yielding event and the subsequent force rebalancing step. This is not the same as simply not implementing force balance, because the propagation of stress during force balance changes the distribution of local strains.
At any strain step, our numerical algorithm is as follows. First, we query how much elastic strain ∆γ must be applied in order to take the least stable element just above its yielding threshold. This amount of strain is then applied to every element on the lattice, l i → l i +∆γ, and the global strain variable incremented by the same amount, γ → γ + ∆γ. The least unstable element, now just above threshold, is then yielded. Force balance is then reimposed across the lattice. This may lead to other elements then exceeding their local thresholds. Any such elements are then yielded, and force balance is reimposed again. This process is repeated iteratively until no elements are left above yield after rebalancing. We then proceed to the next strain step. If, in any strain step, the amount of elastic strain needed to take the least stable element above threshold falls below ∆γ = 0.001, the strain increment is set to this minimal value.
Throughout we use units in which the elastic constant k = 1, and adopt the same yield energy threshold E for for all elements, chosen such that the yield strain of each element l c ≡ 2E/k = 1, as noted above. The physical parameters that remain to be explored are then l 0 , which sets the degree of initial sample annealing prior to shear (recall that small values of l 0 correspond to highly annealed samples), and the size of the sample of material being sheared, as set by the linear size N of the two dimensional N × N lattice.
The macroscopic stress is defined as the average across the lattice of the local elemental stresses:

III. RESULTS
The basic physics that we shall consider is summarised in Fig. 1, which shows a collection of state snapshots, each pertaining to a system of size N = 512 for the model with full 2D force balance. Each snapshot shows as a black dot the lattice sites that have undergone at least one plastic yielding event since shearing commenced. Unyielded sites are shown in white. Rows downwards correspond to progressively less well annealed samples. Columns left to right show the system's state after an increasing number of plastic events have occurred since shearing started.
At early strains after shearing starts, only a few plastic events occur, scattered throughout the sample in a largely uncorrelated way. These lead to a low level of strain heterogeneity slowly accumulating across the sample. We shall call this the "pre-failure" regime. This regime then terminates (for large samples or strong initial annealing) in a catastrophic brittle failure event, in which plastic events spread quickly across the sample in a highly correlated way, in the form of a shear band, and the macroscopic stress signal drops discontinuously.
A. Yielding is always brittle for large system sizes We shall start by confirming that the scenario of Refs. [107,108] also holds in the mean field version of the elastoplastic model simulated here, before then showing that it does not apply beyond mean field. See the stress-strain curves in Fig. 2a). These show a discontinuous stress drop for highly annealed samples (small values of l 0 ), but a continuous stress drop for poor annealing (larger l 0 ). (For the largest values of l 0 , the stress rises monotonically with strain.) As sketched in Fig. 2b) of Ref. [108], the discontinuous stress drop for strong annealing can be understood (in mean field) to stem from an underlying stress-strain curve that has an overhang, with the stress falling discontinuously from the curve's upper to lower branch once the overhang is reached.
The distinction between the regime of discontinuous stress drop for highly annealed samples and of continuous stress drop for poorly annealed samples is further explored in Fig. 3. We define the net stress drop during a shear simulation, for any given level of initial sample annealing, to be the difference between the maximum stress (maximised over all strains) and the (statistically) steady state stress approached as strain γ → ∞. We report this net stress drop as a function of the annealing parameter l 0 via the solid lines in Fig. 3a). We now further consider only the discontinuous part of this stress drop, defined as the maximum stress drop in any single strain increment. We report this quantity via the dashed lines in Fig. 3a). The difference between these two drops is plotted in Fig. 3c). As can be seen, a discontinuous stress drop exists for high levels of sample annealing, with l 0 < l * 0 ≈ 0.1. Its amplitude then falls continuously to zero as l 0 rises towards l * 0 . For poorly annealed samples (l 0 > l * 0 ), there is no discontinuous stress drop, at least for large systems. This behaviour furthermore shows only a rather weak dependence on system size, N .
So far, then, we have confirmed the scenario reported in earlier mean field studies of elastoplastic models [108], with brittle yielding for strong annealing (l 0 < l * 0 ), ductile yielding for poor annealing (l 0 > l * 0 ), and a "random critical point" separating these at l 0 = l * 0 . We now move beyond mean field and propose that the route to reconciling the apparent contradiction between Refs. [107,108] (brittle yielding for high annealing and ductile yielding for poor annealing), and Ref. [109] (brittle yielding for all levels of annealing, provided a stress drop indeed occurs), as discussed in Sec. I above, lies in a careful study of the size of the sample being sheared.
To demonstrate this, we shall simulate our 2D elastoplastic model of N × N elastoplastic elements with Eshelby force balance, for a previously unprecedented range of system sizes from N = 64 to 4096. For small N we shall recover the observation of Refs. [107,108], of brittle yielding for high annealing and ductile yielding for poor annealing. This distinction is however not governed by a critical point in our simulations. For large N we shall recover the scenario of Ref. [109], with brittle yielding for all levels of annealing, provided a stress drop indeed exists. In this way, we suggest that the scenario found in the earlier studies of 2D elastoplastic models [107,108], which apparently confirmed the random critical point found in mean field, may be a finite size effect.
The reconciliation just described is substantiated in Figs. 2 and 3, as follows. Fig. 2b) shows stress-strain curves computed within the 2D elastoplastic model with full force balance, for a large system size N × N = 4096 × 4096. Curves with decreasing maximum stress values correspond to decreasing levels of annealing. As can be seen, the stress drop is discontinuous not only for the highly annealed samples, but also (or very nearly so, for this N = 4096) even for poorly annealed samples. In contrast, simulating the same annealing levels for a small system size N = 256 in Fig. 2c) shows a convincingly discontinuous single stress drop only for strong annealing. In poorly annealed samples, the stress instead declines continuously, when averaged over many runs with different random number seeds. (In any single run, it decreases via a cumulative series of smaller discontinuous drops.) Fig. 2d) shows two collections of curves: the upper collection for strong annealing, and the lower collection for weak annealing, with each collection showing a full sequence of system sizes, N = 64, 128, 256 · · · 4096. For the highly annealed case, the stress drop is discontinuous even for small N . For the poorly annealed case, the stress drop is gradual and continuous for small N , but becomes increasingly steeper with increasing N . This is explored further in Fig. 3, in which panel b) shows as solid lines the net stress drop (from the maximum stress value to the steady state) as a function of the annealing parameter l 0 , averaged over many runs at each l 0 . The discontinuous part of the stress drop, also averaged over many runs at each l 0 , is shown as dashed lines. As can be seen, the discontinuous stress drop persists across the full range of l 0 in this non-mean field system, even for poorly annealed samples (high l 0 ). This is to be contrasted with the vanishing of the discontinuous stress drop at the critical point l 0 = l * 0 ≈ 0.1 in the mean field model. We therefore conclude that the critical point of the mean field model does not inform the yielding behaviour of the non-mean field model.
Subtracting the discontinuous part of the stress drop from the net stress drop for the non-mean field model in panel d), we see that the difference between the net and the discontinuous part of the stress drops shifts progressively to zero as the limit of infinite system size N → ∞ is approached. In the limit N → ∞, therefore, we expect the entire stress drop to occur discontinuously in a single strain increment, however poorly annealed the sample, corresponding to brittle yielding.
For very poorly annealed samples, with a broad initial distribution of local strains, the stress rises monotonically with strain, with no stress overshoot or drop. Indeed, in this case, significant plasticity occurs early in any simulation, because elements in the positive-l wing of the initial strain distribution exceed their local threshold already at small strains. This is true even for large N and clearly would be described as ductile rather than brittle behaviour. However, whether it would be described as yielding is a moot point, because yielding describes the progression with increasing strain from an initially elastic state to a finally plastic one.

Robustness to shape of initial local strain distribution
So far, the elements in our elastoplastic model have all had the same local yield strain, l c , but different initial local strain values drawn from a Gaussian distribution. From this it follows that the values of the initial stress needed to induce a local yielding event -the "stressto-yield" values -are also drawn from a Gaussian distribution. Recent atomistic simulations [119,120] have however indicated that the distribution of the stress-toyield values, approaches zero as a power law, in the limit of small stress-to-yield. In order to test the robustness of our predictions to the shape of the initial distribution of stress-to-yield values, we now consider two alternative initial conditions, aimed at capturing the statistics reported in Refs. [119,120].
The first of these keeps the yield stress values fixed across sites, at l c = 1, but replaces the Gaussian distribution of initial strain values with a β-distribution. This  has a probability density function given by for −1 ≤ l ≤ 1, with f (l) = 0 outside this region. Here Γ is the standard Gamma function. The β−distribution has zero mean, is symmetric about this mean, and decays (for α > 1) towards the yield strain l c = 1 as a power law with exponent α−1. This results in a power law decay of the distribution of stress-to-yield values, close to yield, as seen in the particle simulations [119,120]. After drawing the initial local strain values from this distribution we then impose force balance across the material before shear is applied. The other aspects of our simulation remain as described in Section 2.
As usual, we model the degree of annealing by the variance l 2 0 of the distribution of initial local strain values: a well-annealed sample has a small l 0 , and a poorly annealed sample a large l 0 . The exponent α of the β− distribution is related to its variance via: The condition α > 1, needed for the distribution to decay near threshold, holds as long as l 0 < 1/ √ 3 ≈ 0.5774. Stress-strain curves for different levels of annealing for this beta distribution, with l 0 values chosen to match the variances considered for a Gaussian distribution in Fig. 2(a,b), are shown in Fig. 4(a,b). Consistent with our observations in Fig. 2, we find a discontinuous stress drop for highly annealed samples for all sample sizes. For poorly annealed samples, we instead find a continuous stress drop for small and moderate system sizes, with the drop then tending to a discontinuous one in the limit of large system size. This is further explored in Fig. 4(c), which shows the total and (separately) the discontinuous part of the stress drop. The difference between these are shown in Fig. 4(d), and found to decay slowly as system size increases. In the limit of infinite system size, therefore, the entire stress drop is predicted to be discontinuous, corresponding to brittle yielding.
In our second alternative to a Gaussian distribution of initial local strains, we remove the randomness from the initial local strains, taking the initial strain of each element to be zero. We instead now adopt a distribution of initial local yield strain values, l c . For these we assume a Weibull distribution, which represents a generalisation of the exponential distribution, and which fits the form of the stress-to-yield curves obtained in Refs. [119,120], near threshold. The probability density function is This has two parameters: a scale factor λ and an exponent k ≥ 1, with the exponential distribution corresponding to k = 1. The scale parameter simply determines the units of stress. We choose to work in units such that the mean of the distribution is equal to 1, which is equivalent to choosing This leaves us with the single parameter, k, which sets the exponent of the power law decay l k−1 of the stressto-yield values, close to zero. It is related to the variance l 2 0 of the distribution by Accordingly large values of k correspond to small l 0 and so to well annealed samples. We consider values of l 0 so as to match the values of the variance used previously for the Gaussian and β distributions. At any l 0 , we invert Eqn. 7 numerically to find the corresponding value of k. To ensure that the final state of steady flow is independent of the initial condition, after any element's first local yielding event we choose its new yield strain from a Weibull distribution of unit mean and fixed variance l 2 w = 0.25. Otherwise our algorithm is as described in Section 2.
The resulting stress-strain curves are shown in Fig. 5(a,b). The net stress drop and discontinuous part of the stress drop are shown in Fig. 5(c,d). All the key features found previously for a Gaussian distribution are preserved, with brittle failure for highly annealed samples, and a slow progression from ductile to brittle failure with increasing system size for poorly annealed samples.
Having shown the same qualitative behaviour for three different shapes of the distribution of initial stress-toyield values, we suggest that our physical conclusions are robust with respect to that shape.
B. "Pre-failure" stage: slow accumulation of strain heterogeneity Returning to Fig. 1, we recall that yielding comprises two regimes. In the first regime, plastic events arise sporadically in a way that is (at least initially) relatively spatially uncorrelated. For highly annealed samples (upper rows in Fig. 1), these plastic events occur very infrequently; for poor annealing (lower rows), they are more frequent. In the second regime, a macroscopic shear band spreads quickly across the sample and the sample fails suddenly. In this section we aim to understand the first regime, before sudden failure occurs.
In this first regime, the occurrence of plastic events throughout the sample leads to a progressive accumulation of strain heterogeneity. This can be understood most easily within a simplified one-dimensional description in which we consider our N × N lattice in x and y to comprise an array of N streamlines arranged across the y direction, along each of which we project down the x coordinate by integrating over the N elements along x. (We could equally instead consider an array of N streamlines across x, integrating over the N elements along y.) It is worth noting that this equivalence between x and y holds on account of the π/2 rotational symmetry of the Eshelby propagator, in tandem with the fact that, in common with most numerical studies of elastoplastic models, we have ignored a term that would advect elements with the shear flow. Had we included advection, with a flow direction x and flow-gradient direction y, we anticipate that the strain heterogeneity discussed in this pre-failure regime would arise across the flow-gradient direction y, with the system's state remaining relatively invariant in x: any heterogeneity with wavevector along the x direction would be advected to finally have its wavevector indeed along y.
Consider then (in our model in which x and y are equivalent) a plastic event having just occurred on a given streamline at some y. This causes a loss of stress on that streamline of 1/N . (The streamline stress is defined as the average of the local elemental stresses along that streamline.) This streamline is therefore now out of force balance with the others: in a balanced system, the stress is the same across all streamlines. In the force balance step that immediately follows, therefore, the streamline on which the plastic event just occurred must strain forwards slightly to recover the same stress as all the other streamlines. Every streamline furthermore then strains backwards (by a lesser amount) to respect the fact that some stress has been lost from the system in the plastic event. The physics just described in words is what indeed results after imposing a full Eshelby 2D propagator, and then projecting down the result to one dimension.
Accordingly, after a given strain γ has been applied (on spatial average) across sample as a whole, any streamline of constant y (or of constant x) in Fig. 1 that has experienced a higher number of plastic events than the average number across the whole sample will have strained forwards a little more than the average γ. Those that have experienced statistically fewer events will have strained forwards a little less than γ.
Let us therefore now define the strain heterogeneity, δγ(γ), as the standard deviation of the strain across the y direction, having first integrated across x. This is plotted as a function of the average strain γ applied to the sample as a whole since shearing commenced in Fig. 6. Results are shown for the mean-field model in panel a) and for the 2D lattice model with force balance in b). The solid lines are the strain heterogeneity measured from our simulations. The dashed lines are the results of the analytical calculation that we shall now present in order to understand these simulation results.
Let us start by considering the distribution of local strains as first initialised, denotedP (l, t = 0). Recall that in most of our numerical results this is a Gaussian of width l 0 . (Robustness against switching to a beta or Weibull distribution having been demonstrated above. Indeed, the analytical arguments that follow below do not depend upon a particular choice for the shape of the distribution.) As described in Sec. II, we then immediately impose force balance. This narrows the width of the Gaussian by a factor equal to the square root of the sum across all lattice sites of the square of the 2D Eshelby propagator. We denote this re-normalised distribution by P (l, t = 0). It represents the state of the system immediately before shearing commences.
In what follows, consider the evolution of the probability distribution P (l, t) as a function of the time t or (almost) equivalently the accumulating strain γ =γt applied to the lattice as a whole. We also consider the associated macroscopic stress Σ = dl l P (l, γ) and the strain variable itself, γ. Our strategy will be to examine how these quantities evolve in a perfectly homogeneous "base state" in which the strain remains everywhere equal to γ, with no strain heterogeneity, and then to elucidate via a linear stability analysis the fate of small amplitude heterogeneous perturbations to this base state.
An evolution equation for the probability density function of local strains, P (l, t), can be written as follows (for the moment suppressing any space dependencies in our notation): ∂ t P (l, t) +γ∂ l P = −r(l)P + Y (t)P w (l).
Here r(l) is the yielding rate of an element with local strain l. Recall that r = 0 for |l| < 1 and r = 1 for |l| > 1. The average rate of plastic yielding across the lattice Y (t) = dl r(l) P . The second term on the LHS of Eqn. 8 accounts for the elastic loading of elements with strain rateγ. The first term on the RHS accounts for their "death" due to plastic yielding, and the second for their "rebirth" into new traps, with a post-yielding local strain selected from a distribution P w (l) of zero mean and small width l w . Pre-multiplying Eqn. 8 across by l and averaging over l then gives an evolution equation for the macroscopic stress Σ:Σ =γ − dl l r(l)P, In moving from the first line to the second, we have recognised that in the limit of quasistatic shear explored in this work,γ → 0, each element will yield instantaneously once it reaches its threshold l = 1. We now switch independent variables from the time t since shearing commenced at t = 0 to the strain γ =γt applied to the sample as a whole since t = 0, withγ the corresponding strain rate. We then have In the early stages of the deformation, before significant numbers of elements have yielded, the distribution P (l, γ) can be written to good approximation as P (l, γ) = P (l − γ, γ = 0) for l < 1, P (l, γ) = 0 for l > 1.
This encodes the fact that the stress distribution shifts affinely upwards with elastic loading, up to the yield threshold l = 1 at which elements "die", and remains relatively unperturbed by the "rebirth" of elements via plastic rearrangements until significant numbers of elements have yielded, in the later stages of deformation. Our equation of motion, now specialised to the early stages of quasistatic shear, can therefore be written: As written, Eqns. 8 to 12 appear to neglect the effect of force rebalancing. This can be included by considering a shear strain γ that is now not just the homogeneous strain averaged over the sample as a whole, but an effective heterogeneous strain that also accounts for any shears that arise internally during force balancing, and which lead also to a heterogeneous distribution P .
Let us now incorporate this in our simplified 1D approach. As noted above, our tactic is to consider an initially homogeneous "base state" in which the strain remains everywhere equal to the globally applied strain γ, with no heterogeneity; and then to study the fate of small heterogeneous perturbations to this base state. Ac-cordingly, we write: Σ(y, t) = Σ 0 (t) + k δΣ k (t) exp(iky), γ(y, t) = γ 0 + k δγ k (t) exp(iky), P (l, y, t) = P 0 (l, t) + k δP k (l, t) exp(iky), (13) in which (Σ 0 , γ 0 , P 0 ) denote the base state and (δΣ k , δγ k , δP k ) the small perturbations. Although we have decomposed the perturbations into Fourier modes, each of wavevector k, the governing equation of motion in fact contains no space-differential operators. (This is a property of the Eshelby propagator projected down to 1D [117].) Accordingly, each k−mode will behave in the same way and we now drop the subscript k.
Substituting the base state plus small perturbations into Eqn. 12, expanding in powers of the small perturbations and retaining only terms of first order then gives a linearised equation of motion for the perturbations: Force balance dictates that the stress must remain uniform across streamlines, and accordingly that δΣ = 0. Inserting this into Eqn. 14, and rearranging, gives: in which and (17) Eqn. 15 constitutes a linearised equation of motion governing the growth of the strain heterogeneity δγ as a function of strain γ applied to the sample as a whole since the inception of shear. Its solution is: Inserting into this the forms of f and g from Eqns. 16 and 17 gives, after some manipulation: The integral on the right hand side of Eqn. 19 represents the seeding of strain heterogeneities from the fact that the distribution of local strains, P , as initialised before shear commences and subsequently advected along by shear, is not quite the same across all streamlines, due to the finite number of elements N on each streamline. It is this phenomenon that is represented by the third term on the RHS of Eqn. 14. It arises from the disorder inherent to the amorphous material, as encoded here in the initial distribution of local strains: any streamline that happens to be initialised with a distribution of strains with a slightly higher average than that across all streamlines will suffer more plastic yielding and so strain more than the other streamlines as the level of shear applied externally to the sample as a whole increases.
To investigate this heterogeneous seeding further, we rewrite this integral in Eqn. 19 as We recognise this to be the standard deviation across streamlines of the cumulative fraction of elements with initial local strains between l = 1 − γ and l = 1, i.e., of elements that yield within the first γ strain units.
To compute this standard deviation, we note that among the N elements on any streamline, a fraction will have strains initially in that interval, and a fraction 1 − c will have strains initially outwith that interval. The probability of n elements having a strain in that interval is therefore given by the binomial distribution: which has normalised standard deviation: Finally, recognising that the cumulative fraction c of elements with initial strain in the interval 1−γ and 1 will be the fraction that yield in the first γ strain units, we can relate c to the macroscopic stress and strain variables, Σ 0 and γ 0 via c = γ 0 − Σ 0 . This follows from the fact that the stress Σ 0 (γ 0 ) after γ 0 strain units will be lower than the strain γ 0 by an amount equal to the fraction of elements that have plastically yielded. This can indeed be proved by a straightforward integration of Eqn. 10. The pre-factor to the integral in Eqn. 19 reflects that the innately seeded heterogeneity just discussed will further be amplified via the rebalancing of force in each strain step. Indeed, any streamline that suffers a little more plastic yielding will strain forward a little more in order to recover force balance, and so be susceptible to even further plastic yielding. Comparing this "amplification" pre-factor of Eqn. 19 with Eqn. 12, we see that it is simply equal to the inverse slope of the loading curve of stress as a function of strain, 1/(dΣ 0 /dγ 0 ).
Putting the above results together gives finally: We have now dropped the subscript 0, because the stress and strain in the homogeneous base state equal the macroscopically measured ones, to excellent approximation, in the early time regime in which the strain heterogeneity remains small and our linear calculation is valid.
We pause to emphasise the significance of this result. It predicts how strain heterogeneity develops within the sample as a function of the strain applied on average to the sample as a whole, in the early "pre-failure" phase of the yielding process. The factor 1/(dΣ/dγ) represents an amplification factor, due to the imposition of force balance at each strain step. This amplifies the innate statistical variation across streamlines in the initial seeding of elemental strains, due to the disorder inherent to the initial condition, which is described by the remaining terms. Importantly, Eqn. 24 has been recast in terms of the macroscopic global strain and stress, γ and Σ. In this way, the degree of strain heterogeneity within the sample is related to the externally measured rheological signals.
In the absence of force balance -i.e., in the mean-field elastoplastic model -the amplification factor described above will be missing, but the disorder innate to the initial condition remains. In mean field, therefore, the strain heterogeneity is predicted to grow instead simply as This is plotted by red dashed lines in Fig. 6a) and shows excellent agreement with our simulation results (black solid lines) for the growth of strain heterogeneity as a function of overall applied strain in the mean field model. Our prediction of Eqn. 24 for the 2D lattice elastoplastic model with force balance is plotted by red dashed lines in Fig. 6b). The corresponding simulation data is shown by black solid lines. Excellent agreement is again found between the two early in the yielding process. It is worth noting, however, that in this regime the difference between the mean field and full models is anyway relatively small, because the spatial correlations arising via force balance have not yet had chance to build. At larger strains, once the heterogeneity is better developed, the analytical prediction and simulation results quantitatively differ from each other. This is to be expected, for two reasons. First, once the heterogeneity becomes appreciable the linear assumption of the above analysis breaks down. Second, once significant numbers of plastic events occur, the base state will differ from the simple one assumed in our calculation, which simply takes the elastically shifted counterpart of the initial distribution. Nonetheless, the amplification predicted by the factor 1/(dΣ/dγ) arising from force balance is indeed present in the simulation data for the 2D lattice model (note the upturning black lines in the top right of Fig. 6b), and absent in mean field (note the nearly flatlining black lines in the top right of Fig. 6a). In the full 2D elastoplastic model, with x band denoting the coordinate (x or y) along which the band spreads, shifted such that the initial plastic event is at the origin. Right) In our simplified 1D elastoplastic model, also with N = 512. As can be seen, in each case the band devlopes via two "fronts" that spread out in opposite directions along x band (or its counterpart x1D).
For very poorly annealed samples, for which there is no stress overshoot as a function of strain, the amplification factor 1/(dΣ/dγ) pertaining to this pre-failure regime still predicts a divergence in strain heterogeneity, at the level of this linear calculation, as the stress Σ approaches as constant at large strains γ → ∞, such that dΣ/dγ → 0. Although nonlinear effects that become relevant once the amplitude of the strain heterogeneity becomes significant will then be expected to regularise any divergence, the (initially) leftmost curve of Fig. 6b) indeed confirms that appreciable strain heterogeneity arises even for these very poorly annealed samples.
We recall finally that the 1/ √ N scaling in Eqns. 24 and 25 stems from the fact that the distribution of local strains is not quite the same across all streamlines, due to the finite number of elements N on each streamline. In the limit N → ∞, this difference is predicted to approach zero, giving zero strain heterogeneity in this prefailure regime (before catastrophic failure then takes over, as explored below). This is a direct consequence of the assumption made almost ubiquitously in studies of elasto- b) Cumulative probability distribution of failure strains, extracted from data in a). Failure strain for each data set defined as strain at which δγ = 0.1. c) Average failure strain γ * as a function of the annealing parameter l0 (recall that small l0 corresponds to strong annealing), for system sizes N = 64, 128, 256, 512, 1024, 2048, 4096 in curves black, red, green, blue, violet, cyan, orange upwards at the right. γ * defined as strain at which P cumul = 0.5.
plastic models: that the origin of any initial heterogeneity indeed lies in the disorder inherent to the material. In deformation geometries with a heterogeneous stress field, as pertains for example to a soft material sheared in a curved rheometer, strain heterogeneity would additionally arise from that stress heterogeneity, and would be expected to exceed that arising from the disorder in the large N limit. Our calculation has assumed a regime in which heterogeneity arising from disorder exceeds that arising from systematically varying stress fields.
C. Catastrophic brittle failure stage

2D lattice elastoplastic model
Now we turn to a discussion of the catastrophic failure event, in which a shear band spreads quickly across the sample. Our particular aims are to predict and under- stand its onset in highly annealed samples The left panels of Fig. 7 show in more detail the initial development of the shear band along each of the top three rows in Fig. 1, for values of the annealing parameter l 0 = 0.025, 0.050, 0.075 respectively, corresponding to highly annealed samples. Each of these left panels in Fig. 7 shows the coordinate along the banding direction, x band , at which plastic events occur at successive iterations n of the model dynamics after the final strain step before banding occurs. (Recall that at each iteration we reimpose force balance, then yield any elements that are taken above yield in consequence.) The coordinate x band = y for l = 0.025 and 0.075 and x band = x for l = 0.050 in Fig. 1. x band is then shifted in Fig. 7 such that the plastic event that initiates the band lies at the origin. In the early stage of banding in highly annealed samples, most plastic events occur at the same value of the coordinate orthogonal to x band , with only a small number at the immediately adjacent coordinate value (not shown). This justifies the 1D treatment of banding that will follow. Fig. 8a) shows the strain heterogeneity δγ across the sample as a function of the spatially averaged strain γ, for a series of 128 runs, each with a different random number seed, again in a highly annealed system. In each run, δγ initially grows slowly in the "pre-failure" regime, then shows a sharp sudden upturn as a band nucleates. For each of the 128 runs, we define the strain at which failure occurs as that at which δγ crosses the given threshold, δγ = 0.1. We plot the cumulative distribution P cumul of these failure strains in panel b). We finally extract from it the average failure strain, γ * , defined as the strain at which the cumulative distribution equals 1/2, and plot γ * as a function of the level of the annealing parameter l 0 in c), for several different system sizes. The same data for the average failure strain, γ * are plotted in Fig. 9 as a function of the system size, for different values of l 0 .
For strongly annealed samples (small l 0 ), the failure strain decreases (noticeably) with decreasing levels of sample annealing (increasing l 0 ), and decreases (very weakly) with increasing sample size N . This can be understood as follows. For these highly annealed samples, the initial distribution of local strains is narrow. As the strain applied to the sample as a whole increases, this distribution shifts upwards. Eventually, the single element at its outmost rightward reach (largest l) fails. Once this happens, all other elements are by now near enough the yielding threshold themselves that the effect of the first element yielding will be to trigger a percolating cascade of further yielding events, resulting in a shear band. Therefore, the failure strain in this regime of strong annealing will simply be that at which the first element yields. This in turn is equal to one (the yielding threshold itself) minus the largest local strain initially present in the sample. This largest initial strain will increase with decreasing initial annealing, because the initial distribution is broader, and increase with increasing system size, because there are more elements in the sample to populate the extreme wings of the distribution. The arguments presented in words in this paragraph are encoded in the term P (l − γ, γ = 0) l=1 dγ in the calculation that follows in Sec. III C 3 below.
For weakly annealed samples, the failure strain increases with decreasing levels of sample annealing (increasing l 0 ) and increases with increasing sample size N . This can be understood as follows. For these poorly annealed samples, the the initial distribution of local strains is broad. As this distribution shifts upwards with increasing applied strain, elements start to yield. As the first few do so, however, an insufficient fraction of other elements are themselves close enough to threshold to provide a percolating onwards path along which further events can be triggered. Indeed, because of the breadth of the distribution in this regime, a large stress needs to accumulate in the sample as a whole before there is a percolating path for a shear band through the system. This becomes more pronounced the larger the distribution breadth (larger l 0 ) and the further the shear band has to propagate (larger N ).

Simplified 1D elastoplastic model
Our aims now are to understand the mechanism of shear band propagation in Figs. 1 and 7, to predict the probability distribution of applied strains γ at which banding sets in, as reported in Fig. 8b), and thereby finally to predict the average failure strain γ * in Fig. 8c). To this end, we shall first simulate a simplified elastoplas- tic model comprising a single streamline of elements arranged along the direction along which the band spreads.
For strong annealing, we will find this simplified 1D simulation to predict the banding dynamics of the 2D model actually rather well. We shall then perform a simple 1D analytical calculation, inspired by this 1D simulation, aimed at capturing the observed behaviour.
Our simplified 1D elastoplastic model comprises N elastoplastic elements arranged along a single line in the direction x 1D along which shear banding occurs. (x 1D is therefore the counterpart of x band in the 2D model.) The dynamics of the 1D model are chosen closely to mirror those of the 2D model, as follows. Prior to shear, each element is initialised with a local strain drawn from the same distribution as pertained immediately prior to shear in the 2D model. We then apply shear as follows. At each strain step, we first enquire which is the least stable element, and how much strain must be added to take it just over threshold. We then add that much strain to all elements, and to the global strain variable. That least stable element is now just above threshold, and is yielded. All elements then have their strain adjusted according to a 1D stress propagator that we discuss in the next paragraph, centred on the element that has just yielded. This process of stress propagation may then take some other elements above threshold. Those elements are then yielded, and the strain of all elements again adjusted via a superposition of our 1D stress propagators, each centred on one of the newly yielded elements. This process is repeated iteratively until no elements are left above threshold. We then proceed to the next strain step.
The 1D propagator to which we referred in the previ-ous paragraph is calculated from the full Eshelby propagator of the 2D lattice by taking the values of this function along one dimension only, at the origin of the other dimension. Therefore, if the stress propagator following an event at x = x 0 , y = y 0 on the 2D lattice is E(x−x 0 , y −y 0 ), the stress propagator following an event at y = y 0 on the 1D lattice isẼ(0, y − y 0 ) ≡ E(y − y 0 ). We observe numerically that the function E(y) is well fit to the form 0.309y −2 for y > 2 and y N , up to smaller variations on the spacing of the lattice site. For values of y approaching N , E departs from the power law decay to obey the periodic boundary conditions. E(1) = 0.079.
In the regime of strong annealing (small l 0 ), each run of this 1D model shows a failure event that closely resembles that seen in the full 2D simulation with matched l 0 , at least qualitatively. This can be seen by comparing each left panel of Fig. 7 for the 2D model with its counterpart in the right panel for the 1D model. For any l 0 , we run this 1D code N times, each with a different random number seed, to account for the fact that there are N streamlines in the 2D model. The minimum failure strain across these N separate runs is then taken as the prediction of the 1D model for the strain at which failure would occur in the 2D model. Running that N -fold simulation M times, we construct finally the cumulative probability distribution of failure strains P cumul (γ), across those M runs. From this we calculate finally the average failure strain γ * .
The cumulative distribution of failure strains is shown for three different levels of annealing in Fig. 10a). Solid lines show results for the 2D model, and dashed lines for the 1D model. For the most highly annealed samples, the correspondence between the 1D and 2D models is excellent. For less well annealed samples, we find progressively less good agreement. The average failure strain is shown as a function of the annealing parameter l 0 in Fig. 10b), by the solid and dashed line for the 2D and 1D models respectively. As can be seen, the average failure strain of the 1D model agrees reasonably with that for the 2D model even for samples as poorly annealed as l 0 ≈ 0.1, with full agreement in the strongly annealed limit, l 0 → 0.

Simplified 1D analytical calculation
We seek finally to understand the onset of shear banding via an analytical calculation that is inspired by our simplified 1D elastoplastic model, as follows.
In any strain increment γ → γ + dγ since the inception of shear, the probability of any given element locally yielding is P (l, γ) l=1 dγ. In the early stages of shear, recall that P (l, γ) = P (l − γ, γ = 0) to good approximation, because the initial distribution that pertained before shearing commenced has simply been shifted rightward by elastic loading, with only small changes due to plastic relaxation, which are neglected here. In the limit of small strain increment dγ → 0, the probability of ex-actly one out of N 2 elements yielding, P 1 (γ), is (The N 2 elements can be viewed either as the elements in N runs of our 1D code with N elements along a line; or as the N × N elements on the counterpart 2D lattice. The probability of two elements yielding is O(dγ 2 ) and can be neglected in the limit dγ → 0.) Let us now denote by r(γ) the probability that, when an element yields in the strain interval γ → γ + dγ, it leads to a runaway shear band and therefore sample failure. The probability that a shear band arises in γ → γ + dγ is then N 2 r(γ)P (l − γ, γ = 0) l=1 . Denoting by I the probability that no shear band has arisen by a strain of γ, we have This can be integrated to find I. The probability 1 − I that a shear band will have arisen by a strain γ, i.e., the cumulative probability distribution of failure strains, is then: (28) It remains finally for us to consider the probability r(γ) that, when a single element yields in the interval γ → γ + dγ, it leads to the development of a runaway shear band and therefore to sample failure.
Imagine a shear band that has already grown to comprise n contiguous events, centred on the plastic event that triggered it initially. What is the probability of an event occurring subsequently a distance m lattice sites away, on one side of this band? Such a site would, at the time the band first nucleated, have had a strain γ. It now has a strain due to the effects of stress propagation from each of the n sites that have already yielded along the band. Here E is the stress propagator used in the 1D elastoplastic model, as described previously. The probability of this site now yielding, p nm , is the probability that its initial strain before shearing commenced had been in the interval 1−γ to 1 − γ, which can be written as Here we use the same definition of the cumulative distribution as in Eqn. 20. The probability of zero elements on either side of the band of length n then yielding is in which the square accounts for the fact that there are two sides, one to the left of the pre-existing band, and one to the right. The probability that at least one event occurs and the band keeps spreading is then The probability of the band developing to a large number of events N events , giving a runaway instability, is then Inserting this expression 33 for the probability that a single plastic event occurring at strain γ triggers a runaway banding instability into Eqn. 28, and calculating this integral numerically, we find the cumulative distribution of failure strains given by the dotted lines in Fig. 10, left). (We have taken the limit of large N events and M , but in fact already found good convergence to this limit by N events = 16 and M = 2.) As can be seen, it shows excellent agreement with the the cumulative distribution of failure strains extracted from the simplified 1D elastoplastic model, which in turn shows good agreement with the corresponding distribution for the full 2D elastoplastic model across the range of l 0 values shown, with excellent agreement for the smallest l 0 . The correspondingly predicted averaged failure strain γ * is plotted as a function of the annealing parameter l 0 by the dotted line in Fig. 10, right), along with that from the 1D model simulation (dashed line) and the 2D model simulation (solid line).
While we believe the derivation up to Eqn. 28 to be rigorous, up to corrections that we have discussed, the subsequent arguments used to obtain r(γ) contain more severe assumptions. For example, Eqn. 29 assumes a proto-band of fully contiguous plastic events along a line, whereas in practice, successive events can skip over some lattice sites to leave them unyielded. (Recall Fig. 7.) Nonetheless, as discussed in the previous paragraph, this simple calculation gives convincing agreement with the results of our 1D simulations, at least for values of the annealing parameter up to l 0 = 0.075, and which in turn agree well with our full 2D simulations for the most strongly annealed samples.
Although the expressions 28 and 33 that determine the probability distribution of failure strains look rather complicated, the quantities upon which they depend, p nm and P (l − γ, γ = 0) l=1 , are together determined completely by the parameter l 0 that describes the level of disorder in the sample, as determined by the degree of sample annealing prior to shear, and the system size, N .
It is worth pausing to emphasise the significance of this result, which applies for materials described by the elastoplastic model studied here, in the limit of strong annealing and in quasistatic shear. It predicts the probability distribution of applied strain values at which catastrophic sample failure occurs, in terms of the disorder inherent in the sample prior to shear, as determined by the degree of sample annealing, and the system size N . It therefore indicates a possible way to predict catastrophic material failure before it occurs.

IV. DISCUSSION: IMPLICATIONS FOR EXPERIMENT AND MOLECULAR SIMULATIONS
In this section, we discuss the implications of our findings for experiment and molecular simulation.
Amorphous materials can be broadly subdivided into two (limiting, idealised) subcategories: "athermal" and "thermal". Athermal materials comprise substructures (foam bubbles, emulsion droplets, etc) large enough that Brownian effects can be discarded, with energy barriers impeding bubble/droplet rearrangement that greatly exceed k B T . Commonly studied experimental examples include dry granular packings, dense granular suspensions, foams, and emulsions. In thermal materials, on the other hand, the constituent substructures and/or local energy barriers are small enough that thermal fluctuations play an important role. Colloidal and polymeric glasses, colloidal gels, molecular and metallic glasses typically fall into this second class, although deep within the glass phase, far below the glass transition temperature, energy barriers will be large on the scale of k B T . The athermal model studied here is accordingly expected to apply primarily to materials with substructures large enough that Brownian motion can be neglected upfront, as well as molecular and metallic glasses deep within their glass phase.
Indeed, at any finite temperature T , the timescale for spontaneous particle rearrangement in the absence of shear, τ thermal , is expected to scale as τ 0 exp(E/k B T ), with E the typical energy barrier to rearrangement, and τ 0 the microscopic timescale for local intracage particle motion. In shear, groups of particles are expected to reach the threshold for local yielding on a timescale τ shear that scales as E/k/γ, i.e., as the time required to strain of order E/k at a strain rateγ. (Recall that k is the local modulus.) Once the threshold is reached, yielding is then expected to take place on the microscopic timescale τ 0 . The protocol of athermal quasistatic shear studied here pertains to the regime of well separated rates τ −1 thermal τ −1 shear τ −1 0 . In fully thermal materials, treating τ −1 thermal and τ −1 shear as well separated will not be valid. Indeed, as one moves way from the athermal limit and towards the thermal regime, activated processes will play an increasing role, avalanches associated with yielding are likely to start to overlap, and the yielding transition may be somewhat smoothed compared with that explored here. This would be an interesting topic for future study.
In the context of hard materials, the trend towards increasing fracture toughness with decreasing system size that we find for moderately to well annealed samples is consistent with experiments on nanoglass and metallic glasses [55], where it was suggested to arise from the reduced number of sites for the nucleation of shear bands as system size decreases. Increasing ductility with decreasing sample size was also reported experimentally in Ref. [121], with smaller samples being argued to require a larger stress for the nucleation of shear bands. A trend towards increasing ductility with decreasing sample size in samples subject to compressive loading was also discussed in [54], along with an important difference between failure in compressive and tensile loading, which remains an open challenge for future theoretical work. Decreasing fracture toughness with decreasing "fictive temperature" (increasing annealing) [122] has recently been reported in metallic glasses [123].
In the context of soft materials, our theoretical predictions are consistent with and indeed potentially explain the widespread experimental observation that these materials typically yield in a less abrupt way than hard materials. A key finding of our work has been that abrupt yielding tends to occur, for modestly annealed samples, in the limit of infinite system size (compared with the typical size of a material's constituent molecular or mesoscopic particles) and slow shear rate (compared with the inverse of the intrinsic material timescale). In comparison with hard materials, soft materials typically have much larger constituent mesostructures: the droplets in an emulsion are much larger than the constituent particles of a metallic or molecular glass glass, for example. Accordingly, their associated timescales are typically much slower. The crossover to the abrupt yielding that we predict for large system sizes and slow shear rates is accordingly likely to be accessed only at larger system sizes and slower shear rates for soft materials than hard materials, quite possibly outside the window explored in many existing experimental works.
Nonetheless, a growing body of experimental data increasingly shows that soft materials can also show brittle failure [38][39][40][41][42][43]. Such observations are emerging rapidly in the light of new experimental techniques that can now directly access the microscopic precursors to sudden failure [104][105][106]. Our results predicting the way a microscopic precursor -an initially localised plastic eventleads to a cooperatively growing shear band are particularly pertinent in the context of these new techniques. This connection with be explored in further detail in future work.
We have also made detailed quantitative predictions for the slow initial growth of strain heterogeneity within a material after shearing first commences, prior to any regime of final catastrophic material failure. Of particular relevance to experiment is our elucidation of the way in which this strain heterogeneity within a sample is predicted via the macroscopically measurable external rheological signals of strain and stress. This prediction is directly testable in flow velocimetry experiments of the kind that have been pioneered in studies of shear banding in soft materials such as emulsions and gels in recent years [124]. Specifically, such experiments are now capable of measuring the degree of strain rate heterogeneity across the sample as a function of time, and timeintegrating this signal to obtain the growing degree of strain heterogeneity. This can then be directly compared with our prediction in Fig. 6b and Eqn. 24 above by performing in-tandem measurements of the macroscopic stress-strain curve.
Our findings also have important implications for the interpretation of molecular simulations (in silico experiments) of low temperature glasses, for example as reported in Ref. [112]. In particular, the data in Fig. 3a of that work are consistent with the scenario reported here for moderately or poorly annealed samples, with ductile yielding for small system sizes and a trend towards brittle yielding for larger system sizes. Our work predicts that the same simulation curves reported for even more strongly annealed samples should show a crossover towards brittle yielding at smaller system sizes. This predition is directly testable, and is indeed consistent with some very recent molecular simulation results [125].
We have already made detailed discussion of our work in relation to the molecular simulations of Ref. [115] in Sec. III A above, and we shall not repeat this here.

V. CONCLUSIONS
In summary, we have studied theoretically the yielding of slowly sheared athermal amorphous materials. Our principal contributions have been to predict the conditions under which yielding will be ductile or brittle; to understand the initial growth of strain heterogeneity that is a precursor to material failure; to elucidate the way in which the nucleation and propagation of a shear band leads finally to catastrophic failure; and to predict the distribution of applied shear strains at which catastrophic failure will occur, in terms of the disorder inherent in the sample, as determined by the degree of sample annealing, and the system size.
For highly annealed samples, we have found yielding to be brittle for all samples sizes. In contrast, poorly annealed samples show an important dependence on the size of the sample of material being sheared, with apparently ductile yielding for small samples, and brittle yielding for large samples. This has important implications for experiment, in predicting a tendency towards increasing brittleness (for a fixed sample size) with increasing annealing; and that materials subject to a given level of annealing will show a different mode of failure, depending on the size of the sample being deformed.
Our work has shown that yielding comprises two distinct stages: a "pre-failure" stage, in which small levels of strain heterogeneity slowly accumulate within the material, followed by a catastrophic brittle failure event, in which a shear band quickly propagates across the sample via a cooperating line of plastic events.
In the pre-failure regime, we have provided an exact analytical expression for the slowly growing level of strain heterogeneity, expressed in terms of the macroscopically measurable stress-strain curve and the sample size, and in excellent agreement with our simulation results.
We have further elucidated the mechanism of subsequent catastrophic material failure, in which a shear band nucleates and spreads quickly across the sample. For highly annealed samples, our simulations have shown that that a single element first yields plastically and, in elastically propagating its stress to other elements via the Eshelby propagator, creates further nearby yielding events, leading quickly in turn to a percolating chain of events along a line [113].
In its use of periodic boundary conditions, our work necessarily pertains to the homogeneous nucleation of a shear band, triggered by plastic events that arise within the body of the disordered material, and which determines the ultimate strength of a material in the absence of surface imperfections. Clearly, it is imperative to understand this simpler case of homogeneous nucleation in order to elucidate the ultimate toughness of a material without surface imperfections; and also as a foundation on which to build an understanding of heterogeneous nucleation. In future, it would be interesting to simulate samples with external borders that have imperfections and indentations [122], to study the heterogeneous nucleation of shear bands arising in regions of concentrated stress. We have also limited ourselves to shear deformations, and it would be interesting in future to consider other deformation and loading protocols, such as planar extension.
It remains an open challenge to discriminate between or unify the picture of yielding elucidated here and those put forward in earlier studies, for example within a replica field theory [90][91][92]; as a directed percolation transition [64,65]; within a random first order transition theory for the glass transition [94]; as a Gardner transition [95]; as a spinodal point [96]; and within particle simulations that seed initial weak spots [97].