Elastoplasticity Mediates Dynamical Heterogeneity Below the Mode-Coupling Temperature

As liquids approach the glass transition temperature, dynamical heterogeneity emerges as a crucial universal feature of their behavior. Dynamic facilitation, where local motion triggers further motion nearby, plays a major role in this phenomenon. Here we show that long-range, elastically-mediated facilitation appears below the mode-coupling temperature, adding to the short-range component present at all temperatures. Our results suggest deep connections between the supercooled liquid and glass states, and pave the way for a deeper understanding of dynamical heterogeneity in glassy systems.

As liquids approach the glass transition temperature, dynamical heterogeneity emerges as a crucial universal feature of their behavior. Dynamic facilitation, where local motion triggers further motion nearby, plays a major role in this phenomenon. Here we show that long-range, elastically-mediated facilitation appears below the mode-coupling temperature, adding to the short-range component present at all temperatures. Our results suggest deep connections between the supercooled liquid and glass states, and pave the way for a deeper understanding of dynamical heterogeneity in glassy systems.
Dynamic facilitation is the process whereby relaxation of a local region in a glassy system enables another local region to subsequently move and relax. Facilitation plays a major role in the spatiotemporal pattern of correlated relaxation events characteristic of supercooled liquids, referred to as "dynamical heterogeneity," which is perhaps the most striking hallmark of glassy dynamics [1]. The underlying structural and microscopic origins of dynamic facilitation (also called kinetic facilitation) and dynamical heterogeneity are still widely debated. At high temperatures facilitation appears to be a local, shortranged process [2][3][4][5][6]. However it has been appreciated for decades that some relaxation processes in glassy media involve long-ranged correlations. Given that dynamic facilitation in supercooled liquids arises initially at short time scales, where the system should behave effectively as a solid, it is natural to wonder whether long-ranged elastic processes may influence facilitation and dynamical heterogeneity.
An extreme example of the influence of elasticity is the low-temperature relaxation of tunneling defects in glasses, which is mediated by long-ranged (dipolar) phonon exchange [7][8][9]. Another example is provided by the plastic behavior of amorphous solids. Local plastic deformations in elastic media induce signature quadrupolar perturbations to the stress and strain fields in the surrounding material which decay as power laws in space [10]. The role of such fields in triggering relaxation elsewhere when an amorphous solid is mechanically deformed has been the focus of intense research activity during the last decade [11][12][13][14][15][16][17]. This triggering by stress is a type of "elastically-mediated" dynamic facilitation that has been thoroughly studied in the context of the rheology of amorphous solids by means of elastoplas-tic models [16], in which each rearrangement perturbs the stress of the surrounding material, thereby triggering new rearrangements when the perturbing stress surpasses some local threshold value. It should be noted, however, that elastoplastic models do not incorporate the type of thermal motion that occurs in equilibrium liquids. Recent reports of anisotropic spatial correlations in stress or strain in liquids [18][19][20][21][22][23] have suggested that elasticity might play a role even at high temperatures, but it is now understood that such stress correlations must arise for isotropic systems in mechanical equilibrium and do not require elasticity [24,25].
Here we study by numerical simulation the effects of elasticity on dynamic facilitation in supercooled liquids at temperatures above and below the mode-coupling temperature, T MCT . We examine both the strain response to a rearranging particle as well as the pair correlation function of rearrangements, g hop (r), characterizing the probability of finding a rearrangement at a relative position r within a short time interval following a rearrangement at the origin. At all temperatures studied we find that the strain response of the neighborhood of a rearranging particle is anisotropic and qualitatively consistent with the response of a linearly elastic medium. The pair correlation function of rearrangements, however, reveals a change of behavior with temperature. At high temperatures, dynamic facilitation is local, i.e. it acts only on nearby regions, as indeed assumed in simple models [3,26]. However, dynamic facilitation becomes progressively longer-ranged and elastically-mediated with decreasing temperature. Remarkably, the emergence of long-ranged facilitation takes place near T MCT , strengthening the interpretation of the mode-coupling crossover as the temperature at which the system starts to display arXiv:2103.01852v2 [cond-mat.soft] 5 Mar 2021 solid-like behaviors [17,27,28].
We conduct NVE molecular dynamics simulations, using LAMMPS [29], of a d = 2 dimensional polydisperse system of particles first introduced in [30]. The swap Monte Carlo algorithm [31,32], allows for the equilibration of this system even well below the mode-coupling theory temperature T MCT [30]. The particles have random diameters σ ∈ [σ min , σ max ] with probability density function ∝ σ −3 , where σ min and σ max are determined by the meanσ and coefficient of variation c σ of σ. Particles interact via a pair potential V (r) = GΘ(r 0 −r) r −12 + c 0 + c 2r 2 + c 4r 4 . Here, Θ is the Heaviside step function and c 0 , c 2 , and c 4 are chosen such that V (r 0 ) = 0 and V (r) is continuously differentiable. The argumentr = r/σ is the interparticle separation r between particles i and j with diameters σ i and . Following [30], we chooser 0 = 2.5, = 0.2, and c σ = 0.23. We set G,σ, the (uniform) particle mass m, and the Boltzmann constant k B to unity, thus defining our energy, length, mass, and temperature units. We study N = 10 4 particles in an L × L periodic square box of side length L = 100. For this system the onset temperature is T 0 ≈ 0.236 and the mode-coupling temperature is T MCT ≈ 0.116 [30].
We study rearrangements in 400 trajectories with independent initial velocities chosen from the Boltzmann distribution for each of 202 separate equilibrium configurations, for eight temperatures T = 0.100, 0.105, 0.110, 0.115, 0.120, 0.130, 0.140, and 0.150, doubling our statistics by averaging over trajectories in both directions in time [33]. We identify rearranging particles by focusing on quenched trajectories between inherent structures, and by measuring the squared displacement ∆r 2 between inherent structures [20,33]. Our criterion for rearrangement is ∆r 2 > 5a 2 across a time interval ∆t = 10 2 , where a is the plateau height of the root mean squared displacement for particles of similar size in the unquenched system at the given temperature [33]. The time scale ∆t = 10 2 matches the time at which dynamic facilitation is at play, identified as the one at which individual rearrangements cluster together in space [33]. For reference, the ballistic time is τ ballistic ≈ 10 0 in these units, and the relaxation time is τ α = 3 × 10 2 and τ α = 2 × 10 6 at the highest and lowest temperatures, respectively. For further characterization and details on the choice of ∆t, see the Supplementary Material [33].
Our first observation is that the rearrangement of a given particle causes elastic-like displacements of other particles; i.e. the system, despite being a liquid, displays a solid-like response at short times. To show this, we proceed as follows. For each rearrangement in the system, we obtain the local best-fit linear strain tensor E = E(r) for all particles from the symmetrized local quenched displacement gradient tensor [19,33,34]. We then extract strain invariants, such as the isotropic strain and (e,f) γext = E : eext ⊗ eext across a rearrangement interval ∆t = 10 2 as a function of the position r relative to a rearranging particle that was at r = 0 at the start of the rearrangement interval for systems at temperatures T = 0.100 (a,c,e) and T = 0.150 (b,d,f). Here, E is the best-fit linear strain tensor at r, E is detraced E, and eext is a unit eigenvector corresponding to the largest eigenvalue of E at the origin. The radial axes are plotted on a log scale. γ iso = TrE and the deviatoric strain . Crucially, E also gives us local extensional and compressional axes from its orthonormal eigenvectors: we denote as e ext (e com ) the eigenvector corresponding to the largest (smallest) eigenvalue of E. We then average over neighborhoods of rearranging particles translated such that the rearranging particle is at the origin at the start of the rearrangement interval (of duration ∆t) and rotated such that the extensional axis of its local strain tensor is e ext = (1, 0). We also calculate the strain in this direction, γ ext ≡ E : e ext ⊗ e ext , where e ext and e com are eigenvectors of the local strain tensor E computed at the rearranging particle i, whereas E is the local deviatoric strain tensor of the neighbors j of i. The dipolarity of γ iso , quadrupolarity of γ ext , and the isotropic behavior of γ dev , shown in Fig. 1 for temperatures T = 0.100 (left) and T = 0.150 (right) spanning the full temperature range of our study, are the expected signatures of the response of an elastic solid [35]. We now assess the temperature dependence of the range of the elastic-like response by studying the angular Fourier series coefficients, which characterize the strength of the observed symmetries reported above [33]. More specifically, we consider the zeroth and second order coefficientsγ iso,0 andγ iso,2 , the zeroth order coefficient γ dev,0 and the fourth order coefficientγ ext,4 . As shown in Fig. 2,γ iso,2 andγ ext,4 (Fig. 2(c,d)) decay for T ≤ 0.12 as r −2 from r = 3 up to a cutoff that increases slightly with decreasing T but is roughly r 10. This power-law decay is expected for an elastic response [35]. The radial dependence ofγ dev,0 ( Fig. 2(b)) deviates from the elastic response behavior at higher temperatures but appears to approach the expected behavior (dashed line) with decreasing temperature. This can be explained by the decreasing density of rearrangements per unit time with decreasing T , which allows the elastic displacement field from rearrangements to extend further into the material without being obscured by other rearrangements [33]. Finally, for the zeroth mode of γ iso (Fig. 2(a)), steadystate elasticity predictsγ iso,0 ≡ 0. Instead, we see a rapidly-decaying radial wave of isotropic expansion and contraction, with expansion peaks (compression troughs) roughly coinciding with the peaks (troughs) of g (r) [33]. This is consistent with the results of [35] in the context of athermal quasistatic shear: rearrangements expand the structure at the shortest distances from the rearrangement, thus locally hardening the structure, but compress the structure slightly further out, with this hardeningsoftening pattern repeating out to large distances.
We now focus directly on dynamic facilitation and analyze its relation with the elastic responses characterized above. We use the pair distribution of rearrangements, g hop (r) = g hop (r, θ) (the average density of rearrangements at r given a rearrangement at the origin, divided by the bulk rearrangement density [33]), to investigate the extent to which facilitation of rearrangements is due to elasticity. In Fig. 3(a,b) we displaỹ g hop (r, θ) = 2πg hop (r, θ)/ 2π 0 g hop (r, θ)dθ (the rearrangement pair distribution g hop (r, θ) normalized at each r by its angular mean 1 2π 2π 0 g hop (r, θ) dθ) for rearrangements over the time interval [0, ∆t] [36] at distances r = 1, 3, 10, 30, for temperatures T = 0.100 and T = 0.150. Well below T MCT , at T = 0.100 (Fig. 3(a)), the spatial distribution of subsequent rearrangements remains anisotropic out to r = 30. The angular Fourier series expansion coefficients ofg hop [33] show that this anisotropy comes from the second-and fourth-order Fourier modes ofg hop , and persists up to at least r ≈ 20. Above T MCT , at T = 0.150 ( Fig. 3(b)), no anisostropic structure can be discerned beyond the second neighbor shell, i.e. beyond r > 3.
To quantify the amount of anisotropy in the spatial distribution of rearrangements we introduce the mean squared anisotropy (deviation from isotropy) [37] α(r) = 1 2π 2π 0 (g hop (r, θ) − 1) 2 dθ, as shown in Fig. 3(c). High α at small r shows that displacements at this scale, though highly nonlinear, are measurably influenced by the orientation of the linear strain. This short-ranged influenceindicative of short-ranged dynamic facilitation-persists over the entire temperature range studied. As T is lowered, the large-r tail of α increases markedly, spanning two orders of magnitude, implying a significant difference in the role played by long-range elastic interactions between the lowest and highest temperatures considered here. At the same time, Fig. 3(d) shows that the pair correlation for hops, g hop (r), becomes progressively longerranged with decreasing temperature, indicating that rearrangements are triggered further and further away.
To quantify the extent of elastically-mediated dynamic facilitation, we introduce I α := 40 2.5 α(r)dr, which integrates α over distances beyond the second neighbor shell at r 2.5 (to filter out the effect of short-range facilitation) and up to the limit beyond which finitesize effects appear. Fig. 4(a) displays its evolution with 1/T . We observe a crossover at T MCT from a slow to a rapid growth of I α as a function of 1/T . Above T MCT , the strain-correlated spatial organization of the rearrangements is largely obscured, beyond the secondneighbor shell, by other rearrangements [33]. The competition between these other rearrangements and emergent elasticity at longer distances leads to the observed crossover in I α . Note that the observed crossover at T MCT is robust to choices of lower bound of the integral, rearrangement interval ∆t, strain length ξ FL , and starting point for the integration (r 0 ≥ 2.5) [33]. We similarly define an integrated excess density of rearrangements (relative to the bulk rearrangement density), I hop := 40 2.5 2πr (g hop (r) − 1) dr, where g hop (r) = 1 2π 2π 0 g hop (r, θ)dθ. Fig. 4(b) shows that I hop increases with decreasing T , showing that hops trigger increasingly distant hops. These results for I α and I hop provide direct and explicit evidence that the elastic response to a rearrangement triggers other, distant rearrangements, demonstrating the emergence of elastic dynamic facilitation when crossing T MCT .
Taken together, our findings suggest a crossover from solely short-ranged dynamic facilitation for T > T MCT to additional long-ranged, elastically-mediated dynamic facilitation of avalanches at T < T MCT . This finding provides an intriguing link between mean-field models, where T MCT marks the location where supercooled liquids first exhibit a transient solid-like behavior, and real-space dynamical heterogeneity. We note that the crossover at T MCT is not visible in the strain field induced by a rearrangement; for example, the isotropic strain shows dipolar angular behavior at all temperatures considered. The crossover is only evident in the spatial and angular distribution of rearrangements that follow soon after the original rearrangement at the origin. The anisotropic spatial distribution of hops evidenced in g hop , Fig. 3(a,b), and the increase in range of g hop (r), Fig. 3(d), as well as the crossover at T = T MCT seen in I α , Fig. 4(a), can only be due to the strain field created by the rearrangement at the origin. At high T strain does not facilitate rearrangements because the rearrangements are so abundant and closely spaced that their strain fields effectively interfere [33]. We note that at temperatures below T MCT , dynamical heterogeneity still involves short-ranged facilitation, whereby the local relaxation of clusters coalesces on longer length scales to form avalanches. Nevertheless, the emergence of prominent long-ranged elastic effects must influence the spatiotemporal correlations between such avalanches. Approximately ten years ago it was noticed that rare, long-ranged "surges" of dynamical heterogeneity connected to strain deformation occur in supercooled liquids [38], and that such behavior occurs in some coarsegrained kinetically-constrained models [3,5]. Our findings potentially pinpoint the microscopic underpinnings of such behavior, which could lead to greater understanding of the building blocks of phenomenological models of supercooled liquids. It would be certainly worthwhile to generalize elastoplastic models [16] to describe elastic dynamic facilitation within equilibrium thermal dynamics along the lines of [39,40]. More generally, our results highlight the importance of long-ranged elastic processes in mediating the organization of dynamical heterogeneities on long length and time scales, and the need to incorporate such processes in simplified models of the glass transition.
We thank M. Ozawa and L. Berthier for providing us with equilibrated configurations and the code to generate them, C. Scalliet for providing the swap potential files for use with LAMMPS, and E. Corwin for sharing computational resources with us in the early stages of this project. RNC thanks S. Ridout, G. Zhang, and I. Tah for discussions. This work received funding from the Simons Collaboration "Cracking the glass problem" via 454935 (G. Biroli), 454945 (A. J. Liu), 454951 (D. R. Reichman) and 348126 (R. N. Chacko), and from a Simons Investigator grant (327939 to A. J. Liu). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) [41], which is supported by National Science Foundation grant number ACI-1548562.