The role of fluctuations on the yielding transition of two-dimensional glasses

We numerically study yielding in two-dimensional glasses which are generated with a very wide range of stabilities by swap Monte-Carlo simulations and then slowly deformed at zero temperature. We provide strong numerical evidence that stable glasses yield via a nonequilibrium discontinuous transition in the thermodynamic limit. A critical point separates this brittle yielding from the ductile one observed in less stable glasses. We find that two-dimensional glasses yield similarly to their three-dimensional counterparts but display larger sample-to-sample disorder-induced fluctuations, stronger finite-size effects, and rougher spatial wandering of the observed shear bands. These findings strongly constrain effective theories of yielding.

We numerically study yielding in two-dimensional glasses which are generated with a very wide range of stabilities by swap Monte-Carlo simulations and then slowly deformed at zero temperature. We provide strong numerical evidence that stable glasses yield via a nonequilibrium discontinuous transition in the thermodynamic limit. A critical point separates this brittle yielding from the ductile one observed in less stable glasses. We find that two-dimensional glasses yield similarly to their three-dimensional counterparts but display larger sample-to-sample disorder-induced fluctuations, stronger finite-size effects, and rougher spatial wandering of the observed shear bands. These findings strongly constrain effective theories of yielding.
Amorphous solids encompass a wide variety of systems ranging from molecular and metallic glasses to granular media, also including foams, pastes, emulsions, and colloidal glasses. Their mechanical response to a slowly applied deformation exhibits features such as localized plastic rearrangements, avalanche-type of motion, the emergence of strain localization and shear bands [1][2][3][4][5]. The universality of these phenomena suggests that a unified description may be possible. When slowly deformed at low temperature from an initial quiescent state, amorphous solids yield beyond some finite level of applied strain and reach a steady state characterized by plastic flow. Understanding yielding is a central issue in materials science, where one would like to avoid the unwanted sudden failure of deformed glass samples [6]. It is also a challenging problem in nonequilibrium statistical physics [1].
The ways in which amorphous materials yield can be classified in two main categories: the "brittle" extreme where the sample catastrophically breaks into pieces and show macroscopic shear bands [6], as often observed in molecular and metallic glasses, and the "ductile" behavior in which plastic deformation increases progressively [4], commonly found in soft-matter glassy systems.
The key question is whether the observed variety of yielding behaviors should be described by (i) completely distinct approaches, with, e.g., ductile yielding being describable by soft glassy rheology models [7] while brittle yielding falls into the realm of the theory of fracture [8-10]; (ii) as a unique phenomenon, taken as the ubiquitous limit of stability of a strained solid in the form of a critical spinodal [11][12][13][14], or, (iii) as we have recently argued on the basis of mean-field elasto-plastic models and simulations of a three-dimensional atomic glass [15], within a unique theoretical framework but with the nature of yielding depending on the degree of effective disorder that is controlled by the preparation of the amorphous solid. In the latter case, and in analogy with an athermally driven random-field Ising model (RFIM) [16], yielding evolves from a mere crossover (for poorly annealed samples) to a nonequilibrium discontinuous transition past a spinodal point (for well-annealed, very stable samples); the transition between these two regimes is marked by a critical point that takes place for a specific glass preparation [17] (here we focus on a uniform shear deformation but a similar scenario would apply to an athermal quasi-static oscillatory shear protocol [18,19]). Strictly speaking, these sharp transitions can only be observed at zero temperature in strain-controlled quasi-static protocols [20]. However, temperature is likely to play a minor role in realistic situations given the large energy scales at play. In fact, there is experimental evidence that a given material may indeed show brittle or ductile yielding depending on preparation history of the sample [21-25].
If yielding is a bona fide (albeit nonequilibrium) discontinuous phase transition ending in a critical point akin to that of a RFIM, one should wonder about its universality class and its dependence on space dimension. By default or with the assumption that the phenomenology is qualitatively unchanged when changing dimension, many of the numerical studies of yielding in model amorphous solids have been carried out in two dimensions (2D) in various approaches [26][27][28]. Whereas this may be legitimate when focusing on the flowing steady state or on very ductile behavior, one should be more cautious about the role of spatial fluctuations on the nature, and even the existence, of the yielding transition itself as one decreases the dimension of space. Fluctuations of the order parameter are expected to change the values of the exponents of the critical point as dimension is decreased below an upper critical dimension at which the mean-field description becomes qualitatively valid. More importantly, they smear the transition below a lower critical dimension. In the standard RFIM with ferromagnetic short-range interactions and short-range correlated random fields, this lower critical dimension has been proven to be D = 2 arXiv:1912.06021v1 [cond-mat.dis-nn] 12 Dec 2019 for the equilibrium behavior [29] and, although still debated [30,31], appears to also be D = 2 for the out-ofequilibrium situation of the quasi-statically driven RFIM at zero temperature. However, we expect that the relevant RFIM providing an effective theory for yielding is not the standard one. It is indeed known that elastic interactions in an amorphous solids are long-ranged and anisotropic [32, 33] instead of short-ranged ferromagnetic, as indeed shown by the appearance of strong anisotropic strain localization in the form of shear bands.
Therefore, a careful study of yielding in 2D model atomic glasses as a function of preparation is both of fundamental interest and relevant to two-dimensional physical materials, such as dry foams [23], grains [34], or silica glasses [35]. This is what we report in this letter, where we consider glass samples that are prepared by optimized swap Monte-Carlo simulations [36,37] in a wide range of stabilities from poorly annealed glasses to very stable glasses and that are sheared through an athermal quasi-static protocol. We provide strong evidence that strained 2D stable glasses yield through a sharp discontinuous stress drop, which from finite-size scaling analysis survives in the thermodynamic limit, as in 3D. As the stability of the glass decreases, brittleness decreases and below a critical point, which is characterized by a diverging susceptibility, yielding becomes smooth. Compared to the 3D case, we find that the 2D systems are subject to larger sample-to-sample fluctuations, stronger finite-size effects, and rougher spatial wandering of the shear bands.
We focus on a glass model formed by size-polydisperse particles interacting with a soft repulsive potential [37]. Glass samples are prepared by first equilibrating liquid configurations at a finite temperature, T ini , and then performing a rapid quench to T = 0. The samples are then subsequently deformed at zero temperature. The preparation temperature T ini (sometimes referred to as the "fictive temperature") then uniquely controls the glass stability; we consider a wide range of preparation temperatures, T ini = 0.035 − 0.200 (see the Supplementary Information (SI) [38]). To obtain wellannealed systems, we use the swap Monte Carlo algorithm, which allows equilibration at extremely low temperatures [36]. The considered range of T ini covers poorly annealed glasses (T ini = 0.150 − 0.200), ordinary computer glasses (T ini = 0.100 − 0.120), as well as ultrastable glasses (T ini = 0.035 − 0.070). We then perform a straincontrolled athermal quasi-static shear deformation [39], using Lees-Edwards boundary conditions. To carefully analyze finite-size effects, we vary the number of particles N over a very large range, N = 1000 − 128000, and for a proper averaging we consider several hundreds of samples for each size and T ini . More details about model and methods are given in the SI [38]. ferent types of behavior depending on the initial stability: monotonic crossover for poorly annealed samples (T ini = 0.150 − 0.200), mild stress overshoot for ordinary computer glass samples (T ini = 0.100 − 0.120), and a sharp discontinuous stress drop for very stable samples (T ini = 0.035 − 0.070). This plot is qualitatively similar to that found in 3D [15]: ductile yielding is observed for higher T ini and appears to continuously transform into brittle yielding below T ini ≈ 0.10 − 0.12.
We first give evidence through finite-size scaling analysis that brittle yielding persists in 2D as a nonequilibrium first-order (or discontinuous) transition. For the most stable glass considered (T ini = 0.035), we show in Figs. 1(b) and (c) the stress-strain curves after averaging over many samples and the so-called "disconnected" susceptibility [29], χ dis = N ( σ 2 − σ 2 ) ( · · · denotes average over samples). As N is increased, the slope of the drop following the stress overshoot becomes steeper, suggesting that the averaged stress-strain curve shows a discontinuous jump as N → ∞. As shown below, this is due to the sudden appearance of shear bands [40]. Concomitantly, the disconnected susceptibility χ dis grows with N . We plot its peak values as well as that of the socalled "connected" susceptibility [29], χ con = −d σ /dγ, in Fig. 1(d). We find that both susceptibilities increase with N , χ peak dis ∝ N and χ peak con ∝ N 0.3 , which is a signature of a nonequilibrium first-order (i.e., discontinuous) transition. The stronger divergence of χ peak dis indicates the predominant role of disorder fluctuations, as generically found in the RFIM. The observed finite-size scaling of χ peak dis is the same in 2D and 3D, and reflects the discontinuous nature of the transition where sample-to-sample stress fluctuations at a fixed yield strain are of O(1). On the other hand, the scaling of χ peak con is different from 3D, for which we found χ peak con ∝ N 0.5 [15]. The dominant effect explaining this difference comes from the scaling of the width of the distribution of γ Y at which the largest stress drop takes place. Sample-to-sample fluctuations seem to lead to a standard N −0.5 behavior in 3D but to a broader distribution with a width decaying only as N −0.3 in 2D (see the SI). Within a RFIM perspective [41], this entails that the variance of the effective random field at the transition scales with the linear system size L = N 1/D as ∆ L = χ dis /χ 2 con ∼ L ρ with ρ ≈ 0.8 in 2D and ρ ≈ 0 in 3D. This in turn implies that the random field at yielding has long-range correlations decaying with distance as r −(D−ρ) with ρ > 0 [42, 43] in 2D, a feature that seems absent in 3D yielding.
To further illustrate the strong sample-to-sample fluctuations present in 2D, even in the case of very stable glasses, we show in Fig. 2(a) a zoomed-in plot of the stress-strain curves for a few chosen samples. One can see that in addition to typical samples that display a single sharp, large stress drop (see also Fig. 1(a)) there are samples that yield in multiple stress drops. These samples, which we refer to as "anomalous", display shear bands at yielding that tend to strongly wander and splinter in space (see Fig. 2(d)), whereas typical samples yield via the appearance of a well-defined system-spanning shear band (see Fig. 2(c)) [44]. Anomalous samples lead to very large fluctuations but their fraction, f N , decreases as N increases. To quantify this effect, we have identified these samples from individual stress-strain curves, using the conditions ∆σ max ≤ 0.15, where ∆σ max is the maximum stress drop observed in the strain window γ ∈ [0, γ max ] for each sample (see SI for details [38]). In Fig. 2 we observe that f N decreases with N , in an apparently exponential manner, and appears to vanish as N → ∞ (hence the terminology "anomalous" versus "typical"). Repeating the same analysis for 3D glass samples of a similar stability, we find that f N virtually vanishes above N = 12000 (not shown), which indicates much weaker finite-size effects than in 2D.
Next, we investigate the spatial characteristics of the shear bands. We carefully analyze the wandering of the shear bands in space, comparing 2D and 3D at the same lengthscale and for similar glass stability [45] (excluding anomalous samples). In Figs. 3(a) and (b), we show snapshots right after yielding for 3D (γ = 0.13) and 2D (γ = 0.07), respectively. The 2D shear band appears thicker than the 3D one. More importantly, it seems to wander more or, said otherwise, to be "rougher". This roughness can be quantified by the height-height correlation function [8, [46][47][48]: where Y com (x) is the average height of the shear band and · · · x denotes a spatial average. For the 3D case we use an analogous expression which also takes into account the average in the additional z coordinate (see SI for a detailed explanation). A manifold is rough on large scales if the height-height correlation function scales with distance as C(∆x) ∝ ∆x ζ , where ζ > 0 is the roughness exponent. Therefore, the log-log plot of C(∆x) versus the distance ∆x provides a way to assess the roughness of the shear bands. We show such a plot in Fig. 3(c). The data in 3D show no convincing effect over the (limited) covered range but the results in 2D point to a nontrivial intermediate regime (limited on the longest lengthscales by a saturation due to the system size) where an effective roughness exponent ζ ≈ 0.59 can be seen. It is also clear that the overall magnitude of C(∆x) in 2D is much larger than the one in 3D, and this difference is expected to grow even larger in larger sampels. This again reflects the presence of larger spatial fluctuations in 2D. Because we have found strong evidence that yielding in 2D is a genuine nonequilibrium discontinuous transition for very stable glasses and that poorly annealed samples clearly show a continuous ductile behavior (see Fig. 1(a)), it is tempting to look for signatures of a critical point separating brittle and ductile yielding as one varies the preparation temperature T ini of the glass samples. In Ref.
[15], we showed that the difference between the stress before (σ 1 ) and after (σ 2 ) the largest stress drop, ∆σ max = σ 1 − σ 2 , plays the role of the order parameter distinguishing brittle from ductile yielding in 3D. In particular, we demonstrated that the critical point at T ini,c can be identified by the divergence of the variance of this order parameter, N ( ∆σ 2 max − ∆σ max 2 ). In 2D however, as seen in Fig. 1(a), strong fluctuations seem to also affect the plastic steady state, irrespective of the presence of a critical point. Even for typical samples, this blurs the determination of the largest stress drop when the latter becomes small as one approaches the putative critical point. As an operational procedure to remove this effect, we have therefore defined the order parameter as ∆σ max ≡ σ 1 − σ 2 . In this way the fluctuations of σ 2 that we tentatively attribute to the plastic steadystate regime are explicitly removed. When applied to 3D this new definition captures the critical point even more sharply without affecting the results.
The mean value ∆σ max is shown in Fig. 4(a) (we have removed a trivial offset at high T ini which vanishes in the large-N limit, see the inset); ∆σ max appears rather flat at high T ini and starts to grow and to develop a significant dependence on system size below T ini ≈ 0.1, showing a similar trend as the 3D case. The variance of ∆σ max is displayed in Fig. 4(b). It shows a peak that grows and shifts toward higher T ini with N . The data is not sufficient to allow for a proper determination of critical exponents but give nonetheless support to the existence of a brittle-to-ductile critical point around T ini,c ≈ 0.1.
In summary, we have given strong numerical evidence that 2D yielding of very stable glasses under athermal quasi-static shear remains a nonequilibrium first-order (discontinuous) transition that survives in the thermodynamic limit, with a dominance of the disorder-induced, i.e., sample-to-sample, fluctuations. Furthermore, the transition to ductile yielding is signalled by a critical point. The scenario found in 2D is therefore analogous to the one found in 3D, but with stronger fluctuation effects. On the one hand, this suggests that the brittle-to-ductile transition as a function of sample preparation could be also experimentally observed in 2D or quasi-2D amorphous materials. On the other hand, it confirms that if indeed the effective theory describing the yielding transition is an athermally quasi-statically driven RFIM, the basic features of the model are necessarily modified by the presence of long-range anisotropic Eshelby-like interactions and, in 2D possibly, by long-range correlations in the effective random field. The long-range and quadrupolar nature of the elastic interactions accounts for the appearance of a shear band at the spinodal point marking the start of brittle yielding in stable glasses. In this modified RFIM, the nature of the spinodal and its potential critical character [14] still needs to be investigated.
We acknowledge support from the Simons Foundation (#454933, L. Berthier, #454935, G. Biroli).     [49]. We have checked that the properties of horizontal and vertical shear bands are virtually the same.

Details on the simulation methods
Simulation model. The glass-forming model consists of particles with purely repulsive interactions and a continuous size polydispersity. Particle diameters, d i , are randomly drawn from a distribution of the form: where v 0 sets the unit of energy (and of temperature with the Boltzmann constant k B ≡ 1) and = 0.2 quantifies the degree of nonadditivity of particle diameters. We introduce > 0 in the model to suppress fractionation and thus enhance the glass-forming ability. The constants c 0 , c 1 and c 2 enforce a vanishing potential and continuity of its first-and second-order derivatives at the cut-off distance r cut = 1.25d ij . We simulate a system with N particles within a square cell of area V = L 2 , where L is the linear box length, under periodic boundary conditions, at a number density ρ = N/V = 1. We also compare the results with the corresponding results of the 3D system studied in Ref. [15].
Glass preparation. Glass samples have been prepared by first equilibrating liquid configurations at a finite temperature, T ini (which is sometimes referred to as the fictive temperature of the glass sample) and then performing a rapid quench to T = 0, the temperature at which the samples are subsequently deformed. We prepare equilibrium configurations for the polydisperse disks using swap Monte-Carlo simulations [36]. With probability P swap = 0.2, we perform a swap move where we pick two particles at random and attempt to exchange their diameters, and with probability 1 − P swap = 0.8, we perform conventional Monte-Carlo translational moves. To perform the quench from the obtained equilibrium configurations at T ini down to zero temperature, we use the conjugate-gradient method [50].
The preparation temperature T ini then uniquely controls the stability of glass. We consider a wide range of preparation temperatures, from T ini = 0.035 to T ini = 0.200. To better characterize this temperature span, we give some empirically determined representative temperatures of the model: Onset of slow dynamics [51] takes place at T onset ≈ 0.23, the dynamical mode-coupling crossover [52] at T mct ≈ 0.11, and the estimated experimental glass transition temperature, obtained from extrapolation of the relaxation time [36], at T g ≈ 0.068. Note that these values are slightly different from the ones presented in Ref. [37] due to a small difference in the number density ρ. Our range of fictive temperature T ini therefore covers from slightly below the onset temperature to significantly below the estimated experimental glass-transition temperature.

Mechanical loading.
We have performed straincontrolled athermal quasi-static shear (AQS) deformation using Lees-Edwards boundary conditions [39]. The AQS shear method consists of a succession of tiny uniform shear deformation with ∆γ = 10 −4 , followed by energy minimization via the conjugate-gradient method. The AQS deformation is performed along the x-direction up to the maximum strain γ max = 0.2. Note that during the AQS deformation, the system is always located in a potential energy minimum (except of course during the transient conjugate-gradient minimization), i.e., it stays at T = 0.
Nonaffine displacement. We consider the local nonaffine displacement of a given particle relative to its nearest neighbor particles, D 2 min [53]. D 2 min is always measured between the origin (γ = 0) and a given strain γ. We define nearest neighbors by using the cut-off radius of the interaction range, R cut = 3.0d. We determine the nearest neighbors of a particle from the configuration at γ = 0.

Supporting data
Computation of the susceptibilities. To numerically compute the susceptibilities, we perform a smoothing procedure by averaging over 10 adjacent data points, as described in Ref. [15]. Figure 5 shows the parametric plot of the logarithms of the peak values of χ dis and χ con . The data points roughly follow in a straight line in the putative brittle yielding phase (T ini 0.10), with a slope of about 3.
Finite-size scaling of the susceptibilities for a discontinuous transition in the presence of disorder. We consider the case of a very stable glass (appropriate, e.g., for T ini = 0.035), for which each typical finite-size sam- ple yields through a discontinuous stress drop, ∆σ max , of O(1) at a strain value γ Y . Both the size of the stress drop and the yield strain are sample dependent. It is easily realized that provided the mean value of ∆σ max is strictly positive and of O(1), the fluctuations of ∆σ max lead only to subdominant contributions to the finite-size scaling of the susceptibilities, at least well below T ini,c . The fluctuations of the yield strain on the other hand are crucial. They are likely to regress with the system size N but possibly in a nontrivial fashion. We assume that the values of γ Y are distributed around the peak position, γ * Y , of the distribution according to some probability function scaling with N as with δ > 0 and P(0) > 0. We have computed this distribution of γ Y for 2D samples prepared at T ini = 0.035 and the result is shown in Fig. 6(a). We then perform exercises of scaling collapse assuming Eq. (4) with various δ in Figs. 6(b,c). We find that the conventional value, δ = 0.5 [54], does not work, while a smaller value, δ = 0.3, provides a good scaling collapse. In the following, we will relate the obtained value, δ = 0.3, with the scaling of the susceptibilities.
In the vicinity of yielding, one may describe the stress in each sample as given by where θ(x) is the Heaviside step function and σ 2 is the stress right after the stress drop. As discussed in the main text, this value may also fluctuate, but just as for the fluctuations of ∆σ max this leads to only subdominant corrections to the leading finite-size scaling in the regime where a strong discontinuous transition is present. We therefore assume from now on that neither ∆σ max nor σ 2 fluctuate from sample to sample. One then easily derive that the first two cumulants of σ(γ) close to the yielding transition are expressed as and (7) From the above expressions it is easy to derive the connected susceptibility, and the disconnected one, (9) By taking now into account the scaling form of the distribution P N (γ Y ) in Eq. (4), one immediately obtains that the maximum of the susceptibilities scales as Putting δ = 0.3 obtained numerically by the scaling analysis in Fig. 6, we show a comparison between the above predictions and the directly determined N dependence of χ peak dis and χ peak dis in Fig. 1(d) of the main text. We find excellent agreements.
Identification of the anomalous samples. Figure 7 shows the probability distribution function of the maximum stress drop ∆σ max for 2D (a,b,c) and 3D (d) for low T ini 's. For the most stable case in 2D, T ini = 0.035, there are two peaks for the smaller values of N . The right peak correspond to samples with a single large discontinuous stress drop, which we call typical samples, and the left peak to samples with multiple stress drops at yielding, which we call anomalous samples. For N = 8000 the peak at smaller ∆σ max is dominant, which means that the majority of samples show multiple stress drops. However, the peak at higher ∆σ max grows with increasing N , and for large enough system size, most of the samples show a single large discontinuous stress drop at yielding (hence the denomination of typical samples); yet there remains a tail at smaller ∆σ max corresponding to the anomalous samples. We observe the same trend up to T ini = 0.050, but the peak positions shift toward smaller ∆σ max with increasing T ini . Above this T ini = 0.050, we do not find any hint of two separate peaks, which forbids any sensible distinction of typical and anomalous samples. To separate anomalous samples from typical samples for T ini = 0.035 − 0.050, we choose a cutoff ∆σ max = 0.15 which seems to reasonably distinguish anomalous samples (∆σ max ≤ 0.15) from typical ones (∆σ max > 0.15): see Fig. 7 (a,b,c). There is some leeway in defining this cutoff value, but the conclusion in the main text does not change if we slightly change the value. In contrast to 2D, 3D systems do not show a clear bimodal distribution, as seen in Fig. 7 (d). Besides, the tail at smaller ∆σ max is significantly suppressed compared to 2D. To nonetheless make an attempt to quantify the fraction of anomalous samples, we have chosen a cutoff at ∆σ max = 0.2.
Roughness analysis. We present some details on how we have computed the height-height correlation function C(∆x) for the shear bands from the molecular simulation data.
Particles are considered as part of the shear band if their nonaffine (squared) displacement D 2 min is large enough. Figure 8 shows the probability distribution of D 2 min , P (D 2 min ), for several values of the strain γ covering the regimes before and after yielding. Before yielding (γ = 0.03 − 0.05 for 2D and γ = 0.09 − 0.11 for 3D), P (D 2 min ) is localized near the origin, which reflects the fact that most of the particles show a purely affine deformation and that only a very small fraction of particles undergo nonaffine displacements in the pre-yield regime. After yielding (γ = 0.07−0.09 for 2D and γ = 0.13−0.15 for 3D) on the other hand, a significant tail suddenly appears due to strain localization in the form of shear bands, and this tail grows with increasing γ. By introducing the thresholds shown in the vertical arrows in Fig. 8, we separate particles belonging to the shear band (with a nonaffine displacement above threshold, D 2 min > 0.4 for 2D and D 2 min > 0.59 for 3D) from particles undergoing affine displacement characteristic of a purely elastic solid. An illustration is given in Fig. 9(a) for a 2D sample.
We compute the average location of the shear band (line in 2D or surface in 3D) by discretizing the base space (x for a horizontal shear band in 2D, and (x,z) for a horizontal band in 3D). More specifically, we compute the y-coordinate of the center of mass, Y com (x), of the particles belonging to the shear band and located within the bin specified by the position x (or x and z in 3D).
The output is illustrated in Fig. 9(b). For the bin width we have used 2.53 in 2D and 2.28 in 3D.
The height-height correlation functions C(∆x) is finally defined as C(∆x) = (Y com (x + ∆x) − Y com (x)) 2 1/2 x (11) in 2D. In 3D where the z-axis has to be taken into account, the averaging procedure for C(∆x) is also performed along the z-direction, according to C(∆x) = (Y com (x + ∆x, z) − Y com (x, z)) 2 1/2 x,z .