# Two-step devitrification of ultrastable glasses

Cecilia Herrero, Camille Scalliet, M. Ediger, Ludovic Berthier

► To cite this version:

Cecilia Herrero, Camille Scalliet, M. Ediger, Ludovic Berthier. Two-step devitrification of ultrastable glasses. Proceedings of the National Academy of Sciences of the United States of America, 2023, 120 (16), 10.1073/pnas.2220824120. hal-04117704

HAL Id: hal-04117704

https://hal.science/hal-04117704

Submitted on 5 Jun 2023

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
Two-step devitrification of ultrastable glasses

Cecilia Herrero,¹ Camille Scalliet,² M. D. Ediger,³ and Ludovic Berthier¹,⁴

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

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

³Department of Chemistry, University of Wisconsin-Madison, Madison, Wisconsin 53706, USA

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

(Dated: January 5, 2023)

The discovery of ultrastable glasses has raised novel challenges about glassy systems. Recent experiments studied the macroscopic devitrification of ultrastable glasses into liquids upon heating but lacked microscopic resolution. We use molecular dynamics simulations to analyse the kinetics of this transformation. In the most stable systems, devitrification occurs after a very large time, but the liquid emerges in two steps. At short times, we observe the rare nucleation and slow growth of isolated droplets containing a liquid maintained under pressure by the rigidity of the surrounding glass. At large times, pressure is released after the droplets coalesce into large domains, which accelerates devitrification. This two-step process produces pronounced deviations from the classical Avrami kinetics and explains the emergence of a giant lengthscale characterising the devitrification of bulk ultrastable glasses. Our study elucidates the nonequilibrium kinetics of glasses following a large temperature jump, which differs from both equilibrium relaxation and aging dynamics, and will guide future experimental studies.

Glasses play a critical role in modern technology, in applications as diverse as optical fibers and organic light emitting diode (OLED) displays [1, 2]. However, in spite of their solid-like appearance, all glasses spontaneously evolve in a process known as physical aging [3]. Aging is an inherent feature of the glassy state that results from its nonequilibrium nature, and it is an important challenge because the resulting structural evolution can negatively influence material properties. As a result of recent progress, aging is reasonably well understood if the temperature is lowered by only a few Kelvins near equilibrium. In this case, knowledge of equilibrium response functions allows predictions for the aging dynamics [4]. Larger jumps to low temperature are more challenging, but their qualitative features are well understood [3-7].

By contrast, the opposite case of a large positive temperature jump is much less understood, and enhanced understanding could provide important conceptual advances regarding amorphous materials. Qualitatively different phenomena are observed in comparison to the down jumps and near-equilibrium jumps, with structural relaxation at first occurring very slowly and then accelerating [3, 8]. The recent development of ultrastable glasses, accessed via new experimental [9-11] and simulation [12-14] techniques, provides a new tool to understand large up-jump experiments. Because these materials are prepared in low energy states with high kinetic barriers, extreme up-jump experiments can now be performed with a wider variety of experimental techniques.

In the last decade, positive temperature jumps performed from ultrastable states have revealed previously hidden features of amorphous materials. For ultrastable films, the devitrification from the glass to the liquid occurs heterogeneously, starting from the free surface and propagating towards the interior over large distances (> 1 micron). This process has been studied extensively in

experimental [15-19] and theoretical [20-24] work.

Fewer studies have investigated the homogeneous devitrification of ultrastable glasses in the bulk [25-28]. Surprisingly, macroscopic measurements are best interpreted using an analogy with the nucleation and growth kinetics traditionally observed across first-order phase transitions [15, 29]. However, quantitative analysis requires a number of assumptions [15, 28] regarding the underlying physics, which have not been tested. Previous simulations used relatively small systems [30, 31] or simplified models [23, 29, 32]. Therefore, an understanding of microscopic processes, of the characteristic time scales and length scales and their evolution with glass stability is lacking. At the fundamental level, confirming a deep analogy between crystal melting and devitrification would provide experimental evidence of a thermodynamic interpretation of the glass transition [1, 29].

We report results from molecular dynamics simulations which paint a complete microscopic description of the transformation kinetics of ultrastable glasses following a positive temperature jump. We prepare ultrastable configurations at very low initial temperature $T_i$ using the swap Monte Carlo algorithm [13, 33]. We instantaneously heat the system to the annealing temperature $T_a > T_i$ and observe how the system devitrifies from its initial low-enthalphy state to the equilibrium liquid at $T_a$. Our results generally confirm a nucleation and growth kinetics, but they deviate from the classical Avrami description [34] and existing interpretations on several aspects. In particular, the liquid droplets that nucleate initially grow much more slowly than the large domains found at larger times, because they contain a liquid that is compressed by the surrounding glass. As a result, the large length scale characterising the bulk devitrification of ultrastable glasses is a novel physical quantity, distinct from the crossover length scale discussed for thin films.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Transformation of an ultrastable glass into an equilibrium liquid. (a) Visualising the transformation of an ultrastable glass prepared at $T_i = 0.035$ annealed at $T_a = 0.11$ with $N = 256000$ at times $t$ indicated on the figure. The bond-breaking correlation distinguishes the glass with $C_B^i = 1$ (blue) from the liquid with $C_B^i = 0$ (red). Rare isolated liquid droplets nucleate, grow slowly, then coalesce to transform the entire system. (b) Time evolution of number density $\rho$ for $T_i = 0.035$ and different annealing temperatures $T_a$. The initial glass (I) rapidly expands to the glass state (G), and transforms over a much slower time scale into the equilibrium liquid (L). (c) Time evolution of the liquid fraction $x_l(t)$ for the same parameters as (b). The time axis is normalised by the transformation time $t_{tr}$. Dashed lines: Avrami functions $1 - \exp(-Kt^n)$ with different exponents $n$. (d) Evolution of the effective exponent $n$ as a function of $x_l(t)$. The dashed classical Avrami value $n = 3$ does not adequately describe our results.</figcaption></figure>

We perform molecular dynamics (MD) simulations of a two dimensional size-polydisperse mixture of $N = 64000$ soft repulsive particles, known to be a reliable glass-former [35]. We prepare ultrastable configurations at number density $\rho_i = 1$, temperature $T_i = 0.035$ and pressure $P_i = 2.31$ using the swap Monte Carlo algorithm [13, 33] and reduced numerical units. We then perform classical MD simulations in the NPT ensemble at pressure $P_i$ for a broad range of annealing temperatures $T_a > T_i$. More information about the model, simulations and reduced units are provided in the Methods.

We start in Fig. 1a with a global description of the transformation process for the particular case $T_a = 0.11$. To visualise the transformed regions, we introduce a local structural correlation function $C_B^i(t)$ which records the fraction of nearest neighbors lost by particle $i$ since $t = 0$ [36]. Below, we count a particle as ‘liquid’ whenever $C_B^i(t) \le 0.5$, but we checked equivalence with several other definitions. We observe that devitrification starts from rare regions distributed throughout the system. As time increases, the size of these regions grows and new regions keep appearing. This is consistent with nucleation and growth, but distinct from the pre-existing nuclei picture proposed recently [28]. At larger times, the

growing droplets coalesce and then percolate throughout the system. At very late times, a few glass regions survive, which are eventually invaded by the surrounding liquid. We define the characteristic times $t_{tr}$ and $t_{1/2}$ when respectively 99% and 50% of the particles have become liquid.

The time evolution of the density is shown in Fig. 1b for several annealing temperatures $T_a$. By construction, all curves start from the initial density $\rho_i = 1$, state (I). As a result of the temperature jump from $T_i$ to $T_a$, the system quickly expands without structural rearrangement to reach the glass state (G) in Fig. 1b. Its density is smaller, $\rho_G < \rho_i$. At much larger times, the glass transforms into the final liquid (L) which has an even lower density $\rho_L < \rho_G$. Both $\rho_G$ and $\rho_L$ depend on $T_a$. As found in smaller temperature jumps [3, 8], the dynamics is initially very slow but the glass transformation suddenly accelerates and takes place abruptly. The transformation time $t_{tr}$ increases by 5 orders of magnitude with decreasing the annealing temperature $T_a$.

The transformation kinetics can be quantified using the fraction $x_l(t)$ of liquid particles at time $t$, which goes from 0 at $t = 0$ to 1 at long times, see Fig. 1c. To compare different $T_a$, we rescale the time axis by $t_{tr}$. As $T_a$ de-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Understanding the early stages of the transformation by simulating isolated droplets. (a) Snapshot of a single liquid droplet of radius $R_0 = 80$ in a large $N = 576000$ sample. Color code provided in Fig. 1a. (b) Radial displacement field in the glass after mechanical equilibrium is reached for $R_0 = 30$. (c) Liquid and glass equations of state showing initial state (I), the expanded glass (G) and final liquid (L) for $T_i = 0.035$ and $T_a = 0.11$. (D) indicates the isolated droplet density. (d) The time evolution of the droplet radius allows us to measure the growth velocity $v_1$ of isolated droplets containing the compressed liquid.</figcaption></figure>

creases, the time evolution of $x_l(t)$ becomes sharper and eventually collapses onto a master curve for the lowest annealing temperatures.

Following previous studies [15, 23, 28], we fit $x_l(t)$ to the Avrami equation [34]

$$x_l(t) = 1 - \exp(-Kt^n). \quad (1)$$

In the classical Avrami picture, a nucleation rate describes the initiation of droplets, and a constant velocity captures their growth. Both parameters enter the definition of $K$, and the exponent is predicted to be $n = d+1$, with $d$ the spatial dimension. The classical Avrami prediction is thus $n = 3$ for our two-dimensional simulations. We show in Fig. 1c that this exponent does not describe our measurements accurately. To better quantify these deviations, we follow the experimental literature and report the parametric evolution of the effective exponent $n$ as a function of $x_l$ in Fig. 1d. While $n$ is not well-defined at high $T_a$, it consistently takes a surprisingly large value $n \approx 4.5$ for low $T_a$, indicating a breakdown of the classical nucleation and growth picture. In the following we deploy original numerical strategies to explain this behaviour.

In the early stages of the transformation, the system is composed of isolated liquid droplets confined by the surrounding glass matrix. We analyse an idealised version of this geometry with a single liquid droplet of radius $R_0 \gg 1$ immersed in a very large glass matrix of linear size $L \gg R_0$, see Fig. 2a. We used $L \approx 770$ and found identical results for $R_0 = 30, 80, 50$ and $100$. Our protocol, detailed in the Methods, is the following. We start from the glassy state (G) in Fig. 1b and perform high-temperature dynamics in a spherical domain of radius $R_0$ to transform this region, followed by thermalization at $T_a$. This is done while keeping the particles outside the cavity immobile. At this stage, the density in the liquid droplet is thus equal to that of the glass (G), implying that the liquid pressure has increased above $P_i$. Then,

we simulate the dynamics of the whole sample at temperature $T_a$ and constant volume. At the very beginning of the simulation, the compressed liquid exerts pressure on the surrounding glassy particles to expand, but the expansion is opposed by the rigidity of the glass, until a mechanical equilibrium is found. We have computed the radial displacement in the glass right after mechanical equilibrium is reached. Its angular-averaged amplitude $\delta r$, shown in Fig. 2b, is maximal at the edge of the droplet $r/R_0 = 1$ and decays algebraically to zero at large distances. Consequently, the droplet expands from $R_0$ to $R_0 + \Delta R_0$ so that its density decreases by $\Delta \rho \approx 2\rho_G\Delta R/R_0$. The corresponding droplet (D) density $\rho_D = \rho_G - \Delta \rho$ is shown in Fig. 2c together with the other state points. Far from the droplet, the glass pressure is equal to $P_i$, but it is larger in the liquid droplet which contains a liquid at density $\rho_D$ and pressure $P_D > P_i$ which is therefore distinct from the final equilibrium liquid (L) which is instead at conditions $(P_i, T_a)$.

At much longer times, structural relaxation occurs and the liquid droplet grows slowly. We monitor the time evolution of the liquid fraction $x_l(t)$, which is directly connected to the growing radius $R(t)$ since $x_l(t) \propto R^2(t)$. In Fig. 2d we show that the droplet size increases linearly with time, $R(t) = R_0 + v_1t$, with a velocity $v_1$ that can be measured numerically. We have confirmed that the same value $v_1$ describes the growth of isolated droplets randomly selected in the real process illustrated in Fig. 1. This demonstrates that the idealised geometry depicted in Fig. 2 faithfully represents the initial stages of the bulk devitrification.

To better characterize the late times at which very large liquid domains percolate throughout the system, we simulate the idealised geometry shown in Fig. 3a in which we create a single percolating liquid domain [30] of width $2W_0$. The numerical strategy is similar to that for the droplet geometry, see Methods. Once the vertical liquid domain is prepared, we run NPT simulations at
<figure><img src="image_4.png" /><figcaption>FIG. 3. Understanding the late stages of the transformation by simulating a macroscopic interface. (a): Snapshot of a two-front system with $N = 64000$. Color code as in Fig. 1a. (b) The time evolution of the interface position allows us to measure the growth velocity $v_2$ of large fronts.</figcaption></figure>

temperature $T_a$ and pressure $P_i$ for the whole sample. In this case, the pressure is homogeneous across the system and the liquid density is $\rho_L$. As time increases, we observe that the liquid propagates towards the glass as a front. We extract its velocity by monitoring the time evolution of $x_l(t)$, which is directly proportional the displacement of the interface. This is shown in Fig. 3b for $T_a = 0.11$. The front propagation is again ballistic, with $x_l(t) \propto v_2t$ which defines the front velocity $v_2$. Again, we have checked that the same value $v_2$ describes the growth of large interfaces randomly selected in the real transformation process shown in Fig. 1, confirming the relevance of the idealised geometry studied in Fig. 3.

By gathering velocity measurements for different annealing temperatures and the two idealised geometries shown in Figs. 2, 3 we find that the growth of large interfaces is systematically faster than the growth of isolated droplets, $v_1 < v_2$. This is because the liquid inside droplets is under pressure $P_D > P_i$ and its relaxation time is thus larger than the one of the equilibrium liquid (L). For a given $T_i$, we observe that the difference $v_2 - v_1$ increases with increasing $T_a$, because the separation between the points (G) and (L) increases in Fig. 2c. For a given $T_a$, we also find a larger difference in the velocities for lower $T_i$ because the glass becomes stiffer and compresses the liquid droplet more at short times.

The quantitative comparison between the kinetics of isolated droplets and large fronts demonstrates that a central assumption leading to Avrami kinetics does not hold for the bulk devitrification of ultrastable glasses, as there is not a unique velocity describing the entire process. We argue that this more complex two-step process is responsible for the anomalously large value of the exponent $n$ shown in Fig. 1e. First, we note that at early times $x_l \ll 1$ when only isolated droplets exist, the classical value $n = 3$ is observed. However, as the liquid fraction grows and droplets start to coalesce, $n$ grows to a much larger value $n \approx 4.5$. Qualitatively, this implies that the transformation process accelerates further, as expected if the velocity increases from $v_1$ to $v_2$. In

the Supplementary Information, we show that adding the simplest hypothesis of a smooth transition between two values $v_1$ and $v_2 > v_1$ in the analytic Avrami description generically produces values $n \approx 5$, which supports our physical interpretation of the numerical results. Note that the pressure effect, which controls kinetic stability and the devitrification melting, is absent in several simulation studies [23, 29, 30, 32].

Finally, we collect in Fig. 4 the relevant time scales and length scales characterising the devitrification of ultrastable glasses. Starting with times in Fig. 4a, the slowest time scale is the transformation time $t_{tr}$, which is much larger than any time scale characterising the equilibrium relaxation of the liquid in the same conditions. We display the decay time of the bond-breaking correlation $C_B(t)$, $\tau_B$, and of the self-intermediate scattering function, $\tau_\alpha$, defined in the Methods section. The kinetic stability, encoded in the ratio between $t_{tr}$ and these equilibrium relaxations reaches a large value $10^3 - 10^5$ at the lowest $T_a$, and would increase further if longer timescales could be simulated allowing access to even lower $T_a$. This large kinetic stability is comparable to values previously reported in both vapor deposited films [17, 19, 37] and simulations performed using the swap Monte Carlo algorithm [24, 31].

Turning to length scales, we report in Fig. 4b the evolution of the characteristic transformation length $\ell_{tr}$ suggested by Fig. 4c where snapshots of the system at various $T_a$ are shown at time $t_{1/2}$ when half of the glass has transformed. These images reveal that the transformation process is highly heterogeneous in space, with a characteristic length scale $\ell_{tr}$ which grows as kinetic stability increases. To characterize this length, we employ the method using chord length distributions, following the same recipe as in [38, 39]. We define $\ell_{tr}$ as the first moment of the chord distributions. As shown in Fig. 4b, $\ell_{tr}$ grows with decreasing $T_a$ reaching $\ell_{tr} \approx 60$ for $T_a = 0.10$ and a kinetic stability of about $4 \times 10^5$. In these conditions, the transformed domains barely fit in our simulation box of linear size $L \approx 250$, see the rightmost snapshot in Fig. 4c.

In experimental studies of thin films transforming via front propagation, a large 'crossover' length scale $\ell_c$ has been reported [15, 19]. It is defined as the distance travelled by the liquid front over the transformation time of the corresponding bulk system. Using our notations, we have $\ell_c = v_2 t_{tr}$, which can be determined from our independent measurements of $v_2$ and $t_{tr}$. Previous work [15, 19, 28] assumed that the two length scales $\ell_{tr}$ and $\ell_c$ are equal. This is indeed correct when the classical Avrami picture applies [29]. Our results however imply that this picture does not describe the transformation of ultrastable glasses. In Fig. 4b we confirm that the two characteristic length scales are different, and evolve differently with $T_a$. Again, the existence of two distinct velocities rationalises this finding. A smaller value $v_1 < v_2$ in the first stages of devitrification implies that $t_{tr}$, and thus $\ell_c$, are larger than if $v_1 = v_2$. Instead, a longer $t_{tr}$ al-
<figure><img src="image_5.png" /><figcaption>FIG. 4. Characteristic time scales and length scales of ultrastable glass transformation. (a) Temperature dependence of the transformation time $t_{tr}$ and two characteristic times of the equilibrium liquid: $\tau_B$ defined from $C_B$ and $\tau_\alpha$ from $F_s(q, t)$. (b) Characteristic lengths of the process: $\ell_{tr}$ is the characteristic size of the domains at time $t_{1/2}$, and the crossover length $\ell_c = v_2t_{tr}$. (c) Snapshots taken at time $t_{1/2}$ for $N = 64000$ and different annealing temperatures $T_a = 0.20, 0.17, 0.13, 0.10$ (from left to right).</figcaption></figure>

lows more time for droplets to nucleate so that $\ell_{tr}$, which controlled by the typical distance between droplets, decreases. The opposite influence of the velocity contrast $v_2/v_1$ on physical length scales implies that in general $\ell_{tr}$ is distinct from, and smaller than, $\ell_c$.

Our computational study of the devitrification of ultrastable glasses reveals several unexpected results that had not been predicted by previous theoretical models, and have either not been observed experimentally or contradict published results. We now discuss these points.

Our first novel observation is the nucleation at short times of liquid droplets containing a liquid that is under pressure. This stems from the large density difference between the ultrastable glass and the equilibrium liquid, combined to the rigidity of the surrounding glass, see Fig. 2. Although these two effects may quantitatively differ in different materials, it is difficult to imagine a physical relaxation process that would entirely release this pressure. This finding has several important consequences. First, it implies that the liquid that first appears during the transformation has a relaxation dynamics that is much slower than the equilibrium liquid that forms at long times. Published works probing these dynamics have not reported any evidence for such a slow process [17, 27, 28], and future experiments should revisit this point. A related consequence is that the growth

of these compressed liquid droplets is much slower than the growth of the fronts propagating from a free surface. This directly contradicts an important assumption made in Ref. [15, 28] to interpret experimental results. To our knowledge, the growth of isolated droplets has not yet been studied experimentally.

A final consequence is that the growth velocity of the liquid phase is not constant during the transformation, which directly explains the emergence of an anomalously large Avrami exponent describing the transformation kinetics. This finding constitutes our second key novel observation. It contrasts with the recent report of a smaller exponent extracted from calorimetric measurements [28]. These recent results were explained by hypothesizing a population of pre-existing defects initiating the transformation. We have not observed such defects in our simulations where the appearance of new liquid droplets is instead slow and randomly distributed in time. We can only invoke two important differences between simulations and experiments that could explain this discrepancy. We directly determine the liquid fraction $x_l(t)$ in *situ* whereas in experiments it is indirectly inferred from calorimetric measurements involving a complicated thermal treatment and data analysis. Another difference is the devitrification time window analysed in simulations, which is of order 10 ms in simulations, but about $10^3$
times slower in most experiments [28].

Our third key observation is that a large length scale characterizes the transformation process, which becomes larger as the kinetic stability of the glass increases. Despite our numerical constraints, we directly measure a length scale approaching $ℓ_{tr} \approx 100\sigma$ which would presumably grow larger if we could simulate longer time scales. We are not aware of any bulk relaxation process in equilibrium supercooled liquids or in aging glasses that can reach such large values. Qualitatively, the large length scale $ℓ_{tr}$ reflects the equally-large distance between nucleated liquid droplets in the early stages of the transformation. As $T_a$ decreases, droplets become more sparse as their nucleation rate decreases, and $ℓ_{tr}$ increases. Therefore, it is the large barrier to nucleating liquid droplets which is responsible for the emergence of a large length scale. This interpretation explains why such large length scales have never been found in conventional glasses, as these less stable systems would have smaller barriers.

There has been extensive experimental work deducing that a large crossover length scale $ℓ_c$ appears in thin films transformation via the propagation of a front. We have shown that these two length scales are conceptually and quantitatively distinct. In particular, $ℓ_c$ simply represents a crossover between distinct physical regimes during front propagation, but does not characterise any kind of spatial correlations, contrary to $ℓ_{tr}$. The length scale $ℓ_{tr}$ has recently been indirectly inferred from experiments, but this assumes hypotheses which do not hold in our simulations. A more direct experimental measurement of the large length scale $ℓ_{tr}$ is an important target for future work.

More fundamentally, the observation that the transformation

of an ultrastable glass into a liquid proceeds via a physical mechanism traditionally observed for first-order transitions raises important questions. It was argued in Ref. [29] that the existence of a thermodynamic Kauzmann transition provides the needed ingredients to interpret this analogy. In the alternative approach of Ref. [23] based on a kinetically constrained model, the additional introduction of some thermally activated nucleation sites was needed to reproduce the observed phenomenology, with the growth processes controlled by dynamic facilitation. Future work should therefore concentrate on the first stages of the transformation to understand better their physical origin. It would also be very interesting to understand whether the regions that first transform are related to specific structural features of ultrastable glassy states.

## ACKNOWLEDGMENTS

This work was publicly funded through ANR (the French National Research Agency) under the Investissements d'avenir programme with the reference ANR-16-IDEX-0006. It was also supported by a grant from the Simons Foundation (#454933, LB), the European Research Council under the EU's Horizon 2020 Program, Grant No. 740269 (CS), the U.S. National Science Foundation, CHE-2153944 (MDE), and by a Visiting Professorship from the Leverhulme Trust (VP1-2019-029, LB). CS acknowledges support from a Herchel Smith Fellowship, University of Cambridge, and a Ramon Jenkins Research Fellowship from Sidney Sussex College, Cambridge.

[1] L. Berthier and M. D. Ediger, Facets of glass physics, Physics Today 69, 40 (2016).

[2] J. Ràfols-Ribé, P.-A. Will, C. Hänisch, M. Gonzalez-Silveira, S. Lenk, J. Rodríguez-Viejo, and S. Reineke, High-performance organic light-emitting diodes comprising ultrastable glass layers, Science advances 4, eaar8332 (2018).

[3] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, and S. W. Martin, Relaxation in glassforming liquids and amorphous solids, Journal of Applied Physics 88, 3113 (2000).

[4] B. Riechers, L. A. Roed, S. Mehri, T. S. Ingebrigtsen, T. Hecksher, J. C. Dyre, and K. Niss, Predicting nonlinear physical aging of glasses from equilibrium relaxation via the material time, Science Advances 8, eaal9809 (2022).

[5] G. W. Scherer, Theories of relaxation, Journal of Non-Crystalline Solids 123, 75 (1990).

[6] I. Hodge, Enthalpy relaxation and recovery in amorphous materials, Journal of Non-Crystalline Solids 169, 211 (1994).

[7] S. Simon, J.-Y. Park, and G. McKenna, Enthalpy recovery of a glass-forming liquid constrained in a nanoporous matrix: Negative pressure effects, The European Physical Journal E 8, 209 (2002).

[8] A. J. Kovacs, Transition vitreuse dans les polymères amorphes. étude phénoménologique, in Fortschritte der hochpolymeren-forschung (Springer, 1964) pp. 394–507.

[9] S. F. Swallen, K. L. Kearns, M. K. Mapes, Y. S. Kim, R. J. McMahon, M. D. Ediger, T. Wu, L. Yu, and S. Satija, Organic glasses with exceptional thermodynamic and kinetic stability, Science 315, 353 (2007).

[10] M. D. Ediger, Perspective: Highly stable vapor-deposited glasses, The Journal of Chemical Physics 147, 210901 (2017).

[11] C. Rodriguez-Tinoco, M. Gonzalez-Silveira, M. A. Ramos, and J. Rodriguez-Viejo, Ultrastable glasses: new perspectives for an old problem, La Rivista del Nuovo Cimento, 1 (2022).

[12] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Equilibrium sampling of hard spheres up to the jamming density and beyond, Phys. Rev. Lett. 116, 238002 (2016).

[13] A. Ninarello, L. Berthier, and D. Coslovich, Models and algorithms for the next generation of glass transition studies, Phys. Rev. X 7, 021039 (2017).
[14] A. D. S. Parmar, M. Ozawa, and L. Berthier, Ultrastable metallic glasses in silico, Phys. Rev. Lett. 125, 085505 (2020).

[15] K. L. Kearns, M. D. Ediger, H. Huth, and C. Schick, One micrometer length scale controls kinetic stability of low-energy glasses, The Journal of Physical Chemistry Letters 1, 388 (2010).

[16] D. M. Walters, R. Richert, and M. D. Ediger, Thermal stability of vapor-deposited stable glasses of an organic semiconductor, The Journal of Chemical Physics 142, 134504 (2015).

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

[18] 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 vapour-deposited glasses of toluene, Physical Chemistry Chemical Physics 19, 11089 (2017).

[19] C. Rodríguez-Tinoco, M. Gonzalez-Silveira, J. Ràfols-Ribé, A. Vila-Costa, J. C. Martínez-García, 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).

[20] P. G. Wolynes, Spatiotemporal structures in aging and rejuvenating glasses, Proceedings of the National Academy of Sciences 106, 1353 (2009).

[21] S. Léonard and P. Harrowell, Macroscopic facilitation of glassy relaxation kinetics: Ultrastable glass films with frontlike thermal response, The Journal of Chemical Physics 133, 244502 (2010).

[22] A. Wisitsorasak and P. G. Wolynes, Dynamical heterogeneity of the glassy state, The Journal of Physical Chemistry B 118, 7835 (2014).

[23] R. Gutiérrez and J. P. Garrahan, Front propagation versus bulk relaxation in the annealing dynamics of a kinetically constrained model of ultrastable glasses, Journal of Statistical Mechanics: Theory and Experiment 2016, 074005 (2016).

[24] E. Flenner, L. Berthier, P. Charbonneau, and C. J. Fullerton, Front-mediated melting of isotropic ultrastable glasses, Phys. Rev. Lett. 123, 175501 (2019).

[25] A. Sepúlveda, S. F. Swallen, and M. D. Ediger, Manipulating the properties of stable organic glasses using kinetic facilitation, The Journal of Chemical Physics 138, 12A517 (2013).

[26] C. Rodríguez-Tinoco, J. Ràfols-Ribé, M. González-Silveira, and J. Rodríguez-Viejo, Relaxation dynamics of glasses along a wide stability and temperature range, Scientific reports 6, 1 (2016).

[27] J. Ràfols-Ribé, A. Vila-Costa, C. Rodríguez-Tinoco, A. F. Lopeandía, J. Rodríguez-Viejo, and M. Gonzalez-Silveira, Kinetic arrest of front transformation to gain access to the bulk glass transition in ultrathin films of vapour-deposited glasses, Physical Chemistry Chemical Physics 20, 29989 (2018).

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

[29] R. L. Jack and L. Berthier, The melting of stable glasses is governed by nucleation-and-growth dynamics, The Journal of Chemical Physics 144, 244506 (2016).

[30] G. M. Hocky, L. Berthier, and D. R. Reichman, Equilibrium ultrastable glasses produced by random pinning, The Journal of Chemical Physics 141, 224503 (2014).

[31] C. J. Fullerton and L. Berthier, Density controls the kinetic stability of ultrastable glasses, EPL (Europhysics Letters) 119, 36003 (2017).

[32] M. Lulli, C.-S. Lee, H.-Y. Deng, C.-T. Yip, and C.-H. Lam, Spatial heterogeneities in structural temperature cause kovacs' expansion gap paradox in aging of glasses, Phys. Rev. Lett. 124, 095501 (2020).

[33] L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids, Journal of Statistical Mechanics: Theory and Experiment 2019, 064004 (2019).

[34] M. Avrami, Kinetics of phase change. in general theory, The Journal of Chemical Physics 7, 1103 (1939).

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

[36] C. Scalliet, B. Guiselin, and L. Berthier, Thirty milliseconds in the life of a supercooled liquid, arXiv preprint arXiv:2207.00491 (2022).

[37] K. Dawson, L. A. Kopff, L. Zhu, R. J. McManon, L. Yu, R. Richert, and M. D. Ediger, Molecular packing in highly stable glasses of vapor-deposited tris-naphthylbenzene isomers, The Journal of Chemical Physics 136, 094505 (2012).

[38] V. Testard, L. Berthier, and W. Kob, Influence of the glass transition on the liquid-gas spinodal decomposition, Phys. Rev. Lett. 106, 125702 (2011).

[39] V. Testard, L. Berthier, and W. Kob, Intermittent dynamics and logarithmic domain growth during the spinodal decomposition of a glass-forming liquid, The Journal of Chemical Physics 140, 164502 (2014).

[40] B. Illing, S. Fritschi, H. Kaiser, C. L. Klix, G. Maret, and P. Keim, Mermin-wagner fluctuations in 2d amorphous solids, Proceedings of the National Academy of Sciences 114, 1856 (2017).
METHODS

**Size-polydisperse model glass-former.** We perform molecular dynamics (MD) simulations of a two dimensional size-polydisperse mixture, which is well characterised [35]. The system is composed of soft repulsive spheres whose diameters $\sigma_i$ follow the probability distribution $P(\sigma_i) = A\sigma_i^{-3}$, with A a normalization constant, $\sigma_i \in [\sigma_{min}, \sigma_{max}]$, $\sigma_{min}/\sigma_{max} = 0.45$, and $\sigma_{max} = 1.62$. The pair interaction potential is given by

$$V_{ij}(r) = \varepsilon \left(\frac{\sigma_{ij}}{r}\right)^{12} + c_0 + c_2 \left(\frac{r}{\sigma_{ij}}\right)^2 + c_4 \left(\frac{r}{\sigma_{ij}}\right)^4; \quad (2)$$

where $r = |\mathbf{r}_i - \mathbf{r}_j|$ ($\mathbf{r}_i$ being the position of particle $i$), and non-additive interactions $\sigma_{ij} = 0.5(\sigma_i + \sigma_j)(1 - 0.2|\sigma_i - \sigma_j|)$. We use reduced units based on the particle mass $m$, the energy scale $\varepsilon$, and a microscopic length $\sigma$ defined as the average particle diameter. The time unit is $\tau_{LJ} = \sigma\sqrt{m/\varepsilon}$. The parameters $c_0 = -28\varepsilon/r_c^{12}$, $c_2 = 48\varepsilon/r_c^{14}$, and $c_4 = -21\varepsilon/r_c^{16}$ are introduced in order to ensure that the potential $V_{ij}$ is continuous up to its second derivative at the cutoff distance $r_c = 1.25\sigma_{ij}$.

**Simulating the transformation of an ultrastable glass into a liquid.** The initial configuration is composed of $N = 64000$ particles at a number density $\rho_i = N/L^2 = 1$, in a square box of linear size $L$ and periodic boundary conditions. The initial state is equilibrated using the swap Monte Carlo algorithm [33], at a temperature $T_i = 0.035 \approx T_g/2$, with $T_g$ the experimental glass transition temperature. The corresponding pressure is $P_i = 2.31$. The MD simulations are performed using a Nosé-Hoover thermostat at temperatures $T_a \in [0.1, 0.2]$ and a barostat at pressure $P_i$. The discretisation timestep is set to 0.01. For each $T_a$ we perform MD simulations starting from 8 initial independent configurations, with 3 different sets of initial velocities, therefore gathering statistics for 24 different runs per $T_a$. The errorbars were computed from the statistical error within 95 % of confidence level. In order to test for possible finite size effects, we performed equivalent MD simulations for a system composed of $N = 256000$ particles, running 3 different configurations for $T_a = 0.11, 0.12, 0.13$ and obtained equivalent results.

**Simulations of isolated droplets.** The key technical aspect of these simulations is to allow the system to relax mechanically and possess an inhomogeneous pressure profile. We have found that the following strategy satisfies our needs. We start from the same initial configurations as in the full phase-transformation simulations using an NPT simulation at $T = T_a$ and $P_i$. The expanded glass (G) is reached after a short time, shown in Fig. 2a. Immediately after this rapid expansion, we replicate the system three times in both dimensions, resulting in a glass with $N = 576000$ particles with pressure $P_i$. Such a large system size is needed to ensure that the mechanical response of the glass matrix is not affected by

finite size effects, or, equivalently to ensure that the liquid droplet does not interact mechanically with its replicated images. We insert a liquid droplet in this large glass matrix by performing high-temperature dynamics, at $T = 0.8$ during $10^2 \tau_{LJ}$ in a circular cavity of radius $R_0$, keeping all other particles fixed. This is followed by additional simulations inside the cavity at temperature $T_a$ for $10^4 \tau_{LJ}$ to ensure equilibration inside the cavity. We then perform NVT simulations of the entire system at $T_a$. In the first stages of these simulations, the liquid droplet expands into the glass matrix, which pushes back until mechanical equilibrium is reached, state (D) in Fig. 3c. In this state, we have measured the pressure $P_D$ inside the droplet using different methods. First, we used direct computation from the diagonal components of the pressure tensor in a large portion of the liquid droplet. We confirmed that the potential energy and structural relaxation time correspond to the pressure $P_D$. We also measured the expansion of the droplet by following the displacement field in the radial direction. The resulting density $\rho_D$ agrees with the measured pressure $P_D$. At a given $T_a$, simulations were performed for three different droplet sizes $R_0 = 50, 80, 100$ which gave identical results. We also checked that our results do not depend on the specific initial configuration.

**Simulations of macroscopic fronts.** For large fronts we again start from the expanded glass state (G) but we do not need to replicate the system in both directions. We directly perform high-temperature simulations to transform a slit of width $2W_0$ into a liquid, as shown in Fig. 3a. This is followed by some dynamics at $T_a$. We then perform NPT simulations at pressure $P_i$ and temperature $T_a$. At a given $T_a$, simulations were performed for three different widths $W_0 = 50, 80, 100$. We also checked that our results do not depend on the specific initial configuration.

**Bond-breaking correlation function and phase definition.** We distinguish between the glass and liquid states using the bond-breaking correlation function given by

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

where $n_i(t)$ is the number of neighbours of particle $i$ at time $t$, and $n_i(t|0)$ is the number of neighbors of particle $i$ at time which were also neighbor at $t = 0$. At $t = 0$, the neighbors of particle $i$ are defined as all particles $j \neq i$ whose interparticle distance is smaller than a threshold, $r_{ij}/\sigma_{ij} < 1.35$, which corresponds to the first minima of the rescaled radial distribution function $g(r/\sigma_{ij})$. We define liquid particles as those with $C_B^i(t) \le 0.5$. The bond-breaking correlation function is defined from the ensemble average:

$$C_B(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} C_B^i(t) \right\rangle, \quad (4)$$

where $N$ is the number of particles. One can define from the bond-breaking correlation a characteristic time for
the structural relaxation of the equilibrium bulk system as $C_B(t = \tau_B) = 0.5$.

**Self-intermediate scattering function.** We use an- other common correlation function to characterize the bulk dynamics, namely the self-intermediate scattering function, given by

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

where $\delta\mathbf{r}_i = \mathbf{r}_i(t) - \mathbf{r}_i(0)$ and wave vectors with $|\mathbf{q}| = 6.9$, corresponding to the position of the global maximum in the static structure factor $S(\mathbf{q})$. We define the $\alpha$-relaxation time as $F_s(t = \tau_\alpha) = 1/e$. For large 2d systems $F_s(t)$ can be affected by Mermin-Wagner fluctuations, which depend logarithmically on the system size [40]. To avoid this effect, we performed bulk MD equilibrium sim- ulations of a relatively small system with $N = 2000$.

**Chord length distribution and characteristic length scale.** We employ the chord length distribu- tion [38, 39] to characterize the size of the liquid do- mains. Briefly, we discretize the spatial coordinates on a grids with cells of linear size 1.47, so that the cell is large enough to ensure that there is at least one particle in the cell [36]. Then, we average $C_B^i$ over all the particles in a given cell. We binarize the result so that the cell is either in state 0 (liquid) when $C_B^i \le 0.5$ or 1 (glass) otherwise. The chord length $\ell$ is then defined from the intersection of segments in the $x$ and $y$ directions with the liquid domains. By averaging over all segments in both directions, we get the chord length distribution $P(\ell)$ and we extract a characteristic length scale as the first moment of this distribution.
## SUPPLEMENTARY INFORMATION

The transformation of a glass into a liquid has been previously described using a classical nucleation-and-growth Avrami dynamics [29], which leads to the Avrami equation [34]:

$$x_l(t) = 1 - \exp(-Kt^n), \quad (6)$$

where $x_l$ is the liquid fraction at time $t$, K and $n$ are parameters dependent on the transformation characteristics.

A more general expression for the phase transformation is obtained as follows. The number of nuclei $dN$ that appear during a time interval $d\tau$ is given by $dN = IVd\tau$, with $I$ the nucleation rate per unit volume, and $V$ is the volume of the system. One can suppose that the volume of each nucleus grows isotropically in time, with a radius $r(t)$. Then, the extended volume of the new phase resulting from the nucleation and growth of the nuclei that appeared in $d\tau$ is $dV_{\text{ext}} = g_d r^d dN$, where $d$ is the system's dimension and $g_d$ a geometrical factor (e.g. $g_d = 4\pi/3$ in 3d and $\pi$ in 2d). This extended transformed volume defines the extended fraction $dx_{\text{ext}} = dV_{\text{ext}}/V$. However, part of this extended fraction covers already transformed material. Therefore, the actual liquid fraction that forms during an increment $d\tau$ will be proportional to the fraction of the untransformed phase: $dx_l = dx_{\text{ext}}(1 - x_l)$. Integrating this equation one obtains

$$x_l = 1 - \exp(-x_{\text{ext}}). \quad (7)$$

The Avrami equation is a particular solution of Eq. (7) in the case where the growth velocity $v$ and the nucleation rate $I$ are constant. From these assumptions, one can show that, for constant $I > 0$ then $n = d + 1$ (with $d$ the spatial dimension), while for $I = 0$ and all nucleation sites are present at $t = 0$ one gets $n = d$. Specifically, in the case $d = 2$, a constant $I$ implies that circular domains grow with a linear size $r(t) \sim vt$, so that

$$x_{\text{ext}} = \int_{0}^{t} \pi I r^{2}(\tau) d\tau = \frac{\pi}{3} I v^{2} t^{3},$$

yielding the Avrami exponent $n = d + 1 = 3$.

To understand the anomalous Avrami exponent ($n \approx 4.5$) obtained in our simulations, we introduce a simple modification to the above picture using insights from our MD simulations where we observe a non-constant growth velocity, which goes from $v_1$ observed for isolated droplets, to $v_2 > v_1$ observed for large fronts. We introduce the dimensionless parameter $\alpha = v_2/v_1 > 1$ to quantify this effect.

We model the evolution of the growth velocity so that it smoothly evolves from $v_1$ at short times to $v_2$ at long times. To do so with minimal ingredients, we introduce a crossover length scale $R^*$, or, equivalently a crossover time scale $t^* = R^*/v_1$ controlling the velocity evolution.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Minimal model for time-dependent velocity. We introduce a simple model where the velocity smoothly evolves from $v_1$ to $v_2 > v_1$ at a typical timescale $t^*$, using a sigmoidal function, Eq. (8), which is well approximated by a linear time dependence near the inflection point. Here, $v_1 = 5$, $v_2 = 6.25$, $I = 8 \times 10^{-7}$ and $\Delta t = 2$.</figcaption></figure>

We estimate $R^*$ as the average distance between nucleated droplets [29], namely

$$R^* = \frac{(v_1/I)^{1/(d+1)}}{2},$$

so that $R^*$ is related to the typical droplet size when they start to merge to form large fronts. We also introduce a short duration $\Delta t$ for the crossover to occur leading to a sigmoidal functional form,

$$v(t) = v_1 \frac{1 + \alpha \exp\left(\frac{t-t'}{\Delta t}\right)}{1 + \exp\left(\frac{t-t'}{\Delta t}\right)}; \quad (8)$$

where $t' = t^* + 2\Delta t$, and $t$ is the simulation time. A specific example is shown in Fig. 5.

In the intermediate time regime, when $t \sim t'$, one can perform the Taylor series of Eq. (8) to get

$$v(t \sim t') = v_1 \left[ \frac{\alpha + 1}{2} + \frac{1}{4}(\alpha - 1) \frac{t - t'}{\Delta t} + \mathcal{O}(t^3) \right],$$

allowing us to linearise the velocity in the transition regime as $v(t) \sim v_1(at + b)$ with $a = \frac{\alpha - 1}{4\Delta t}$ and $b = \frac{\alpha + 1}{2} - at'$. In this regime, one gets

$$x_{\text{ext}} = \int_{0}^{t} \pi I r^{2}(\tau) d\tau \\ = \pi I v_{1}^{2} \left[ a^{2} \frac{t^{5}}{5} + 2ab \frac{t^{4}}{4} + b^{2} \frac{t^{3}}{3} \right], \quad (9)$$

suggesting that an effective exponent larger than 3 and close to $n \approx 5$ holds.

To confirm these analytical results, we have performed Monte Carlo simulations of this kinetics, as illustrated in Fig. 6a. In these simulations, we introduce a discrete lattice. The liquid nucleates stochastically with a constant
<figure><img src="image_7.png" /><figcaption>FIG. 6. Monte Carlo simulations. (a) Snapshots of Monte Carlo simulations of the nucleation and growth dynamics with time-dependent velocity. The color codes for glass (blue) and liquid (red). The velocity time evolution is the same as in Fig. 5, yielding $a \sim 0.031$. (b) Liquid fraction versus normalized time for different values of a, which characterizes the steepness of the velocity change, along with different Avrami functions. (c) Effective Avrami exponent n from the data in (b).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Characteristic length scale in Monte Carlo simulations. Snapshots of two different Monte Carlo simulations for $I = 5 \times 10^{-7}$ and $\Delta t = 5$, taken at $t = t_{1/2}$. The characteristic distance between nucleation events, and thus $\ell_{tr}$, is smaller when $v_1 < v_2$ (b) compared to the standard Avrami case where $v_1 = v_2$ (a).</figcaption></figure>

nucleation rate $I$ on the sites of the lattice. The nucleated sites then grow circular domains with a velocity given by Eq. (8), independently of other liquid domains. The model has four parameters: the nucleation rate $I$, the two velocities $v_1$ and $v_2 = \alpha v_1$ and the timescale $\Delta t$ controlling the sharpness of the change between these two values. We find, however, that it is the parameter $a = (\alpha - 1)/\Delta t$, which characterizes the steepness of the velocity change, which mostly controls the obtained kinetics.

We can extract the time evolution of the fraction of liquid $x_l(t)$ and the effective Avrami exponent $n$ from these Monte Carlo simulations for different values of $a$, as shown in Fig. 6b. For $a=0$, the velocity is constant throughout the simulation, we recover the classical Avrami kinetics, yielding $n=3$, see Fig. 6c. However, as soon as the velocity changes, one finds a larger Avrami exponent, going up to $n=5$ for $a=0.042$.

Finally, it is interesting to use this minimal model to explore the two characteristic lengths discussed in the main text, $\ell_{tr}$ and $\ell_c$. When $v_1 = v_2$, we expect that $\ell_{tr} = \ell_c$. However, by decreasing the value of $v_1$ to explore the regime $v_1 < v_2$ (keeping all other parameters unchanged), the transformation time becomes longer and thus $\ell_c = v_2 t_{tr}$ increases. Nevertheless, as discussed in the main text, larger transformation times at constant nucleation rate implies the presence of a larger number of nuclei at time $t_{1/2}$, decreasing the distance between them and therefore decreasing $\ell_{tr}$. The generic expectation that $\ell_{tr} < \ell_c$ is consistent with our MC simulations.

In Fig. 7 we show two different simulations with the same nucleation frequency but distinct velocity evolution. One can directly observe how decreasing $v_1$ results in a significantly smaller characteristic length $\ell_{tr}$.

This minimal analytic model where a single modification to the Avrami picture is introduced is therefore able to reproduce the results obtained from the MD simulations described in the main text.