# Emergent scales and spatial correlations at the yielding transition of glassy materials

Stefano Aime, Domenico Truzzolillo

► To cite this version:

Stefano Aime, Domenico Truzzolillo. Emergent scales and spatial correlations at the yielding transition of glassy materials. Physical Review Materials, 2025, GDR_MFA, 9 (5), pp.055602. 10.1103/Phys-RevMaterials.9.055602. hal-05080626

HAL Id: hal-05080626

https://hal.science/hal-05080626v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Emergent scales and spatial correlations at the yielding transition of glassy materials

Stefano Aime

Laboratoire Matière Molle et Chimie, (UMR 7167, ESPCI-CNRS) ESPCI ParisTech, 10 rue Vauquelin, 75005 Paris, France

Domenico Truzzolillo

Laboratoire Charles Coulomb (L2C), UMR 5221, CNRS-Université de Montpellier, France

(*domenico.truzzolillo@umontpellier.fr)

(*stefano.aime@espci.fr)

(Dated: May 22, 2025)

Glassy materials yield under large external mechanical solicitations. Under oscillatory shear, yielding shows a well-known rheological fingerprint, common to samples with widely different microstructures. At the microscale, this corresponds to a transition between slow, solid-like dynamics and faster liquid-like dynamics, which can coexist at yielding in a finite range of strain amplitudes. Here, we capture this phenomenology in a lattice model with two main parameters: glassiness and disorder, describing the average coupling between adjacent lattice sites, and their variance, respectively. In absence of disorder, our model yields a law of correspondent states equivalent to trajectories on a cusp catastrophe manifold, a well-known class of problems including equilibrium liquid-vapour phase transitions. Introducing a finite disorder in our model entails a qualitative change, to a continuous and rounded transition, whose extent is controlled by the magnitude of the disorder. We show that a spatial correlation length $\xi$ emerges spontaneously from the coupling between disorder and bifurcating dynamics. With vanishing disorder, $\xi$ diverges and yielding becomes discontinuous, suggesting that the abruptness of yielding can be rationalized in terms of a lengthscale of dynamic heterogeneities.$^a$

# I. INTRODUCTION

Understanding the flow of glassy materials is a scientific challenge of great importance for many practical purposes that has been animating the research activity of a wide scientific community. Particular effort has been devoted in the last decades to understanding the onset of flow, from the solid-like response characterizing small deformations to progressive fluidization, with a response dominated by plastic dissipation as the deformation is augmented beyond the linear regime. This transition is called yielding, and exhibits similar rheological features in samples of interest for several research fields and applications, ranging from food processing and cosmetics to oil extraction, concrete formation and pharmaceutics [1, 2]. To capture this common phenomenology, many attempts have been made to find a unifying description of yielding, and possibly predict the conditions under which industrially-relevant samples start to flow. The research of the basic ingredients needed to reproduce the observed yielding behavior has been helped by elastoplastic models [1, 3], mode-coupling theories [4, 5], as well as other theoretical frameworks, such as the soft glassy rheology [6, 7] and fluidity models [8, 9]. Starting from a minimal description of the dynamics common to all arrested systems, these models capture the evolution of the viscoelastic moduli as a function of the strain amplitude $\gamma_0$. This evolution is often gradual, and it takes place over an extended range of $\gamma_0$, apparently with no univoquely well-defined yield strain. Yet, the extent of this range may vary, depending not only on the sample studied, but also on many experimental parameters such as sample age, previous mechanical history, deformation geometry, boundary conditions [10]. Controlling

the transition between gradual and abrupt yielding, sometimes referred to as ductile and brittle yielding in analogy to failure of hard materials, is of paramount importance in many industrial contexts. However, the physics underlying these different behaviors is still poorly understood. Significant breakthrough came recently thanks to an advanced analysis of the rheological signal, showing that even gradual yielding corresponds to the well-defined onset of non-recoverable deformation, which increases with $\gamma_0$ and becomes predominant in the regime of large deformations [11, 12].

The notion of yielding as a transition to irreversible, or unrecoverable, deformation has been particularly explored by probing the microscopic dynamics induced by shear in rheo-optical experiments and particle-based simulations. The large susceptibility typical of soft materials makes microscopic dynamics a very sensitive probe of nonlinear deformations: their acceleration with increasing $\gamma_0$ revealed that, at the microscale, yielding corresponds to a well-defined transition to irreversible particle motion [13–16]. Numerical simulations confirmed this scenario [17–23], suggesting that it can be described in the broader framework of absorbing phase transitions [24–37]. Yet, the nature of the yielding transition has been the object of a longstanding debate due to apparently contrasting results: on the one hand, some experiments reveal features typical of a discontinuous (first-order) transition, such as the abrupt jump of a microscopic order parameter [15, 18, 30, 38–41], the coexistence of states with dissimilar dynamics [30], and hysteretic behavior [42, 43], while on the other hand, other experiments probe features typical of a continuous (second-order) transition, such as the gradual evolution of the order parameter [13, 14, 17], sluggish relaxations towards steady viscoelastic moduli or microscopic dynamics after preshear [14, 44], and emerging large fluctuations [14, 15, 45].

This apparent contradiction has been recently rationalized

<footnote>$^a$Accepted for publication in Physical Review Materials, DOI:10.1103/PhysRevMaterials.9.055602</footnote>
in terms of an equation of state relating the microscopic dy-
namics to the applied deformation [46, 47]: in this framework,
yielding is regarded as a dynamic transition that is discon-
tinuous for samples deeply quenched into the glassy state,
and that becomes continuous at the glass transition. Imple-
mented on a lattice, this model predicts the coexistence of
solid-like and fluid-like relaxation modes at yielding, giving
rise to dynamic heterogeneities whose magnitude and char-
acteristic length scale depend on the glassiness of the sam-
ple. This scenario, confirmed experimentally by dynamic light
scattering and microscopy on a variety of soft glassy samples
[46], suggests that gradual yielding observed at the macro
scale may result from locally discontinuous processes: the
relationship between microscopic dynamics and macroscopic
rheology emerges therefore as a major knowledge gap towards
a comprehensive understanding of yielding in soft materials.

A route towards filling this gap is to inspect and to ratio-
nalize the appearance of spatial correlation of the dynamics at
the rheological yielding. In the recent few years, experiments
and simulations have unveiled the emergence of a dynamic
length scale in yielding soft glasses [9, 14, 17, 48–60] that
has been related to the extent of localized plastic events. This
length scale, also observed upon melting [61–64], may be sig-
nificantly larger than that of dynamic heterogeneities at rest
[65–67], and may even exceed the size of the whole system
as it flows in tight geometries, giving rise to size-dependent
rheological behavior [48]. Similarly, abrupt yielding in soft
amorphous materials has been observed as individual plastic
events organize in system-spanning shear bands [68–77]. This
behavior exhibits strong analogies with quasi-brittle failure
in disordered hard materials, such as wood, concrete, com-
posites, rocks and metallic alloys [78–85], which has been
described as a critical-like phenomenon with a characteris-
tic length scale that grows and eventually diverges as macro-
scopic failure is approached [86–91]. Similar ideas guide the
research of precursory signals anticipating granular or fric-
tional slip [92–95] and even rockfalls and earthquakes [96–98]
Conversely, when failure occurs through sparse, uncorrelated
plastic events, ductile behavior is observed [69, 99–102]: this
further confirms that the length scale characterizing the spatial
correlation of plastic events is key to determine the yielding
behavior of amorphous solids.

Yet, despite its relevance, little is known about the origin
of this length scale, its connection to the properties of the
sample, and its emergence under a given applied deforma-
tion. Structural heterogeneity has been shown to suppress
spatial correlations and to enhance ductility [99, 103–105],
but even without significant structural differences, the same
effect can be obtained through mechanical history and sample
equilibration [20, 55, 57, 69, 106–110]. Recent experimen-
tal works on soft and metallic glasses [108, 111–113], high-
lighed the role of structural and mechanical disorder for the
sharpness of mechanical yielding. Theoretical and numerical
works [69, 104, 114–117] have rationalized this result, show-
ing that more disordered samples yield through a larger num-
ber of smaller plastic events, resulting in the smooth evolu-
tion of the stress. By contrast, better-annealed glassy samples,
characterized by a reduced disorder in local mechanical prop-

erties [114], yield more abruptly, and tend to develop shear
bands or fractures [69, 115, 118, 119].

Here, we study the effects of disorder on the sharpness of
the dynamic transition observed at the microscale upon yield-
ing. Building upon the on-lattice model developed to repro-
duce the stroboscopic dynamics measured by rheo-DLS [46],
here we focus on the role that coupling disorder plays in de-
termining both the nature of the transition (continuous vs. dis-
continuous) and the extent of dynamic heterogeneities, which
are aspects and consequences of our model that are not dis-
cussed in [46]. In absence of disorder, our model, yields
a law of correspondent states equivalent to trajectories on a
cusp catastrophe manifold, a well-known class of problems in-
cluding equilibrium liquid-vapour phase transitions. By gen-
eralizing the model proposed in [46] to account for an arbi-
trary range of coupling between shear-induced relaxations we
briefly discuss in mean field approximation the impact of this
range on the yielding transition. We show how introducing a
finite disorder in our model entails a qualitative change, to a
continuous and rounded transition, whose extent is radically
affected by the magnitude of the disorder. Finally, and more
importantly, in the present work we show that a spatial cor-
relation length ξ emerges spontaneously from the coupling
between disorder and bifurcating dynamics. With vanishing
disorder, ξ diverges and yielding becomes discontinuous, sug-
gesting that the abruptness of yielding can be rationalized in
terms of a diverging lengthscale of dynamic heterogeneities.
This is what we consider the most significant result of the
present work, which substantially enriches the physics emerg-
ing from the experiments presented in [46]

The reminder of this paper is structured as follows. In sec-
tion II we highlight the experimental grounds on which our
model is built. In section III we formulate the model in its
general form and we detail its physical bases. In section IV
we discuss the results obtained in mean field approximation,
and we formulate a law of corresponding states allowing to
map distinct fluids on the same equation of state. In section V
we discuss the role of disorder in our model by implementing
it on a 2D lattice. Simulations indicate that with increasing
disorder, yielding undergoes a sharp-to-gradual transition that
exhibits strong system size effects. This transition is rational-
ized in section VI by recasting our model as an elementary
catastrophe, and showing that sharp yielding is a manifesta-
tion of metastability, quantified in terms of an effective po-
tential. In section VII we show that yielding is associated to
a dynamic correlation length that grows as coupling disorder
is reduced, and eventually exceeds the system size, entailing
the transition to abrupt failure, which therefore emerges as a
finite-size effect. Finally, in section VIII we make some con-
cluding remarks, we summarize the key results and we high-
light the potential perspectives of our work.

II. EXPERIMENTAL FINDINGS

Dynamic Light Scattering (DLS) quantifies the microscopic
dynamics via the intensity correlation function g₂(τ) − 1,
which decays from 1 to 0 as microscopic displacements grow
<figure><img src="image_2.png" /><figcaption>FIG. 1. Yielding as a dynamic transition. a) Sketch of the parallel plate geometry employed in rheo-DLS experiments. b) Schematic of the stroboscopic collection of speckle images during one oscillatory experiment: One image is recorded every oscillation period T. The lag time is then a multiple of a period, τ = nT. c): Intensity autocorrelation functions under oscillatory shear, for concentrated Ludox nanoparticles (φ = 0.45), plotted vs the normalized time delay ωτ, with 2π/ω the period of the oscillations. The strain amplitude γ₀ increases from blue to red shades, spanning the rheological yielding transition. Symbols: experimental data. Lines: fits using Eq. 1. The compressed exponent βs of the slow mode is 1.9±0.05, fixed by the spontaneous dynamics measured on the samples at rest. The stretching exponent βf is 0.9±0.05. The fitting parameters for the same samples are shown in panel d). Left axis: γ₀ dependence of the rates Γs (slow mode, blue squares) and Γf (fast mode, red circles). Right axes: relative amplitude χ of the slow mode (stars). Lines: numerical results of the general model for the yielding transition discussed in the text. The best fit parameters are reported in [46, 47]. In panels d) and e), the gray shaded rectangles indicate the range of γ₀ from the onset of the third harmonic in stress response to the strain crossover where G' = G''.</figcaption></figure>

beyond a length scale π/q ≈ (0.1 − 1) µm, set by the scattering vector q [120]. To investigate the microscopic dynamics under deformation through rheo-DLS, we apply an oscillatory shear with angular frequency ω and strain amplitude γ₀, and stroboscopically measure g₂(τ) − 1, for lag times τ that are multiples of the oscillation period, sketched in Figure 1-a,b. Here we report only a small subset of the full dataset presented in Ref. [46], where we investigated a variety of systems—including microgels, Ludox particles, and emulsions—and explicitly explored the role of the oscillatory strain frequency. The data shown in Figure 1 have been obtained for a glassy sample of Ludox particles at ω = 3.14 rad/s. We find that for small enough γ₀ the dynamics are independent of the applied shear, and that dynamics accelerate with growing strain at larger γ₀. Concomitantly, the shape of g₂ − 1 evolves from a compressed exponential at low γ₀ to a stretched exponential at large γ₀, with an intermediate regime where correlation functions exhibit a two-step decay well described by a linear combination of the two modes, as shown in Figure 1-c. To model this evolution, correlation functions measured at all strain amplitudes are fitted by the following expression:

$$g_2(\tau) - 1 = \left\{ \chi \exp \left[ -(\Gamma_s \tau)^{\beta_s} \right] + (1 - \chi) \exp \left[ -(\Gamma_f \tau)^{\beta_f} \right] \right\}^2, \quad (1)$$

where Γs and Γf are the relaxation rates of the slow and fast decay modes, respectively; χ is the relative weight of the slow decay mode, and βf and βs are the stretching and compressing exponents of the two modes, shared by all correlation functions. The fit results suggest that, as γ₀ is increased, three distinct regimes can be identified:

i For γ₀ < 5%, the correlation function decays as a single compressed exponential (χ = 1), with an exponent βs ≈ 1.9 and a slow relaxation rate Γs, hardly dependent on γ₀. We find that Γs is approximately equal to Γ₀ = 10⁻⁴ s⁻¹, the relaxation rate measured in absence of shear. Therefore, we refer to this slow, compressed-exponential relaxation as the solid-like dynamics.

ii For 5% < γ₀ < 9%, a second, faster mode, characterized by a stretched exponential relaxation, with βf ≈ 0.9, adds to the spontaneous relaxation mode, and becomes increasingly dominant with increasing γ₀, as χ decays from 1 to 0.

iii Finally, for γ₀ > 9%, the correlation functions decay as single stretched exponentials (χ = 0), with a relaxation rate Γf growing as a power law of γ₀, with a sample-dependent exponent n ≈ 3. Because such stretched-exponential dynamics are similar to those found in dense suspensions just below the glass transition [121],
we refer to this faster relaxation mode as fluid-like dynamics.

Comparison with the rheological response indicates that these three regimes correspond to: (i) the linear regime, characterized by a purely harmonic stress response and strain-independent viscoelastic moduli, with $G' > G''$; (ii) the yielding regime, with an increasingly inharmonic stress response and an overshoot in the first-harmonic $G''(\gamma_0)$; (iii) the fluidized regime, with $G'' > G'$ and both nonlinear moduli decaying as power laws of $\gamma_0$, as shown in Figure 1-e.

Comparing the experimental results reported in Figure 1 to the result of analogous experiments on different soft glassy materials [46], we find that this scenario does not depend on sample details such as the nature of the microscopic constituents, their soft repulsive interactions or their concentration, as long as the sample exhibits glassy dynamics at rest. The most remarkable feature changing from sample to sample is the abruptness of the transition, quantified by the extent of regime (ii). For instance, we find that, for a given class of samples, say, charged nanoparticles as in Fig. 1, more concentrated samples exhibit more abrupt yielding. A similar trend is found as one given sample is probed at increasingly large angular frequency $\omega$, or if it is aged for longer time prior to being tested. These findings suggests that a deeper understanding of the dynamic transition observed at the microscale can reveal the very nature of yielding: this motivates the search for a general model able to describe this behavior.

### III. LATTICE MODEL WITH FACILITATED ADVECTION

To model the microscopic relaxation rate $\Gamma$, we assume that it can be written as the sum of two independent processes: the spontaneous dynamics and the shear-induced dynamics, characterized by relaxation rates $\Gamma_0$ and $\Gamma_{sh}(\gamma_0)$, respectively. Previous works [122, 123] have postulated a generic power-law dependence of $\Gamma_{sh}$ on the macroscopic shear rate. Here, we adapt this ansatz to oscillatory rheology, taking $\dot{\gamma}_0 = \omega\gamma_0$ as a representative shear rate, and using the experimental observation that $\Gamma_{sh} \propto \omega$ to properly capture the dependence on $\omega$ [46]. We obtain:

$$ \Gamma = \Gamma_0 + K\omega\gamma_0^l. \qquad (2) $$

This result captures well the limits of small and large $\gamma_0$, but fails describing the well-defined transition found in experiments, rather predicting a continuous, smooth acceleration of $\Gamma$ with increasing $\gamma_0$, as shown by the black dashed line in Fig. 1-d. To improve the model predictions in the yielding region, we consider the role of dynamic and structural heterogeneities in the sample [124–126] arising from dynamic correlations between nearby sample regions. Indeed, glassy dynamics at rest are known to be governed by dynamic facilitation [125, 126]: portions of the system with faster dynamics promote structural relaxation of nearby regions. By analogy, we introduce a facilitated advection (FAd) model that couples the shear-induced dynamics of nearby sample regions in an

Ising-like fashion. We implement this model on a lattice and we assume that the shear-induced relaxation of a generic site $i$ does not depend only on the applied macroscopic shear field but also on the local environment, that screens the imposed global deformation and delays local rearrangements:

$$ \Gamma_i = \Gamma_0 + K[\tau_{sh}(\gamma_0) + \tau_F(i, j)]^{-1}. \qquad (3) $$

Here $\tau_{sh}(\gamma_0) = (\omega\gamma_0^l)^{-1}$, whereas $\tau_F(i, j)$ is a characteristic time that depends on the dynamics of two adjacent lattice sites $i$ and $j$, and represents the leading-order expansion of the perturbation of the shear-imposed time scale due to local dynamics. $\Gamma_i$ is therefore the local relaxation of the $i$-th lattice cell representing a mesoscopic region of the sample, that reduces to the rate given by Eq. 2 only when the dynamics are uncorrelated. This is the case of ergodic fluids with no dynamical heterogeneities and very poor spatial correlation of local dynamics, characterized by $\tau_F(i, j) = 0$. More generally, we model the effect of FAd by introducing a term, $\tau_F(i, j)$, that penalizes strong differences in the dynamics of nearby cells, while preserving the small and large amplitude limits of Eq. 2:

$$ \Gamma_i = \Gamma_0 + \omega K \left( \frac{1}{\gamma_0^l} + \frac{\omega^2}{N_z} \sum_{j=1}^{N_z} \frac{\alpha_{ij}}{\Gamma_i \Gamma_j} \right)^{-1} \qquad (4) $$

Here $\alpha_{ij}$ are positive constants coupling the dynamics of the site $i$ to the ones of its $N_z$ nearest neighbors. In the limit of $\alpha_{ij} = 0$, Eq. 4 reduces to the original ansatz of Eq. 2 with $\Gamma_i = \Gamma$ for all lattice sites. With $\alpha_{ij} > 0$, the local rates vary between lattice sites: lower-than-average $\Gamma_j$ result in larger FAd term for its neighbors, decreasing their rate relative to sites in faster-relaxing neighborhoods. Therefore, the FAd term introduced in Eq. 4 suppresses differences between neighbor lattice sites, like in dynamic facilitation [127, 128].

### IV. MEAN FIELD APPROXIMATION

To illustrate the new features introduced by FAd, we first set all coupling constants to the same value $\alpha$, and we study the resulting mean-field equation of state (EOS) for the average relaxation rate $\Gamma(\gamma_0)$:

$$ \Gamma = \Gamma_0 + \omega K \gamma_0^l \left( 1 + \frac{\alpha \omega^2 \gamma_0^l}{\Gamma^2} \right)^{-1} \qquad (5) $$

For any given set of model parameters $\{\Gamma_0, \omega, K, n, \alpha\}$, this expression reduces to Eq. 2 both in the limit of small strains, $\gamma_0^l \ll \Gamma_0^2/\alpha\omega^2$, and of large strains, $\gamma_0^l \gg \alpha/K^2$, where $\Gamma^2 \sim \omega^2 K^2 \gamma_0^{2n} \gg \alpha\omega^2\gamma_0^l$. However, for $K \lesssim \alpha\omega/\Gamma_0$, there exists an intermediate range of strain amplitudes where the model with FAd exhibits a qualitatively new behavior. To show this, we recast Eq. 5 in a form that is formally identical to the Van der Waals (VdW) equation of state for real gases, with $\Gamma = 1/\gamma_0^l$ and $\omega K$ playing the role of volume, pressure and temperature, respectively [129]:
$$ (\Gamma - \Gamma_0) \left( \frac{1}{\gamma_0^l} + \frac{\alpha \omega^2}{\Gamma^2} \right) = \omega K \quad (6) $$

Guided by this analogy, we represent the mean-field equation of state (EOS) as curves in the $1/\gamma_0$ vs $\Gamma$ plane, analogous to isotherms in $p-V$ diagrams for real gases. These are curves defined by the equation:

$$ \frac{1}{\gamma_0} = \left[ \frac{\omega K}{\Gamma - \Gamma_0} - \alpha \left( \frac{\omega}{\Gamma} \right)^{\beta} \right]^{1/n}, \quad (7) $$

where the exponent $\beta > 1$ is introduced to generalize Eq. 4 (in which $\beta = 2$) to longer-range dynamic coupling, extending to $\beta - 1$ lattice sites (see Appendix A for the expression on lattice). Like in real gases, the EOS of Eq. 7 describes $\gamma_0^{-1}(\Gamma)$ curves that decrease monotonically for sufficiently large $K$, and that become nonmonotonic as $K$ is reduced below a critical value, as shown in Fig. 2.

In a typical oscillatory strain sweep, the system is prepared in an aged and stationary glassy state at rest, and then it is subject to oscillatory deformation at increasingly large strain amplitude. We model the evolution of the microscopic dynamics by fixing $K$ and $\alpha$, and following the increase in $\Gamma$ as $\gamma_0$ is increased in Eq. 7. This corresponds to following a path going from the top-left to the bottom-right corner of the $\gamma_0^{-1} - \Gamma$ diagram in Fig. 2. For small $\gamma_0$, corresponding to the linear regime probed by rheology, $\Gamma(\gamma_0)$ is nearly constant as the EOS exhibits a very steep, solid-like branch analogous to that of nearly-incompressible VdW liquids. The presence of a local minimum in Eq. 7 sets an upper bound, $\gamma_{th}^+$, for the strain amplitude that can be attained along this section of the EOS. Beyond $\gamma_{th}^+$, $\Gamma$ can be no longer increased continuously with increasing $\gamma_0$. Instead, a small increase of $\gamma_0$ causes $\Gamma$ to jump discontinuously to the fluid branch of the EOS, $\Gamma_f(\gamma_0)$, shown as a red line in Fig. 2. This discontinuous jump represents an abrupt yielding event. Conversely, if $\gamma_0$ is decreased starting from the fluidized state in a reverse strain sweep, the system will follow the fluid branch until the local maximum at a second threshold strain amplitude, $\gamma_{th}^-$, where it discontinuously jumps back on the solid branch. Therefore, in mean field, our model predicts an abrupt, discontinuous yielding with hysteresis for a range of model parameters such that Eq. 7 gives nonmonotonic $\gamma_0^{-1}(\Gamma)$.

While still unable to fully capture the coexistence of solid-like and liquid-like dynamics for a finite range of $\gamma_0$, this mean field model substantially improves Eq. 2 by introducing a well-defined dynamic transition at yielding. Remarkably, this is achieved with the introduction of one single additional parameter, the coupling constant $\alpha$, which can be adjusted to quantitatively reproduce the magnitude of the dynamic acceleration measured upon yielding for each sample. In an attempt to rationalize the impact of model parameters on the predicted dynamic acceleration, we refer once more to the formal analogy with VdW gases: there, the liquid-vapor phase transition of different gases at different temperatures can be described in a universal manner by the so-called law of corresponding

<figure><img src="image_3.png" /><figcaption>FIG. 2. Inverse strain amplitude $\gamma_0^{-1}(\Gamma)$ obtained according to equation 7. Transitions following an ascending ramp ($\Delta\gamma_0 > 0$, increasing $\gamma_0$ from the solid branch) and a descending ramp ($\Delta\gamma_0 < 0$, decreasing $\gamma_0$ from the fluid branch) occur at the two extrema of the curve for $\gamma_0 = \gamma_{th}^+$ and $\gamma_0 = \gamma_{th}^-$, as indicated by the horizontal blue and red arrows respectively. Here $K = 300$, $n = 3$, $\beta = 2$, $\alpha = 0.17$, $\omega = 0.66 \text{ rad/s}$ and $\Gamma_0 = 10^{-4} \text{ s}^{-1}$. The gray dashed lines show as reference the EOS for $\alpha = 0$ with the other parameters left unaltered. Inset: Reduced equation of state $\gamma_r^{-1}(\Gamma_r)$ obtained for $n = 3$, $\beta = 2$ and different values of $K_r$ as indicated in the figure. The horizontal dotted lines indicate the predicted yielding for ascending strain amplitude ramps, $\Delta\gamma_0 > 0$ (decreasing $\gamma_r^{-1}$), starting from a sample at rest.</figcaption></figure>

states, stating that the behavior of a real gas in a given thermodynamic state ($p, V, T$) only depends on the distance from the critical point. In an analogous way, here we seek for a normalized description in terms of rescaled variables {$\gamma_r^{-1} = \gamma_c/\gamma_0$, $\Gamma_r = \Gamma/\Gamma_c$, $K_r = K/K_c$}, where ($\gamma_c^{-1}, \Gamma_c, K_c$) defines the critical point, i.e. the horizontal inflection point of the EOS, satisfying the condition: $\partial\gamma_0^{-1}/\partial\Gamma = \partial^2\gamma_0^{-1}/\partial\Gamma^2 = 0$. Imposing this condition in Eq. 7 we obtain:

$$ \Gamma_c = \Gamma_0 \frac{\beta + 1}{\beta - 1} \quad (8.a) $$

$$ \gamma_c^{-1} = \left[ \alpha \left( \frac{\omega}{\Gamma_0} \right)^{\beta} \left( \frac{\beta - 1}{\beta + 1} \right)^{\beta+1} \right]^{1/n} \quad (8.b) $$

$$ K_c = 4\alpha \left( \frac{\omega}{\Gamma_0} \right)^{\beta-1} \frac{\beta(\beta-1)^{\beta-1}}{(\beta+1)^{\beta+1}} \quad (8.c) $$

Rewriting Eq. 6 in terms of these rescaled variables, we obtain a reduced equation of state:

$$ \left( \Gamma_r - \frac{1}{\mu} \right) \left( \gamma_r^{-n} + \frac{\mu}{\Gamma_r^\beta} \right) = \frac{4\beta}{\beta^2 - 1} K_r \quad (9) $$
where $\mu = (\beta + 1)/(\beta - 1)$, and $K_r$ plays the role of a reduced temperature: $K_r < 1$ defines the range of model parameters for which "iso-$K_r$" $\gamma^{-1}(\Gamma)$ curves are nonmonotonic and therefore describe yielding soft solids. Quantitative differences in the yielding behavior of different materials are incorporated in the value of $K_r$: for $K_r \lesssim 1$, $\Gamma(\gamma_0)$ undergoes a sizeable increase along the solid branch prior to yielding, upon which it experiences a relatively small jump to reach the fluid branch. This behavior mimics the gradual yielding observed for soft materials close to the glass transition [15]. Conversely, for smaller $K_r$, $\Gamma(\gamma_0)$ is nearly constant in the solid branch, and the discontinuity at $\gamma_{th}^+$ grows, mimicking the larger dynamic acceleration upon yielding that is observed in experiments [15] and numerical simulations [69] on samples that have been well equilibrated deep in the glassy regime. Therefore, we refer to $1 - K_r$ as the sample *glassiness*. From Eq. 8.c, we find that $K_r$ decreases with both increasing $\alpha$, indicating that stronger FAd coupling enhances dynamic acceleration upon yielding, and with increasing $\omega$, indicating that yielding is time-dependent, like the glassy state itself, which is characterized by a finite relaxation time $\Gamma_0$. As longer timescales are probed by reducing the Deborah number $De = \omega/\Gamma_0$, the sample will be effectively less glassy, with a progressively smaller jump of the relaxation rate $\Gamma$ at yielding, until $\omega/\Gamma_0 = K/4\alpha$, a threshold beyond which the model predicts fluid-like behavior at all strain amplitudes, as $K_r > 1$.

The amount of dynamic acceleration prior to yielding and the jump of relaxation rate at the transition also depend on $\beta$, which represents the spatial range of FAd coupling. In particular, for increasing $\beta$ we find that the extent of the discontinuity of $\Gamma$ at yielding increases, and that the value of $\Gamma$ just before yielding decreases, approaching $\Gamma_0$, as detailed in Appendix A. The result of an augmented range of FAd coupling is in line with enhanced fluidization and stress drops that emerge in numerical simulations of deeply quenched glasses and percolating networks, where long-ranged stress correlations determine rigidity at rest and possibly brittleness under shear [69, 130, 131].

On a final note, we observe that as $K_r$ is decreased, the normalized threshold strain $\gamma_{r,th}^+$ increases. Interestingly, we find that $\gamma_{r,th}^+$ diverges for a critical, finite value of $K_r$: $K_{rm} = (1 + \beta^{-1})^{\beta+1}/4$. This singularity cannot be properly discussed in the framework of the mean field model, but hints at a second transition, to a different, more catastrophic failure mechanism.

## V. EFFECT OF DISORDER: NUMERICAL SIMULATIONS

The mean-field model described by Eqs. 5-9 accounts for both the linear and the fully-fluidized regimes, predicting an abrupt transition between the two states. This is in contrast with the experimental observation that yielding may manifest as a gradual process, with fast and slow relaxation modes coexisting in a finite range of strain amplitudes [46]. To account for this coexistence, we discuss model predictions beyond mean field, by considering the effect of structural disorder

on the local heterogeneity of shear-induced relaxations in glassy materials [132–134]. The disorder in the coupling strength therefore reflects the underlying structural spatial variability in a glass, such as polydispersity and spatial packing heterogeneity, which has a huge impact on dynamical heterogeneities (e.g., different local relaxation behaviors), contributing to the observed facilitation-like interactions between neighboring regions. Hereafter, without loss of generality, we will focus on the particular case $\beta = 2$, describing short-range FAd coupling.

<figure><img src="image_4.png" /><figcaption>FIG. 3. **Effect of disorder in numerical simulations.** (a) $\gamma_r^{-1}(\Gamma_r)$ diagrams obtained by averaging local rates $\Gamma_i$ in simulations with $K_r = 0.86$ and disorder $d$ increasing from $10^{-5}$ to 0.175 (green to blue to red shades). Dashed lines represent simulations where yielding was abrupt due to finite size effects. Dotted line: mean-field equation of state, eq. 9, with $\beta = 2$ and $K_r = 0.86$. (b) Histogram of local rates $\Gamma_i$, with $K_r = 0.86$ and $\gamma_r$ increasing from top to bottom, as denoted by the colored arrows: $1/\gamma_r = 0.69$ (purple), 0.43 (gold), 0.275 (cyan), 0.07 (brown). Vertical dashed black lines represent mean field values $\Gamma^s$, $\Gamma^f$. Bottom panel: snapshots of 2D configuration of local rates. Rate increases from blue ($\Gamma = 1$) to red shades ($\Gamma = 100$) in a logarithmic fashion, as shown by the color bar. Colored arrows mark strain amplitudes shown in panel b, dashed and dotted black vertical lines mark the reference threshold strain amplitudes predicted by mean-field, $\gamma_{th}^+$, and by the Maxwell rule, $\gamma_{th}^{(M)}$, respectively.</figcaption></figure>

In presence of disorder, the local rates $\Gamma_i$ vary from site
to site, and to find their values, one needs to solve N coupled nonlinear equations, each of which may have multiple solutions. Tackling this problem analytically goes beyond the scope of this paper, though a simple and preliminary analysis based on the central limit theorem [135] helps in setting the scenario, as discussed in Appendix B. Here, to get insights into the new behavior arising from disorder, we study it numerically, on a $N \times N$ square lattice with periodic boundary conditions. Each site, $i$, is assigned a local relaxation rate, $\Gamma_i$, and each couple of neighbor sites is assigned a coupling constant, $\alpha_{ij}$, randomly drawn from a probability distribution, $P(\alpha_{ij})$. For a given strain amplitude $\gamma_0$, solving our model numerically consists in finding a configuration of local rates $\{\Gamma_i\}$ that satisfies Eq. 4 for all sites. In our implementation, $\{\Gamma_i\}$ is approached iteratively: starting from an initial configuration of local rates $\{\tilde{\Gamma}_i^{(0)}\}$, we compute a set of target rates $\{\tilde{\Gamma}_i^{(1)}\}$ through Eq. 4 using $\{\tilde{\Gamma}_i^{(0)}\}$ to evaluate the FAd coupling term. If all $\tilde{\Gamma}_i^{(0)}$ satisfied Eq. 4, we would obtain $\{\tilde{\Gamma}_i^{(1)}\} = \{\tilde{\Gamma}_i^{(0)}\} = \{\Gamma_i\}$. If this is not the case, we use $\{\tilde{\Gamma}_i^{(1)}\}$ to compute a new set of target rates $\{\tilde{\Gamma}_i^{(2)}\}$, and keep iterating until the algorithm converges on a stationary set of rates $\{\tilde{\Gamma}_i^{(s)}\}$, for which the loss function $\mathcal{L}(s) = \sum_i \left[ \frac{\log \tilde{\Gamma}_i^{(s+1)} - \log \tilde{\Gamma}_i^{(s)}}{\log \tilde{\Gamma}_i^{(s)}} \right]^2$ vanishes: we then take $\{\Gamma_i\} = \{\tilde{\Gamma}_i^{(s)}\}$. Although close to yielding many iterations may be needed to find this solution, we find that this algorithm always converges, as detailed in Appendix C.

Numerical simulations with disorder can quantitatively reproduce the gradual yielding observed in experiments, as shown by the solid lines in Figure 1d. To highlight the role of glassiness and disorder on the yielding transition, in the following we fix $n = \omega = \Gamma_0 = K = 1$, and we draw the coupling parameters $\alpha_{ij}$ from a Gaussian distribution with average $\bar{\alpha}$ and variance $\sigma_{\bar{\alpha}}^2$, which we use to define the disorder parameter $d = \sigma_{\bar{\alpha}}^2/\bar{\alpha}^2$. To avoid unphysical negative coupling constants, we truncate $P(\alpha_{ij})$ to ensure that $P(\alpha_{ij} \le 0) = 0$. To mimic the effect of a strain sweep applied on a sample initially at rest, we start from a solid-like condition, with $\tilde{\Gamma}_i^{(0)} = \Gamma_0$, and solve the numerical model with a set of strain amplitudes $\gamma_0$ spanning the physically interesting range across the yielding transition. For each $\gamma_0$, we then compute the average relaxation rate $\bar{\Gamma} = \langle \Gamma_i \rangle$, and we compare it to the mean-field EOS. In the limit of $d \to 0$, we find that our simulations reproduce the mean field results: $\bar{\Gamma}$ closely follows Eq. 6, with $\alpha = \bar{\alpha}$ and $\beta = 2$, and it jumps abruptly as $\gamma_0$ reaches a threshold value $\gamma_y$ very close to $\gamma_{th}^+$, the local minimum of Eq. 7, as shown by the light green dashed lines in Figure 3a.

Increasing $d$, we find that $\gamma_y$ decreases, seemingly approaching a second characteristic strain, marked as $\gamma_{th}^{(M)}$ in Fig. 3a. This suggests that it exists a regime of strain amplitudes smaller than $\gamma_{th}^+$ for which the solid branch of the mean-field EOS is actually metastable, similarly to what occurs in first order thermodynamic phase transitions [129] between binodal and spinodal lines. Hence, a discontinuous jump of dynamics—reminiscent of a first-order transition—does appear in finite systems, provided that the coupling disorder is

sufficiently small. This behavior aligns with the sharp transitions observed in finite-size athermal quasi static simulations of highly annealed systems [41]. As $d$ is further increased, simulations show a qualitative change, to a continuous, gradual increase of $\bar{\Gamma}(\gamma_0)$, as shown by solid lines in Figure 3a. To understand the origin of this qualitative change, we analyze the distribution of local relaxation rates, $P(\Gamma_i)$. For small $d$, we find that $P(\Gamma_i)$ is sharply peaked around $\bar{\Gamma}$, with a variance proportional to $d$, reflecting the heterogeneity of local couplings. As $d$ increases and the discontinuity in $\Gamma(\gamma_0)$ disappears, we find that, within a finite interval of strain amplitudes around $\gamma_s$, $P(\Gamma_i)$ becomes bimodal: a second peak appears at larger $\Gamma_i$, roughly located at $\Gamma^f(\gamma_0)$, the fluidized branch of the mean-field EOS, as shown in Fig. 3b. In this regime, two distinct populations of local rates can be identified: slower ones, following the solid branch of the EOS, and faster ones, which jump to the fluidized branch without entailing the transition of the entire lattice.

As $\gamma_0$ is increased, we find that the two peaks in $P(\Gamma_i)$ change their relative amplitude: the fluid peak increases in amplitude and slightly shifts towards larger rates, while the solid peak gradually disappears without significant shift in its position. This indicates that the continuous evolution of $\gamma^{-1}(\bar{\Gamma})$ reported in Fig. 3a actually reflects a spatially heterogeneous process, which is locally discontinuous. This behavior recalls the dynamic coexistence of fast and slow relaxation modes observed experimentally in Fig. 1, therefore we refer to the range of model parameters for which $P(\Gamma_i)$ is bimodal as the coexistence region. As $\gamma_0$ is increased within this region, an increasing fraction of local rates $\Gamma_i$ jumps from the solid to the fluidized branch of the EOS, reproducing the continuous decay of $\chi$ in Figure 1d.

For the most disordered configurations, this coexistence extends over a wide range of strain amplitudes, as liquid-like local rates appear at small $\gamma_0$, while solid-like ones persist until larger $\gamma_0$: this entails a more gradual increase of $\bar{\Gamma}(\gamma_0)$ as shown by the red solid lines in Fig. 3. In this sense, we observe that increasing disorder in our model produces more gradual, "ductile-like" yielding. Its ability to control the abruptness of yielding through one single model parameter makes our model a convenient playground to study the impact of disorder on the yielding transition, which has been the focus of recent theoretical work [69, 104, 114]. Comparing dashed and solid lines in Fig. 3a, it is tempting to conclude that our model predicts a transition from continuous to discontinuous yielding at a well-defined critical disorder $d^*$. To better understand this transition, we look at how the local rates $\Gamma_i$ are distributed in space, with particular focus on the coexistence region. We find that slow and fast-relaxing sites are not distributed randomly, following the spatially-uncorrelated coupling constants, but they are organized in finite-sized domains separated by relatively sharp interfaces, as shown in the color maps of Fig. 3. In particular, we find that the size of these domains grows with decreasing $d$, and becomes comparable with the size of the whole lattice for $d \approx d^*$, at the transition to discontinuous yielding.

We interpret this result in analogy with nucleation and growth in first-order phase transitions: the FAd term cou-
ples nearby lattice sites, penalizing gradients in the local rate, therefore it favors the formation of domains with homogeneous local dynamics. In absence of disorder, this penalty produces the hysteresis shown in Figure 2, analogous to the behavior observed in metastable systems close to a first-order transition such as supercooled liquids [136] or ferromagnetic materials under an external magnetic field [137], and also reported for colloidal glasses under oscillatory shear [138]. With disorder, we observe that solid-to-liquid transition is first nucleated at sites where the local coupling is weak, which then act as nucleation seeds. If the disorder is small, this nucleation occurs when $\gamma_0$ is in the metastable region, and triggers the abrupt transition of the entire lattice. Conversely, for large disorders, liquid domains form at $\gamma_0 < \gamma_y$, and coexist with solid-like domains pinned to sites with strong enough local coupling. As disorder increases, $P(\alpha_{ij})$ gets broader, and the number density of sites susceptible of pinning solid or liquid domains increases. This yields a larger number of smaller-sized domains, coexisting for a wider range of strain amplitudes. Conversely, for smaller disorders, the reduced variance of local couplings produces larger-sized domains, only observable in a narrower range of $\gamma_0$. In the extreme case of very small disorder, we interpret the transition to discontinuous yielding as a result of the scarcity of sites with weak enough coupling to seed the nucleation of liquid domains at $\gamma_0 < \gamma_y$. This interpretation is analogous to the one emerging in numerical simulations of dense 2D soft disks, where brittleness is related to the density of "softer" elements [104].

To test this hypothesis in our model, we select a representative disorder, $d = 0.02$, and we run a series of simulations on increasingly smaller subsets of a 1024 × 1024 lattice with fixed $\alpha_{ij}$. We find that in full-scale simulations $\bar{\Gamma}(\gamma_0)$ increases gradually, and that slow and fast domains coexist in a finite range of strains around $\gamma_0 \approx 17$. We then extract a subset of the lattice with $N = 512$ and run the simulation again. The resulting $\bar{\Gamma}(\gamma_0)$ is fully consistent with the one obtained on the larger lattice, and slow and fast domains are reproduced with the same spatial patterns, as shown in the bottom rows of Figure 4a. The same result holds for $N = 256$, proving that in this range of $N$ finite size effects are negligible. However, we find that, as $N$ if reduced below a threshold value $N^* \sim 128$, roughly corresponding to the size of the largest domains, $\bar{\Gamma}(\gamma_0)$ becomes discontinuous, yielding is delayed to larger $\gamma_0$, and coexistence is no longer observed. This corroborates the hypothesis that the coexistence of solid-like and fluidized domains is intrinsic to our model with disorder, and that the discontinuous jump observed for smaller $N$ or $d$ is an effect of the finite size of the simulation lattice.

To better assess these finite size effects, we repeat this series of simulations for various disorders, and for each disorder we run many independent simulations, as detailed in Appendix D. Close to the transition between abrupt and gradual yielding, we observe that, depending on the randomly generated configuration of $\alpha_{ij}$, some simulations exhibit gradual yielding with coexistence, and some do not. For each $N$ and $d$, we then compute the fraction, $\Psi$, of simulations exhibiting gradual yielding. We find that $\Psi$ increases with both $N$ and $d$, reflecting the higher probability to find weakly-coupled

<figure><img src="image_5.png" /><figcaption>FIG. 4. **Effect of lattice size.** (a) snapshots of 2D simulations on increasingly small subset of an original 1024x1024 lattice (bottom row). Local rate increases from blue to red shades. Simulation parameters are: $K_r = 0.86$ and $d = 0.02$. The white dashed square in each left column snapshot represent the subset of the lattice that has been used for the smaller scale simulations reported in the row above. (b) Symbols: fraction of simulations exhibiting gradual yielding, $\Psi$, as a function of quenched disorder and for different lattice sites, increasing from $N = 10$ (purple) to $N = 512$ (red). Lines: Weibull cumulative distributions fitted to the data. (c) Weibull fit parameters: shape function, $p$ (red, right axis) and threshold disorder, $d^*$ (black, left axis), as a function of simulation lattice size.</figcaption></figure>

sites able to nucleate fast domains before the solid state becomes metastable. For any given $N$, we find that $\Psi_N(d)$ has a sigmoidal shape well-described by a cumulative Weibull function, $1 - \exp[-(d/d^*)^p]$, further confirming that the onset of gradual yielding with coexistence depends on the low-$\alpha$ tail of the distribution of coupling constants, which makes it an extreme value problem, analogous to rupture nucleated by defects [119, 139–141]. Fitting $\Psi_N(d)$, we extract the characteristic disorder $d^*$ and the exponent $p$. We find that the larger variability of smaller-scale simulations reflects in broader Weibull distributions, characterized by a lower exponent $p$. Moreover, we find that $d^*$ decreases with increasing $N$, confirming that for smaller disorders, larger-scale simulations are needed to capture the larger domains coexisting
at the yield point, as shown in Fig. 4c. We conclude that in the thermodynamic limit, $N \to \infty$, our model predicts gradual yielding with coexistence for all finite disorders, and that simulations showing discontinuous yielding, marked by dashed lines in Figure 3a, are biased by finite size effects.

## VI. ELEMENTARY CATASTROPHE AND MAXWELL RULE

Restricting our analysis to gradually-yielding simulations, we observe that decreasing disorder (red to blue solid lines in Figure 3a) highlights the existence of a well-defined yield strain, $\gamma_{th}^{(M)}$, distinct from the local extrema of the mean field EOS, $\gamma_{th}^{\pm}$. Although our simulation data can only track gradual yielding down to $d = 0.02$, our results strongly suggest that in the limit of infinitely large simulation lattices, $N \to \infty$, a discontinuity of $\bar{\Gamma}(\gamma_0)$ at $\gamma_0 = \gamma_{th}^{(M)}$ arises in the limit of $d \to 0^+$. The presence of this second yield strain was not predicted by the mean-field model discussed above, as it can only be observed with a finite disorder. However, its persistence for $d \to 0^+$ suggests that $\gamma_{th}^{(M)}$ is somehow related to some hidden property of the mean-field model itself. Once more exploiting the formal analogy with VdW gases, we speculate that this yield strain, laying in between $\gamma_{th}^{-}$ and $\gamma_{th}^{+}$, may arise from an optimum criterion analogous to the Maxwell rule, which selects the exact phase transition pressure within the interval where multiple solutions are analytically admitted, based on the minimization of the free energy. Indeed, we find that the value of $\gamma_{th}^{(M)}$ estimated using numerical simulations closely matches the one computed applying the Maxwell rule to Eq. 7: for the model parameters used in Fig. 4 we find $\gamma_{th}^{(M)} = 15.6$, closely matching simulation results with small disorders as shown in Fig. 3 and as further documented in Appendix D. This correspondence corroborates our hypothesis, suggesting the existence of an underlying potential function playing the role of a free energy. This is particularly surprising because, in contrast with VdW, here Equation 9 describes a dynamic transition in systems out of equilibrium, and is not derived from energy minimization.

To illustrate how something similar to a free energy can be derived in our model, we take $\beta = 2$, and we recast Eq. 9 in the form: $\gamma_r^{-n}\Gamma_r^3 - (\gamma_r^{-n} + 8K_r)\Gamma_r^2/3 + 3\Gamma_r - 1 = 0$. This expression takes an even simpler form in terms of a new variable, $\rho$, defined such that the critical point is in $\rho = 0$:

$$ \rho^3 + a\rho + b = 0 \qquad (10) $$

with:

$$ \rho = \frac{1}{\Gamma_r} - 1 \qquad (11.a) $$

$$ a = \frac{8K_r}{3} + \frac{1}{3}\gamma_r^{-n} - 3 \qquad (11.b) $$

$$ b = \frac{8K_r}{3} - \frac{2}{3}\gamma_r^{-n} - 2 \qquad (11.c) $$

Eq. 10 defines a manifold unifying all trajectories followed by samples at various $K_r$ subject to varying strain amplitudes. Its functional form is well-known in statistical physics and bifurcation theory, as it represents the manifold of a so-called cusp catastrophe, the same describing the liquid-vapour transition of gases [142].

To highlight its features, we show it in Figure 5 as a function of the new set of variables $\{\rho, a, b\}$. The non-monotonic EOS characterizing glassy samples defined by $K_r < 1$ correspond to paths crossing the region where the manifold folds over itself, such that for a given strain, fixing the couple $(a, b)$, multiple values of $\rho$, thus of the microscopic rate $\Gamma$, are allowed. In this region, iso-$K_r$ paths have local extrema, where the manifold is locally perpendicular to planes of constant $\rho$. Differentiating Eq. 10, we find that this occurs for $3\rho^2 + a = 0$, which defines the so-called bifurcation set, $\Sigma$ [142]. The projection of $\Sigma$ on the $(a, b)$ plane is defined by $a = -3(b/2)^{3/2}$, which has a cusp in the origin, as shown in Fig. 5b. Hence, in our model yielding falls into the broader class of critical phenomena described as a cusp catastrophe [142, 143]. Projected on the $(a, b)$ plane, strain sweeps are straight lines defined by $a = 4(K_r - 1) - b/2$, with $\gamma$ increasing with $b$. For $K_r < 1$, these lines intercept $\Sigma$ in two points corresponding to $\gamma_{th}^{\pm}$ and $\gamma_{th}^{-}$. Between these two points lays the critical strain resulting from simulations with disorder, $\gamma_{th}^{(M)}$. Unlike $\gamma_{th}^{\pm}$, this new critical strain does not correspond to a critical point of the manifold. To relate it to the mean field model, we integrate Eq. 10 to obtain a potential function $V(\rho; a, b)$ such that $\partial V / \partial \rho = 0$ yields the equation of state. For Van der Waals gases, $V$ would correspond to Helmholtz free energy. Here, we obtain:

$$ V(\rho; a, b) = \frac{1}{4}\rho^4 + \frac{a}{2}\rho^2 + b\rho \qquad (12) $$

The shape of this function is dictated by the parameters $a$ and $b$. In particular, we observe that $V$ is convex for $a > 0$, whereas for $a < 0$ its second derivative becomes negative for $a/3 < \rho^2 < -a/3$, and the function develops two minima, corresponding to the two branches of the critical set $\Sigma$. This evolution of $V$ recalls the behavior of the free energy functional in Landau theory for magnetic phase transitions in the presence of an external field [129], with $\rho$ as the order parameter. In that framework, the requirement that the free energy should be minimized provides a stability criterion that selects the state of the system when the equation of state has multiple solutions. Here, for $a < 0$, this criterion corresponds to selecting the deepest of the two minima, which depends on the sign of $b$, the coefficient of the only odd term in $V$. For $b = 0$, the potential is an even function of $\rho$ and the two minima are equivalent (Fig. 5d): this criterion defines the so-called Maxwell set, which is shown as a dashed line in Fig. 5b. This second critical set is profoundly different from $\gamma_{th}^{\pm}$, as it can only be highlighted by tracking the potential functional $V(\rho; a, b)$ across yielding. When yielding occurs at $\gamma_{th}^{(M)}$, $V$ evolves continuously, as the two minima are equally deep, when $\rho$ jumps from one to the other. This is what characterizes first-order phase transitions, in which the free energy is continuous, although its derivative is not. On the contrary, when yielding
<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) Cusp catastrophe manifold and trajectories. Three trajectories ($K_r = 0.915, 1.000, 1.428$) are shown and represent ascending strain sweeps $\Delta\gamma_0 > 0$ for a glassy (blue), critical (purple) and liquid (red). The trajectory of the glassy system undergoes a discontinuous jump as it meets the critical set $\Sigma$. (b): Bifurcation set (black solid line) and three representative trajectories $a(b)$ for $K_r > 1$ (red, liquid), $K_r = 1$ (purple, critical iso-K) and $K_r < 1$ (blue, glass) in the control parameter space $(a,b)$ defined by eqs (11.b) and (11.c). The cusp point coincides to the critical point $(\gamma_r, \Gamma_r, K_r)=(1, 1, 1)$. The two intersections between the curve at $K_r < 1$ and the bifurcation set correspond to the yield points in the absence of coupling disorder for ascending (right) and descending (left) strain ramps. (c-e) Potential $V(\rho)$ for $K_r = 0.88$, $n = 3$ and the three critical strains corresponding to the critical set $\Sigma$ and zero-order transition (c and e), with a jump occurring from the vanishing minimum to the other (indicated by the arrows), and the Maxwell set (d) where the transition is reversible and first-order.</figcaption></figure>

occurs at $\gamma_{th}^{+}$, $V$ jumps discontinuously to the deeper minimum as its second local minimum disappears. This defines a different class of transitions, of the zeroth order.

In the case of thermodynamic phase transitions, the Maxwell set corresponds to the phase boundary. Here, $V(\rho;a,b)$ was introduced by formal analogy to thermodynamic free energy: although by construction its local minimization yields the equation of state, its global minimization does not have a straightforward physical relevance. Yet, surprisingly, simulations indicate that the Maxwell set, which we can write from Eq. 11.c as $\gamma_{th}^{(M)} = \gamma_c(4K_r-3)^{-1/n}$, does correspond to the yield point in the limit of large-scale simulations and for $d \to 0$. This result suggests that $V$ has an unanticipated

physical relevance, which is revealed by a finite, albeit small, disorder. The profound physical meaning of this functional is yet to be unveiled, and goes beyond the scope of this paper.

Here, we emphasize that disorder, which was initially introduced to account for the gradual onset of yielding as observed in experiments, plays a key role in determining the nature of the yielding transition itself. According to our FAd model, yielding is thus a rounded and continuous first-order transition with an emergent criticality at the yield strain, sharing features with both pure first- and second-order transitions. The yield point is set by a Maxwell rule, describing the strain at which solid-like dynamics become metastable and fluid-like domains can nucleate at weakly-coupled lattice sites, playing the role of soft spots [144]. The density of such soft spots sets a characteristic length scale that diverges at $d = 0$, where yielding becomes discontinuous.

## VII. EMERGING DYNAMIC LENGTH SCALE

This characteristic length scale emerges spontaneously from our on-lattice simulations, and dictates whether yielding manifests as an abrupt or a continuous transition in finite systems. To characterize this length scale, we compute the connected spatial correlation function of local rates: $c(r) \propto \langle\log\Gamma_i\log\Gamma_j\rangle - \langle\log\Gamma_i\rangle\langle\log\Gamma_j\rangle$, by averaging on all pairs of lattice sites at a given Euclidean distance $r$. We normalize $c(r)$, such that all correlations decay from 1 to 0, and we define the correlation length as its integral average $\xi = \int c(r)dr$. Operationally, we compute the integral by fitting $c(r)$ to Kohlrausch-Williams-Watts (KWW) functions $c(r) = \exp[-(\lambda r)^p]$, and expressing the integral as a function of the fit parameters $\xi = \lambda^{-1}\Gamma_e(p^{-1}+1)$, where $\Gamma_e(x)$ denotes Euler's Gamma function. Figure 6a shows representative $c(r)$ for $K_r = 0.86$ and $d = 0.03$ for strains spanning the coexistence region. For simulations exhibiting abrupt yielding, we find that $c(r)$ decays very fast for all values of $\gamma_0$: $\xi$ is as small as a single lattice spacing, reflecting the lack of spatial correlations in the coupling constants. By contrast, for simulations exhibiting gradual yielding, $\xi$ depends strongly on $\gamma_0$, and is significantly larger than the lattice spacing in the coexistence region. It peaks to a maximum value, $\xi^*$, quantifying the growth and coalescence of fluid domains, as shown in Figure 6b. We find that $\xi^*$ increases with decreasing $d$, reflecting the larger size of slow and fast domains coexisting at the transition, as shown previously in the snapshots of Fig. 3. Moreover, we observe an increase of the same peak height with decreasing reduced temperature $K_r$ as shown in figure 6c, that is coherent with the experimental evidence that the more a system is quenched in the glassy phase, the more it will be characterized by long-ranged correlated dynamics and sharp yielding.

This increase of $\xi^*$ with decreasing $d$ requires increasingly large simulation lattices to be properly captured. To find a more direct connection between $\xi^*$ and $N$, we invert $d^*(N)$ from Figure 4, to obtain a threshold lattice size for a given disorder, $N^*(d)$. We observe that the growth of $N^*$ with decreasing $d$ follows closely that of the largest lengthscale, $\xi^*$.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Spatial heterogeneities (a) symbols: spatial correlation of local rates for $K_r = 0.86$ and $d = 0.03$, for $\gamma_r$ across the yielding transition: $\gamma_r = 2.15$ (black), 2.25 (red), 2.41 (green), 2.53 (purple), 2.59 (gold). Lines: stretched exponential fits. (b) characteristic correlation length, $\xi$, as a function strain amplitudes in the coexistence region, for disorder $d$ increasing from 0.02 (blue) to 0.2 (red). (c) Symbols: maximum correlation length, $\xi^*$, as a function of disorder for $K_r$ decreasing from 0.98 (purple) to 0.86 (black), for $N = 512$. Full data-points are averaged on 10 independent simulations, all of which exhibited gradual yielding. Open black symbols denote disorders for which one or more simulations exhibited abrupt yielding: these simulations were excluded from the averaged dataset. Lines and crosses: threshold simulation lattice size, $N^*$, extracted for each disorder by inverting $d^*(N)$ data from Figure 4c.</figcaption></figure>

and that this holds for all the investigated temperatures, as shown by solid lines in Fig. 6c. We thus conclude that the transition from abrupt yielding to gradual yielding is dictated by the competition of two lengthscales: the system size, $N$, and the characteristic lengthscale of shear-induced dynamic heterogeneities, $\xi$. For $N > \xi^*$, the system yields gradually through the nucleation and growth of fast-relaxing domains, whereas for smaller $N$ the scarcity of weakly-coupled lattice sites delays the transition to larger strain amplitudes, where the system is metastable, and eventually yields abruptly.

This change emerges in our simulations as a finite size effect, rather than an intrinsic feature of our model itself, which predicts gradual yielding, independent of the sample size, for any $d > 0$. However, this effect highlights that a diverging correlation length may have tangible consequences on the manifestation of yielding at the macro scale. Exploiting our model to investigate the conditions under which a soft glassy sample may nucleate system-spanning heterogeneities such as shear bands animates our current research effort.

VIII. CONCLUSIONS

In this work we have presented a new model reproducing the key features of microscopic yielding in glassy systems under oscillatory shear, as reported by recent rheo-optical experiments [46]. We model the microscopic relaxation rate as a function of the imposed strain amplitude, by accounting for both the ultraslow, solid-like aging dynamics found in samples at rest or under small-amplitude deformations, and the faster, liquid-like dynamics measured stroboscopically in samples subject to large-amplitude oscillatory shear. We show that these two behaviors are captured by established models postulating a generic power-law dependence of shear-induced dynamics on the macroscopic shear rate, which however fail predicting the well-defined transition observed upon yielding. To properly describe this transition, we borrow ideas from dynamic facilitation of quiescent glasses, and we introduce a new dynamic coupling, which we called facilitated advection (FAd), between neighboring sample regions. The mean field solution of this model yields an equation of state predicting an abrupt transition between glassy and fluidized states as a function of strain amplitude. This equation of state shares the same cusp catastrophe manifold as Van der Waals gases, and predicts the existence of a critical point with a genuine second order transition for a set of model parameters corresponding to samples at the glass transition. We are currently working to test this prediction in experiments on colloidal suspensions at different concentrations, across the glass transition.

To describe gradual yielding as observed in experiments, we performed numerical simulations introducing disorder in the FAd coupling. We showed that a finite disorder changes qualitatively the nature of the transition: in presence of disorder, we find that solid and liquid-like domains coexist in a disorder-dependent range of strain amplitudes. Within this range, we measure a dynamic correlation length, $\xi$, that peaks at yielding and increases with decreasing $d$, diverging for $d \to 0$. The system size therefore plays a key role in determining the nature of the transition: for large enough systems, yielding is a rounded and continuous transition sharing features with both first-order and second-order transitions, akin to what has been observed in spin systems with coupling disorder [145]. Conversely, as $\xi$ exceeds the system size, we observe a qualitative change to abrupt yielding. This change emerges in our simulations as a finite size effect, rather than a feature of our model itself. However, it highlights that a diverging correlation length may have tangible consequences on the manifestation of yielding at the macro scale. We believe that these findings, which are entirely absent from [46], mark a substantial advancement in the understanding of the interplay between disorder and bifurcating dynamics in yielding phenomena, and help unify the observations of sharp yielding behavior in finite-size simulations [41] with a broader theoretical framework that accounts for disorder and spatial correlations.

To explore the relationship between local dynamics and rheological behavior, our ongoing effort is devoted to integrate our model with constitutive equations defining the time evolution of the stress within each mesoscopic portion of sheared
materials. This will allow us to test the generality of the FAd model to other rheological tests beyond oscillatory shear. In addition to that, we are investigating whether introducing terms including gradients of local rates ($\nabla\Gamma_i$) could lead to coalescence of fluidized (yielded) regions, giving rise potentially to steady shear-bands.

Finally, although our model has been developed to describe yielding in soft glasses, its formulation in terms of a strain-dependent local relaxation rate allows its further generalization to systems with more complex structure and rheology, such as gels, attractive glasses or anisotropic systems, in a broader attempt to achieve a unified description of yielding in soft amorphous materials.

## ACKNOWLEDGMENTS

This work was financially supported by CNES. SA acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-22-CE06-0025. We thank L. Cipelletti and W. Kob for illuminating discussions.

## Appendix A: Model with arbitrary range of facilitated advection coupling

Equation 4 of the main text takes a more general form when the dynamics of each site is coupled via bonds comprised within the $\beta-1$ shells centered on the site $i$. More in detail: The rate of each site $i$ is coupled to the all combinations of rates characterizing the nodes included in the shortest path going from the site $i$ to the site $j_{\beta-1}$, the latter being distant ($\beta-1$) lattice spacings from the former, as sketched in Fig. 7. Therefore the dynamics of each cell is coupled to $N_z(N_z-1)^{\beta-2}$ neighboring $j$-cells via positive coupling constants $\alpha_{i_{j_1}...j_{\beta-1}}$, where $2 \le \beta \le \infty$ dictates the range of the interaction. The most general form of this model then reads:

$$ \Gamma_i - \Gamma_0 = \left[ \frac{\omega K}{\frac{1}{\gamma_0^i} + \omega^\beta \frac{1}{N_z(N_z-1)^{\beta-2}} \sum_{j_1=1}^{N_z} \sum_{j_2=1}^{N_z-1} \cdots \sum_{j_{\beta-1}=1}^{N_z-1} \frac{\alpha_{i_{j_1}...j_{\beta-1}}}{\Gamma_i \prod_{k=1}^{\beta-1} \Gamma_{j_k}}} \right] \quad (A1) $$

where $\Gamma_{j_k} = 1/\tau_{j_k}$ is the relaxation rate characterizing the site $j_k$. The term containing the sums running over the $j_k$ indices represents the general facilitated advection term, modeling the average coupling between adjacent mesoscopic regions of soft glasses relaxing at different rates.

The mean field approximation of equation A1 is Equation 7 of the main text and it is obtained by imposing $\Gamma_i = \Gamma_{j_k} = \Gamma$ and $\alpha_{i_{j_1}...j_{\beta-1}} = \alpha$ for each site of index $i$. The impact of the coupling range in mean-field approximation is shown in Fig. 8.

<figure><img src="image_8.png" /><figcaption>FIG. 7. Sketch of the interaction between sites on a hexagonal lattice for $\beta = 4$. The site $i$ is coupled via the coupling constants $\alpha_{i_{j_1}j_2j_3}$ with all sites comprised in the first three shells centered on the site $i$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Effect of the spatial range of facilitated advection (FAd) on the mean field equation of state. a) Four selected equations of state with $\alpha = 0.17$, $K = 300$, $\omega = 0.66$ Hz and different values of $\beta$ as indicated in the figure. The values of the relaxation rate at the local minimum of the EOS $\Gamma(\gamma_{th}^+)$ are shown in the inset. Here $\beta$ is such we obtain only positive $\gamma_{th}^+$. Horizontal dashed lines show two representative jump of relaxation rate at the yielding transition. b) $\Gamma(\gamma_{th}^+)$ in function of $\beta$. Full points are the same shown in the inset of panel a) Empty points are obtained for larger values of $\beta$ for which $(\gamma_{th}^+)^{-n}$ is negative. The full horizontal line marks the limiting case $\Gamma(\gamma_{th}^+) = \Gamma_0$</figcaption></figure>

## Appendix B: Role of disorder: general considerations using the central limit theorem

As for classical first order transition the question we ask is whether the macroscopic phases survive in the presence of disorder or the system forms domains (droplets containing $N_D$ sites) that follow the local value of the random-interaction. To this aim we first make some general considerations on the effect of the disorder on the equation of state and consider a virtual single domain of one phase (say fluid) embedded in the other phase (say solid). We take equation A1 with $\beta = 2$ and
average over $N_D$ sites forming a fluid droplet:

$$ \frac{1}{N_D} \sum_{i=1}^{N_D} (\Gamma_i - \Gamma_0) = \frac{1}{N_D} \sum_{i=1}^{N_D} \left[ \frac{\omega K}{\left( \frac{1}{\gamma_0^\eta} + \frac{\omega^2}{N_z} \sum_{j=1}^{N_z} \frac{\alpha_{ij}}{\Gamma_i \Gamma_j} \right)} \right] = \frac{1}{N_D} \sum_{i=1}^{N_D} Y_i \quad (B1) $$

If $\alpha_{ij}$ is extracted from a distribution defined on a bounded support $[0, \alpha_M]$, then $Y = \{Y_i\}$ has also a pdf $P(Y)$ defined on a bounded support $\frac{\omega K}{(\gamma_0^n + (\frac{\omega}{\omega_0})^2 N_z \alpha_M)} \le Y_i \le K \omega \gamma_0^n$. Hence, the

variance of $P(Y)$ ($\sigma_Y^2$) exists and is finite if $P(Y)$ calculated at the two bounds is finite. For $N_D \gg 1$ we can apply the central limit theorem for weakly dependent random variables [135] to the expression (B1) and write

$$ \Gamma_D - \Gamma_0 = \frac{1}{N_D} [N_D E[Y] \pm \sigma_Y \sqrt{N_D}] = E[Y] \pm \frac{\sigma_Y}{\sqrt{N_D}} \sim E[Y] \pm \frac{\sigma_Y}{L^{d/2}} \quad (B2) $$

where $\Gamma_D$ is the average relaxation rate in the droplet, $E[Y]$ is the expectation value of $Y$ and $d$ is the dimension of the space. The expectation value $E[Y]$ is equal to the arithmetic mean in the limit $N_D \to \infty$.

Equation B2 can be viewed as a local MF equation for the droplet where now the effective spontaneous dynamics is characterized by a relaxation rate $(\Gamma_0 \pm \frac{\sigma_Y}{L^{d/2}})$ depending on the disorder. It is important to point out that the expectation value $E[Y]$ is not given by the mean field value obtained by fixing all coupling constants to one value $\alpha_{ij} = \alpha$, since it depends on the specific form of the probability density $P(Y)$ which in turn depends on the chosen pdf of the coupling constants $p(\alpha_{ij})$. Equation B2 shows that for $\sigma_Y = 0$ the value of the average relaxation rate of the droplet does not depend on its size: the dynamics coincides with that predicted by the mean field (MF) approximation. On the contrary, when $\sigma_Y > 0$ the system dynamics may split into fluid and solid droplets depending on whether locally $Y$ undergoes positive or negative fluctuations.

Equation B2, however, is not enough to inspect the effect of disorder and the occurrence of fluid nucleation and heterogeneous dynamics. For this purpose we re-write equation B1 with the random coupling constants $\alpha_{ij} = \langle \alpha \rangle + \epsilon_{ij}$, where $\epsilon_{ij}$ are random deviations with zero mean that here for simplicity we assume gaussian: $\epsilon_{ij} = N(0, \sigma_\alpha)$.

$$ \Gamma_D - \Gamma_0 = \frac{1}{N_D} \sum_{i=1}^{N_D} \left[ \frac{\omega K(\varphi)}{\left( \frac{1}{\gamma_0^\eta} + \frac{\omega^2}{N_z} \sum_{j=1}^{N_z} \frac{\langle \alpha \rangle + \epsilon_{ij}}{\Gamma_i \Gamma_j} \right)} \right] \quad (B3) $$

Since the system will split into fluid-like and solid-like droplets we further simplify the problem and we seek for a solution of equation B3 by adopting the approximation of homogeneous droplets, namely $\langle \Gamma_i \Gamma_j \rangle = \langle \Gamma \rangle^2 = \Gamma_D^2$. By doing

so and expanding for small $\epsilon_{ij}$, we obtain

$$ \Gamma_D - \Gamma_0 \approx \frac{K \omega}{\gamma_0^{-n} + \frac{\omega^2 \langle \alpha \rangle}{\Gamma_D^2}} \left[ 1 - \frac{1}{N_D} \sum_{i=1}^{N_D} \frac{\frac{\omega^2}{\Gamma_D^2 N_z} \sum_{j=1}^{N_z} \epsilon_{ij}}{\gamma_0^{-n} + \frac{\omega^2 \langle \alpha \rangle}{\Gamma_D^2}} \right] \quad (B4) $$

We can now apply the central limit theorem, knowing that the mean of the deviations over the nearest neighbor sites is still normally distributed with zero mean and standard deviation $\sigma_\alpha/\sqrt{N_z}$. Equation B4 becomes

$$ \Gamma_D - \Gamma_0 \approx \frac{K \omega}{\gamma_0^{-n} + \frac{\omega^2 \langle \alpha \rangle}{\Gamma_D^2}} \left[ 1 \pm \frac{\frac{\omega^2 \sigma_\alpha}{\sqrt{N_z N_D}}}{\Gamma_D^2 \gamma_0^{-n} + \frac{\omega^2 \langle \alpha \rangle}{\Gamma_D^2}} \right] \quad (B5) $$

From Equation B5 therefore we reach the following conclusions.

i) By increasing the number of nearest neighbors $N_z$ or the lattice sites $N_D$ belonging to a droplet the effect of disorder weakens. The mean field equation becomes thus exact in presence of disorder in infinite dimensions.

ii) In the presence of disorder ($\sigma_\alpha > 0$) the equation of state differs from the meanfield one and the system splits in solid and fluid regions, corresponding respectively to positive and negative local fluctuations of the coupling constants $\alpha_{ij}$ over droplets of size $N_D$. Fluid nucleation occurs for strain amplitudes $\gamma_0$ lower than $\gamma_{th}^+$, correspondent to the local minimum of the equation of state obtained for $\sigma_\alpha = 0$. Similarly, for $\gamma_0 > \gamma_{th}^+$ we will still have solid regions, since part of the system will follow equation B5 with the negative sign in its r.h.s. This simplifying approach, although it is nothing but a two-state model for nucleation as it considers only the average fluctuations via the central limit theorem, it allows to draw some important conclusions more quantitatively. Figure 9 shows different equations of state for zero disorder and both positive and negative fluctuations of the coupling constants in the droplets. By increasing progressively the strain we move downwards along the solid branch of the curves and reach the strain $\gamma_0 = \gamma_0^{(1)}$. In absence of disorder, this point would correspond to a homogeneous solid-like state, since the system did not reach the local minimum. In the presence of quenched disorder however small droplets can nucleate and reach, on average, a typical size equal to $N_D = 20$ for the case shown in figure 9, while other regions will stay solid with an average dynamics slower than that predicted by the MF equation. On the other hand some regions of the system will stay solid even beyond the yield strain predicted for $\sigma_\alpha = 0$, since they have not yet reached the critical set, namely the minimum of the (blue) mean field curves in figure 9. This critical strain, we repeat for clarity, depends on the droplet size. Such a condition, for example, occurs at $\gamma_0 = \gamma_0^{(2)}$, where only solid regions smaller than $N_D = 20$ survive. This two-state approximation therefore confirms that i) a yielding transition starts at lower strain with respect to the MF prediction when disorder is taken into account, and ii) that solid-like domains persist well beyond the yielding threshold that the system would experience in absence of disorder. This scenario is confirmed
<figure><img src="image_10.png" /><figcaption>FIG. 9. Equation of state (Eq. B5) in the presence of finite disorder with positive and negative average fluctuations of coupling constants in droplets of size $N_D$ as indicated in the figure for a representative set of parameters.</figcaption></figure>

by simulations discussed in section (V) and (VII) of the main manuscript.

## Appendix C: Convergence kinetics

Figure 10 shows the evolution of the average relaxation rate $\langle\Gamma\rangle$, the fraction of fast-relaxing sites $\chi$ and the value of the loss function $H$, as a function of the number of iterations. Figure 11 shows the number of iterations needed to reach convergence as a function of strain amplitude, $\gamma_0$.

## Appendix D: Statistics at fixed disorder (Figure 4 main text)

Figure 12 reports the analysis of multiple independent $N \times N$ simulations at a fixed disorder, to quantify the variability and assess finite size effects.

<figure><img src="image_11.png" /><figcaption>FIG. 10. Left column (a-c): evolution of the average relaxation rate $\langle\Gamma\rangle$ (a), the fraction of fast-relaxing sites $1-\chi$ (b) and the value of the loss function $\mathcal{L}$ (c), as a function of the number of iterations, for different values of the strain amplitude, increasing from blue to red shades. Black circles denote the convergence to the numerical solution. Right column (d-f): magnified view of the last iterations for the first strain amplitude above yielding. The probed strain amplitudes range from $10^{-2}$ up to 100, $\tau_0 = 10^4$s, $K=1$, $\beta=2$, $n=2$, $\alpha = 1.26 \cdot 10^{-4}$. The resulting yield strain is $(5.43536 \pm 0.00001) \cdot 10^{-2}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11. (a) Black lines, left axis: average relaxation rate as a function of $\gamma_0$. Solid, dashed and dotted lines correspond to three system sizes, $N = 10^2, 10^3, 10^4$, and $d = 10^{-9}$. Blue to red symbols, right axis: iterations required to reach convergence. Symbol color changes with $\gamma_0$, and correspond to line colors in Figure 10 (b). To highlight the power-law divergence at yielding, the same data from panel a are plotted in logarithmic scale, as a function of the distance from yielding: $\gamma_y/\gamma - 1$ and $\gamma/\gamma_y - 1$ for data below and above yielding, respectively.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 12. (a) Average rate, $\Gamma$, as a function of $\gamma_0$, for a selection of independent realizations of $N \times N$ simulation lattices, with $N = 8, 32, 128$ for black, blue and purple lines, respectively. (b) Solid lines: average on many independent simulations for different $N$, as specified in the legend. Shaded areas represent the dispersion of individual datapoints. (c) Histograms of threshold strains, $\gamma_{th}$, defined at $\chi = 0.5$, for many independent $N \times N$ simulations, with $N$ increasing from 8 to 512 as specified in the label. Dash-dotted line: $\gamma_v^{(M)}$. (d) average threshold strain, $\langle\gamma_{th}\rangle$, as a function of simulation lattice size, $N$, for various $d$ increasing from $10^{-5}$ to 0.15, as specified in the label. Error bars represent the standard variation of the threshold strain distribution. Dotted black line: mean-field threshold strain. Dash-dotted line: $\gamma_v^{(M)}$. (e) Red solid line: $\langle\gamma_{th}\rangle$ measured with $N = 512$ as a function of $d$. Red datapoints: threshold strain measured at $\chi = 0.1$. Open symbols: simulations with abrupt yielding. Red shaded area: span between threshold strain defined at $\chi = 0.1$ and $\chi = 0.9$. Dotted black line: $\gamma_{th}^+$, dash-dotted black line: $\gamma_v^{(M)}$. Blue line, symbols and shaded area: threshold strain measured in a descending strain amplitude ramp. Dashed black line: mean-field prediction for the threshold strain in a descending ramp.</figcaption></figure>
[1] A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Deformation and flow of amorphous solids: Insights from elasto-plastic models, Rev. Mod. Phys. **90**, 045006 (2018).

[2] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield Stress Materials in Soft Condensed Matter, Rev. Mod. Phys. **89**, 035005 (2017), arXiv:1502.05281.

[3] J. T. Parley, S. Sastry, and P. Sollich, Mean-Field Theory of Yielding under Oscillatory Shear, Phys. Rev. Lett. **128**, 198001 (2022).

[4] J. M. Brader, M. Siebenbürger, M. Ballauff, K. Reinheimer, M. Wilhelm, S. J. Frey, F. Weysser, and M. Fuchs, Nonlinear response of dense colloidal suspensions under oscillatory shear: Mode-coupling theory and Fourier transform rheology experiments, Phys. Rev. E **82**, 061401 (2010).

[5] T. Voigtmann, Nonlinear glassy rheology, Current Opinion in Colloid & Interface Science **19**, 549 (2014).

[6] C. Liu, K. Martens, and J.-L. Barrat, Mean-Field Scenario for the Athermal Creep Dynamics of Yield-Stress Fluids, Phys. Rev. Lett. **120**, 028004 (2018).

[7] P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Rheology of soft glassy materials, Phys Rev. Lett. **78**, 2020 (1997).

[8] G. Picard, A. Ajdari, L. Bocquet, and F. Lequeux, Simple model for heterogeneous flows of yield stress fluids, Physical Review E **66**, 10.1103/PhysRevE.66.051501 (2002).

[9] R. Benzi, T. Divoux, C. Barentin, S. Manneville, M. Sbragaglia, and F. Toschi, Unified Theoretical and Experimental View on Transient Shear Banding, Phys. Rev. Lett. **123**, 248001 (2019).

[10] T. Divoux, E. Agoritsas, S. Aime, C. Barentin, J.-L. Barrat, R. Benzi, L. Berthier, D. Bi, G. Biroli, D. Bonn, P. Bourrianne, M. Bouzid, E. Del Gado, H. Delanoë-Ayari, K. Farain, S. M. Fielding, M. Fuchs, J. Van Der Gucht, S. Henkes, M. Jalaal, Y. M. Joshi, A. Lemaître, R. L. Leheny, S. Manneville, K. Martens, W. Poon, M. Popovic, I. Procaccia, L. Ramos, J. A. Richards, S. A. Rogers, S. Rossi, M. Sbragaglia, G. Tarjus, F. Toschi, V. Trappe, J. Vermant, M. Wyart, F. Zamponi, and D. Zare, Ductile-to-brittle transition and yielding in soft amorphous materials: Perspectives and open questions, Soft Matter **20**, 6868-6888 (2024).

[11] G. J. Donley, P. K. Singh, A. Shetty, and S. A. Rogers, Elucidating the $G''$ overshoot in soft materials with a yield transition via a time-resolved experimental strain decomposition, Proc Natl Acad Sci U.S.A. **117**, 21945 (2020).

[12] K. Kamani, G. J. Donley, and S. A. Rogers, Unification of the Rheological Physics of Yield Stress Fluids, Phys. Rev. Lett. **126**, 218002 (2021).

[13] D. J. Pine, J. P. Gollub, J. F. Brady, and A. M. Leshansky, Chaos and threshold for irreversibility in sheared suspensions, Nature **438**, 997 (2005).

[14] K. Hima Nagamanasa, S. Gokhale, A. K. Sood, and R. Ganapathy, Experimental signatures of a nonequilibrium phase transition governing the yielding of a soft glass, Phys. Rev. E **89**, 062308 (2014).

[15] E. D. Knowlton, D. J. Pine, and L. Cipelletti, A microscopic view of the yielding transition in concentrated emulsions, Soft Matter **10**, 6931 (2014), arXiv:1403.4433.

[16] A. D. Gopal and D. J. Durian, Nonlinear Bubble Dynamics in a Slowly Driven Foam, Phys. Rev. Lett. **75**, 2610 (1995).

[17] D. Fiocco, G. Foffi, and S. Sastry, Oscillatory athermal quasistatic deformation of a model glass, Phys. Rev. E **88**, 020301 (2013).

[18] T. Kawasaki and L. Berthier, Macroscopic yielding in jammed solids is accompanied by a nonequilibrium first-order transition in particle trajectories, Phys. Rev. E **94**, 022615 (2016).

[19] I. Regev, T. Lookman, and C. Reichhardt, Onset of irreversibility and chaos in amorphous solids under periodic shear, Phys. Rev. E **88**, 062401 (2013).

[20] H. Bhaumik, G. Foffi, and S. Sastry, The role of annealing in determining the yielding behavior of glasses under cyclic shear deformation, Proc Natl Acad Sci USA **118**, e2100227118 (2021).

[21] C. Ness and M. E. Cates, Absorbing-State Transitions in Granular Materials Close to Jamming, Phys. Rev. Lett. **124**, 088004 (2020).

[22] R. Mari, E. Bertin, and C. Nardini, Absorbing phase transitions in systems with mediated interactions, Phys. Rev. E **105**, L032602 (2022).

[23] M. Mungan and T. A. Witten, Cyclic annealing as an iterated random map, Phys. Rev. E **99**, 052132 (2019).

[24] H. Hinrichsen, Non-equilibrium critical phenomena and phase transitions into absorbing states, Advances in Physics **49**, 815 (2000).

[25] S. Lübeck, Universal scaling behavior of non-equilibrium phase transitions, Int. J. Mod. Phys. B **18**, 3977 (2004).

[26] R. Dickman, A. Vespignani, and S. Zapperi, Self-organized criticality as an absorbing-state phase transition, Phys. Rev. E **57**, 5095 (1998).

[27] R. Maire, A. Plati, M. Stockinger, E. Trizac, F. Smallenburg, and G. Foffi, Interplay between an Absorbing Phase Transition and Synchronization in a Driven Granular System, Phys. Rev. Lett. **132**, 238202 (2024).

[28] R. M. Ziff, E. Gulari, and Y. Barshad, Kinetic Phase Transitions in an Irreversible Surface-Reaction Model, Phys. Rev. Lett. **56**, 2553 (1986).

[29] L. Corté, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Random organization in periodically driven systems, Nature Phys **4**, 420 (2008).

[30] R. Jeanneret and D. Bartolo, Geometrically protected reversibility in hydrodynamic Loschmidt-echo experiments, Nat. Commun. **5**, 3474 (2014).

[31] N. C. Keim and P. E. Arratia, Mechanical and Microscopic Properties of the Reversible Plastic Regime in a 2D Jammed Material, Phys. Rev. Lett. **112**, 028302 (2014).

[32] E. Tjhung and L. Berthier, Hyperuniform Density Fluctuations and Diverging Dynamic Correlations in Periodically Driven Colloidal Suspensions, Phys. Rev. Lett. **114**, 148301 (2015).

[33] J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo, Emergent Hyperuniformity in Periodically Driven Emulsions, Phys. Rev. Lett. **115**, 108301 (2015).

[34] K. A. Takeuchi, M. Kuroda, H. Chaté, and M. Sano, Directed Percolation Criticality in Turbulent Liquid Crystals, Phys. Rev. Lett. **99**, 234503 (2007).

[35] R. Möbius and C. Heussinger, (Ir)reversibility in dense granular systems driven by oscillating forces, Soft Matter **10**, 4806.

[36] B. Néel, I. Rondini, A. Turzillo, N. Mujica, and R. Soto, Dynamics of a first-order transition to an absorbing state, Phys. Rev. E **89**, 042206 (2014).

[37] M. Chantry, L. S. Tuckerman, and D. Barkley, Universal continuous transition to turbulence in a planar shear flow, J. Fluid Mech. **824**, R1 (2017).

[38] M. C. Rogers, K. Chen, M. J. Pagenkopp, T. G. Mason, S. Narayanan, J. L. Harden, and R. L. Leheny, Microscopic
signatures of yielding in concentrated nanoemulsions under large-amplitude oscillatory shear, Physical Review Materials **2**, 095601 (2018).

[39] P. K. Jaiswal, I. Procaccia, C. Rainone, and M. Singh, Mechanical Yield in Amorphous Solids: A First-Order Phase Transition, Phys. Rev. Lett. **116**, 085501 (2016).

[40] S. Karmakar, E. Lerner, and I. Procaccia, Statistical physics of the yielding transition in amorphous solids, Phys. Rev. E **82**, 055103 (2010).

[41] P. Leishangthem, A. D. S. Parmar, and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nat Commun **8**, 14653 (2017).

[42] T. Divoux, V. Grenard, and S. Manneville, Rheological Hysteresis in Soft Glassy Materials, Phys. Rev. Lett. **110**, 018304 (2013).

[43] M. Mungan, S. Sastry, K. Dahmen, and I. Regev, Networks and Hierarchies: How Amorphous Materials Learn to Remember, Phys. Rev. Lett. **123**, 178002.

[44] N. C. Keim and P. E. Arratia, Yielding and microstructure in a 2D jammed material under shear deformation, Soft Matter **9**, 6222 (2013).

[45] K. N. Nordstrom, J. P. Gollub, and D. J. Durian, Dynamical heterogeneity in soft-particle suspensions under shear, Phys. Rev. E **84**, 021403 (2011).

[46] S. Aime, D. Truzzolillo, D. J. Pine, L. Ramos, and L. Cipelletti, A unified state diagram for the yielding transition of soft colloids, Nat. Phys. 10.1038/s41567-023-02153-w (2023).

[47] Data of Ref [46] publicly available at: https://springernature.figshare.com/articles/dataset/A_unified_state_diagram_for_the_yielding_transition_ofSoft_colloids/22674334?file=40244296? (2023).

[48] J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, and L. Bocquet, Spatial cooperativity in soft glassy flows, Nature **454**, 84 (2008).

[49] C. Bonnoit, J. Lanuza, A. Lindner, and E. Clement, Mesoscopic Length Scale Controls the Rheology of Dense Suspensions, Phys. Rev. Lett. **105**, 108302 (2010).

[50] T. Sentjabrskaja, P. Chaudhuri, M. Hermes, W. C. K. Poon, J. Horbach, S. U. Egelhaaf, and M. Laurati, Creep and flow of glasses: Strain response linked to the spatial distribution of dynamical heterogeneities, Sci Rep **5**, 11884 (2015).

[51] S. Mukherji, N. Kandula, A. K. Sood, and R. Ganapathy, Strength of Mechanical Memories is Maximal at the Yield Point of a Soft Glass, Phys. Rev. Lett. **122**, 158001 (2019).

[52] P. Jop, V. Mansard, P. Chaudhuri, L. Bocquet, and A. Colin, Microscale Rheology of a Soft Glassy Material Close to Yielding, Phys. Rev. Lett. **108**, 148301 (2012).

[53] H. A. Vinutha, M. Marchand, M. Caggioni, V. V. Vasisht, E. Del Gado, and V. Trappe, Memory of shear flow in soft jammed materials, PNAS Nexus **3**, pgae441 (2024).

[54] K. Martens, L. Bocquet, and J.-L. Barrat, Spontaneous formation of permanent shear bands in a mesoscopic model of flowing disordered matter, Soft Matter **8**, 4197 (2012).

[55] N. C. Keim, J. Hass, B. Kroger, and D. Wieker, Global memory from local hysteresis in an amorphous solid, Phys. Rev. Res. **2**, 012004 (2020).

[56] H. Bhaumik, G. Foffi, and S. Sastry, Avalanches, Clusters, and Structural Change in Cyclically Sheared Silica Glass, Phys. Rev. Lett. **128**, 098001 (2022).

[57] H. Bhaumik, G. Foffi, and S. Sastry, Yielding transition of a two dimensional glass former under athermal cyclic shear deformation, J. Chem. Phys. **156**, 064502 (2022).

[58] E. Del Gado, P. Ilg, M. Kröger, and H. C. Öttinger, Non-

affine Deformation of Inherent Structure as a Static Signature of Cooperativity in Supercooled Liquids, Phys. Rev. Lett. **101**, 095501 (2008).

[59] M. L. Manning and A. J. Liu, Vibrational Modes Identify Soft Spots in a Sheared Disordered Packing, Phys. Rev. Lett. **107**, 108302 (2011).

[60] P. Schall, D. A. Weitz, and F. Spaepen, Structural Rearrangements That Govern Flow in Colloidal Glasses, Science **318**, 1895 (2007).

[61] R. L. Jack and L. Berthier, The melting of stable glasses is governed by nucleation-and-growth dynamics, The Journal of Chemical Physics **144**, 244506 (2016).

[62] K. L. Kearns, M. D. Ediger, H. Huth, and C. Schick, One Micrometer Length Scale Controls Kinetic Stability of Low-Energy Glasses, J. Phys. Chem. Lett. **1**, 388 (2010).

[63] S. F. Swallen, K. Traynor, R. J. McMahon, M. D. Ediger, and T. E. Mates, Stable Glass Transformation to Supercooled Liquid via Surface-Initiated Growth Front, Phys. Rev. Lett. **102**, 065503 (2009).

[64] C. Herrero, C. Scalliet, M. D. Ediger, and L. Berthier, Two-step devitrification of ultrastable glasses, Proc. Natl. Acad. Sci. U.S.A. **120**, e2220824120 (2023).

[65] A. S. Keys, L. O. Hedges, J. P. Garrahan, S. C. Glotzer, and D. Chandler, Excitations Are Localized and Relaxation Is Hierarchical in Glass-Forming Liquids, Phys. Rev. X **1**, 021013 (2011).

[66] M. D. Ediger, Spatially Heterogeneous Dynamics in Supercooled Liquids, Annu. Rev. Phys. Chem **51**, 99 (2000).

[67] E. Donth, H. Huth, and M. Beiner, Characteristic length of the glass transition, J. Phys.: Condens. Matter **13**, L451 (2001).

[68] P. Edera, M. Briziolı, M. Madani, E. Ngouamba, P. Cousset, V. Trappe, G. Petekidis, F. Giavazzi, and R. Cerbino, Yielding under the microscope: A multi-scale perspective on brittle and ductile behaviors in oscillatory shear (2024), arXiv:2402.00221 [cond-mat].

[69] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Random critical point separates brittle and ductile yielding transitions in amorphous materials, Proc. Natl. Acad. Sci. U.S.A. **115**, 6656 (2018).

[70] J. Pollard and S. M. Fielding, Yielding, shear banding, and brittle failure of amorphous materials, Phys. Rev. Research **4**, 043037 (2022).

[71] A. Pommella, L. Cipelletti, and L. Ramos, Role of Normal Stress in the Creep Dynamics and Failure of a Biopolymer Gel, Phys. Rev. Lett. **125**, 268006 (2020).

[72] A. Le Bouil, A. Amon, S. McNamara, and J. Crassous, Emergence of Cooperativity in Plasticity of Soft Glassy Materials, Phys. Rev. Lett. **112**, 246001 (2014).

[73] D. Houdoux, T. B. Nguyen, A. Amon, and J. Crassous, Plastic flow and localization in an amorphous material: Experimental interpretation of the fluidity, Phys. Rev. E **98**, 022905 (2018).

[74] A. S. Argon and H. Y. Kuo, Plastic flow in a disordered bubble raft (an analog of a metallic glass), Mater. Sci. Eng. **39**, 101 (1979).

[75] C. Maloney and A. Lemaître, Universal Breakdown of Elasticity at the Onset of Material Failure, Phys. Rev. Lett. **93**, 195501 (2004).

[76] R. Dasgupta, H. G. E. Hentschel, and I. Procaccia, Microscopic Mechanism of Shear Bands in Amorphous Solids, Phys. Rev. Lett. **109**, 255502 (2012).

[77] J. Antonaglia, W. J. Wright, X. Gu, R. R. Byer, T. C. Hufnagel, M. LeBlanc, J. T. Uhl, and K. A. Dahmen, Bulk Metallic Glasses Deform via Slip Avalanches, Phys. Rev. Lett. **112**, 155501 (2014).
[78] J. M. Alegre, I. I. Cuesta, and M. Lorenzo, An Extension of the Monkman-Grant Model for the Prediction of the Creep Rupture Time Using Small Punch Tests, Exp Mech **54**, 1441 (2014).

[79] J. Koivisto, M. Ovaska, A. Miksic, L. Laurson, and M. J. Alava, Predicting sample lifetimes in creep fracture of heterogeneous materials, Physical Review E **94**, 23002 (2016).

[80] A. Garcimartín, A. Guarino, L. Bellon, and S. Ciliberto, Statistical Properties of Fracture Precursors, Phys. Rev. Lett. **79**, 3202 (1997).

[81] N. Kandula, B. Cordonnier, E. Boller, J. Weiss, D. K. Dysthe, and F. Renard, Dynamics of Microscale Precursors During Brittle Compressive Failure in Carrara Marble, J. Geophys. Res. Solid Earth **124**, 6121 (2019).

[82] F. Renard, B. Cordonnier, M. Kobchenko, N. Kandula, J. Weiss, and W. Zhu, Microscale characterization of rupture nucleation unravels precursors to faulting in rocks, Earth Planet. Sci. Lett. **476**, 69 (2017).

[83] D. A. Lockner, J. D. Byerlee, V. Kuksenko, A. Ponomarev, and A. Sidorin, Quasi-static fault growth and shear fracture energy in granite, Nature **350**, 39 (1991).

[84] D. Lockner, The role of acoustic emission in the study of rock fracture, Int. J. Rock Mech. Min. Sci. **30**, 883 (1993).

[85] C. A. Schuh and A. C. Lund, Atomistic basis for the plastic yield criterion of metallic glass, Nature Mater **2**, 449 (2003).

[86] G. Zöller, S. Hainzl, and J. Kurths, Observation of growing correlation length as an indicator for critical point behavior prior to large earthquakes, J. Geophys. Res. **106**, 2167 (2001).

[87] S. M. Potirakis, A. Karadimitrakis, and K. Eftaxias, Natural time analysis of critical phenomena: The case of pre-fracture electromagnetic emissions, Chaos **23**, 023117 (2013).

[88] L. Girard, J. Weiss, and D. Amitrano, Damage-Cluster Distributions and Size Effect on Strength in Compressive Failure, Phys. Rev. Lett. **108**, 225502 (2012).

[89] D. Sornette, Predictability of catastrophic events: Material rupture, earthquakes, turbulence, financial crashes, and human birth, Proc. Natl. Acad. Sci. U.S.A. **99**, 2522 (2002).

[90] D. Sornette and J. Andersen, Scaling with respect to disorder in time-to-failure, Eur. Phys. J. B **1**, 353 (1998).

[91] S. Zapperi, A. Vespignani, and H. E. Stanley, Plasticity and avalanche behaviour in microfracturing phenomena, Nature **388**, 658 (1997).

[92] P. A. Johnson, B. Ferdowsi, B. M. Kaproth, M. Scuderi, M. Griffa, J. Carmeliet, R. A. Guyer, P.-Y. Le Bas, D. T. Trugman, and C. Marone, Acoustic emission and microslip precursors to stick-slip failure in sheared granular material, Geophys. Res. Lett. **40**, 5627 (2013).

[93] D. Houdoux, A. Amon, D. Marsan, J. Weiss, and J. Crassous, Micro-slips in an experimental granular shear band replicate the spatiotemporal characteristics of natural earthquakes, Commun Earth Environ **2**, 1 (2021).

[94] A. A. Ostapchuk and K. G. Morozova, On the Mechanism of Laboratory Earthquake Nucleation Highlighted by Acoustic Emission, Sci Rep **10**, 7245 (2020).

[95] S. Gvirtzman and J. Fineberg, Nucleation fronts ignite the interface rupture that initiates frictional motion, Nat. Phys. **17**, 1037 (2021).

[96] R. A. Kromer, D. J. Hutchinson, M. J. Lato, D. Gauthier, and T. Edwards, Identifying rock slope failure precursors using LiDAR for transportation corridor hazard management, Eng. Geol. **195**, 93 (2015).

[97] M. M. Scuderi, C. Marone, E. Tinti, G. Di Stefano, and C. Collettini, Precursory changes in seismic velocity for the spectrum of earthquake failure modes, Nature Geosci **9**, 695 (2016).

[98] T. H. Dixon, Y. Jiang, R. Malservisi, R. McCaffrey, N. Voss, M. Protti, and V. Gonzalez, Earthquake and tsunami forecasts: Relation of slow slip events to subsequent earthquake rupture, Proceedings of the National Academy of Sciences **111**, 17039 (2014).

[99] J. Vasseur, F. B. Wadsworth, Y. Lavallée, A. F. Bell, I. G. Main, and D. B. Dingwell, Heterogeneity: The key to failure forecasting, Sci Rep **5**, 13259 (2015).

[100] V. Kádár, G. Pál, and F. Kun, Record statistics of bursts signals the onset of acceleration towards failure, Sci. Rep. **10**, 2508 (2020).

[101] J. Yu, A. Datye, Z. Chen, C. Zhou, O. E. Dagdeviren, J. Schroers, and U. D. Schwarz, Atomic-scale homogeneous plastic flow beyond near-theoretical yield stress in a metallic glass, Commun Mater **2**, 1 (2021).

[102] B. A. Sun, H. B. Yu, W. Jiao, H. Y. Bai, D. Q. Zhao, and W. H. Wang, Plasticity of Ductile Metallic Glasses: A Self-Organized Critical State, Phys. Rev. Lett. **105**, 035501 (2010).

[103] H. Idrissi, M. Ghidelli, A. Béché, S. Turner, S. Gravier, J.-J. Blandin, J.-P. Raskin, D. Schryvers, and T. Pardoen, Atomic-scale viscoplasticity mechanisms revealed in high ductility metallic glass films, Sci Rep **9**, 13426 (2019).

[104] D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J.-L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner, and M. L. Manning, Predicting plasticity in disordered solids from structural indicators, Phys. Rev. Materials **4**, 113609 (2020).

[105] S. Rossi and G. Tarjus, Emergence of a random field at the yielding transition of a mean-field elasto-plastic model, J. Stat. Mech. **2022**, 093301 (2022).

[106] D. Rodney, A. Tanguy, and D. Vandembroucq, Modeling the mechanics of amorphous solids at different length scale and time scale, Modelling Simul. Mater. Sci. Eng. **19**, 083001 (2011).

[107] B. F. Di Dio, F. Khabaz, R. T. Bonnecaze, and M. Cloitre, Transient dynamics of soft particle glasses in startup shear flow. Part II: Memory and aging, J. Rheol. **66**, 717 (2022).

[108] P. Edera, M. Bantawa, S. Aime, R. T. Bonnecaze, and M. Cloitre, Tuning residual stress, directional memory and aging in soft glassy materials (2024), arXiv:2402.08293 [cond-mat].

[109] M. Vasoya, C. H. Rycroft, and E. Bouchbinder, Notch Fracture Toughness of Glasses: Dependence on Rate, Age, and Geometry, Phys. Rev. Appl. **6**, 024008 (2016).

[110] M. Fan, M. Wang, K. Zhang, Y. Liu, J. Schroers, M. D. Shat-tuck, and C. S. O'Hern, Effects of cooling rate on particle rearrangement statistics: Rapidly cooled glasses are more ductile and less reversible, Phys. Rev. E **95**, 022611 (2017).

[111] G. Kumar, P. Neibecker, Y. H. Liu, and J. Schroers, Critical fictive temperature for plasticity in metallic glasses, Nat. Commun. **4**, 1536 (2013).

[112] D. Şopu, F. Spieckermann, X. Bian, S. Fellner, J. Wright, M. Cordill, C. Gammer, G. Wang, M. Stoica, and J. Eckert, Rejuvenation engineering in metallic glasses by complementary stress and structure modulation, NPG Asia Mater **15**, 1 (2023).

[113] X. Bian, D. Şopu, G. Wang, B. Sun, J. Bednarčik, C. Gammer, Q. Zhai, and J. Eckert, Signature of local stress states in the deformation behavior of metallic glasses, NPG Asia Mater **12**, 1 (2020).

[114] S. Patinet, D. Vandembroucq, and M. L. Falk, Connecting Local Yield Stresses with Plastic Activity in Amorphous Solids, Phys. Rev. Lett. **117**, 045501 (2016).
[115] S. Rossi, G. Biroli, M. Ozawa, and G. Tarjus, Far-from-equilibrium criticality in the random-field Ising model with Eshelby interactions, Phys. Rev. B **108**, L220202 (2023).

[116] J. P. Sethna, K. A. Dahmen, and C. R. Myers, Crackling noise, Nature **410**, 242 (2001).

[117] S. Rossi, G. Biroli, M. Ozawa, G. Tarjus, and F. Zamponi, Finite-Disorder Critical Point in the Yielding Transition of Elastoplastic Models, Phys. Rev. Lett. **129**, 228002 (2022).

[118] A. D. S. Parmar, S. Kumar, and S. Sastry, Strain Localization Above the Yielding Point in Cyclically Deformed Glasses, Phys. Rev. X **9**, 021018 (2019).

[119] M. Ozawa, L. Berthier, G. Biroli, and G. Tarjus, Rare events and disorder control the brittle yielding of well-annealed amorphous solids, Phys. Rev. Res. **4**, 023227 (2022).

[120] B. J. Berne and R. Pecora, *Dynamic Light Scattering with Application to Chemistry, Biology, and Physics*. (John Wiley and Sons, INC, 1990).

[121] A.-M. Philippe, D. Truzzolillo, J. Galvan-Myoshi, P. Dieudonné-George, V. Trappe, L. Berthier, and L. Cipelletti, Glass transition of soft colloids, Phys. Rev. E **97**, 040601 (2018).

[122] H. M. Wyss, K. Miyazaki, J. Mattsson, Z. Hu, D. R. Reichman, and D. A. Weitz, Strain-Rate Frequency Superposition: A Rheological Probe of Structural Relaxation in Soft Materials, Phys. Rev. Lett. **98**, 238303 (2007).

[123] A. Hess and N. Aksel, Yielding and structural relaxation in soft materials: Evaluation of strain-rate frequency superposition data by the stress decomposition method, Phys. Rev. E **84**, 051502 (2011).

[124] J. E. Hallett, F. Turci, and C. P. Royall, Local structure in deeply supercooled liquids exhibits growing lengthscales and dynamical correlations, Nat Commun **9**, 3272 (2018).

[125] S. Gokhale, K. Hima Nagamanasa, R. Ganapathy, and A. K. Sood, Growing dynamical facilitation on approaching the random pinning colloidal glass transition, Nat Commun **5**, 4685 (2014).

[126] C. K. Mishra, K. Hima Nagamanasa, R. Ganapathy, A. K. Sood, and S. Gokhale, Dynamical facilitation governs glassy dynamics in suspensions of colloidal ellipsoids, Proc. Natl. Acad. Sci. U.S.A. **111**, 15362 (2014).

[127] D. Chandler and J. P. Garrahan, Dynamics on the Way to Forming Glass: Bubbles in Space-Time, Annu. Rev. Phys. Chem **61**, 191 (2010).

[128] T. Speck, Dynamic facilitation theory: A statistical mechanics approach to dynamic arrest, J. Stat. Mech. **2019**, 084015 (2019).

[129] L. D. Landau, E. M. Lifšic, L. P. Pitaevskij, and L. D. Landau, *Statistical Physics. 1: By E. M. Lifshitz and L. P. Pitaevskii*, 3rd ed., Course of Theoretical Physics / L. D. Landau and E. M. Lifshitz No. 5 (Elsevier Butterworth Heinemann, Amsterdam Heidelberg, 2011).

[130] H. Tong, S. Sengupta, and H. Tanaka, Emergent solidity of amorphous materials as a consequence of mechanical self-organisation, Nat Commun **11**, 4863 (2020).

[131] J. Colombo, A. Widmer-Cooper, and E. Del Gado, Microscopic Picture of Cooperative Processes in Restructuring Gel Networks, Phys. Rev. Lett. **110**, 198301 (2013).

[132] L. Ortlieb, T. S. Ingebrigtsen, J. E. Hallett, F. Turci, and C. P. Royall, Probing excitations and cooperatively rearranging regions in deeply supercooled liquids, Nat Commun **14**, 2621 (2023).

[133] E. D. Cubuk, R. J. S. Ivancic, S. S. Schoenholz, D. J. Strickland, A. Basu, Z. S. Davidson, J. Fontaine, J. L. Hor, Y.-R. Huang, Y. Jiang, N. C. Keim, K. D. Koshigan, J. A. Lefever, T. Liu, X.-G. Ma, D. J. Magagnosc, E. Morrow, C. P. Ortiz, J. M. Rieser, A. Shavit, T. Still, Y. Xu, Y. Zhang, K. N. Nordstrom, P. E. Arratia, R. W. Carpick, D. J. Durian, Z. Fakhraai, D. J. Jerolmack, D. Lee, J. Li, R. Riggleman, K. T. Turner, A. G. Yodh, D. S. Gianola, and A. J. Liu, Structure-property relationships from universal signatures of plasticity in disordered solids, Science **358**, 1033 (2017).

[134] H. Tong and H. Tanaka, Structural order as a genuine control parameter of dynamics in simple glass formers, Nat Commun **10**, 5596 (2019).

[135] M. Fleermann and W. Kirsch, The central limit theorem for weakly dependent random variables by the moment method, arXiv:2202.04717, 17 (2022).

[136] O. Petrov and I. Furó, A study of freezing–melting hysteresis of water in different porous materials. Part II: Surfactant-templated silicas, Phys. Chem. Chem. Phys. **13**, 16358 (2011).

[137] S. Chikazumi, C. D. Graham, and S. Chikazumi, *Physics of Ferromagnetism*, 2nd ed., The International Series of Monographs on Physics No. 94 (Clarendon Press ; Oxford University Press, Oxford : New York : Warszawa, 1997).

[138] M. T. Dang, D. Denisov, B. Struth, A. Zaccone, and P. Schall, Reversibility and hysteresis of the sharp yielding transition of a colloidal glass under oscillatory shear, Eur. Phys. J. E **39**, 44 (2016).

[139] W. Weibull, *A statistical theory of the strength of materials*, (Generalstabens litografiska anstälts förlag, Stockholm, 1939).

[140] B. D. Coleman, Statistics and Time Dependence of Mechanical Breakdown in Fibers, J. Appl. Phys. **29**, 968 (1958).

[141] S. Aime, N. D. Eisenmenger, and T. A. P. Engels, A model for failure in thermoplastic elastomers based on Eyring kinetics and network connectivity, J. Rheol. **61**, 1329 (2017).

[142] A. Okniński, *Catastrophe Theory*, Comprehensive Chemical Kinetics No. v. 33 (Elsevier ; PWN–Polish Scientific Publishers, Amsterdam ; New York : Warszawa, 1992).

[143] M. Demazure, *Bifurcations and Catastrophes* (Springer, Berlin, 2000).

[144] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, Irreversible reorganization in a supercooled liquid originates from localized soft modes, Nature Phys **4**, 711 (2008).

[145] A. Bellafard, S. Chakravarty, M. Troyer, and H. G. Katzgraber, The effect of quenched bond disorder on first-order phase transitions, Annals of Physics **357**, 66 (2015).