Ali Tahaei

Giulio Biroli

Misaki Ozawa

Marko Popović

Matthieu Wyart

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

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

As a liquid is cooled, the relaxation time τ α 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 [START_REF] Anderson | Through the glass lightly[END_REF][START_REF] Mark D Ediger | Supercooled liquids and glasses[END_REF][START_REF] Debenedetti | Supercooled liquids and the glass transition[END_REF][START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF]. In this regime the effective activation energy associated to τ α grows for many liquids, leading to a super-Arrhenius behavior. Approaching T g , dynamics also becomes heterogeneous on a growing correlation length scale ξ [START_REF] Kob | Dynamical heterogeneities in a supercooled lennard-jones liquid[END_REF][START_REF] Yamamoto | Dynamics of highly supercooled liquids: Heterogeneity, rheology, and diffusion[END_REF][START_REF] Dalle-Ferrier | Spatial correlations in the dynamics of glassforming liquids: Experimental determination of their temperature dependence[END_REF][START_REF] Karmakar | Growing length scales and their relation to timescales in glass-forming liquids[END_REF]. 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 [START_REF] Theodore R Kirkpatrick | Scaling concepts for the dynamics of viscous liquids near an ideal glassy state[END_REF][START_REF] Lubchenko | Theory of structural glasses and supercooled liquids[END_REF][START_REF] Biroli | The random first-order transition theory of glasses: a critical assessment[END_REF] 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 [START_REF] Ritort | Glassy dynamics of kinetically constrained models[END_REF][START_REF] Berthier | Dynamical heterogeneities in glasses, colloids, and granular media[END_REF]. This lead to a theory [START_REF] Juan | Geometrical explanation and scaling of dynamical heterogeneities in glass forming systems[END_REF][START_REF] Juan P Garrahan | Dynamical first-order phase transition in kinetically constrained models of glasses[END_REF][START_REF] Hedges | Dynamic order-disorder in atomistic models of structural glass formers[END_REF] 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 ξ. Free volume [START_REF] Turnbull | Free-volume model of the amorphous phase: glass transition[END_REF] or elastic [START_REF] Dyre | Colloquium: The glass transition and elastic models of glass-forming liquids[END_REF][START_REF] Rainone | Pinching a glass reveals key properties of its soft spots[END_REF][START_REF] Kapteijns | Does mesoscopic elasticity control viscous slowing down in glassforming liquids?[END_REF] 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 [START_REF] Pica Ciamarra | The energy cost of local rearrangements, not cooperative effects, makes glasses solid[END_REF]. 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 [START_REF] Dyre | Colloquium: The glass transition and elastic models of glass-forming liquids[END_REF][START_REF] Thomas | Solid-like mean-square displacement in glass-forming liquids[END_REF]: a rearrangement corresponds to a local (plastic) strain, that affects stress away from it. Lemaitre [START_REF] Lemaître | Structural Relaxation is a Scale-Free Process[END_REF] 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 [START_REF] Chowdhury | Long range stress correlations in the inherent structures of liquids at rest[END_REF][START_REF] Tong | Emergent solidity of amorphous materials as a consequence of mechanical self-organisation[END_REF][START_REF] Wu | Anisotropic stress correlations in two-dimensional liquids[END_REF][START_REF] Maier | Emergence of long-ranged stress correlations at the liquid to glass transition[END_REF][START_REF] Steffen | Molecular simulations and hydrodynamic theory of nonlocal shear-stress correlations in supercooled fluids[END_REF][START_REF] Flenner | Long-range spatial correlations of particle displacements and the emergence of elasticity[END_REF][START_REF] Klochko | Theory of length-scale dependent relaxation moduli and stress fluctuations in glass-forming and viscoelastic liquids[END_REF]. Recently, it was shown that indeed elastic interactions play an important role in dynamical facilitation [START_REF] Rahul N Chacko | Elastoplasticity mediates dynamical heterogeneity below the mode coupling temperature[END_REF]. Moreover, Ref. [START_REF] Lerbinger | Relevance of shear transformations in the relaxation of supercooled liquids[END_REF] showed that heterogeneous relaxations take place in the form of plastic rearrangements that are called shear transformations [START_REF] Michael | Dynamics of viscoplastic deformation in amorphous solids[END_REF]. Following these molecular simulations studies, Ref. [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF] demonstrated with elasto-plastic models (traditionally used to study the plasticity of amorphous solids under loading [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF][START_REF] Nicolas | Deformation and flow of amorphous solids: Insights from elastoplastic models[END_REF][START_REF] Vandembroucq | Universal depinning force fluctuations of an elastic line: Application to finite temperature behavior[END_REF][START_REF] Lin | Scaling description of the yielding transition in soft amorphous solids at zero temperature[END_REF][START_REF] Rossi | Finite-disorder critical point in the yielding transition of elastoplastic models[END_REF]) 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 [START_REF] Berthier | Dynamical heterogeneities in glasses, colloids, and granular media[END_REF][START_REF] Berthier | Direct experimental evidence of a growing length scale accompanying the glass transition[END_REF]. Furthermore, such models also capture the emergence of a gap in the distribution of local barriers under cooling, which controls the dynamics [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF][START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF]. 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 [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF] 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 [START_REF] Candelier | Spatiotemporal hierarchy of relaxation events, dynamical heterogeneities, and structural reorganization in a supercooled liquid[END_REF][START_REF] Keys | Excitations are localized and relaxation is hierarchical in glassforming liquids[END_REF][START_REF] Yanagishima | Common mechanism of thermodynamic and mechanical origin for ageing and crystallization of glasses[END_REF]. Our analysis builds a link with systems that crackle such as disordered magnets, granular materials, and earthquakes [START_REF] James P Sethna | Crackling noise[END_REF][START_REF] Rosso | Avalanches and deformation in glasses and disordered systems[END_REF] even at finite temperatures [START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF][START_REF] Víctor | Creep and thermal rounding close to the elastic depinning threshold[END_REF][START_REF] Yao | Thermal vestiges of avalanches in the driven random field ising model[END_REF][START_REF] Korchinski | Dynamic phase diagram of plastically deformed amorphous solids at finite temperature[END_REF], revealing that dynamical heterogeneities are controlled by a critical point at zero temperature, with a diverging length scale ξ ∼ T -ν and correlation volume χ * 4 ∼ T -γ . We provide a scaling argument expressing ν and γ 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 [START_REF] Scalliet | Excess wings and asymmetric relaxation spectra in a facilitated trap model[END_REF][START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF]. Based on the properties of the zerotemperature fixed point governing thermal avalanches, we also show that the decoupling Dτ α between particle diffusion D and relaxation time τ α (the so-called Stoke-Einstein violation- [START_REF] Tarjus | Breakdown of the stokes-einstein relation in supercooled liquids[END_REF][START_REF] Mark | Spatially heterogeneous dynamics in supercooled liquids[END_REF][START_REF] Sengupta | Breakdown of the stokes-einstein relation in two, three, and four dimensions[END_REF][START_REF] Charbonneau | Hopping and the stokes-einstein relation breakdown in simple glass formers[END_REF][START_REF] Kawasaki | Identifying time scales for violation/preservation of stokes-einstein relation in supercooled water[END_REF]) diverges as Dτ α ∼ 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 τ α [START_REF] Jung | Excitation lines and the breakdown of stokeseinstein relations in supercooled liquids[END_REF][START_REF] Berthier | Length scale for the onset of fickian diffusion in supercooled liquids[END_REF][START_REF] Lester O Hedges | Decoupling of exchange and persistence times in atomistic models of glass formers[END_REF][START_REF] Chaudhuri | Universal nature of particle displacements close to glass and jamming transitions[END_REF][START_REF] Pastore | Breakdown of the stokes-einstein relation in supercooled liquids: A cage-jump perspective[END_REF]. 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 lowtemperature glass-forming liquids, we employ a generalization of elasto-plastic models [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF][START_REF] Nicolas | Deformation and flow of amorphous solids: Insights from elastoplastic models[END_REF] 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 [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF][START_REF] Nicolas | Deformation and flow of amorphous solids: Insights from elastoplastic models[END_REF][START_REF] Lin | Scaling description of the yielding transition in soft amorphous solids at zero temperature[END_REF][START_REF] Rossi | Finite-disorder critical point in the yielding transition of elastoplastic models[END_REF] (some exceptions e.g., Ref. [START_REF] Vv Bulatov | A stochastic model for continuum elasto-plastic behavior. ii. a study of the glass transition and structural relaxation[END_REF]) and they were generalized to take into account thermal fluctuations [START_REF] Ferrero | Relaxation in yield stress systems through elastically interacting activated events[END_REF] and then to study glass-forming liquids [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF]. 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 [START_REF] Dyre | Colloquium: The glass transition and elastic models of glass-forming liquids[END_REF][START_REF] Lemaître | Structural Relaxation is a Scale-Free Process[END_REF], 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 [START_REF] Eduardo | Tensorial description of the plasticity of amorphous composites[END_REF] and a scalar elasto-plastic model in which the shear stress is represented by a scalar variable [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF] (see Appendix A for details). In both models, thermal fluctuations are implemented through a probability rate τ -1 0 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 τ 0 is a microscopic relaxation timescale [START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF][START_REF] Ferrero | Relaxation in yield stress systems through elastically interacting activated events[END_REF][START_REF] Ferrero | Yielding of amorphous solids at finite temperatures[END_REF][START_REF] Rodriguez-Lopez | Temperature dependence of fast relaxation processes in amorphous materials[END_REF]. 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 α , which is suggested by recent elasto-plastic models and molecular simulations [START_REF] Lerbinger | Relevance of shear transformations in the relaxation of supercooled liquids[END_REF][START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF][START_REF] Ferrero | Yielding of amorphous solids at finite temperatures[END_REF][START_REF] Rodriguez-Lopez | Temperature dependence of fast relaxation processes in amorphous materials[END_REF][START_REF] Yue Fan | How thermally activated deformation starts in metallic glass[END_REF]. In this study, we set τ 0 = 1, c = 1, and the value α = 3/2 corresponds to a smooth microscopic potential [START_REF] Aguirre | Critical exponents of the yielding transition of amorphous solids[END_REF]. As stress is relaxed locally at the site of the plastic event, stress in the system changes according to the Eshelby kernel [START_REF] Picard | Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids[END_REF]. 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 [START_REF] Craig | Amorphous systems in athermal, quasistatic shear[END_REF]. 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 [START_REF]Full isotropy is slightly broken by the choice of biperiodic boundary conditions[END_REF]. 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, ⟨π(t)⟩ time , which has been used widely in the context of kinetically constrained models [START_REF] Ritort | Glassy dynamics of kinetically constrained models[END_REF]. The observable π(t) is defined by π(t) = 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 ⟨• • • ⟩ 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 ⟨π(t)⟩ time at different temperatures, as shown in Fig. 1(a). In Fig. 1(b), we find that the associated relaxation time τ α defined by ⟨π(τ α )⟩ time = 1/2 increases in a Arrhenius way. Figure 2 snapshot of local persistence for τ /τ α ≃ 0.53, demonstrating spatially heterogeneous dynamics. To quantify the magnitude of dynamical heterogeneity, we compute the four-point correlations function [START_REF] Donati | Theory of non-linear susceptibility and correlation length in glasses and liquids[END_REF], χ 4 (t), defined by

χ 4 (t) = N ⟨π 2 (t)⟩ time -⟨π(t)⟩ 2 time . (1)
χ 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 [START_REF] Berthier | Dynamical heterogeneities in glasses, colloids, and granular media[END_REF], as it has been estimated in real experiments [START_REF] Dalle-Ferrier | Spatial correlations in the dynamics of glassforming liquids: Experimental determination of their temperature dependence[END_REF][START_REF] Berthier | Direct experimental evidence of a growing length scale accompanying the glass transition[END_REF][START_REF] Capaccioli | Dynamically correlated regions and configurational entropy in supercooled liquids[END_REF][START_REF] Dauchot | The glass transition in molecules, colloids and grains: universality and specificity[END_REF] as well as molecular simulations [START_REF] Donati | Theory of non-linear susceptibility and correlation length in glasses and liquids[END_REF][START_REF] Karmakar | Growing length and time scales in glass-forming liquids[END_REF][START_REF] Coslovich | Dynamic and thermodynamic crossover scenarios in the kob-andersen mixture: Insights from multi-cpu and multi-gpu simulations[END_REF]. Figure 2(b) shows the time and temperature evolution of χ 4 (t). It takes a peak near the relaxation timescale τ α , 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 χ 4 (t), denoted χ * 4 , versus temperature T for several system sizes in Fig. 3(a). The system size dependence, akin to molec- ular simulations [START_REF] Karmakar | Growing length and time scales in glass-forming liquids[END_REF][START_REF] Chakrabarty | Block analysis for the calculation of dynamic and static length scales in glass-forming liquids[END_REF], allowed us to perform finite size scaling. We find that for increasing the system size L, χ * 4 follows a scaling form, χ * 4 ∼ T -γ with a critical point at T = 0 and the associated exponent γ. We obtain a scaling collapse for χ * 4 (L, T ) in Fig 3(b) (see also molecular simulation studies [START_REF] Karmakar | Growing length and time scales in glass-forming liquids[END_REF][START_REF] Chakrabarty | Block analysis for the calculation of dynamic and static length scales in glass-forming liquids[END_REF]), indicating that the dynamics is governed by a diverging correlation length toward T = 0, i.e., ξ ∼ T -ν with an exponent ν. The corresponding data for the scalar model are presented in Appendix B. Besides, we obtained a consistent result in terms of ν by directly measuring the dynamical correlation lengthscale based on the four-point structure factor [START_REF] Lačević | Spatially heterogeneous dynamics investigated via a time-dependent four-point density correlation function[END_REF]. 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. 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 3/2
min , 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 selforganized-criticality [START_REF] Paczuski | Avalanche dynamics in evolution, growth, and depinning models[END_REF] and some disordered materials under quasi-static driving [START_REF] Víctor | Creep and thermal rounding close to the elastic depinning threshold[END_REF][START_REF] Baret | Extremal model for amorphous media plasticity[END_REF], including annealed glasses [START_REF] Kumar | Mapping out the glassy landscape of a mesoscopic elastoplastic model[END_REF].

In the thermodynamic limit, L → ∞, the extremal dynamics leads to an absorbing condition for E ≤ E c , where E c is a critical energy barrier as found (for T → 0) in Ref. [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF]. As a result, the distribution of local (activation) energy barriers P (E) vanishes for E ≤ E c (see the sketch in Fig. 9(a)). This implies that lim L→∞ P (x) = 0 for x ≤ x c = E 2/3 c , where P (x) is the distribution of x. Thus only a sub-extensive number of sites are found for x ≤ 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 [START_REF] Paczuski | Avalanche dynamics in evolution, growth, and depinning models[END_REF]. 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 [START_REF] Víctor | Creep and thermal rounding close to the elastic depinning threshold[END_REF][START_REF] Paczuski | Avalanche dynamics in evolution, growth, and depinning models[END_REF]. For a finite L, we fix a chosen threshold stability value x 0 (≤ 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 (eventbased avalanche size), and the total number of sites S that relaxed at least once during an avalanche (site-based avalanche size). By construction, S ≤ S and S ≤ L d . We show snapshots having S and the corresponding 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 [START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF]), whereas the sitebased avalanche size 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 χ 4 ). Sizes of avalanches, S and 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 [START_REF] Víctor | Creep and thermal rounding close to the elastic depinning threshold[END_REF][START_REF] Paczuski | Avalanche dynamics in evolution, growth, and depinning models[END_REF][START_REF] Han | Critical behavior of a generalized bak-sneppen model[END_REF]:

P (S) ∼ S -τ g(S/S c ), (2)
where g(z) is a scaling function and S c is a cutoff size which takes the form:

S c ∼ (x c -x 0 ) -1/σ f L d f (x c -x 0 ) -1/σ , (3)
where 1/σ and d f are critical exponents and f (z) = 1 for z ≫ 1 and f (z) = z for z → 0. Thus,

S c ∼ L d f when x 0 → x c , whereas S c ∼ (x c -x 0 ) -1/σ when L → ∞.
The same expressions are expected to hold for the sitebased avalanche size S, defining exponents τ , 1/σ, and df . To estimate S c we use the fact that for 1 < τ < 3 the ratio ⟨S 3 ⟩/⟨S 2 ⟩ is proportional to the cutoff value S c , where

⟨• • • ⟩ = ∞ 0 dSP (S)(• • • ).
As we are interested in scaling of S c with system size, the numerical constant is irrelevant and we define S c ≡ ⟨S 3 ⟩/⟨S 2 ⟩. The same expression is used mutatis mutandis for the site-based avalanche size, S. More details about avalanche statistics can be found in Appendix C.

We determine x c and the exponents 1/σ, 1/σ, d f and df by measuring S c and Sc 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 Sc , as shown in Figs. 5 (a,b). We obtain x c = 0.560 ± 0.001. The values of the critical exponents are presented in Table II. Fig- ures 5 (c,d) display S c and Sc as a function of x c -x 0 , for various system size L. They indeed show that S c and Sc 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 ( 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 ∼ L d f and the previously obtained values of d f and df , see Figs 6 (c,d). From the data collapse we determine values of avalanche exponents τ and τ (see Table II). The successful collapse of the data confirms the validity of the scaling ansatz as well as the value of critical exponents. 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 [START_REF] James P Sethna | Crackling noise[END_REF][START_REF] Rosso | Avalanches and deformation in glasses and disordered systems[END_REF].

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) ∼ (x -x c ) θ , defining an exponent θ [START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF]. This behavior also occurs near the yielding transition of amorphous solids under shearing (with x c = 0) [START_REF] Lin | On the density of shear transformations in amorphous solids[END_REF], where it affects the scaling of flow properties [START_REF] Lin | Scaling description of the yielding transition in soft amorphous solids at zero temperature[END_REF] and plasticity [START_REF] Karmakar | Statistical physics of the yielding transition in amorphous solids[END_REF], and more generally in glassy systems with longrange interactions [START_REF] Müller | Marginal stability in structural, spin, and electron glasses[END_REF]. 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 θ = 1. Yet the kicks are much more broadly distributed, and in higher dimensions, this random process is akin to a Levy flight [START_REF] Lemaître | Plastic Response of a 2D Amorphous Solid to Quasi-Static Shear : II -Dynamical Noise and Avalanches in a Mean Field Model[END_REF], which can be shown to imply θ MF = 1/2 [START_REF] Lin | Mean-field description of plastic flow in amorphous solids[END_REF]. 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) ∼ (x -x c ) θ .
The following two features of P (x) connect the extremal and finite temperature dynamics. First, the point x c or the corresponding energy scale

E c = x 3/2 c
controls the effective energy barrier associated to τ α . Indeed, at small temperature the dynamics proceeds by relaxing the sites with smallest barriers, i.e., sites having barriers close to E c [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF][START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF]. This is in agreement with the relaxation time τ α observed in Fig. 1(b), which scales as τ α ∼ e Ec/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, ⟨x second -x min ⟩ ∼ N -δ . 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 ∼ x second -x min , we conclude

⟨E second -E min ⟩ ∼ N -δ . ( 4
)
Thus the exponent δ 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 δ is reported in Table II. Extreme value statistics argument would suggest δ = 1/(1 + θ) [START_REF] Lin | Scaling description of the yielding transition in soft amorphous solids at zero temperature[END_REF] (although near the yielding point deviations from this relation have been reported in finite dimension [START_REF] Ezequiel | Properties of the density of shear transformations in driven amorphous solids[END_REF], and explained in [START_REF] Korchinski | Signatures of the spatial extent of plastic events in the yielding transition in amorphous solids[END_REF]). Within the meanfield theory [START_REF] Lin | Mean-field description of plastic flow in amorphous solids[END_REF], one finds δ MF = 1/(1 + θ 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 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 S L -d f 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.

P( S)

L = 8 L = 16 L = 32 L = 64 L = 128 L = 256 (a) (b) (c) (d)

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 II. Critical exponents obtained from extremal dynamics simulations at T = 0 + for the scalar and tensorial elastoplastic 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.

Scalar model

Tensorial model δ = 0.64 ± 0.01 δ = 0.62 ± 0.01 τ = 1.25 ± 0.05 τ = 1.30 ± 0.05

d f = 2.3 ± 0.1 d f = 2.3 ± 0.1 1/σ = 2.2 ± 0.1 1/σ = 1.95 ± 0.03 τ = 1.25 ± 0.05 τ = 1.35 ± 0.03 df = 2.00 ± 0.02 df = 1.95 ± 0.05 1/σ = 1.9 ± 0.1 1/σ = 1.75 ± 0.03
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 ξ 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 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 min ≈ E c that occurs on a timescale τ α ∼ e Ec/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 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 powerlaw 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, ⟨E second -E min ⟩, otherwise sites with higher energy barriers (such as E second ) might relax and dynamics cease to be extremal. According to Eq. ( 4), ⟨E second -E min ⟩ depends on the system size N . Thus, the condition for the extremal dynamics which interplays T and N is given by T

≪ ⟨E second -E min ⟩ ∼ N -δ , i.e., N ≪ N T ∼ T -1/δ .
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 independent. 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. Moreover, it also directly predicts that the size of dynamically correlated regions characterized by χ * 4 , which is given by the cross-over size above which the condition for the extremal dynamics breaks down: χ * 4 ∼ ξ df leading to χ * 4 ∼ T -df /dδ . To conclude, we derive two scaling relations for thermal avalanches:

γ = df dδ , (5)
ν = 1 dδ . (6)
These scaling relations connect dynamical heterogeneities at finite T (characterized by γ and ν) and the distribution of local energy barriers at T = 0 + (characterized by δ) together with the morphology of avalanches during the extremal dynamics (by df ). We thus predict γ and ν by using δ and df measured in the extremal dynamics in Sec. IV. These predictions are in good agreement 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 df = d = 2 lead to similar values. This argument is expected to break down if large spatial correlations characterize the structure of the system, causing the locations where avalanches start to be correlated. Indeed, as is always the case in disordered materials at zero temperature, any critical threshold 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 spatial 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 ν = ν ′ , corresponding to ∆E c ∼ E second -E min . ν ′ can be related to previously introduced exponents, 1/σ and df , as follows. If a finite system displays a system spanning avalanche S ∼ Sc ∼ L df 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 ∼ Sc ∼ L df is associated with a characteristic change of the critical threshold ∆x

c = x c -x 0 ∼ L -σ df (when z ≈ 1 in f (z)), implying that ξ ∼ ∆x -1/(σ df ) c
. Since ∆x c ∼ ∆E c and ∆E c ∼ ξ -1/ν ′ , we obtain ν ′ = 1/(σ df ). Using numerical values for σ and df , this expression leads to ν ′ ≈ 0.95 and ν ′ ≈ 0.9, respectively, for the scalar and tensorial models. We thus have in the present system ν ≈ ν ′ , supporting 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 0avalanches 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 3/2 0 . Hence, it carries information both about the typical time-scale τ (x 0 ) ∼ e E(x0)/T and the avalanche (cut-off) size Sc (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, Sc , follows Sc ∼ (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(x0)/T . We thus obtain a relation connecting the duration of avalanches, τ (x 0 ) (measured in the unit of the relaxation time τ α ∼ e Ec/T ), with the cutoff size Sc as S-σ c ∼ (x c -x 0 ) ∼ E c -E(x 0 ) ∼ T ln(τ α /τ (x 0 )). Therefore, 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 Sc (t), which is given by

Sc (t) ∼ [T ln(τ α /t)] -1/σ . (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 [START_REF] Candelier | Spatiotemporal hierarchy of relaxation events, dynamical heterogeneities, and structural reorganization in a supercooled liquid[END_REF][START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF]. In Ref. [START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF], the time evolution of a chord length ⟨ℓ⟩ characterizing the linear size of mobile domains has been measured. Since we expect Sc ∼ ⟨ℓ⟩ df , our prediction for ⟨ℓ⟩ is given by

⟨ℓ⟩ Pred = A(T ) [T ln(τ α /t)] -1/(σ df ) , (8)
where A(T ) is a function with a finite limiting value as T → 0. In Fig. 8, we compare the molecular simulation FIG. 8. Comparison between the molecular simulation data (empty points) in Ref. [START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF] and our theoretical prediction (solid curves) in Eq. ( 8) for the time evolution of a linear avalanche size ⟨ℓ⟩. We set 1/σ = 1.8 and df = 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.

1 0 -5 1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 0 1 0 1 T =
data for a two-dimensional polydisperse mixture [START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF] 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 Sc (t) and χ 4 (t) explicitly. This allows us to predict the time evolution of χ 4 (t) for times 1 ≪ t ≪ τ α -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 [START_REF] Tarjus | Breakdown of the stokes-einstein relation in supercooled liquids[END_REF][START_REF] Mark | Spatially heterogeneous dynamics in supercooled liquids[END_REF][START_REF] Sengupta | Breakdown of the stokes-einstein relation in two, three, and four dimensions[END_REF][START_REF] Charbonneau | Hopping and the stokes-einstein relation breakdown in simple glass formers[END_REF][START_REF] Kawasaki | Identifying time scales for violation/preservation of stokes-einstein relation in supercooled water[END_REF].

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 [START_REF] Jung | Excitation lines and the breakdown of stokeseinstein relations in supercooled liquids[END_REF][START_REF] Berthier | Length scale for the onset of fickian diffusion in supercooled liquids[END_REF][START_REF] Lester O Hedges | Decoupling of exchange and persistence times in atomistic models of glass formers[END_REF] (or the caging time in molecular simulations [START_REF] Pastore | Breakdown of the stokes-einstein relation in supercooled liquids: A cage-jump perspective[END_REF][START_REF] Pica Ciamarra | Particle jumps in structural glasses[END_REF]), which has been argued to be the key ingredient of decoupling between diffusion and relaxation in super-cooled liquids [START_REF] Jung | Excitation lines and the breakdown of stokeseinstein relations in supercooled liquids[END_REF][START_REF] Berthier | Length scale for the onset of fickian diffusion in supercooled liquids[END_REF][START_REF] Lester O Hedges | Decoupling of exchange and persistence times in atomistic models of glass formers[END_REF][START_REF] Chaudhuri | Universal nature of particle displacements close to glass and jamming transitions[END_REF][START_REF] Pastore | Breakdown of the stokes-einstein relation in supercooled liquids: A cage-jump perspective[END_REF].

In our case, this effect originates from the difference between the event-based and site-based avalanche sizes characterized by d f > df . 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 τ α and liner size is ξ. During the formation, a single site relaxes order of ξ d f /ξ df 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, ξ (d f -df ) /τ α , leading to a Stoke-Einstein breakdown Dτ α of order Dτ α ∼ ξ (d f -df ) ∼ T -h that diverges at vanishing temperature with h = ν(d f -df ). Conventionally, the Stokes-Einstein violation is considered a consequence of spatially heterogeneous dynamics [START_REF] Mark | Spatially heterogeneous dynamics in supercooled liquids[END_REF]. We directly connect the former and the length scale of dynamical heterogeneity ξ. 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 ∼ τ -ζ α , or Dτ α ∼ τ 1-ζ α , but instead Dτ α ∼ (log τ α ) 1-ζ ′ where 1 -ζ ′ = ν(d f -df ).
The former is the fit that is usually conjectured from experimental data [START_REF] Stephen F Swallen | Self-diffusion of trisnaphthylbenzene near the glass transition tempera-ture[END_REF][START_REF] Mallamace | Transport properties of glassforming liquids suggest that dynamic crossover temperature is as important as the glass transition temperature[END_REF] with ζ ≈ 0.8. However, given the small value of 1 -ζ 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τ α 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 [START_REF] Jung | Excitation lines and the breakdown of stokeseinstein relations in supercooled liquids[END_REF][START_REF] Berthier | Length scale for the onset of fickian diffusion in supercooled liquids[END_REF] (see Appendix B for details). Dτ α increases with decreasing T , following the scaling prediction Dτ α ∼ T -h with h = ν(d f -df ) 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 [START_REF] Mark | Spatially heterogeneous dynamics in supercooled liquids[END_REF][START_REF] Ozawa | Tuning pairwise potential can control the fragility of glass-forming liquids: From a tetrahedral network to isotropic soft sphere models[END_REF]. 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 [START_REF] Dyre | Colloquium: The glass transition and elastic models of glass-forming liquids[END_REF][START_REF] Pica Ciamarra | The energy cost of local rearrangements, not cooperative effects, makes glasses solid[END_REF][START_REF] Lerbinger | Relevance of shear transformations in the relaxation of supercooled liquids[END_REF][START_REF] Sivert | Dielectric relaxation of isoamyl bromide[END_REF][START_REF] Muhammad | A theory of localized excitations in supercooled liquids[END_REF]. 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 powerlaw relations, rather than in terms of the relaxation time scale.

One important aspect of the models we studied is that 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 [START_REF] Ritort | Glassy dynamics of kinetically constrained models[END_REF].

Although local excitations identified in the dynamic facilitation scenario [START_REF] Keys | Excitations are localized and relaxation is hierarchical in glassforming liquids[END_REF][START_REF] Isobe | Applicability of dynamic facilitation theory to binary hard disk systems[END_REF] 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, 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 [START_REF] Yael S Elmatad | Finite-temperature critical point of a glass transition[END_REF][START_REF] Turci | Nonequilibrium phase transition in an atomistic glassformer: The connection to thermodynamics[END_REF]. 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 χ * 4 increases by about two decades as the glass transition is reached [START_REF] Dalle-Ferrier | Spatial correlations in the dynamics of glassforming liquids: Experimental determination of their temperature dependence[END_REF][START_REF] Dauchot | The glass transition in molecules, colloids and grains: universality and specificity[END_REF]), whereas the magnitude of the Stoke-Einstein breakdown is comparable to that of rather strong glass-forming liquids [START_REF] Mark | Spatially heterogeneous dynamics in supercooled liquids[END_REF]. According to our scaling argument, such a breakdown will increase if the system shows more intensive accumulations of multiple relaxations characterized by larger ν(d f -df ). 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 [START_REF] Agoritsas | On the relevance of disorder in athermal amorphous materials under shear[END_REF], 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 [START_REF] Tanaka | Critical-like behaviour of glass-forming liquids[END_REF][START_REF] Royall | The role of local structure in dynamical arrest[END_REF][START_REF] Tanaka | Revealing key structural features hidden in liquids and glasses[END_REF][START_REF] Paret | Assessing the structural heterogeneity of supercooled liquids through community inference[END_REF]. Such generalization would allow us to study the structure-dynamics relationship [START_REF] Widmer-Cooper | How reproducible are dynamic heterogeneities in a supercooled liquid?[END_REF][START_REF] Widmer-Cooper | Irreversible reorganization in a supercooled liquid originates from localized soft modes[END_REF][START_REF] Hocky | Correlation of local order with particle mobility in supercooled liquids is highly system dependent[END_REF] in elastoplastic models. Further improvements may be achieved by adding fluctuations and non-linearities to the propagator, which are present at short-range [START_REF] Lemaître | Anomalous elasticity and plastic screening in amorphous solids[END_REF][START_REF] Lerner | Breakdown of continuum elasticity in amorphous solids[END_REF]. An interesting line of research to develop quantitative models is a mapping from a molecular simulation to an elastoplastic model [START_REF] Liu | Elastoplastic approach based on microscopic insights for the steady state and transient dynamics of sheared disordered solids[END_REF][START_REF] Fernández Castellanos | Insights from the quantitative calibration of an elasto-plastic model from a lennard-jones atomic glass[END_REF][START_REF] David F Castellanos | History dependent plasticity of glass: A mapping between atomistic and elasto-plastic models[END_REF] or the one pursued in Refs. [START_REF] Tah | Fragility in glassy liquids: A structural approach based on machine learning[END_REF][START_REF] Xiao | Machine learning-informed structuro-elastoplasticity predicts ductility of disordered solids[END_REF], 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 [START_REF] Mark D Ediger | Supercooled liquids and glasses[END_REF]. 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 [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF][START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF]. 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 ∞ (T ) of the material, the activation energy E c will be proportional to G ∞ (T ), a correlation known to exist in some glass-forming liquids [START_REF] Dyre | Colloquium: The glass transition and elastic models of glass-forming liquids[END_REF][START_REF] Hecksher | A review of experiments testing the shoving model[END_REF]. More recent works relate this energy scale of local rearrangements to local (instead of global) elasticity [START_REF] Kapteijns | Does mesoscopic elasticity control viscous slowing down in glassforming liquids?[END_REF], plasticity [START_REF] Lerbinger | Relevance of shear transformations in the relaxation of supercooled liquids[END_REF], or alternatively to the varying geometry of elementary rearrangements under cooling [START_REF] Ji | Mean-field description for the architecture of low-energy excitations in glasses[END_REF]. Thus, it is worthwhile to reveal the relationship between the activation energy barriers in liquids and lowenergy excitations in glasses [START_REF] Scalliet | Nature of excitations and defects in structural glasses[END_REF][START_REF] Mizuno | Intermittent rearrangements accompanying thermal fluctuations distinguish glasses from crystals[END_REF][START_REF] Richard | Detecting low-energy quasilocalized excitations in computer glasses[END_REF]. Local energy barriers would also be related to locally-favored structures in microscopic perspectives [START_REF] Royall | The role of local structure in dynamical arrest[END_REF][START_REF] Coslovich | Understanding fragility in supercooled lennard-jones mixtures. i. locally preferred structures[END_REF]. 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 [START_REF] Kapteijns | Does mesoscopic elasticity control viscous slowing down in glassforming liquids?[END_REF][START_REF] Pica Ciamarra | The energy cost of local rearrangements, not cooperative effects, makes glasses solid[END_REF][START_REF] Lerbinger | Relevance of shear transformations in the relaxation of supercooled liquids[END_REF] 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 [START_REF] Scalliet | Thirty milliseconds in the life of a supercooled liquid[END_REF][START_REF] Biroli | The rfot theory of glasses: Recent progress and open issues[END_REF]. 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 [START_REF] Fernandez | Avalanche behavior in creep failure of disordered materials[END_REF][START_REF] Bauer | Collective rearrangement at the onset of flow of a polycrystalline hexagonal columnar phase[END_REF][START_REF] Caton | Plastic behavior of some yield stress fluids: from creep to long-time yield[END_REF][START_REF] Divoux | From stress-induced fluidization processes to herschel-bulkley behaviour in simple yield stress fluids[END_REF][START_REF] Siebenbürger | Creep in colloidal glasses[END_REF][START_REF] Grenard | Timescales in creep and yielding of attractive gels[END_REF][START_REF] Leocmach | Creep and fracture of a protein gel under stress[END_REF] or that of pinned elastic interfaces [START_REF] Víctor | Creep and thermal rounding close to the elastic depinning threshold[END_REF][START_REF] Bustingorry | Thermal rounding of the depinning transition[END_REF][START_REF] Kolton | Creep motion of an elastic string in a random potential[END_REF] below the threshold force where they spontaneously flow. 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.

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 σi . Note that symbols σ and σ 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: σxx,i and σ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 θ Y with respect to the σxx axis in shear stress space, centered at zero shear stress and separated by 2σ Y (see Fig 11). The local yield surface is assigned for each site i. We initiate θ Y,i with a uniformly distributed random number in [0, 2π).

When a site i becomes unstable it undergoes a plastic event over a timescale τ 0 : σxx,i → σxx,i -δσ xx,i and σxy,i → σxy,i -δσ xy,i , where the amount of stress drops δσ xx,i and δσ xy,i are given by

δσ xx,i = -(z -x) sin(θ Y,i ) sgn(σ xx,i ), (A2) δσ xy,i = -(z -x) cos(θ Y,i ) sgn(σ xy,i ), (A3)
respectively. sgn(x) is the sign function and z is a random number drawn from an exponential distribution p(z) = e -z/z0 /z 0 with z 0 = 1. The duration of a plastic event τ 0 is accounted for by triggering the relaxation with a probability per unit time 1/τ 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 /τ 0 , where E(x) = x 3/2 is the local energy barrier [START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF]. After each plastic event, we redraw the angle of the yield surface θ Y,i from a uniform random distribution. Note that to simulate such a dynamics at low temperatures we implement a Gillespie type of algorithm [START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF], which operates as follows. Consider an event occurring at some time t. Following it, a relaxation time τ i for each site i is chosen with an exponential distribution of mean e -E(xi)/T /τ 0 . The next event corresponds to the smallest τ i , leading to a plastic event on the corresponding site, which occurs at a time t + τ i . At that point, stresses are computed once again, the variables τ i are sampled from 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σY and each line makes an angle of θY with σxx axis. The distance from the yield surface is shown by x.

. 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,

G xx,xx (q) = - (q 2 x -q 2 y ) 2 (q 2 x + q 2 y ) 2 , (A4)
G xy,xy (q) = - 4q 2 x q 2 y (q 2 x + q 2 y ) 2 , (A5) G xx,xy (q) = G xy,xx (q) = - 2q x q y (q 2 x -q 2 y ) (q 2 x + q 2 y ) 2 , (A6)
where 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 2 x = 2 -2 cos (2πn x /L) , q 2 y = 2 -2 cos (2πn y /L) and q x q y = 2 sin (2πn x /L) sin (2πn y /L) where

n α = -L/2 + {1, • • • , L} with α = x, y.
text, respectively. We then measure a four-point correlations function, χ 4 (t), defined by

χ 4 (t) = L d ⟨π 2 (t)⟩ time -⟨π(t)⟩ 2 time . (B1)
χ 4 (t) for the scalar and tensorial models are presented in Fig. 3 of Ref. [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF] and Fig. 2(b) in the main text, respectively. χ 4 (t) quantifies the size of the dynamically correlated region, because one can rewrite it as

χ 4 (t) = 1 L d i,j ⟨p i (t)p j (t)⟩ time -⟨π(t)⟩ 2 time (B2) = 1 L d i k ⟨ϕ i (t)ϕ i+k (t)⟩ time , (B3)
where ϕ i (t) = p i (t) -⟨π(t)⟩ time . For example, ϕ i (τ α ) = ±1/2. Therefore, χ 4 (t) is proportional to the average number of sites correlated dynamically. Therefore, its peak value, χ * 4 , contains essentially the same information as S, in particular, χ * 4 ∼ Sc at T = 0 + . Figure 12(a) shows χ * 4 versus T for several system sizes L for the scalar model. One can see a scaling regime, χ * 4 ∼ T -γ at lower T and larger L. A scaling collapse is obtained in Fig. 12(b), which determines another critical exponent ν associated with a lengthscale of dynamical heterogeneity. The obtained values for the critical exponents are reported in Table II.

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) [START_REF] Lačević | Spatially heterogeneous dynamics investigated via a time-dependent four-point density correlation function[END_REF], defined by

S 4 (q, t) = 1 L d ij ⟨p i (t)p j (t)⟩ time -⟨π(t)⟩ 2 time e iq•(ri-rj ) , (B4)
where q = |q|. In Fig. 13(a), we show S 4 (q, t) for the scalar model at t = τ * when χ 4 (t) takes the peak value, χ * 4 = χ 4 (τ * ). We find that τ * is close to τ α , and thus S 4 (q, τ * ) encodes heterogeneity associated with structural relaxation. We note that at the long wave-length limit, S 4 (q, t) converges to χ 4 (t), namely, lim q→0 S 4 (q, t) = χ 4 (t). We then assume the Ornstein-Zernike form at lower q, S 4 (q, τ * ) =

χ * 4 1 + (qξ 4 ) a , ( B5
)
where ξ 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 ξ 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 ξ 4 ∼ T -ν with ν = 0.9, which is consistent with the one estimated from the finite size scaling in Fig. 12. Moreover, the χ * 4 versus ξ 4 plot [START_REF] Flenner | Universal features of dynamic heterogeneity in supercooled liquids[END_REF][START_REF] Kim | Multiple length and time scales of dynamic heterogeneities in model glass-forming liquids: A systematic analysis of multi-point and multitime correlations[END_REF] in Fig. 13(d) provides us with the fractal dimensions, df = 2, which is also consistent with the one measured in the extremal dynamics in Fig. 18.

Prediction for the four-point correlation function

We connect the cutoff size for the site-based avalanche size, Sc (t), in a given time interval t, and the time evolution of the four-point correlation function, χ 4 (t).

We first compute the site-based avalanche size by S(t) = N i=1 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 S(t) measured by the time average are given by respectively. Consider a correlation volume of linear extension ξ. On this length scale, χ 4 (t) crosses-over toward its value for an infinite system. It is also the largest length for which extremal dynamics applies, implying that S(t) is distributed in a power-law fashion. Thus, Sc (t) can be estimated by ⟨ S2 (t)⟩ time /⟨ S(t)⟩ time , as given by Sc (t) = i,j {⟨p i (t)p j (t)⟩ time -(2⟨π(t)⟩ time -1)} N (1 -⟨π(t)⟩ time ) .

⟨ S(t)⟩ time = N (1 -⟨π(t)⟩ time ), ⟨ S2 (t)⟩ time = i,j {⟨p i (t)p j (t)⟩ time -(2⟨π(t)⟩ time -1)} ,
(B6) In general, one can expect that ⟨π(t)⟩ time follows the (stretched) exponential decay, ⟨π(t)⟩ time ≃ e -(t/τα) β , where β is an exponent. Typically, 0 < β ≤ 1 for equilibrium supercooled liquids. Our elastoplastic models show nearly exponential relaxation with β ≃ 1. We now consider the early time stage, where ⟨π(t)⟩ time can be approximated by ⟨π(t)⟩ time ≃ 1 -(t/τ α ) β . Under such a circumstance, Eq. (B6) and χ 4 (t) defined in Eq. (B2) suggest that Sc (t) ≃ (t/τ α ) -β χ 4 (t). Together with Eq. ( 7) in the main text, we predict the time evolution of χ 4 (t) as

χ 4 (t) ≃ A(t/τ α ) β [T ln(τ α /t)] -1/σ , ( B7
)
where A is a constant which does not depend on t and T . Figure 14(a) shows a Log-Log plot for χ 4 (t) measured at finite temperature simulations for the scalar model. The initial growth can be fitted effectively by a powerlaw, χ 4 (t) ∼ t b [START_REF] Biroli | Dynamical heterogeneity in glass-forming liquids[END_REF][START_REF] Flenner | Dynamic heterogeneity in two-dimensional supercooled liquids: Comparison of bond-breaking and bond-orientational correlations[END_REF], with b ≃ 1.4. Instead, our argument in Eq. (B7) predicts a linear growth with a log- arithmic correction. In Fig. 14(b), we show a parametric plot to numerically test Eq. (B7) with β = 1, A = 0.5, and 1/σ = 1.9 (see Table II). We find that the simulated χ 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.

1 0 0 1 0 1 1 0 2 1 0 3 1 0 4 1 0 5 1 0 6 1 0 7 1 0 -5 1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 2 1 0 3 ( a ) 1 0 -5 1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 1 1 0 2 1 0 3 1 0 -5 1 0 -3 1 0 -1

Tracer particles

We monitor the diffusion of tracer particles [START_REF] Jung | Excitation lines and the breakdown of stokeseinstein relations in supercooled liquids[END_REF][START_REF] Berthier | Length scale for the onset of fickian diffusion in supercooled liquids[END_REF] 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 τ α , 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 ∆ 2 (t) of the tracers by

∆ 2 (t) = 1 N t Nt k=1 ⟨∆r 2 k (t)⟩ time , (B8)
where ∆r 2 and N t is the number of tracer particles. In Fig. 15(b) we show ∆ 2 (t) for different temperatures. We find a diffusive behavior at a larger time, ∆ 2 (t) = Dt, from which we extract the diffusion coefficient D for each temperature.

2 k (t) = (x k (t) -x k (0)) 2 + (y k (t) -y k (0))
Appendix C: Avalanches at T = 0 +

Extremal dynamics

We explain the extremal dynamics at T = 0 + . In the finite temperature simulations described in Sec. A, we 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 ≤ x ≤ 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 3/2 min , 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 [START_REF] Víctor | Creep and thermal rounding close to the elastic depinning threshold[END_REF][START_REF] Paczuski | Avalanche dynamics in evolution, growth, and depinning models[END_REF][START_REF] Baret | Extremal model for amorphous media plasticity[END_REF].

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. [START_REF] Paczuski | Avalanche dynamics in evolution, growth, and depinning models[END_REF] 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 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 , ..., 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, S, by the number of sites participating in a single avalanche. By construction, S ≤ S. The distinction between S and 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 S depend on the threshold value x 0 . As x 0 is increased, the size of avalanches, S and S, increases. As we will discuss further below, avalanches become a system spanning when x 0 → x c , where x c is the critical value associated with the critical energy gap E c = x 3/2 c . We will vary x 0 systematically to probe the critical behavior associated with x c (see below).

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 , ..., S M } and { S1 , S2 , ..., SM }, respectively, where M is the number of data points. In this paper, we analyze S and S in parallel. Below we will explain how to analyze the data using S, but the same procedures are applied for S. We first define the m-th moments of avalanche distribution where P (S) is the distribution of avalanches. One expects that P (S) follows a power law distribution, P (S) ∼ S -τ g(S/S c ), (C2)

where τ is a critical exponent, S c is a cutoff size, and g(z) is a scaling function. Figure 17 shows P (S) and P ( S) for several x 0 for the scalar model. These plots demonstrate that the size of avalanches (both S and 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 < τ < 2, one obtains

⟨S m ⟩ ∼ ∞ 0 dSS m-τ g(S/S c ) ∼ S m+1-τ c , (C3)
which implies S c ∼ ⟨S m+1 ⟩/⟨S m ⟩. Thus, in practice, we define S c by S c = ⟨S 3 ⟩/⟨S 2 ⟩. Following Ref. [START_REF] Han | Critical behavior of a generalized bak-sneppen model[END_REF], we assume

S c ∼ (x c -x 0 ) -1/σ f L d f (x c -x 0 ) -1/σ , ( C4
)
where f (z) = 1 for z ≫ 1 and f (z) = z for z → 0. Thus, S c ∼ L d f when x 0 → x c and S c ∼ (x c -x 0 ) -1/σ when L → ∞. Figures 18(a,b) show S c and Sc approaching x 0 → x c for the scalar model. Both S c and Sc increase with increasing x 0 with an eventual saturation due to a finite-size effect. One can see the expectated behavior, S c ∼ (x c -x 0 ) -1/σ at larger L ( Sc 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/σ, 1/σ, d f , and df , 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 ( S) at x 0 = x c , as shown in Figs. 19(a,b) for the scalar model. They show characteristic power-law behavior with cut-off S c and Sc due to a finite size effect, which scales as S c ∼ L d f and Sc ∼ L df , respectively. The power-law behavior determines τ and τ 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 configuration right before (or after) each avalanche defined by x 0 = x c starts (or ends). Thus we exclude configurations during each avalanche and focus only on stable configurations expected to hold at strictly T = 0, where dynamics is not allowed. In Fig. 20(a), we show the measured P (x) for the scalar model together with P (x) obtained from the finite temperature simulations. As T is lowered, P (x) for finite T converges to P (x) obtained from the extremal dynamics, whose functional form is consistent with P (x) ∼ (x -x c ) θ , expected from other disordered systems. We then measure ⟨E second -E min ⟩ ∼ N -δ in Fig. 20(b), which encodes an important feature of the distribution of activation energy barriers at T = 0 (see Sec. IV for detail discussions). The obtained exponent δ is reported in Table II.

FIG. 1 .FIG. 2 .

 12 FIG. 1. Dynamics of the tensorial model in finite temperature simulations. (a): Mean persistence correlation function, ⟨π(t)⟩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 τα versus 1/T . The red dashed-line corresponds to τα ∼ e Ec/T with Ec = x 3/2 c = 0.42 measured independently in Sec. IV.

FIG. 3 .

 3 FIG. 3. (a): The peak of four-point correlation function, χ * 4 , versus T for several L for the tensorial model. The red dashed-line follows χ * 4 ∼ T -γ . (b): The corresponding scaling collapse.

FIG. 4 .

 4 FIG. 4. Snapshots of an avalanche formation for the extremal dynamics of the tensorial model with L = 256. Event-based avalanche size is S ≃ 2.7 × 10 5 , while the site-based avalanche size is S ≃ 3 × 10 4 . Purple shows immobile sites (zero event), and the colorbar shows the number of relaxation events in mobile sites.

FIG. 5 .

 5 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 Sc = ⟨S 3 ⟩/⟨S 2 ⟩ based on Eq. (3), for various L for the event-based (a) and site-based (b) avalanche sizes, which determines the critical threshold xc and critical exponents, 1/σ, 1/σ, d f , and df . (c, d): Sc versus xc-x0. The red dashed line corresponds to Sc ∼ (xc-x0) -1/σ in (c) and Sc ∼ (xc -x0) -1/σ in (d).

FIG. 6 .

 6 FIG. 6. (a, b): Distribution of avalanche size P (S) (a) and P ( S) (b) for a stability threshold x0 = xc for the tensorial model, with varying the system size L. The dashed lines follow P (S) ∼ S -τ (a) and P ( S) ∼ S-τ (b). (c, d): The corresponding scaling collapse following Eq. (2) with Sc ∼ L d f (c). The same expressions are used for P ( S) in (d).

FIG. 7 .

 7 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 xc. (b): Average of E second -Emin obtained from the extremal dynamics for various N = L d . The red dashed line corresponds to ⟨E second -Emin⟩ ∼ N -δ .

FIG. 9 .

 9 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 Ec. On a time scale τα ∼ e Ec/T , a site with a barrier near Ec relaxes (red arrow), which we call a "nucleation event". As a result, due to elastic interactions, other sites may display lower barriers E < Ec (blue arrows), which we call "facilitated sites". They relax on a time scale much faster than τα, 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 ξ corresponds to the maximum size for which extremal dynamics applies, which defines the length of dynamical heterogeneity at finite temperature.

FIG. 10 .

 10 FIG. 10. The Stoke-Einstein decoupling Dτα, where D is the diffusion coefficient of tracers, for the tensorial elasto-plastic model with L = 128. (a): Dτα versus τα could suggest an effective power-law behaviour (the fractional Stokes-Einstein violation), with two apparent exponents in the range τα < 10 4 and τα > 10 4 . (b): Our prediction suggests instead a powerlaw behaviour in terms of T . The red dashed line represents Dτα ∼ T -ν(d f -df ) .

FIG. 12 .

 12 FIG. 12. The four-point correlation function χ4 for the scalar model. (a): The peak value χ * 4 versus T for several L in a Log-Log plot. The dashed line follows χ * 4 ∼ T -γ . (d): Scaling collapse of χ * 4 (L, T ), which determines ν.

FIG. 13 .

 13 FIG. 13. The four-point structure factor S4(q, t) and the associated correlation length ξ4 for the scalar model. (a): S4(q, τ * ) 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 ξ4 versus T . The dashed line follows ξ4 ∼ T -ν . (b): χ * 4 versus ξ4. The dashed line follows χ * 4 ∼ ξ df 4 .

FIG. 14 .

 14 FIG. 14. (a): Log-Log plot for χ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.

FIG. 15 .

 15 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 (τα = 4.14 × 10 12 ). (b): Meansquared displacement ∆ 2 (t) of tracer particles are shown with points. The solid lines follow ∆ 2 (t) = Dt, from which we extract the diffusion coefficient D.

3 S i m u l a t i o n s t e p , FIG. 16 .

 3,16 FIG.[START_REF] Hedges | Dynamic order-disorder in atomistic models of structural glass formers[END_REF]. An example of the evolution of xmin 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, S1, S2, ..., are presented based on the threshold value x0 = 0.2. xc = 0.281 is indicated by the horizontal dashed line.

FIG. 17 .

 17 FIG. 17. Avalanche distributions P (S) (a) and P ( S) (b) for several x0 approaching the critical point xc = x0 = 0.281 for the scalar model with L = 256. The dashed lines in (a) and (b) follow P (S) ∼ S -τ and P ( S) ∼ S -τ , respectively.

FIG. 18 .FIG. 19 .

 1819 FIG. 18. Cut-off in avalanche distributions for the scalar model. (a, b): Sc = ⟨S 3 ⟩/⟨S 2 ⟩ (a) and Sc = ⟨ S3 ⟩/⟨ S2 ⟩ (b) versus xc -x0 for various system sizes. The dashed lines in (a) and (b) follow Sc ∼ (xc -x0) -1/σ and Sc ∼ (xc -x0) -1/σ , respectively. (c,d): Scaling collapse assuming Eq. (C4). The critical point xc, exponents σ, σ, d f , and df are determined by these scaling plots.

FIG. 20 .

 20 FIG. 20. (a): Distribution P (x) for the scalar model with L = 256, obtained from the finite temperature simulations and the extremal dynamics. The vertical arrow locates the critical threshold xc. (b): ⟨E second -Emin⟩, where the average is taken over configurations right before (or after) each avalanche formation, i.e., for which xmin > x0 = xc. The dashed line follows ⟨E second -Emin⟩ ∼ N -δ .

TABLE I .

 I Critical exponents (γ, ν) obtained from finite T simulations in the scalar and tensorial elastoplastic models in two-dimensions, compared with their predicted values proposed in Section V.

		Scalar model Tensorial model Pred. d = 2 Pred. MF
	γ	1.8 ± 0.1	1.7 ± 0.1	1.60 ± 0.05	1.5
	ν 0.85 ± 0.05	0.80 ± 0.05	0.80 ± 0.03	0.75
		IV. CRITICAL POINT AND EXTREMAL
			DYNAMICS AT T = 0 +

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

Scalar model

We study a scalar elastoplastic model [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF] in a twodimensional lattice whose linear box length is L using the lattice constant as the unit of length. For each site, we assign local shear stress σ i (scalar variable) at a position r i .

The dynamical rule for the simulation model is akin to Monte-Carlo dynamics [START_REF] Berthier | The monte carlo dynamics of a binary lennard-jones glass-forming mixture[END_REF]. We pick a site, say i, up randomly among L 2 sites. If σ i is greater (or lower) than or equal to a threshold σ Y > 0 (or -σ Y < 0), namely, |σ i | ≥ σ Y , this site shows a plastic event: σ i → σ i -δσ i , where δσ i is the local stress drop. We use an uniform threshold, σ Y = 1. Instead, if |σ i | < σ Y , with probability e -E(σi)/T , where E(σ i ) is a local energy barrier and T is the temperature, this site shows a plastic event: σ i → σ i -δσ i . This corresponds to a plastic rearragement induced by a local thermal activation. We employ [START_REF] Craig | Energy barrier scalings in driven systems[END_REF]. By introducing the local stress distance to threshold, x i = σ Y -|σ i |, we can rewrite E(x) = x 3/2 . This specific form of the local energy barrier is suggested by molecular simulation studies [START_REF] Lerbinger | Relevance of shear transformations in the relaxation of supercooled liquids[END_REF][START_REF] Yue Fan | How thermally activated deformation starts in metallic glass[END_REF] and previous elastoplastic models under shear [START_REF] Popović | Thermally activated flow in models of amorphous solids[END_REF][START_REF] Ferrero | Yielding of amorphous solids at finite temperatures[END_REF]. The stress drop δσ i associated with a plastic event is a stochastic variable. In this paper, we use

, where sgn(x) is the sign function and z > 0 is a random number drawn by an exponential distribution, p(z) = 1 z0 e -z/z0 . z 0 is the mean value and we set z 0 = 1. This exponential distribution would be realistic according to molecular simulations in Ref. [START_REF] Barbot | Local yield stress statistics in model amorphous solids[END_REF].

A local plastic event at site i influences all other sites (∀j ̸ = i) as

where r ji = r j -r i and ψ i ∈ [0, π/2) is a random orientation of the Eshelby kernel G ψ r . Numerical implementation of G ψ r is described in Ref. [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF]. 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 σ i (∀i) randomly while keeping the force balance, i.e., the sum of stresses along each row and column of lattice sites is strictly zero [START_REF] Popović | Elastoplastic description of sudden failure in athermal amorphous materials during quasistatic loading[END_REF][START_REF] Pollard | Yielding, shear banding, and brittle failure of amorphous materials[END_REF]. To study dynamical properties We explain how to analyze correlation functions for finite temperature simulations. We first consider the persistence two-point time correlation function, ⟨π(t)⟩ time , which is defined by π(t) = 1 L d 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. ⟨• • • ⟩ time denotes the time average at the stationary state. ⟨π(t)⟩ time for the scalar and tensorial models are presented in Fig. 1 of Ref. [START_REF] Ozawa | Elasticity, facilitation, and dynamic heterogeneity in glass-forming liquids[END_REF] and Fig. 1(a) in the main