# Thirty milliseconds in the life of a supercooled liquid

Camille Scalliet, Benjamin Guiselin, Ludovic Berthier

► To cite this version:

Camille Scalliet, Benjamin Guiselin, Ludovic Berthier. Thirty milliseconds in the life of a supercooled liquid. Annual Review of Physical Chemistry, 2022, 12 (4), pp.041028. 10.1103/PhysRevX.12.041028. hal-03915196

HAL Id: hal-03915196

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

Submitted on 13 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Thirty Milliseconds in the Life of a Supercooled Liquid

Camille Scalliet,¹ Benjamin Guiselin,² and Ludovic Berthier³,⁴,*

¹Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, United Kingdom

²ENSL, CNRS, Laboratoire de physique, F-69342 Lyon, France

³Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, 34095 Montpellier, France

⁴Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom

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

(Received 21 July 2022; revised 25 September 2022; accepted 12 October 2022; published 12 December 2022)

We combine the swap Monte Carlo algorithm to long multi-CPU molecular dynamics simulations to analyze the equilibrium relaxation dynamics of model supercooled liquids over a time window covering 10 orders of magnitude for temperatures down to the experimental glass transition temperature $T_g$. The analysis of several time correlation functions coupled to spatiotemporal resolution of particle motion allow us to elucidate the nature of the equilibrium dynamics in deeply supercooled liquids. We find that structural relaxation starts at early times in rare localized regions characterized by a waiting-time distribution that develops a power law near $T_g$. At longer times, relaxation events accumulate with increasing probability in these regions as $T_g$ is approached. This accumulation leads to a power-law growth of the linear extension of relaxed domains with time with a large, temperature-dependent dynamic exponent. Past the average relaxation time, unrelaxed domains slowly shrink with time due to relaxation events happening at their boundaries. Our results provide a complete microscopic description of the particle motion responsible for key experimental signatures of glassy dynamics, from the shape and temperature evolution of relaxation spectra to the core features of dynamic heterogeneity. They also provide a microscopic basis to understand the emergence of dynamic facilitation in deeply supercooled liquids and allow us to critically reassess theoretical descriptions of the glass transition.

DOI: 10.1103/PhysRevX.12.041028

Subject Areas: Chemical Physics, Statistical Physics

## I. INTRODUCTION

There exists a large corpus of experimental studies analyzing the physical properties of supercooled liquids undergoing a glass transition [1–3]. Glassy systems exhibit well-established signatures characterizing their thermodynamic, rheological, and dynamic properties. An important goal of this experimental quest is to develop a sufficiently precise understanding of the physical behavior of liquids undergoing a glass transition to guide and constrain theoretical developments. A successful theoretical framework should explain the observed behaviors with precise assumptions that can be directly tested by experiments [4,5]. This program is not yet complete, and different theoretical explanations remain able to account for experimental results using hypotheses that can be difficult to validate experimentally [6–10].

Computer simulations have an important role to play in this endeavor, as they offer by construction a complete spatiotemporal resolution of glassy dynamics and the possibility to measure observables which are difficult or impossible to access experimentally [11,12]. For a long time, a major obstacle was the inability to study realistic models of glassy liquids in the temperature regime relevant to experiments. The situation changed radically five years ago when the swap Monte Carlo algorithm [13] was optimized and novel glass models were developed [14–16]. The swap Monte Carlo algorithm employs unphysical particle moves to accelerate the equilibration of supercooled liquids and can reach equilibrium states down to the experimental glass transition temperature $T_g$ or even below. This algorithmic development allowed progress regarding the analysis of structural and thermodynamic properties of liquid states [17–21], as well as characterization of the glass below $T_g$ [22–28]. However, because it employs unphysical particle motion, the swap Monte Carlo algorithm cannot be used to directly analyze the dynamics of supercooled liquids near $T_g$. To date, simulations of the dynamics were mostly performed in a relatively high-temperature regime corresponding, at best, to an average relaxation time up to 1 µs when converted into experimental units [29–31]. This

<footnote>*Corresponding author.
ludovic.berthier@umontpellier.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
corresponds to a temperature scale near the mode-coupling temperature crossover $T_{MCT} > T_g$. Fortuitously, experiments performed with colloidal particles cover a similar dynamic range [32,33]. As a result, a detailed microscopic characterization of the dynamics in the mode-coupling regime $T \ge T_{MCT}$ is available [34,35]. There exists ample evidence that physics in the regime $T < T_{MCT}$ may be of a different nature [36–40], but a thorough numerical exploration of glassy dynamics deep in this temperature regime is currently lacking.

Here we show that month-long multi-CPU molecular dynamics simulations, started from configurations that are first equilibrated using the swap Monte Carlo algorithm, open a novel window to analyze the equilibrium dynamics of deeply supercooled liquids in the regime $T_g < T < T_{MCT}$. In practice, we follow the equilibrium dynamics over ten decades in time with particle-scale resolution at temperatures down to $T_g$. Converted to experimental units, this approach allows us to follow the entire structural relaxation at temperatures well below the mode-coupling crossover $T < T_{MCT}$, up to relaxation times of $\tau_\alpha \approx 10$ ms. For even lower temperatures $T \approx T_g$, we study the first ten decades of the structural relaxation up to a maximal timescale of about 30 ms for our longest simulations.

We explore the relaxation dynamics of model supercooled liquids in a regime that was too difficult or impossible to access before. This important numerical effort has led to a previous work [41], where we concentrated on the emergence of excess wings in dynamic spectra obtained at low temperatures. Our ambition and focus here are very different as we provide a complete view of all relaxation processes from microscopic times up to timescales several times longer than the relaxation time.

Our philosophy in this article is to first report our numerical observations and quantify them with as little interference from theoretical models as possible. We pay special attention to the regime below $T_{MCT}$ that has not been accessed before and contrast our findings with earlier work at higher temperatures. In the final part of the manuscript only, we critically compare our results to existing theoretical frameworks. For the two models studied, a clear physical picture of the structural relaxation emerges at low temperatures, which stems from early relaxation events that are broadly distributed followed by increasingly correlated motions in space and time. These observations account for the emergence of dynamic facilitation in the dynamics of deeply supercooled states.

Our manuscript is organized as follows. In Sec. II, we define our computer models and numerical strategy. In Sec. III, we present results concerning ensemble-averaged time correlation functions. In Sec. IV, we offer visualization of the relaxation dynamics over a broad range of timescales, length scales, and temperatures. In Sec. V, we analyze in more detail the early times of the relaxation. In Sec. VI, we explain how the structural relaxation unfolds

from early to large times. In Sec. VII, we provide a discussion of our results, comparing them with earlier numerical work and theoretical views.

## II. COMPUTER MODELS AND METHODS

### A. Glass-forming models

We study size-polydisperse mixtures of soft repulsive spheres in two and three spatial dimensions $d = 2, 3$. These two models have been shown to be representative computational glass formers [17,18] and extensively studied before [23,26,27,42]. The particle diameters $\{\sigma_i\}$ are drawn from the probability distribution $P(\sigma) = A/\sigma^3$ with A a normalization constant with bounds $\sigma_{min} \le \sigma_i \le \sigma_{max}$. Two particles i and j at positions $\mathbf{r}_i$ and $\mathbf{r}_j$ and separated by a distance $r_{ij} = |\mathbf{r}_{ij}| = |\mathbf{r}_i - \mathbf{r}_j|$ interact within a cutoff $r_{ij}/\sigma_{ij} < x_c = 1.25$ with a repulsive interaction potential:

$$v(r_{ij}) = \epsilon \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} + c_0 + c_2 \left(\frac{r_{ij}}{\sigma_{ij}}\right)^2 + c_4 \left(\frac{r_{ij}}{\sigma_{ij}}\right)^4, \quad (1)$$

where the constants $c_0 = -28\epsilon/x_c^{12}$, $c_2 = 48\epsilon/x_c^{14}$, and $c_4 = -21\epsilon/x_c^{16}$ ensure continuity of the potential and of its first two derivatives at the cutoff $x_c$. We employ a non-additive mixing rule $\sigma_{ij} = 0.5(\sigma_i + \sigma_j)(1 - \eta|\sigma_i - \sigma_j|)$ to avoid fractionation and crystallization at low temperatures [14]. All particles have an equal mass $m$. We use the average diameter $\bar{\sigma}$ as unit length, $\epsilon$ as unit energy with the Boltzmann constant set to unity, and $\sqrt{m\bar{\sigma}^2/\epsilon}$ as unit time. In these units, we choose $\sigma_{min} = 0.73$, $\sigma_{max} = 1.62$, and $\eta = 0.2$.

We simulate N particles at number density $\rho = N/L^d = 1$ in a cubic or square box of linear size L with periodic boundary conditions. The results reported for the 3D model are obtained for systems of N = 1200 particles. We perform simulations of a larger system with N = $10^4$ particles to check for finite-size effects in the dynamics and for visualization purposes. Measurements in the 2D model are reported for N = 2000 particles. Some 2D simulations are performed with N = $10^4$ particles as well, in particular to generate snapshots.

### B. Monte Carlo and molecular dynamics simulations

We use a hybrid scheme combining Monte Carlo (MC) particle-swap moves and molecular dynamics (MD) in order to generate equilibrium configurations [15]. The algorithm alternates between ordinary MD simulation sequences during which the particle positions and velocities evolve with a fixed particle diameter, and Monte Carlo sequences during which particle diameters are swapped at fixed positions and velocities. The MD sequences last $t_{MD}$ and take place at constant temperature T imposed by a Nosé-Hoover thermostat. To perform MC sequences, the positions and velocities are frozen and
$n_{swap}N$ swap moves are attempted. During a swap move, two particles are randomly selected and an exchange of their diameters is proposed. After computing the change in potential energy, the move is accepted following the Metropolis acceptance rule at temperature $T$, which ensures equilibrium sampling. Following Ref. [15], we optimize $t_{MD}$ and $n_{swap}$ for maximal efficiency. Using this hybrid algorithm, we prepare a large number of statistically independent configurations $N \in [100, 500]$ for temperatures ranging from the onset of glassy dynamics $T_o$ down to the experimental glass transition temperature $T_g$. These characteristic temperature scales are defined and numerically determined below.

We take these equilibrium configurations as initial conditions for standard MD simulations without swap MC moves. We perform microcanonical NVE simulations in 3D and canonical NVT simulations with a Nosé-Hoover thermostat in 2D, using in both cases a time step $dt = 0.01$ for the numerical integration of the equations of motion. These different choices are made for convenience only, as the statistical ensemble does not influence the local dynamics. The duration of most simulations is $10^7$, corresponding to $10^9$ MD steps and a simulation time of about one week. Ensemble-averaged observables are computed by averaging over $N$ independent trajectories. For selected state points in 2D, we increase the simulation time to $8.2 \times 10^8$ for $N = 10^4$ particles, which was reached by using a parallel code running over 24 processors (CPUs) for two and a half months.

### C. Why 30 ms?

The temperatures and timescales over which equilibrium dynamics can be probed with MD simulations are limited, as recently reviewed in Ref. [30]. Molecular dynamics simulations last at most $10^{10}$ MD steps for a simple pair potential simulated over several weeks and a reasonable system size ($N \sim 10^3-10^4$). This translates into a maximum simulation time $t_{max} \sim 10^8$, considering a typical discretization time step of order $10^{-2}$. Without the swap Monte Carlo algorithm, the simulated time is necessarily split between equilibration and production runs. To achieve both a proper equilibration and a significant exploration of the configuration space in equilibrium conditions, the MD equilibration run must last at least about 100 times the averaged relaxation time $\tau_\alpha$. This sets an upper limit $\tau_\alpha \le 10^6$ on the accessible relaxation times even using extensive simulations and an efficient implementation of the MD. Most numerical studies have therefore investigated the dynamics of supercooled liquids with $\tau_\alpha \le 10^5$.

Our numerical strategy completely circumvents the need for time-consuming equilibration runs by exploiting the huge equilibration speedup afforded by the swap Monte Carlo algorithm. This idea is in its infancy [31,41,43]. In this approach, preparing multiple equilibrated

configurations down to $T_g$ is a simple task, and we can thus ensure both equilibration and statistical accuracy of the results. This implies that the computational time can be entirely spent on the production runs to simulate the dynamics of low-temperature states over a time window limited to $t_{max} \approx 10^8$. Crucially [31], this time window is available even at temperatures where $\tau_\alpha \gg t_{max}$, which were inaccessible in equilibrium conditions in previous work.

Following earlier works [14,31,42], we translate simulation timescales into experimental ones using the structural relaxation time $\tau_o$ at the onset temperature $T_o$ as a reference timescale. For a broad range of molecular liquids [44], one measures $\tau_o \approx 10^{-10}$ s. For the computer models studied here, one finds $\tau_o \approx 3$ in simulation units. The production runs over which ensemble-averaged quantities are measured last $10^7 \approx 3 \times 10^6 \tau_o$ in simulation units, which translates into $3 \times 10^6 \tau_o \approx 0.3$ ms in experimental units. Below, we also present the results of very long simulations which last $t_{max} = 8.2 \times 10^8$. In experimental units, this corresponds to observing the relaxation dynamics of a supercooled liquid over about 30 ms, as announced.

### III. EQUILIBRIUM RELAXATION DYNAMICS DOWN TO THE EXPERIMENTAL GLASS TRANSITION

#### A. Time correlation functions

We first characterize the equilibrium dynamics using ensemble-averaged time correlation functions. In 3D, we use the self-intermediate scattering function

$$F_s(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} \cos[\mathbf{q} \cdot \delta\mathbf{r}_i(t)] \right\rangle, \quad (2)$$

where $\delta\mathbf{r}_i(t) = \mathbf{r}_i(t) - \mathbf{r}_i(0)$, and the brackets indicate an average over the $N$ independent runs at temperature $T$. We also perform an angular average over wave vectors with $|\mathbf{q}| = 6.9$ corresponding to the first peak in the total structure factor $S(q)$. We define the $\alpha$-relaxation time $\tau_\alpha(T)$ as $F_s(\tau_\alpha) = e^{-1}$.

In Fig. 1(a), we present equilibrium results for $F_s(t)$ in 3D at several temperatures. We concentrate on data in the unexplored low-temperature regime where $\tau_\alpha > 10^4$. The highest temperature shown, $T = 0.095$, roughly corresponds to the lowest temperature explored in the majority of earlier studies, while the lowest temperature shown is close to the experimental glass transition temperature $T_g$. This broad temperature range concretely demonstrates how our numerical strategy opens a novel temperature window to study the dynamical slowdown in supercooled liquids.

All correlation functions display a fast decay at very short time $t \approx 1$ due to thermal motion within the arrested amorphous structure. At very low temperature, elasticity and inertial dynamics give rise to a weak oscillatory decay
toward the plateau. After an extended plateau regime which becomes more extended at lower temperature, full decorrelation is eventually observed taking the form of a stretched exponential. As the temperature decreases, structural relaxation shifts to longer times, and it is no longer observed for temperatures $T \lesssim 0.075$. At the lowest temperature $T \lesssim 0.059$, we observe an extended plateau covering about seven decades in time. More careful inspection reveals that the plateau is not strictly constant but decays extremely weakly with time, as we illustrate more clearly below. Recall that despite the absence of decorrelation in the numerical time window, all correlation functions are representative of equilibrium dynamics.

In 2D, collective long-ranged fluctuations give rise to particle displacements which affect the behavior of $F_s(t)$ in an unwanted way [45,46]. We measure the bond-orientational correlation function $C_\psi(t)$ which is sensitive to changes in the local environment of particles [47,48]. We define the sixfold order parameter of particle i,

$$ \Psi_i(t) = \frac{1}{n_i(t)} \sum_{j=1}^{n_i(t)} e^{i6\theta_{ij}(t)}, \quad (3) $$

where $n_i(t)$ is the number of neighbors it has at time t defined as particles j with $r_{ij} < 1.45$. The cutoff corresponds to the first minimum in the total pair distribution function, and we checked that alternative definitions of neighbors via a Voronoi tessellation or a solid-angle-based method [49] lead to similar results. The angle $\theta_{ij}(t)$ is defined between the x axis (without loss of generality) and $r_{ij}(t)$. The bond-orientational correlation function then reads

$$ C_{\psi}(t) = \left\langle \frac{\sum_{i} \Psi_{i}(t) \Psi_{i}^{*}(0)}{\sum_{i} |\Psi_{i}(0)|^{2}} \right\rangle, \quad (4) $$

where the star denotes the conjugate complex and the sums run over all particles $i = 1, \dots, N$. We extract the $\alpha$-relaxation time $\tau_\alpha(T)$ as $C_\psi(\tau_\alpha) = e^{-1}$.

In Fig. 1(b), we present equilibrium results for $C_\psi(t)$ in 2D at several temperatures selected as in 3D below the mode-coupling crossover. Clearly, the time and temperature evolution of the bond-orientational correlation function is qualitatively similar to that of the self-intermediate scattering function in 3D shown in Fig. 1(a).

## B. Relaxation times and temperature scales

We determine three temperature scales relevant to describe the dynamic slowdown in supercooled liquids: the onset temperature $T_o$, the mode-coupling crossover temperature $T_{MCT}$, and the experimental glass transition temperature $T_g$. We first provide the value of these temperature scales and detail below how we estimate them. We obtain $T_o = 0.2$, $T_{MCT} = 0.095$, and $T_g = 0.056$ in 3D,

<figure><img src="image_3.png" /><figcaption>FIG. 1. Equilibrium time correlation functions characterizing relaxation dynamics of deeply supercooled liquids. (a) Self-intermediate scattering function $F_s$ in 3D. (b) Bond-orientational correlation function $C_\psi$ in 2D. In both panels, we show only temperatures between the mode-coupling crossover and the experimental glass transition temperature, $T_g \le T \le T_{MCT}$.</figcaption></figure>

and $T_o = 0.2$, $T_{MCT} = 0.12$, and $T_g = 0.07$ in 2D. These characteristic temperatures will be useful to interpret and organize our results.

In Figs. 2(a) and 2(b), we present the relaxation time $\tau_\alpha$ for the 3D and 2D models, respectively. At the onset temperature of glassy dynamics $T_o$, where the relaxation time equals $\tau_o$, the relaxation time departs from its high-temperature Arrhenius dependence, $\tau_\alpha(T) \propto e^{E_\infty/T}$. The energy scale $E_\infty$ is equal to $E_\infty = 0.23$ (in 3D) and $E_\infty = 0.7$ (in 2D), and $\tau_o \approx 3$ (in $d = 2, 3$). In Fig. 2, we focus on the supercooled regime and report $\tau_\alpha$ normalized by $\tau_o$ as a function of the inverse temperature. In this representation, Arrhenius behavior translates into a straight line. In the simulations, we directly measure the relaxation in the range $\tau_\alpha < 10^7$, equivalently, $\log_{10}(\tau_\alpha/\tau_o) < 6.5$. These measurements are labeled as “Data” in Fig. 2.

Then, we locate the mode-coupling crossover temperature $T_{MCT}$ by using a power-law fit $\tau_\alpha(T) \sim (T - T_{MCT})^{-\gamma}$ in the regime $0 \le \log_{10}(\tau_\alpha/\tau_o) \le 3$ [50] using the exponents $\gamma = 2.7$ (in 2D) and $\gamma = 2.5$ (in 3D). We find $\tau_\alpha(T_{MCT}) \sim 10^4 \tau_o$ in 2D and 3D.

Direct measurements of $\tau_\alpha$ stop at temperatures where the correlation functions do not reach $e^{-1}$ in the numerical time window. This occurs roughly five decades
<figure><img src="image_4.png" /><figcaption>FIG. 2. Averaged relaxation times in 3D and 2D glass-forming models. Relaxation time $\tau_{\alpha}$ rescaled by its value $\tau_{o}$ at the onset temperature, as a function of the inverse temperature in 3D (a) and 2D (b). The data points are directly measured from Fig. 1, and extended over 1.5 decades using time-temperature superposition (TTS). An Arrhenius fit (full line) extrapolates $\tau_{\alpha}/\tau_{o}$ to the value $10^{12}$ to locate the experimental glass transition temperature $T_{g}$ (dashed line). We concentrate on the regime between $T_{MCT}$ (dash-dotted line) and $T_{g}$.</figcaption></figure>

before the extrapolated glass transition $T_{g}$, which is defined by $\log_{10}[\tau_{\alpha}(T_{g})/\tau_{o}] = 12$. Still, at temperatures where $\tau_{\alpha} < 10^{7}$, the final decay is well fitted by a stretched exponential form $F_{0}e^{-(t/\tau_{\alpha})^{\beta}}$. The stretching exponents $\beta \approx 0.56$ (in 3D) and $\beta \approx 0.6$ (in 2D) are nearly temperature independent within statistical uncertainty, and the amplitude $F_{0}$ increases weakly with decreasing temperature. This indicates that time-temperature superposition (TTS) is well obeyed in our systems. Fixing $\beta$ to a constant value, we can use TTS to estimate $\tau_{\alpha}$ when it falls outside our numerical window [31]. In practice, we can safely extend our measurements of $\tau_{\alpha}(T)$ by an additional 1.5 decades, as indicated by the “TTS” points in Fig. 2, which then reach $\log_{10}(\tau_{\alpha}/\tau_{o}) \approx 8$ (i.e., about 10 ms).

To estimate the experimental glass transition temperature $T_{g}$, we must extrapolate the relaxation time data by 4 additional orders of magnitude. To do so, we describe $\tau_{\alpha}$ over the remaining decades using an Arrhenius fit $\tau_{\alpha}(T) \propto e^{E_{A}/T}$, with $E_{A} = 2.67$ (in 3D) and $E_{A} = 2.97$ (in 2D), and locate $T_{g}$ where $\log_{10}[\tau_{\alpha}(T_{g})/\tau_{o}] = 12$. Although the temperature dependence of $\tau_{\alpha}$ is not purely Arrhenius over the entire numerical window, our choice of an Arrhenius extrapolation which neglects fragility provides a lower bound to the correct value of $T_{g}$.

The swap Monte Carlo algorithm allows us to easily prepare equilibrium configurations even at temperatures close to the determined $T_{g}$. We can thus safely claim that we analyze the equilibrium dynamics of liquids down to the experimental glass temperature $T_{g}$. The curves at the lowest temperatures shown in Fig. 1 correspond to the dynamics at, or very close to, the experimental glass transition

temperature where we can access the first 30 ms of the relaxation dynamics in our longest simulations.

### C. Mean-squared displacements: E pur si muove

To characterize the average motion of the particles, we compute the mean-squared displacement (MSD):

$$ \Delta(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} |\mathbf{r}_{i}(t) - \mathbf{r}_{i}(0)|^{2} \right\rangle. \qquad (5) $$

The full lines in Fig. 3 present the time dependence of the MSD at various temperatures in 3D. The data are not presented in 2D because of the large and spurious collective fluctuations mentioned above. When the temperature is not too low, the curves exhibit the usual time dependence with a ballistic regime at early times followed by a long plateau at intermediate times before finally entering a diffusive regime at very long times which we can observe down to $T = 0.09$. As the temperature decreases, the intermediate plateau extends over longer timescales, and the diffusive regime eventually shifts outside the numerical time window. Close to $T_{g}$, we can observe only the very beginning of the escape from the plateau, mirroring the behavior of the self-intermediate scattering function in Fig. 1(a). The MSD at long times reaches the small value $\Delta \approx 0.03$ corresponding to an average particle displacement of about one fifth of the average particle diameter. In other words, the system appears totally frozen on the very long timescale explored by the simulations. And yet, as we see, particles move.

Indeed, a more subtle picture emerges by filtering thermal motion [51,52]. We use a conjugate gradient algorithm to bring the configurations explored dynamically $\{\mathbf{r}_{i}(t)\}_{i=1,...,N}$ to their nearest potential energy minimum, also called inherent structure (IS) $\{\mathbf{r}_{i}^{\text{IS}}(t)\}_{i=1,...,N}$. From these, we compute a version of the MSD $\Delta^{\text{IS}}(t)$, where the erratic thermal motion of particles exploring their cages no

<figure><img src="image_5.png" /><figcaption>FIG. 3. Mean-squared displacement in 3D deeply supercooled states. From $T_{MCT}$ to $T_{g}$ (left to right) in the normal dynamics $\Delta(t)$ (full lines) and in the inherent structures (IS) $\Delta^{\text{IS}}(t)$ (symbols). The legend is as in Fig. 1(a). The dash-dotted line corresponds to $\Delta = 6 \times 10^{-3}$ used in Fig. 7 below.</figcaption></figure>
longer contributes. When particles simply fluctuate around their average position, this procedure returns a vanishing value for $\Delta^{\text{IS}}$, while $\Delta$ plateaus at the Debye-Waller factor.

The data for $\Delta^{\text{IS}}$ shown in Fig. 3 paint a completely different picture. The two MSD curves coincide only in the long-time diffusive regime for large particle displacements, where the small contribution due to thermal vibrations becomes negligible. At early times in the ballistic regime, we have $\Delta^{\text{IS}} \ll \Delta$ as previously reported [52]. Surprisingly, we observe that close to $T_g$, $\Delta^{\text{IS}}(t)$ develops a nontrivial time dependence as soon as $t > 1$. Although the average structural relaxation time is close to $\tau_\alpha = 10^{12}$, there are already nontrivial particle motion and relaxation taking place at times that are many orders of magnitude shorter. We conclude that the extended plateau in the MSD in fact masks early nontrivial relaxation events which take place much before the average structural relaxation time. These motions would be difficult to detect when using the self-intermediate scattering function, as pioneered in Ref. [51]. The nontrivial time evolution of the MSD when measured in inherent structures also shows that a large number of ISs are explored over the time window preceding the structural relaxation, which raises questions about the relevance of the potential energy landscape to describe the $\alpha$-relaxation dynamics of deeply supercooled liquids.

These observations constitute a central theme of this work: By resolving particle motion in equilibrium conditions, we reveal how the structural relaxation unfolds over an extended period of time between $t \approx 1$ and $t \gg \tau_\alpha$ at very low temperatures. In the rest of the article, we analyze the corresponding particle motion and develop a physical understanding of the structural relaxation in deeply supercooled liquids.

# IV. QUALITATIVE OVERVIEW OF THE RELAXATION DYNAMICS

## A. Direct visualization of the relaxation

To resolve relaxation dynamics at the particle scale, we use the bond-breaking correlation function $C_B^i(t)$ [53,54]. By definition, this quantity measures the fraction of the neighbors of particle $i$ defined at time $t = 0$ that are still neighbors later at time $t$,

$$C_B^i(t) = \frac{n_i(t|0)}{n_i(0)}, \qquad (6)$$

where $n_i(0)$ is the number of neighbors of particle $i$ at time $t = 0$, and $n_i(t|0)$ the number of those particles that are still neighbors after time $t$. At $t = 0$, neighbors are defined via the criterion $r_{ij}/\sigma_{ij} < 1.485$ (in 3D) or $r_{ij}/\sigma_{ij} < 1.3$ (in 2D), which corresponds to the first minimum in the rescaled pair distribution function constructed using distances $r_{ij}$ rescaled by $\sigma_{ij}$. At $t > 0$, we define neighbors via $r_{ij}/\sigma_{ij} < 1.7$ (in $d = 2, 3$), which is below the second peak

in the rescaled pair distribution function. This slightly larger cutoff avoids spurious noise due to particles frequently exiting or entering the neighbor shell due to thermal vibrations [48].

With these choices, $C_B^i(t)$ decays only whenever a rearrangement involving particle $i$ takes place, and is thus an excellent indicator of local relaxation events [53], although not being a single-particle quantity. Thresholding single-particle displacements (physical displacements or the ones computed from the visited inherent states) yields qualitatively similar pictures, but the definition of local relaxation may become blurred, especially at the highest studied temperatures. The correlation $C_B^i(t)$ thus tracks the cumulative effect of rearrangements taking place around particle $i$ over time. Particles which have lost half of their initial neighbors $C_B^i \le 0.5$ are called “mobile” [55]. Particles with $C_B^i = 1$ have not rearranged and are called “quiescent.” We define

$$C_B(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} C_B^i(t) \right\rangle \qquad (7)$$

as the ensemble average of $C_B^i(t)$.

We illustrate the spatiotemporal evolution of the relaxation dynamics in a two-dimensional liquid in Fig. 4. The corresponding videos are provided in Supplemental Material [56]. The temperature is $T = 0.09$, well below $T_{\text{MCT}} = 0.12$. This temperature is the lowest one at which we can observe full decorrelation within the simulation time $t_{\text{max}} = 8.2 \times 10^8$ equivalent to 30 ms, since the relaxation time defined by $C_\psi(\tau_\alpha) = e^{-1}$ is $\tau_\alpha = 1.5 \times 10^8$. On a logarithmic scale, this temperature is thus roughly halfway between $T_{\text{MCT}}(\tau_\alpha/\tau_o \approx 10^4)$ and $T_g(\tau_\alpha/\tau_o \approx 10^{12})$. We show snapshots of the liquid at six times along the trajectory $t_1, \dots, t_6 = 2 \times 10^6, 10^7, 4.6 \times 10^7, 1.1 \times 10^8, 2.3 \times 10^8, 8.3 \times 10^8$. The particles are colored according to the value of the bond-breaking correlation at time $t_j$ using a blue to red code for $C_B^i(t_j) = 1$ to 0.

At time $t_1$, two decades before $\tau_\alpha$, the average correlation $C_B(t_1) = 0.97$ is still very close to 1. While most particles are quiescent (blue), we distinguish a few mobile particles (white or red) clustered into sparse, compact, and localized domains. These clusters typically contain a few red particles, which have $C_B^i(t_1) = 0$. These particles have undergone a series of nontrivial rearrangements which lead to a complete renewal of their local environment. Particles at the boundary of these clusters typically appear white, indicating $C_B^i(t_1) = 0.5$: They have lost half of their initial neighbors (those inside the mobile cluster). The images display the raw data for $C_B^i$, with no spatial averaging. The sharp separation observed between mobile and quiescent particles is thus a genuine physical effect. When the average correlation is $C_B = 0.5$ near time $t_4$, we see that in real space most particles have either $C_B^i = 1$ or $C_B^i = 0$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Thirty milliseconds in the life of a supercooled liquid. Spatiotemporal evolution (left to right, top to bottom) of the structural relaxation in a 2D deeply supercooled liquid at $T = 0.09$, where $\tau_{\alpha} = 1.5 \times 10^8$. Particle colors indicate the bond-breaking correlation $C_B^i(t_j)$ from blue (quiescent) to red (mobile) following the color bar indicated top left. The snapshots are taken at times $t_1, \dots, t_6 = 2 \times 10^6, 10^7, 4.6 \times 10^7, 1.1 \times 10^8, 2.3 \times 10^8, 8.2 \times 10^8$, at which the average bond-breaking correlation equals 0.97, 0.86, 0.7, 0.54, 0.33, and 0.1. The linear size of the box is $L = 100$.</figcaption></figure>

with equal probability. This important property becomes
very clear at this very low temperature, as we demonstrate
shortly.

As time increases in Fig. 4, we observe that relaxation
proceeds via two distinct mechanisms. From one snapshot
to the next, we first observe the appearance of new relaxed
clusters in regions that are quiescent in the preceding
snapshot. Clearly, the emergence of such new clusters is a
stochastic process characterized by a time distribution that
we analyze in Sec. V D.

The second physical process leading to the appearance of
newly relaxed particles is via the coarsening of mobile
clusters which exist at time $t_j$ and get larger at time $t_{j+1}$.
The time dependence of the domain size is studied in
Sec. V D. At large times, the growth of relaxed regions
leads to their coalescence, so that the system exhibits a

bicontinuous structure of mobile or quiescent domains near
$\tau_{\alpha}$. This bicontinuous structure is obviously characterized
by a distribution of domain sizes that we investigate in
Sec. V D.

At very long times in Fig. 4, we observe small regions
that are still quiescent after $5.5\tau_{\alpha}$. They represent the most
stable regions of the structure at time $t=0$. Compared to
the substantial past and present research focusing on fast
relaxing regions and defects, very little is known about
these slow domains. Inspection of the last two snapshots
reveals that the slow domains are slowly invaded by the red
mobile regions and slowly shrink in size as time increases.
In Sec. V I E, we discuss the physics associated with these
slow domains and their relation with the lifetime of
dynamic heterogeneities.

B. Evolution with temperature

The physical picture of the structural relaxation
revealed by the images in Fig. 4 appears relatively simple
because there is a very sharp separation between mobile
and quiescent particles. However, such a clear distinction
emerges only when the temperature is low enough, and it
is not present in the temperature scale traditionally
studied in computer simulations. To illustrate this point,
we compare the pattern of structural relaxation at differ-
ent temperatures. The corresponding videos are provided
in Supplemental Material [56]. We keep the value of the
average correlation function $C_B(t)$ equal, as this corre-
sponds to a similar degree of decorrelation from the
initial structure.

We first investigate the very initial stages of the relax-
ation process at $t \ll \tau_\alpha$. We show in Fig. 5 three snapshots
of 2D liquids at a time $t$ where the average correlation is
$C_B(t) = 0.978$. The temperatures are equal, or close, to
three relevant temperatures: $T_o = 0.2$ (a), $T_{MCT} = 0.12$ (b),
and $T = 0.075 = 1.07T_g$ (c). Clearly, the picture changes
with the temperature. At the onset temperature $T_o$, many
particles turn to have lost a few neighbors and appear in
light blue. These particles are scattered throughout the
sample and do not form well-defined clusters. Moving to
the mode-coupling crossover temperature $T_{MCT}$, we dis-
tinguish relatively extended regions containing only qui-
escent particles. The particles which have undergone some
rearrangements start to cluster. Still, the clusters remain
rough, numerous, with ill-defined boundaries. At the lowest
temperature $T = 0.075$ closest to $T_g = 0.07$, the early stage
of relaxation looks qualitatively different. We observe a
small number of very mobile (red) particles localized into
compact clusters with clear boundaries. These particles
have undergone several relaxation events leading to
$C_B^i = 0$. The particles in the rest of the sample are all
quiescent and have not undergone a single rearrangement.

Turning to larger times where structural relaxation occurs
$C_B(t) = 0.5$, we again observe a dramatic temperature
dependence. The dynamics is almost spatially homogeneous
<figure><img src="image_7.png" /><figcaption>FIG. 5. Temperature evolution of dynamic heterogeneity in 2D. Particles are colored following the bond-breaking correlation $C_B^i$; see legend in Fig. 4. (a)–(c) Relaxation at early times $t \ll \tau_\alpha$ where $C_B(t) = 0.978$ at different temperatures. Snapshots are taken at $t = 6, 2 \times 10^3, 3.2 \times 10^8$ from (a) to (c). The circles indicate the point-to-set length scale $\xi_{PTS}$ discussed in Fig. 21. (d)–(f) Relaxation close to the $\alpha$ relaxation, where $C_B(t) = 0.5$. Snapshots are taken at $t = 500, 1.5 \times 10^5, 1.7 \times 10^8$ from (d) to (f). The linear size of the box is $L = 100$.</figcaption></figure>

at the onset temperature and becomes spatially correlated near $T_{MCT}$ where relatively well-defined domains have emerged, yet with large fluctuations within them and rough boundaries. As the temperature is lowered further, we observe again a much clearer contrast between relaxed and unrelaxed domains which become more compact and less irregular at their boundaries. Therefore, even though the average correlation $C_B(t)$ is the same within the top and bottom panels of Fig. 5, the distribution of local values $C_B^i(t)$ evolves from a single peak around the average at high temperatures to a nearly bimodal distribution, peaking at 0 and 1 at low temperature. We also note that the characteristic size of the heterogeneity shown in these plots changes a lot between $T_o$ and $T_{MCT}$ but appears to evolve more slowly below $T_{MCT}$. We also document this point below; see Fig. 21.

In Fig. 6, we show that similar conclusions are supported by our simulations in 3D, although snapshot rendering is more intricate. To ease visualization, not only particle colors but also diameters follow the bond-breaking correlation $C_B^i$: At short (resp, large) times, mobile (resp, quiescent) particles are shown with a larger diameter. In 3D, the emergence at very low temperatures of well-defined mobile clusters at short times is very clear, whereas a very diffuse pattern is observed at higher temperatures. At large times and low temperature, extended clusters of quiescent particles (blue) are visible. We instead notice that it is difficult to observe the dynamic heterogeneity field near the mode-coupling crossover in 3D.

<figure><img src="image_8.png" /><figcaption>FIG. 6. Temperature evolution of dynamic heterogeneity in 3D. (a)–(c) Relaxation at early times $t \ll \tau_\alpha$ where $C_B(t) = 0.978$, at different temperatures taken at $t = 6.5, 475$, and $1.8 \times 10^7$ from (a) to (c). (d)–(f) Relaxation close to the $\alpha$ relaxation, where $C_B(t) = 0.5$ at different temperatures taken at $t = 45, 2.7 \times 10^4, 1.2 \times 10^7$ from (d) to (f). The particles are colored following their bond-breaking correlation (see Fig. 4) and drawn with a diameter proportional to either $1 - C_B^i(t)$ (top) or $C_B^i(t)$ (bottom). The linear size of the box is $L = 21.5$.</figcaption></figure>

# V. EARLY TIMES OF THE STRUCTURAL RELAXATION

## A. van Hove distribution functions

To analyze particle motion at very short times, we first record the self-part of the van Hove distribution function

$$G_s(r, t) = \frac{1}{4\pi r^2} \langle \delta(r - |\mathbf{r}_i(t) - \mathbf{r}_i(0)|) \rangle \quad (8)$$

in the three-dimensional liquid; see Fig. 7. To follow the evolution with the temperature of this quantity, we select

<figure><img src="image_9.png" /><figcaption>FIG. 7. Heterogeneous dynamics at early times. van Hove distribution of particle displacements $G_s(r, t)$ in 3D at the time when $\Delta^{\text{IS}}(t) = 6 \times 10^{-3}$ (dash-dotted line in Fig. 3) for $T = 0.095, 0.08$, and 0.059 [same legend as in Fig. 1(a)]. The dash-dotted line is an exponential fit of the tail of the distribution at the lowest temperature. The inset compares the van Hove distributions computed at $T = 0.059$ either from real displacements or from the visited inherent structures (IS).</figcaption></figure>
times where the MSD in the inherent structure is constant $\Delta^{\text{IS}}(t) = 6 \times 10^{-3}$, which is indicated as a horizontal line in Fig. 3. This threshold corresponds to the value of $\Delta^{\text{IS}}$ as the MSD $\Delta$ departs from its plateau. As such, we probe the distribution of particle displacements at extremely early times compared to $\tau_{\alpha}$.

At all temperatures, the distribution displays a peak around $r=0$ caused by particles vibrating around their position at $t=0$. Although thermal vibrations occur on a length scale $r \approx 0.15$, the van Hove distribution extends to much broader values, with a tail which is well described by an exponential decay [57]. The exponential tail encompasses displacements that exceed the particle size and extends to an increasing range of particle displacements as the temperature lowers. When $T \le 0.07$, a peak emerges near $r \approx 0.9$, which coincides with the position of the first peak of the pair distribution function $g(r)$. Therefore, the emerging peak in $G_s(r, t)$ is due to particles hopping to the position previously occupied by one of their neighbors.

This behavior means that dynamic heterogeneity becomes more pronounced as $T$ decreases, even at $t \ll \tau_{\alpha}$. On these short timescales, a large majority of particles are quiescent. These particles have very small displacements due to thermal vibrations, and they populate the broad peak near $r=0$ in $G_s(r, t)$. However, a very small fraction of particles are already quite mobile and display displacements that are broadly distributed and reach several particle diameters long. These mobile particles give rise to the exponential tail in $G_s(r, t)$. In real space, they correspond to the rare clusters of relaxed particles detected at short times in Fig. 4. This small population of mobile particles is responsible for the slow growth of the MSD computed within inherent states in Fig. 3. This is confirmed in the inset of Fig. 7 which demonstrates that the exponential tail is unaffected by removing thermal fluctuations. Finally, the temperature evolution of the van Hove distribution also confirms that the contrast between mobile and quiescent particles becomes stronger at lower temperatures, which leads to the stronger contrast observed in the relaxation snapshots in Sec. IV.

## B. Complex particle motion inside isolated clusters

The van Hove distribution at short times and low temperatures shows that a small population of highly mobile particles coexists with a majority of particles which simply undergo thermal vibrations. Yet, the distribution itself provides no information as to how these mobile particles are organized in space and what type of rearrangements they undergo.

In Fig. 8, we highlight these mobile particles and their displacements for $T = T_{\text{MCT}}$, and for $T = 0.075$ near $T_g$ in the two-dimensional system. We represent with arrows the displacement vectors of the particles $\delta\mathbf{r}_i(t) = \mathbf{r}_i(t) - \mathbf{r}_i(0)$, at $t = 2 \times 10^3$, $3.2 \times 10^8$ for Figs. 8(a) and 8(b), respectively. These two displacement fields correspond to the

<figure><img src="image_10.png" /><figcaption>FIG. 8. Displacement field at early times. Displacement fields at (a) $T = T_{\text{MCT}} = 0.12$ and (b) $T = 0.075$ corresponding to the images in Figs. 5(b) and 5(c). Arrows show the particle displacement vectors $\delta\mathbf{r}_i(t) = \mathbf{r}_i(t) - \mathbf{r}_i(0)$ at (a) $t = 2 \times 10^3$, (b) $3.2 \times 10^8$. The linear size of the box is $L = 100$. Right: enlargement of the squared region; displacement vectors are colored by their magnitude.</figcaption></figure>

images in Figs. 5(b) and 5(c). The vector magnitudes are color coded from dark to bright with a maximum magnitude clearly depending on the temperature. We observe a clear spatial heterogeneity in particle displacements, with a coexistence of mobile particles characterized by large displacement vectors and quiescent particles with very small displacements. The displacement field in the quiescent regions qualitatively resembles the structure of the long-wavelength modes among the eigenstates of the Hessian matrix of amorphous glassy states [58–60].

Yet, a careful look reveals that the physics is qualitatively different at the two temperatures. These differences are most appreciated on the right panels, where we magnify some mobile regions. Near $T_{\text{MCT}}$, the mobile cluster contains particles which have moved at most $|\delta\mathbf{r}_i(t)| \approx 1$. The displacements are coherent, extending over distances up to ten particles. This explains the partial decorrelation of the structure, as shown in Fig. 5(b). These displacements resemble stringlike motion [61].

In comparison, mobile clusters close to $T_g$ are composed of particles that have moved significantly larger distances and have fully decorrelated the initial structure. We still observe modelike displacements in the unrelaxed regions, although with a much smaller amplitude. The mobile
particles are clearly identified, and they form clusters of various sizes and geometries, which are quite compact with sharp boundaries. The inset reveals that mobile particles have undergone a large number of rearrangements leading to displacement vectors that are very entangled.

At very low temperatures, the physical picture is that of highly localized rearrangements taking place in an otherwise elastic matrix responding as an amorphous solid. It would be interesting to measure spatial correlations of these small displacements for several temperatures, as they could indicate a possible crossover from delocalized modelike motion at high temperature to elasticlike displacements at lower temperatures.

Inside the clusters of mobile particles, the particle displacements appear complicated and do not form simple patterns such as strings [61,62], loops, or simple swaps. Anticipating Sec. VI, we understand that these complex displacement patterns result, in fact, from the accumulation of a large number of elementary relaxation events which take place inside these localized clusters.

### C. Cluster analysis: Statistical properties

One way to quantify the statistical properties of the localized clusters where structural relaxation occurs at early times is to perform a cluster analysis. We first threshold particle mobility into two families and then consider that two mobile particles belong to the same cluster if their relative distance is smaller than 1.5, corresponding to the first minimum in the total pair distribution function. With these definitions, we can group relaxed particles into independent clusters and perform an analysis of the statistical properties of the clusters observed in the simulations. A similar analysis has been performed before in the regime above $T_{MCT}$ [61,63–65]. Back then, the procedure required relatively arbitrary thresholding in particle mobility, which instead becomes more physical at the low temperatures investigated here. We implement several thresholding procedures, which all yield the same qualitative picture.

In Fig. 9(a), we show the average number of distinct clusters $n_c(t)$ as a function of time for different temperatures in 3D. This quantity is extensive at short times, but its time and temperature evolution should not depend on $N$. At all temperatures, $n_c$ increases with time before eventually reaching a maximum and decreasing at longer times. We compute the pair distribution function of clusters (not reported), which does not show any peculiar spatial dependence at short times. This suggests that the clusters are, at least initially, randomly distributed in space and uncorrelated. The time dependence of $n_c(t)$ is interesting as the growth toward its maximum is well described by a power law which extends over a broader range as $T$ decreases. The corresponding exponent slowly decreases with decreasing $T$ from $n_c(t) \sim t^{0.5}$ at $T_{MCT}$ to $n_c(t) \sim t^{0.38}$ close to $T_g$.

<figure><img src="image_11.png" /><figcaption>FIG. 9. Statistical analysis of relaxed clusters in 3D. (a) Average number of clusters $n_c(t)$ of mobile particles for temperatures as in Fig. 1(a). Dashed lines indicate the time $t_c$ at which $n_c$ is maximum. At low temperature, the number of clusters grows as a power law (line) with exponent 0.38. (b) Average number of mobile particles $n(t)$ per cluster, with a slow power-law growth indicated as a guide to the eye. Inset: average number $n_c(t)$ and size $n(t)$ of clusters at $t = t_c$ as a function of the inverse temperature.</figcaption></figure>

We display in Fig. 9(b) the average number of particles $n(t)$ per mobile cluster. At early times, when the number of clusters starts to grow, their average size $n$ remains modest (about two to three particles), but the cluster size is distributed (see below). This analysis confirms the intuition obtained from Fig. 4: At early times, the dynamics is very heterogeneous in space and time, and structural relaxation is initiated in randomly distributed localized clusters containing only a few particles, which relax in an extended “sea” of blue, quiescent particles. As time increases, the typical size of these relaxed clusters increases, and mobile regions slowly coarsen with time. At the lowest temperature studied, the growth of $n(t)$ at times much shorter than $\tau_\alpha$ is very slow, as revealed by the power law $n(t) \sim t^{0.08}$ indicated in Fig. 9(b).

At later times, the number $n_c$ of mobile clusters reaches a maximum at a time $t = t_c$ before decreasing toward 1 at very long times, when all particles are mobile and the system forms a single mobile cluster [62]. Therefore, we can generically expect the emergence of a maximum number of clusters. At $t_c$, independent clusters start to merge either because two clusters appear at nearby locations or because two nearby growing clusters merge. Either way, the emergence and growth of mobile clusters should naturally induce a form of mobility percolation as time increases [62].
The timescale $t_c$ where $n_c$ is maximal grows with decreasing temperature. From $T = 0.097$ to 0.0715, the ratio $t_c/\tau_\alpha$ decreases from 0.3 to 0.02, while $\tau_\alpha$ grows by about 4 orders of magnitude. Thus, $t_c$ is mostly driven by the evolution of $\tau_\alpha$ but both timescales are not precisely proportional. The time $t_c$ results from the competition between the emergence of new relaxed clusters and the merging of preexisting ones. As such, we see no reason to interpret $t_c$ as an important timescale to describe the structural relaxation, in contrast to Ref. [62]. In particular, $t_c$ does not lead to any observable signature in time correlation functions or relaxation spectra.

Our interpretation for the physical origin of a maximum $n_c$ is supported by the steep evolution observed at $t > t_c$ for the average cluster size $n(t)$. Finally the distribution of cluster sizes $P(n,t)$ at $t = t_c$ is also consistent with the physics of percolation, as it shows a power-law decay with a temperature-independent exponent $P(n,t_c) \sim n^{-1.8}$ (data not shown). Similar distributions have been obtained before [61,64]. The absence of temperature evolution suggests that the strong thresholding performed at relatively high temperatures seriously weakens the physical evolution observed in Fig. 4 and may artificially render the high-temperature regime more heterogeneous and clustered than it really is.

Two characteristic quantities exhibit an interesting temperature evolution, as shown in the inset of Fig. 9. As the temperature decreases, the average number of clusters $n_c(t_c)$ decreases, while their typical size $n(t_c)$ increases. This echoes the observation made from Figs. 5(d)–5(f) close to the structural relaxation time: Relaxed domains become larger and more compact as the temperature decreases. The compactification of the domains at the moment they percolate is confirmed by the temperature evolution of the product $n(t_c) \times n_c(t_c)$ which increases by about 50% (from 80 to 120) in the temperature regime shown in Fig. 9.

## D. Waiting-time distribution of cluster relaxation: Emergence of a power-law tail

The power-law growth of the number of relaxed clusters $n_c(t)$ suggests that the appearance of a new cluster of mobile particles is a stochastic process. Here we define $\tau$ as the first time $t$ at which a new cluster of mobile particles appears. We measure the distribution $\Pi(\log \tau)$ of waiting times $\tau$ for new clusters of relaxed particles in both 2D and 3D; see Fig. 10. To avoid statistical noise, we represent only data for temporal bins which have accumulated more than six events among more than $5 \times 10^3$ for each temperature. At the lowest temperatures, it is no longer possible to normalize the distributions, as we cannot observe the entire range of waiting times for the appearance of relaxed clusters. In the log-log representation of Fig. 10, the low-temperature distributions are defined up to an arbitrary vertical shift. Finally, to compare the distributions at

<figure><img src="image_12.png" /><figcaption>FIG. 10. Emergence of a power-law tail in the waiting-time distribution of relaxed clusters. Waiting-time distribution $\Pi(\log \tau)$ for the appearance of new clusters of relaxed particles in (a) 3D and (b) 2D. Close to $T_g$, the distribution is well fitted at short times $\tau/\tau_\alpha \ll 1$ by a power law with an exponent $b \approx 0.38$ (in 3D) and 0.45 (in 2D). Legend as in Fig. 1.</figcaption></figure>

different temperatures, we normalize the time axis by the average relaxation time $\tau_\alpha(T)$ at each temperature.

For temperatures near $T_{MCT}$ and slightly below, we measure the entire distribution of waiting times $\Pi(\log \tau)$. Strikingly, even at these relatively high temperatures (from an experimental viewpoint), the distributions are already quite broad since they extend over more than four decades of waiting times. Compared to an average relaxation time $\tau_\alpha \approx 10^4$ around $T_{MCT}$, this implies that the width of the distribution is comparable to its average. In other words, there is dynamics happening at all timescales between the microscopic time $t \approx 1 \approx \tau_o$ and the structural relaxation time $\tau_\alpha$.

Exploring the temperature regime below $T_{MCT}$ toward $T_g$, we find that this trend becomes even more prominent. The waiting-time distributions broaden as $T$ decreases. Interestingly, the functional form of the distribution changes when the temperature gets closer to $T_g$ in the sense that it develops a clear power-law tail at $\tau \ll \tau_\alpha$,

$$ \Pi(\log \tau) \sim \tau^b. \qquad (9) $$

The exponent $b(T)$ depends very weakly on the temperatures, as it decreases slowly with decreasing $T$. In the vicinity of $T_g$, we measure $b \approx 0.38$ (in 3D) and $b \approx 0.45$
(in 2D). This power-law behavior, which emerges only at $T \ll T_{MCT}$ is a novel feature revealed by our numerical strategy which probes equilibrium dynamics near $T_g$ for the first time.

The similar power-law behavior found for $\Pi(\log \tau)$ and for $n_c(t)$ in 3D should not come as a surprise, as one generally expects that

$$n_c(t) = \int_{-\infty}^{\log t} d \log \tau \Pi(\log \tau). \quad (10)$$

Plugging Eq. (9) into the above equation and performing a change of variable directly yields $n_c(t) \sim t^b$, as indeed observed in Fig. 9(a).

### E. Short-time dynamics in the frequency domain: Emergence of excess wings

In experiments, the relaxation dynamics of supercooled liquids is often probed thanks to spectroscopy techniques (dielectric, mechanical, light scattering, NMR) which measure linear response functions in the frequency domain [66–68]. To compare our results with experimental measurements, we define the analog of a susceptibility spectrum in the frequency domain. We assume that the relaxation dynamics stems from a distribution of timescales $G(\log \tau)$, such that a correlation function $C(t)$ which monitors structural relaxation can be written as [69]

$$C(t) = \int_{-\infty}^{+\infty} G(\log \tau) e^{-t/\tau} d \log \tau. \quad (11)$$

This expression amounts to decomposing the broad spectrum of relaxation times characterizing supercooled liquids into a series of elementary exponential processes. In the Fourier domain, this yields

$$\chi(\omega) = - \int_{-\infty}^{+\infty} G(\log \tau) \frac{1}{1 + i\omega\tau} d \log \tau. \quad (12)$$

We define the relaxation spectrum $\chi''(\omega)$ as the imaginary part of the susceptibility $\chi(\omega)$. Following earlier work, we express the distribution of timescales $G$ from the measured time correlation function through the relation $G \approx -dC/d \log t$ [70]. This allows us to compute numerically an analog of the out-of-phase susceptibility spectrum $\chi''(\omega)$ [71]. In practice, we choose $C(t)$ to be the self-intermediate scattering function in 3D and the bond-breaking correlation in 2D. While the quantitative detail of the spectra may change for other choices of correlation functions (mainly the relative amplitude of the various processes), the main features are quite robust. This reflects analogous findings of observable dependence in experimental studies [44,66–68].

We report the relaxation spectra in Fig. 11 for $d = 2, 3$. The spectra exhibit a peak at high frequency $\omega \approx 1$ reflecting the rapid decay of the time correlation functions

<figure><img src="image_13.png" /><figcaption>FIG. 11. Emergence of excess wings in relaxation spectra near $T_g$. Equilibrium relaxation spectra $\chi''$ as a function of the frequency $\omega$ in (a) 3D and (b) 2D. The black lines are power-law fits $\chi''(\omega) \sim \omega^{-\sigma}$ with $\sigma \approx 0.38$ (in 3D) and $\sigma \approx 0.45$ (in 2D). Legend as in Fig. 1.</figcaption></figure>

on a microscopic timescale corresponding to thermal motion within a well-defined cage. The lower amplitude of the microscopic peak in 2D compared to 3D is easily explained. The definition of $C_B$, with a larger cutoff employed to define neighbors at later times, indeed makes it less sensitive than $F_s$ to such motion.

The second peak found at lower frequency corresponds to the $\alpha$-relaxation peak and is typically located at frequency $\omega_\alpha \approx 1/\tau_\alpha$. At the lowest temperatures, the $\alpha$ peak exits the numerically accessible frequency window, and we are simply left with the high-frequency flank of the structural relaxation peak.

In both 2D and 3D, the spectra display a power-law behavior for the lowest temperatures near $T_g$,

$$\chi''(\omega) \sim \omega^{-\sigma}, \quad (13)$$

where $\sigma(T)$ is an exponent that appears to depend very weakly on the temperature and slowly decreases as $T$ decreases. The power law in Eq. (13) emerges in the range $\omega \in [10^{-6}, 10^{-3}]$, with an exponent $\sigma \approx 0.38$ in 3D and $\sigma \approx 0.45$ in 2D. The exponents $\sigma$ appearing in the spectra at high frequencies are much smaller than the stretching exponents $\beta = 0.56$ (in 3D) and 0.6 (in 2D) obtained from fitting the corresponding time correlation functions in the $\alpha$-relaxation regime. Therefore, the power law revealed by
Eq. (13) cannot be understood as the high-frequency limit of the $\alpha$ peak. In Ref. [41], we demonstrated that the power laws correspond to the excess wings appearing in the spectra of many molecular liquids. We emphasize that the excess wings become clearly visible only at temperatures $T \gtrsim T_g$, explaining why they have not been reported before in numerical simulations. We also notice that they do not take the form of a “$\beta$ peak” separated from the $\alpha$ peak, as argued recently in Ref. [62], and do not require an internal degree of freedom.

Assuming that the distribution of relaxation times $G(\log \tau)$ in Eq. (12) is dominated by the appearance of new clusters of mobile particles at short times, the power law $\Pi(\log \tau) \sim \tau^b$ would translate into $G(\log \tau) \sim \tau^b$. Plugging this expression into Eq. (12) and performing a change of variable leads to $\chi''(\omega) \sim \omega^{-b}$. This suggests that the exponent $\sigma$ describing $\chi''(\omega)$ is equal to the exponent $b$ describing the excess wing. For the two systems studied here, we indeed get a good agreement with $b \approx \sigma \approx 0.38$ in 3D and $b \approx \sigma \approx 0.45$ in 2D. The quantitative agreement between the measured exponent of the excess wings and the waiting-time distribution of the relaxed cluster shows that the latter process provides a microscopic explanation for this characteristic spectroscopic signature in supercooled liquids near $T_g$.

## F. Suppression of excess wings in ultrastable glasses

To characterize the physical processes responsible for excess wings (or, more generally, of secondary processes) in molecular liquids, relaxation spectra have been measured in nonequilibrium situations [72,73]. For instance, aging studies report a slow decrease of the amplitude of secondary processes with the waiting time $t_w$ after a rapid quench below $T_g$. This suggests that secondary processes become weaker as the glass is annealed toward more stable states. This effect was directly evidenced in recent experiments using physical vapor deposition as a way to efficiently synthesize glasses with profoundly different degrees of stability [73]. This would correspond to varying the waiting time in aging experiments, albeit over inaccessible long timescales. The study unequivocally confirmed the suppression of secondary processes in extremely stable glasses.

To assess whether a similar behavior is observed numerically, we investigate the dynamics of glasses with varying stability at temperature $T = 0.062$ in 3D. We vary glass stability by first preparing equilibrium liquids at three temperatures $T_i = 0.07, 0.08$, and $0.10$ using the swap Monte Carlo algorithm. Using conventional molecular dynamics with a Nosé-Hoover thermostat, we then quench them suddenly to $T = 0.062$. The resulting glasses are aged at $T = 0.062$ during $t_w = 2 \times 10^6$. We check that our measurements do not sensitively depend on $t_w$: It is large enough to ensure that all important structural aging

<figure><img src="image_14.png" /><figcaption>FIG. 12. Suppression of excess wings in ultrastable glasses in 3D. (a) Spectra $\chi''(\omega)$ measured at $T = 0.062$ in glasses first equilibrated at $T_i$, then quenched at $T = 0.062$ where they age during $t_w = 2 \times 10^6$. We report the equilibrium spectrum labeled by $T_i = 0.062$. (b) Average number of clusters $n_c(t)$ of mobile particles with the same definition as in Fig. 9.</figcaption></figure>

processes have shifted outside the observation time window. The temperature $T_i$, which may be interpreted as a fictive temperature [74], mainly controls the degree of stability of the resulting glasses. A similar strategy was adopted in Refs. [23,25,26]. Using the relaxation time $\tau_\alpha(T_i)$ as a quantitative measure of stability, the range of temperatures $T_i$ translates into a variation in glass preparation times between $10^3$ and $10^{12}$, which represents a significant dynamic range.

In Fig. 12(a), we report the relaxation spectra measured from $F_s(t)$ with conventional MD simulations starting from the samples aged during $t_w$. All spectra exhibit a similar frequency dependence with a power law following Eq. (13) and an exponent $\sigma$ that appears almost independent of the glass stability. Still, the glass stability strongly affects the amplitude of the relaxation so that

$$ \chi''(\omega) \approx A(T_i)\omega^{-\sigma} \quad (14) $$

with a prefactor $A(T_i)$ which decreases by more than 1 order of magnitude from $T_i = 0.1$ to $0.062$. Our measurements directly reveal that secondary processes are suppressed in glasses of increasing stability, in excellent agreement with experimental observations [73] where a suppression factor of 3 has been achieved.
<figure><img src="image_15.png" /><figcaption>FIG. 13. Relaxed clusters at equal times in poorly aged glasses and equilibrium liquids in 2D. Snapshots of two systems simulated for $t = 10^7$ at temperature $T = 0.075$. Initial conditions are (a) sample equilibrated at $T_i = T_{MCT} = 0.12$, then aged for $t_w = 4 \times 10^7$ at $T = 0.075$, and (b) equilibrium sample at $T_i = 0.075$. Particles are colored following the bond-breaking correlation, as in Fig. 4. The linear size of the box is $L = 100$.</figcaption></figure>

Our simulations allow us to uncover the physical origin of this suppression. Following the cluster analysis of Sec. V, we present in Fig. 12(b) the time dependence of the number $n_c(t)$ of clusters of mobile particles in aged glasses. All curves display the same time dependence, essentially a power law, with a prefactor that changes by 1 order of magnitude depending on $T_i$. This observation extends the link established above between the relaxed clusters observed at short times and the excess wing observed in relaxation spectra to nonequilibrium glasses.

Finally, we illustrate this point thanks to simulation snapshots drawn from 2D simulations. We repeat a protocol similar to 3D, starting from an equilibrium configuration at $T_i = 0.12 = T_{MCT}$, which is suddenly quenched to $T = 0.075$ where it ages for a long waiting time $t_w = 4 \times 10^7$ using conventional canonical molecular dynamics. We then run the dynamics over a duration $t = 10^7$ and record particle motion. We compare the results to the equilibrium dynamics obtained over the same duration and temperature $T = 0.075$. In Fig. 13, we compare the corresponding mobility fields measured in the poorly annealed system ($T_i = 0.12$) and the very stable one ($T_i = 0.075$). These images directly confirm that the number and size of the relaxed clusters in the two systems are very different. The depletion of short-time excitations as stability increases explains the suppression of excess wings and of secondary relaxations observed in ultrastable glasses. This finding adds a new item to the growing list of glassy excitations that get depleted as glass stability is varied [22–26].

The overall physical conclusion is that the population of localized mobile clusters observed at the early times during the equilibrium relaxation of supercooled liquids becomes stability dependent when studying the nonequilibrium dynamics of glasses. Many more relaxation events are observed in less stable glasses (a factor of 10 in Fig. 12). The nature of these additional events is different from

equilibrium ones, as they necessarily correspond to irreversible events which slowly drive the system toward equilibrium. Our measurements for the amplitude $A(T_i)$ imply that 90% of the secondary processes at play in poorly stable glasses prepared using conventional numerical methods are actually not present in equilibrium conditions [62,75] and great care must be paid regarding the physical interpretation of such studies.

# VI. HOW STRUCTURAL RELAXATION UNFOLDS FROM EARLY TO LATE TIMES

## A. Maps of local relaxation time

We first illustrate how structural relaxation takes place in space and time starting from the early mobile clusters characterized above. We generate a spatial map of local relaxation times for the long trajectory $t_{max} = 8.2 \times 10^8$ at $T = 0.09$ of Fig. 4. Such a representation of the spatially heterogeneous dynamics was used long ago [76–78]. To this end, we attach to each particle its relaxation time $\tau_\alpha^i$ defined as the first time at which particle $i$ becomes mobile, keeping our definition of mobility from the inequality $C_B^i(\tau_\alpha^i) \le 0.5$.

We show in Fig. 14 two representative maps of the local relaxation time obtained at $T = 0.12 \approx T_{MCT}$ [Fig. 14(a)] and a temperature $T = 0.09$ halfway between $T_{MCT}$ and $T_g$ [Fig. 14(b)]. The latter corresponds to the long trajectory of Fig. 4 which lasts $t_{max} = 8 \times 10^8$. The gray regions in Fig. 14(b) contain the particles that have not yet relaxed after $t_{max}$. The color code uses a logarithmic scale from bright (yellow) at short times to dark (black) at long times. As found in Fig. 10, there are particles which become mobile at extremely short times. Note that the shortest $\tau_\alpha^i$ are bounded by the time delay $\Delta t$ used to store configurations along the trajectory [we use $\Delta t = 10^3$ and $10^5$ for

<figure><img src="image_16.png" /><figcaption>FIG. 14. Maps of local relaxation time reveal the emergence of dynamic facilitation. Map of a local relaxation time $\tau_\alpha^i$ in the two-dimensional system defined as the first time at which $C_B^i \le 0.5$ and (a) $T = 0.12 = T_{MCT}$ and (b) $T = 0.09$. The color codes for the value of $\tau_\alpha^i$ are on a logarithmic scale. In (b), the arrow highlights the emergence of dynamic facilitation, and gray regions have not yet relaxed at $t_{max} = 8 \times 10^8$. The box size is $L = 100$.</figcaption></figure>
Figs. 14(a) and 14(b), respectively]. The upper bound is set by the longest local relaxation time in Fig. 14(a) (about $20\tau_\alpha$) or by the simulation time $t_{max}$ in Fig. 14(b).

The images in Fig. 14 illustrate the very broad waiting-time distributions for the appearance of new mobile particles in the sample. We observe, in particular, many zones where particles become mobile orders of magnitude before the average structural relaxation time $\tau_\alpha$. They appear as bright yellow zones, and coexist with regions that instead have a relaxation time larger than $\tau_\alpha$. We identify many more fast zones near $T_{MCT}$ than at $T = 0.09$. The slow coarsening of relaxed regions observed in Fig. 4 translates into a smooth spatial spreading of the color in these relaxation-time maps. Mobility appearing in a localized region at early time extends at longer times to its neighboring regions which thus appear darker. We also observe that the typical size of the correlated regions seems to increase on a logarithmic scale: Bright yellow domains are very small while darker domains are larger. These images thus confirm that the relaxation starting at localized regions at early times slowly spreads to neighboring regions over logarithmically increasing times.

This very clear spatiotemporal view of the structural relaxation becomes actually much more complicated at higher temperatures. Around the mode-coupling temperature in Fig. 14(a), the relaxation starts at early times (bright) at numerous locations that are rather fuzzy. Similarly, the spread of mobility to neighboring regions is much harder to identify. This observation suggests that dynamic facilitation is much more clearly defined close to $T_g$ than in the higher-temperature regime explored by conventional simulations.

We also remark that mobility does not spread isotropically. On the bottom left of Fig. 14(b), one can indeed distinguish a domain which relaxes very early, surrounded in one direction by a facilitated region identified by a smooth color gradient to the south (indicated by an arrow), while to the north, one distinguishes a gray region in which no relaxation takes place over the subsequent four decades. Therefore, dynamic facilitation acts isotropically on average, but it appears to be locally anisotropic. This suggests that the local disorder plays an important role and controls how mobility can propagate in space. This observation seems in harmony with the anisotropic kinetic constraints introduced in certain kinetically constrained models [70,79,80].

In three-dimensional liquids, we measure the local relaxation time and display how the dynamics evolves in space and time. To ease visualization, we show in Fig. 15 a selected relaxed cluster of particles at very low temperature $T = 0.07$. Following Fig. 14, the color code indicates the local relaxation time. As for 2D systems at low temperature, mobility starts at very early times in small localized regions and spreads to the neighboring particles over time to form a larger cluster of mobile particles at larger times.

<figure><img src="image_17.png" /><figcaption>FIG. 15. Dynamic facilitation in the 3D model. A connected cluster containing 163 mobile particles at $T = 0.07$ inside a 3D system composed of $N = 10^4$ particles (linear size $L = 21.5$) at time $t = 4.1 \times 10^6$. The particles are represented at their position at $t = 0$ using a color which codes for $\tau_\alpha^i$.</figcaption></figure>

Qualitatively, dynamic facilitation thus appears to also play an important role in 3D.

## B. Space-time trajectories

In constructing maps of the local relaxation time, each particle is shown only once with an indication about the first time it becomes mobile. These maps say nothing about what happens after this first relaxation event. To address this point, we construct a $d+1$ space-time trajectory representing the evolution of the mobility field along a time axis for $d=2$ liquids [8]. To define mobility on a given time interval of duration $\Delta t$, we use the usual criterion $C_B^i(t + \Delta t, t) < 0.5$. By stacking the mobility field over consecutive time slices of duration $\Delta t$, we can visualize how mobility propagates in space and time.

In Fig. 16, we show a representative space-time mobility trajectory at $T = 0.09$ obtained from the same data as in Figs. 4 and 14(b). The trajectory of duration $t_{max} = 8 \times 10^8$ is split into 400 time slices of width $\Delta t = 2 \times 10^6$. The choice of $\Delta t$ results from a compromise. It must be small enough compared to $\tau_\alpha$ so that we can resolve how the structural relaxation unfolds from early times, but not too small so the fraction of mobile particles within each time interval is nonvanishing. With our choice, about 5% of the particles are mobile in each slice. Only particles that are mobile with a time slice are rendered. To ease visualization, we show only particles with $0 < x < L/2$. We highlight the main features of the space-time trajectory by constructing a surface representation of the set of mobile particles, as implemented in the OVITO software [81].

The most striking novel observation deduced from Fig. 16 is the emergence of “tubes” in space-time. Physically, this corresponds to localized regions in space which keep relaxing in many successive time slices. We conclude that a large number of relaxation events accumulates in some regions, even long after the first relaxation event. In Sec. VI C, we quantify this accumulation effect by introducing an appropriate statistical tool.
<figure><img src="image_18.png" /><figcaption>FIG. 16. Space-time trajectory in the 2D model. A trajectory starting at $t=0$ (left) finishing at $t_{max} = 8 \times 10^8$ (right) at $T=0.09$ is split into 400 time slices of duration $\Delta t = 2 \times 10^6$. We show all particles that are mobile within a given time slice, together with a surface representation which draws a contour of the coarse-grained mobility field. To ease $d+1$ visualization, the color linearly codes for the time dimension, and one half of the system in the $x$ direction is removed.</figcaption></figure>

In addition, a careful examination of Fig. 16 reveals that these tubes are not perfectly aligned along the time direction. They can bend, widen, or merge with neighboring tubes. In terms of particle motion, this means that the relaxation events which accumulate in a localized region are not simply the repetition of the exact same event (for instance, a complex but reversible particle motion), which means that relaxation can spread to neighboring particles. As a result, mobility propagates in space from one time slice to another. This observation accounts for the slow coarsening of the relaxed regions observed in Fig. 4. In Sec. VI D below, we quantify this coarsening process. The propagation of mobility from one region to a neighboring one via the accumulation of complex localized relaxation events provides the microscopic origin of dynamic facilitation.

### C. Accumulation of localized relaxation events

To quantify the accumulation of mobility in the same regions over time, we denote as $m_1(\Delta t)$ the fraction of particles that are mobile (defined again via $C_B^i$) in the time slice $[t, t + \Delta t]$. We then identify the particles that are mobile over the next time slice $[t + \Delta t, t + 2\Delta t]$. Finally, we define $m_2(\Delta t)$ as the fraction of particles that are mobile in two consecutive time intervals, i.e., mobile in both $[t, t + \Delta t]$ and $[t + \Delta t, t + 2\Delta t]$.

By construction, $m_2 \le m_1$, the equality being obtained when the dynamics repeats itself exactly from one time

slice to the next. In the opposite case of completely uncorrelated dynamics between two consecutive intervals, one gets $m_2 = m_1^2 < m_1$. To quantify the mobility correlation between consecutive time slices, we introduce the quantity $c(\Delta t)$,

$$c(\Delta t) = \frac{m_2 - m_1^2}{m_1 - m_1^2}, \qquad (15)$$

whose normalization ensures that $0 \le c(\Delta t) \le 1$. By definition, we have $c(\Delta t) = 1$ if the dynamics is perfectly correlated over two consecutive intervals, while $c(\Delta t) = 0$ if it is completely uncorrelated. For a given $\Delta t$, a large value $c(\Delta t)$ implies that a large fraction of mobile particles are identical in successive time intervals of duration $\Delta t$. This is illustrated in Fig. 17(a).

We investigate the evolution of $c(\Delta t)$ as a function of $\Delta t$ and temperature $T$ in both 2D and 3D; see Fig. 17. To

<figure><img src="image_19.png" /><figcaption>FIG. 17. The accumulation of relaxation events becomes more probable toward $T_g$. (a) The quantity $c(\Delta t)$ in Eq. (15) quantifies the probability for particles to be mobile within two consecutive intervals of duration $\Delta t$. Data for $c(\Delta t)$ in 3D [(b), legend provided in Fig. 1(a)] and 2D (c) liquids at various temperatures. The time axis is rescaled by $\tau_\alpha$ measured from $F_s(t)$ and $C_B(t)$, respectively.</figcaption></figure>
compare different temperatures, we rescale $\Delta t$ by the
average relaxation time $\tau_\alpha$. The latter is measured from
$F_s(t)$ in 3D and $C_B(t)$ in 2D. All curves exhibit a
qualitatively similar time dependence with $c(\Delta t) \approx 0$ at
both microscopic times $\Delta t \approx 1$ and at very long times
$\Delta t \gg \tau_\alpha$, where the dynamics in successive frames is
obviously uncorrelated. The probability $c(\Delta t)$ thus exhibits
a maximum at a time which scales almost like $\tau_\alpha$. The value
at the maximum increases with decreasing the temperature.
We also see that at any fixed value of $\Delta t/\tau_\alpha$, the value of
$c(\Delta t)$ grows when $T$ decreases.

These measurements quantitatively confirm the tube
interpretation of the dynamics detected in the space-time
trajectory of Fig. 16. At low temperature, the mobility field
measured in a given time frame becomes increasingly
correlated with the mobility field measured in the con-
secutive time frame: Relaxation events accumulate at
identical locations over long times to form tubes in
d + 1 dimensions.

The existence of a growing maximum at $\Delta t \approx \tau_\alpha$ shows that the spatial structure of dynamic heterogeneities characterizing the structural relaxation becomes increasingly similar from one slice of duration $\tau_\alpha$ to the next as the temperature decreases. For example, at $T = 0.075$ in 3D, we find that $c(\tau_\alpha) \approx 0.5$, which implies that half of the particles that relax over one relaxation time $\tau_\alpha$ also relax over the next relaxation time. One can physically anticipate that such a mechanism leads to an increased lifetime of dynamic heterogeneities at low temperature.

D. Slow temperature-dependent coarsening
of relaxed domains

We now investigate how mobility spreads in space
over larger timescales. Natural quantities to measure the
growth of spatially correlated regions are four-point
spatial correlation functions of the mobility field, either
in real or in Fourier space [82–84]. Both functions are
related to the four-point dynamic susceptibility that has
been the subject of a large number of studies [85–88]. It is,
however, well known that collecting good statistics for
these functions requires very large systems [89–92]. This
represents a numerical effort on its own, which we leave
for future work.

We characterize growing dynamic length scales using
an alternative method based on the chord length distri-
bution. This method was introduced to analyze porous and
bicontinuous structures [93–95]. We note that the method
does not faithfully characterize dynamic heterogeneities at
temperatures above T<sub>MCT</sub>, where relaxed domains have
fuzzy shapes and ill-defined boundaries. It was, however,
shown to efficiently determine characteristic length
scales in compact bicontinuous structures [94], which
are observed at T<sub>g</sub> < T < T<sub>MCT</sub>. Finally, this method does
not suffer from considerations related to the choice of
statistical ensembles [90,91].

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

FIG. 18. Slow coarsening of relaxed domains measured by chord length distributions in 2D. (a) Snapshots of a portion of the two-dimensional liquid with a mobility field at $T = 0.095$ (color code from Fig. 4). (b) Discretized version of the mobility field used to measure the length $\ell$ of chords drawn vertically and horizontally which intersect relaxed domains (white). (c) Chord length distribution $P(\ell, t)$ at $T = 0.095$ and increasing times. (d) Growth of the average chord length $\langle \ell \rangle$ with time. A slow algebraic growth develops at early times and low temperatures.

First, we identify relaxed regions via the bond-breaking
correlation; see Fig. 18(a). To determine the spatial exten-
sion of relaxed domains, it is convenient to perform a
coarse-graining and binarize the $C_B^i$ field. To do so, we
discretize the configuration into a grid of cells of linear
size 1 and compute their average $C_B^i$, particles being
weighted by their area overlap with each cell, which is
then thresholded at 0.5 between mobile and quiescent cells.
We illustrate in Fig. 18 how our procedure transforms a
particle configuration [Fig. 18(a)] into a discrete lattice of
mobile and quiescent cells [Fig. 18(b)]. The binary grid is
then used to measure the distribution of chord length $\ell$
defined as horizontal and vertical segments intersecting
mobile domains, as shown in Fig. 18(b).

After averaging, the method yields the probability
distribution of chord lengths $P(\ell, t)$ at any temperature $T \le$
$T_{\text{MCT}}$ and time $t$. As an example, we show in Fig. 18(c) the
chord length distribution measured at $T = 0.095$ in 2D for
three different times. All distributions have a maximum at
very small $\ell \approx 1$, presumably due to the roughness of the
domains. At larger $\ell$, the distribution is smooth and decays over a typical distance which grows with increasing time. The decay at large $\ell$ is roughly exponential, as seen in other systems [94].

To extract a time-dependent length scale characterizing the growing size of the domains of relaxed particles, we compute the first moment of the chord length distribution,

$$ \langle \ell \rangle = \int_{0}^{\infty} d\ell P(\ell, t)\ell. \qquad (16) $$

The time dependence of $\langle \ell \rangle$ is shown in Fig. 18(d) for a broad range of temperatures.

The time dependence of the typical chord length qualitatively changes with the temperature. Close to $T_{MCT} = 0.12$ and above, the average relaxation time $\tau_{\alpha}$ is not yet very large. Many small domains appear spontaneously across the system at early times, and they rapidly merge with one another. The merging of independent clusters is the mechanism that controls the growth of the typical domain size in the time regime near $\tau_{\alpha}$. In other words, dynamic facilitation is present above $T_{MCT}$, but it does not seem to directly control the time dependence of spatially correlated dynamic domains.

At lower temperature instead, we observe a sparse population of regions that appear at early times and coarsen independently over a large time window before different domains start to merge. In this early-time regime, which exists only at temperatures much lower than $T_{MCT}$, the typical domain size appears to grow as a power law

$$ \langle \ell \rangle (t) \sim t^{1/z(T)} \qquad (17) $$

with a dynamic exponent $z(T)$ which increases slowly as $T$ decreases. For $T = 0.0853$, we find, for instance, $z(T) \approx 5.9$ (i.e., $1/z \approx 0.17$). This very large value of the dynamic exponent $z(T)$ implies that the growth of relaxed regions is strongly subdiffusive in the time regime $t \ll \tau_{\alpha}$. The decrease of this exponent with temperature demonstrates that, while the dynamics repeats itself with increasing probability, it also takes longer for mobility to spread to neighboring regions. As a word of caution, we note that our numerical determination of $z(T)$ is performed over ranges of timescales and temperatures which are highly constrained, leaving open the possibility that the domain growth is actually better described by a logarithmic increase, for instance. Finally, the very slow growth of $\langle \ell \rangle$ in the 2D model echoes the similarly slow growth of the average cluster size $n(t)$ measured in the 3D model in Fig. 9.

## E. How the story ends: Lifetime of dynamic heterogeneity

We finally describe how structural relaxation takes place at very long times, beyond the structural relaxation time. Typically, $\tau_{\alpha}$ is defined from the decay of a time correlation function $C(t)$ below an arbitrary level, for instance,

<figure><img src="image_21.png" /><figcaption>FIG. 19. Relaxation of very slow domains in the 2D model. $x-y$ projection of the space-time mobility trajectory shown in Fig. 16. The color code indicates, in a linear scale, the time at which particles become mobile for the first time.</figcaption></figure>

$C(\tau_{\alpha}) = e^{-1}$. When $C(t)$ is the self-intermediate function, this means that particles have, on average, moved by a distance $2\pi/q$ after $\tau_{\alpha}$. Of course, a similar value is reached if a fraction of the particles has moved a lot, while the rest is still immobile. Using the bond-breaking correlation function at half value, as we frequently do above, means that either all particles have, on average, lost half of their neighbors, or that one half of the particles have lost all of their neighbors while the others did nothing. In other words, the definition of an “average” relaxation time $\tau_{\alpha}$ does not imply that all particles have completely relaxed after a time $\tau_{\alpha}$. In fact, the snapshots shown above clearly demonstrate that, at very low temperatures, a significant fraction of the particles have not relaxed at all after $\tau_{\alpha}$ while the others have fully relaxed. This paragraph merely reformulates more explicitly the well-established statement that the dynamics of supercooled liquids is spatially heterogeneous [96].

We push the simulations at $T = 0.09$ up to 5 times longer than the $\alpha$-relaxation time and represent the outcome in Fig. 19. We show the $x-y$ projection of the space-time representation of Fig. 16, the time axis now pointing inside the plane, with a linear color code from short (white) to long (dark) times. When accumulated over such a large timescale, we observe that a significant fraction of the system remains completely unrelaxed (gray domains). We also observe that the latest stages of the dynamics where some slow regions relax for the first time (red regions) spread from regions that have already relaxed. We make similar observations in multiple trajectories: The slowest regions to relax are always invaded and relaxed from their boundaries. In other words, the regions that have already relaxed at $\tau_{\alpha}$ are progressively
invading the ones that relax slower than average. This
suggests that the tail of the correlation function at very
long times is physically controlled by the slow coarsening
process described in the previous section.

This has two important consequences.

First, this implies that the gray unrelaxed regions in
Fig. 19 are extremely stable and would tend to relax over a
timescale that is much larger than the average relaxation
time τα. However, their environment is not as stable, and
the relaxation which takes place there facilitates or accel-
erates the relaxation of these slow domains. In Ref. [97], we
argued that this phenomenon provides a physical explan-
ation for the asymmetric shape of relaxation spectra in
deeply supercooled liquids in terms of an underlying
distribution of “natural” relaxation times that has a very
broad tail that is cut off by the facilitated relaxation of the
slowest regions. This mechanism has been discussed in
various contexts [98–102].

Second, the timescale over which the slowest regions
relax in the system controls the so-called lifetime τ<sub>DH</sub> of the
dynamic heterogeneity. Physically, a bicontinuous pattern
of the mobility field with fast and slow particles emerges
near τ<sub>α</sub>, as shown, for instance, in Fig. 5. If the relaxation
time in the slow regions in such plots were given by a
timescale τ<sub>DH</sub> ≫ τ<sub>α</sub>, then the same mobility pattern would
be found over many consecutive time intervals of duration
τ<sub>α</sub> and would start to change significantly only after a time
τ<sub>DH</sub> [103]. The study of this lifetime has been the subject of
intense experimental [96,104–107] and numerical [43,108–
110] studies.

Our simulations show that, as a result of dynamic
facilitation, the lifetime of the slow regions is controlled
by a combination of two factors: (i) the relaxation time of the
fast regions and (ii) the coarsening of mobility from one
relaxed region to the next. Because the latter factor depends
very weakly on the temperature, we expect that τ<sub>DH</sub> may
increase slightly more slowly than τ<sub>α</sub> but should not be
strongly decoupled from it. In simulations where complete
decorrelation can be observed, we find that after 20–40 times
τ<sub>α</sub>, the entire system has completely relaxed and information
about dynamic heterogeneity is then totally lost. This
conclusion appears consistent with recent experiments
[106,107] and simulations [43]. Such a long timescale also
justifies our statement about equilibration in supercooled
liquids, which cannot be reached before at least 20–40 τ<sub>α</sub> as
otherwise, some particles would not have relaxed at all.

VII. DISCUSSION

A. Summary of main novel results

Let us first briefly recapitulate the main new results
arising from studying dynamics over long timescales at
very low temperatures.

First, we observe that in the very low-temperature regime
the dynamics starts at sparse locations and takes the form of
localized but quite complex relaxation events. Detailed

analysis reveals the emergence of power laws characteriz-
ing the time dependence of three quantities: the number of
independent clusters of mobile particles, the distribution of
waiting time for their appearance, and the high-frequency
dependence of ensemble-averaged relaxation spectra
(excess wing). We demonstrate that these power laws are
intimately connected, and their amplitude strongly depends
on the stability of the system.

Second, we describe how these early relaxation
events induce an accumulation of relaxation events repeat-
edly taking place at roughly the same location leading to
the slow coarsening of the relaxed domains. These obser-
vations account for the emergence, at a coarse-grained level,
of dynamic facilitation. In particular, we characterize a
dynamic exponent z(T) relating timescales and length
scales, which we extract by introducing a chord length
analysis of the dynamic heterogeneity. The corresponding
growth of relaxed regions is very slow and strongly sub-
diffusive with a large value of z(T), possibly logarithmic.

Third, we observe a qualitative evolution of the nature of
the spatially heterogeneous dynamics at all timescales as
the temperature is decreased much below the mode-
coupling crossover. In particular, we find an increasing
segregation between mobile and immobile regions result-
ing in spatial fluctuations of the mobility that become more
compact with smoother boundaries at low temperatures.

Fourth, we find that at low temperature the regions that
relax the slowest at timescales much larger than τα seem so
stable that the fastest mechanism to relax them is via the
propagation of mobility from faster regions surrounding
them, that is, via dynamic facilitation. Our simulations thus
demonstrate that the lifetime of the dynamic heterogeneity
is controlled by dynamic facilitation.

B. Comparison with previous simulations

Computational research regarding the dynamics of
supercooled liquids is immense and covers several decades
of work [11]. Because we explore a temperature regime that
was not accessible to simulations before, it is useful to
contrast our main results with earlier studies. This is
organized in three broad topics.

First, we discuss the emergence of excess wings
associated with short-time localized relaxation events
characterized by a broad waiting-time distribution, as
reported in Sec. V. This finding more broadly refers to
the topic of secondary relaxations, or β processes, in the
dynamics of supercooled liquids [2,111]. We recall that
excess wings appear in experiments in a time window
roughly between 1 µs and 1 s for temperatures much
lower than T<sub>MCT</sub> [112–114]. Simulations performed with-
out the swap Monte Carlo algorithm cannot reach equi-
librium at temperatures where excess wings appear. As a
consequence, earlier attempts to explore similar time
and temperature windows necessarily dealt with non-
equilibrium glasses obtained by crossing a computer glass
transition at a temperature much higher than $T_g$, typically near $T_{MCT}$. In Sec. V F, we show that relaxation events detected in such poorly annealed glasses are not representative of those observed in equilibrium materials at the same temperature. In a recent series of simulations using computer models not very different from ours [62,75], these events were mathematically described as an additive secondary peak in relaxation spectra. Other works have reported similar conclusions regarding the existence of a $\beta$ process [115]. Our results show that this description may not be adequate, as the high-frequency shoulder observed in these works in fact turns into extended power laws in equilibrium at lower temperatures.

Our work also sheds light on an important issue regarding secondary relaxations. A number of numerical studies have reported the existence of a $\beta$ process taking the form of a peak in the frequency domain located at a frequency fully decoupled from $\tau_\alpha$. These studies analyzed the relaxation behavior of particles that are more complex than the point particles studied here, for instance, polymeric systems [116,117] or particles with shapes [118–120]. These studies suggest that intramolecular degrees of freedom can make relaxation spectra even more complex than the ones we report. However, our observation that excess wings taking the form of extended power laws exist in our much simpler model shows that secondary processes do exist in the absence of intramolecular degrees of freedom.

A second major theme in computer studies is the idea that a crossover temperature separates two physical regimes in the dynamics of supercooled liquids [35]. This is theoretically rationalized by the existence of a critical temperature controlling the dynamics of supercooled liquids in the context of mean-field approaches [50,121,122]. This critical temperature is similar to the one predicted from mode-coupling theory, although its interpretation as an artifact of mean-field approximations has only recently been fully clarified [123,124]. It is generally agreed that such mean-field approaches may usefully describe the first few decades of the dynamic slowdown, which is precisely the temperature regime studied over the years by computer simulations [35]. Having access to considerably lower temperatures confirms that relaxation dynamics keeps changing with decreasing temperatures. However, the present results are in fact more naturally described as a progressive evolution toward low-temperature physics, with no sharp change at a specific temperature. Importantly, most of the features that we report become prominent and unambiguous only at the lowest temperatures that we can analyze. We suggest that a considerable amount of work is needed to reassess conclusions drawn in previous works at higher temperatures.

Third, our results show how dynamic facilitation emerges at very low temperature and plays an important role in both the approach to $\tau_\alpha$ and at longer times. The

concept of dynamic facilitation has a long history in the glass literature [98,125], to the point that it serves as a basis for the construction of a large family of kinetic glass models [126], which can display many phenomena observed in more realistic particle systems and have been, in particular, very useful to characterize dynamic heterogeneity [127].

Direct signatures of dynamic facilitation in molecular dynamics simulations have been discussed before [128–132]. The most direct studies have, in particular, followed similar principles: First, use some thresholding procedure to detect the time at which a given particle is moving (i.e., performing some jump), and then search for the enhanced probability that a nearby particle will relax not too far in the future [129]. Our philosophy is conceptually different. We do not perform single-particle analysis of jump dynamics but rather quantify how the observed collective relaxation events effectively give rise to dynamic facilitation when observed over some large timescales and length scales, as in the space-time representation of Fig. 16 or via the emergence of a dynamic exponent $z(T)$. Our simulations suggest that dynamic facilitation becomes more prominent and should be more easily characterized at lower temperatures, which seems to contradict earlier findings [133].

Finally, a recent numerical study suggests that dynamic facilitation may become long-ranged below $T_{MCT}$ and may result from elastic interactions [134]. We do not repeat the detailed analysis proposed in this work at the much lower temperatures analyzed here, but the gradual coarsening of relaxed domains suggests that facilitation does not act over very large distances and can be quantitatively described, at a coarse-grained level, by a dynamic exponent similar to the one emerging from nearest-neighbor facilitation in kinetically constrained models. Reconciling these two observations represents a worthwhile topic for future work.

### C. Consequences for glass transition physics

What do we learn about the physics of supercooled liquids, and how does this knowledge impact existing theories of the glass transition?

An important lesson drawn from our results is that equilibrium relaxation dynamics at low temperatures appears qualitatively different from observations made in the usual temperature regime covered by molecular dynamics studies. Whereas particle motion results from collective modelike mobility patterns involving most particles at high $T$, we observe instead localized particle motion that completely relaxes the structure of mobile particles leaving the rest of the system essentially unaffected at low $T$. These observations capture, in real space, the commonly accepted idea that dynamics smoothly transitions from a nonactivated relaxation regime that efficiently relaxes the system at high temperature to a regime dominated by activated processes involving a finite number of particles at low
<figure><img src="image_22.png" /><figcaption>FIG. 20. Convergence of relaxation times at low temperature in 3D (a) and 2D (b). Relaxation times $\tau_\alpha$ defined from the self-intermediate scattering function $F_s(t)$ (a) and bond-orientational correlation $C_\psi(t)$ (b) compared to the relaxation time $\tau_\alpha^B$ defined from the bond-breaking correlation $C_B(t)$. The very different timescales measured at high temperatures become identical at low temperatures, much below $T_{MCT}$ (dotted dashed line).</figcaption></figure>

temperature embedded in an otherwise frozen elastic matrix.

In practice, this crossover happens quite gradually as the relaxation time increases by about 7 orders of magnitude, and it is difficult to claim that a particular temperature scale, such as the mode-coupling crossover $T_{MCT}$, marks a sharp change in the physics. Previous studies offered examples of both sharp [135–137] and smooth [37,39,40,138,139] changes between the two types of physics, and our analysis of the relaxation dynamics is in better agreement with the latter family of studies.

To corroborate this point, we present in Fig. 20 the temperature evolution of the average relaxation times deduced from two types of correlation functions. Previous simulations performed in the nonactivated high-temperature regime have reported a strong decoupling [48,54,140,141] between the bond-breaking correlation $C_B(t)$ and more conventional time correlation functions, such as the self-intermediate scattering function $F_s(t)$ or the bond-orientational correlation $C_\psi(t)$. This has been rationalized by the presence of vibrational modelike mobility patterns such as the one shown in Fig. 8, that are able to decorrelate $F_s(t)$ and $C_\psi(t)$, but would leave $C_B(t)$ essentially unaffected [54]. The results shown in Fig. 20 confirm this finding at high temperatures. However, pushing the analysis to much lower temperatures reveals that these two timescales eventually converge. At the lowest temperature, localized and complex activated relaxation processes result in large displacements that are able to simultaneously decorrelate $F_s(t)$ or $C_\psi(t)$ and $C_B(t)$. This qualitative change in the dynamics toward localized activated relaxation events emerges gradually, and the two relaxation times start to coincide at temperatures that are in fact much lower than $T_{MCT}$. In practice, this implies that the dynamic heterogeneity at these low temperatures becomes

fairly independent of the particular observable and specific choices made to detect particle mobility.

The idea of a crossover toward activated dynamics at low temperatures has a long history and is, in particular, at the core of the random first-order-transition (RFOT) theory [50,142,143]. The theoretical explanation is that the dynamic transition predicted in the mean-field limit must be avoided in finite dimensions, as the lifetime of glassy states necessarily becomes finite [122]. A crossover between nonactivated high-temperature dynamics and activated localized events also emerges naturally in kinetically constrained models [144], but this occurs without invoking an avoided dynamic transition. Therefore, the observation of a crossover, in itself, does not confirm any particular theoretical approach.

More interestingly perhaps, the direct observation of localized activated events reported here suggests that one should now be able to understand and characterize their nature better. Our preliminary investigations show that relaxation dynamics at low temperatures involves an extremely large number of transitions between inherent structures, for instance, implying that the potential energy landscape cannot be used to gather information about relaxation dynamics as already noticed [145]. It remains numerically challenging to unambiguously group inherent structures into larger metabasins [138,146–149] or to quantify free-energy barriers [150], as envisioned within RFOT theory. Revisiting these earlier attempts is an important goal to characterize the nature of activated dynamics in the regime where it is effectively present.

We repeatedly argue that dynamic facilitation emerges in the relaxation dynamics at low temperatures. This observation is a crucial new piece of evidence in the long-lasting debate between dynamic and static explanations of the glass transition [6]. Clearly, the thermodynamic RFOT theory picture of a liquid broken into a mosaic of droplets undergoing collective activated relaxation events [142,151] is invalidated by our results, at least in the simplest picture of independently relaxing droplets.

There have been several attempts to introduce some degree of dynamic facilitation within RFOT theory, with the argument that fast relaxing droplets could affect the dynamics of slower ones in their neighborhood, as described phenomenologically in Refs. [99,102], or more formally in Refs. [101,152]. Our results show that such a combination would be needed to account for the dynamics in laboratory experiments, and quantitatively explain the temperature evolution of central physical quantities, such as the $\alpha$-relaxation time or the length scale of dynamic heterogeneity [153].

In Fig. 21, we provide a decisive quantitative support for this conclusion in 2D. We directly compare the temperature evolution of the static point-to-set length $\xi_{PTS}$ [154] that sets the typical size of the droplets in the mosaic picture [142,151], and the typical length scale of the dynamic
<figure><img src="image_23.png" /><figcaption>FIG. 21. Decoupling of static and dynamic length scales in 2D supercooled liquids. The dynamic length scale is the average chord length $\langle l \rangle$ of dynamic heterogeneities at $C_B = 0.5$. The point-to-set length scale $\xi_{PTS}$ is equal to $2R_c$, with $R_c$ the crossover radius from Ref. [18]. Mode-coupling crossover temperature (dashed-dotted line) and glass transition temperature (dashed) are indicated.</figcaption></figure>

heterogeneity at the structural relaxation time deduced from the chord length analysis described in Sec. VI D. The data for $\xi_{PTS}$ are taken from Ref. [18] for our 2D model. A strong decoupling between static and dynamic length scales is observed. This has been reported before for $T > T_{MCT}$ [37,155] in a regime where activated dynamics should not be relevant. This had left open the possibility that the static point-to-set length could catch up with the dynamic one in a low-temperature regime dominated by activated processes, as discussed in Ref. [156]. Our results invalidate this hypothesis, and demonstrate that, at least for our 2D system, the characteristic size of the dynamic heterogeneity is not controlled by $\xi_{PTS}$, even at temperatures near the experimental glass transition $T_g$. This is a major conclusion that emerges from the present study: The simplest formulation of RFOT theory does not describe our data well without invoking an important facilitation component [99,101,102].

Does this conclusion invalidate thermodynamic theories of the glass transition? Clearly not. First, future work should consolidate our findings across a broader range of models, and, more importantly, in 3D where the measurement and comparison of various length scales require a massive computational effort. Establishing whether 3D models behave as in Fig. 21 is an important task. Second, it could still be that the emergence of clusters at short times that triggers further relaxation over growing domains is due to the type of activated dynamics envisioned by RFOT theory, and these events could very well be controlled by the growing static length scale $\xi_{PTS}$. As an illustration of this idea, we mark in Fig. 5 the spatial extent of the point-to-set length scale, which compares reasonably well with the spatial extension of the domains which first relax. If correct, this interpretation would imply that structural relaxation emerges from a combination of physical

ingredients, requiring descriptive tools stemming from both dynamic and static approaches [101]. More broadly, future work should reconcile the present findings regarding the dynamics of supercooled liquids to the recent detailed characterization of global [19,20] and local [43,157] thermodynamic fluctuations of the overlap order parameter [158,159] which are at the root of RFOT theory. Third, the accumulation of relaxation events over long periods of time in certain locations suggests that dynamic facilitation necessarily has a static origin [160,161]. In addition, explaining how relaxation at one place triggers relaxation in the neighborhood should also require linking statics to dynamics [132]. Therefore, the observation of dynamic facilitation does not automatically imply that static considerations and structural and thermodynamic information become irrelevant. Future work should explore how to best describe or incorporate dynamic facilitation effects in the context of RFOT theory.

Does the observation of dynamic facilitation automatically validate all predictions stemming from analogies with kinetically constrained models [8]? Clearly not, and a lot remains to be done. A key result we obtain in this regard is the direct measurement of a large dynamic exponent $z(T)$. The measured value invalidates entire families of kinetically constrained models, including the simplest version of the Fredrickson-Andersen facilitated model [162], but also the large family of collective kinetically constrained models [125,163–166], where nontrivial defects display diffusive motion [167] corresponding to $z = 2$. Clearly, our results are much closer to the behavior found in the family of anisotropic east [70,79] and arrow [80] models where defects are subdiffusive [127,168,169]. However, further work should establish, for instance, whether the power laws characterizing the short-time dynamics revealed in Sec. V are similar to the ones reported in Ref. [70] in the 3D east model which also give rise to excess wings. It would also be useful to understand whether and how the concept of kinetic constraints emerges from microscopic atomic motion and particle interaction, and whether the exponent $z(T)$ varies quantitatively as predicted in the east model.

## D. Further perspectives

Finally we describe four important research directions suggested by our results.

The first important task ahead of us is to test whether the results obtained here can be generalized to other models by changing the type of particle interactions and possibly generalizing to more complicated particle shapes and molecules. This clearly requires improving the swap Monte Carlo algorithm even further [15,16] as well as developing molecular dynamics software and hardware that are even more efficient [170,171]. It would also be useful to vary the spatial dimensions over a broader spectrum [172–174] to see whether glassy dynamics changes qualitatively by increasing $d$ over a larger range.
A second important task concerns a more detailed characterization of the spatially heterogeneous dynamics at the very low temperatures that can now be accessed. Decades of work have led to the development of an arsenal of tools consisting of multipoint correlation functions and dynamic susceptibilities [87,88]. These quantities are well known and their physical content fully understood. Technically, however, there are subtleties and potential artifacts that need to be carefully considered [90]. In addition, a proper measurement of correlation length scales, their time and temperature evolution, as well as the characterization of the geometry of these correlations require significant statistics along with simulation boxes that are large enough [89,92]. Such measurements would allow a quantitative characterization of the growth of characteristic length scales for dynamic heterogeneity as well as the evolution of the geometry of correlated regions for which quantitative scenarios exist [70,175]. This is a clear task for future work, but it necessitates a large investment of computational resources.

A third research effort should aim at a better characterization of dynamic facilitation, which becomes more prominent at low temperatures than it is above $T_{MCT}$. It could be useful to explore the behavior of the various tools introduced before [130,132,149] to correlate single-particle motion in space and time at much lower temperatures. Starting with Ref. [132], several papers have developed computational tools to identify the analog of the excitations that define kinetically constrained models in particle-based simulations [160,176–178]. Would these algorithms be able to detect the sparse population of relaxing clusters that drives the relaxation near $T_g$? Methods should also be developed to more directly and more precisely measure the dynamic exponent $z(T)$ relating timescales and length scales, as this is an important emerging consequence of dynamic facilitation. Finally, we merely observe the emergence of dynamic facilitation, but we do not provide a microscopic understanding of its origin in terms of structure, thermodynamics, or geometry of the supercooled liquid. This ambitious task would be required to provide a fully microscopic picture of structural relaxation near $T_g$.

A fourth ambitious goal is the development of more efficient molecular dynamics schemes that would allow us to fill the remaining gap between the 30 ms studied here and the 100 s timescale that characterize most experiments. Our results show that the physics keeps changing qualitatively when shifting the “glass ceiling” by several orders of magnitude. For the types of models studied here, the remaining gap in timescales is about $10^4$. Such a colossal gap cannot be filled by hardware improvements or by performing longer simulations and will require creative development and invention of molecular dynamics techniques, maybe coupled to smart Monte Carlo moves. Filling this gap appears as the last obstacle before reaching the Holy Grail: simulating 100 s of the life of a supercooled liquid at the experimental glass transition.

## ACKNOWLEDGMENTS

We would like to thank Jean-Philippe Bouchaud for his support, questions, and constant interest over the last two years. Many members of the Simons Collaboration “Cracking the glass problem” have challenged our numerical results and interpretations over the last two years, which helped us improve our manuscript. We thank G. Jung, D. S. Simmons, and P. G. Wolynes for comments on our manuscript. This work is supported by a grant from the Simons Foundation (Grant No. 454933, L. B.) and the European Research Council under the EU’s Horizon 2020 Program, Grant No. 740269. C. S. acknowledges support from the Fondation l’Oréal for a L’Oréal-UNESCO for Women in Science Fellowship, a Herchel Smith Fellowship, University of Cambridge, and a Ramon Jenkins Research Fellowship from Sidney Sussex College, Cambridge. B. G. acknowledges support by Capital Fund Management—Fondation pour la Recherche.

C. Scalliet and B. Guiselin contributed equally to this work.

[1] M. D. Ediger, C. A. Angell, and S. R. Nagel, Supercooled Liquids and Glasses, J. Phys. Chem. **100**, 13200 (1996).

[2] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, and S. W. Martin, Relaxation in Glassforming Liquids and Amorphous Solids, J. Appl. Phys. **88**, 3113 (2000).

[3] C. A. Angell, Formation of Glasses from Liquids and Biopolymers, Science **267**, 1924 (1995).

[4] L. Berthier and G. Biroli, Theoretical Perspective on the Glass Transition and Amorphous Materials, Rev. Mod. Phys. **83**, 587 (2011).

[5] P. G. Debenedetti and F. H. Stillinger, Supercooled Liquids and the Glass Transition, Nature (London) **410**, 259 (2001).

[6] G. Tarjus, An Overview of the Theories of the Glass Transition, in Dynamical Heterogeneities in Glasses, Colloids, and Granular Media (Oxford University Press, New York, 2011), Vol. 150, pp. 39–67.

[7] V. Lubchenko and P. G. Wolynes, Theory of Structural Glasses and Supercooled Liquids, Annu. Rev. Phys. Chem. **58**, 235 (2007).

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

[9] G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot, The Frustration-Based Approach of Supercooled Liquids and the Glass Transition: A Review and Critical Assessment, J. Phys. Condens. Matter **17**, R1143 (2005).

[10] C. P. Royall and S. R. Williams, The Role of Local Structure in Dynamical Arrest, Phys. Rep. **560**, 1 (2015).

[11] J.-L. Barrat and L. Berthier, Computer Simulations of the Glass Transition and Glassy Materials, arXiv:2206.01013.
[12] L. Berthier and D. R. Reichman, *Modern Computational Studies of the Glass Transition*, arXiv:2208.02206.

[13] T. S. Grigera and G. Parisi, *Fast Monte Carlo Algorithm for Supercooled Soft Spheres*, Phys. Rev. E **63**, 045102(R) (2001).

[14] A. Ninarello, L. Berthier, and D. Coslovich, *Models and Algorithms for the Next Generation of Glass Transition Studies*, Phys. Rev. X **7**, 021039 (2017).

[15] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, *Efficient Swap Algorithms for Molecular Dynamics Simulations of Equilibrium Supercooled Liquids*, J. Stat. Mech. (2019) 064004.

[16] A. D. S. Parmar, B. Guiselin, and L. Berthier, *Stable Glassy Configurations of the Kob-Andersen Model Using Swap Monte Carlo*, J. Chem. Phys. **153**, 134505 (2020).

[17] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, *Configurational Entropy Measurements in Extremely Supercooled Liquids That Break the Glass Ceiling*, Proc. Natl. Acad. Sci. U.S.A. **114**, 11356 (2017).

[18] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, *Zero-Temperature Glass Transition in Two Dimensions*, Nat. Commun. **10**, 1508 (2019).

[19] B. Guiselin, L. Berthier, and G. Tarjus, *Random-Field Ising Model Criticality in a Glass-Forming Liquid*, Phys. Rev. E **102**, 042129 (2020).

[20] B. Guiselin, L. Berthier, and G. Tarjus, *Statistical Mechanics of Coupled Supercooled Liquids in Finite Dimensions*, SciPost Phys. **12**, 091 (2022).

[21] Y. Nishikawa, M. Ozawa, A. Ikeda, P. Chaudhuri, and L. Berthier, *Relaxation Dynamics in the Energy Landscape of Glass-Forming Liquids*, Phys. Rev. X **12**, 021001 (2022).

[22] C. Scalliet, L. Berthier, and F. Zamponi, *Absence of Marginal Stability in a Structural Glass*, Phys. Rev. Lett. **119**, 205501 (2017).

[23] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, *Random Critical Point Separates Brittle and Ductile Yielding Transitions in Amorphous Materials*, Proc. Natl. Acad. Sci. U.S.A. **115**, 6656 (2018).

[24] C. Scalliet, L. Berthier, and F. Zamponi, *Nature of Excitations and Defects in Structural Glasses*, Nat. Commun. **10**, 5102 (2019).

[25] D. Khomenko, C. Scalliet, L. Berthier, D. R. Reichman, and F. Zamponi, *Depletion of Two-Level Systems in Ultrastable Computer-Generated Glasses*, Phys. Rev. Lett. **124**, 225901 (2020).

[26] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel, and E. Flenner, *Low-Frequency Vibrational Modes of Stable Glasses*, Nat. Commun. **10**, 26 (2019).

[27] W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, *Glass Stability Changes the Nature of Yielding under Oscillatory Shear*, Phys. Rev. Lett. **124**, 225502 (2020).

[28] L. Wang, L. Berthier, E. Flenner, P. Guan, and G. Szamel, *Sound Attenuation in Stable Glasses*, Soft Matter **15**, 7018 (2019).

[29] M. P. Eastwood, T. Chitra, J. M. Jumper, K. Palmo, A. C. Pan, and D. E. Shaw, *Rotational Relaxation in Ortho-Terphenyl: Using Atomistic Simulations to Bridge Theory and Experiment*, J. Phys. Chem. B **117**, 12898 (2013).

[30] J.-H. Hung, T. K. Patra, V. Meenakshisundaram, J. Hari Mangalara, and D. S. Simmons, *Universal Localization Transition Accompanying Glass Formation: Insights from Efficient Molecular Dynamics Simulations of Diverse Supercooled Liquids*, Soft Matter **15**, 1223 (2019).

[31] L. Berthier and M. D. Ediger, *How to “Measure” a Structural Relaxation Time That Is Too Long to Be Measured?*, J. Chem. Phys. **153**, 044501 (2020).

[32] G. L. Hunter and E. R. Weeks, *The Physics of the Colloidal Glass Transition*, Rep. Prog. Phys. **75**, 066501 (2012).

[33] S. Gokhale, A. K. Sood, and R. Ganapathy, *Deconstructing the Glass Transition through Critical Experiments on Colloids*, Adv. Phys. **65**, 363 (2016).

[34] W. Kob and H. C. Andersen, *Testing Mode-Coupling Theory for a Supercooled Binary Lennard-Jones Mixture I: The van Hove Correlation Function*, Phys. Rev. E **51**, 4626 (1995).

[35] W. Kob, *Computer Simulations of Supercooled Liquids and Glasses*, J. Phys. Condens. Matter **11**, R85 (1999).

[36] G. Brambilla, D. El Masri, M. Pierno, L. Berthier, L. Cipelletti, G. Petekidis, and A. B. Schofield, *Probing the Equilibrium Dynamics of Colloidal Hard Spheres above the Mode-Coupling Glass Transition*, Phys. Rev. Lett. **102**, 085703 (2009).

[37] W. Kob, S. Roldán-Vargas, and L. Berthier, *Non-Monotonic Temperature Evolution of Dynamic Correlations in Glass-Forming Liquids*, Nat. Phys. **8**, 164 (2012).

[38] C. Dalle-Ferrier, C. Thibierge, C. Alba-Simionesco, L. Berthier, G. Biroli, J.-P. Bouchaud, F. Ladieu, D. L'Hôte, and G. Tarjus, *Spatial Correlations in the Dynamics of Glassforming Liquids: Experimental Determination of Their Temperature Dependence*, Phys. Rev. E **76**, 041510 (2007).

[39] P. Das and S. Sastry, *Crossover in Dynamics in the Kob-Andersen Binary Mixture Glass-Forming Liquid*, J. Non-Cryst. Solids **14**, 100098 (2022).

[40] D. Coslovich, M. Ozawa, and W. Kob, *Dynamic and Thermodynamic Crossover Scenarios in the Kob-Andersen Mixture: Insights from Multi-CPU and Multi-GPU Simulations*, Eur. Phys. J. E **41**, 62 (2018).

[41] B. Guiselin, C. Scalliet, and L. Berthier, *Microscopic Origin of Excess Wings in Relaxation Spectra of Supercooled Liquids*, Nat. Phys. **18**, 468 (2022).

[42] M. Ozawa, C. Scalliet, A. Ninarello, and L. Berthier, *Does the Adam-Gibbs Relation Hold in Simulated Supercooled Liquids?*, J. Chem. Phys. **151**, 084504 (2019).

[43] L. Berthier, *Self-Induced Heterogeneity in Deeply Supercooled Liquids*, Phys. Rev. Lett. **127**, 088002 (2021).

[44] B. Schmidtke, N. Petzold, R. Kahlau, M. Hofmann, and E. A. Rössler, *From Boiling Point to Glass Transition Temperature: Transport Coefficients in Molecular Liquids Follow Three-Parameter Scaling*, Phys. Rev. E **86**, 041507 (2012).

[45] B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret, and P. Keim, *Mermin-Wagner Fluctuations in 2D Amorphous Solids*, Proc. Natl. Acad. Sci. U.S.A. **114**, 1856 (2017).

[46] S. Vivek, C. P. Kelleher, P. M. Chaikin, and E. R. Weeks, *Long-Wavelength Fluctuations and the Theory and Experiment*, J. Phys. Chem. B **117**, 12898 (2013).
Glass Transition in Two Dimensions and Three Dimensions, Proc. Natl. Acad. Sci. U.S.A. **114**, 1850 (2017).

[47] E. Flenner and G. Szamel, Fundamental Differences between Glassy Dynamics in Two and Three Dimensions, Nat. Commun. **6**, 7392 (2015).

[48] Z. Flenner and G. Szamel, Viscoelastic Shear Stress Relaxation in Two-Dimensional Glass-Forming Liquids, Proc. Natl. Acad. Sci. U.S.A. **116**, 2015 (2019).

[49] J. A. van Meel, L. Filion, C. Valeriani, and D. Frenkel, A Parameter-Free, Solid-Angle Based, Nearest-Neighbor Algorithm, J. Chem. Phys. **136**, 234107 (2012).

[50] W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory (Oxford University Press, New York, 2009), Vol. 143.

[51] T. B. Schrøder, S. Sastry, J. C. Dyre, and S. C. Glotzer, Crossover to Potential Energy Landscape Dominated Dynamics in a Model Glass-Forming Liquid, J. Chem. Phys. **112**, 9834 (2000).

[52] T. B. Schrøder and J. C. Dyre, Solid-like Mean-Square Displacement in Glass-Forming Liquids, J. Chem. Phys. **152**, 141101 (2020).

[53] R. Yamamoto and A. Onuki, Heterogeneous Diffusion in Highly Supercooled Liquids, Phys. Rev. Lett. **81**, 4915 (1998).

[54] H. Shiba, T. Kawasaki, and A. Onuki, Relationship between Bond-Breakage Correlations and Four-Point Correlations in Heterogeneous Glassy Dynamics: Configuration Changes and Vibration Modes, Phys. Rev. E **86**, 041504 (2012).

[55] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, Localized Soft Modes and the Supercooled Liquid's Irreversible Passage through Its Configuration Space, J. Chem. Phys. **131**, 194508 (2009).

[56] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.12.041028 for videos corresponding to Figs. 4 and 5.

[57] P. Chaudhuri, L. Berthier, and W. Kob, Universal Nature of Particle Displacements Close to Glass and Jamming Transitions, Phys. Rev. Lett. **99**, 060604 (2007).

[58] A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, Irreversible Reorganization in a Supercooled Liquid Originates from Localized Soft Modes, Nat. Phys. **4**, 711 (2008).

[59] G. M. Hocky and D. R. Reichman, A Small Subset of Normal Modes Mimics the Properties of Dynamical Heterogeneity in a Model Supercooled Liquid, J. Chem. Phys. **138**, 12A537 (2013).

[60] C. Brito and M. Wyart, Geometric Interpretation of Previtrification in Hard Sphere Liquids, J. Chem. Phys. **131**, 024504 (2009).

[61] C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Stringlike Cooperative Motion in a Supercooled Liquid, Phys. Rev. Lett. **80**, 2338 (1998).

[62] H.-B. Yu, R. Richert, and K. Samwer, Structural Rearrangements Governing Johari-Goldstein Relaxations in Metallic Glasses, Sci. Adv. **3**, e1701577 (2017).

[63] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Dynamical Heterogeneities in a Supercooled Lennard-Jones Liquid, Phys. Rev. Lett. **79**, 2827 (1997).

[64] F. W. Starr, J. F. Douglas, and S. Sastry, The Relationship of Dynamical Heterogeneity to the Adam-Gibbs and Random First-Order Transition Theories of Glass Formation, J. Chem. Phys. **138**, 12A541 (2013).

[65] P. Chaudhuri, S. Sastry, and W. Kob, Tracking Heterogeneous Dynamics during the α Relaxation of a Simple Glass Former, Phys. Rev. Lett. **101**, 190601 (2008).

[66] P. Lunkenheimer, U. Schneider, R. Brand, and A. Loid, Glassy Dynamics, Contemp. Phys. **41**, 15 (2000).

[67] T. Körber, R. Stäglich, C. Gainaru, R. Böhmer, and E. A. Rössler, Systematic Differences in the Relaxation Stretching of Polar Molecular Liquids Probed by Dielectric vs Magnetic Resonance and Photon Correlation Spectroscopy, J. Chem. Phys. **153**, 124510 (2020).

[68] M. Flämiq, M. Hofmann, N. Fatkullin, and E. A. Rössler, NMR Relaxometry: The Canonical Case Glycerol, J. Phys. Chem. B **124**, 1557 (2020).

[69] R. Böhmer, R. V. Chamberlin, G. Diezemann, B. Geil, A. Heuer, G. Hinze, S. C. Kuebler, R. Richert, B. Schiener, H. Sillescu, H. W. Spiess, U. Tracht, and M. Wilhelm, Nature of the Non-Exponential Primary Relaxation in Structural Glass-Formers Probed by Dynamically Selective Experiments, J. Non-Cryst. Solids **235**, 1 (1998).

[70] L. Berthier and J. P. Garrahan, Numerical Study of a Fragile Three-Dimensional Kinetically Constrained Model, J. Phys. Chem. B **109**, 3578 (2005).

[71] T. Blochowicz, C. Tschirwitz, S. Benkhof, and E. A. Rössler, Susceptibility Functions for Slow Relaxation Processes in Supercooled Liquids and the Search for Universal Relaxation Patterns, J. Chem. Phys. **118**, 7544 (2003).

[72] U. Schneider, R. Brand, P. Lunkenheimer, and A. Loidl, Excess Wing in the Dielectric Loss of Glass Formers: A Johari-Goldstein β Relaxation?, Phys. Rev. Lett. **84**, 5560 (2000).

[73] H.-B. Yu, M. Tylinski, A. Guiseppi-Elie, M. D. Ediger, and R. Richert, Suppression of β Relaxation in Vapor-Deposited Ultrastable Glasses, Phys. Rev. Lett. **115**, 185501 (2015).

[74] A. Q. Tool, Relation between Inelastic Deformability and Thermal Expansion of Glass in Its Annealing Range, J. Am. Ceram. Soc. **29**, 240 (1946).

[75] H.-B. Yu, M.-H. Yang, Y. Sun, F. Zhang, J.-B. Liu, C.-Z. Wang, K.-M. Ho, R. Richert, and K. Samwer, Fundamental Link between β Relaxation, Excess Wings, and Cage-Breaking in Metallic Glasses, J. Phys. Chem. Lett. **9**, 5877 (2018).

[76] M. Foley and P. Harrowell, The Spatial Distribution of Relaxation Times in a Model Glass, J. Chem. Phys. **98**, 5069 (1993).

[77] P. Harrowell, Visualizing the Collective Motions Responsible for the α and β Relaxations in a Model Glass, Phys. Rev. E **48**, 4359 (1993).

[78] M. M. Hurley and P. Harrowell, Kinetic Structure of a Two-Dimensional Liquid, Phys. Rev. E **52**, 1694 (1995).

[79] J. Jäckle and S. Eisinger, A Hierarchically Constrained Kinetic Ising Model, Z. Phys. B **84**, 115 (1991).
[80] J. P. Garrahan and D. Chandler, *Coarse-Grained Microscopic Model of Glass Formers*, Proc. Natl. Acad. Sci. U.S.A. **100**, 9710 (2003).

[81] A. Stukowski, *Visualization and Analysis of Atomistic Simulation Data with OVITO–The Open Visualization Tool*, Model. Simul. Mater. Sci. Eng. **18**, 015012 (2010).

[82] N. Lačević, F. W. Starr, T. B. Schröder, and S. C. Glotzer, *Spatially Heterogeneous Dynamics Investigated via a Time-Dependent Four-Point Density Correlation Function*, J. Chem. Phys. **119**, 7372 (2003).

[83] C. Bennemann, C. Donati, J. Baschnagel, and S. C. Glotzer, *Growing Range of Correlated Motion in a Polymer Melt on Cooling towards the Glass Transition*, Nature (London) **399**, 246 (1999).

[84] C. Donati, S. C. Glotzer, and P. H. Poole, *Growing Spatial Correlations of Particle Displacements in a Simulated Liquid on Cooling toward the Glass Transition*, Phys. Rev. Lett. **82**, 5064 (1999).

[85] S. C. Glotzer, V. N. Novikov, and T. B. Schröder, *Time-Dependent, Four-Point Density Correlation Function Description of Dynamical Heterogeneity and Decoupling in Supercooled Liquids*, J. Chem. Phys. **112**, 509 (2000).

[86] C. Donati, S. Franz, S. C. Glotzer, and G. Parisi, *Theory of Non-Linear Susceptibility and Correlation Length in Glasses and Liquids*, J. Non-Cryst. Solids **307**, 215 (2002).

[87] C. Toninelli, M. Wyart, L. Berthier, G. Biroli, and J.-P. Bouchaud, *Dynamical Susceptibility of Glass Formers: Contrasting the Predictions of Theoretical Scenarios*, Phys. Rev. E **71**, 041505 (2005).

[88] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos, *Dynamical Heterogeneities in Glasses, Colloids, and Granular Media* (Oxford University Press, New York, 2011), Vol. 150.

[89] S. Karmakar, C. Dasgupta, and S. Sastry, *Analysis of Dynamic Heterogeneity in a Glass Former from the Spatial Correlations of Mobility*, Phys. Rev. Lett. **105**, 015701 (2010).

[90] L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob, K. Miyazaki, and D. R. Reichman, *Spontaneous and Induced Dynamic Fluctuations in Glass Formers. I. General Results and Dependence on Ensemble and Dynamics*, J. Chem. Phys. **126**, 184503 (2007).

[91] L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob, K. Miyazaki, and D. R. Reichman, *Spontaneous and Induced Dynamic Correlations in Glass Formers. II. Model Calculations and Comparison to Numerical Simulations*, J. Chem. Phys. **126**, 184504 (2007).

[92] E. Flenner and G. Szamel, *Dynamic Heterogeneities above and below the Mode-Coupling Temperature: Evidence of a Dynamic Crossover*, J. Chem. Phys. **138**, 12A523 (2013).

[93] P. Levitz, *Off-Lattice Reconstruction of Porous Media: Critical Evaluation, Geometrical Confinement and Molecular Transport*, Adv. Colloid Interface Sci. **76**, 71 (1998).

[94] V. Testard, L. Berthier, and W. Kob, *Intermittent Dynamics and Logarithmic Domain Growth during the Spinodal Decomposition of a Glass-Forming Liquid*, J. Chem. Phys. **140**, 164502 (2014).

[95] V. Testard, L. Berthier, and W. Kob, *Influence of the Glass Transition on the Liquid-Gas Spinodal Decomposition*, Phys. Rev. Lett. **106**, 125702 (2011).

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

[97] C. Scalliet, B. Guiselin, and L. Berthier, *Excess Wings and Asymmetric Relaxation Spectra in a Facilitated Trap Model*, J. Chem. Phys. **155**, 064505 (2021).

[98] R. G. Palmer, D. L. Stein, E. Abrahams, and P. W. Anderson, *Models of Hierarchically Constrained Dynamics for Glassy Relaxation*, Phys. Rev. Lett. **53**, 958 (1984).

[99] X. Xia and P. G. Wolynes, *Microscopic Theory of Heterogeneity and Nonexponential Relaxations in Supercooled Liquids*, Phys. Rev. Lett. **86**, 5526 (2001).

[100] C. Rehwald and A. Heuer, *How Coupled Elementary Units Determine the Dynamics of Macroscopic Glass-Forming Systems*, Phys. Rev. E **86**, 051504 (2012).

[101] S. M. Bhattacharyya, B. Bagchi, and P. G. Wolynes, *Facilitation, Complexity Growth, Mode Coupling, and Activated Dynamics in Supercooled Liquids*, Proc. Natl. Acad. Sci. U.S.A. **105**, 16077 (2008).

[102] L. Berthier, G. Biroli, J.-P. Bouchaud, and G. Tarjus, *Can the Glass Transition Be Explained without a Growing Static Length Scale?*, J. Chem. Phys. **150**, 094501 (2019).

[103] A. Heuer, *Information Content of Multitime Correlation Functions for the Interpretation of Structural Relaxation in Glass-Forming Systems*, Phys. Rev. E **56**, 730 (1997).

[104] K. Schmidt-Rohr and H. W. Spiess, *Nature of Nonexponential Loss of Correlation above the Glass Transition Investigated by Multidimensional NMR*, Phys. Rev. Lett. **66**, 3020 (1991).

[105] C.-Y. Wang and M. D. Ediger, *How Long Do Regions of Different Dynamics Persist in Supercooled o-Terphenyl?*, J. Phys. Chem. B **103**, 4177 (1999).

[106] K. Paeng, H. Park, D. T. Hoang, and L. J. Kaufman, *Ideal Probe Single-Molecule Experiments Reveal the Intrinsic Dynamic Heterogeneity of a Supercooled Liquid*, Proc. Natl. Acad. Sci. U.S.A. **112**, 4952 (2015).

[107] K. Paeng and L. J. Kaufman, *Single Molecule Experiments Reveal the Dynamic Heterogeneity and Exchange Time Scales of Polystyrene near the Glass Transition*, Macromolecules **49**, 2876 (2016).

[108] S. Léonard and L. Berthier, *Lifetime of Dynamic Heterogeneity in Strong and Fragile Kinetically Constrained Spin Models*, J. Phys. Condens. Matter **17**, S3571 (2005).

[109] K. Kim and S. Saito, *Multi-Time Density Correlation Functions in Glass-Forming Liquids: Probing Dynamical Heterogeneity and Its Lifetime*, J. Chem. Phys. **133**, 044511 (2010).

[110] K. Kim and S. Saito, *Multiple Length and Time Scales of Dynamic Heterogeneities in Model Glass-Forming Liquids: A Systematic Analysis of Multi-Point and Multi-Time Correlations*, J. Chem. Phys. **138**, 12A506 (2013).

[111] J. D. Stevenson and P. G. Wolynes, *A Universal Origin for Secondary Relaxations in Supercooled Liquids and Structural Glasses*, Nat. Phys. **6**, 62 (2010).

[112] P. K. Dixon, L. Wu, S. R. Nagel, B. D. Williams, and J. P. Carini, *Scaling in the Relaxation of Supercooled Liquids*, Phys. Rev. Lett. **65**, 1108 (1990).
[113] N. Menon and S. R. Nagel, Evidence for a Divergent Susceptibility at the Glass Transition, Phys. Rev. Lett. **74**, 1230 (1995).

[114] A. Kudlik, C. Tschirwitz, T. Blochowicz, S. Benkhof, and E. Rössler, Slow Secondary Relaxation in Simple Glass Formers, J. Non-Cryst. Solids **235**, 406 (1998).

[115] H. Zhang, X. Wang, H.-B. Yu, and J. F. Douglas, Fast Dynamics in a Model Metallic Glass-Forming Material, J. Chem. Phys. **154**, 084505 (2021).

[116] D. Bedrov and G. D. Smith, Molecular Dynamics Simulation Study of the α- and β-Relaxation Processes in a Realistic Model Polymer, Phys. Rev. E **71**, 050801(R) (2005).

[117] Dmitry Bedrov and Grant D Smith, Secondary Johari-Goldstein Relaxation in Linear Polymer Melts Represented by a Simple Bead-Necklace Model, J. Non-Cryst. Solids **357**, 258 (2011).

[118] D. Fragiadakis and C. M. Roland, Molecular Dynamics Simulation of the Johari-Goldstein Relaxation in a Molecular Liquid, Phys. Rev. E **86**, 020501(R) (2012).

[119] D. Fragiadakis and C. M. Roland, Cooperativity of Short-Time Dynamics Revisited, J. Non-Cryst. Solids **500**, 18 (2018).

[120] K. Shiraishi, H. Mizuno, and A. Ikeda, Johari-Goldstein β Relaxation in Glassy Dynamics Originates from Two-Scale Energy Landscape, arXiv:2209.01726.

[121] T. R. Kirkpatrick and D. Thirumalai, p-Spin-Interaction Spin-Glass Models: Connections with the Structural Glass Problem, Phys. Rev. B **36**, 5388 (1987).

[122] G. Parisi, P. Urbani, and F. Zamponi, Theory of Simple Glasses: Exact Solutions in Infinite Dimensions (Cambridge University Press, Cambridge, England, 2020).

[123] A. Ikeda and K. Miyazaki, Mode-Coupling Theory as a Mean-Field Description of the Glass Transition, Phys. Rev. Lett. **104**, 255704 (2010).

[124] T. Maimbourg, J. Kurchan, and F. Zamponi, Solution of the Dynamics of Liquids in the Large-Dimensional Limit, Phys. Rev. Lett. **116**, 015902 (2016).

[125] G. H. Fredrickson and H. C. Andersen, Kinetic Ising Model of the Glass Transition, Phys. Rev. Lett. **53**, 1244 (1984).

[126] F. Ritort and P. Sollich, Glassy Dynamics of Kinetically Constrained Models, Adv. Phys. **52**, 219 (2003).

[127] J. P. Garrahan, P. Sollich, and C. Toninelli, Kinetically Constrained Models, in Dynamical Heterogeneities in Glasses, Colloids, and Granular Media (Oxford University Press, New York, 2011), Vol. 150, pp. 111–137.

[128] Y. Gebremichael, M. Vogel, and S. C. Glotzer, Particle Dynamics and the Development of String-like Motion in a Simulated Monoatomic Supercooled Liquid, J. Chem. Phys. **120**, 4415 (2004).

[129] M. Vogel and S. C. Glotzer, Spatially Heterogeneous Dynamics and Dynamic Facilitation in a Model of Viscous Silica, Phys. Rev. Lett. **92**, 255901 (2004).

[130] M. N. J. Bergroth, M. Vogel, and S. C. Glotzer, Examination of Dynamic Facilitation in Molecular Dynamics Simulations of Glass-Forming Liquids, J. Phys. Chem. B **109**, 6748 (2005).

[131] R. Candelier, A. Widmer-Cooper, J. K. Kummerfeld, O. Dauchot, G. Biroli, P. Harrowell, and D. R. Reichman, Spatiotemporal Hierarchy of Relaxation Events, Dynamical Heterogeneities, and Structural Reorganization in a Supercooled Liquid, Phys. Rev. Lett. **105**, 135702 (2010).

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

[133] R. Candelier, O. Dauchot, and G. Biroli, Building Blocks of Dynamical Heterogeneities in Dense Granular Media, Phys. Rev. Lett. **102**, 088001 (2009).

[134] R. N. Chacko, F. P. Landes, G. Biroli, O. Dauchot, A. J. Liu, and D. R. Reichman, Elastoplasticity Mediates Dynamical Heterogeneity below the Mode Coupling Temperature, Phys. Rev. Lett. **127**, 048002 (2021).

[135] K. Broderix, K. K. Bhattacharya, A. Cavagna, A. Zippelius, and I. Giardina, Energy Landscape of a Lennard-Jones Liquid: Statistics of Stationary Points, Phys. Rev. Lett. **85**, 5360 (2000).

[136] T. S. Grigera, A. Cavagna, I. Giardina, and G. Parisi, Geometric Approach to the Dynamic Glass Transition, Phys. Rev. Lett. **88**, 055502 (2002).

[137] D. Coslovich, A. Ninarello, and L. Berthier, A Localization Transition Underlies the Mode-Coupling Crossover of Glasses, SciPost Phys. **7**, 077 (2019).

[138] B. Doliwa and A. Heuer, What Does the Potential Energy Landscape Tell Us about the Dynamics of Supercooled Liquids and Glasses?, Phys. Rev. Lett. **91**, 235501 (2003).

[139] L. Berthier, G. Biroli, D. Coslovich, W. Kob, and C. Toninelli, Finite-Size Effects in the Dynamics of Glass-Forming Liquids, Phys. Rev. E **86**, 031502 (2012).

[140] R. Yamamoto and A. Onuki, Dynamics of Highly Supercooled Liquids: Heterogeneity, Rheology, and Diffusion, Phys. Rev. E **58**, 3515 (1998).

[141] T. Kawasaki and K. Kim, Identifying Time Scales for Violation/Preservation of Stokes-Einstein Relation in Supercooled Water, Sci. Adv. **3**, e1700399 (2017).

[142] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Scaling Concepts for the Dynamics of Viscous Liquids near an Ideal Glassy State, Phys. Rev. A **40**, 1045 (1989).

[143] V. Lubchenko and P. G. Wolynes, Barrier Softening near the Onset of Nonactivated Transport in Supercooled Liquids: Implications for Establishing Detailed Connection between Thermodynamic and Kinetic Anomalies in Supercooled Liquids, J. Chem. Phys. **119**, 9088 (2003).

[144] L. Berthier and J. P. Garrahan, Real Space Origin of Temperature Crossovers in Supercooled Liquids, Phys. Rev. E **68**, 041201 (2003).

[145] A. Heuer, Exploring the Potential Energy Landscape of Glass-Forming Systems: From Inherent Structures via Metabasins to Macroscopic Transport, J. Phys. Condens. Matter **20**, 373101 (2008).

[146] S. Büchner and A. Heuer, Metastable States as a Key to the Dynamics of Supercooled Liquids, Phys. Rev. Lett. **84**, 2168 (2000).

[147] B. Doliwa and A. Heuer, Hopping in a Supercooled Lennard-Jones Liquid: Metabasins, Waiting Time Distribution, and Diffusion, Phys. Rev. E **67**, 030501(R) (2003).
[148] B. Doliwa and A. Heuer, *Energy Barriers and Activated Dynamics in a Supercooled Lennard-Jones Liquid*, Phys. Rev. E **67**, 031506 (2003).

[149] M. Vogel, B. Doliwa, A. Heuer, and S. C. Glotzer, *Particle Rearrangements during Transitions between Local Minima of the Potential Energy Landscape of a Binary Lennard-Jones Liquid*, J. Chem. Phys. **120**, 4404 (2004).

[150] M. Baity-Jesi, G. Biroli, and D. R. Reichman, *Revisiting the Concept of Activation in Supercooled Liquids*, Eur. Phys. J. E **44**, 77 (2021).

[151] J.-P. Bouchaud and G. Biroli, *On the Adam-Gibbs-Kirkpatrick-Thirumalai-Wolynes Scenario for the Viscosity Increase in Glasses*, J. Chem. Phys. **121**, 7347 (2004).

[152] S. M. Bhattacharyya, B. Bagchi, and P. G. Wolynes, *Bridging the Gap between the Mode Coupling and the Random First Order Transition Theories of Structural Relaxation in Liquids*, Phys. Rev. E **72**, 031509 (2005).

[153] G. Biroli and J.-P. Bouchaud, *The RFOT Theory of Glasses: Recent Progress and Open Issues*, arXiv:2208.05866.

[154] A. Montanari and G. Semerjian, *Rigorous Inequalities between Length and Time Scales in Glassy Systems*, J. Stat. Phys. **125**, 23 (2006).

[155] P. Charbonneau and G. Tarjus, *Decorrelation of the Static and Dynamic Length Scales in Hard-Sphere Glass Formers*, Phys. Rev. E **87**, 042305 (2013).

[156] W. Kob, S. Roldan-Vargas, and L. Berthier, *Spatial Correlations in Glass-Forming Liquids across the Mode-Coupling Crossover*, Phys. Procedia **34**, 70 (2012).

[157] B. Guiselin, G. Tarjus, and L. Berthier, *Static Self-Induced Heterogeneity in Glass-Forming Liquids: Overlap as a Microscope*, J. Chem. Phys. **156**, 194503 (2022).

[158] B. Guiselin, G. Tarjus, and L. Berthier, *On the Overlap between Configurations in Glassy Liquids*, J. Chem. Phys. **153**, 224502 (2020).

[159] B. Guiselin, G. Tarjus, and L. Berthier, *Is Glass a State of Matter?*, arXiv:2207.14204.

[160] L. Ortlieb, T. S. Ingebrigtsen, J. E. Hallett, F. Turci, and C. P. Royall, *Relaxation Mechanisms in Supercooled Liquids Past the Mode-Coupling Crossover: Cooperatively Rearranging Regions vs Excitations*, arXiv:2103.08060.

[161] D. Ganapathi, A. K. Sood, and R. Ganapathy, *Structural Origin of Excitations in a Colloidal Glass-Former*, J. Chem. Phys. **156**, 214502 (2022).

[162] S. Whitelam, L. Berthier, and J. P. Garrahan, *Renormalization Group Study of a Kinetically Constrained Model for Strong Glasses*, Phys. Rev. E **71**, 026128 (2005).

[163] G. H. Fredrickson and H. C. Andersen, *Facilitated Kinetic Ising Models and the Glass Transition*, J. Chem. Phys. **83**, 5822 (1985).

[164] W. Kob and H. C. Andersen, *Kinetic Lattice-Gas Model of Cage Effects in High-Density Liquids and a Test of Mode-Coupling Theory of the Ideal-Glass Transition*, Phys. Rev. E **48**, 4364 (1993).

[165] C. Toninelli, G. Biroli, and D. S. Fisher, *Jamming Percolation and Glass Transitions in Lattice Models*, Phys. Rev. Lett. **96**, 035702 (2006).

[166] C. Toninelli and G. Biroli, *A New Class of Cellular Automata with a Discontinuous Glass Transition*, J. Stat. Phys. **130**, 83 (2008).

[167] C. Toninelli, G. Biroli, and D. S. Fisher, *Spatial Structures and Dynamics of Kinetically Constrained Models of Glasses*, Phys. Rev. Lett. **92**, 185504 (2004).

[168] P. Sollich and M. R. Evans, *Glassy Time-Scale Divergence and Anomalous Coarsening in a Kinetically Constrained Spin Chain*, Phys. Rev. Lett. **83**, 3238 (1999).

[169] P. Sollich and M. R. Evans, *Glassy Dynamics in the Asymmetrically Constrained Kinetic Ising Chain*, Phys. Rev. E **68**, 031504 (2003).

[170] S. Plimpton, *Fast Parallel Algorithms for Short-Range Molecular Dynamics*, J. Comput. Phys. **117**, 1 (1995).

[171] N. Bailey, T. Ingebrigtsen, J. S. Hansen, A. Veldhorst, L. Bøhling, C. Lemarchand, A. Olsen, A. Bacher, L. Costigliola, U. Pedersen, H. Larsen, J. C. Dyre, and T. B. Schrøder, *RUMD: A General Purpose Molecular Dynamics Package Optimized to Utilize GPU Hardware Down to a Few Thousand Particles*, SciPost Phys. **3**, 038 (2017).

[172] J. D. Eaves and D. R. Reichman, *Spatial Dimension and the Dynamics of Supercooled Liquids*, Proc. Natl. Acad. Sci. U.S.A. **106**, 15171 (2009).

[173] M. Adhikari, S. Karmakar, and S. Sastry, *Spatial Dimensionality Dependence of Heterogeneity, Breakdown of the Stokes-Einstein Relation, and Fragility of a Model Glass-Forming Liquid*, J. Phys. Chem. B **125**, 10232 (2021).

[174] L. Berthier, P. Charbonneau, and J. Kundu, *Finite Dimensional Vestige of Spinodal Criticality above the Dynamical Glass Transition*, Phys. Rev. Lett. **125**, 108001 (2020).

[175] J. D. Stevenson, J. Schmalian, and P. G. Wolynes, *The Shapes of Cooperatively Rearranging Regions in Glass-Forming Liquids*, Nat. Phys. **2**, 268 (2006).

[176] T. Speck and D. Chandler, *Constrained Dynamics of Localized Excitations Causes a Non-Equilibrium Phase Transition in an Atomistic Model of Glass Formers*, J. Chem. Phys. **136**, 184509 (2012).

[177] M. Isobe, A. S. Keys, D. Chandler, and J. P. Garrahan, *Applicability of Dynamic Facilitation Theory to Binary Hard Disk Systems*, Phys. Rev. Lett. **117**, 145701 (2016).

[178] M. Campo and T. Speck, *Dynamical Coexistence in Moderately Polydisperse Hard-Sphere Glasses*, J. Chem. Phys. **152**, 014501 (2020).