# Finite-dimensional vestige of spinodal criticality above the dynamical glass transition

Ludovic Berthier, Patrick Charbonneau, Joyjit Kundu

► To cite this version:

Ludovic Berthier, Patrick Charbonneau, Joyjit Kundu. Finite-dimensional vestige of spinodal criticality above the dynamical glass transition. Physical Review Letters, 2020, pp.108001. 10.1103/Phys-RevLett.125.108001. hal-02986305

HAL Id: hal-02986305

https://hal.science/hal-02986305

Submitted on 2 Jun 2021

**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
Finite-dimensional vestige of spinodal criticality above the dynamical glass transition

Ludovic Berthier,¹,² Patrick Charbonneau,³,⁴ and Joyjit Kundu³

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

²Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge CB2 1EW, United Kingdom

³Department of Chemistry, Duke University, Durham, North Carolina 27708, USA

⁴Department of Physics, Duke University, Durham, North Carolina 27708, USA

(Dated: December 30, 2019)

Finite-dimensional signatures of spinodal criticality are notoriously difficult to come by. The dynamical transition of glass-forming liquids, first described by mode-coupling theory, is a spinodal instability preempted by thermally activated processes that also limit how close the instability can be approached. We combine numerical tools to directly observe vestiges of the spinodal criticality in finite-dimensional glass formers. We use the swap Monte Carlo algorithm to efficiently thermalise configurations beyond the mode-coupling crossover, and analyze their dynamics using a scheme to screen out activated processes, in spatial dimensions ranging from $d = 3$ to $d = 9$. We observe a strong softening of the mean-field square-root singularity in $d = 3$ that is progressively restored as $d$ increases above $d = 8$, in surprisingly good agreement with perturbation theory.

# I. INTRODUCTION

Spinodals predicted by mean-field theories do not exist in finite-dimensional systems because thermal (or other) fluctuations destabilize the precursor metastable state before the critical point can ever be reached [1–4]. An Ising system prepared with the metastable magnetization, for instance, grows nuclei (instantons [5, 6]) of the opposite magnetization well ahead of the formal instability region. As the single-flip relaxation dynamics slows down critically, that of nucleation accelerates [7]; cluster moves only worsen the imbalance. Numerical studies have nevertheless observed convincing hints of the Ising spinodal pseudocriticality in high enough spatial dimension $d$. Nucleation kinetics being exponentially suppressed as $d$ increases, the pseudocritical power-law scaling can then be made sufficiently extended [7].

More theoretically enticing than the Ising spinodal is that of models with disorder, which capture the essence of systems ranging from magnetic [8] to mesoporous [9, 10] materials and also appear in social sciences and economics [11]. These models exhibit a rich set of activated processes, such as avalanches [8, 12] and hopping [13, 14], in addition to nucleation, which makes their criticality especially challenging to scrutinize. Hence, although the Ginzburg criterion for the corresponding cubic field theory without activation gives an upper critical dimension $d_u = 8$ [15–17], it is unclear how relevant the associated pseudocriticality might be in any given system. Below $d_u$, perturbative expansions relying on dimensional reduction [16, 18] or direct expansion [19] further do not occur. Even if they did, one may expect non-perturbative fluctuations to also contribute [20]. Although these fluctuations should limit the relevance of perturbative approaches and question the very existence of an upper critical dimension, very little is known about their effects [21, 22], emphasizing further the need for quantitative results in finite-dimensional systems.

This theoretical haze has not held back the use of

mode-coupling theory (MCT) to describe the dynamics of supercooled liquids [23], which has been the subject of countless numerical and experimental tests [23, 24]. Although unclear in its initial derivations, it is now understood that MCT [25] as well as the mean-field $d \to \infty$ description of liquids [26] indeed correspond to the limit of stability of the glass phase upon heating (or decompressing). The dynamical glass transition should thus be described as a spinodal (thermodynamic) instability in the presence of quenched disorder [20, 27]. Unfortunately, the spinodal is found to be totally hidden by finite-dimensional effects in direct free energy calculations [28, 29].

The associated critical scaling laws of the structural relaxation time $\tau_\alpha$ and of time correlation functions upon approaching the avoided MCT (dynamical) transition from the equilibrium liquid are instead much more frequently examined [23, 24, 26]. These quantities are straightforwardly measured in both simulations and experiments, but they are also non-universal, i.e., model dependent. Theoretical predictions for the associated critical exponents are not only sensitive to the spatial dimension (even above $d_u$) and to activated processes, but also to fine details of the liquid structure and pair interactions. These predictions are thus typically of limited quantitative validity in finite $d$, but this limitation did not prevent MCT from making valid predictive statements regarding the glassy dynamics of a variety of materials [23, 24].

In addition to non-universal scaling laws, the spinodal criticality is associated to a few universal signatures. In particular, a square root singularity of the Edwards-Anderson parameter directly follows [23, 30], which is dynamically accessible as the plateau height in time correlation functions or the typical cage size in particle displacements. Treating fluctuations beyond mean-field leads also to universal predictions regarding the behaviour of four-point susceptibilities [25, 31, 32]. Yet, because of the computational difficulty of equilibrating
liquids beyond the avoided dynamical transition, and of the lack of experimental methods to screen out activated processes, it remains difficult to assess signatures of the spinodal, even in simple glass-forming liquids. As a result, the validity of the square root singularity remains a debated issue [33].

In this work, we exploit a recent implementation of swap Monte Carlo (SWAP) for continuously polydisperse systems, which bypasses the sluggishness associated with approaching the avoided dynamical transition [34–36], in order to probe the spinodal criticality beyond the MCT crossover in equilibrium conditions. By carefully controlling for activated processes, our analysis manages to extract a sufficiently long scaling regime of the typical cage size to estimate effective critical exponents and the putative presence of the square-root singularity across several space dimensions from the experimentally relevant $d = 3$, where we conclude that the singularity is considerably softened, up to $d = 9$ where a nearly perfect square-root scaling can be convincingly observed.

## II. SIMULATION DETAILS

We consider a continuously polydisperse system of $N$ hard spheres under periodic boundary conditions in a simulation box of constant volume $V$. A hypercubic box is used in $d = 3, 4, 5, 6$, and $8$, while in $d = 7$ and in $d = 9$ we use the Wigner-Seitz cell of the checkerboard lattice that decreases the number of simulated particles while preserving the same effective box size. The particle-size distribution is $P(\sigma) = K/\sigma^3$, with normalization constant $K$ for diameters $\sigma \in \{\sigma_{max}, \sigma_{min}\}$ where $\sigma_{max}, \sigma_{min}$ are the maximum and minimum diameters for a given polydispersity. The average diameter $\bar{\sigma}$ sets the unit of length, the degree of polydispersity is defined by the standard deviation of the diameter distribution, and the packing fraction is $\varphi = \rho\bar{V}_d$ for a number density $\rho = N/V$ and average volume of a $d$-dimensional hypersphere $\bar{V}_d$. The degree of polydispersity is chosen to be the minimum needed for the SWAP efficiency to saturate, i.e., 23% for $d = 3$, 10% for $d = 4–8$, and 8% for $d = 9$ [36], and suitably optimized SWAP sampling is used to equilibrate initial configurations. This approach ensures fast structural relaxation without crystallization or fractionation. Structural equilibration is notably validated by the complete decay of the self-part of the particle-scale overlap function

$$Q(t) = \frac{1}{N} \sum_{i=1}^{N} \Theta(a - |\mathbf{r}_i(t) - \mathbf{r}_i(0)|), \quad (1)$$

where $\Theta$ is the Heaviside function and $a = 0.3\bar{\sigma}$ is about the typical particle cage size [29, 35–37]. The associated structural relaxation time, $\tau_\alpha$, is defined as $Q(\tau_\alpha) = 1/e$. From these initial equilibrium configurations, multiple simulations are then run with a purely local Monte Carlo dynamics. This computational scheme achieves equilibration at densities 3–8% above the estimated dynamical

<figure><img src="image_2.png" /><figcaption>FIG. 1. Time evolution of the MSD $\langle\Delta r^2(t)\rangle$ and of the non-Gaussian parameter $\alpha_2(t)$ for different densities in $d = 3$. Vertical lines denote the start and end times of the window over which the cage size is measured for different densities.</figcaption></figure>

transition, depending on $d$ [35, 36]. This strategy opens a comfortable regime to study glassy dynamics approaching the mode-coupling crossover from the arrested phase, unavailable to previous computational work.

## III. TYPICAL CAGE SIZE

We first consider the size of the typical cage, $\hat{\Delta}$, which in the MCT and the mean-field description of hard spheres is expected to scale as

$$\hat{\Delta}(\varphi) = \hat{\Delta}_d - A_d(\varphi - \varphi_d)^{1/\delta}, \quad (2)$$

with $1/\delta = 1/2$ and $\hat{\Delta}_d = \hat{\Delta}(\varphi_d)$ for densities beyond the dynamical transition, $\varphi > \varphi_d$. In the high-dimensional limit, $\hat{\Delta}$ could be extracted from the long-time plateau of the mean square displacement (MSD) of an individual particle, $\hat{\Delta} = \Delta r_i^2(t) = \int r^2 G_s(\mathbf{r}, t) d\mathbf{r}$, where $G_s(\mathbf{r}, t)$ is the self-part of the van Hove function. In finite dimensions, two difficulties arise. First, caging is heterogeneous, and hence the full distribution of cages must be considered. Second, a sharp plateau in the MSD can only be identified much beyond the (avoided) dynamical transition, i.e. too far beyond the regime of interest. We consider the second effect first. Because activated processes interfere with the formation of the MSD plateau, the size of the transient cage must be extracted over a finite time window. In order to identify the mean-field-like dynamical caging regime, we rely on the non-Gaussian parameter, $\alpha_2(t) = \frac{d}{d+2} \frac{\langle r^4(t) \rangle}{\langle r^2(t) \rangle^2} - 1$, as illustrated in Fig. 1. More specifically, the upper bound of the time window is set at 20% of peak non-Gaussianity,
<table><caption>TABLE I. Fit parameters $A_d$ and $\hat{\Delta}_d$ for Eq. (3) for different dimensions. The results for $\varphi_d$ are obtained by standard MCT dynamical scaling [26]. The direct evaluation of the typical cage size at the estimated $\varphi_d$ validates the value of the fitted quantity. Error bars on the fit parameters are determined from quality of the fit ($R^2$) to Eq. (2).</caption><thead><tr><th>d</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><td>φ<sub>d</sub></td><td>0.600(2)</td><td>0.410(2)</td><td>0.277(1)</td><td>0.1808(8)</td><td>0.1147(6)</td><td>0.0716(3)</td><td>0.0426(2)</td></tr><tr><td>Δ̂<sub>d</sub></td><td>0.015(2)</td><td>0.021(2)</td><td>0.020(2)</td><td>0.018(1)</td><td>0.016(1)</td><td>0.010(1)</td><td>0.013(1)</td></tr><tr><td>Δ(φ<sub>d</sub>)</td><td>0.0148(5)</td><td>0.0208(5)</td><td>0.0196(5)</td><td>0.0175(6)</td><td>0.0161(6)</td><td>0.0101(4)</td><td>0.0129(5)</td></tr><tr><td>1/δ</td><td>0.74(9)</td><td>0.68(9)</td><td>0.65(8)</td><td>0.61(9)</td><td>0.59(9)</td><td>0.53(9)</td><td>0.55(10)</td></tr><tr><td>A<sub>d</sub></td><td>0.11(2)</td><td>0.16(2)</td><td>0.15(2)</td><td>0.12(2)</td><td>0.11(2)</td><td>0.04(1)</td><td>0.09(3)</td></tr></tbody></table>

i.e., the maximum of $\alpha_2(t)$, and the lower bound is set to the short-time plateau of $\alpha_2(t)$. Because the lower bound is much smaller than the upper bound, results are insensitive to its precise value, while the upper bound only weakly affects the subsequent analysis as long as it is chosen consistently. Note that for packing fractions much above $\varphi_d$ (in $d = 3$, $\varphi \gtrsim 0.635$), the system does not relax within the simulation timescale (with standard dynamics), which results in a clear plateau in the MSD (Fig. 1). The upper bound is then even tighter. Note also that generalizing the MSD to further suppress the contribution of activated processes as in Ref. [13] markedly flattens the MSD, but does not quantitatively affect the subsequent analysis (see Appendix A).

Once the timescale for identifying the cage dynamics is set, we can measure the cage size for each particle, $\Delta_i$,

<figure><img src="image_3.png" /><figcaption>FIG. 2. Cage size distributions, $P(\Delta_i)$, in (a) $d = 3$ for $\varphi = 0.6005, 0.6032, 0.6111, 0.6208, 0.6309$, and $0.6414$ and (b) $d = 6$ for $\varphi = 0.1810, 0.1829, 0.1842, 0.1865, 0.1892$, and $0.1916$. Note that $\varphi_d = 0.600(2)$ and $0.1808(8)$ in $d = 3$ and $6$ respectively. Fat tails at large displacements persist over the whole density regime accessible in simulations. The estimator $\hat{\Delta} = \operatorname{argmax}_{\Delta_i} P(\Delta_i)$ (vertical dashed lines) nonetheless monotonically shifts to smaller values as $\varphi$ increases.</figcaption></figure>

and average over samples and particles to define the cage size distribution $P(\Delta_i)$. The results in Fig. 2 show that cages generally tighten as density increases in all dimensions, but that fat tails at large displacements persist for all densities. These tails deviate significantly from the log-normal forms reported in some mean-field models [13]. Although relatively little theoretical guidance is available as to what the proper functional form for $P(\Delta_i)$ should be [38], our observations suggest that activated processes are not fully eliminated from the MSD analysis. To further sidestep this issue, we use as estima-

<figure><img src="image_4.png" /><figcaption>FIG. 3. Critical scaling of the typical cage size in $d = 3...9$ as well as (inset) fitted power-law exponents, $1/\delta$. For visual clarity, data is vertically shifted by a factor of $3^{d-3}$. Error bars reflect the measurement uncertainty of $\Delta(\varphi)$ only with $\Delta_d$ and $\varphi_d$ here chosen to optimize the quality ($R^2$) of the fit. Inset: the critical exponent $1/\delta$ versus $d$. A clear change in critical scaling is seen around the upper critical dimension, $d_u = 8$, predicted from perturbative approaches. For $d > d_u$, the results are consistent with the mean-field prediction ($1/\delta = 1/2$, dashed line); for $d < d_u$ the deviation from the mean-field result is approximately linear $1/\delta - 1/2 = B(d_u - d)$ with $B = 0.049(3)$.</figcaption></figure>
tor of the typical cage size the mode of the distribution, $\hat{\Delta} = \operatorname{argmax}_{\Delta_i} P(\Delta_i)$, which is much less sensitive to the activated processes that appear in the fat tail of the distribution than the mean cage size, but converges to the same quantity as $d \to \infty$.

In order to assess the critical scaling of the typical cage size, $\hat{\Delta}(\varphi)$ is fitted to Eq. (2) using $\delta$, $\hat{\Delta}_d$ and $A_d$ as parameters, while $\varphi_d$ is obtained independently from the growth of the relaxation time $\tau_\alpha(\varphi)$ [36] (Fig. 3, TableI). Because of the uncertainty on $\varphi_d$, the cage size at $\varphi_d$, is not directly measured; but the fitted value is consistent with the direct estimate, which validates our approach. The values of the fit parameters and direct measurements are listed in Table I. Note that both the fit error on $\hat{\Delta}_d$ at a given $\varphi_d$ and the propagated uncertainty from $\varphi_d$ are then included. In contrast to earlier (cruder) estimates [39], we find that for a given polydispersity, $\hat{\Delta}_d$ decreases monotonically with increasing dimension [30]. More significantly, we also find that $1/\delta$ decreases monotonically, and nearly linearly, as dimension increases, from 0.74(9) in $d=3$ down to values numerically indistinguishable from the mean-field prediction, $1/\delta = 1/2$, around $d \ge 8$. Although the precise numerical estimates of $\delta$ are fairly robust to the details of the above analysis, they may still be fragile to the overall scheme. The detection of the upper critical dimension in the vicinity of $d=8$, and the systematic softening of the square-root singularity below $d_u = 8$ are nonetheless numerically robust (see Appendix A and B). Interestingly, the latter is in sharp contrast with the prediction from dimensional reduction [16, 18].

## IV. CONCLUSION

The finite-dimensional vestige of the spinodal criticality associated with the dynamical transition of glass forming liquids has here been characterized by numerical simulations using the SWAP algorithm and a careful screening of activated processes across a broad range of spatial dimensions. Our simulations reveal that the square-root singularity, often used to describe experimental measurements in molecular glass formers, is dramatically softened in $d=3$ for hard spheres, a canonical model for testing MCT predictions. The measured effective exponent, $1/\delta \approx 0.75$, can still be considered as indirect evidence for an underlying avoided singularity, because $1/\delta = 1$ would be trivially expected for a featureless evolution of the cage size. Further, the slow variation of $1/\delta$ towards $1/2$ with spatial dimension $d$ suggests that strong deviations from mean-field criticality exist even in large spatial dimensions. It takes simulations in dimensions $d \ge 8$ to observe direct signatures of the square-root scaling that underlies the mean-field dynamical glass transition. The agreement between our results with perturbative approaches is nevertheless surprising, given the expected role of nonperturbative physics in the avoidance of the dynamical glass transition [20], which can in principle

<figure><img src="image_5.png" /><figcaption>FIG. 4. Generalized mean square displacement for different $n$ and different $\varphi$ in $d=3$.</figcaption></figure>

persist even above $d=8$. A possible explanation might be the relative insensitivity of our specific estimators to these effects.

Because this spinodal critical point is part of a broad universality class [20], we expect our results to apply to a variety of other systems, for which the interplay between activation and criticality might be harder to control. Most crucially, these results further motivate the use of the dynamical criticality and of the mean-field description in describing the behavior of finite-dimensional liquid glass formers. It was recently shown that deviations from the dynamic glass transition can be studied considering the degree of localization of unstable modes in the potential energy landscape [40]. Localized excitations are indeed expected to disappear as $d$ increases, and our study thus also motivates verifying this prediction in larger dimensions.

## ACKNOWLEDGMENTS

We thank G. Biroli, G. Carra, S. Franz, G. Tarjus, and F. Zamponi for stimulating discussions, as well as Tom Milledge for technical support. We acknowledge support from the Simons Foundation grant (#454933, LB, # 454937, PC). Simulations were performed at Duke Compute Cluster (DCC), and at the Extreme Science and Engineering Discovery Environment (XSEDE), supported by National Science Foundation grant number ACI-1548562.

## Appendix A: Generalized mean square displacement

In order to further suppress the contribution of activated processes, we can generalize the definition of the mean square displacement as in Ref. [13]

$$ \Delta r^2(t) = \lim_{n \to 0} \left\{ \int r^{2n} G_s(r,t) S_{d-1}(r) dr \right\}^{1/n}, \quad (A1) $$
<figure><img src="image_6.png" /><figcaption>FIG. 5. Critical scaling of the typical cage size in $d = 4$ calculated from the generalized MSD with packing fraction for different values of $n$. Inset: the critical exponent $\delta$ is robust against the choice of $n$. All data are for $N = 1000$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Cage size distribution for $d = 4$ at $\varphi = 0.411, 0.416$, and $0.431$ for different system sizes. The typical cage size is found to be robust against changes to the system size.</figcaption></figure>

where $S_{d-1}(r) = 2\pi^{d/2}r^{d-1}/\Gamma(d/2)$ is the surface area of a $d$-dimensional hypersphere, and $G_s(r, t)$ is the self-part of the van Hove function. Note, however, that a direct use of this expression is numerically impractical, because it requires perfect sampling around the maximum of the van Hove function. In practice, we thus consider the convergence of the analysis upon reducing $n$ from the standard MSD definition with $n = 1$. Figure 4 shows that reducing $n$ markedly enhances plateau formation, as expected. Thus, for smaller $n$, the cage size distribution can be measured over a longer time window, which is particularly helpful as $\varphi$ approaches $\varphi_d$. For instance, with $n = 0.07$ the upper bound of the time window can be set at 25% (instead of 20%) of the peak position of $\alpha_2(t)$. The critical exponent extracted from fitting results for different values of $n$ to Eq. 2 is robust. As an illustration, results for $d = 4$ are reported in Fig. 5.

## Appendix B: Finite-size Effects

To ensure the robustness of the critical scaling with system size, we consider the cage size distributions at three different densities above $\varphi_d$ for two system sizes in Fig. 6. The typical cage size or the mode of the distributions at a given density are here indistinguishable. The critical scaling of the typical cage is thus insensitive to the system size, in the regime considered here.
[1] J. D. Gunton and M. C. Yalabik, Phys. Rev. B **18**, 6199 (1978).

[2] J. D. Gunton and M. Droz, *Intro to the Theory of Metastable and Unstable States* (Springer, 1983).

[3] T. S. Ray and W. Klein, J. Stat. Phys. **61**, 891 (1990).

[4] F. H. Stillinger, *Energy Landscapes, Inherent Structures, and Condensed-Matter Phenomena* (Princeton University Press, 2015).

[5] J. S. Langer, Ann. Phys. **41**, 108 (1967).

[6] J. S. Langer, Ann. Phys. **54**, 258 (1969).

[7] T. S. Ray, J. Stat. Phys **62**, 463 (1991).

[8] J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Phys. Rev. Lett. **70**, 3347 (1993).

[9] G. J. Aubry, F. Bonnet, M. Melich, L. Guyon, P. Spathis, F. Despetis, and P.-E. Wolf, Phys. Rev. Lett. **113**, 085301 (2014).

[10] F. Detcheverry, E. Kierlik, M. L. Rosinberg, and G. Tarjus, Langmuir **20**, 8006 (2004).

[11] J.-P. Bouchaud, J. Stat. Phys. **151**, 567 (2013).

[12] O. Perković, K. A. Dahmen, and J. P. Sethna, Phys. Rev. B **59**, 6106 (1999).

[13] P. Charbonneau, A. Ikeda, G. Parisi, and F. Zamponi, Proc. Natl. Acad. Sci. U.S.A. **109**, 13939 (2012).

[14] The impact of thermal fluctuations can be eliminated by looking at the zero temperature spinodal of the random-field Ising model (RFIM), which exhibits system-spanning avalanches [8]. However, in this case, rare athermal nonperturbative fluctuations drastically change the nature of this spinodal as well as its location, making its study much less controlled than otherwise imagined [20].

[15] G. Biroli and J.-P. Bouchaud, J. Phys.: Condens. Matter **19**, 205101 (2007).

[16] S. Franz, G. Parisi, F. Ricci-Tersenghi, and T. Rizzo, Eur. Phys. J. E **34**, 102 (2011).

[17] S. Franz, H. Jacquin, G. Parisi, P. Urbani, and F. Zamponi, Proc. Natl. Acad. Sci. U.S.A. **109**, 18725 (2012).

[18] F. Zhong, Front. Phys. **12**, 126402 (2017).

[19] T. Rizzo, Phys. Rev. B **94**, 014202 (2016).

[20] S. K. Nandi, G. Biroli, and G. Tarjus, Phys. Rev. Lett. **116**, 145701 (2016).

[21] V. Lubchenko and P. G. Wolynes, J. Chem. Phys. **119**, 9088 (2003).

[22] S. M. Bhattacharyya, B. Bagchi, and P. G. Wolynes, Proc. Nat. Acad. Sci. U.S.A. **105**, 16077 (2008).

[23] W. Götze, *Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory* (Oxford University Press, 2012).

[24] W. Götze, J. Phys.: Condens. Matter **11**, A1 (1999).

[25] A. Andreanov, G. Biroli, and J.-P. Bouchaud, Europhys. Lett. **88**, 16001 (2009).

[26] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, Annu. Rev. Condens. Matter Phys. **8**, 265 (2017).

[27] S. Franz and G. Parisi, Phys. Rev. Lett. **79**, 2486 (1997).

[28] L. Berthier, Phys. Rev. E **88**, 022313 (2013).

[29] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, Proc. Natl. Acad. Sci. U.S.A. **114**, 11356 (2017).

[30] G. Parisi and F. Zamponi, Rev. Mod. Phys. **82**, 789 (2010).

[31] C. Donati, S. Franz, S. C. Glotzer, and G. Parisi, J. Non-Cryst. Solids **307**, 215 (2002).

[32] L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob, K. Miyazaki, and D. R. Reichman, J. Chem. Phys. **126**, 184503 (2007).

[33] L. Berthier, Phys. Rev. E **76**, 011507 (2007).

[34] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Phys. Rev. Lett. **116**, 238002 (2016).

[35] A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev. X **7**, 021039 (2017).

[36] L. Berthier, P. Charbonneau, and J. Kundu, Phys. Rev. E **99**, 031301(R) (2019).

[37] L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa, and S. Yaida, Nat. Comm. **10**, 1508 (2019).

[38] P. Chaudhuri, L. Berthier, and W. Kob, Phys. Rev. Lett. **99**, 060604 (2007).

[39] P. Charbonneau, Y. Jin, G. Parisi, and F. Zamponi, Proc. Natl. Acad. Sci. U.S.A. **111**, 15025 (2014).

[40] D. Coslovich, A. Ninarello, and L. Berthier, Scipost **7**, 077 (2019).