# Dynamical Facilitation Governs the Equilibration Dynamics of Glasses

Rahul Chacko, François Landes, Giulio Biroli, Olivier Dauchot, Andrea Liu, David Reichman

► To cite this version:

Rahul Chacko, François Landes, Giulio Biroli, Olivier Dauchot, Andrea Liu, et al. Dynamical Facil- itation Governs the Equilibration Dynamics of Glasses. Physical Review X, 2024, 14 (3), pp.031012. 10.1103/PhysRevX.14.031012. hal-04749267

HAL Id: hal-04749267

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

Submitted on 8 Nov 2024

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
# Dynamical Facilitation Governs the Equilibration Dynamics of Glasses

Rahul N. Chacko<sup>1</sup>,<sup>†</sup> François P. Landes<sup>2</sup>,<sup>†</sup> Giulio Biroli,<sup>3</sup> Olivier Dauchot<sup>4</sup>,<sup>†</sup> Andrea J. Liu<sup>1</sup>,<sup>†</sup> and David R. Reichman<sup>5</sup>

<sup>1</sup>Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, Pennsylvania 19104, USA

<sup>2</sup>Université Paris-Saclay, CNRS, INRIA, Laboratoire Interdisciplinaire des Sciences du Numérique,
TAU team, 91190 Gif-sur-Yvette, France

<sup>3</sup>Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne
Université, Université Paris-Diderot, Sorbonne Paris Cité, 75005 Paris, France

<sup>4</sup>UMR Gulliver 7083 CNRS, ESPCI, PSL Research University, 10 rue Vauquelin, 75005 Paris, France

<sup>5</sup>Department of Chemistry, Columbia University, 3000 Broadway, New York, New York 10027, USA

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

(Received 1 January 2024; revised 18 March 2024; accepted 7 May 2024; published 19 July 2024)

Convincing evidence of domain growth in the heating of ultrastable glasses suggests that the equilibration dynamics of supercooled liquids could be driven by a nucleation and growth mechanism. We investigate this possibility by simulating the equilibration dynamics of a model glass during both heating and cooling between poorly and well-annealed states. Though we do observe the growth of domains during heating, we find that domains are absent during cooling. This absence is inconsistent with classical nucleation theory. By comparing the equilibration dynamics of our glass with that of two models with kinetic constraints, we demonstrate that dynamical facilitation generically leads to heating driven by domain growth and cooling without domains. Our results provide strong evidence that dynamical facilitation, not nucleation and interfacial-tension-driven domain growth, is the driving mechanism for the equilibration dynamics of glass formers.

DOI: 10.1103/PhysRevX.14.031012

Subject Areas: Computational Physics, Soft Matter,
Statistical Physics

## I. INTRODUCTION

The phenomenology of glasses and supercooled liquids (i.e., glasses in metastable equilibrium such that detailed balance is satisfied on all relevant timescales) bears a striking resemblance to that of crystalline solids [1–8]. As with transformations of Eshelby disclinations in crystals [9], local plastic events in these systems perturb their neighborhoods with long-ranged quadrupolar displacement fields [10,11]. A large body of experimental evidence also finds that ultrastable glasses evolve upon heating in a manner consistent with the Avrami framework of classical nucleation theory [12]: Fronts of high-temperature supercooled liquid advance at constant speed into the bulk ultrastable glass [5,13–29], with the fraction of material retaining ultrastable glass structure decaying with time according to a compressed exponential law [15,16,18,30]. By contrast, poorly annealed (“ordinary”) glasses are found to transform homogeneously into the supercooled liquid upon heating above the glass

transition temperature [20,27]. This difference of behavior, also observed in pressure-controlled molecular dynamics (MD) simulations [7,8,30,31], suggests a process analogous to crystal freezing and melting [6–8,32] and raises the possibility that a thermodynamic first-order transition mediated by nucleation and domain growth separates poorly annealed glass from well-annealed glass [6,8,33].

Indeed, as with crystals, a notion of structural order also applies to supercooled liquids. In these amorphous systems, local equilibrium structure constrains [34] and is constrained by [35] neighboring structure up to a finite distance ξ<sub>PTS</sub> (the “point-to-set” length) away [36]. Random first-order transition theory (RFOT), known to be exact in infinite dimension [37], predicts systems below a critical temperature T<sub>MCT</sub> to be separated into a mosaic of ξ<sub>PTS</sub>-sized domains of mutually incompatible amorphous order [1,36,38]. As demonstrated by randomly pinned glasses [39] and supercooled liquids coupled to a reference configuration [3], amorphous order can in principle lead to first-order transitions in amorphous systems, supporting the possibility of a first-order transition in supercooled liquids [40].

Nucleation and domain growth is not the only available explanation for the equilibration dynamics of glass formers, however. Particles in supercooled liquids are trapped into “cages” by their nearest neighbors, rattling around for long

<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>
periods between cage-escape events in which they hop out of one cage and into another [41,42]. As a local plastic event within an otherwise elastic material, a cage escape weakens its neighborhood, creating the conditions for more cage escapes [11]. Such a process, in which rearrangements propagate mobility to their local neighborhoods, is known as dynamical facilitation [43,44] and plays an important role in glassy liquids [45]. This self-propagation of mobility results in avalanches of plastic events correlated in space and time [42]. On long timescales, these avalanches appear as large clusters of spatially correlated mobility that grow as a function of the timescale [43,46], resembling domain growth. Dynamical facilitation plays an important role in the equilibrium dynamics of super-cooled liquids, and it could also be the key mechanism at play in equilibration dynamics, i.e., relaxation toward equilibrium. Indeed, many of the phenomena observed in glasses upon heating or cooling are also seen in kinetically constrained models (KCMs), spin-lattice models in which spins interact only via facilitation [47] upon heating [2,48,49] or cooling [50]. Dynamical facilitation therefore offers a possible alternative explanation for the behavior of glasses during heating [20,51] and cooling [50]. In this paper, we show that it is indeed dynamical facilitation that governs the equilibration dynamics of glasses, not nucleation followed by domain growth.

We first address the question of nucleation interpreted in the classical sense [52] of a competition between bulk and interfacial free energies which sets a critical size above which randomly nucleated domains will grow and below which they will shrink. We distinguish nucleation, which involves a thermodynamic drive toward domain growth as small domains (of linear size smaller than the nucleation length) are suppressed, from the evolution of structure via random, thermally activated rearrangements, with no suppression of small domains. In the latter case, an additional mechanism is required to explain domain growth. In the second part of this paper, we show that dynamical facilitation fills this role.

We conduct MD simulations of a polydisperse glass at fixed temperature $T_{\text{eq}}$ and pressure $P$ for two sets of initial conditions. In simulation runs using the first set of initial conditions, which we refer to as homogeneous simulations, the initial state corresponds to equilibrium at an initial temperature $T_0$ different from the thermostat temperature $T_{\text{eq}}$. In simulation runs using the second type of initial condition, which we refer to as slab simulations, we introduce a system-percolating slab of $T_{\text{eq}}$ equilibrium structure lying in the $x=0$ plane that interrupts the $T_0$ structure of the surrounding system [see Figs. 4(a) and 4(b)]. In this geometry, the slab can expand while keeping the size and shape of its boundary fixed, avoiding any putative interfacial free-energy cost as it expands. In a classical nucleation theory scenario, the slab represents a nucleus above the critical nucleus size. We fix

$T_{\text{low}} = \min\{T_0, T_{\text{eq}}\}$ and $T_{\text{high}} = \max\{T_0, T_{\text{eq}}\}$ so that any critical point separating $T_0$ and $T_{\text{eq}}$ during heating will also do so during cooling.

These simulations provide a number of tests for nucleation and domain growth:

(1) We can characterize the local structure and dynamics, and simply look (by eye) for growing domains of $T_{\text{eq}}$ structure in individual realizations of the homogeneous simulation.

(2) We can also look for evidence of a nucleation time in these simulations.

(3) We can quantitatively track domain growth via the displacement of the slab boundary in slab simulations.

(4) We can assess the roughness of domains, both qualitatively in individual realizations of the homogeneous simulation and quantitatively (via the width of the slab boundary) in the slab simulation.

As we will show, all four tests argue against nucleation, with the most striking evidence provided by the absence of domain growth during cooling, as we shall discuss.

We validate these tests by first applying them to a crystallizing monodisperse system with temperatures $T_{\text{low}}$ and $T_{\text{high}}$ straddling its melting point, verifying our expectations for the case of a genuine first-order transition. We then apply this methodology to our polydisperse model glass, obtaining results inconsistent with the nucleation and growth mechanism. Finally, we support our claim that dynamical facilitation is responsible for the equilibration dynamics of glasses by applying the same four tests to two plaquette models, the triangular plaquette model (TPM) [50,53,54] and the square pyramidal plaquette model (SPPM) [55].

Plaquette models are spin-lattice models with a dual representation in terms of plaquettes with trivial thermodynamics and kinetically constrained dynamics [56] corresponding to KCMs in the plaquette representation. The TPM and SPPM are similar by design, and the mechanisms underlying dynamical facilitation in these two systems differ greatly from the elastoplastic mechanism for facilitation in supercooled liquids and glasses [11]. We find that just like our polydisperse model glass, these plaquette models exhibit domain growth during heating and none during cooling. For all three systems, domain growth during heating proceeds at a constant speed that persists even as the system relaxes. The only notable difference between the systems is the evolution of domain roughness, which is different for all three systems and therefore not a test for facilitation-driven dynamics.

In investigating nucleation and domain growth and the role of dynamical facilitation as above, we also establish two side results of particular interest. First, we provide direct support for a two-state equilibration scenario in which glasses evolve from their initial state into the target equilibrium state (set by the thermostat) via direct
transformation of structure into the target equilibrium structure, as opposed to via one or more intermediate structural states. This scenario has long been the experimental consensus for heating [6,14,22,24,28,29,57,58]. However, we notably also find it to hold for cooling. Second, to support dynamical facilitation as the governing mechanism for the equilibration dynamics of glasses, we show that the equilibration dynamics of the TPM in the slab geometry, which closely resembles that of the supercooled liquid, can be captured using the same phenomenological model for both heating and cooling.

The structure of this paper is as follows. In Sec. II, we discuss in detail our model crystal and glass systems, and the way in which our simulation protocols are applied to them. In Secs. III and IV, we apply our four tests to the model crystal and glass systems. We defer a broader discussion of the significance of these results until Sec. VI. In Sec. V, we introduce our two plaquette models and show the results of applying the four tests to them. Finally, in Sec. VI, we synthesize the results for the four different systems, discussing the conclusions that can be drawn from comparing these results with each other. Readers may find it helpful to refer to Sec. VI after they read about each of the four model systems (two particulate systems and two plaquette models) in turn.

## II. SYSTEM

Our three-dimensional monodisperse crystal and two-dimensional polydisperse glass systems comprise particles of uniform mass $m$ interacting via the pair potential

$$V(\tilde{r}) = \begin{cases} V_0(\tilde{r}^{-12} + c_0 + c_2\tilde{r}^2 + c_4\tilde{r}^4), & \tilde{r} \le 1.25, \\ 0, & \tilde{r} > 1.25, \end{cases} \quad (1)$$

where

$$\tilde{r} = 2r/[(\sigma_i + \sigma_j)(1 - 0.2|\sigma_i - \sigma_j|)].$$

Here, $r$ is the separation between particles $i$ and $j$ with sizes $\sigma_i$ and $\sigma_j$, and the constants $c_0$, $c_2$, and $c_4$ are chosen such that $V(1.25) = V'(1.25) = V''(1.25) = 0$ (Appendix A 1). Following Ref. [59], for the polydisperse system, $\sigma_i$ for each particle is drawn independently from a distribution with probability density $\propto \sigma^{-3}$ in the interval $[\sigma_{min}, \sigma_{max}]$, with $\sigma_{min}$ and $\sigma_{max}$ chosen such that the coefficient of variation of $\sigma$ is $c_\sigma = 0.23$. This system, introduced in Ref. [59], is designed to allow equilibration to temperatures well below $T_{MCT}$ using swap Monte Carlo simulations. We take $V_0$, $\bar{\sigma}$, $m$, and $V_0/k_B$ to be our energy, length, mass, and temperature units, where $k_B$ is the Boltzmann constant and where $\bar{\sigma}$ is the expected value of the particle size distribution. The system coordinates are centered at the origin, and periodic boundary conditions are imposed along each axis of the system.

We prepare our initial states using constant-volume Monte Carlo simulations implemented as described in Appendix A 2. For the monodisperse case, the homogeneous simulations have $N = 27648$ particles and initial lengths $48, 12\sqrt{3}$, and $8\sqrt{6}$ along the $x$, $y$, and $z$ axes, respectively, corresponding to number density $\rho = \sqrt{2}$. These system dimensions are compatible with a perfect fcc crystal, the ground state for this system (Appendix A 7). For melting, we heat the system from an initial configuration equilibrated to $T_{low} = 1.90$ using a thermostat at $T_{high} = 2.35$, while for freezing, we cool the system from an initial configuration equilibrated to $T_{high} = 2.10$ using a thermostat at $T_{low} = 1.48$.

For the slab simulations in the monodisperse system, we fix $T_{low} = 2.0$ and $T_{high} = 2.1$, which straddle the melting point. We preserve the compatibility of the system dimensions with a perfect fcc crystal by not attempting to eliminate the initial pressure gradient at the slab boundary in slab simulations, other than by our choice of narrowly separated $T_{low}$ and $T_{high}$. Instead, we double the number of particles and the system size along the $x$ axis, as compared to the systems used for homogeneous simulations. This allows for time between the early-time pressure wave dissipating and the left and right slab fronts colliding at the periodic boundary of the system. For freezing, we define our frozen slab to be the region with $x$ coordinate satisfying $|x| < 6$. For melting, our preparation protocol (Appendix A 2) requires us to specify the size of the frozen nonslab region instead, leaving the size of the melted slab region to be an emergent property of our preparation protocol dependent on the size of the interface. We define the nonslab region to be that with $x$ coordinate satisfying $|x| \ge 8$.

For the polydisperse case, we simulate $N = 20000$ particles and fix $T_{low} = 0.06$ and $T_{high} = 0.14$ for both homogeneous and slab simulations. The homogeneous heating simulations have initial lengths 200.00 and 100.00 along the $x$ and $y$ axes, respectively, corresponding to number density $\rho = 1$. At this density, $T_{MCT} \approx 0.12$ [46], such that (accounting for density changes [60]) $T_{low}$ and $T_{high}$ straddle $T_{MCT}$. The homogeneous heating simulations follow these systems to steady state at $T_{high}$, corresponding to (average) system lengths 204.70 and 102.35 ($\rho \approx 0.955$) along the $x$ and $y$ axes. These $T_{high}$ equilibrium configurations are then used as the initial $T_0 = T_{high}$ states for the homogeneous cooling simulations.

Our choice of $T_{low} = 0.06$ and $T_{high} = 0.14$ comes from practical considerations. We require $T_{low} < T_{MCT}$ to ensure that $T_{low}$ is low enough for facilitation to play an important role in the equilibrium dynamics [43] and for amorphous order to set in [1,38]. The lower we set $T_{low}$, the stronger we expect dynamical facilitation to be and the slower we expect the dynamics to be. We choose $T_{low} = 0.06$ as a compromise between strong facilitation and dynamics fast
enough that we are able to see appreciable structural evolution on the timescale of a week of run-time. While this value of $T_{low}$ corresponds to a less stable ultrastable glass than is typically encountered in the experimental literature [25,32], it was found in Ref. [8] that choosing a much cooler $T_{low}$ ($T_{low}$ = 0.035) only quantitatively changes the behavior of fronts observed during heating. This suggests that $T_{low}$ = 0.06 is low enough to capture the physics of extremely stable systems.

At number density $\rho \approx 0.955$ the onset temperature is $T_{onset} \approx 0.15$ [46]. We check that increasing $T_{high}$, even above $T_{onset}$, does not affect our results for heating (in terms of the four tests). However, choosing $T_{high}$ = 0.14, below $T_{onset}$, ensures that we are probing a potential ultrastable-to-poorly-annealed glass transition, rather than a transition between an ultrastable glass and a high-temperature liquid.

For the slab simulations in the polydisperse system, we define the slab to be the region satisfying $|x| < b$, where $b = 25.00$ for heating and 25.56 for cooling. In slab simulations, we greatly reduce the initial pressure gradient at the slab boundary by first equilibrating the system in its entirety to $T_{eq}$, with system dimensions 200.00 and 100.00 if $T_{eq} = T_{low}$ or 204.70 and 102.35 if $T_{eq} = T_{high}$, then shrinking or stretching the nonslab region along the x axis such that $\rho \approx 0.955$ in this region if $T_0 = T_{high}$ or $\rho = 1$ if $T_0 = T_{low}$, before finally evolving this nonslab region with Monte Carlo moves at temperature $T_0$ until it reaches steady state.

Our MD simulations are performed in LAMMPS [61] (see Appendix A 3), adopting a Nosé-Hoover thermostat and a Berendsen barostat that keeps the relative sizes of the system along each axis constant. In order to minimize the perturbation to the system due to the start-up of the barostat, we fix the pressure such that the virial contribution to the pressure at the beginning of the MD simulation matches that of the target $T_{eq}$ equilibrium state. For any given MD simulation, we average our results over 100 realizations, except when explicitly considering individual realizations or when performing the homogeneous simulations on the polydisperse glass, in which case we average over only ten realizations.

We define time $t = 0$ as the time at which the MD simulations begin. For our analysis, it is necessary to measure the distributions of structural and dynamical quantities in the initial, untransformed state during homogeneous simulations, as well as at equilibrium at the target temperature $T_{eq}$. Our methodology for obtaining these distributions can be found in Appendix A 4.

### III. CRYSTAL CASE STUDY

We first validate our tests by applying them to the first-order melting and freezing transitions of a monodisperse system of particles, in which the growth of nuclei is driven by bulk free-energy differences [62], with structural

frustration at the boundary of the nucleus suppressing roughness and promoting domain growth [63,64]. We track the evolution of structure in this system using Tong and Tanaka's $\Theta$ order parameter [65], a local measure of packing inefficiency that vanishes in the limit of perfect fcc structure. We also track the corresponding evolution of the mobility $\mu(t)$ defined here for each particle as the magnitude of its inherent-state displacement between times $t$ and $t + 0.1$.

**Test 1:** Here we simply look for the emergence and growth of $T_{eq}$-structure domains in individual realizations of the homogeneous simulation. Figures 1(a) for melting and 1(e) for freezing (see also Appendix B 5) both show a domain of $T_{eq}$ structure and dynamics which grows to span the system along the narrow horizontal axis. At this point, the domain can expand without gaining additional interfacial area, making this width an upper bound for the critical nucleus size, and we see the domain quickly grow to envelop the entire system.

**Test 2:** Here we look for evidence of a waiting time for nucleation. In Fig. 1(a), we see that the time between forming a system-spanning nucleus and melting the entire system is of order $t \sim 10^0$, while it takes a time of order $t \sim 10^2$ for this system-spanning nucleus to appear. For cooling [Fig. 1(e)], these timescales are $t \sim 10^2$ and $t \sim 10^3$, respectively. This scale separation between the nucleation and domain growth times holds true across realizations (Sec. VI).

**Two-state scenario:** Our third and fourth tests involve slab simulations. To ensure that the results from these are relevant to systems without an artificial slab, we first demonstrate that the structural evolution observed in the homogeneous simulations corresponds to a two-state scenario, such that at any time $t$, a fraction $f$ of the system has local structure distributed according to equilibrium at the thermostat temperature $T_{eq}$ (the target state), while the remaining fraction $1-f$ retains its initial distribution of structure, corresponding to equilibrium at temperature $T_0$ (the initial state). In this scenario, widely expected to hold for the heating of ultrastable glasses [6,14,22,24,28,29,57,58], the distribution $P(\phi, t)$ of a local structural indicator $\phi$ with initial-state distribution $P_0(\phi)$ and target-state distribution $P_{eq}(\phi)$ will evolve according to

$$P(\phi, t) = (1 - f(t))P_0(\phi) + f(t)P_{eq}(\phi), \quad (2)$$

interpolating between $P_0(\phi)$ and $P_{eq}(\phi)$ as the transformed fraction $f$ grows with time. Given distributions $P_0(\phi)$, $P_{eq}(\phi)$, and $P(\phi, t)$ for a structural quantity $\phi$, we can extract the transformed fraction $f$ from $P(\phi, t)$ under the assumption that Eq. (2) holds (see Appendix A 8). Alternatively, Eq. (2) implies that
<figure><img src="image_3.png" /><figcaption>FIG. 1. Crystal melting (a)–(d) and freezing (e)–(h) transitions in a monodisperse system. (a),(e) Snapshots in the $z=0$ plane of (top row) the $\Theta$ order parameter and (bottom row) the mobility $\mu$ as a nucleus appears and grows in an individual trajectory. (b),(c),(f),(g) Distributions (b),(f) $P(\log_{10} \Theta)$ and (c),(g) $P(\log_{10} \mu)$ at times (b),(c) $t = 100, 200, \dots, 600$ and (f),(g) $t = 1000, 2000, \dots, 7000$ (later time, lighter color). The red curves in (b) and (c) are equilibrium distributions. (d),(h) Untransformed fraction $1-f$ for $f$ calculated from $P(\log_{10} \Theta)$ (dashed red), $P(\log_{10} \mu)$ (dashed green), potential energy $\bar{V}$ (solid orange), inherent-state potential energy $V_{\text{IS}}$ (solid blue), $\bar{\Theta}$ (solid red), $\overline{\log_{10} \Theta}$ (dotted red), $\bar{\mu}$ (solid green), and $\overline{\log_{10} \mu}$ (dotted green).</figcaption></figure>
$$1 - f(t) = \frac{\bar{\phi}_{\text{eq}} - \bar{\phi}(t)}{\bar{\phi}_{\text{eq}} - \bar{\phi}_0}, \qquad (3)$$

where $\bar{\phi}_0$, $\bar{\phi}_{\text{eq}}$, and $\bar{\phi}(t)$ are the mean values of the quantity $\phi$ distributed according to $P_0(\phi)$, $P_{\text{eq}}(\phi)$, and $P(\phi, t)$, respectively, allowing us to calculate $f$ from the evolving mean of $\phi$.

If the two-state scenario holds, we should obtain the same value for the transformed fraction $f$, regardless of whether we use the distribution $P(\phi, t)$ or the mean $\bar{\phi}(t)$ to calculate it, and regardless of which local structural variable $\phi$ we choose. We can therefore test the validity of the two-state scenario by calculating $f$ from both $P(\phi, t)$ and $\bar{\phi}$ for a variety of different structural quantities $\phi$ (Appendix B 2). In practice, we only calculate $f$ from $P(\phi, t)$ for a single structural quantity $\phi$, since this is more complicated to do (Appendix A 8). Though it is not a requirement of the two-state equilibration scenario, we also check if we can obtain $f$ from $P(\phi, t)$ and $\bar{\phi}$ when $\phi$ is a dynamical, rather than structural, quantity, since the coupling of structure and dynamics is a question of physical interest.

In Figs. 1(b), 1(c), 1(f), and 1(g), we see that the $T_0$ and $T_{\text{eq}}$ distributions of $\log_{10} \Theta$ and $\log_{10} \mu$ are well separated [66], making the validity of the two-state scenario immediately apparent. We calculate $1-f$ from these distributions and plot the resulting curves against $1-f$ calculated from $\bar{\phi}$ for $\phi$ corresponding to $\Theta$, $\log_{10} \Theta$, $\mu$, $\log_{10} \mu$, the potential energy $V$ associated with each particle, and the inherent-state potential energy $V_{\text{IS}}$ associated with each particle. As shown in Figs. 1(d) and 1(h), these different methods of calculating $1-f$ agree well with one another, confirming the validity of the two-state scenario and therefore the soundness of the slab simulation as a test for domain growth, and showing that structure and dynamics are coupled together in crystal freezing and melting.

Having validated the slab simulations, we calculate the (ensemble-averaged) $\Theta$ profile along the $x$ axis shown as transparent curves in Figs. 2(c) and 2(d). We identify the interface at the boundary between a $T_{\text{eq}}$-structure slab and the surrounding $T_0$-structure region by fitting this profile with a logistic function

$$\Theta(\tilde{x}) = \frac{\Theta_{\text{eq}} + \Theta_0}{2} + \frac{\Theta_{\text{eq}} - \Theta_0}{2} \tanh \frac{\tilde{x} - x_*}{2\lambda}, \qquad (4)$$

where the argument $\tilde{x} = 2|x|/L_x$ accounts for changes in $L_x$ (the system size along $x$) due to the barostat and takes advantage of the $x \mapsto -x$ symmetry of our slab geometry [67]. The parameters $\Theta_{\text{eq}}$, $\Theta_0$, $x_*$, and $\lambda$ are fit at each snapshot of time $t$, representing, respectively, the plateau value of $\Theta$ within the slab of $T_{\text{eq}}$ structure, the plateau value of $\Theta$ outside this slab, the interface center, and the interface width. The order parameter $\Theta$ is invariant with respect to isotropic expansion, and is therefore insensitive to the

<figure><img src="image_4.png" /><figcaption>FIG. 2. Dynamics of melting and freezing for temperatures $T_{\text{low}} = 2.0$ and $T_{\text{high}} = 2.1$ tightly straddling the critical point, given an artificial target-state slab at time $t = 0$. (a),(b) Slice along $z = 0$ of the initial configuration, colored according to $\Theta$ on a logarithmic scale from $10^{-3}$ to $10^{-1}$ (higher $\Theta$, lighter color). (c),(d) Average $\Theta(\tilde{x})$ at times $t = 0, 50, \dots, 300$ for melting, and $t = 0, 40, \dots, 400$ for freezing (later time, lighter color). Transparent curves correspond to raw data, solid curves to fits. Here, $\tilde{x} = 2|x|/L_x$, where $L_x$ is the time-dependent system size along the $x$ axis. (e)-(h) Change, relative to their initial value, of the (e),(f) interface center $x_*$ and (g),(h) relative width $\lambda/\lambda_0$ obtained when fitting the profiles to Eq. (4). Here, the initial interface width $\lambda_0$ has values $\lambda_0 = 0.0213$ and $0.0159$ in (e) and (f), respectively. The dashed lines in (e) and (f) are linear fits.</figcaption></figure>

transient effect of the pressure gradient at time $t = 0$, making it an especially suitable choice for our analysis.

Test 3: Here we check for domain growth in the slab simulation. As we see in Figs. 2(c) and 2(d), the slab of $T_{\text{eq}}$ structure indeed grows in both melting and freezing simulations. Figures 2(e) and 2(f) show that after an initial transient of around 150 time units, which we associate with the dissipation of the pressure gradient present in the initial configuration, the displacement $\Delta x_* = x_*(t) - x_*(0)$ of the interface center increases at constant speed. Both this constant front speed and the two-state equilibration scenario accord with the description of domain growth within the Avrami formalism [12] for crystallization and melting.

Test 4: In our final test, we look for evidence that interfacial area is being minimized. In Fig. 1(a), we see a transition from the domain boundary being preferentially aligned along a crystal lattice vector to being aligned horizontally, while in Fig. 1(e), the domain boundary
has no particular initial orientation (the surrounding struc-
ture is amorphous), but also quickly becomes horizontal.
This is evidence of an interfacial energy penalty promoting
a reduction of interfacial area. For the slab simulation, the
initial width $\lambda_0$ of the interface has to do with the ability of
the $T_0$ and $T_{\text{eq}}$ structure to constrain other structure nearby,
but the two-state equilibration scenario means we can
associate further growth in $\lambda$ with roughening of the
interface. In Figs. 2(g) and 2(h), after the initial pres-
sure-driven transient ($t \ge 150$), fluctuations in $\Delta\lambda = \lambda(t) -
\lambda_0$ are comparable in magnitude to the size of the noise
floor. This suppression of roughening is further evidence of
an interfacial energy penalty promoting a reduction of
interfacial area.

In summary, the results of our four tests match our
expectations for a nucleation and domain growth process,
validating our methodology.

IV. GLASSES

By adjusting our methodology for calculating Θ (Appendix A 5) and redefining the mobility μ(t) to be the inherent-state displacement of a given particle between times t and t + 100 (accounting for Mermin-Wagner fluctuations as described in Appendix A 6), we can use Θ and μ to characterize the local structure and dynamics of our two-dimensional polydisperse glass.

Test 1: Here we look for growing $T_{\text{eq}}$-structure domains in an individual realization of the homogeneous simulation (Appendix B 5). In Fig. 3(a), we see growing high-Θ and matching high-μ domains, representing clear evidence that heating is mediated by domain growth. We do not, however, see any visible low-Θ or low-μ domains during cooling in Fig. 3(e), suggesting a lack of domain growth in this case. This is a first and important difference from the crystal case.

Two-state scenario: Before applying our second test, which will benefit from the results of the slab simulations, we confirm the validity of the two-state scenario. We find good agreement between $1 - f$ calculated from the different structural quantities $P(\Theta)$, $\bar{V}$, $\bar{V}_{\text{IS}}$, $\bar{\Theta}$, and $\overline{\log_{10} \Theta}$ in Figs. 3(d) and 3(h), establishing the validity of the two-state scenario for both heating and cooling. These figures also highlight a major difference between heating and cooling: the coupling of structure and dynamics. For heating, this coupling is evident from the matching high-Θ and high-μ domains and the good agreement between $1 - f$ calculated using any of our structural quantities with $1 - f$ calculated from $P(\log_{10} \mu)$ or $\bar{\mu}$ [Fig. 3(d)] demonstrates this yet further. [We do see a deviation from this trend of $1 - f$ calculated from $\overline{\log_{10} \mu}$ due to the rightward shift of the $P_0(\log_{10} \mu)$ mode with time in Fig. 3(c); see also Appendix B 1.]

For cooling, however, $1 - f$ calculated from $\bar{\mu}$ or $\overline{\log_{10} \mu}$
decays much faster than $1 - f$ calculated from structural

quantities [Fig. 3(h)], and while $P(\Theta)$ evolves in accor-
dance with the two-state scenario of Eq. (2) [Fig. 3(f); see
also Appendix A 8], the evolution of $P(\log_{10} \mu, t)$, com-
prising a continuous shift of the initial distribution
$P_0(\log_{10} \mu)$ toward lower mobility [Fig. 3(g)], does not.
We therefore see that structure and dynamics are decoupled
during cooling.

Test 2: Our second test is to look for evidence of a nucleation time. With no evidence of domain growth during cooling, we focus here on heating. In a classical nucleation theory scenario, small domains are suppressed, so the nucleation time should be much larger than the timescale of individual rearrangements (successful nucleation should require many attempts), as was the case for the crystal (Fig. 1). We see in Fig. 3(a) that $T_{\text{high}}$-structure domains first appear at time $t \lesssim 10^4$. As can be seen in Fig. 3(d) [see also Fig. 8(c)], this timescale is no larger than the timescale $t \sim 10^4$ for the relaxation of the system as the domains grow.

This absence of a significant waiting time (and hence,
putative nucleation time) is reinforced by the results of the
slab simulations. In these, we find [Fig. 4(e)] that $T_{\text{high}}$-
structure fronts advance with speed $1.57 \times 10^{-3}$ into the
bulk. A domain growing isotropically at this speed would
take $6.37 \times 10^4$ time units to span the $x$ axis ($L_x \approx 200$).
This is the time required for a growing, initially point-sized
domain to fully relax the system into the $T_{\text{high}}$ state, without
allowing for any waiting time. We see in Fig. 3(d) that this
is a good quantitative match with the actual time taken for
the system to fully relax into the $T_{\text{high}}$ state, implying that
the waiting time is small on the timescale of domain
growth. The high density of disjoint growing domains
visible at $t = 20000$ in Fig. 3(a) attests to this too.

Test 3: Here we see if the $T_{\text{eq}}$-structure slabs grow in the
slab simulation. As we have already mentioned,
$T_{\text{eq}} = T_{\text{high}}$-structure slabs indeed grow in heating simu-
lations, with the center of the slab interface $x_*$ advancing
with constant speed up to the system size limit [Figs. 4(c)
and 4(e)]. In the case of cooling, however, we see that
slabs of $T_{\text{eq}} = T_{\text{low}}$ structure do not grow at all [Figs. 4(d)
and 4(f)]. This shows that the reason no domains of $T_{\text{low}}$
structure are visible in Fig. 3(e) is that such domains do not
grow in the first place. In a classical nucleation theory
scenario, the competition between bulk and interfacial free
energies would apply to both heating and cooling, with the
free-energy difference between $T_{\text{low}}$ and $T_{\text{high}}$ structure in
fact being larger during cooling given the lower thermostat
temperature $T_{\text{eq}} = T_{\text{low}}$. This asymmetry and lack of
domain growth during cooling is therefore a strong piece
of evidence against this scenario.

Test 4: Here we investigate whether the size of the
interface is being minimized. The high-mobility domains
visible in Fig. 3(a) clearly have rough edges, and we see no
alignment along the horizontal axis to minimize the length
of the interface. Indeed, spanning the system appears to
<figure><img src="image_5.png" /><figcaption>FIG. 3. Heating (a)–(d) and cooling (e)–(h) of a two-dimensional polydisperse supercooled liquid between temperatures 0.06 and 0.14. (a),(d) Snapshots of the mobility $\mu$ as the system evolves, exhibiting domain growth during heating and quadrupoles during cooling, once rearrangements are sufficiently rare. (b),(c),(f),(g) Evolution of the (b),(f) $\Theta$ and (c),(g) $\log_{10} \mu$ distributions at times (b),(c) $t = 6000, 12000, \dots, 60000$ and (f),(g) $t = 10^{2.0}, 10^{2.5}, \dots, 10^{6.0}$ (later time, lighter color). The red curves in (f) and (g) are equilibrium distributions obtained from equilibrium trajectories at the target temperature and pressure. (d),(h) Untransformed fraction $1-f$ for $f$ calculated from $P(\log_{10} \Theta)$ (dashed red), $P(\log_{10} \mu)$ (dashed green), $\bar{V}$ (solid orange), $\bar{V}_{\text{IS}}$ (solid blue), $\bar{\Theta}$ (solid red), $\overline{\log_{10} \Theta}$ (dotted red), $\bar{\mu}$ (solid green), and $\overline{\log_{10} \mu}$ (dotted green).</figcaption></figure>

have no impact on the growth of domains in Fig. 3(a), unlike in the crystal case, where spanning the system was associated with an acceleration of domain growth and alignment of the domain boundaries with the short axis of the system [Figs. 1(a) and 1(e)]. The growth of the interface

width $\lambda$ in slab simulations provides further evidence against a thermodynamic drive to minimize the size of the interface. After an early-time transient ($t \lesssim 10^3$) of diffusive growth ($\Delta\lambda \sim t^{1/2}$), we see that the growth in $\lambda$ during heating [Fig. 4(g)] is superdiffusive until the time
<figure><img src="image_6.png" /><figcaption>FIG. 4. Dynamics of (left column) heating and (right column) cooling supercooled liquids between temperatures $T_{low} = 0.06$ and $T_{high} = 0.14$ straddling $T_{MCT}$, given an artificial target-state slab at time $t = 0$. (a),(b) Initial configuration for (a) heating and (b) cooling, colored according to $\Theta$ on a linear scale from 0.0 to 0.1 (higher $\Theta$, lighter color). (c),(d) Average $\Theta(\tilde{x})$ at times (c) $t = 0, 2000, \dots, 20000$ and (d) $t = 10^{2.0}, 10^{2.5}, \dots, 10^{6.0}$ (later time, lighter color). Transparent curves correspond to raw data, solid curves to fits. As before, $\tilde{x} = 2|x|/L_x$, where $L_x$ is the time-dependent size of the system along the $x$ axis. (e)–(h) Change relative to their initial value of the (e),(f) interface center $x_*$ and (g),(h) relative width $\lambda/\lambda_0$ obtained when fitting the profiles to Eq. (4). Here, the initial interface width $\lambda_0$ has values $\lambda_0 = 0.00428$ and $0.00414$ in (g) and (h), respectively. The dashed line in (e) is a linear fit, while the dashed line in (g) represents diffusive growth of approximately $t_i^2$, showing that $\lambda$ is growing superdiffusively. The dotted curve in (c) and vertical lines in (e) and (g) correspond to time $t = 12000$, above which we judge the left and right boundaries of the slab to interact with one another at the $\tilde{x} = \pm 1$ periodic boundary of the system.</figcaption></figure>

$t \approx 12000$ at which the opposite sides of the slab meet at the periodic boundary. This indicates significant interfacial roughening, faster even than the naïve diffusive expectation for an unsuppressed, freely diffusing interface.

In summary, the glass behaves very differently from the crystal upon equilibration. All four tests reveal discrepancies between observations and expectations based on a nucleation and growth picture.

## V. PLAQUETTE MODELS

Here we show that the very same findings that are inconsistent with a nucleation and growth picture are explained naturally with dynamical facilitation, understood

in its broadest sense as the self-propagation of mobility, independent of the mechanisms through which this propagation is achieved. To this end, we apply our four tests to the TPM and SPPM. These models have been shown to display a diverging (at zero temperature) point-to-set length, but the interfacial tension between amorphous metastable states in these systems is zero, and their equilibrium dynamics is completely driven by dynamical facilitation, not by RFOT dynamics [68].

### A. Model details

The TPM and SPPM are, respectively, two- and three-dimensional lattice models comprising spins $s_i \in \{-1, 1\}$ [i being the lattice coordinate $(i, j)$ or $(i, j, k)$] interacting via Hamiltonians

$$H = -\frac{J}{2} \sum_{i} p_i, \qquad (5)$$

where $J$ is the interaction energy, and the plaquette $p_i \in \{-1, 1\}$ is the product of the triplet of spins $s_{i-1,j}$, $s_{i-1,j+1}$, and $s_{i,j}$ in the TPM and the quintuplet of spins $s_{i-1,j-1,k+1}$, $s_{i-1,j,k+1}$, $s_{i,j-1,k+1}$, $s_{i,j,k+1}$, and $s_{i,j,k}$ in the SPPM. We evolve these systems according to Glauber dynamics [69]: Nodes are randomly selected, one at a time, to undergo a trial flip $s_i \mapsto -s_i$ with probability

$$P(\text{flip}) = \frac{1}{1 + e^{\beta \Delta E}} \qquad (6)$$

of acceptance, where $\Delta E$ is the change in energy that would result from this spin flip and $\beta = 1/k_B T$ for temperature $T$. Given periodic boundary conditions and a number of particles along at least one axis that is a power of 2, there is a one-to-one correspondence between spins and plaquettes [50], and it is in the plaquette representation that the TPM and SPPM are seen to be KCMs, with noninteracting Hamiltonians [Eq. (5)] and facilitation resulting from the kinetic constraint that trial moves correspond to a triplet of plaquettes $p_{i,j}$, $p_{i+1,j-1}$, and $p_{i+1,j}$ flipping simultaneously in the TPM and a quintuplet of simultaneously flipping plaquettes $p_{i,j,k}$, $p_{i,j,k-1}$, $p_{i,j+1,k-1}$, $p_{i+1,j,k-1}$, and $p_{i+1,j+1,k-1}$ in the SPPM.

We take the indices $i, j$, and $k$ to vary along the $x$, $y$, and $z$ axes, respectively, with bond length $b$ separating adjacent nodes along an axis, and adopt $b, J$, and $J/k_B$ as our length, energy, and temperature units. We assume a constant rate of $N$ trial flips per unit of time, where $N$ is the number of lattice points in the system. For both the TPM and SPPM, we simulate systems with 128 lattice points along each axis.

We once again conduct homogeneous and slab simulations. For the slab simulations, we define our slab to be the set of points with axis coordinate $\alpha$ satisfying $|\alpha| < 12$, where $\alpha$ is the $x$ coordinate in the TPM case and the $z$
coordinate in the SPPM case. Since plaquettes do not interact in these systems [Eq. (5)], we equilibrate the slab and nonslab regions by randomly sampling the excitation state of each plaquette site with equilibrium probability $1/(1 + e^{1/T})$ of being excited, where $T = T_{eq}$ in the slab and $T_0$ outside the slab. Because the TPM is not invariant under $x \mapsto -x$ and the SPPM is not invariant under $z \mapsto -z$, we treat the interfaces on the left ($\alpha < 0$) and right ($\alpha > 0$) sides of the slab separately.

## B. Test results

Test 1: Here we look for growing domains in an individual realization of the homogeneous simulation. As we see in Fig. 5, there is domain growth in both the TPM and SPPM during heating but not during cooling. When present, the domains grow with a characteristic triangular shape set by the kinetic constraint, highlighting the fact that this domain growth is due to dynamical facilitation. This asymmetry of having domain growth during heating but not

<figure><img src="image_7.png" /><figcaption>FIG. 5. Structural evolution of the TPM and SPPM during heating and cooling. Systems are initialized at equilibrium at temperature $T_0 = 0.0$ for heating (a),(b) and 1.0 for cooling (c),(d), and the system is then evolved with Glauber dynamics at temperature $T_{eq} = 0.3$ (a), 0.5 (b), or 0.2 (c),(d). Markers correspond to the locations of excited plaquettes in the given snapshot of the system. In (b), we show a slice along the y plane containing the largest number of excited plaquettes at time $t = 15\,000$, so as to capture the birth of the domain. In (d), we instead show a slice along the y plane containing the smallest number of excited plaquettes at time $t = 10^6$, such that if a low-excitation domain were present at this time, we would capture it.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Dynamics of (left columns) heating and (right columns) cooling the TPM and SPPM from temperature $T_0$ to $T_{eq}$ for ($T_0$, $T_{eq}$) = (0.0, 0.3) (first column), (0.0,0.5) (second column), (0.3,0.2) (third column), and (0.5,0.4) (last column), given an artificial target-state slab at time $t=0$. (a)–(d) Initial configuration for (a) TPM heating, (b) SPPM heating, (c) TPM cooling, and (d) SPPM cooling, with markers corresponding to excited plaquettes (in the $y=0$ plane for the SPPM). The hotter region is shaded red, while the colder region is shaded blue. (e)–(h) Normalized concentration $\bar{c} = \frac{1}{2} + \{[c - (c_{eq} + c_0)/2]/|c_{eq} - c_0|\}$ of excited plaquettes within bins of constant $x$ (TPM) or $z$ (SPPM) at times (e) $t=0, 10\,000, \dots, 100\,000$, (f) $t=0, 1000, \dots, 10\,000$, (g) $t=0, 10^3, 10^4, 10^5$, and $10^6$, and (h) $t=0, 10^1, 10^2, 10^3, 10^4$, and $10^5$ (later time, lighter color). Here, $c$ is the unnormalized concentration, and $c_0$ and $c_{eq}$ are the equilibrium concentrations $1/(1+e^{1/T})$ for $T=T_0$ and $T_{eq}$, respectively. Transparent curves show raw data; solid curves are fits to Eq. (4). The solid curves fully overlap with the transparent curves, apart from in (g). Bottom rows: change in the fit parameters (i)–(l) $x_*$ and (m)–(p) $\lambda$ relative to their values in the fit to the $t=0$ profile. The orange curves correspond to the slab interface initially at negative $x$, the blue curves to the slab interface initially at positive $x$. The dashed lines in (i)–(l) are linear fits; the dashed lines in (m)–(p) are power-law fits. There is no blue curve in (o) because it is below the noise floor. The dotted curve in (e) and vertical lines in (i) and (m) correspond to time $t=50\,000$, above which we judge the left and right boundaries of the slab to interact with one another.</figcaption></figure>

during cooling is consistent with the results of Sec. IV for the polydisperse glass, and different from the results of Sec. III for a genuine first-order phase transition in the monodisperse system.

Test 2: Here we look for evidence of a nucleation time. Ironically, given that there is by design no nucleation in these systems (the plaquettes do not interact, so there can be no interfacial tension), the high-energy barrier to thermal activation in an excitation-free system means that we do see a significant waiting time of order $t \sim 10^4$ before the emergence of domains, which grow to envelop the system on

a time of order $t \sim 10^3$. This waiting time arises due to the need to escape the large potential energy basin of the defect-free initial system in order for domain growth to begin.

Two-state scenario: The TPM and SPPM comprise noninteracting plaquettes in one of two states, so we can safely assume that the results of slab simulations will relate to the behavior seen in the homogeneous simulation. There is therefore no need to test the two-state scenario as we did for the glass and crystal.

Test 3: Here we track domain growth quantitatively in the slab simulation. As we show in Figs. 6(e), 6(f), 6(i),
and 6(j), the result of the third test is that slabs in the TPM and SPPM indeed grow upon heating, with the displacement of the interface center from the slab center increasing at a constant rate. The speed at which the center of the left interface advances is different from the corresponding speed for the right interface, as expected given the asymmetric kinetic constraints. By contrast, we see in Figs. 6(g), 6(h), 6(k), and 6(l) that $T_{eq}$-structure slabs in the TPM and SPPM do not grow during cooling, but sometimes even shrink slightly. These results are, once again, consistent with the results of the polydisperse glass.

Test 4: Here we look for evidence of constraints on the roughness of growing domains. During heating, the growing TPM domains in Fig. 5(a) have rough boundaries, but we see little roughness in the growing SPPM domains in Fig. 5(b). There is also no evidence of accelerated domain growth once the system spans the box width, as seen in Figs. 1(a) and 1(e) for the monodisperse crystal. This is confirmed by the growth of the interface width $\lambda$ in the slab simulations. As shown in Fig. 6(m), $\lambda$ grows approximately diffusively during heating in the TPM, while, as seen in Fig. 6(n), the growth of $\lambda$ in the heated SPPM is strongly subdiffusive. Despite the similar designs of the TPM and SPPM, the roughness of growing domains during heating has very different dynamics in the two cases. We note that even within a given model, the power law for the growth of $\lambda$ for the left interface differs from the power law for the right interface. This shows that the asymmetric constraint is affecting the rate of roughening, and highlights the sensitivity of roughening to the mechanism underlying facilitation. The fact that neither the TPM nor the SPPM see the superdiffusive growth of $\lambda$ seen in the heated polydisperse glass is consistent with this.

We see that the main phenomena associated with the equilibration dynamics of glasses—Avrami-like domain growth during heating and an absence of domains during cooling—can be obtained purely via dynamical facilitation, even with trivial thermodynamics, as in the KCMs studied here, suggesting that these two-state scenario behaviors are generic to dynamically facilitated systems. The one phenomenon seen in the glass but not recovered in the KCMs is superdiffusive roughening during heating, but this roughening behavior is highly sensitive to the specific details of the system, varying not only between KCMs, but even between the left and right boundaries of a single slab.

### C. Phenomenological model

From the shared phenomenology of the polydisperse glass and the KCMs, we can build a minimal, phenomenological model of slab growth,

$$ \partial_t c = (c_{\text{eq}} - c)/\tau + s|\partial_x c| + D\partial_x^2 c, \quad (7) $$

where c is the excitation state of the system (e.g., $\Theta$ or the plaquette concentration c), $\tau$ is the local relaxation time of

the system, s is the speed of mobility propagation (the absolute value ensuring that transport is always directed from high to low mobility), and D is a diffusion constant accounting for the roughening of the system. This allows us to discuss the observed dynamics of these systems in terms of specific dynamical ingredients. For simplicity, we focus on slab growth in the TPM, identifying the contributions to the rate of change $\partial_t c$ of the plaquette concentration c, and noting system-dependent considerations when they arise.

The first contribution we consider is the mobility propagation $s|\partial_x c|$. In the ultrastable glass heating literature, experiments [24,28] and a KCM study [51] find front speeds $s = l_s/\tau_{\text{eq}}(T_{\text{eq}})^\gamma$, where $\gamma \approx 1$. Comparing Figs. 4(c) and 6(f) with Figs. 4(i) and 6(j), we see that in both the polydisperse glass and the TPM, the velocity of the interface center remains constant even after significant

<figure><img src="image_9.png" /><figcaption>FIG. 7. Dynamics of the phenomenological model governed by Eq. (7) for constants $l_s = 0.18$ and $l_D = 0.25$ given an artificial target-state slab at time $t = 0$. Heating (left column) is performed between temperatures $T_0 = 0.0$ and $0.3$ using $\tau_0 = 0.549$ and $a = 0.253$, and cooling (right column) is performed between temperatures $T_0 = 0.3$ and $0.2$ using $\tau_0 = 0.705$ and $a = 1.00$. Parameters are chosen so as to yield behavior similar to the TPM, establishing the reasonableness of this model, but no explicit fit to the TPM is performed. (a),(b) Normalized excitation state $\tilde{c}$ defined as in Fig. 6(e) at times (a) $t = 0, 1000, \dots, 10000$ and (b) $t = 0, 10^3, 10^4, 10^5$, and $10^6$ (later time, lighter color). Transparent curves correspond to numerically solved Eq. (7), solid curves to fits to the $x \ge 0$ part of this solution, noting that $c(-x) = c(x)$. (c)–(f) Change, relative to its initial value, of the interface center position $x_*(c),(d)$, and width $\lambda$ (e),(f) obtained when fitting the profiles to Eq. (4). The dashed lines in (c) and (d) are linear fits, while the dashed lines in (e) and (f) represent diffusive growth of approximately $t^{1/2}$.</figcaption></figure>
relaxation of the nonslab region toward equilibrium, so we take $l_s$ to be a constant.

The specific form of the relaxation time $\tau$ is system dependent, but any sufficiently fragile system would yield the speeds $s \approx 0$ for low $T_{eq}$ and $s > 0$ for high $T_{eq}$ seen in Figs. 4 and 6. From the equilibrium behavior $\tau \sim e^{1/\log(3)T^2}$ [68] of the TPM, we obtain

$$ \tau = \tau_0 \exp \left( \frac{1/\log 3}{[(1-a)T_{eq} + aT_f]^2} \right) \quad (8) $$

by modifying the Tool, Narayanaswami, and Moynihan model [70] to yield the TPM relaxation in equilibrium and to work in the limit of vanishing fictive temperature $T_f = 1/\log(c^{-1} - 1)$. Here, $\tau_0$ is a constant timescale, and $a$ is a fit parameter (see Appendix B 4).

Finally, the rate of roughening is system dependent, as we have discussed. In the TPM, this roughening yields an approximately diffusive growth of the slab width $\lambda \sim t^{1/2}$, which can be obtained using a diffusive term $D\partial_x^2 c$. From Figs. 4 and 6, we see that when the interface speed is small, so is the diffusion rate, so we assume proportionality between the two: $D = l_D^2/\tau_{eq}$, where $l_D$ is a constant.

Comparing Figs. 4 and 7, we see that this phenomenological model of relaxation dynamics based on dynamical facilitation captures the key features of slab evolution, both during heating and during cooling, justifying the separation of dynamics into distinct local relaxation, constant-speed drive, and roughening pieces.

## VI. DISCUSSION

We have investigated whether the relaxation of a glass toward equilibrium during heating and cooling proceeds via nucleation and domain growth in the sense of classical nucleation theory. By applying the same four tests to the heating and cooling of a three-dimensional monodisperse crystal, a two-dimensional polydisperse glass, and the TPM and SPPM, we found that the equilibration dynamics of our glass does not resemble that of the melting or freezing crystal, but instead resembles that of the dynamically facilitated plaquette models in every discriminating aspect.

This is made immediately apparent in Fig. 8, which shows the evolution of $1-f$ calculated from Eq. (3) applied to $\Theta$ [Figs. 8(a)–8(d)] or the plaquette concentration $c$ [Figs. 8(e)–8(h)] for ten individual homogeneous simulation trajectories (cf. the ensemble-averaged data in Figs. 1, 3, and 5). For crystal melting and freezing in our monodisperse system, we see a nearly discontinuous change of $1-f$ in Figs. 8(a) and 8(b) once a nucleus is formed, though freezing trajectories are seen to get trapped in a number of metastable states in Fig. 8(b). By contrast, the evolution of $1-f$ in both our polydisperse glass and the two plaquette models [Figs. 8(c)–8(h)] is much smoother and more gradual, reflecting the absence of a nucleation time in these systems.

<figure><img src="image_10.png" /><figcaption>FIG. 8. (a)–(d) Evolution of $1-f$ in homogeneous simulations calculated from Eq. (3) using the mean value of $\Theta$ for ten individual realizations contributing to the averaged results shown in (a) Fig. 1(d), (b) Fig. 1(h), (c) Fig. 3(d), and (d) Fig. 3(h). (e)–(h) Evolution of $1-f$ calculated from Eq. (3) using the concentration $c$ of excited plaquettes for ten individual realizations in the same regimes as those shown in Fig. 5.</figcaption></figure>

For cooling [Figs. 8(d), 8(f), and 8(h)], the lack of variation between individual trajectories reflects the lack of domain growth in these systems, which results in better averaging within individual realizations.

Figure 8 highlights the fact that while the absence of a nucleation time and the asymmetry between heating being driven by domain growth and cooling being a homogeneous process is not consistent with a classical nucleation theory scenario, the equilibration dynamics of our glass is in fact generic to dynamically facilitated systems. Individual trajectories of the glass resemble those of the plaquette models and not the crystal because dynamical facilitation, not nucleation and domain growth, governs the equilibration dynamics of glasses.

### A. Nucleation and domain growth

While we do find that the equilibration of our glass during heating is mediated by domain growth (Figs. 3, 4,
and 8), we do not find evidence of nucleation for either heating or cooling. We find no evidence of a nucleation time for the emergence of $T_{\text{high}}$ domains during heating, which appear on a similar timescale to the timescale over which the front of a growing $T_{\text{high}}$ domain advances by one length unit. We also find no evidence of any thermodynamic drive to minimize the total interface between the domains of $T_{\text{low}}$ structure and $T_{\text{high}}$ structure during this heating process. Unlike in the crystal case [Figs. 1(a) and 1(e)], there is no alignment of the growing $T_{\text{high}}$ domain in Fig. 3(a) along the horizontal axis to minimize its size, and in slab simulations, the slab boundary exhibits superdiffusive roughening. Cooling ultimately provides the most striking evidence against a drive to minimize the interface between domains of $T_{\text{high}}$ and $T_{\text{low}}$ structure, since the absence of domain growth combined with the validity of the two-state equilibration scenario implies that small domains of $T_{\text{low}}$ structure surrounded by $T_{\text{high}}$ structure proliferate during the cooling process.

This result concerns only nucleation in the thermodynamic sense of classical nucleation theory. First-order transitions in systems coupled to a reference configuration 3]] and dynamical phase transitions [71] are found in KCMs, and an explanation of relaxation dynamics based on those transitions may be possible [4,55]. As we have shown, however, this explanation cannot involve traditional nucleation and growth, and we believe that it would ultimately boil down to a complementary view on how dynamical facilitation influences relaxation dynamics.

Note, moreover, that our results also do not contradict the existence of a finite point-to-set length $\xi_{\text{PTS}}$ or the possibility of an RFOT mosaic. In fact, the asymmetry between heating and cooling, including the decoupling of structure and mobility during the latter process, can be derived via the microscopic RFOT theory [72,73]. As shown in Appendix B 3, individual rearrangements perturb the surrounding structure over a length scale $\xi_{\text{perturb}} \sim 10^0$ similar to the point-to-set length $\xi_{\text{PTS}} \sim 10^0$ for our temperature range [35], with the typical size of perturbations to $\Theta$ comparable to the size of $\Theta$ itself. This raises the possibility that individual rearrangements create individual domains of $T_{\text{low}}$ structure within a mosaic as the system cools.

## B. Dynamical facilitation

We have established that dynamical facilitation drives the domain growth observed when heating the polydisperse glass. Despite very different mechanisms for dynamical facilitation, the interface center $x_*$ evolves in the same manner in both the polydisperse glass and our two plaquette models. In all cases, its behavior is well captured by a constant-speed driving term $s|\partial_x c|$ in Eq. (7).

In recent work, Herrero *et al.* [8] conducted very similar slab simulations to ours using the same polydisperse

system and found, suggestively, that it is individual rearrangements near the slab boundary that advance the transformation front. In light of the two-state equilibration scenario, this explains why the significant relaxation of the nonslab region at times $t \ge 10^3$ in the slab heating simulations [Fig. 4(e)] does not affect the constant front speed. As this region relaxes, domains of $T_{\text{high}}$ structure grow which contribute to raising the height of the $T_{\text{low-}}$ region plateau in Fig. 4(c) [parameter $\Theta_0$ in Eq. (4)], but which do not advance the front. Instead, the front advances via the transformation of $T_{\text{low}}$ structure adjacent to the front facilitated by rearrangements near the interface. The same scenario explains observations of the same behavior in the TPM [Fig. 6(i)].

The interface width $\lambda$, for its part, provides system-specific information. The initial interface widths $\lambda_0$ relate to the ability of local structure to constrain structure nearby, hence, $\lambda_0$ in the crystal (Fig. 2) being larger than in the glass (Fig. 4), while $\lambda_0 \approx 0$ in the noninteracting plaquette systems (Fig. 4). Given the two-state equilibration scenario, the evolution of roughness is captured by the evolution of $\lambda$ in the ensemble- and spatially averaged profile $\Theta(\bar{x})$ [Eq. (4)] in slab simulations. As we have seen, this evolution is strongly system dependent, to the point where the widths of the interfaces on the left and right sides of the TPM and SPPM slabs grow as power laws $\lambda \sim t^\nu$ with different power-law exponents $\nu$. Among our glass and two plaquette models, we find exponents $\nu$ spanning values well above [Fig. 3(g)] to well below [Fig. 5(n)] the naïve diffusive exponent $\nu = 1/2$ for roughness arising from random fluctuations in the local front speed.

It was shown in Ref. [8] that slab roughness in the slab-simulation heating of very cold systems at very long times reaches a steady state with characteristic length scale set by the dynamic correlation length $\xi_4$ [74]. Future work may similarly uncover the physical origin of the power-law exponent $\nu$.

## C. Technical advances

Our investigations of nucleation and domain growth and dynamical facilitation are aided by two key technical advances with implications of their own. The first is the use of Eqs. (2) and (3) to confirm the two-state equilibration scenario. While, on the basis of heat capacity measurements, this two-state scenario has long been expected to hold for the heating of ultrastable glasses [6, 14, 22, 24, 28, 29, 57, 58], our work here contains the most direct confirmation of the two-state scenario to date.

Remarkably, we found that the two-state scenario holds for cooling as well as heating, despite the lack of domain growth in the former. This has encouraging implications for the modeling of glasses out of equilibrium, showing that local structure after a rearrangement can simply be sampled directly and independently from the equilibrium
distribution of structure at the target equilibrium temperature, as, for instance, assumed in trap models [75].

A second technical advance of our work is the introduction of the phenomenological model of slab growth, Eq. (7). This model shows how slab growth during heating can be reconciled with homogeneous relaxation during cooling given the same set of dynamical ingredients. By identifying the distinct components of the dynamics of equilibrating glasses, this model should also aid the development of future models for the equilibration of glasses in nonslab geometries.

Data sets and code for the work in this paper can be found at Ref. [76].

## ACKNOWLEDGMENTS

The authors thank C. Scalliet for sharing datasets from Ref. [46] with them. R. N. C. is grateful to L. Berthier, H.-H. Boltz, R. C. Dennis, Z. Fakhraai, M. A. Galvani Cunha, C. Herrero, P. Luo, P. D. Olmsted, M. Ozawa, and A. G. Yodh for helpful discussions concerning this project. The design of Figs. 1 and 3 in this paper was heavily inspired by the figures in Ref. [77]. This work was supported by the Simons Foundation via the “Cracking the glass problem” collaboration (Grant No. 454945, R. N. C. and A. J. L.). A. J. L. also thanks the Simons Foundation for support via Grant No. 327939 as well as for the hospitality of the Center for Computational Biology at the Flatiron Institute, and the Isaac Newton Institute for Mathematical Sciences at Cambridge University (EPSRC Grant No. EP/R014601/1) for support and hospitality. This work used the Anvil Supercomputer at the Rosen Center for Advanced Computing [78] through Allocation No. PHY200101 from the Extreme Science and Engineering Discovery Environment [79], which was supported by National Science Foundation Grant No. 1548562. G. B. acknowledges support from the Simons Foundation Grant (No. 454935 Giulio Biroli).

## APPENDIX A: METHODOLOGICAL DETAILS

### 1. Swap potential

For convenience, we note here that in Eq. (1),

$$c_0 = -28\tilde{r}_{\text{cut}}^{-12}. \qquad (A1)$$

$$c_1 = 48\tilde{r}_{\text{cut}}^{-14}. \qquad (A2)$$

$$c_2 = -21\tilde{r}_{\text{cut}}^{-16}, \qquad (A3)$$

where $\tilde{r}_{\text{cut}} = 1.25$. We also note that $A$, $\sigma_{\text{min}}$, and $\sigma_{\text{max}}$ in the probability density function

$$P(\sigma) = \begin{cases} A\sigma^{-3} & \sigma \in [\sigma_{\text{min}}, \sigma_{\text{max}}], \\ 0 & \text{otherwise,} \end{cases} \qquad (A4)$$

for the particle size $\sigma$ can be found by solving

$$2A \tanh\left(\frac{\bar{\sigma}^2(1 + c_{\sigma}^2)}{2A}\right) = \bar{\sigma}^2 \qquad (A5)$$

numerically for $A$ and noting that

$$\frac{\sigma_{\text{min}}}{\sigma_{\text{max}}} = \exp\left(-\frac{\bar{\sigma}^2(1 + c_{\sigma}^2)}{A}\right) \qquad (A6)$$

and

$$\bar{\sigma} = A(\sigma_{\text{min}}^{-1} - \sigma_{\text{max}}^{-1}), \qquad (A7)$$

ultimately yielding

$$A \approx 1.3145271918, \qquad (A8)$$

$$\sigma_{\text{min}} \approx 0.7244461244, \qquad (A9)$$

$$\sigma_{\text{max}} \approx 1.6138530488 \qquad (A10)$$

for $\bar{\sigma} = 1$ and $c_{\sigma} = 0.23$.

### 2. Preparation protocols

To homogeneously equilibrate to low temperatures, we conduct Monte Carlo simulations at constant volume and temperature until steady state, randomly selecting particles and conducting trial displacements with components randomly sampled from the interval (-0.05, 0.05). For polydisperse systems, we implement a swap Monte Carlo methodology originally developed by Ozawa for Ref. [35]. With probability 0.2, a trial move corresponds to a trial swap instead of a trial displacement. If a particle with size $\sigma_i$ is selected for a swap move, a second particle with size $\sigma_j$ satisfying $|\sigma_j - \sigma_i| < 0.2$ is randomly selected as its swap partner. We start our Monte Carlo simulations from a configuration with random initial positions, except when initializing a crystal in the monodisperse system, in which case we start from a perfect fcc crystal in order to accelerate the equilibration process.

When preparing polydisperse systems containing an artificial target-state slab, we use swap Monte Carlo to equilibrate the system to the target temperature $T_{\text{eq}}$. We then designate a slab-shaped region corresponding to points with $x$ coordinate satisfying $|x| > a$ for some $a > 0$ to be outside the target-state slab. We choose a new, off-target temperature $T_0$ for this region and find the number density for a homogeneous system at equilibrium at temperature $T_0$ and with virial pressure matching that of the target equilibrium system. We stretch or compress the off-slab region along the $x$ axis such that it matches this number density, thus reducing the pressure gradient that would otherwise exist between the slab and the nonslab region. MD simulations are then run at the temperature and pressure corresponding to the target equilibrium system.
When preparing monodisperse systems containing an artificial target-state slab, we choose crystal and liquid temperatures that are similar in magnitude, such that the pressure gradient is less pronounced. It is not possible to stretch or compress the off-slab region along the x axis to reduce the pressure gradient, as we do in the poly-disperse case, as this would prevent us from guaranteeing system dimensions compatible with a perfect fcc crystal. We first homogeneously equilibrate the system to the target crystal state, then conduct Monte Carlo simulations at high temperature, sampling only from the designated liquid region. Once the liquid region has melted, we conduct Monte Carlo simulations at the target, lower temperature for the liquid region, sampling only from this region. MD simulations in this case are run at constant temperature and a constant pressure chosen such that the virial pressure at equilibrium matches the mean virial pressure in the initial configuration.

### 3. Molecular dynamics

We conduct our MD simulations in LAMMPS [61] (June 23, 2022, Update 2) at constant temperature and pressure implemented via LAMMPS' Nosé-Hoover thermostat (FIX NVT) with damping parameter Tdamp = 1 and Berendsen barostat (FIX PRESS/BERENDSEN) with damping parameter Pdamp = 1. If a system has been initialized to temperature T₀ using Monte Carlo simulations, we find the virial contribution to the pressure at this temperature, and choose a pressure such that this contribution is maintained once the temperature changes to the target temperature Tₐₐ. We do this so as to minimize changes in the system size upon changing the temperature.

### 4. Initial and steady-state distributions

Because the start-up of the thermostat and barostat perturbs the system at early times, we must wait until a time t₀ > 0 before we can measure initial distributions for the homogeneous simulations. We use t₀ = 55 for simulations of the monodisperse system and t₀ = 100 for the polydisperse system. For homogeneous melting simulations in the monodisperse system, equilibrium distributions are obtained from t = 600 data, restricting to realizations with bulk average $\bar{\Theta}(t = 600) > 0.091$. For the homogeneous freezing simulations, “equilibrium” distributions are obtained from t = 10 000 data. This does not correspond to equilibrium at the target temperature Tₐₐ, but rather a collection of metastable states that each realization is trapped in at that time (see Sec. VI). For the polydisperse system, equilibrium distributions at temperature Tₐₐ = Tₐₐₗ are obtained from homogeneous heating simulation data for t ≥ 80 000, while equilibrium distributions at temperature Tₐₐ = Tₐₐₗ are obtained from simulations at that temperature in the microcanonical ensemble.

### 5. Θ

The main structural indicator used in this work is Tong and Tanaka's Θ order parameter [65], which measures how inefficiently packed the structure local to each particle is. For each particle i in the system, we iterate over pairs of nearest neighbors j and k such that j and k are also nearest neighbors of each other. We consider them perfectly efficiently packed if the pair potentials Vᵢⱼ, Vⱼₖ, and Vᵢₖ are all equal. If these particles have corresponding sizes σᵢ, σⱼ, and σₖ, this means their centers form a triangle similar to the triangle with sides of length lᵢⱼ = 1/2(σᵢ + σⱼ), lᵢₖ = 1/2(σᵢ + σₖ), and lⱼₖ = 1/2(σⱼ + σₖ). At i, this triangle makes an angle θᵢⱼₖ such that

$$ \cos \theta_{ijk}^* = \frac{l_{ij}^2 + l_{ik}^2 - l_{jk}^2}{2l_{ij}l_{ik}}. \qquad (A11) $$

Letting θᵢⱼₖ be defined such that

$$ \cos \theta_{ijk} = \frac{r_{ij}^2 + r_{ik}^2 - r_{jk}^2}{2r_{ij}r_{ik}}, \qquad (A12) $$

where rᵢⱼ is the distance separating particles i and j, with rⱼₖ and rᵢₖ defined similarly, we define Θ at i to be

$$ \Theta_i = \langle |\theta_{ijk} - \theta_{ijk}^*| \rangle_{j,k}, \qquad (A13) $$

where the average is over nearest-neighbor pairs j and k that are nearest neighbors of i as well as each other. The difference in angle |θᵢⱼₖ - θᵢⱼₖ| is calculated modulo π.

For the purposes of calculating Θᵢ in our two-dimensional polydisperse system, we consider i and j to be nearest neighbors if they have normalized separation $\tilde{r}_{ij} < \tilde{r}_{NN}(\sigma_i)$, where

$$ \tilde{r}_{NN}(\sigma_i) = 2.025 - 0.515\sigma_i. \qquad (A14) $$

This corresponds to the first minimum after the nearest-neighbor peak of g(\tilde{r}), where g(\tilde{r}) is proportional to the probability of finding a particle at normalized separation \tilde{r} given a particle at the origin with size σ, as shown in Fig. 9. For our three-dimensional monodisperse system, the corresponding quantity is \tilde{r} = 1.3. When a particle i does not have any pair of nearest neighbors j and k that are nearest neighbors of each other, we leave its Θᵢ undefined and exclude it when calculating averaged quantities.

### 6. Mermin-Wagner fluctuations

Our two-dimensional glass is influenced by Mermin-Wagner fluctuations [80–82]. We account for this by subtracting from any displacement vector Δpᵢ of a particle i the mean displacement vector ⟨Δpⱼ⟩ⱼ averaged over neighboring particles j ≠ i within a distance ξ of particle
<figure><img src="image_11.png" /><figcaption>FIG. 9. (a) Pair distribution $g(\tilde{r})$ for the polydisperse system at temperature $T = 0.14$ calculated separately for particles at the origin of size $\sigma$ within each of ten logarithmically spaced bins plotted as separate curves (larger $\sigma$, darker curve). This pair distribution is calculated using the same normalization as the usual pair distribution $g(r)$, but replacing the raw interparticle separation $r$ with $\tilde{r}$, hence why $g(\tilde{r}) \to 1$ as $\tilde{r} \to \infty$. (b) $\tilde{r}$ position $\tilde{r}_{NN}$ of the trough following the first-nearest-neighbor peak of $g(\tilde{r})$ as a function of the size $\sigma$ of the particle at the origin for temperatures $T = 0.06$ (blue crosses) and $T = 0.14$ (orange circles). The dashed line is the linear fit yielding Eq. (A14).</figcaption></figure>

i. We choose the value of $\xi$ to be that appropriate to equilibrium at temperature $T_{eq}$ determined as in Ref. [11].

### 7. fcc ground state

To determine that the ground state of the monodisperse system is an fcc crystal, we reason that the ground state of an isotropic repulsive pair interaction decreasing monotonically with distance must correspond to the distance-maximizing fcc or hcp (hexagonal close packing) structures. As shown in Fig. 10, hcp packings at temperature $T < 2.7$ (thermostat implemented using Monte Carlo trial steps at constant volume) transiently increase in potential energy before relaxing into what we find to be an fcc crystal, implying that the fcc state is more stable. (For $T \ge 2.7$, the packings melt into a liquid instead.)

### 8. Interpolating distributions

In the main text, we calculate the transformed fraction $f$ from the evolving distribution $P(\phi, t)$ for different quantities $\phi$ characterizing the local structure or mobility of

<figure><img src="image_12.png" /><figcaption>FIG. 10. Evolution in the canonical ensemble of the mean potential energy of an initially perfect hcp crystal with 22 particles along each basis vector. Curves correspond to temperatures $T = 1.5, 1.6, ..., 5.0$ (lighter color, higher temperature). Time is in units of number of trial Monte Carlo steps per particle.</figcaption></figure>

the system, plotting them as dashed curves in Figs. 1(d), 1(h), 3(d), and 3(h). We achieve this using one of two approaches, depending on the behavior of the growing $P_{eq}(\phi)$ mode. When fluctuations in the distribution near the $P_{eq}(\phi)$ mode are not too severe, we use

$$f(t) = \frac{P(\phi_{*}, t) - P_{0}(\phi_{*})}{P_{\text{eq}}(\phi_{*}) - P_{0}(\phi_{*})}, \quad (A15)$$

where $\phi_{*} = \underset{\phi}{\operatorname{argmax}} P_{\text{eq}}(\phi)$. This is the case for $1 - f$ calculated from $P(\log_{10} \Theta)$ and $P(\log_{10} \mu)$ in Fig. 1(d) (dashed red curve and dashed green curve, respectively), and for $1 - f$ calculated from $P(\Theta)$ in Figs. 3(d) and 3(h) (dashed red curves). However, where the growing $P_{\text{eq}}(\phi)$ mode is sensitive to noise, or where it shifts slightly as it grows [Fig. 3(c)], we use

$$f(t) = \frac{\max_{\phi \in I} P(\phi, t) - \max_{\phi \in I} P_0(\phi)}{\max_{\phi} P_{\text{eq}}(\phi_{*}) - \max_{\phi \in I} P_0(\phi)}, \quad (A16)$$

which is insensitive to these fluctuations. Here, $I$ is an appropriately chosen interval. This approach is used to obtain the red and green dashed curves in Fig. 1(d) [$I = (-\infty, -1.75)$ (red) and $(-\infty, -2)$ (green)], as well as the green dashed curve in Fig. 3(d) [$I = (-0.5, \infty)$].

We also calculate $f$ from Eq. (3) using the averages $\bar{\phi}$ of various quantities $\phi$ and plot them as solid or dotted curves in Figs. 1(d), 1(h), 3(d), and 3(h). To demonstrate that our different methods of calculating $f$, which collapse onto one another, successfully yield the interpolation parameter $f$ in Eq. (2), we compare in Fig. 11 the evolutions $P(\phi, t)$ shown in Figs. 1(d), 1(h), 3(d), and 3(h) to the evolution obtained from Eq. (2) where $P_0(\phi)$ and $P_{\text{eq}}(\phi)$ are taken as input from simulation data and $f$ is obtained from Eq. (3) using $\phi = \log_{10} \Theta$ (for the crystal) or $\Theta$ (for the glass).

## APPENDIX B: SUPPORTING DETAILS

### 1. Perturbed crystal

In Secs. III and IV, we considered systems well described by the simple interpolation between initial and target states expressed in Eq. (2). We stressed that this implies a clean two-state equilibration scenario, as the systems evolve directly from the initial state into the final state. A good test of this idea is provided by systems in which the initial state is perturbed during the heating process, violating Eq. (2). In Fig. 3, for instance, there is a visible rightward shift of the $P_0(\log_{10} \mu)$ peak in Fig. 3(c) as target-state domains grow and a corresponding deviation of the decay of the mean of $\log_{10} \mu$ from that of the fraction $f$ of untransformed structure in Fig. 3(d).

Here, we show in Fig. 12 that slightly increasing the temperature at which we heat our crystal ($T_{\text{high}} = 1.40$, as
<figure><img src="image_13.png" /><figcaption>FIG. 11. Validation of Eq. (2). Each row corresponds to data from the main text: (a),(b), Fig. 1(b); (c),(d) Fig. 1(f); (e),(f) Fig. 3(b); (g),(h) Fig. 3(f). Left-hand side: same data as in the corresponding figure from the main text. Right-hand side: output of Eq. (2) [where $\phi = \log_{10} \Theta$ (a)–(d) or $\Theta$ (e)–(h)] with $f$ calculated from Eq. (3) applied to the quantity $\Theta$.</figcaption></figure>

compared to $T_{high} = 1.35$ for melting in Fig. 1) leads to large rightward shifts of the $P_0(\log_{10} \Theta)$ and $P_0(\log_{10} \mu)$ peaks as the melting proceeds [Figs. 12(a) and 12(b)]. This is accompanied by deviations of $1-f$ calculated from Eq. (3) $\overline{\log_{10} \Theta}$ and $\overline{\log_{10} \mu}$ from the correct value of $1-f$ obtained from the heights of the $P_{eq}(\log_{10} \Theta)$ and $P_{eq}(\log_{10} \mu)$ peaks [Fig. 12(c)]. Inspection of individual trajectories (not shown) suggests that this deviation is due to a large defect density, with pairs of defects imposing strain on crystalline domains trapped between them [note the visible interface of length approximately equal to four particle diameters at the edge of the growing domain in Fig. 1(a), top row].

By contrast, we find that increasing $T_{high}$ to 0.17 in the homogeneous heating of the polydisperse glass, far above $T_{onset}$, does not lead to a perturbation of the untransformed structure and resultant violation of the two-state equilibra- tion scenario beyond that seen in Fig. 3 for $\log_{10} \mu$. This is likely the result of the highly localized nature of facilitation at low temperatures, as seen in Fig. 3(d) of Ref. [11], and

<figure><img src="image_14.png" /><figcaption>FIG. 12. Melting transition in a monodisperse system heated from $T_0 = 1.90$ to $T_{eq} = 2.40$ [cf. $T_{eq} = 2.35$ in Figs. 1(a)–1(d)]. Distributions (a) $P(\log_{10} \Theta)$ and (b) $P(\log_{10} \mu)$ at times $t = 0.0, 3.2, \dots, 48.0$ (later time, lighter curve). (c) Untransformed fraction $1-f$ for $f$ calculated from (dashed red) $P(\log_{10} \Theta)$ and (dashed green) $P(\log_{10} \mu)$ plotted against the mean (solid orange) potential energy $V$, (solid red) $\Theta$, (dotted red) $\log_{10} \Theta$, (solid green) $\mu$, and (dotted green) $\log_{10} \mu$ affinely transformed accord- ing to Eq. (3) with $t_0 = 2.0$. On these short timescales, start-up of the barostat leads to erratic behavior of the inherent-state energy, which is therefore not shown.</figcaption></figure>

the fact that domains of $T_{high}$ and $T_{low}$ do not noticeably statically interact, as we establish in the main text.

## 2. Two-state scenario test: Three-state case

Our two-state scenario test allows us to distinguish multistep relaxation processes from the direct transforma- tion of material from initial-state structure into target-state structure. To illustrate this, we consider a system that fully transforms from the initial state into an intermediate state before transforming into the target equilibrium state. In this case, we have

$$ \bar{\phi} = (1 - f_0)\bar{\phi}_0 + (f_0 - f_1)\bar{\phi}_1 + f_1\bar{\phi}_{eq}, \quad (B1) $$

where $f_0$ is the fraction of the system that has transformed from the initial state into the intermediate state, $f_1$ (which vanishes unless $f_0 = 1$) is the fraction of the system that has transformed from the intermediate state into the target state, and $\bar{\phi}_1$ is the mean of $\phi$ in the intermediate state. It is clear from Eq. (B1) that in this case, $1-f$ calculated from Eq. (3) will be different for different quantities $\phi$.

To confirm this, we adapt the homogeneous simulations of Sec. IV, introducing an intermediate temperature
<figure><img src="image_15.png" /><figcaption>FIG. 13. (a) Heating from temperature $T_0 = 0.06$ to $T_1 = 0.12$ until time $t = 2 \times 10^6$ followed by heating at $T_{eq} = 0.14$. (b) Cooling from temperature $T_0 = 0.14$ to $T_1 = 0.12$ until time $t = 2 \times 10^5$ followed by cooling at $T_{eq} = 0.06$. Plots show $1-f$ calculated from potential energy $\bar{V}$ (solid orange), inherent-state potential energy $\bar{V}_{IS}$ (solid blue), $\bar{\Theta}$ (solid red), $\overline{\log_{10} \Theta}$ (dotted red), $\bar{\mu}$ (solid green), and $\overline{\log_{10} \bar{\mu}}$ (dotted green). Inset to (a): enlargement of times during which $1-f$ evolves rapidly.</figcaption></figure>

$T_1 = 0.12$. We evolve the glass from its initial $T_0$ equilibrium state until steady state with a thermostat at temperature $T_1$, then change the thermostat temperature to $T_{eq}$ and evolve the system further. We fix the pressure based on $T_{eq}$, matching the pressures used in the homogeneous simulations of Sec. IV.

As shown in Fig. 13, the curve collapses seen in Figs. 3(d) and 3(h) are not recovered in this multistage relaxation scenario. For both heating and cooling, the structural quantities $V$, $V_{IS}$, $\Theta$, and (to a lesser extent) $\log_{10} \Theta$ evolve differently from one another, both during the initial evolution at the intermediate temperature $T_1$, and after the thermostat is changed to the target temperature $T_{eq}$. This validates calculating $1-f$ from Eq. (3) for different structural quantities $\phi$ as a test of the two-state scenario.

### 3. Structural perturbation due to rearrangements

As the glass evolves during cooling in Fig. 3(e), it becomes possible to distinguish individual rearrangement events via the characteristic quadrupolar displacement fields they induce, as seen in the snapshots at times $t \ge 10^3$. In Fig. 14, we leverage this to identify the effect of such individual rearrangements on the local structure, identifying rearrangements with particles that have mobility $\mu > 0.1$ and have a larger mobility than any other particle within eight units of distance. We see that changes in $\Theta$ near the rearrangement are of the same magnitude, $|\Delta\Theta| \sim 10^{-2}$, as $\Theta$ itself [Figs. 3(b) and 3(f)] up to a distance $r \approx 5$ from the rearrangement, decaying exponentially with decay length $\xi_{perturb} = 2.34$.

<figure><img src="image_16.png" /><figcaption>FIG. 14. Perturbation to the local $\Theta$ field due to rearrangements in the homogeneously cooled glass studied in Figs. 3(e)–3(h) at time $t = 10^6$. (a) Mean absolute change in $\Theta$ between times $10^6$ and $10^6 + 10^2$ for particles a distance $r$ from a rearranging particle, averaged over ten independent realizations. The inset shows these data on a log scale after subtracting the global mean absolute change in $\Theta$. The dashed line corresponds to a fit to an exponential decay with decay length $\xi_{perturb} = 2.43$. (b) Change in $\Theta$ across the same time interval as (a) for the realization shown in Fig. 3(e).</figcaption></figure>

### 4. Local relaxation in the TPM

Appropriate values for parameters $\tau_0$ and $a$ in Eq. (8) are obtained by fitting Eq. (7) for a homogeneous system to the decay of $\tilde{c}(x = \pm 64)$ toward equilibrium in Figs. 6(e) and 6(g). We show these fits in Fig. 15.

### 5. Supplemental Material videos

In the Supplemental Material [83], we include videos for the homogeneous simulations found in the main text, including many more snapshots than could be put in the main text. We show multiple realizations of each simulation, with the realization shown in the main text included as the leftmost column. We pair the snapshots with the corresponding evolution of $1-f$ for the given realization calculated as in Fig. 8. Supplemental Material Video 1 corresponds to Figs. 1(a) and 8(a), Supplemental Material Video 2 corresponds to Figs. 1(e) and 8(b), Supplemental Material Video 3 corresponds to Figs. 3(a) and 8(c), and Supplemental Material Video 4 corresponds to Figs. 3(e)

<figure><img src="image_17.png" /><figcaption>FIG. 15. Fits (dashed lines) to the decay of $\tilde{c}(x = \pm 64)$ toward equilibrium in the TPM during (a) heating [Fig. 6(e)] and (b) cooling [Fig. 6(g)]. In these plots, $\tilde{c}$ is transformed according to Eq. (3) to make the decay clear.</figcaption></figure>
and 8(d). Supplemental Material Videos 5 and 6 are copies of Supplemental Material Videos 3 and 4, respectively, but with the Θ plot replaced with the change ΔΘ in Θ between times t and t + 100.

[1] 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).

[2] I. Douglass and P. Harrowell, *Can a stable glass be superheated? Modelling the kinetic stability of coated glassy films*, J. Chem. Phys. **138**, 12A516 (2013).

[3] L. Berthier and R. L. Jack, *Evidence for a disordered critical point in a glass-forming liquid*, Phys. Rev. Lett. **114**, 205701 (2015).

[4] R. L. Jack and L. Berthier, *The melting of stable glasses is governed by nucleation-and-growth dynamics*, J. Chem. Phys. **144**, 244506 (2016).

[5] U. S. Cubeta and V. Sadtchenko, *Glass softening kinetics in the limit of high heating rates*, J. Chem. Phys. **150**, 094508 (2019).

[6] A. Vila-Costa, J. Ràfols-Ribé, M. González-Silveira, A. F. Lopeandía, L. Abad-Muñoz, and J. Rodríguez-Viejo, *Nucleation and growth of the supercooled liquid phase control glass transition in bulk ultrastable glasses*, Phys. Rev. Lett. **124**, 076002 (2020).

[7] C. Herrero, C. Scalliet, M. D. Ediger, and L. Berthier, *Two-step devitrification of ultrastable glasses*, Proc. Natl. Acad. Sci. U.S.A. **120**, e2220824120 (2023).

[8] C. Herrero, M. D. Ediger, and L. Berthier, *Front propagation in ultrastable glasses is dynamically heterogeneous*, J. Chem. Phys. **159**, 114504 (2023).

[9] J. D. Eshelby, *The determination of the elastic field of an ellipsoidal inclusion, and related problems*, Proc. R. Soc. A **241**, 376 (1957).

[10] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet, *Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids*, Eur. Phys. J. E **15**, 371 (2004).

[11] 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).

[12] M. Fanfoni and M. Tomellini, *The Johnson-Mehl-Avrami-Kohnogorov model: A brief review*, Nuovo Cimento Soc. Ital. Fis. **20D**, 1171 (1998).

[13] S. F. Swallen, K. Traynor, R. J. McMahon, M. D. Ediger, and T. E. Mates, *Stable glass transformation to supercooled liquid via surface-initiated growth front*, Phys. Rev. Lett. **102**, 065503 (2009).

[14] K. L. Kearns, K. R. Whitaker, M. D. Ediger, H. Huth, and C. Schick, *Observation of low heat capacities for vapor-deposited glasses of indomethacin as determined by AC nanocalorimetry*, J. Chem. Phys. **133**, 014702 (2010).

[15] K. L. Kearns, M. D. Ediger, H. Huth, and C. Schick, *One micrometer length scale controls kinetic stability of low-energy glasses*, J. Phys. Chem. Lett. **1**, 388 (2010).

[16] K. J. Dawson, L. Zhu, L. Yu, and M. D. Ediger, *Anisotropic structure and transformation kinetics of vapor-deposited indomethacin glasses*, J. Phys. Chem. B **115**, 455 (2011).

[17] K. R. Whitaker, M. Ahrenberg, C. Schick, and M. D. Ediger, *Vapor-deposited α, α, β-tris-naphthylbenzene glasses with low heat capacity and high kinetic stability*, J. Chem. Phys. **137**, 154502 (2012).

[18] A. Sepúlveda, E. Leon-Gutierrez, M. Gonzalez-Silveira, M. T. Clavaguera-Mora, and J. Rodríguez-Viejo, *Anomalous transformation of vapor-deposited highly stable glasses of toluene into mixed glassy states by annealing above Tg*, J. Phys. Chem. Lett. **3**, 919 (2012).

[19] Z. Chen, A. Sepúlveda, M. D. Ediger, and R. Richert, *Dynamics of glass-forming liquids. XVI. Observation of ultrastable glass transformation via dielectric spectroscopy*, J. Chem. Phys. **138**, 12A519 (2013).

[20] A. Sepúlveda, S. F. Swallen, and M. D. Ediger, *Manipulating the properties of stable organic glasses using kinetic facilitation*, J. Chem. Phys. **138**, 12A517 (2013).

[21] A. Sepúlveda, M. Tylinski, A. Guiseppi-Elie, R. Richert, and M. D. Ediger, *Role of fragility in the formation of highly stable organic glasses*, Phys. Rev. Lett. **113**, 045901 (2014).

[22] C. Rodríguez-Tinoco, M. Gonzalez-Silveira, J. Ràfols-Ribé, A. F. Lopeandía, M. T. Clavaguera-Mora, and J. Rodríguez-Viejo, *Evaluation of growth front velocity in ultrastable glasses of indomethacin over a wide temperature interval*, J. Phys. Chem. B **118**, 10795 (2014).

[23] D. Bhattacharya and V. Sadtchenko, *Enthalpy and high temperature relaxation kinetics of stable vapor-deposited glasses of toluene*, J. Chem. Phys. **141**, 094502 (2014).

[24] C. Rodríguez-Tinoco, M. Gonzalez-Silveira, J. Ràfols-Ribé, A. F. Lopeandía, and J. Rodríguez-Viejo, *Transformation kinetics of vapor-deposited thin film organic glasses: The role of stability and molecular packing anisotropy*, Phys. Chem. Chem. Phys. **17**, 31195 (2015).

[25] S. S. Dalal and M. D. Ediger, *Influence of substrate temperature on the transformation front velocities that determine thermal stability of vapor-deposited glasses*, J. Phys. Chem. B **119**, 3875 (2015).

[26] M. Tylinski, A. Sepúlveda, D. M. Walters, Y. Z. Chua, C. Schick, and M. D. Ediger, *Vapor-deposited glasses of methyl-m-toluate: How uniform is stable glass transformation?, J. Chem. Phys. **143**, 244509 (2015).

[27] D. M. Walters, R. Richert, and M. D. Ediger, *Thermal stability of vapor-deposited stable glasses of an organic semiconductor*, J. Chem. Phys. **142**, 134504 (2015).

[28] J. Ràfols-Ribé, M. Gonzalez-Silveira, C. Rodríguez-Tinoco, and J. Rodríguez-Viejo, *The role of thermodynamic stability in the characteristics of the devitrification front of vapor-deposited glasses of toluene*, Phys. Chem. Chem. Phys. **19**, 11089 (2017).

[29] U. Cubeta, D. Bhattacharya, and V. Sadtchenko, *Communication: Surface-facilitated softening of ordinary and vapor-deposited glasses*, J. Chem. Phys. **147**, 071101 (2017).

[30] C. J. Fullerton and L. Berthier, *Density controls the kinetic stability of ultrastable glasses*, Europhys. Lett. **119**, 36003 (2017).
[31] E. Flenner, L. Berthier, P. Charbonneau, and C. J. Fullerton, *Front-mediated melting of isotropic ultrastable glasses*, Phys. Rev. Lett. **123**, 175501 (2019).

[32] M. Ruiz-Ruiz, A. Vila-Costa, T. Bar, C. Rodríguez-Tinoco, M. Gonzalez-Silveira, J. A. Plaza, J. Alcalá, J. Fraxedas, and J. Rodriguez-Viejo, *Real-time microscopy of the relaxation of a glass*, Nat. Phys. **19**, 1509 (2023).

[33] C. A. Angell, *On the uncertain distinction between fast landscape exploration and second amorphous phase (ideal glass) interpretations of the ultrastable glass phenomenon*, J. Non-Cryst. Solids **407**, 246 (2015).

[34] R. Das, B. P. Bhowmik, A. B. Puthirath, T. N. Narayanan, and S. Karmakar, *Soft pinning: Experimental validation of static correlations in supercooled molecular glass-forming liquids*, Proc. Natl. Acad. Sci. U.S.A. **2**, pgad277 (2023).

[35] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, *Zero-temperature glass transition in two dimensions*, Nat. Commun. **10**, 1508 (2019).

[36] G. Biroli and J. Bouchaud, *The random first-order transition theory of glasses: A critical assessment*, in *Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications*, edited by P. G. Wolynes and V. Lubchenko (Wiley, New York, 2012), pp. 31–113.

[37] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, *Exact theory of dense amorphous hard spheres in high dimension. III. The full replica symmetry breaking solution*, J. Stat. Mech. (2014) P10009.

[38] 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).

[39] G. M. Hocky, L. Berthier, and D. R. Reichman, *Equilibrium ultrastable glasses produced by random pinning*, J. Chem. Phys. **141**, 224503 (2014).

[40] C. Cammarota, *A general approach to systems with randomly pinned particles: Unfolding and clarifying the random pinning glass transition*, Europhys. Lett. **101**, 56001 (2013).

[41] E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, *Three-dimensional direct imaging of structural relaxation near the colloidal glass transition*, Science **287**, 627 (2000).

[42] 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).

[43] 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).

[44] G. Tarjus, *An overview of the theories of the glass transition*, in *Dynamical Heterogeneities in Glasses, Colloids, and Granular Media*, edited by L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos (Oxford University Press, New York, 2011), pp. 39–67.

[45] O. Dauchot, D. J. Durian, and M. van Hecke, *Dynamical heterogeneities in grains and foams*, in *Dynamical Heterogeneities in Glasses, Colloids, and Granular Media*, edited by L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and

W. van Saarloos (Oxford University Press, New York, 2011), pp. 203–228.

[46] C. Scalliet, B. Guiselin, and L. Berthier, *Thirty milliseconds in the life of a supercooled liquid*, Phys. Rev. X **12**, 041028 (2022).

[47] J. P. Garrahan, P. Sollich, and C. Toninelli, *Kinetically constrained models*, in *Dynamical Heterogeneities in Glasses, Colloids, and Granular Media*, edited by L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. von Sarloos (Oxford University Press, New York, 2011), pp. 341–369.

[48] S. Butler and P. Harrowell, *Glassy relaxation at surfaces: The correlation length of cooperative motion in the facilitated kinetic Ising model*, J. Chem. Phys. **95**, 4466 (1991).

[49] R. Gutiérrez and J. P. Garrahan, *Front propagation versus bulk relaxation in the annealing dynamics of a kinetically constrained model of ultrastable glasses*, J. Stat. Mech. (2016) 074005.

[50] J. P. Garrahan and M. E. J. Newman, *Glassiness and constrained dynamics of a short-range nondisordered spin model*, Phys. Rev. E **62**, 7670 (2000).

[51] S. Léonard and P. Harrowell, *Macroscopic facilitation of glassy relaxation kinetics: Ultrastable glass films with frontlike thermal response*, J. Chem. Phys. **133**, 244502 (2010).

[52] D. Kashchiev, *Nucleation: Basic Theory with Applications* (Butterworth Heinemann, Washington, DC, 2000).

[53] J. P. Garrahan, *Glassiness through the emergence of effective dynamical constraints in interacting systems*, J. Phys. Condens. Matter **14**, 1571 (2002).

[54] J. P. Garrahan, *Transition in coupled replicas may not imply a finite-temperature ideal glass transition in glass-forming systems*, Phys. Rev. E **89**, 030301(R) (2014).

[55] R. M. Turner, R. L. Jack, and J. P. Garrahan, *Overlap and activity glass transitions in plaquette spin models with hierarchical dynamics*, Phys. Rev. E **92**, 022115 (2015).

[56] F. Ritort and P. Sollich, *Glassy dynamics of kinetically constrained models*, Adv. Phys. **52**, 219 (2003).

[57] C. Rodríguez-Tinoco, M. Gonzalez-Silveira, J. Ràfols-Ribé, A. Vila-Costa, J. C. Martinez-Garcia, and J. Rodríguez-Viejo, *Surface-bulk interplay in vapor-deposited glasses: Crossover length and the origin of front transformation*, Phys. Rev. Lett. **123**, 155501 (2019).

[58] A. Vila-Costa, M. Gonzalez-Silveira, C. Rodríguez-Tinoco, M. Rodríguez-López, and J. Rodriguez-Viejo, *Emergence of equilibrated liquid regions within the glass*, Nat. Phys. **19**, 114 (2023).

[59] A. Ninarello, L. Berthier, and D. Coslovich, *Models and algorithms for the next generation of glass transition studies*, Phys. Rev. X **7**, 021039 (2017).

[60] S. Sengupta, T. B. Schröder, and S. Sastry, *Density-temperature scaling of the fragility in a model glass-former*, Eur. Phys. J. E **36**, 141 (2013).

[61] S. Plimpton, *Fast parallel algorithms for short-range molecular dynamics*, J. Comput. Phys. **117**, 1 (1995).

[62] Y. Peng, W. Li, T. Still, A. G. Yodh, and Y. Han, *In situ observation of coalescence of nuclei in colloidal crystal-crystal transitions*, Nat. Commun. **14**, 4905 (2023).
[63] R. L. Penn and J. F. Banfield, *Imperfect oriented attachment: Dislocation generation in defect-free nanocrystals*, Science **281**, 969 (1998).

[64] R. L. Penn and J. F. Banfield, *Oriented attachment and growth, twinning, polytypism, and formation of metastable phases; insights from nanocrystalline TiO₂*, Am. Mineral. **83**, 1077 (1998).

[65] H. Tong and H. Tanaka, *Revealing hidden structural order controlling both fast and slow glassy dynamics in supercooled liquids*, Phys. Rev. X **8**, 011041 (2018).

[66] Note that $P_0(\log_{10}\Theta)$ for heating is not the same as $P_{eq}(\log_{10}\Theta)$ for cooling, since a randomly nucleated crystal is unlikely to be oriented in a direction compatible with the system axes, so equilibrium at the perfect fcc crystal structure is not generally attainable upon cooling.

[67] Symmetry-breaking sample-to-sample fluctuations are suppressed upon ensemble averaging in practice, and in the limit of infinite size along the y and z axes in principle.

[68] R. L. Jack and J. P. Garrahan, *Caging and mosaic length scales in plaquette spin models of glasses*, J. Chem. Phys. **123**, 164508 (2005).

[69] R. J. Glauber, *Time-dependent statistics of the Ising model*, J. Math. Phys. (N.Y.) **4**, 294 (1963).

[70] C. T. Moynihan, A. J. Easteal, M. A. De Bolt, and J. Tucker, *Dependence of the fictive temperature of glass on cooling rate*, J. Am. Ceram. Soc. **59**, 12 (1976).

[71] L. O. Hedges, R. L. Jack, J. P. Garrahan, and D. Chandler, *Dynamic order-disorder in atomistic models of structural glass formers*, Science **323**, 1309 (2009).

[72] V. Lubchenko and P. G. Wolynes, *Theory of structural glasses and supercooled liquids*, Annu. Rev. Phys. Chem. **58**, 235 (2007).

[73] P. G. Wolynes, *Spatiotemporal structures in aging and rejuvenating glasses*, Proc. Natl. Acad. Sci. U.S.A. **106**, 1353 (2009).

[74] 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).

[75] C. Monthus and J.-P. Bouchaud, *Models of traps and glass phenomenology*, J. Phys. A **29**, 3847 (1996).

[76] *EquilibrationDynamics_repository*: https://bitbucket.org/chackorn/equilibrationdynamics_repository; R. Chacko, F. Landes, G. Biroli, O. Dauchot, A.Liu, and D. R. Reichman, *Data sets for: Dynamical Facilitation Governs the Equilibration Dynamics of Glasses*. figshare (2024), 10.6084/m9.fishshare.25971439.v1.

[77] A. Chardac, L. A. Hoffmann, Y. Poupart, L. Giomi, and D. Bartolo, *Topology-driven ordering of flocking matter*, Phys. Rev. X **11**, 031069 (2021).

[78] X. Carol Song, Preston Smith, Rajesh Kalyanam, Xiao Zhu, Eric Adams, Kevin Colby, Patrick Finnegan, Erik Gough, Elizabeth Hillery, Rick Irvine, Amiya Maji, and Jason St. John, *Anvil—System Architecture and Experiences from Deployment and Early User Operations. In Practice and Experience in Advanced Research Computing (PEARC '22)* (Association for Computing Machinery, New York, NY, USA, 2022), Article 23, 1–9, 10.1145/3491418.3530766.

[79] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, *XSEDE: Accelerating scientific discovery*, Comput. Sci. Eng. **16**, 62 (2014).

[80] E. Flenner and G. Szamel, *Fundamental differences between glassy dynamics in two and three dimensions*, Nat. Commun. **6**, 7392 (2015).

[81] S. Vivek, C. P. Kelleher, P. M. Chaikin, and E. R. Weeks, *Long-wavelength fluctuations and the glass transition in two dimensions and three dimensions*, Proc. Natl. Acad. Sci. U.S.A. **114**, 1850 (2017).

[82] 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).

[83] See *Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.14.031012* for six videos.