Scaling Description of Dynamical Heterogeneity and Avalanches of Relaxation in Glass-Forming Liquids

Ali Tahaei,¹ Giulio Biroli,² Misaki Ozawa,³ Marko Popović,¹,⁴ and Matthieu Wyart⁵

¹Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str.38, 01187 Dresden, Germany

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

³Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France

⁴Center for Systems Biology Dresden, Pfotenhauerstrasse 108, 01307 Dresden, Germany

⁵Institute of Physics, EPFL, Lausanne, Switzerland

We provide a theoretical description of dynamical heterogeneities in glass-forming liquids, based on the premise that relaxation occurs via local rearrangements coupled by elasticity. In our framework, the growth of the dynamical correlation length $\xi$ and of the correlation volume $\chi_4$ are controlled by a zero-temperature fixed point. We connect this critical behavior to the properties of the distribution of local energy barriers at zero temperature. Our description makes a direct connection between dynamical heterogeneities and avalanche-type relaxation associated to dynamic facilitation, allowing us to relate the size distribution of heterogeneities to their time evolution. Within an avalanche, a local region relaxes multiple times, the more the larger is the avalanche. This property, related to the nature of the zero-temperature fixed point, directly leads to decoupling of particle diffusion and relaxation time (the so-called Stokes-Einstein violation). Our most salient predictions are tested and confirmed by numerical simulations of scalar and tensorial thermal elasto-plastic models.

# I. INTRODUCTION

As a liquid is cooled, the relaxation time $\tau_\alpha$ below which it acts as a solid – before displaying flow – grows from picoseconds at high temperatures up to minutes at the glass transition temperature $T_g$ [1–4]. In this regime the effective activation energy associated to $\tau_\alpha$ grows for many liquids, leading to a super-Arrhenius behavior. Approaching $T_g$, dynamics also becomes heterogeneous on a growing correlation length scale $\xi$ [5–8]. The underlying causes for these observations are still debated. In some views, activation is cooperative: the slowing down of the dynamics is governed by complex motion taking place on an increasingly large static length-scale. In particular, cooperativity is central in the Random First Order Theory [9–11] and due to the growth of amorphous order. Another approach focuses on dynamical facilitation, the phenomenon by which a region's relaxation is made much more likely by a relaxation nearby. In kinetically constrained models, such as the East model, kinetic rules induce dynamic facilitation and growth of dynamical correlations [12, 13]. This lead to a theory [14–16] in which thermodynamics plays almost no role, but dynamics is heterogeneous (due to kinetic constraints) and the super-Arrhenius behavior is due to non-local rearrangements taking place over $\xi$. Free volume [17] or elastic [18–20] models assume that the activation energy is not controlled by a static growing length scale: it is governed by the energy barrier of elementary rearrangements, or excitations. Recently, measurements indicated that the distribution of local energy barriers shifts to higher energy under cooling, opening up a gap at low energies where excitations are nearly absent [21]. This shift accounts quantitatively for the dynamical slowing down of the liquid, supporting that local energy barriers may indeed control the dynamics. Yet in these views, what causes

the existence of a growing length scale is unclear. An intuitive resolution of this paradox stems from the fact that on relatively short time scales, a super-cooled liquid acts as a solid [18, 22]: a rearrangement corresponds to a local (plastic) strain, that affects stress away from it. Lemaitre [23] was the first in stressing the possible relevance of long-range stress correlations in the dynamics of super-cooled liquids, in line with various theoretical and numerical studies [24–30]. Recently, it was shown that indeed elastic interactions play an important role in dynamical facilitation [31]. Moreover, Ref. [32] showed that heterogeneous relaxations take place in the form of plastic rearrangements that are called shear transformations [33]. Following these molecular simulations studies, Ref. [34] demonstrated with elasto-plastic models (traditionally used to study the plasticity of amorphous solids under loading [35–39]) that while being controlled by local energy barriers, the dynamics can at the same time display growing dynamical correlations similar to observations in experiments and molecular simulations [13, 40]. Furthermore, such models also capture the emergence of a gap in the distribution of local barriers under cooling, which controls the dynamics [34, 41]. We expect that the elasto-plastic description discussed above becomes more and more relevant with decreasing temperature. This paper focuses on such a lower temperature regime.

In this work we propose a scaling description of dynamical heterogeneities in glass-forming liquids, modeled as undergoing local irreversible rearrangements coupled by elasticity. We show that the dynamical correlations observed in elasto-plastic models of equilibrium glassy dynamics [34] are due to "thermal avalanches", where rare nucleated events are followed by a pulse of faster (or facilitated) events. These avalanches are very reminiscent of the ones found in molecular simulations of
super-cooled liquids [42–44]. Our analysis builds a link with systems that crackle such as disordered magnets, granular materials, and earthquakes [45, 46] even at finite temperatures [41, 47–49], revealing that dynamical heterogeneities are controlled by a critical point at zero temperature, with a diverging length scale $\xi \sim T^{-\nu}$ and correlation volume $\chi_4^* \sim T^{-\gamma}$. We provide a scaling argument expressing $\nu$ and $\gamma$ in terms of the distribution of energy barriers at $T = 0$. Our analysis makes predictions on the power-law distribution of the size of dynamical heterogeneities, which could be tested in numerical simulations of super-cooled liquids thanks to the recent advances in the characterization of dynamical correlations [50, 51]. Based on the properties of the zero-temperature fixed point governing thermal avalanches, we also show that the decoupling $D\tau_\alpha$ between particle diffusion $D$ and relaxation time $\tau_\alpha$ (the so-called Stoke-Einstein violation—[52–56]) diverges as $D\tau_\alpha \sim T^{-h}$, where $h$ can be expressed in terms of avalanche properties at vanishing temperature. The key physical mechanism behind the Stoke-Einstein violation is the intensive accumulation of rearrangements in the mobile region that quantifies a larger diffusion constant $D$, relative to the structural relaxation time $\tau_\alpha$ [57–61]. We perform numerical simulations and analyze the results by finite size scaling in order to test our salient predictions. Our findings agree well with the scaling theory both for scalar and tensorial thermal elasto-plastic models.

## II. THERMAL ELASTO-PLASTIC MODELS

To study dynamical heterogeneities in low-temperature glass-forming liquids, we employ a generalization of elasto-plastic models [35, 36] which are a class of mesoscopic models designed to capture the essential features of localized plastic events (shear transformations) with stress relaxation accompanied by long-range elastic interactions. An elasto-plastic model contains $N = L^d$ mesoscopic sites that are arranged on a regular grid whose linear size is $L$, and $d$ is the spacial dimension. Each site exhibits a plastic event when the magnitude of the local shear stress becomes sufficiently large, which leads to local stress relaxation and stress redistribution in the rest of the system via the form of Eshelby fields. Elasto-plastic models were originally introduced to study the flow of amorphous solids under external loading [35, 36, 38, 39] (some exceptions e.g., Ref. [62]) and they were generalized to take into account thermal fluctuations [63] and then to study glass-forming liquids [34]. Here, we study their low-temperature relaxation dynamics in the absence of loading. Following the idea that supercooled liquids can be viewed as solids that flow [18, 23], we develop a physical scenario based on the assumption that local rearrangements are elastically coupled, and we model this physical mechanism by elasto-plastic models.

In this paper, we use a tensorial elasto-plastic model

that accounts for the shear stress tensor [64] and a scalar elasto-plastic model in which the shear stress is represented by a scalar variable [34] (see Appendix A for details). In both models, thermal fluctuations are implemented through a probability rate $\tau_0^{-1}e^{-E/T}$ to trigger a plastic (mobile) event in an otherwise stable elastic (immobile) site, where $E$ is the local activation energy barrier and $\tau_0$ is a microscopic relaxation timescale [41, 63, 65, 66]. The energy barrier $E$ can be related to the minimal amount of additional shear stress $x_i$ required to destabilize a site $i$. In particular, we consider $E(x) = cx^\alpha$, which is suggested by recent elasto-plastic models and molecular simulations [32, 41, 65–67]. In this study, we set $\tau_0 = 1$, $c = 1$, and the value $\alpha = 3/2$ corresponds to a smooth microscopic potential [68]. As stress is relaxed locally at the site of the plastic event, stress in the system changes according to the Eshelby kernel [35]. Specifically, in the scalar model the kernel is randomly oriented at each relaxation event. This random orientation of the Eshelby kernel is a crucial to describe for isotropic supercooled liquids in a quiescent state, in contrast to amorphous solids under shear, where Eshelby fields align along the shear direction [69]. Note that in the tensorial model, after each relaxation the yielding surface is reoriented uniformly at random and one uses the full Eshelby propagator, hence no additional symmetrization is required [70]. Besides, our elastoplastic models have zero total (macroscopic) stress, as anticipated in quiescent supercooled liquids. Detailed description of the tensorial and scalar models are found in Appendix A. We perform numerical simulations in a two-dimensional periodic lattice ($d = 2$), whereas we develop scaling arguments in general $d$ dimensions. Since we find that both scalar and tensor models show qualitatively and quantitatively similar results (e.g., critical exponents), we focus on the tensorial model in the main text and present the scalar one in Appendix.

## III. DYNAMICAL HETEROGENEITIES AT FINITE TEMPERATURES

We first perform finite temperature simulations and characterize dynamical properties. To this end, we consider the persistence two-point time correlation function, $\langle\pi(t)\rangle_{\text{time}}$, which has been used widely in the context of kinetically constrained models [12]. The observable $\pi(t)$ is defined by $\pi(t) = \sum_i p_i(t)/N$, where $p_i(t) = 1$ if the site $i$ did not exhibit a plastic event until time $t$ and remained immobile, and $p_i(t) = 0$ for mobile sites that relaxed at least once. The notation $\langle\cdots\rangle_{\text{time}}$ denotes the time average at the stationary state at temperature $T$ reached after a long enough equilibration time. We have verified that the model does show dynamical slowing down by measuring $\langle\pi(t)\rangle_{\text{time}}$ at different temperatures, as shown in Fig. 1(a). In Fig. 1(b), we find that the associated relaxation time $\tau_\alpha$ defined by $\langle\pi(\tau_\alpha)\rangle_{\text{time}} = 1/2$ increases in a Arrhenius way. Figure 2(a) represents a
<figure><img src="image_1.png" /><figcaption>FIG. 1. Dynamics of the tensorial model in finite temperature simulations. (a): Mean persistence correlation function, $\langle\pi(t)\rangle_{\text{time}}$, for $L = 128$ and $T = 0.050, 0.040, 0.030, 0.025, 0.020, 0.018$, and 0.015 (from left to right). (b): The relaxation time $\tau_{\alpha}$ versus $1/T$. The red dashed-line corresponds to $\tau_{\alpha} \sim e^{E_c/T}$ with $E_c = x_c^{3/2} = 0.42$ measured independently in Sec. IV.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>FIG. 2. Dynamical heterogeneity of the tensorial model in finite temperature simulations. (a): Snapshot characterized by local persistence, $p_i(\tau)$, when $\tau/\tau_{\alpha} \approx 0.53$ at $T = 0.015$. Red and purple sites correspond to mobile ($p_i(\tau) = 0$) and immobile ($p_i(\tau) = 1$) sites, respectively. The system size is $L = 128$. (b): Four-point correlation function, $\chi_4(t)$, for $L = 128$ and $T = 0.050, 0.040, 0.030, 0.025, 0.020, 0.018$, and 0.015 (from left to right)..</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 3. (a): The peak of four-point correlation function, $\chi_4^*$, versus $T$ for several $L$ for the tensorial model. The red dashed-line follows $\chi_4^* \sim T^{-\gamma}$. (b): The corresponding scaling collapse.</figcaption></figure>

snapshot of local persistence for $\tau/\tau_{\alpha} \approx 0.53$, demonstrating spatially heterogeneous dynamics. To quantify the magnitude of dynamical heterogeneity, we compute the four-point correlations function [71], $\chi_4(t)$, defined by

$$\chi_4(t) = N (\langle \pi^2(t) \rangle_{\text{time}} - \langle \pi(t) \rangle_{\text{time}}^2). \quad (1)$$

$\chi_4(t)$ is proportional to the size of the dynamically correlated region (see Appendix B for details), which is the central observable characterizing dynamical heterogeneity approaching the glass transition [13], as it has been estimated in real experiments [7, 40, 72, 73] as well as molecular simulations [71, 74, 75].

Figure 2(b) shows the time and temperature evolution of $\chi_4(t)$. It takes a peak near the relaxation timescale $\tau_{\alpha}$, and the peak grows with decreasing temperature, which is the hallmark of dynamical heterogeneity in glassy dynamics. We have then plotted the peak value of $\chi_4(t)$, denoted $\chi_4^*$, versus temperature $T$ for several system sizes in Fig. 3(a). The system size dependence, akin to molec-

ular simulations [74, 76], allowed us to perform finite size scaling. We find that for increasing the system size $L$, $\chi_4^*$ follows a scaling form, $\chi_4^* \sim T^{-\gamma}$ with a critical point at $T = 0$ and the associated exponent $\gamma$. We obtain a scaling collapse for $\chi_4^*(L, T)$ in Fig 3(b) (see also molecular simulation studies [74, 76]), indicating that the dynamics is governed by a diverging correlation length toward $T = 0$, i.e., $\xi \sim T^{-\nu}$ with an exponent $\nu$. The corresponding data for the scalar model are presented in Appendix B. Besides, we obtained a consistent result in terms of $\nu$ by directly measuring the dynamical correlation lengthscale based on the four-point structure factor [77]. The observed critical exponents and predictions (see below) are summarized in Table I. The main goal of this paper is to provide a scaling theory for these critical behavior, connecting the dynamics at finite temperature to a zero-temperature critical point.
<table><caption>TABLE I. Critical exponents ($\gamma$, $\nu$) obtained from finite T simulations in the scalar and tensorial elastoplastic models in two-dimensions, compared with their predicted values proposed in Section V.</caption><thead><tr><th></th><th>Scalar model</th><th>Tensorial model</th><th>Pred. d = 2</th><th>Pred. MF</th></tr></thead><tbody><tr><th>γ</th><td>1.8 ± 0.1</td><td>1.7 ± 0.1</td><td>1.60 ± 0.05</td><td>1.5</td></tr><tr><th>ν</th><td>0.85 ± 0.05</td><td>0.80 ± 0.05</td><td>0.80 ± 0.03</td><td>0.75</td></tr></tbody></table>

## IV. CRITICAL POINT AND EXTREMAL DYNAMICS AT $T = 0^+$

We now consider dynamics at vanishing temperature, $T = 0^+$, and show that it is related to a critical point. At $T = 0^+$ the site with the smallest energy barrier $E_{min} = x_{min}^{3/2}$, the weakest site, always yields first [78], where $x_{min}$ is the corresponding stress required to destabilise the site. Therefore, one can simulate dynamics at $T = 0^+$ by relaxing always the weakest site, instead of relaxing a random site weighted by the relaxation rate $e^{-E(x)/T}$. This algorithm allows us to access dynamical information even at vanishing temperature (see details in Appendix C). It is an example of extremal dynamics, which is well-studied in the context of self-organized-criticality [79] and some disordered materials under quasi-static driving [47, 80], including annealed glasses [81].

In the thermodynamic limit, $L \to \infty$, the extremal dynamics leads to an absorbing condition for $E \le E_c$, where $E_c$ is a critical energy barrier as found (for $T \to 0$) in Ref. [34]. As a result, the distribution of local (activation) energy barriers $P(E)$ vanishes for $E \le E_c$ (see the sketch in Fig. 9(a)). This implies that $\lim_{L \to \infty} P(x) = 0$ for $x \le x_c = E_c^{2/3}$, where $P(x)$ is the distribution of $x$. Thus only a sub-extensive number of sites are found for $x \le x_c$ at a finite $L$. In this paper, we use $P(E)$ and $P(x)$ interchangeably since they carry essentially the same information.

The extremal dynamics consists of a succession of avalanches. In fact, relaxation at a site can change local energy barriers $E$ at different sites by elastic interactions. As long as those are below $E_c$ the corresponding sites belong to the same avalanche [79]. This is a vivid realization of the phenomenon of dynamic facilitation. To characterize the avalanches, we follow the method introduced in the study of extremal dynamics [47, 79]. For a finite $L$, we fix a chosen threshold stability value $x_0$ ($\le x_c$), and define $x_0$-avalanches as the sequences of events for which $x_{min} < x_0$ (see Fig. 16 in Appendix C). Two useful characterizations of avalanche size are the total number of relaxation events $S$ in a given sequence (event-based avalanche size), and the total number of sites $\tilde{S}$ that relaxed at least once during an avalanche (site-based avalanche size). By construction, $\tilde{S} \le S$ and $\tilde{S} \le L^d$. We show snapshots having $S$ and the corresponding $\tilde{S}$ in Fig. 4. The two can be (and, as we shall show, are)

different as a site can relax multiple times within the same avalanche. Thus, the event-based avalanche size $S$ will allows us to quantify the accumulation of relaxation events in the mobile region (as emphasized in a recent molecular simulation study [51]), whereas the site-based avalanche size $\tilde{S}$ will be associated to the spatial extent of dynamically correlated regions (as it contains the essentially same information as the persistence map in Fig. 2(a) and hence $\chi_4$). Sizes of avalanches, $S$ and $\tilde{S}$, depend on the threshold $x_0$. They grow with increasing $x_0$ and diverge at $x_c$.

By systematically exploring different values of $x_0$ one can determine the critical point $x_c$. In fact, one expects that the avalanche distribution $P(S)$ during the extremal dynamics follows a power-law with a scaling form [47, 79, 82]:

$$P(S) \sim S^{-\tau} g(S/S_c), \quad (2)$$

where $g(z)$ is a scaling function and $S_c$ is a cutoff size which takes the form:

$$S_c \sim (x_c - x_0)^{-1/\sigma} f \left( \frac{L^{d_f}}{(x_c - x_0)^{-1/\sigma}} \right), \quad (3)$$

where $1/\sigma$ and $d_f$ are critical exponents and $f(z) = 1$ for $z \gg 1$ and $f(z) = z$ for $z \to 0$. Thus, $S_c \sim L^{d_f}$ when $x_0 \to x_c$, whereas $S_c \sim (x_c - x_0)^{-1/\sigma}$ when $L \to \infty$. The same expressions are expected to hold for the site-based avalanche size $\tilde{S}$, defining exponents $\tilde{\tau}$, $1/\tilde{\sigma}$, and $\tilde{d}_f$. To estimate $S_c$ we use the fact that for $1 < \tau < 3$ the ratio $\langle S^3 \rangle / \langle S^2 \rangle$ is proportional to the cutoff value $S_c$, where $\langle \cdots \rangle = \int_0^\infty dS P(S) (\cdots)$. As we are interested in scaling of $S_c$ with system size, the numerical constant is irrelevant and we define $S_c \equiv \langle S^3 \rangle / \langle S^2 \rangle$. The same expression is used mutatis mutandis for the site-based avalanche size, $\tilde{S}$. More details about avalanche statistics can be found in Appendix C.

We determine $x_c$ and the exponents $1/\sigma$, $1/\tilde{\sigma}$, $d_f$ and $\tilde{d}_f$ by measuring $S_c$ and $\tilde{S}_c$ for different values of threshold $x_0$ and system size $L$ and collapsing them using the scaling form in Eq. (3) for both $S_c$ and $\tilde{S}_c$, as shown in Figs. 5 (a, b). We obtain $x_c = 0.560 \pm 0.001$. The values of the critical exponents are presented in Table II. Figures 5 (c, d) display $S_c$ and $\tilde{S}_c$ as a function of $x_c - x_0$, for various system size $L$. They indeed show that $S_c$ and $\tilde{S}_c$ grow with $x_0$ and saturate due to finite size effects.

Once $x_c$ is determined, we can study the statistics of the system-spanning avalanches relevant to the thermodynamic limit. Thus, we fix $x_0 = x_c$ and measure the distribution $P(S)$ and $P(\tilde{S})$ of the avalanche sizes, see Figs 6 (a,b). We find that avalanche sizes are power-law distributed with eventual cut-off, consistent with the general assumption in Eq. (2). The scaling form in Eq. (2) collapses the data using $S_c \sim L^{d_f}$ and the previously obtained values of $d_f$ and $\tilde{d}_f$, see Figs 6 (c,d). From the data collapse we determine values of avalanche exponents $\tau$ and $\tilde{\tau}$ (see Table II). The successful collapse of the data confirms the validity of the scaling ansatz as well as the value of critical exponents.
<figure><img src="image_4.png" /><figcaption>FIG. 4. Snapshots of an avalanche formation for the extremal dynamics of the tensorial model with $L = 256$. Event-based avalanche size is $S \simeq 2.7 \times 10^5$, while the site-based avalanche size is $\tilde{S} \simeq 3 \times 10^4$. Purple shows immobile sites (zero event), and the colorbar shows the number of relaxation events in mobile sites.</figcaption></figure>

The above analysis reveals that the extremal dynamics of our model of glass forming-liquid displays, scale-free, avalanche-type dynamics akin to the ones of other disordered systems under external loading [45, 46].

Coming back to the distribution of energy barriers $P(E)$ and the corresponding distribution $P(x)$, one would expect that it continuously vanishes at $x_c$ as $P(x) \sim (x - x_c)^\theta$, defining an exponent $\theta$ [41]. This behavior also occurs near the yielding transition of amorphous solids under shearing (with $x_c = 0$) [83], where it affects the scaling of flow properties [38] and plasticity [84], and more generally in glassy systems with long-range interactions [85]. The underlying reason is that at each step of the extremal dynamics, each site receives a stress kick, such that $x_i$ of a site $i$ follows a random process, with an effective absorbing condition at $x_c$. If this random process was a Brownian motion, then one would obtain $\theta = 1$. Yet the kicks are much more broadly distributed, and in higher dimensions, this random process is akin to a Levy flight [86], which can be shown to imply $\theta_{MF} = 1/2$ [87]. In Fig. 7(a), we plot $P(x)$ for the extremal dynamics, measured from configurations just before (or after) each avalanche defined as $x_{min} > x_0 = x_c$, together with $P(x)$ obtained from finite temperature simulations studied in Sec. III. At higher $T$, $P(x)$ shows a broader distribution with a smooth decay. As $T$ is decreased, $P(x)$ converges to the one obtained by the extremal dynamics at $T = 0^+$, whose shape is consistent with $P(x) \sim (x - x_c)^\theta$.

The following two features of $P(x)$ connect the extremal and finite temperature dynamics. First, the point

<figure><img src="image_5.png" /><figcaption>FIG. 5. Statistical properties of avalanches during the extremal dynamics at $T = 0^+$ for the tensorial model. (a, b): Scaling collapse for the cutoff size $S_c = \langle S^3 \rangle / \langle S^2 \rangle$ based on Eq. (3), for various $L$ for the event-based (a) and site-based (b) avalanche sizes, which determines the critical threshold $x_c$ and critical exponents, $1/\sigma$, $1/\tilde{\sigma}$, $d_f$, and $\tilde{d}_f$. (c, d): $S_c$ versus $x_c - x_0$. The red dashed line corresponds to $S_c \sim (x_c - x_0)^{-1/\sigma}$ in (c) and $\tilde{S}_c \sim (x_c - x_0)^{-1/\tilde{\sigma}}$ in (d).</figcaption></figure>

$x_c$ or the corresponding energy scale $E_c = x_c^{3/2}$ controls the effective energy barrier associated to $\tau_\alpha$. Indeed, at small temperature the dynamics proceeds by relaxing the sites with smallest barriers, i.e., sites having barriers close to $E_c$ [34, 41]. This is in agreement with the relaxation time $\tau_\alpha$ observed in Fig. 1(b), which scales as $\tau_\alpha \sim e^{E_c/T}$.

Second, more importantly for what follows, in a finite system of size $N = L^d$, the typical scale of $x_{min}$ and the second lowest $x$, denoted as $x_{second}$, is expected to follow a power law, $\langle x_{second} - x_{min} \rangle \sim N^{-\delta}$. As we shall show, this plays a key role in the characterization of thermal avalanches. Since the energy difference between the lowest and second lowest activation energies is given by $E_{second} - E_{min} \sim x_{second} - x_{min}$, we conclude

$$ \langle E_{\text{second}} - E_{\text{min}} \rangle \sim N^{-\delta}. \quad (4) $$

Thus the exponent $\delta$ characterizes an important feature of the energy barrier relevant for low temperature dynamics. We numerically confirmed this scaling in Fig. 7(b), and the obtained value of $\delta$ is reported in Table II. Extreme value statistics argument would suggest $\delta = 1/(1 + \theta)$ [38] (although near the yielding point deviations from this relation have been reported in finite dimension [88], and explained in [89]). Within the mean-field theory [87], one finds $\delta_{MF} = 1/(1 + \theta_{MF}) = 2/3$, a value close to what we observe in Fig. 7(b).

The results presented in this section show that the extremal dynamics of the tensorial elastoplastic model of super-cooled liquids is governed by system-spanning avalanches. We have fully characterized the associated
<figure><img src="image_6.png" /><figcaption>FIG. 6. (a, b): Distribution of avalanche size $P(S)$ (a) and $P(\tilde{S})$ (b) for a stability threshold $x_0 = x_c$ for the tensorial model, with varying the system size $L$. The dashed lines follow $P(S) \sim S^{-\tau}$ (a) and $P(\tilde{S}) \sim \tilde{S}^{-\tilde{\tau}}$ (b). (c, d): The corresponding scaling collapse following Eq. (2) with $S_c \sim L^{d_f}$ (c). The same expressions are used for $P(\tilde{S})$ in (d).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 7. (a): $P(x)$ from finite $T$ simulations and extremal dynamics at $T = 0^+$ for the tensorial model with $L = 256$. The red arrow indicates $x_c$. (b): Average of $E_{\text{second}} - E_{\text{min}}$ obtained from the extremal dynamics for various $N = L^d$. The red dashed line corresponds to $\langle E_{\text{second}} - E_{\text{min}} \rangle \sim N^{-\delta}$.</figcaption></figure>

zero temperature critical point by obtaining all relevant exponents, summarized in Table II. We also obtained qualitatively and quantitatively (e.g., critical exponents) similar results in the scalar model, which are presented in Appendix C. Remarkably, we find that the scalar and tensorial models display very close values of the critical exponents and likely correspond to the same universality class.

## V. SCALING THEORETICAL ARGUMENTS

In the following, we develop a scaling theory that connects dynamical heterogeneities observed in finite temperature simulations (Sec. III) and the zero-temperature critical point, and associated avalanches, of the extremal

<table><caption>TABLE II. Critical exponents obtained from extremal dynamics simulations at $T = 0^+$ for the scalar and tensorial elasto-plastic models in two-dimensions. The reported error for the measured exponents corresponds to the range of parameters over which the power law behaviour successfully collapses the data.</caption><thead><tr><th>Scalar model</th><th>Tensorial model</th></tr></thead><tbody><tr><td>$\delta = 0.64 \pm 0.01$</td><td>$\delta = 0.62 \pm 0.01$</td></tr><tr><td>$\tau = 1.25 \pm 0.05$</td><td>$\tau = 1.30 \pm 0.05$</td></tr><tr><td>$d_f = 2.3 \pm 0.1$</td><td>$d_f = 2.3 \pm 0.1$</td></tr><tr><td>$1/\sigma = 2.2 \pm 0.1$</td><td>$1/\sigma = 1.95 \pm 0.03$</td></tr><tr><td>$\tilde{\tau} = 1.25 \pm 0.05$</td><td>$\tilde{\tau} = 1.35 \pm 0.03$</td></tr><tr><td>$\tilde{d}_f = 2.00 \pm 0.02$</td><td>$\tilde{d}_f = 1.95 \pm 0.05$</td></tr><tr><td>$1/\tilde{\sigma} = 1.9 \pm 0.1$</td><td>$1/\tilde{\sigma} = 1.75 \pm 0.03$</td></tr></tbody></table>

dynamics (Sec. IV). We will also discuss other important physical consequences, such as the time evolution of avalanche sizes and the Stokes-Einstein violation.

## A. Length scale of dynamical heterogeneity

We consider the effect of finite temperature $T$ on the extremal dynamics. As we will discuss below, the breakdown of the condition for the extremal dynamics naturally leads to the lengthscale $\xi$ of dynamical heterogeneity at finite $T$. The proposed picture is schematically depicted in Fig. 9.

At finite but low $T$, the dynamics is expected to be extremal on large but finite lengthscales. To understand the underlying mechanism, let us consider a finite size system. For a fixed system size, if $T$ is small enough, the site with the lowest activation energy $E_{\text{min}}$ typically relaxes first like for $T = 0^+$ (the probability to relax the site with the second lowest energy is negligible). This relaxation corresponds to a slow "nucleation" event when $E_{\text{min}} \approx E_c$ that occurs on a timescale $\tau_{\alpha} \sim e^{E_c/T}$, as schematically shown in Fig. 9(a). The nucleation event lowers energy barriers of some other sites of the system, which we call "facilitated" sites, due to elastic interactions, and causing a sequence of faster events when $E_{\text{min}} < E_c$, as shown in Fig. 9(b). This cascade process forms an avalanche, which eventually stops. A new avalanche starts again once another nucleation event occurs. The dynamics is thus intermittent with the power-law avalanches discussed in Sec. IV. Clearly, there is a typical size $N_T$ above which this picture breaks down. Indeed, the above description holds when the thermal energy $T$ is much smaller than the typical energy difference between the lowest and second lowest activation energies, $\langle E_{\text{second}} - E_{\text{min}} \rangle$, otherwise sites with higher energy barriers (such as $E_{\text{second}}$) might relax and dynamics cease to be extremal. According to Eq. (4), $\langle E_{\text{second}} - E_{\text{min}} \rangle$ depends on the system size $N$. Thus, the condition for the extremal dynamics which interplays $T$ and $N$ is given by $T \ll \langle E_{\text{second}} - E_{\text{min}} \rangle \sim N^{-\delta}$, i.e., $N \ll N_T \sim T^{-1/\delta}$.

In consequence, when the system size $N$ is too large
at fixed T, the above description cannot hold. In this case, in particular in the thermodynamic limit, multiple nucleation events followed by avalanches will take place (independently) in parallel in the system, as sketched in Fig. 9(c). The cut-off length scale ξ encompassing a single avalanche is not limited by system size but by other avalanches in the system, which can be estimated assuming that the locations of these avalanches are in- dependent. Such an approximation will be accurate if structural spatial correlations are not prepunderant in this system, as discussed below. Assuming finite size scaling, ξ^d must then corresponds to the largest system size N_T for which extremal dynamics at finite T holds, i.e., ξ^d ~ N_T. This provides the link between finite size zero-temperature avalanches and thermal ones. More- over, it also directly predicts that the size of dynam- ically correlated regions characterized by χ_4^*, which is given by the cross-over size above which the condition for the extremal dynamics breaks down: χ_4^* ~ ξ^d_f lead- ing to χ_4^* ~ T^-(d_f/d^*). To conclude, we derive two scaling relations for thermal avalanches:

$$ \gamma = \frac{\tilde{d}_f}{d\delta}, \qquad (5) $$

$$ \nu = \frac{1}{d\delta}. \qquad (6) $$

These scaling relations connect dynamical hetero- geneities at finite T (characterized by γ and ν) and the distribution of local energy barriers at T = 0^+ (charac- terized by δ) together with the morphology of avalanches during the extremal dynamics (by d̃_f). We thus predict γ and ν by using δ and d̃_f measured in the extremal dy- namics in Sec. IV. These predictions are in good agree- ment with finite temperature simulations in Sec. III, as summarized in Table I. We also find that the mean-field predictions using δ_MF = 2/3 and d̃_f = d = 2 lead to similar values.

This argument is expected to break down if large spa- tial correlations characterize the structure of the sys- tem, causing the locations where avalanches start to be correlated. Indeed, as is always the case in disor- dered materials at zero temperature, any critical thresh- old such as x_c or E_c must display finite-size fluctuations Δx_c ~ ΔE_c. These fluctuations are described by some exponent ν', such that in a system of finite size L, one has ΔE_c ~ L^-(1/ν'). The value of ν' is affected by spa- tial correlations. In general one must have ν' ≥ ν, since the fluctuations of Δx_c or ΔE_c must be at least as large as the typical distance between most unstable sites in a quiescent system E_second - E_min ~ L^-(1/ν): indeed, E_c cannot be more precisely defined than this difference. We expect our argument above to hold when ν = ν', corre- sponding to ΔE_c ~ E_second - E_min. ν' can be related to previously introduced exponents, 1/δ and d̃_f, as follows. If a finite system displays a system spanning avalanche ˜S ~ ˜S_c ~ L^d_f and entirely rearranges, then the change of x_c, or equivalently that of E_c, must be of order Δx_c.

According to Eq. (3), an avalanche of size ˜S ~ ˜S_c ~ L^d_f is associated with a characteristic change of the critical threshold Δx_c = x_c - x_0 ~ L^-(δ̃d_f) (when z ≈ 1 in f(z)), implying that ξ ~ Δx_c^-(δ̃d_f). Since Δx_c ~ ΔE_c and ΔE_c ~ ξ^-(1/ν'), we obtain ν' = 1/(δ̃d_f). Using numerical values for δ̃ and d̃_f, this expression leads to ν' ≈ 0.95 and ν' ≈ 0.9, respectively, for the scalar and tensorial models. We thus have in the present system ν ≈ ν', sup- porting our assumption of independent avalanches. Note however that such an equality does not need to hold in general, especially in large d.

## B. Time evolution of thermal avalanches

We now focus on the time evolution of the size of thermal avalanches, based on the scaling results for the extremal dynamics. In Sec. IV, we introduced x_0- avalanches to more generally probe the critical behavior at x_0 = x_c. This turns out to be useful also to work out the relationship between time and length scales for thermal avalanches. In fact, x_0 is associated to a typical energy scale E(x_0) = x_0^3/2. Hence, it carries information both about the typical time-scale τ(x_0) ~ e^E(x_0)/T and the avalanche (cut-off) size ˜S_c(x_0).

According to the scaling argument discussed before, even at a finite temperature T, the extremal dynamics can be applied to a finite system, in particular, smaller than N_T. The cutoff size of x_0-avalanches, ˜S_c, follows ˜S_c ~ (x_c - x_0)^-(1/δ̃), see Eq. (3) and Fig. 5(b). During an extremal dynamics, the duration of an avalanche will be dominated by the relaxation of the most stable site it involves, of characteristic timescale τ(x_0) ~ e^E(x_0)/T. We thus obtain a relation connecting the duration of avalanches, τ(x_0) (measured in the unit of the relax- ation time τ_α ~ e^E_c/T), with the cutoff size ˜S_c as ˜S_c^-(δ̃) ~ (x_c - x_0) ~ E_c - E(x_0) ~ T ln(τ_α/τ(x_0)). There- fore, avalanches in such a finite system and over a finite time interval t (smaller than the relaxation time scale τ_α) are intermittent. Namely, the system is quiescent most of the time, yet when it is not, its avalanche size ˜S(t) is power-law distributed with a cutoff size ˜S_c(t), which is given by

$$ \tilde{S}_c(t) \sim [T \ln(\tau_\alpha/t)]^{-1/\delta}. \qquad (7) $$

Thus, we find that the size of avalanches grows very slowly – only logarithmically – with time.

The above prediction can be tested experimentally or in molecular dynamics simulations [42, 51]. In Ref. [51], the time evolution of a chord length ⟨ℓ⟩ characterizing the linear size of mobile domains has been measured. Since we expect ˜S_c ~ ⟨ℓ⟩^d_f, our prediction for ⟨ℓ⟩ is given by

$$ \langle l \rangle^{\text{Pred}} = A(T) [T \ln(\tau_{\alpha}/t)]^{-1/(\delta \tilde{d}_f)}, \qquad (8) $$

where A(T) is a function with a finite limiting value as T → 0. In Fig. 8, we compare the molecular simulation
<figure><img src="image_8.png" /><figcaption>FIG. 8. Comparison between the molecular simulation data (empty points) in Ref. [51] and our theoretical prediction (solid curves) in Eq. (8) for the time evolution of a linear avalanche size $\langle l \rangle$. We set $1/\bar{\sigma} = 1.8$ and $\tilde{d}_f = 2$, based on Table II. We use $A = 0.6, 0.7, 0.9, 1.1, 1.2$ and $1.3$ for $T = 0.15, 0.125, 0.11, 0.095, 0.0853$ and $0.0775$, respectively.</figcaption></figure>

data for a two-dimensional polydisperse mixture [51] and our prediction. The prediction is very good, in particular, at lower temperatures where the elastoplastic description is supposed to work well.

In Appendix B, we also connect $\tilde{S}_c(t)$ and $\chi_4(t)$ explicitly. This allows us to predict the time evolution of $\chi_4(t)$ for times $1 \ll t \ll \tau_\alpha$ – a prediction that agrees with the numerics.

### C. Decoupling between diffusion and relaxation

We now show that the scaling theory developed above directly implies decoupling between diffusion and relaxation as found in super-cooled liquids, the so-called Stokes-Einstein violation [52–56].

One of the remarkable aspects found in the numerical simulations in Sec. IV is that sites relax *multiple times* within the same avalanche. In terms of a dynamical trajectory, a site waits a long time (remains immobile) before relaxing, but once it relaxes, it redoes multiple times within a short period of time. This is characterized by a difference between the so-called persistence time and exchange time [57–59] (or the caging time in molecular simulations [61, 90]), which has been argued to be the key ingredient of decoupling between diffusion and relaxation in super-cooled liquids [57–61].

In our case, this effect originates from the difference between the event-based and site-based avalanche sizes characterized by $d_f > \tilde{d}_f$. To connect it to the zero temperature critical point, let us focus on the formation of a thermal avalanche whose time-scale is the order of $\tau_\alpha$ and liner size is $\xi$. During the formation, a single

site relaxes order of $\xi^{d_f}/\xi^{\tilde{d}_f}$ times. Assuming that each relaxation gives a random displacement to particles in the neighborhood of this site, their diffusion constant $D$ must be proportional to the rate for the relaxation event, $\xi^{(d_f-\tilde{d}_f)/\tau_\alpha}$, leading to a Stoke-Einstein breakdown $D\tau_\alpha$ of order $D\tau_\alpha \sim \xi^{(d_f-\tilde{d}_f)} \sim T^{-h}$ that diverges at vanishing temperature with $h = \nu(d_f - \tilde{d}_f)$. Conventionally, the Stokes-Einstein violation is considered a consequence of spatially heterogeneous dynamics [53]. We directly connect the former and the length scale of dynamical heterogeneity $\xi$. On top of that, our scaling theory emphasizes accumulations of multiple events inside a mobile region as the microscopic mechanism leading to the Stokes-Einstein violation.

Note that our scaling argument does not predict a fractional Stokes-Einstein violation in which $D \sim \tau_\alpha^{-\zeta}$, or $D\tau_\alpha \sim \tau_\alpha^{1-\zeta}$, but instead $D\tau_\alpha \sim (\log \tau_\alpha)^{1-\zeta'}$ where $1-\zeta' = \nu(d_f - \tilde{d}_f)$. The former is the fit that is usually conjectured from experimental data [91, 92] with $\zeta \approx 0.8$. However, given the small value of $1-\zeta$ the latter fit is also a viable option [93].

We have numerically tested our prediction for the Stokes-Einstein violation in Fig. 10, showing the product $D\tau_\alpha$ in the tensorial model. In this plot, $D$ is estimated numerically using tracer particles that jump randomly by one lattice spacing each time relaxation occurs in their current site, similarly to what was originally done for kinetic constrained models [57, 58] (see Appendix B for details). $D\tau_\alpha$ increases with decreasing $T$, following the scaling prediction $D\tau_\alpha \sim T^{-h}$ with $h = \nu(d_f - \tilde{d}_f)$ measured independently in Secs. III and IV. The observed amount of the violation we find is not large and hence more representative of strong glass-forming liquids than fragile ones since it has been reported that the magnitude of the violation and fragility are correlated [53, 94]. We will come back to this point in Sec. VI.

## VI. CONCLUSION AND DISCUSSION

We provided a theoretical description of dynamical heterogeneities in super-cooled liquids based on the assumption that local rearrangements are elastically coupled. The elasto-plastic models we studied offer a quantitative solution for how dynamical correlations can emerge even in cases in which local barriers control the dynamics [18, 21, 32, 95, 96]. Our main result is the theoretical explanation of dynamical correlations in terms of a zero-temperature critical point with the associated scaling relations. This leads to quantitative predictions on the power-law statistics of thermal avalanches testable in more realistic systems. Our study suggests that dynamical heterogeneities in super-cooled liquids should be investigated in terms of the temperature $T$ to seek power-law relations, rather than in terms of the relaxation time scale.

One important aspect of the models we studied is that
<figure><img src="image_9.png" /><figcaption>FIG. 9. Proposed picture for dynamical heterogeneities in glass-forming liquids. (a): At low temperatures, the distribution of activation energy barriers $P(E)$ presents a gap below some energy $E_c$. On a time scale $\tau_\alpha \sim e^{E_c/T}$, a site with a barrier near $E_c$ relaxes (red arrow), which we call a "nucleation event". As a result, due to elastic interactions, other sites may display lower barriers $E < E_c$ (blue arrows), which we call "facilitated sites". They relax on a time scale much faster than $\tau_\alpha$, leading to a rapid sequence of events, forming a thermal avalanche. The corresponding real space picture is shown in (b, c). (b): A nucleation event (red circle) triggers facilitated sites (blue circles) by elastic interactions (wavy arrows). These induced events again induce other sites, forming an avalanche growth. (c): Avalanche growth is cut off due to other avalanches originating from different nucleation events taking place simultaneously in the system. The cutoff length $\xi$ corresponds to the maximum size for which extremal dynamics applies, which defines the length of dynamical heterogeneity at finite temperature.</figcaption></figure>

they encode in a very simple and natural way the coupling of local relaxation and elastic interaction. Their simplicity, combined with the richness of the dynamical behavior – in particular, the emergence of facilitation and dynamical correlations – is a remarkable aspect, as it shows which salient facts one can obtain with minimal physical ingredients. Kinetically constrained models have instead abstract kinetic rules and show a large variety of behaviors depending on the kinetic constraints [12]. Although local excitations identified in the dynamic facilitation scenario [43, 97] would have connections with local activations in our framework, the physical interpretation of kinetic rules in kinetic constrained models is still an open and crucial challenge. Along this line of thought, it would be interesting to devise a kinetic constraint rule effectively incorporating elastoplasticity. Nevertheless,

<figure><img src="image_10.png" /><figcaption>FIG. 10. The Stoke-Einstein decoupling $D\tau_\alpha$, where $D$ is the diffusion coefficient of tracers, for the tensorial elasto-plastic model with $L = 128$. (a): $D\tau_\alpha$ versus $\tau_\alpha$ could suggest an effective power-law behaviour (the fractional Stokes-Einstein violation), with two apparent exponents in the range $\tau_\alpha < 10^4$ and $\tau_\alpha > 10^4$. (b): Our prediction suggests instead a power-law behaviour in terms of $T$. The red dashed line represents $D\tau_\alpha \sim T^{-\nu(d_f - \tilde{d}_f)}$.</figcaption></figure>

various concepts and theoretical tools developed in the study of kinetically constrained models provide important guidelines and, in fact, played an essential role for our analysis of thermal elasto-plastic models. It has been demonstrated in the dynamic facilitation theory that a critical point in an extended non-equilibrium phase diagram influences glassy dynamics [98, 99]. It would also be interesting to study whether such a critical point exists in elastoplastic models.

Quantitatively, the magnitude of dynamical heterogeneities we simulated is comparable to most supercooled liquids (as the estimated $\chi_4^*$ increases by about two decades as the glass transition is reached [7, 73]), whereas the magnitude of the Stoke-Einstein breakdown is comparable to that of rather strong glass-forming liquids [53]. According to our scaling argument, such a breakdown will increase if the system shows more intensive accumulations of multiple relaxations characterized by larger $\nu(d_f - \tilde{d}_f)$. This effect could be achieved by imposing that some of the model parameters (such as the local values of the yield stress, or how sites are coupled to the elastic field) are randomly distributed [100], instead of being single-valued as assumed here for simplicity. These effects are expected in glass-forming liquids due to the presence of structural heterogeneity of local orders [101–104]. Such generalization would allow us to study the structure-dynamics relationship [105–107] in elastoplastic models. Further improvements may be achieved by adding fluctuations and non-linearities to the propagator, which are present at short-range [108, 109]. An interesting line of research to develop quantitative models is a mapping from a molecular simulation to an elastoplastic model [110–112] or the one pursued in Refs. [113, 114], which uses machine learning methods and the so-called softness field to obtain quantitative effective models.

Our results also underline important themes to study in the future, including the nature of local rearrangements in glass-forming liquids, and their connection to
fragility [2]. Concerning the latter, the current elastoplastic models correspond to strong glass-formers with Arrhenius behavior, with activation energy given by the magnitude of the gap $E_c$ entering the distribution of local barriers [34, 41]. This point results from the simplifying assumption that the energy scale of local rearrangements does not vary with temperature. In an improved (still simplified) model where all elastic energies follow the high-frequency elastic modulus $G_\infty(T)$ of the material, the activation energy $E_c$ will be proportional to $G_\infty(T)$, a correlation known to exist in some glass-forming liquids [18, 115]. More recent works relate this energy scale of local rearrangements to local (instead of global) elasticity [20], plasticity [32], or alternatively to the varying geometry of elementary rearrangements under cooling [116]. Thus, it is worthwhile to reveal the relationship between the activation energy barriers in liquids and low-energy excitations in glasses [117–119]. Local energy barriers would also be related to locally-favored structures in microscopic perspectives [102, 120]. More progress along those lines will be instrumental to understand what controls fragility in glass-forming liquids.

Note that although we mainly focused our attention on the theoretical scenario based on local barriers driving the dynamics [20, 21, 32] together with facilitation and avalanches, the physical phenomena discussed in this work are more general. For example, they can also apply to cases in which cooperative rearrangements take place. In these perspectives, in particular within Random First Order Transition theory, the local relaxation event would correspond to a cooperative rearrangement [51, 121]. Obviously, our current elastoplastic models do not take into account growing cooperativity as a static correlation, and they have a singularity only at $T = 0$ in contrast to RFOT with a finite temperature singularity at $T_K > 0$. One could incorporate a growing static correlation in the models by increasing the number of sites involving a thermal activation process. It would be interesting to work out precise predictions in this case.

Finally, we expect our scaling theory which connects spatial correlations at finite temperature to extremal dynamics at zero temperature, to be relevant for a broader class of problems beyond the context of the glass transition studied here. Phenomena in which the implications of these arguments could be studied include the creep flow of disordered materials [122–128] or that of pinned elastic interfaces [47, 129, 130] below the threshold force where they spontaneously flow.

## ACKNOWLEDGEMENT

We thank the Simons collaboration for discussions, in particular J. Baron, L. Berthier, C. Brito, C. Gavazzoni, E. Lerner, C. Liu, D. R. Reichman, and C. Scalliet. We thank S. Patinet for insightful discussions. We appreciate C. Scalliet for sharing the data in Ref. [51]. GB thanks J.P. Bouchaud for discussions, in particular on the

Bak-Sneppen model and models of dynamic facilitation. MW thanks A. Rosso and M. Muller for discussions on avalanche-type response over the years, and T. de Geus, W. Ji and M. Pica Ciamarra for discussions. MW acknowledges support from the Simons Foundation Grant (No. 454953 Matthieu Wyart) and from the SNSF under Grant No. 200021-165509. GB acknowledges support from the Simons Foundation Grant (No. 454935 Giulio Biroli).

# Appendix A: Thermal elastoplastic models

## 1. Scalar model

We study a scalar elastoplastic model [34] in a two-dimensional lattice whose linear box length is $L$ using the lattice constant as the unit of length. For each site, we assign local shear stress $\sigma_i$ (scalar variable) at a position $\mathbf{r}_i$.

The dynamical rule for the simulation model is akin to Monte-Carlo dynamics [131]. We pick a site, say $i$, up randomly among $L^2$ sites. If $\sigma_i$ is greater (or lower) than or equal to a threshold $\sigma_Y > 0$ (or $-\sigma_Y < 0$), namely, $|\sigma_i| \ge \sigma_Y$, this site shows a plastic event: $\sigma_i \to \sigma_i - \delta\sigma_i$, where $\delta\sigma_i$ is the local stress drop. We use an uniform threshold, $\sigma_Y = 1$. Instead, if $|\sigma_i| < \sigma_Y$, with probability $e^{-E(\sigma_i)/T}$, where $E(\sigma_i)$ is a local energy barrier and $T$ is the temperature, this site shows a plastic event: $\sigma_i \to \sigma_i - \delta\sigma_i$. This corresponds to a plastic rearrangement induced by a local thermal activation. We employ $E(\sigma_i) = (\sigma_Y - |\sigma_i|)^\alpha$ with $\alpha = 3/2$ [132]. By introducing the local stress distance to threshold, $x_i = \sigma_Y - |\sigma_i|$, we can rewrite $E(x) = x^{3/2}$. This specific form of the local energy barrier is suggested by molecular simulation studies [32, 67] and previous elastoplastic models under shear [41, 65]. The stress drop $\delta\sigma_i$ associated with a plastic event is a stochastic variable. In this paper, we use $\delta\sigma_i = (z + |\sigma_i| - \sigma_Y)\text{sgn}(\sigma_i)$, where $\text{sgn}(x)$ is the sign function and $z > 0$ is a random number drawn by an exponential distribution, $p(z) = \frac{1}{z_0}e^{-z/z_0}$. $z_0$ is the mean value and we set $z_0 = 1$. This exponential distribution would be realistic according to molecular simulations in Ref. [133].

A local plastic event at site $i$ influences all other sites ($\forall j \ne i$) as

$$ \sigma_j \to \sigma_j + G_{\mathbf{r}_{ji}}^{\psi_i} \delta\sigma_i, \quad (A1) $$

where $\mathbf{r}_{ji} = \mathbf{r}_j - \mathbf{r}_i$ and $\psi_i \in [0, \pi/2)$ is a random orientation of the Eshelby kernel $G_r^\psi$. Numerical implementation of $G_r^\psi$ is described in Ref. [34].

Similar to the Monte-Carlo dynamics, we repeat the above attempt $L^2$ times, which corresponds to unit time.

For the initial condition, we draw the local stress $\sigma_i$ ($\forall i$) randomly while keeping the force balance, i.e., the sum of stresses along each row and column of lattice sites is strictly zero [134, 135]. To study dynamical properties
at the steady-state, we monitor the waiting time dependence of observables, and we report them only at the steady-state, discarding the initial transient part.

## 2. Tensorial model

We have implemented a two-dimensional elastoplastic model in which we account for the tensorial nature of the shear stress tensor. In this tensorial version of the elastoplastic model, the state of each site $i$ is described by its local shear stress tensor $\tilde{\sigma}_i$. Note that symbols $\sigma$ and $\tilde{\sigma}$ are also used as critical exponents in the main text, not to be confused with local shear stress defined here. The shear stress tensor is traceless and symmetric and hence in two dimensions it is defined by two independent components: $\tilde{\sigma}_{xx,i}$ and $\tilde{\sigma}_{xy,i}$.

The local yield stress is defined by a surface in the shear stress space, with the region inside and outside the surface corresponding to mechanically stable (elastic, immobile) and unstable (plastic, mobile) states, respectively. The minimum amount of shear stress required to make a site unstable is the distance to the yield surface, and we denote its magnitude by $x$, as schematically shown in Fig 11. We choose the local yield surface to consist of two parallel lines at an angle $\theta_Y$ with respect to the $\tilde{\sigma}_{xx}$ axis in shear stress space, centered at zero shear stress and separated by $2\sigma_Y$ (see Fig 11). The local yield surface is assigned for each site $i$. We initiate $\theta_{Y,i}$ with a uniformly distributed random number in $[0, 2\pi)$.

When a site $i$ becomes unstable it undergoes a plastic event over a timescale $\tau_0$: $\tilde{\sigma}_{xx,i} \rightarrow \tilde{\sigma}_{xx,i} - \delta\tilde{\sigma}_{xx,i}$ and $\tilde{\sigma}_{xy,i} \rightarrow \tilde{\sigma}_{xy,i} - \delta\tilde{\sigma}_{xy,i}$, where the amount of stress drops $\delta\tilde{\sigma}_{xx,i}$ and $\delta\tilde{\sigma}_{xy,i}$ are given by

$$ \delta\tilde{\sigma}_{xx,i} = -(z-x) \sin(\theta_{Y,i}) \operatorname{sgn}(\tilde{\sigma}_{xx,i}), \quad (A2) $$

$$ \delta\tilde{\sigma}_{xy,i} = -(z-x) \cos(\theta_{Y,i}) \operatorname{sgn}(\tilde{\sigma}_{xy,i}), \quad (A3) $$

respectively. $\operatorname{sgn}(x)$ is the sign function and $z$ is a random number drawn from an exponential distribution $p(z) = e^{-z/z_0}/z_0$ with $z_0 = 1$. The duration of a plastic event $\tau_0$ is accounted for by triggering the relaxation with a probability per unit time $1/\tau_0$ whenever the site is unstable. In an athermal system sites can only relax by first becoming unstable. At finite temperature $T$ stable sites undergo relaxation at the rate $e^{-E(x)/T}/\tau_0$, where $E(x) = x^{3/2}$ is the local energy barrier [41]. After each plastic event, we redraw the angle of the yield surface $\theta_{Y,i}$ from a uniform random distribution.

Note that to simulate such a dynamics at low temperatures we implement a Gillespie type of algorithm [41], which operates as follows. Consider an event occurring at some time $t$. Following it, a relaxation time $\tau_i$ for each site $i$ is chosen with an exponential distribution of mean $e^{-E(x_i)/T}/\tau_0$. The next event corresponds to the smallest $\tau_i$, leading to a plastic event on the corresponding site, which occurs at a time $t + \tau_i$. At that point, stresses are computed once again, the variables $\tau_i$ are sampled from

<figure><img src="image_11.png" /><figcaption>FIG. 11. A schematic plot of the yield surface for the tensorial model. The blue dot shows the state of a site in the shear stress space. The yield surface is described by two parallel lines separated by $2\sigma_Y$ and each line makes an angle of $\theta_Y$ with $\tilde{\sigma}_{xx}$ axis. The distance from the yield surface is shown by $x$.</figcaption></figure>

their new distributions. this algorithm is then repeated iteratively.

To maintain the force balance, a plastic event at site $i$ redistributes the shear stress field in the other sites following the force dipole propagator. In Fourier space the elements of elastic kernel are given by,

$$ \hat{G}_{xx,xx}(\mathbf{q}) = - \frac{(q_x^2 - q_y^2)^2}{(q_x^2 + q_y^2)^2}, \quad (A4) $$

$$ \hat{G}_{xy,xy}(\mathbf{q}) = - \frac{4q_x^2 q_y^2}{(q_x^2 + q_y^2)^2}, \quad (A5) $$

$$ \hat{G}_{xx,xy}(\mathbf{q}) = \hat{G}_{xy,xx}(\mathbf{q}) = - \frac{2q_x q_y (q_x^2 - q_y^2)}{(q_x^2 + q_y^2)^2}, \quad (A6) $$

where $\mathbf{q} = (q_x, q_y)$ is the Fourier vector. For a discrete system with periodic boundary conditions, we introduce a correction term to the Fourier modes, given by $q_x^2 = 2 - 2 \cos(2\pi n_x/L)$, $q_y^2 = 2 - 2 \cos(2\pi n_y/L)$ and $q_x q_y = 2 \sin(2\pi n_x/L) \sin(2\pi n_y/L)$ where $n_\alpha = -L/2 + \{1, \dots, L\}$ with $\alpha = x, y$.

# Appendix B: Characterizations of dynamical heterogeneities for finite temperature simulations

## 1. Four-point correlation function

We explain how to analyze correlation functions for finite temperature simulations. We first consider the persistence two-point time correlation function, $\langle\pi(t)\rangle_{\text{time}}$, which is defined by $\pi(t) = \frac{1}{L^d} \sum_i p_i(t)$, where $p_i(t) = 1$ (immobile) if the site $i$ did not show a plastic event until time $t$ from $t = 0$, and $p_i(t) = 0$ (mobile) otherwise. $\langle\cdots\rangle_{\text{time}}$ denotes the time average at the stationary state. $\langle\pi(t)\rangle_{\text{time}}$ for the scalar and tensorial models are presented in Fig. 1 of Ref. [34] and Fig. 1(a) in the main
text, respectively. We then measure a four-point correlations function, $\chi_4(t)$, defined by

$$ \chi_4(t) = L^d \left( \langle \pi^2(t) \rangle_{\text{time}} - \langle \pi(t) \rangle_{\text{time}}^2 \right). \quad (B1) $$

$\chi_4(t)$ for the scalar and tensorial models are presented in Fig. 3 of Ref. [34] and Fig. 2(b) in the main text, respectively. $\chi_4(t)$ quantifies the size of the dynamically correlated region, because one can rewrite it as

$$ \chi_4(t) = \frac{1}{L^d} \sum_{i,j} \left( \langle p_i(t) p_j(t) \rangle_{\text{time}} - \langle \pi(t) \rangle_{\text{time}}^2 \right) \quad (B2) $$

$$ = \frac{1}{L^d} \sum_i \sum_k \langle \phi_i(t) \phi_{i+k}(t) \rangle_{\text{time}}, \quad (B3) $$

where $\phi_i(t) = p_i(t) - \langle \pi(t) \rangle_{\text{time}}$. For example, $\phi_i(\tau_\alpha) = \pm 1/2$. Therefore, $\chi_4(t)$ is proportional to the average number of sites correlated dynamically. Therefore, its peak value, $\chi_4^*$, contains essentially the same information as $\tilde{S}$, in particular, $\chi_4^* \sim \tilde{S}_c$ at $T = 0^+$.

Figure 12(a) shows $\chi_4^*$ versus $T$ for several system sizes $L$ for the scalar model. One can see a scaling regime, $\chi_4^* \sim T^{-\gamma}$ at lower $T$ and larger $L$. A scaling collapse is obtained in Fig. 12(b), which determines another critical exponent $\nu$ associated with a lengthscale of dynamical heterogeneity. The obtained values for the critical exponents are reported in Table II.

## 2. Dynamical correlation lengthscale

We consider extracting a correlation length directly instead of performing finite-size scaling. To this end, we measure the spatial dependence of the four-point structure factor, $S_4(q, t)$ [77], defined by

$$ S_4(q,t) = \frac{1}{L^d} \sum_{ij} \left( \langle p_i(t) p_j(t) \rangle_{\text{time}} - \langle \pi(t) \rangle_{\text{time}}^2 \right) e^{i \mathbf{q} \cdot (\mathbf{r}_i - \mathbf{r}_j)}, \quad (B4) $$

where $q = |\mathbf{q}|$. In Fig. 13(a), we show $S_4(q, t)$ for the scalar model at $t = \tau^*$ when $\chi_4(t)$ takes the peak value, $\chi_4^* = \chi_4(\tau^*)$. We find that $\tau^*$ is close to $\tau_\alpha$, and thus $S_4(q, \tau^*)$ encodes heterogeneity associated with structural relaxation. We note that at the long wave-length limit, $S_4(q, t)$ converges to $\chi_4(t)$, namely, $\lim_{q \to 0} S_4(q, t) = \chi_4(t)$. We then assume the Ornstein-Zernike form at lower $q$,

$$ S_4(q, \tau^*) = \frac{\chi_4^*}{1 + (q\xi_4)^a}, \quad (B5) $$

where $\xi_4$ is the dynamical correlation length extracted and $a$ is an exponet. As shown in Fig. 13(b), we find this scaling with $a = 2.2$. From this plot, we extract $\xi_4$ and present its temperature dependence in Fig. 13(c). This analysis can be done only for larger systems, $L = 32$, 64, and 128, since the scaling regime cannot be reached in smaller systems within our simulations. We find

<figure><img src="image_12.png" /><figcaption>FIG. 12. The four-point correlation function $\chi_4$ for the scalar model. (a): The peak value $\chi_4^*$ versus $T$ for several $L$ in a Log-Log plot. The dashed line follows $\chi_4^* \sim T^{-\gamma}$. (d): Scaling collapse of $\chi_4^*(L, T)$, which determines $\nu$.</figcaption></figure>

$\xi_4 \sim T^{-\nu}$ with $\nu = 0.9$, which is consistent with the one estimated from the finite size scaling in Fig. 12. Moreover, the $\chi_4^*$ versus $\xi_4$ plot [136, 137] in Fig. 13(d) provides us with the fractal dimensions, $\tilde{d}_f = 2$, which is also consistent with the one measured in the extremal dynamics in Fig. 18.

## 3. Prediction for the four-point correlation function

We connect the cutoff size for the site-based avalanche size, $\tilde{S}_c(t)$, in a given time interval $t$, and the time evolution of the four-point correlation function, $\chi_4(t)$.

We first compute the site-based avalanche size by $\tilde{S}(t) = \sum_{i=1}^N n_i(t)$, where $n_i(t) = 0$ if the site $i$ did not exhibit a plastic event until time $t$, and $n_i(t) = 1$ for mobile sites that relaxed at least once. Thus, $n_i(t)$ can be written by $n_i(t) = 1 - p_i(t)$. The first and second moments of $\tilde{S}(t)$ measured by the time average are given by

$$ \begin{aligned} \langle \tilde{S}(t) \rangle_{\text{time}} &= N(1 - \langle \pi(t) \rangle_{\text{time}}), \\ \langle \tilde{S}^2(t) \rangle_{\text{time}} &= \sum_{i,j} \{ \langle p_i(t) p_j(t) \rangle_{\text{time}} - (2 \langle \pi(t) \rangle_{\text{time}} - 1) \}, \end{aligned} $$
<figure><img src="image_13.png" /><figcaption>FIG. 13. The four-point structure factor $S_4(q, t)$ and the associated correlation length $\xi_4$ for the scalar model. (a): $S_4(q, \tau^*)$ for several temperatures for $L = 128$. (b): The corresponding plot for the Ornstein-Zernike form in Eq. (B5). The dashed-line defines a slope corresponding to the exponent $a$. (c): The extracted $\xi_4$ versus $T$. The dashed line follows $\xi_4 \sim T^{-\nu}$. (b): $\chi_4^*$ versus $\xi_4$. The dashed line follows $\chi_4^* \sim \xi_4^{\bar{d}_f}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 14. (a): Log-Log plot for $\chi_4(t)$ for the scalar model with $L = 128$ for $T = 0.050, 0.040, 0.030, 0.025, 0.020, 0.015$, and $0.013$ (from left to right). (b): Parametric plot to test the prediction in Eq. (B7). The dashed line defines the linear relation.</figcaption></figure>

respectively. Consider a correlation volume of linear extension $\xi$. On this length scale, $\chi_4(t)$ crosses-over toward its value for an infinite system. It is also the largest length for which extremal dynamics applies, implying that $\tilde{S}(t)$ is distributed in a power-law fashion. Thus, $\tilde{S}_c(t)$ can be estimated by $\langle \tilde{S}^2(t) \rangle_{\text{time}}/\langle \tilde{S}(t) \rangle_{\text{time}}$, as given by

$$ \tilde{S}_c(t) = \frac{\sum_{i,j} \{\langle p_i(t)p_j(t) \rangle_{\text{time}} - (2\langle\pi(t)\rangle_{\text{time}} - 1)\}}{N(1 - \langle\pi(t)\rangle_{\text{time}})} \quad (B6) $$

In general, one can expect that $\langle \pi(t) \rangle_{\text{time}}$ follows the (stretched) exponential decay, $\langle \pi(t) \rangle_{\text{time}} \approx e^{-(t/\tau_\alpha)^\beta}$, where $\beta$ is an exponent. Typically, $0 < \beta \le 1$ for equilibrium supercooled liquids. Our elastoplastic models show nearly exponential relaxation with $\beta \approx 1$. We now consider the early time stage, where $\langle \pi(t) \rangle_{\text{time}}$ can be approximated by $\langle \pi(t) \rangle_{\text{time}} \approx 1 - (t/\tau_\alpha)^\beta$. Under such a circumstance, Eq. (B6) and $\chi_4(t)$ defined in Eq. (B2) suggest that $\tilde{S}_c(t) \approx (t/\tau_\alpha)^{-\beta}\chi_4(t)$. Together with Eq. (7) in the main text, we predict the time evolution of $\chi_4(t)$ as

$$ \chi_4(t) \approx A(t/\tau_\alpha)^\beta [T \ln(\tau_\alpha/t)]^{-1/\bar{\sigma}}, \quad (B7) $$

where $A$ is a constant which does not depend on $t$ and $T$.

Figure 14(a) shows a Log-Log plot for $\chi_4(t)$ measured at finite temperature simulations for the scalar model. The initial growth can be fitted effectively by a power-law, $\chi_4(t) \sim t^b$ [138, 139], with $b \approx 1.4$. Instead, our argument in Eq. (B7) predicts a linear growth with a log-

arithmetic correction. In Fig. 14(b), we show a parametric plot to numerically test Eq. (B7) with $\beta = 1$, $A = 0.5$, and $1/\bar{\sigma} = 1.9$ (see Table II). We find that the simulated $\chi_4(t)$ for different temperatures follows our prediction at early times. Deviations from the prediction can be observed on very short timescales. This is presumably due to the fact that on such time scales, the corresponding energy scale is too small compared with $E_c$, which violates the assumption underlying the asymptotic argument of Sec. V.

## 4. Tracer particles

We monitor the diffusion of tracer particles [57, 58] due to the local relaxations. We consider one tracer particle in each site of the elastoplastic model; each tracer particle moves randomly to one of the four nearest neighbors (in $d=2$) after a plastic event in that site. The trajectory of the $k$-th tracer particle is specified by $(x_k(t), y_k(t))$. Typical trajectories of the tracer particles are shown in Fig. 15(a) as a function of time. In a time scale comparable to $\tau_\alpha$, the tracer travels over multiple sites in a very short time and spends most of the time without any activity. We then define the mean-squared displacement $\Delta^2(t)$ of the tracers by

$$ \Delta^2(t) = \frac{1}{N_t} \sum_{k=1}^{N_t} \langle \Delta r_k^2(t) \rangle_{\text{time}}, \quad (B8) $$

where $\Delta r_k^2(t) = (x_k(t) - x_k(0))^2 + (y_k(t) - y_k(0))^2$ and $N_t$ is the number of tracer particles. In Fig. 15(b) we show $\Delta^2(t)$ for different temperatures. We find a diffusive behavior at a larger time, $\Delta^2(t) = Dt$, from which we extract the diffusion coefficient $D$ for each temperature.

# Appendix C: Avalanches at $T = 0^+$

## 1. Extremal dynamics

We explain the extremal dynamics at $T = 0^+$. In the finite temperature simulations described in Sec. A, we
<figure><img src="image_15.png" /><figcaption>FIG. 15. Diffusion of tracer particles for the tensorial model with $L = 128$. (a): The y component of typical trajectory of a tracer particle with $T = 0.015$ ($\tau_{\alpha} = 4.14 \times 10^{12}$). (b): Mean-squared displacement $\Delta^2(t)$ of tracer particles are shown with points. The solid lines follow $\Delta^2(t) = Dt$, from which we extract the diffusion coefficient $D$.</figcaption></figure>

take into account local thermal activation for a plastic event based on the probability $e^{-E(x)/T}$, where $E(x) = x^{3/2}$ at $0 \le x \le 1$. At vanishing temperature, $T = 0^{+}$, this probability is extremely small. Therefore, the site with the smallest $x$, denoted as $x_{\min}$, associated with the lowest energy barrier $E_{\min} = x_{\min}^{3/2}$, shows the next plastic event. Thus, in practice, one can choose the weakest site having $x_{\min}$ sequentially, instead of asking $e^{-E(x)/T}$ each time and waiting until it shows an event. This algorithm enormously accelerates dynamics and allows us to access information about plastic activities even at $T = 0^{+}$. This is the so-called extremal dynamics [47, 79, 80].

Finding one $x_{\min}$ corresponds to one simulation step. This is not directly related to physical time (that is why one can simulate it even at $T = 0^{+}$), yet one can associate the simulation step with the size of an avalanche (see below). Simulations start with the same initial condition used in the finite temperature simulations. The system enters the stationary state after passing the initial transient regime. We carefully checked the stationarity by monitoring the waiting time dependence of $P(x)$. We report data taken only from the stationary state.

In Fig. 16, we show a representative trajectory of $x_{\min}$ during an extremal dynamics simulation at the stationary state, which is an analog of FIG. 5 in Ref. [79] for a model for self-organized criticality. Typically, the weakest site with $x_{\min}(s)$ at a simulation step $s$ induces the next weakest site at step $s+1$ with $x_{\min}(s+1)$ at a neighbor region because of elastic interactions. In particular, $x_{\min}(s) > x_{\min}(s+1)$, when the previous weakest site at step $s$ destabilizes the next weakest site at step $s+1$. Therefore, a sequence of the weakest sites is dynamically correlated, forming an avalanche until the last weakest site is found at an uncorrelated place with a higher value of $x_{\min}$. The determination of uncorrelation and hence the termination of an avalanche has some ambiguity. Thus, following previous works, we introduce the threshold $x_0$ below which a sequence of $x_{\min}$ is correlated. In particular, we define the size of event-based avalanche, $S$, by the number of chosen $x_{\min}$ forming a sequence with $x_{\min} < x_0$ (in other words, the duration of

<figure><img src="image_16.png" /><figcaption>FIG. 16. An example of the evolution of $x_{\min}$ during the extremal dynamics at a steady state for the scalar model. The system size is $L = 256$. A series of event-based avalanche sizes, $S_1, S_2, \dots$, are presented based on the threshold value $x_0 = 0.2$. $x_c = 0.281$ is indicated by the horizontal dashed line.</figcaption></figure>

simulation steps with $x_{\min} < x_0$). As shown in Fig. 16, one can extract a series of event-based avalanche sizes, $S_1, S_2, \dots$, from the trajectory of the extremal dynamics simulation. A given site may be chosen as the weakest site several times during one avalanche formation, which all contribute to $S$. Instead, one can define the size of site-based avalanche, $\tilde{S}$, by the number of sites participating in a single avalanche. By construction, $\tilde{S} \le S$. The distinction between $S$ and $\tilde{S}$ provides us with important physical information about the accumulation of multiple relaxation activities, which leads to an argument about the Stokes-Einstein violation, as discussed in Sec. V.

By construction, $S$ and $\tilde{S}$ depend on the threshold value $x_0$. As $x_0$ is increased, the size of avalanches, $S$ and $\tilde{S}$, increases. As we will discuss further below, avalanches become a system spanning when $x_0 \rightarrow x_c$, where $x_c$ is the critical value associated with the critical energy gap $E_c = x_c^{3/2}$. We will vary $x_0$ systematically to probe the critical behavior associated with $x_c$ (see below).

## 2. Avalanche statistics

We describe how to analyze the avalanche data obtained during $T = 0^{+}$ extremal dynamics simulations. During simulations, we record the series of the event and site-based avalanche sizes, given by $\{S_1, S_2, \dots, S_M\}$ and $\{\tilde{S}_1, \tilde{S}_2, \dots, \tilde{S}_M\}$, respectively, where $M$ is the number of data points. In this paper, we analyze $S$ and $\tilde{S}$ in parallel. Below we will explain how to analyze the data using $S$, but the same procedures are applied for $\tilde{S}$. We first define the $m$-th moments of avalanche distribution
<figure><img src="image_17.png" /><figcaption>FIG. 17. Avalanche distributions $P(S)$ (a) and $P(\tilde{S})$ (b) for several $x_0$ approaching the critical point $x_c = x_0 = 0.281$ for the scalar model with $L = 256$. The dashed lines in (a) and (b) follow $P(S) \sim S^{-\tau}$ and $P(\tilde{S}) \sim S^{-\tilde{\tau}}$, respectively.</figcaption></figure>

$(m = 1, 2, ...)$ by

$$ \langle S^m \rangle = \int_0^\infty dS P(S) S^m = \frac{1}{M} \sum_{k=1}^M S_k^m, \quad (C1) $$

where $P(S)$ is the distribution of avalanches. One expects that $P(S)$ follows a power law distribution,

$$ P(S) \sim S^{-\tau} g(S/S_c), \quad (C2) $$

where $\tau$ is a critical exponent, $S_c$ is a cutoff size, and $g(z)$ is a scaling function. Figure 17 shows $P(S)$ and $P(\tilde{S})$ for several $x_0$ for the scalar model. These plots demonstrate that the size of avalanches (both $S$ and $\tilde{S}$) grow with increasing $x_0$, as expected. In particular, a scale-free, power-law behavior (with the eventual cut-off) is being developed by approaching the critical point $x_c$, which proves that $x_0$ is the relevant parameter that dictates the critical behavior of the system.

Assuming Eq. (C2) and $1 < \tau < 2$, one obtains

$$ \langle S^m \rangle \sim \int_0^\infty dS S^{m-\tau} g(S/S_c) \sim S_c^{m+1-\tau}, \quad (C3) $$

which implies $S_c \sim \langle S^{m+1} \rangle / \langle S^m \rangle$. Thus, in practice, we define $S_c$ by $S_c = \langle S^3 \rangle / \langle S^2 \rangle$. Following Ref. [82], we assume

$$ S_c \sim (x_c - x_0)^{-1/\sigma} f \left( \frac{L^{d_f}}{(x_c - x_0)^{-1/\sigma}} \right), \quad (C4) $$

where $f(z) = 1$ for $z \gg 1$ and $f(z) = z$ for $z \to 0$. Thus, $S_c \sim L^{d_f}$ when $x_0 \to x_c$ and $S_c \sim (x_c - x_0)^{-1/\sigma}$ when $L \to \infty$. Figures 18(a,b) show $S_c$ and $\tilde{S}_c$ approaching $x_0 \to x_c$ for the scalar model. Both $S_c$ and $\tilde{S}_c$ increase with increasing $x_0$ with an eventual saturation due to a finite-size effect. One can see the expected behavior, $S_c \sim (x_c - x_0)^{-1/\sigma}$ at larger $L$ ($\tilde{S}_c$ as well). We then perform the scaling collapse in Figs. 18(c,d) following the scaling form in Eq. (C4). These scaling plots determine the critical point $x_c = 0.281$ and the critical exponents, $1/\sigma$, $1/\tilde{\sigma}$, $d_f$, and $\tilde{d}_f$, for the scalar model. The obtained values are reported in Table II.

Once the critical threshold $x_c$ is determined, we measure avalanche distributions $P(S)$ and $P(\tilde{S})$ at $x_0 = x_c$,

<figure><img src="image_18.png" /><figcaption>FIG. 18. Cut-off in avalanche distributions for the scalar model. (a, b): $S_c = \langle S^3 \rangle / \langle S^2 \rangle$ (a) and $\tilde{S}_c = \langle \tilde{S}^3 \rangle / \langle \tilde{S}^2 \rangle$ (b) versus $x_c - x_0$ for various system sizes. The dashed lines in (a) and (b) follow $S_c \sim (x_c - x_0)^{-1/\sigma}$ and $\tilde{S}_c \sim (x_c - x_0)^{-1/\tilde{\sigma}}$, respectively. (c,d): Scaling collapse assuming Eq. (C4). The critical point $x_c$, exponents $\sigma$, $\tilde{\sigma}$, $d_f$, and $\tilde{d}_f$ are determined by these scaling plots.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 19. (a,b): Distribution of avalanche size $P(S)$ (a) and $P(\tilde{S})$ (b) for a stability threshold $x_0 = x_c$ for the scalar model, with varying the system size $L$. The dashed lines follow $P(S) \sim S^{-\tau}$ (a) and $P(\tilde{S}) \sim \tilde{S}^{-\tilde{\tau}}$ (b). (c): The corresponding scaling collapse following Eq. (C2) with $S_c \sim L^{d_f}$. (d): Same for $P(\tilde{S})$.</figcaption></figure>

as shown in Figs. 19(a,b) for the scalar model. They show characteristic power-law behavior with cut-off $S_c$ and $\tilde{S}_c$ due to a finite size effect, which scales as $S_c \sim L^{d_f}$ and $\tilde{S}_c \sim L^{\tilde{d}_f}$, respectively. The power-law behavior determines $\tau$ and $\tilde{\tau}$ whose values are reported in Table II. We then perform scaling collapses in Figs. 19(c,d), following Eq. (C2), which validates the scaling ansatz and measured critical exponents.

Finally, we compute the distribution $P(x)$ from the
inherent structures of liquids at rest,” The Journal of chemical physics **144**, 124508 (2016).

[25] Hua Tong, Shiladitya Sengupta, and Hajime Tanaka, “Emergent solidity of amorphous materials as a consequence of mechanical self-organisation,” Nature communications **11**, 1–10 (2020).

[26] Bin Wu, Takuya Iwashita, and Takeshi Egami, “Anisotropic stress correlations in two-dimensional liquids,” Physical Review E **91**, 032301 (2015).

[27] Manuel Maier, Annette Zippelius, and Matthias Fuchs, “Emergence of long-ranged stress correlations at the liquid to glass transition,” Physical review letters **119**, 265701 (2017).

[28] David Steffen, Ludwig Schneider, Marcus Müller, and Jörg Rottler, “Molecular simulations and hydrodynamic theory of nonlocal shear-stress correlations in supercooled fluids,” The Journal of Chemical Physics **157**, 064501 (2022).

[29] Elijah Flenner and Grzegorz Szamel, “Long-range spatial correlations of particle displacements and the emergence of elasticity,” Physical Review Letters **114**, 025501 (2015).

[30] L Klochko, J Baschnagel, JP Wittmer, H Meyer, O Benzerara, and AN Semenov, “Theory of length-scale dependent relaxation moduli and stress fluctuations in glass-forming and viscoelastic liquids,” The Journal of Chemical Physics **156**, 164505 (2022).

[31] Rahul N Chacko, François P Landes, Giulio Biroli, Olivier Dauchot, Andrea J Liu, and David R Reichman, “Elastoplasticity mediates dynamical heterogeneity below the mode coupling temperature,” Physical Review Letters **127**, 048002 (2021).

[32] Matthias Lerbinger, Armand Barbot, Damien Vandembroucq, and Sylvain Patinet, “Relevance of shear transformations in the relaxation of supercooled liquids,” Physical Review Letters **129**, 195501 (2022).

[33] Michael L Falk and James S Langer, “Dynamics of viscoplastic deformation in amorphous solids,” Physical Review E **57**, 7192 (1998).

[34] Misaki Ozawa and Giulio Biroli, “Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids,” Physical Review Letters **130**, 138201 (2023).

[35] Guillemette Picard, Armand Ajdari, François Lequeux, and Lydéric Bocquet, “Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids,” The European Physical Journal E **15**, 371–381 (2004).

[36] A. Nicolas, E.E. Ferrero, K. Martens, and J.-L. Barrat, “Deformation and flow of amorphous solids: Insights from elastoplastic models,” Rev. Mod. Phys. **90**, 045006 (2018).

[37] D. Vandembroucq, R. Skoe, and S. Roux, “Universal depinning force fluctuations of an elastic line: Application to finite temperature behavior,” Phys. Rev. E **70**, 051101 (2004).

[38] J. Lin, E. Lerner, A. Rosso, and M. Wyart, “Scaling description of the yielding transition in soft amorphous solids at zero temperature,” Proc. Natl. Acad. Sci. **111**, 14382–14387 (2014).

[39] Saverio Rossi, Giulio Biroli, Misaki Ozawa, Gilles Tarjus, and Francesco Zamponi, “Finite-disorder critical point in the yielding transition of elastoplastic models,” Physical Review Letters **129**, 228002 (2022).

[40] Ludovic Berthier, Giulio Biroli, J-P Bouchaud, Luca Cipelletti, D El Masri, Denis L'Hôte, Francois Ladieu, and Matteo Pierno, “Direct experimental evidence of a growing length scale accompanying the glass transition,” Science **310**, 1797–1800 (2005).

[41] M. Popović, T.W.J. de Geus, W. Ji, and M. Wyart, “Thermally activated flow in models of amorphous solids,” Phys. Rev. E **104**, 025010 (2021).

[42] Raphaël Candelier, Asaph Widmer-Cooper, Jonathan K Kummerfeld, Olivier Dauchot, Giulio Biroli, Peter Harrowell, and David R Reichman, “Spatiotemporal hierarchy of relaxation events, dynamical heterogeneities, and structural reorganization in a supercooled liquid,” Physical review letters **105**, 135702 (2010).

[43] Aaron S Keys, Lester O Hedges, Juan P Garrahan, Sharon C Glotzer, and David Chandler, “Excitations are localized and relaxation is hierarchical in glass-forming liquids,” Physical Review X **1**, 021013 (2011).

[44] Taiki Yanagishima, John Russo, and Hajime Tanaka, “Common mechanism of thermodynamic and mechanical origin for ageing and crystallization of glasses,” Nature communications **8**, 1–10 (2017).

[45] James P Sethna, Karin A Dahmen, and Christopher R Myers, “Crackling noise,” Nature **410**, 242–250 (2001).

[46] Alberto Rosso, James P Sethna, and Matthieu Wyart, “Avalanches and deformation in glasses and disordered systems,” arXiv preprint arXiv:2208.04090 (2022).

[47] Víctor Hugo Purrello, Jose Luis Iguain, Alejandro Benedykt Kolton, and Eduardo Alberto Jagla, “Creep and thermal rounding close to the elastic depinning threshold,” Physical Review E **96**, 022112 (2017).

[48] Liheng Yao and Robert L Jack, “Thermal vestiges of avalanches in the driven random field ising model,” Journal of Statistical Mechanics: Theory and Experiment **2023**, 023303 (2023).

[49] Daniel Korchinski and Jörg Rottler, “Dynamic phase diagram of plastically deformed amorphous solids at finite temperature,” Physical Review E **106**, 034103 (2022).

[50] Camille Scalliet, Benjamin Guiselin, and Ludovic Berthier, “Excess wings and asymmetric relaxation spectra in a facilitated trap model,” The Journal of Chemical Physics **155**, 064505 (2021).

[51] Camille Scalliet, Benjamin Guiselin, and Ludovic Berthier, “Thirty milliseconds in the life of a supercooled liquid,” Physical Review X **12**, 041028 (2022).

[52] Gilles Tarjus and Daniel Kivelson, “Breakdown of the stokes–einstein relation in supercooled liquids,” The Journal of chemical physics **103**, 3071–3073 (1995).

[53] Mark D Ediger, “Spatially heterogeneous dynamics in supercooled liquids,” Annual review of physical chemistry **51**, 99–128 (2000).

[54] Shiladitya Sengupta, Smarajit Karmakar, Chandan Dasgupta, and Srikanth Sastry, “Breakdown of the stokes-einstein relation in two, three, and four dimensions,” The Journal of chemical physics **138**, 12A548 (2013).

[55] Patrick Charbonneau, Yuliang Jin, Giorgio Parisi, and Francesco Zamponi, “Hopping and the stokes-einstein relation breakdown in simple glass formers,” Proceedings of the National Academy of Sciences **111**, 15025–15030 (2014).

[56] Takeshi Kawasaki and Kang Kim, “Identifying time scales for violation/preservation of stokes-einstein relation in supercooled water,” Science Advances **3**, e1700399 (2017).
[57] YounJoon Jung, Juan P Garrahan, and David Chandler, “Excitation lines and the breakdown of stokes-einstein relations in supercooled liquids,” Physical Review E 69, 061205 (2004).

[58] Ludovic Berthier, David Chandler, and Juan P Garrahan, “Length scale for the onset of fickian diffusion in supercooled liquids,” Europhysics Letters 69, 320 (2004).

[59] Lester O Hedges, Lutz Maibaum, David Chandler, and Juan P Garrahan, “Decoupling of exchange and persistence times in atomistic models of glass formers,” The Journal of chemical physics 127, 211101 (2007).

[60] Pinaki Chaudhuri, Ludovic Berthier, and Walter Kob, “Universal nature of particle displacements close to glass and jamming transitions,” Physical review letters 99, 060604 (2007).

[61] Raffaele Pastore, Takuma Kikutsuji, Francesco Rusciano, Nobuyuki Matubayasi, Kang Kim, and Francesco Greco, “Breakdown of the stokes-einstein relation in supercooled liquids: A cage-jump perspective,” The Journal of Chemical Physics 155, 114503 (2021).

[62] VV Bulatov and AS Argon, “A stochastic model for continuum elasto-plastic behavior. ii. a study of the glass transition and structural relaxation,” Modelling and Simulation in Materials Science and Engineering 2, 185 (1994).

[63] Ezequiel E Ferrero, Kirsten Martens, and Jean-Louis Barrat, “Relaxation in yield stress systems through elastically interacting activated events,” Physical review letters 113, 248301 (2014).

[64] Eduardo Alberto Jagla, “Tensorial description of the plasticity of amorphous composites,” Physical Review E 101, 043004 (2020).

[65] Ezequiel E Ferrero, Alejandro B Kolton, and Eduardo A Jagla, “Yielding of amorphous solids at finite temperatures,” Physical Review Materials 5, 115602 (2021).

[66] Gieberth Rodriguez-Lopez, Kirsten Martens, and Ezequiel E Ferrero, “Temperature dependence of fast relaxation processes in amorphous materials,” arXiv preprint arXiv:2302.06471 (2023).

[67] Yue Fan, Takuya Iwashita, and Takeshi Egami, “How thermally activated deformation starts in metallic glass,” Nature communications 5, 5083 (2014).

[68] I Fernández Aguirre and Eduardo Alberto Jagla, “Critical exponents of the yielding transition of amorphous solids,” Physical Review E 98, 013002 (2018).

[69] Craig E Maloney and Anael Lemaitre, “Amorphous systems in athermal, quasistatic shear,” Physical Review E 74, 016118 (2006).

[70] Full isotropy is slightly broken by the choice of biperiodic boundary conditions.

[71] Claudio Donati, Silvio Franz, Sharon C Glotzer, and Giorgio Parisi, “Theory of non-linear susceptibility and correlation length in glasses and liquids,” Journal of non-crystalline solids 307, 215–224 (2002).

[72] Simone Capaccioli, Giancarlo Ruocco, and Francesco Zamponi, “Dynamically correlated regions and configurational entropy in supercooled liquids,” The Journal of Physical Chemistry B 112, 10652–10658 (2008).

[73] Olivier Dauchot, François Ladieu, and C Patrick Royall, “The glass transition in molecules, colloids and grains: universality and specificity,” arXiv preprint arXiv:2211.03158 (2022).

[74] Smarajit Karmakar, Chandan Dasgupta, and Srikanth Sastry, “Growing length and time scales in glass-forming liquids,” Proceedings of the National Academy of Sciences 106, 3675–3679 (2009).

[75] Daniele Coslovich, Misaki Ozawa, and Walter Kob, “Dynamic and thermodynamic crossover scenarios in the kob-andersen mixture: Insights from multi-cpu and multi-gpu simulations,” The European Physical Journal E 41, 1–11 (2018).

[76] Saurish Chakrabarty, Indrajit Tah, Smarajit Karmakar, and Chandan Dasgupta, “Block analysis for the calculation of dynamic and static length scales in glass-forming liquids,” Physical Review Letters 119, 205502 (2017).

[77] N Lačević, Francis W Starr, TB Schröder, and Sharon C Glotzer, “Spatially heterogeneous dynamics investigated via a time-dependent four-point density correlation function,” The Journal of chemical physics 119, 7372–7387 (2003).

[78] We are considering either finite systems or the thermodynamic limit taken after the zero-temperature limit.

[79] Maya Paczuski, Sergei Maslov, and Per Bak, “Avalanche dynamics in evolution, growth, and depinning models,” Physical Review E 53, 414 (1996).

[80] J.-C. Baret, D. Vandembroucq, and S. Roux, “Extremal model for amorphous media plasticity,” Phys. Rev. Lett. 89, 195506 (2002).

[81] Dheeraj Kumar, Sylvain Patinet, Craig E Maloney, Ido Regev, Damien Vandembroucq, and Muhittin Mungan, “Mapping out the glassy landscape of a mesoscopic elastoplastic model,” The Journal of Chemical Physics 157, 174504 (2022).

[82] Jihui Han, Wei Li, and Weibing Deng, “Critical behavior of a generalized bak-sneppen model,” in Journal of Physics: Conference Series, Vol. 1113 (IOP Publishing, 2018) p. 012011.

[83] J. Lin, A. Saade, E. Lerner, A. Rosso, and M. Wyart, “On the density of shear transformations in amorphous solids,” Europhys. Lett. 105, 26003 (2014).

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

[85] Markus Müller and Matthieu Wyart, “Marginal stability in structural, spin, and electron glasses,” Annu. Rev. Condens. Matter Phys. 6, 177–200 (2015).

[86] A. Lemaître and C. Caroli, “Plastic Response of a 2D Amorphous Solid to Quasi-Static Shear: II - Dynamical Noise and Avalanches in a Mean Field Model,” arXiv preprint 0705.3122 (2007), 10.48550/arXiv:0705.3122.

[87] Jie Lin and Matthieu Wyart, “Mean-field description of plastic flow in amorphous solids,” Physical review X 6, 011005 (2016).

[88] Ezequiel E Ferrero and Eduardo A Jagla, “Properties of the density of shear transformations in driven amorphous solids,” Journal of Physics: Condensed Matter 33, 124001 (2021).

[89] Daniel Korchinski, Céline Ruscher, and Jörg Rottler, “Signatures of the spatial extent of plastic events in the yielding transition in amorphous solids,” Physical Review E 104, 034603 (2021).

[90] Massimo Pica Ciamarra, Raffaele Pastore, and Antonio Coniglio, “Particle jumps in structural glasses,” Soft matter 12, 358–366 (2016).

[91] Stephen F Swallen, Paul A Bonvallet, Robert J McMahon, and MD Ediger, “Self-diffusion of tris-naphthylbenzene near the glass transition temperature,
ture,” Physical review letters **90**, 015901 (2003).

[92] Francesco Mallamace, Caterina Branca, Carmelo Cor- saro, Nancy Leone, Jeroen Spooren, Sow-Hsin Chen, and H Eugene Stanley, “Transport properties of glass- forming liquids suggest that dynamic crossover temper- ature is as important as the glass transition tempera- ture,” Proceedings of the National Academy of Sciences **107**, 22457–22462 (2010).

[93] In fact, in some kinetically constrained models [140], the fractional Stokes-Einstein violation (with a similar value of $\zeta$ with experiments) was conjectured based on numerical data and physical arguments. Recent rigor- ous results [141] showed that this is not what happens, but there is likely a violation as the one we find. This suggests that the logarithmic Stokes-Einstein violation can be easily misinterpreted as a fractional one with a small exponent $1 - \zeta$.

[94] Misaki Ozawa, Kang Kim, and Kunimasa Miyazaki, “Tuning pairwise potential can control the fragility of glass-forming liquids: From a tetrahedral network to isotropic soft sphere models,” Journal of Statistical Me- chanics: Theory and Experiment **2016**, 074002 (2016).

[95] Sivert H Glarum, “Dielectric relaxation of isoamyl bro- mide,” The Journal of Chemical Physics **33**, 639–643 (1960).

[96] Muhammad R Hasyim and Kranthi K Mandadapu, “A theory of localized excitations in supercooled liquids,” The Journal of chemical physics **155**, 044504 (2021).

[97] Masaharu Isobe, Aaron S Keys, David Chandler, and Juan P Garrahan, “Applicability of dynamic facilitation theory to binary hard disk systems,” Physical review letters **117**, 145701 (2016).

[98] Yael S Elmatad, Robert L Jack, David Chandler, and Juan P Garrahan, “Finite-temperature critical point of a glass transition,” Proceedings of the National Academy of Sciences **107**, 12793–12798 (2010).

[99] Francesco Turci, C Patrick Royall, and Thomas Speck, “Nonequilibrium phase transition in an atomistic glass- former: The connection to thermodynamics,” Physical Review X **7**, 031028 (2017).

[100] Elisabeth Agoritsas, Eric Bertin, Kirsten Martens, and Jean-Louis Barrat, “On the relevance of disorder in athermal amorphous materials under shear,” The Eu- ropean Physical Journal E **38**, 1–22 (2015).

[101] H. Tanaka, T. Kawasaki, H. Shintani, and K. Watan- abe, “Critical-like behaviour of glass-forming liquids,” Nat. Mater. **9**, 324–331 (2010).

[102] C Patrick Royall and Stephen R Williams, “The role of local structure in dynamical arrest,” Physics Reports **560**, 1–75 (2015).

[103] Hajime Tanaka, Hua Tong, Rui Shi, and John Russo, “Revealing key structural features hidden in liquids and glasses,” Nature Reviews Physics **1**, 333–348 (2019).

[104] Joris Paret, Robert L Jack, and Daniele Coslovich, “Assessing the structural heterogeneity of supercooled liquids through community inference,” The Journal of chemical physics **152**, 144502 (2020).

[105] Asaph Widmer-Cooper, Peter Harrowell, and H Fynew- ever, “How reproducible are dynamic heterogeneities in a supercooled liquid?” Physical review letters **93**, 135701 (2004).

[106] A. Widmer-Cooper, H. Perry, P. Harrowell, and D.R. Reichman, “Irreversible reorganization in a supercooled liquid originates from localized soft modes,” Nat. Phys.

4, 711–715 (2008).

[107] Glen M Hocky, Daniele Coslovich, Atsushi Ikeda, and David R Reichman, “Correlation of local order with par- ticle mobility in supercooled liquids is highly system de- pendent,” Physical review letters **113**, 157801 (2014).

[108] Anaël Lemaître, Chandana Mondal, Michael Moshe, Itamar Procaccia, Saikat Roy, and Keren Screiber- Re'em, “Anomalous elasticity and plastic screening in amorphous solids,” Physical Review E **104**, 024904 (2021).

[109] E. Lerner, E. DeGiuli, G. Düring, and M. Wyart, “Breakdown of continuum elasticity in amorphous solids,” Soft Matter **10**, 5085 (2014).

[110] Chen Liu, Suman Dutta, Pinaki Chaudhuri, and Kirsten Martens, “Elastoplastic approach based on mi- croscopic insights for the steady state and transient dy- namics of sheared disordered solids,” Physical Review Letters **126**, 138005 (2021).

[111] David Fernández Castellanos, Stéphane Roux, and Syl- vain Patinet, “Insights from the quantitative calibration of an elasto-plastic model from a lennard-jones atomic glass,” Comptes Rendus. Physique **22**, 1–28 (2021).

[112] David F Castellanos, Stéphane Roux, and Sylvain Patinet, “History dependent plasticity of glass: A map- ping between atomistic and elasto-plastic models,” Acta Materialia **241**, 118405 (2022).

[113] Indrajit Tah, Sean A Ridout, and Andrea J Liu, “Fragility in glassy liquids: A structural approach based on machine learning,” The Journal of Chemical Physics **157**, 124501 (2022).

[114] Hongyi Xiao, Ge Zhang, Entao Yang, Robert JS Ivan- cic, Sean A Ridout, Robert Riggleman, Douglas J Durian, and Andrea J Liu, “Machine learning-informed structuro-elastoplasticity predicts ductility of disor- dered solids,” arXiv preprint arXiv:2303.12486 (2023).

[115] Tina Hecksher and Jeppe C Dyre, “A review of ex- periments testing the shoving model,” Journal of Non- Crystalline Solids **407**, 14–22 (2015).

[116] Wencheng Ji, Tom WJ de Geus, Elisabeth Agoritsas, and Matthieu Wyart, “Mean-field description for the architecture of low-energy excitations in glasses,” Phys- ical Review E **105**, 044601 (2022).

[117] Camille Scalliet, Ludovic Berthier, and Francesco Zam- poni, “Nature of excitations and defects in structural glasses,” Nature communications **10**, 5102 (2019).

[118] Hideyuki Mizuno, Hua Tong, Atsushi Ikeda, and Ste- fano Mossa, “Intermittent rearrangements accompany- ing thermal fluctuations distinguish glasses from crys- tals,” The Journal of Chemical Physics **153**, 154501 (2020).

[119] David Richard, Geert Kapteijns, and Edan Lerner, “Detecting low-energy quasilocalized excitations in computer glasses,” arXiv preprint arXiv:2303.12887 (2023).

[120] Daniele Coslovich and Giorgio Pastore, “Understanding fragility in supercooled lennard-jones mixtures. i. locally preferred structures,” The Journal of chemical physics **127**, 124504 (2007).

[121] Giulio Biroli and Jean-Philippe Bouchaud, “The rfot theory of glasses: Recent progress and open issues,” arXiv preprint arXiv:2208.05866 (2022).

[122] David Fernandez Castellanos and Michael Zaiser, “Avalanche behavior in creep failure of disordered ma- terials,” Physical review letters **121**, 125501 (2018).
[123] T. Bauer, J. Oberdisse, and L. Ramos, “Collective re-
arrangement at the onset of flow of a polycrystalline
hexagonal columnar phase,” *Phys. Rev. Lett.* **97** (2006),
10.1103/PhysRevLett.97.258303.

[124] F. Caton and C. Baravian, “Plastic behavior of some
yield stress fluids: from creep to long-time yield,” *Rheol.
Acta* **47**, 601–607 (2008).

[125] T. Divoux, C. Barentin, and S. Manneville, “From stress-induced fluidization processes to herschel-bulkley behaviour in simple yield stress fluids,” *Soft Matter* **7**, 8409 (2011).

[126] M. Siebenbürger, M. Ballauff, and T. Voigtmann, “Creep in colloidal glasses,” *Phys. Rev. Lett.* **108**, 255701 (2012).

[127] V. Grenard, T. Divoux, N. Taberlet, and S. Manneville, “Timescales in creep and yielding of attractive gels,” *Soft Matter*, **17** (2014).

[128] M. Leocmach, C. Perge, T. Divoux, and S. Man-
neville, “Creep and fracture of a protein gel under
stress,” *Phys. Rev. Lett.* **113** (2014), 10.1103/Phys-
RevLett.113.038303.

[129] Sebastian Bustingorry, AB Kolton, and Thierry Gi-
amarchi, “Thermal rounding of the depinning transi-
tion,” Europhysics Letters **81**, 26005 (2007).

[130] Alejandro B Kolton, Alberto Rosso, and Thierry Gia-
marchi, “Creep motion of an elastic string in a random
potential,” Physical review letters **94**, 047002 (2005).

[131] Ludovic Berthier and Walter Kob, “The monte carlo
dynamics of a binary lennard-jones glass-forming mix-
ture,” *Journal of Physics: Condensed Matter* **19**, 205130
(2007).

[132] Craig E Maloney and Daniel J Lacks, “Energy bar-
rier scalings in driven systems,” Physical Review E **73**,
061106 (2006).

[133] Armand Barbot, Matthias Lerbinger, Anier Hernandez-

Garcia, Reinaldo García-García, Michael L Falk,
Damien Vandembroucq, and Sylvain Patinet, “Local
yield stress statistics in model amorphous solids,” Phys-
ical Review E **97**, 033001 (2018).

[134] Marko Popović, Tom WJ de Geus, and Matthieu Wyart, “Elastoplastic description of sudden failure in athermal amorphous materials during quasistatic loading,” Physical Review E **98**, 040901 (2018).

[135] Joseph Pollard and Suzanne M Fielding, “Yielding, shear banding, and brittle failure of amorphous materials,” Physical Review Research **4**, 043037 (2022).

[136] Elijah Flenner, Hannah Staley, and Grzegorz Szamel, “Universal features of dynamic heterogeneity in supercooled liquids,” Physical review letters **112**, 097801 (2014).

[137] Kang Kim and Shinji Saito, “Multiple length and time scales of dynamic heterogeneities in model glass-forming liquids: A systematic analysis of multi-point and multi-time correlations,” *The Journal of chemical physics* **138**, 12A506 (2013).

[138] Giulio Biroli, Kunimasa Miyazaki, and David R Re-
ichman, “Dynamical heterogeneity in glass-forming liq-
uids,” arXiv preprint arXiv:2209.02825 (2022).

[139] Elijah Flenner and Grzegorz Szamel, “Dynamic heterogeneity in two-dimensional supercooled liquids: Comparison of bond-breaking and bond-orientational correlations,” *Journal of Statistical Mechanics: Theory and Experiment* **2016**, 074008 (2016).

[140] Juan P Garrahan, Peter Sollich, and Cristina Toninelli, “Kinetically constrained models,” Dynamical heterogeneities in glasses, colloids, and granular media **150**, 111–137 (2011).

[141] Oriane Blondel and Cristina Toninelli, “Is there a frac-
tional breakdown of the stokes-einstein relation in ki-
netically constrained models at low temperature?” Eu-
rophysics Letters **107**, 26005 (2014).