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

As a liquid is cooled, the relaxation time $\tau_{\alpha}$ below which it acts as a solid – before displaying flow – grows from picoseconds at high temperatures up to minutes at the glass transition temperature $T_{g}$ Anderson (1995); Ediger et al. (1996); Debenedetti and Stillinger (2001); Berthier and Biroli (2011). In this regime the effective activation energy associated to $\tau_{\alpha}$ grows for many liquids, leading to a super-Arrhenius behavior. Approaching $T_{g}$, dynamics also becomes heterogeneous on a growing correlation length scale $\xi$  Kob et al. (1997); Yamamoto and Onuki (1998); Dalle-Ferrier et al. (2007); Karmakar et al. (2014). 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 Kirkpatrick et al. (1989); Lubchenko and Wolynes (2007); Biroli and Bouchaud (2012) 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 Ritort and Sollich (2003); Berthier et al. (2011). This lead to a theory Garrahan and Chandler (2002); Garrahan et al. (2007); Hedges et al. (2009) in which thermodynamics plays almost no role, but dynamics is heterogeneous (due to kinetic constraints) and the super-Arrhenius behavior is due to non-local rearrangements taking place over $\xi$. Free volume Turnbull and Cohen (1961) or elastic Dyre (2006); Rainone et al. (2020); Kapteijns et al. (2021) 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 Ciamarra et al. (2023). 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 Dyre (2006); Schrøder and Dyre (2020): a rearrangement corresponds to a local (plastic) strain, that affects stress away from it. Lemaitre Lemaître (2014) 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 Chowdhury et al. (2016); Tong et al. (2020); Wu et al. (2015); Maier et al. (2017); Steffen et al. (2022); Flenner and Szamel (2015); Klochko et al. (2022). Recently, it was shown that indeed elastic interactions play an important role in dynamical facilitation Chacko et al. (2021). Moreover, Ref. Lerbinger et al. (2022) showed that heterogeneous relaxations take place in the form of plastic rearrangements that are called shear transformations Falk and Langer (1998). Following these molecular simulations studies, Ref. Ozawa and Biroli (2023) demonstrated with elasto-plastic models (traditionally used to study the plasticity of amorphous solids under loading Picard et al. (2004); Nicolas et al. (2018); Vandembroucq et al. (2004); Lin et al. (2014a); Rossi et al. (2022)) 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 Berthier et al. (2005, 2011). Furthermore, such models also capture the emergence of a gap in the distribution of local barriers under cooling, which controls the dynamics Popović et al. (2021); Ozawa and Biroli (2023). 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 Ozawa and Biroli (2023) 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 Candelier et al. (2010); Keys et al. (2011); Yanagishima et al. (2017). Our analysis builds a link with systems that crackle such as disordered magnets, granular materials, and earthquakes Sethna et al. (2001); Rosso et al. (2022) even at finite temperatures Popović et al. (2021); Purrello et al. (2017); Yao and Jack (2023); Korchinski and Rottler (2022), revealing that dynamical heterogeneities are controlled by a critical point at zero temperature, with a diverging length scale $\xi\sim T^{-\nu}$ and correlation volume $\chi_{4}^{*}\sim T^{-\gamma}$. We provide a scaling argument expressing $\nu$ and $\gamma$ in terms of the distribution of energy barriers at $T=0$. Our analysis makes predictions on the power-law distribution of the size of dynamical heterogeneities, which could be tested in numerical simulations of super-cooled liquids thanks to the recent advances in the characterization of dynamical correlations Scalliet et al. (2021, 2022). Based on the properties of the zero-temperature fixed point governing thermal avalanches, we also show that the decoupling $D\tau_{\alpha}$ between particle diffusion $D$ and relaxation time $\tau_{\alpha}$ (the so-called Stoke-Einstein violation— Tarjus and Kivelson (1995); Ediger (2000); Sengupta et al. (2013); Charbonneau et al. (2014); Kawasaki and Kim (2017)) diverges as $D\tau_{\alpha}\sim T^{-h}$, where $h$ can be expressed in terms of avalanche properties at vanishing temperature. The key physical mechanism behind the Stoke-Einstein violation is the intensive accumulation of rearrangements in the mobile region that quantifies a larger diffusion constant $D$, relative to the structural relaxation time $\tau_{\alpha}$ Jung et al. (2004); Berthier et al. (2004); Hedges et al. (2007); Chaudhuri et al. (2007); Pastore et al. (2021). 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.

To study dynamical heterogeneities in low-temperature glass-forming liquids, we employ a generalization of elasto-plastic models Picard et al. (2004); Nicolas et al. (2018) 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 Picard et al. (2004); Lin et al. (2014a); Nicolas et al. (2018); Rossi et al. (2022) (some exceptions e.g., Ref. Bulatov and Argon (1994)) and they were generalized to take into account thermal fluctuations Ferrero et al. (2014) and then to study glass-forming liquids Ozawa and Biroli (2023). 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 Dyre (2006); Lemaître (2014), 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 Jagla (2020) and a scalar elasto-plastic model in which the shear stress is represented by a scalar variable Ozawa and Biroli (2023) (see Appendix A for details). In both models, thermal fluctuations are implemented through a probability rate $\tau_{0}^{-1}e^{-E/T}$ to trigger a plastic (mobile) event in an otherwise stable elastic (immobile) site, where $E$ is the local activation energy barrier and $\tau_{0}$ is a microscopic relaxation timescale Ferrero et al. (2014); Popović et al. (2021); Ferrero et al. (2021); Rodriguez-Lopez et al. (2023). The energy barrier $E$ can be related to the minimal amount of additional shear stress $x_{i}$ required to destabilize a site $i$. In particular, we consider $E(x)=cx^{\alpha}$, which is suggested by recent elasto-plastic models and molecular simulations Ferrero et al. (2021); Popović et al. (2021); Fan et al. (2014); Lerbinger et al. (2022); Rodriguez-Lopez et al. (2023). In this study, we set $\tau_{0}=1$, $c=1$, and the value $\alpha=3/2$ corresponds to a smooth microscopic potential Aguirre and Jagla (2018). As stress is relaxed locally at the site of the plastic event, stress in the system changes according to the Eshelby kernel Picard et al. (2004). 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 Maloney and Lemaitre (2006). 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 ¹¹1Full isotropy is slightly broken by the choice of bi-periodic boundary conditions.. 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.

We first perform finite temperature simulations and characterize dynamical properties. To this end, we consider the persistence two-point time correlation function, $\langle\pi(t)\rangle_{\rm time}$, which has been used widely in the context of kinetically constrained models Ritort and Sollich (2003). The observable $\pi(t)$ is defined by $\pi(t)=\sum_{i}p_{i}(t)/N$, where $p_{i}(t)=1$ if the site $i$ did not exhibit a plastic event until time $t$ and remained immobile, and $p_{i}(t)=0$ for mobile sites that relaxed at least once. The notation $\langle\cdots\rangle_{\rm time}$ denotes the time average at the stationary state at temperature $T$ reached after a long enough equilibration time. We have verified that the model does show dynamical slowing down by measuring $\langle\pi(t)\rangle_{\rm time}$ at different temperatures, as shown in Fig. 1(a). In Fig. 1(b), we find that the associated relaxation time $\tau_{\alpha}$ defined by $\langle\pi(\tau_{\alpha})\rangle_{\rm time}=1/2$ increases in a Arrhenius way. Figure 2(a) represents a snapshot of local persistence for $\tau/\tau_{\alpha}\simeq 0.53$, demonstrating spatially heterogeneous dynamics. To quantify the magnitude of dynamical heterogeneity, we compute the four-point correlations function Donati et al. (2002), $\chi_{4}(t)$, defined by

$\chi_{4}(t)$ is proportional to the size of the dynamically correlated region (see Appendix B for details), which is the central observable characterizing dynamical heterogeneity approaching the glass transition Berthier et al. (2011), as it has been estimated in real experiments Berthier et al. (2005); Dalle-Ferrier et al. (2007); Capaccioli et al. (2008); Dauchot et al. (2022) as well as molecular simulations Donati et al. (2002); Karmakar et al. (2009); Coslovich et al. (2018).

Figure 2(b) shows the time and temperature evolution of $\chi_{4}(t)$. It takes a peak near the relaxation timescale $\tau_{\alpha}$, and the peak grows with decreasing temperature, which is the hallmark of dynamical heterogeneity in glassy dynamics. We have then plotted the peak value of $\chi_{4}(t)$, denoted $\chi_{4}^{*}$, versus temperature $T$ for several system sizes in Fig. 3(a). The system size dependence, akin to molecular simulations Karmakar et al. (2009); Chakrabarty et al. (2017), allowed us to perform finite size scaling. We find that for increasing the system size $L$, $\chi_{4}^{*}$ follows a scaling form, $\chi_{4}^{*}\sim T^{-\gamma}$ with a critical point at $T=0$ and the associated exponent $\gamma$. We obtain a scaling collapse for $\chi_{4}^{*}(L,T)$ in Fig 3(b) (see also molecular simulation studies Karmakar et al. (2009); Chakrabarty et al. (2017)), indicating that the dynamics is governed by a diverging correlation length toward $T=0$, i.e., $\xi\sim T^{-\nu}$ with an exponent $\nu$. The corresponding data for the scalar model are presented in Appendix B. Besides, we obtained a consistent result in terms of $\nu$ by directly measuring the dynamical correlation lengthscale based on the four-point structure factor Lačević et al. (2003). The observed critical exponents and predictions (see below) are summarized in Table 1. 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_{\rm min}=x_{\rm min}^{3/2}$, the weakest site, always yields first ²²2We are considering either finite systems or the thermodynamic limit taken after the zero-temperature limit., where $x_{\rm min}$ is the corresponding stress required to destabilise the site. Therefore, one can simulate dynamics at $T=0^{+}$ by relaxing always the weakest site, instead of relaxing a random site weighted by the relaxation rate $e^{-E(x)/T}$. This algorithm allows us to access dynamical information even at vanishing temperature (see details in Appendix C). It is an example of extremal dynamics, which is well-studied in the context of self-organized-criticality Paczuski et al. (1996) and some disordered materials under quasi-static driving Baret et al. (2002); Purrello et al. (2017), including annealed glasses Kumar et al. (2022).

In the thermodynamic limit, $L\to\infty$, the extremal dynamics leads to an absorbing condition for $E\leq E_{c}$, where $E_{c}$ is a critical energy barrier as found (for $T\rightarrow 0$) in Ref. Ozawa and Biroli (2023). As a result, the distribution of local (activation) energy barriers $P(E)$ vanishes for $E\leq E_{c}$ (see the sketch in Fig. 9(a)). This implies that $\lim_{L\rightarrow\infty}P(x)=0$ for $x\leq x_{c}=E_{c}^{2/3}$, where $P(x)$ is the distribution of $x$. Thus only a sub-extensive number of sites are found for $x\leq 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 Paczuski et al. (1996). 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 Paczuski et al. (1996); Purrello et al. (2017). For a finite $L$, we fix a chosen threshold stability value $x_{0}\ (\leq x_{c})$, and define $x_{0}$-avalanches as the sequences of events for which $x_{\rm min}<x_{0}$ (see Fig. 16 in Appendix C). Two useful characterizations of avalanche size are the total number of relaxation events $S$ in a given sequence (event-based avalanche size), and the total number of sites $\tilde{S}$ that relaxed at least once during an avalanche (site-based avalanche size). By construction, $\tilde{S}\leq S$ and $\tilde{S}\leq L^{d}$. We show snapshots having $S$ and the corresponding $\tilde{S}$ in Fig. 4. The two can be (and, as we shall show, are) different as a site can relax multiple times within the same avalanche. Thus, the event-based avalanche size $S$ will allows us to quantify the accumulation of relaxation events in the mobile region (as emphasized in a recent molecular simulation study Scalliet et al. (2022)), whereas the site-based avalanche size $\tilde{S}$ will be associated to the spatial extent of dynamically correlated regions (as it contains the essentially same information as the persistence map in Fig. 2(a) and hence $\chi_{4}$). Sizes of avalanches, $S$ and $\tilde{S}$, depend on the threshold $x_{0}$. They grow with increasing $x_{0}$ and diverge at $x_{c}$.

By systematically exploring different values of $x_{0}$ one can determine the critical point $x_{c}$. In fact, one expects that the avalanche distribution $P(S)$ during the extremal dynamics follows a power-law with a scaling form Paczuski et al. (1996); Purrello et al. (2017); Han et al. (2018):

where $g(z)$ is a scaling function and $S_{c}$ is a cutoff size which takes the form:

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

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

Once $x_{c}$ is determined, we can study the statistics of the system-spanning avalanches relevant to the thermodynamic limit. Thus, we fix $x_{0}=x_{c}$ and measure the distribution $P(S)$ and $P(\tilde{S})$ of the avalanche sizes, see Figs 6 (a,b). We find that avalanche sizes are power-law distributed with eventual cut-off, consistent with the general assumption in Eq. (2). The scaling form in Eq. (2) collapses the data using $S_{c}\sim L^{d_{f}}$ and the previously obtained values of $d_{f}$ and $\tilde{d}_{f}$, see Figs 6 (c,d). From the data collapse we determine values of avalanche exponents $\tau$ and $\tilde{\tau}$ (see Table 2). 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 Sethna et al. (2001); Rosso et al. (2022).

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

The following two features of $P(x)$ connect the extremal and finite temperature dynamics. First, the point $x_{c}$ or the corresponding energy scale $E_{c}=x_{c}^{3/2}$ controls the effective energy barrier associated to $\tau_{\alpha}$. Indeed, at small temperature the dynamics proceeds by relaxing the sites with smallest barriers, i.e., sites having barriers close to $E_{c}$  Ozawa and Biroli (2023); Popović et al. (2021). This is in agreement with the relaxation time $\tau_{\alpha}$ observed in Fig. 1(b), which scales as $\tau_{\alpha}\sim e^{E_{c}/T}$.

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

Thus the exponent $\delta$ characterizes an important feature of the energy barrier relevant for low temperature dynamics. We numerically confirmed this scaling in Fig. 7(b), and the obtained value of $\delta$ is reported in Table 2. Extreme value statistics argument would suggest $\delta=1/(1+\theta)$ Lin et al. (2014a) (although near the yielding point deviations from this relation have been reported in finite dimension Ferrero and Jagla (2021), and explained in Korchinski et al. (2021)). Within the mean-field theory Lin and Wyart (2016), one finds $\delta_{\rm MF}=1/(1+\theta_{\rm 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 zero temperature critical point by obtaining all relevant exponents, summarized in Table 2. 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.

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

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 Glarum (1960); Dyre (2006); Lerbinger et al. (2022); Hasyim and Mandadapu (2021); Ciamarra et al. (2023). Our main result is the theoretical explanation of dynamical correlations in terms of a zero-temperature critical point with the associated scaling relations. This leads to quantitative predictions on the power-law statistics of thermal avalanches testable in more realistic systems. Our study suggests that dynamical heterogeneities in super-cooled liquids should be investigated in terms of the temperature $T$ to seek power-law relations, rather than in terms of the relaxation time scale.

One important aspect of the models we studied is that 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 Ritort and Sollich (2003). Although local excitations identified in the dynamic facilitation scenario Isobe et al. (2016); Keys et al. (2011) 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 Elmatad et al. (2010); Turci et al. (2017). 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 super-cooled liquids (as the estimated $\chi_{4}^{*}$ increases by about two decades as the glass transition is reached Dalle-Ferrier et al. (2007); Dauchot et al. (2022)), whereas the magnitude of the Stoke-Einstein breakdown is comparable to that of rather strong glass-forming liquids Ediger (2000). According to our scaling argument, such a breakdown will increase if the system shows more intensive accumulations of multiple relaxations characterized by larger $\nu(d_{f}-\tilde{d}_{f})$. This effect could be achieved by imposing that some of the model parameters (such as the local values of the yield stress, or how sites are coupled to the elastic field) are randomly distributed Agoritsas et al. (2015), 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 Tanaka et al. (2010); Royall and Williams (2015); Tanaka et al. (2019); Paret et al. (2020). Such generalization would allow us to study the structure-dynamics relationship Widmer-Cooper et al. (2004, 2008); Hocky et al. (2014) in elastoplastic models. Further improvements may be achieved by adding fluctuations and non-linearities to the propagator, which are present at short-range Lemaître et al. (2021); Lerner et al. (2014). An interesting line of research to develop quantitative models is a mapping from a molecular simulation to an elastoplastic model Liu et al. (2021); Castellanos et al. (2021, 2022) or the one pursued in Refs. Tah et al. (2022); Xiao et al. (2023), 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 Ediger et al. (1996). 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 Popović et al. (2021); Ozawa and Biroli (2023). This point results from the simplifying assumption that the energy scale of local rearrangements does not vary with temperature. In an improved (still simplified) model where all elastic energies follow the high-frequency elastic modulus $G_{\infty}(T)$ of the material, the activation energy $E_{c}$ will be proportional to $G_{\infty}(T)$, a correlation known to exist in some glass-forming liquids Dyre (2006); Hecksher and Dyre (2015). More recent works relate this energy scale of local rearrangements to local (instead of global) elasticity Kapteijns et al. (2021), plasticity Lerbinger et al. (2022), or alternatively to the varying geometry of elementary rearrangements under cooling Ji et al. (2022). Thus, it is worthwhile to reveal the relationship between the activation energy barriers in liquids and low-energy excitations in glasses Scalliet et al. (2019); Mizuno et al. (2020); Richard et al. (2023). Local energy barriers would also be related to locally-favored structures in microscopic perspectives Coslovich and Pastore (2007); Royall and Williams (2015). 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 Kapteijns et al. (2021); Lerbinger et al. (2022); Ciamarra et al. (2023) 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 Scalliet et al. (2022); Biroli and Bouchaud (2022). 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_{\rm 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 Castellanos and Zaiser (2018); Bauer et al. (2006); Caton and Baravian (2008); Divoux et al. (2011); Siebenbürger et al. (2012); Grenard et al. (2014); Leocmach et al. (2014) or that of pinned elastic interfaces Bustingorry et al. (2007); Purrello et al. (2017); Kolton et al. (2005) below the threshold force where they spontaneously flow.

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. Scalliet et al. (2022). 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).