# Relaxation dynamics in the energy landscape of glass-forming liquids

Yoshihiko Nishikawa, Misaki Ozawa, Atsushi Ikeda, Pinaki Chaudhuri, Ludovic Berthier

► To cite this version:

Yoshihiko Nishikawa, Misaki Ozawa, Atsushi Ikeda, Pinaki Chaudhuri, Ludovic Berthier. Relaxation dynamics in the energy landscape of glass-forming liquids. Physical Review X, 2022, pp.021001. 10.1103/PhysRevX.12.021001. [hal-03636595](https://hal.science/hal-03636595)

HAL Id: hal-03636595

https://hal.science/hal-03636595

Submitted on 31 May 2022

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
# Relaxation Dynamics in the Energy Landscape of Glass-Forming Liquids

Yoshihiko Nishikawa,¹ Misaki Ozawa,² Atsushi Ikeda,³ Pinaki Chaudhuri,⁴ and Ludovic Berthier¹,⁵

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

²Laboratoire de Physique de l'Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, Paris 75005, France

³Graduate School of Arts and Sciences, The University of Tokyo, Tokyo 153-8902, Japan

⁴The Institute of Mathematical Sciences, C.I.T. Campus, Taramani, Chennai 600 113, India

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

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

(Received 5 July 2021; revised 25 October 2021; accepted 2 February 2022; published 1 April 2022)

We numerically study the zero-temperature relaxation dynamics of several glass-forming models to their inherent structures, following quenches from equilibrium configurations sampled across a wide range of initial temperatures. In a mean-field Mari-Kurchan model, we find that relaxation changes from a power law to an exponential decay below a well-defined temperature, consistent with recent findings in mean-field *p*-spin models. By contrast, for finite-dimensional systems, the relaxation is always algebraic, with a nontrivial universal exponent at high temperatures crossing over to a harmonic value at low temperatures. We demonstrate that this apparent evolution is controlled by a temperature-dependent population of localized glassy excitations. Our work unifies several recent lines of studies aiming at a detailed characterization of the complex potential energy landscape of glass formers, and challenges both mean-field and real space descriptions of glasses.

DOI: 10.1103/PhysRevX.12.021001

Subject Areas: Chemical Physics
Computational Physics
Statistical Physics

## I. INTRODUCTION

Many systems of scientific interest are described as “complex,” even though definitions of complexity may vary across scientific fields [1]. For many-body interacting systems, the potential energy landscape $E(r)$, which describes the potential energy $E$ of the system as a function of the complete set of coordinates $\{r\}$ of its constituents, has become a central object of study [2,3]. It serves both empirical goals, for instance, to picture the dynamic evolution of a system in a “rugged” landscape [4], but can also be described mathematically very precisely [2,5,6]. The detailed characterization and dynamic exploration of complex potential energy landscapes are important problems for amorphous materials [4,7,8], optimization problems [9], machine learning algorithms [10,11], and other disordered systems [12].

Since the work of Goldstein [13], the physics of glassy systems is often described in terms of the properties of their potential energy landscapes. The large number of energy minima connected by complex dynamic pathways is

typically invoked in introductory lectures about amorphous media [4], and the sketch of complex energy landscapes very often accompanies the interpretation of experimental measurements [14], which makes this object more than a pure theoretical curiosity. Analytically, the properties of the potential energy landscape of glass-forming models have been studied extensively at the mean-field level through the analysis of fully connected disordered spin models, such as *p*-spin models. In this limit, the phase space can be divided into long-lived metastable states (or, pure states), and both free-energy and energy landscapes can be studied in great detail, thus providing a firm relation between the landscape structure and the thermodynamics and dynamics of the system [5,6,15]. Current efforts in this area concern the analysis of dynamic pathways [16], or corrections to mean field [17].

In finite dimensions, the study of energy minima, or inherent structures, first gained momentum when Stillinger and Weber transformed Goldstein’s ideas into concrete tools to both explore and exploit the potential energy landscape of glasses [4,18]. A key step is the tiling of the equilibrium configuration space, pertinent to describe physical properties, into basins of attraction surrounding energy minima. It is this mapping which putatively connects the thermodynamic and dynamic properties of glass formers to the topography of their potential energy

<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>
landscape, although the relevance of such an approach has often been debated [19,20], because the pure states defined in the mean-field limit do not exist in finite dimensions [21]. The analysis of energy minima has been used to estimate the configurational entropy [7], while saddle points were discussed in connection with the dynamic mode-coupling crossover [22–25]. However, these approaches do not have the same level of rigor as those in *p*-spin models since inherent structures are different from pure states [21,26,27]: inherent structures are configurations that are energetically stable against infinitesimal particle moves, whereas the pure states are defined as free-energy minima. The structure of the potential energy landscape and its precise relationship with dynamics and thermodynamics remain under intense scrutiny [8,28]. In particular, the role of excitations in the potential energy landscape has been discussed in connection with sound propagation [29], specific heat [30], and vibrational [31] and mechanical properties [32].

Virtually all studies of glassy landscapes start by “instantaneously” relaxing configurations to the “nearest” energy minimum, following known numerical recipes [33]. Strangely, however, only a few studies have been dedicated to the physical processes at play during the energy minimization itself [34–40]. In our view, this represents an important vacuum because this relaxation dynamics in fact provides a convenient way to navigate the potential energy landscape and explore its geometry and the nature as well as interactions between excitations that are relevant to describe glassy materials. Suppose, for instance, that the landscape is simple and smooth. Using steepest descent dynamics, the system should then settle in an inherent structure very quickly while, on a rugged landscape, the system meanders and crosses many saddles during relaxation [41]. Similarly, the steepest descent dynamics obtained within kinetically constrained lattice models simply stems from a noninteracting set of excited defects and is therefore unremarkable [19,42]. Thus, in the context of glassy systems, the steepest descent dynamics probes the detailed structure of the potential energy landscape, potentially illuminates its connection to the physical dynamics, and provides novel constraints on physical descriptions of glassy excitations.

Recently, the analysis of steepest descent in mixed mean-field *p*-spin glass models revealed the existence of two important characteristic temperatures [35]. First, starting from initial states prepared at high temperatures *T*, the energy density of the final inherent state is constant for *T* > *T*<sub>onset</sub>, and it decreases with decreasing *T* when *T* ≤ *T*<sub>onset</sub>. This sharp onset temperature does not affect the relaxation dynamics itself which obeys a nontrivial power-law time dependence as long as *T* ≥ *T*<sub>SF</sub>. By contrast, the decay is exponentially fast below *T*<sub>SF</sub> (where SF is state following). This implies that the system is

always close to an energy minimum for *T* ≤ *T*<sub>SF</sub> in which it converges very quickly by steepest descent. The critical temperature *T*<sub>SF</sub> < *T*<sub>onset</sub> also reflects a change in the structure of the potential energy landscape, as inherent states have a marginal density of states above *T*<sub>SF</sub>, which becomes gapped below. Within *p*-spin models, these two characteristic temperatures are unrelated to the equilibrium dynamics, which becomes nonergodic at the mode-coupling temperature *T*<sub>MCT</sub>, distinct from both *T*<sub>SF</sub> and *T*<sub>onset</sub>, showing that even at mean-field level free-energy and energy landscapes are different objects.

In numerical studies, the relaxation dynamics in *D* = 2 and 3 (*D* is the space dimension) harmonic spheres just above jamming was recently studied starting from high temperatures including random configurations at *T* = ∞ [34], and a power-law time decay was found with a nontrivial, dimension-dependent exponent. Another recent work [36] explores the statistics of single particle displacements between initial and final configurations in several three-dimensional models and reports the existence of a crossover temperature separating high-temperature from low-temperature behaviors. These interesting studies do not provide a complete physical picture of the relaxation dynamics toward energy minima, neither do they assess the existence of the critical temperature *T*<sub>SF</sub> found in mean-field approaches. The universality of the power-law time dependence found near jamming across models, and even the effect of spatial dimension and initial temperatures, were not fully elucidated either.

Here, we provide a comprehensive numerical study of the steepest descent dynamics in generic glass-forming liquids. We address its dimensionality, universality, and initial stability dependences by studying a mean-field Mari-Kurchan (MK) model and three finite-dimensional models in two, three, four, and eight dimensions using a wide range of initial states obtained through the swap Monte Carlo algorithm. We numerically detect the predicted mean-field transition at *T*<sub>SF</sub> in the Mari-Kurchan model. However, the transition is absent in all finite-dimensional models, where it is replaced by a smooth temperature evolution between two nontrivial limits that we analyze in detail. We show that this crossover is controlled by a finite population of localized defects where particle rearrangements take place during the minimization, with the overall concentration of these defects decreasing at lower temperatures. Therefore, finite-dimensional glass-forming systems at finite temperatures cannot be seen as inherent structures excited by small thermal fluctuations, since they are described neither by mean-field energy landscapes nor by a simple picture of noninteracting localized defects. Our results provide a complete picture of the relaxation dynamics in glassy landscapes and illuminate the role, nature, and interactions of localized defects in finite-dimensional structural glasses.
## II. RESULTS

### A. Steepest descent dynamics

We numerically solve the equations of motion of steepest descent dynamics,

$$ \zeta \frac{dr_i}{dt} = - \frac{\partial E}{\partial r_i}, \qquad (1) $$

starting at time $t = 0$ from an equilibrium configuration prepared at initial temperature $T$, where $\zeta$ is the damping coefficient and $E$ is the potential energy. The time unit is $\tau_0 = \zeta \ell^2 / v_0$, where $\ell$ is the unit length scale and $v_0$ is the unit energy for particle interactions. In Eq. (1), energy is dissipated via a uniform background. We have not tested more complicated dissipation mechanisms such as used in dense particle suspensions [43]. Note that the dynamics in Eq. (1) is fully athermal (there is no noise term) and the temperature $T$ that we vary only controls the Boltzmann distribution from which initial conditions for the dynamics are drawn.

We monitor the mean energy $\langle E(t) \rangle$ and the root mean squared velocity,

$$ \langle |\mathbf{v}(t)| \rangle = \left\langle \sqrt{\frac{1}{N} \sum_i \left| \frac{dr_i}{dt} \right|^2} \right\rangle, \qquad (2) $$

during the relaxation dynamics, where the brackets represent an average over initial equilibrium configurations, and $N$ is the number of particles. We define an exponent $\beta$ for the time decay [34] as

$$ \langle |\mathbf{v}(t)| \rangle \sim t^{-\beta}. \qquad (3) $$

For the dynamics in Eq. (1), the energy decay is exactly related to the velocity decay as $(1/N)[d/(dt)]\langle E(t) \rangle = -\zeta \langle |\mathbf{v}(t)|^2 \rangle$. As a result, the energy decay can be expressed using the same exponent: $\langle E(t) - E(t \to \infty) \rangle \sim t^{-(2\beta-1)}$. Therefore, we focus on the velocity relaxation and Eq. (3).

We consider several structural glass models in various dimensions and interaction potentials over a wide range of preparation temperatures. We study a soft-sphere version of the mean-field Mari-Kurchan model [44], polydisperse soft-sphere models in two [45] and three dimensions [46], harmonic spheres [47] in two, three, four, and eight dimensions, and the Kob-Andersen model [48] for two and three dimensions. Note that soft-sphere models have a steep repulsive interaction with an $r^{-12}$ core and a short cutoff (we have checked that extremely few rattler particles [49] are found in the corresponding inherent structures), whereas the harmonic potential models have a very soft core, which may affect the $T \to \infty$ limit for initial conditions. The Kob-Andersen model uses the Lennard-Jones potential with a steep repulsive core and attractive forces at larger distances.

To prepare equilibrium configurations in a wide range of temperatures, we use the planting method [50] for the Mari-Kurchan model and the swap Monte Carlo algorithm [46] for some of the finite-dimensional systems, which should allow us to detect any of the putative transitions predicted from mean-field landscapes. Further details about the models and simulation protocols are provided in the Appendix A and Supplemental Material (SM) [51].

### B. Mean-field Mari-Kurchan model

Thanks to its mean-field nature, we can apply the replica liquid theory to the Mari-Kurchan model, as detailed in SM [51], and obtain the dynamical mode-coupling transition temperature: $T_{MCT} \approx 0.0084$. We also studied the equilibrium dynamics using a simple Metropolis algorithm, and find that the theoretical estimate of $T_{MCT}$ describes the numerical data reasonably well. This study allows us to also estimate the onset temperature for slow dynamics: $T_{onset} \approx 0.015$.

We study the steepest descent starting from equilibrium configurations in the range $T \in [0.0001, 0.015]$. Figure 1 shows the velocity decay $\langle |\mathbf{v}(t)| \rangle$ for various temperatures and system sizes. Figure 1(a) shows that the relaxation dynamics strongly depends on initial equilibrium

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Velocity $\langle |\mathbf{v}(t)| \rangle$ as a function of time $t$ in the Mari-Kurchan soft-sphere model at several temperatures with $N = 16384$. (b) Velocity decay at two selected initial equilibrium temperatures and several system sizes. The dashed line indicates $\beta = 0.75$.</figcaption></figure>
temperature. For high temperatures, $T \gtrsim 0.006$, $\langle|v(t)|\rangle$ follows a clear power-law decay with an exponent that we estimate as $\beta \approx 0.75$. In a finite-size system, this power-law decay is interrupted at long times. On the other hand, at low temperatures, an exponential decay occurs. These results suggest that the high- and low-temperature relaxation dynamics are qualitatively different, and are separated by a critical temperature.

To fully confirm the distinct occurrence of power-law and exponential decays, we analyze finite-size effects. In Fig. 1(b), we show the velocity $\langle|v(t)|\rangle$ for several system sizes at two selected temperatures. At high initial temperature $T = 0.015$, $\langle|v|\rangle$ has a strong system-size dependence. Larger systems take longer times to reach energy minima and follow the power-law decay with $\beta \approx 0.75$ over a broader time window. The system-size dependence suggests that in the thermodynamic limit, $N \to \infty$, the velocity decay has a genuine power-law behavior with a diverging timescale. At very low initial temperature, $T = 0.0001$ instead, the velocity decay has almost no system-size dependence, implying that the time to reach energy minima remains finite in the thermodynamic limit, confirming the exponential decay.

Therefore, in the initial temperature regime between 0.0001 and 0.006, the Mari-Kurchan model displays a transition characterizing the nature of the relaxation dynamics akin to the behavior reported at the temperature $T_{\text{SF}}$ discussed in p-spin models. Interestingly, for the MK model, the $T_{\text{SF}}$ is noticeably smaller than the estimated

$T_{\text{MCT}} \approx 0.0084$, confirming that these two temperatures should be distinguished (they are equal in some versions of p-spin models).

While the existence of the transition is compatible with results for the mixed p-spin glass model, the decay exponent $\beta \approx 0.75$ at high temperatures observed for the Mari-Kurchan model differs slightly from the spin glass model where $\beta \approx 0.83$ [60], or the random Lorentz gas in $d \to \infty$ where $\beta \approx 1$ [40]. A broader class of mean-field models and even more extensive numerics should be explored to settle the relevance of this small difference.

### C. Finite-dimensional models

For finite $D$ systems, we first consider the relaxation dynamics starting from the high-temperature limit, $T \to \infty$, in various models and spatial dimensions. Figure 2(a) shows the results for monodisperse harmonic spheres in dimensions $D=2, 3, 4$, and $8$. In all dimensions, we observe a power-law decay, but the exponent $\beta$ depends on $D$. We find $\beta \approx 0.92$ for $D=2$ and $\beta \approx 0.85$ for $D=3$, which are consistent with previous work [34]. For $D=4$ and $8$, on the other hand, the results are described by the same exponent $\beta \approx 0.75$. This suggests that a mean-field value $\beta = 0.75$ in the high-temperature regime is reached for $D \ge 4$. Interestingly, this exponent is close to the one observed in the Mari-Kurchan model at high temperatures. In SM, we discuss the system-size dependence of $\langle|v(t)|\rangle$ [51]. Larger systems always take a longer time to reach

<figure><img src="image_4.png" /><figcaption>FIG. 2. The velocity $\langle|v(t)|\rangle$ as a function of time $t$ for the harmonic spheres in two, three, four, and eight dimensions (a), the three-dimensional soft spheres (b), the Kob-Andersen Lennard-Jones model (KALJ) (c), the two-dimensional harmonic spheres (d), the soft spheres (e), and the Kob-Andersen Lennard-Jones model (f). In (a), the number of particles is $N = 64000$ for $D=2$, $N = 65536$ for $D=3$ and $D=4$, and $N = 16384$ for $D=8$ (data are shifted vertically for clarity). In (b) and (d)–(f), $N = 96000$, 64000, 64000, and 125000, respectively. In (c), $N = 27135$ at $T = 0.37$ and $N = 76800$ for other temperatures.</figcaption></figure>
energy minima, which is consistent with a pure power-law decay in the thermodynamic limit $N \to \infty$ at large times.

To investigate universality, we look at the results at high temperatures for various models. Figures 2(b)–2(f) show the velocity decay for harmonic and soft spheres, and the Kob-Andersen model. Note that for soft spheres and Kob-Andersen models, the influence of the repulsive core can be felt at arbitrarily large temperatures. Nevertheless, all models at higher initial temperatures asymptotically show $\beta \approx 0.92$ and $0.85$ in $D=2$ and $D=3$, respectively. Therefore, we conclude that the value of $\beta$ is universal, irrespective of the details of the interaction potentials, size polydispersity, or the proximity of a jamming transition.

We then study the effect of the initial stability on the relaxation dynamics. We first consider $D=3$ polydisperse soft spheres. Using the swap Monte Carlo algorithm, we vary initial equilibrium temperature quite significantly, $T \in [0.062, 1.0]$, which includes $T_{MCT} \approx 0.1$ and $T_{onset} \approx 0.18$ determined by standard methods [46]. The model reproduces the universal exponent in $D=3$, $\beta = 0.85$, at finite but high temperatures; see Fig. 2(b). With decreasing temperature, the velocity relaxation $\langle |v(t)| \rangle$ becomes faster, as expected from the physical intuition that the system starts closer to an energy minimum in a smoother landscape. However, even at temperatures much below $T_{MCT}$, the velocity relaxation $\langle |v(t)| \rangle$ displays a power-law decay, but with a larger apparent exponent, $\beta \approx 1.25$. The same exponent at low $T$ is found in the $D=3$ Kob-Andersen model [see Fig. 2(c)], suggesting that $\beta \approx 1.25$ is also universal. We vary the initial stability for $D=2$ harmonic spheres, soft spheres, and the Kob-Andersen model in Figs. 2(d)–2(f). The data demonstrate the same trend as in $D=3$ models, yet the low-temperature velocity decay exponent is now $\beta \approx 1$ in all three models, different from the $D=3$ value. Importantly, we do not observe an exponential decay at any studied temperatures in any of the finite-dimensional models, in contrast to the mean-field spin glass and Mari-Kurchan models. Instead, we find that the high- and low-temperature regimes are both characterized by universal power laws, with an exponent which only depends on the spatial dimension.

### D. Harmonic limit

We can rationalize our numerical observations at low temperatures using a harmonic dynamical description [61]. At very low initial temperatures, the initial equilibrium configuration is located nearer to the final inherent state. Thus, it makes sense to approximate the energy during the steepest descent dynamics using a harmonic expansion,

$$E(t) \approx E(t \to \infty) + \frac{1}{2} \Delta \mathbf{r}(t) \cdot \mathbf{H} \cdot \Delta \mathbf{r}(t), \quad (4)$$

with $\Delta \mathbf{r}(t) = \mathbf{r}(t \to \infty) - \mathbf{r}(t)$ and $\mathbf{H}$ the Hessian matrix in the energy minimum. Let us assume that the phononic

modes following the Debye law and quasilocalized modes following the non-Debye quartic law coexist in the low-frequency region of the vibrational density of states [52]. By linearizing the equations of motion, we can relate the time decay of the velocity to the properties of the Hessian matrix and we find that the velocity should decay with an exponent $\beta_{harm} = D/4 + 1/2$, yielding $\beta_{harm} = 1$ and $1.25$ for $D=2$ and $D=3$, respectively. (See Appendix B for a more detailed discussion of the harmonic approximation, which also shows that quasilocalized modes provide a subdominant contribution to the velocity decay for $D < 5$.) These values are fully consistent with our numerical observations, which means that the low-temperature relaxation dynamics appears to be well described, at least over the simulated timescales and system sizes, by a simple harmonic approach, for all types of particle interactions.

### E. Localized defects

The harmonic analysis shows that, because of phonons, the mean-field transition at $T_{SF}$ to gapped energy minima reached exponentially fast cannot exist in finite $D$, and relaxation dynamics is in fact necessarily algebraic, even in a harmonic, “state following” limit. Our numerics is nevertheless compatible with two distinct temperature regimes, with nonharmonic effects becoming predominant at high initial $T$. Is a sharp transition separating these two regimes?

To address this question, we must understand the microscopic relaxation mechanism beyond the harmonic limit. At very low temperatures, we expect that the initial configuration and the final inherent state differ by small displacements which do not affect much the geometry of the particle packing. Figures 3(a)–3(c) show the displacement

<figure><img src="image_5.png" /><figcaption>FIG. 3. Two-dimensional soft spheres. (a)–(c) Nonaffine displacements $D_{min}^2$ and (d)–(f) defects for three different configurations. In (d)–(f), defects are shown in red, other particles are in blue. In all panels, the displacement vectors $\mathbf{r}_i(t=0) - \mathbf{r}_i(t \to \infty)$ (arrows) are amplified by a factor of 3. The temperature is $T = 0.035$ in (a) and (d), 0.2 in (b) and (e), and 0.8 in (c) and (f).</figcaption></figure>
and nonaffine displacement fields [62] between initial and final configurations for the $D=2$ soft-sphere system. For $T=0.035$, where the harmonic description works well (and $\beta = \beta_{\text{harm}} = 1.0$ is measured), particle displacements are indeed very small, implying that most particles interact with the same neighbors in initial and final states; see Fig. 3(a). For $T=0.2$, however, larger displacements are observed, and the most mobile particles are spatially correlated; see Fig. 3(b). Large nonaffine displacements are associated with localized particle rearrangements occurring during the steepest descent, which we call “defects.” As the temperature is increased further, more particles have large nonaffine displacements, see Fig. 3(c), and the initial and final configurations become substantially different.

To quantify particle rearrangements during minimization, we introduce a variable $\phi_i$ for each particle defined such that $\phi_i = 0$ if particle $i$ neither loses nor gains any neighbor during steepest descent, and $\phi_i = 1$ otherwise (see SM for precise definitions [51]), and denote $\phi = (1/N) \sum_i \phi_i$ the concentration of such defects in a given configuration with $N$ particles. The field $\phi_i$ thus identifies the location of particle rearrangements, as shown in Figs. 3(d)–3(f). Red particles with $\phi_i = 1$ are found in high $D_{\text{min}}^2$ regions, which validates the proposed identification of defects. The defects are also observed in the displacement field $|\mathbf{r}_i(t) - \mathbf{r}_i(\infty)|$; see SM for further discussion. In Ref. [34], similar defects were visualized using the nonaffine velocity field.

In Fig. 4, we show the average concentration of defects, $\langle\phi\rangle$, for $D=2$ and $D=3$ soft-sphere models, and the collective susceptibility of defects, $\chi = N(\langle\phi^2\rangle - \langle\phi\rangle^2)$. The average defect density is a smooth function of temperature which seems to remain finite at any initial $T > 0$. The susceptibility shows a well-defined peak, whose shape and location are independent of the system size; see Fig. 4(b). These results indicate that no sharp phase transition (with a vanishing $\langle\phi\rangle$) separates the relaxation dynamics between high and low temperatures. The defect density has a sigmoidal shape as it saturates to unity at large $T$ and decreases very rapidly to small values as $T \to 0$. It displays an inflection point at a temperature $T_{\text{def}}$ that also corresponds to the peak of the susceptibility. Physically, $T_{\text{def}}$ represents the temperature where $\langle\phi\rangle$ varies more strongly with $T$ and has the largest fluctuations, thus separating the high-$T$ regime where $\langle\phi\rangle$ approaches unity, from low $T$ where it is very small. The gradual disappearance of localized rearrangements presumably explains the temperature evolution of the self-part of the van Hove function [36]. The discussion of the harmonic limit in Sec. II D showed that the defects revealed by steepest descent dynamics at lower temperatures do not simply result from the harmonic excitation of the quasilocalized modes populating the low-frequency part of the density of states (which would lead to a different power-law decay), although a more complicated relation could exist.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Soft spheres. (a) Average concentration $\langle\phi\rangle$ and (b) collective susceptibility $\chi$ of defects, as a function of initial equilibrium temperature. The temperature evolution of the defect concentration is smooth, with a maximum variation near $T_{\text{def}} \approx 0.25$ and $T_{\text{def}} \approx 0.15$ for $D=2$ and 3, respectively.</figcaption></figure>

### III. DISCUSSION

We studied the physical dynamics during steepest descent energy minimization for various glass-forming models in spatial dimensions $D=2$ to $D=8$ and also in the mean-field limit, for a wide range of initial conditions. Focusing on the exponent $\beta$ characterizing the algebraic decay of the average velocity, we identified its universal, finite-dimensional features. First, we showed that the mean-field transition at temperature $T_{\text{SF}}$ to an exponential decay cannot exist in finite $D$ due to the presence of phonons. More importantly, we showed that the measured evolution of $\beta$ from its high-temperature universal value toward a larger harmonic value $\beta_{\text{harm}} = D/4 + 1/2$, observed at low temperatures, reflects in fact the gradual suppression of a population of localized defects with decreasing $T$. The relative importance of defects and plane waves explains the observed evolution of $\beta$. Since $\beta_{\text{harm}}$ is larger than its high-$T$ value, we expect the latter exponent to dominate the long-time limit of the velocity decay at any finite temperature in the thermodynamic limit. In this view, the harmonic regime is only a transient which lasts longer at lower temperature when there are less defects. As a result, the mean-field critical temperature $T_{\text{SF}}$ has no analog in finite $D$. This implies that, at finite temperature, an instantaneous configuration of a finite-dimensional glass-forming system can never be seen as an inherent structure
excited by small thermal fluctuations. It would be interesting to explore theoretical models alternative to mean-field glass models, such as elastoplastic models [63], to better account for our numerical observations, in particular, the value of the exponent $\beta$.

Our results have broad physical consequences. First, they imply that the defect dynamics leading to the coarsening of the nonaffine velocity field described in Ref. [34] (see also SM [51]) is actually relevant for generic finite D glass-forming liquids, and is unrelated to the athermal jamming transition. The observed universal exponent $\beta$ implies similarly universal geometrical features of the potential energy landscapes of generic structural glasses. Interestingly, an experimental realization of the steepest descent dynamics has recently been proposed [64]. By perturbing a stable foam configuration in two dimensions, localized defects during the relaxation were also observed. Such experiments could validate our numerical findings, especially the universal exponent $\beta$ found at high initial temperatures. More generally, our observations about defects at lower initial temperatures is another supporting evidence of the existence of localized excitations in stable glasses relevant for metallic and molecular glasses [65].

Second, together with recent analytic and numerical works [35,66], our results shed new light on the connection between equilibrium glassy dynamics and stationary points of the potential energy landscape. The interpretation of the mode-coupling temperature $T_{MCT}$ as a topographic change in the potential energy landscape does not hold in mixed $p$-spin models [35]. Our simulations of the Mari-Kurchan model confirm that the saddle-to-minima transition occurs at a temperature $T_{SF}$ distinct from $T_{MCT}$, already at mean-field level. The emergence of localized defects in finite dimensions found here is consistent with the recent conclusion [66] that the critical transition at $T_{SF}$ is replaced in finite D by a smooth crossover. A similar scenario controlled by noninteracting localized defects was also found in kinetically constrained models [42], thus suggesting a potential connection between the defects revealed by steepest descent dynamics and those discussed in the context of dynamic facilitation [67]. However, the power-law decay revealed by our study cannot result from the deexcitation of a noninteracting gas of isolated defects, and steepest descent dynamics in kinetically constrained models would instead be unremarkable. It is also unclear whether elastoplastic models where relaxation events are coupled by elasticity can account for our findings.

Third, our finding that a finite concentration of defects controls the nonharmonic relaxation from equilibrated configurations to inherent states suggests that the potential energy landscape of glass formers is both rugged and chaotic. To test this idea numerically, we applied a very small random perturbation to the initial configuration and monitored the subsequent steepest descent dynamics. We found that a slight perturbation typically leads to different

inherent structures (not shown), consistent with earlier work [68,69]. The strong chaoticity of the minimization dynamics implies that the energy minimum reached from a given equilibrium configuration in fact strongly depends on the minimization algorithm itself [70]. The steepest descent dynamics we used is just the simplest algorithm for numerical optimization, but there are several other (usually more efficient) ways to reach the bottom of the potential energy landscape, such as conjugate gradient (CG) [71] and fast inertial relaxation engine (FIRE) [53]. Indeed, we find that starting from the same initial configuration, the steepest descent, CG, and FIRE dynamics typically converge to different inherent structures, as quantified by their mutual distances (see SM [51]). The evolution of this distance mirrors the temperature evolution of the defect concentration in Fig. 4(a), higher initial temperatures leading to larger separations. In Fig. 5, we show representative snapshots of the displacement field between two inherent structures obtained by two different algorithms starting from a unique initial configuration. A single localized defect can be seen at low $T$, which naturally gives rise to a quadrupolar Eshelby-like displacement field. Defects proliferate at higher temperature. This shows that mapping an equilibrium liquid state to an inherent structure is a fully dynamical problem, which becomes uniquely defined only after a specific choice for the minimization algorithm is made. Localized defects, which had been used by Stillinger [72] to construct an argument against the existence of a Kauzmann transition (see the discussion in Ref. [73]), instead weaken the thermodynamic significance of a tiling of configuration space directly based on inherent states.

In recent years, localized glassy defects have been reported from the study of harmonic [31] and non-harmonic excitations, in the fields of plasticity [32] and low-temperature transport properties [30], and in connection with secondary relaxations in deeply supercooled liquids [74,75] and dynamic facilitation [67]. Steepest

<figure><img src="image_7.png" /><figcaption>FIG. 5. Displacement fields, $r_i^{CG} - r_i^{FIRE}$, between pairs of minima obtained via conjugate gradient ($\{r_i^{CG}\}$) and FIRE algorithms ($\{r_i^{FIRE}\}$), following a quench from the same initial configuration at $T = 0.07$ (a) and $T = 0.15$ (b), for 2D soft spheres with $N = 64000$. Arrows magnified by a factor 40 and 2 for (a) and (b), respectively.</figcaption></figure>
descent dynamics thus corresponds to another situation where localized excitations control structural rearrangements at the particle scale and reveal that they interact in a nontrivial manner. Future work should establish the similarities and differences between these disparate observations. Ultimately, we expect that a unifying picture of localized defects with specific interactions will soon become available and applicable to a host of different physical situations.

## ACKNOWLEDGMENTS

We thank G. Biroli, R. Chacko, G. Folena, and H. Ikeda for discussions. We also thank A. D. S. Parmar for sharing stable Kob-Andersen configurations. This work was supported by grants from the Simons Foundation (Grant No. 454935, L. B.) and JSPS KAKENHI (Grants No. 18H05225, No. 19H01812, No. 20H01868, and No. 20H00128, A. I.).

# APPENDIX A: METHODS

## 1. Models

We study the steepest descent dynamics of models with three different interaction potentials: soft spheres, harmonic spheres, and Lennard-Jones interactions. The dimensionality dependence, including the mean-field limit, of this dynamics is studied by using the models in two, three, four, eight dimensions, and the mean-field Mari-Kurchan model [44,76].

### a. Soft spheres

The two- and three-dimensional soft-sphere models [45,46] consist of particles with purely repulsive interactions and a continuous size polydispersity. Particle diameters $d_i$ are randomly drawn from a distribution of the form $f(d) = Ad^{-3}$, for $d \in [d_{min}, d_{max}]$, where A is a normalization constant. The size polydispersity is quantified by $\delta = (\bar{d}^2 - \bar{d}^2)^{1/2}/\bar{d}$, where the over line denotes an average over the distribution $f(d)$. Here we choose $\delta = 0.23$ by imposing $d_{min}/d_{max} = 0.449$. The average diameter $\bar{d}$ sets the unit of length. The soft-sphere interactions are pairwise and described by an inverse power-law potential,

$$ 
\begin{aligned}
U_{ij}(r) &= v_0 \left(\frac{d_{ij}}{r}\right)^{12} + c_0 + c_1 \left(\frac{r}{d_{ij}}\right)^2 + c_2 \left(\frac{r}{d_{ij}}\right)^4, \\
d_{ij} &= \frac{(d_i + d_j)}{2} (1 - \epsilon |d_i - d_j|),
\end{aligned} 
\quad (A1) 
$$

where $v_0$ sets the unit of energy (and of temperature with the Boltzmann constant $k_B \equiv 1$) and $\epsilon = 0.2$ quantifies the degree of nonadditivity of particle diameters. We introduce $\epsilon > 0$ in the model to suppress fractionation and thus to enhance the glass-forming ability. The constants $c_0$, $c_1$, and $c_2$ enforce a vanishing potential and continuity of its

first- and second-order derivatives at the cutoff distance $r_{cutoff} = 1.25d_{ij}$. We simulate a system with N particles within a square cell of area (volume) $V = L^2$ ($V = L^3$), where L is the linear box length, under periodic boundary conditions, at a number density $\rho = N/V = 1$ (1.02) for 2D (3D).

We prepare equilibrium configurations using the swap Monte Carlo algorithm [46]. With probability $P_{swap} = 0.2$, we perform a swap move where we randomly pick two particles (i and j) having similar diameters ($|d_i - d_j| < 0.2$) and attempt to exchange their diameters. With probability $1 - P_{swap} = 0.8$, instead, we perform conventional Monte Carlo translational moves, where we pick one particle and displace it within a box with linear length $\delta_{max} = 0.12\bar{d}$.

### b. Harmonic spheres

We study the harmonic sphere model [47,77] in two, three, four, and eight dimensions. The harmonic sphere model has an interaction potential,

$$ U_{ij}(r_{ij}) = \frac{v_0}{2} \left(1 - \frac{r_{ij}}{d_{ij}}\right)^2, \quad (A2) $$

$$ d_{ij} = \frac{(d_i + d_j)}{2} (1 - \epsilon |d_i - d_j|), \quad (A3) $$

where $v_0$ is again the unit of the energy scale. For the two-dimensional model, to avoid crystallization at low temperature, we use the continuously polydisperse nonadditive model with the same distribution of the particle diameters used in the soft-sphere model and $\epsilon = 0.2$ in two dimensions. The unit length scale for the two-dimensional model is $\bar{d}$ as well as the soft-sphere model. We again use the swap Monte Carlo algorithm with the same setting and parameters as for the polydisperse soft spheres to equilibrate down to very low temperatures. In three, four, and eight dimensions, crystallization is highly suppressed, and the simple additive ($\epsilon = 0$) monodisperse model is enough to study the relaxation dynamics to disordered states. Because of the finite range of the interaction, the system has a critical jamming transition at finite density, below which the relaxation dynamics shows an exponentially fast decay toward zero energy states [78]. Since in other models we study the relaxation dynamics toward energy minima with a finite energy, a direct comparison is possible when the inherent structures of harmonic spheres have finite energies as well. We thus set the volume fraction above the jamming transition to $\phi = 1.2, 0.73, 0.5, 0.1$ in two, three, four, and eight dimensions, respectively, so that the final energies are always finite.
### c. Kob-Andersen Lennard-Jones model

For the case of the well-studied Kob-Andersen binary Lennard-Jones (KALJ) model, the interaction between two particles has the following form:

$$U_{ij}(r_{ij}) = 4\epsilon_{ij}\left(\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right), \quad (A4)$$

for $r < r_{cutoff} = 2.5\sigma_{ij}$, and particles i, j can belong to either A or B species which constitute the binary mixture. $r_{cutoff}$ is the cutoff distance at which the potential $U_{ij}(r_{ij})$ is truncated. The different interaction parameters for the binary mixture take the following values: $\epsilon_{AA} = 1.0$, $\epsilon_{AB} = 1.5\epsilon_{AA}$, $\epsilon_{BB} = 0.5\epsilon_{AA}$; $\sigma_{AA} = 1.0$, $\sigma_{AB} = 0.8\sigma_{AA}$, $\sigma_{BB} = 0.88\sigma_{AA}$. The mixture has 80:20 composition in 3D and 65:35 composition in 2D, to optimize glass-forming ability. In $D = 2$, we study a system consisting of 125 000 particles, and for $D = 3$, we study a system of 76 800 particles. Additionally, in $D = 3$, we study a system of 27 135 particles to probe the quench dynamics of states sampled at a low temperature ($T = 0.37$), where configurations are obtained by a swap Monte Carlo scheme developed in Ref. [79].

### d. Mari-Kurchan model

We study the mean-field Mari-Kurchan (MK) model [44,76] in three dimensions with the simple monodisperse soft-sphere interaction in Eq. (A1) with $\epsilon = 0$ and the cutoff length $r_{cutoff} = 4d$, where $d$ is the diameter of particles. The volume fraction is $\phi = 0.5$. The MK model has quenched randomness in the particle distance, and the interaction potential is thus $U(|\mathbf{r}_i - \mathbf{r}_j + \mathbf{A}_{ij}|)$, where $\mathbf{A}_{ij}$ is a three-dimensional vector with each component sampled from the uniform distribution in the interval $[0, L]$ ($L$ the box size). Equilibrium configurations of the MK model are produced by using the planting technique [44,54]. For systems with general isotropic interactions, the cubic shape of the box complicates the direct sampling of the random shifts from the Boltzmann distribution:

$$P(\mathbf{A}_{ij}|\mathbf{r}_{ij}) = \frac{\exp[-\beta U_{ij}(|\mathbf{r}_i - \mathbf{r}_j + \mathbf{A}_{ij}|)]}{\int d\mathbf{A}_{ij} \exp[-\beta U_{ij}(|\mathbf{r}_i - \mathbf{r}_j + \mathbf{A}_{ij}|)]} \quad (A5)$$

We thus use the Markov chain Monte Carlo method to sample the random shifts $\{\mathbf{A}_{ij}\}$ from the distribution Eq. (A5) so that any given particle configuration follows from the Boltzmann distribution. For each pair of particles i and j, we take $\mathbf{A}_{ij}$ as the random shift after 200 Monte Carlo sweeps with the simple Metropolis algorithm starting from uniformly random numbers.

## APPENDIX B: HARMONIC EXPONENT

We discuss the asymptotic decay of the velocity by assuming that the system is perfectly harmonic and the

vibrational density of states follows the Debye law. Let the Hessian matrix $H$ of an inherent structure have eigenvalues $\{\lambda_a\}$ and corresponding eigenvectors $\{\mathbf{x}_a\}$. Since the Hessian matrix is real symmetric, eigenvectors are orthogonal; $\mathbf{x}_a \cdot \mathbf{x}_b = \delta_{ab}$, where $\delta_{ab}$ is the Kronecker's delta. Using the eigenvectors, we have the particle displacement written as

$$\Delta\mathbf{r}(t) = \sum_a c_a(t)\mathbf{x}_a, \quad (B1)$$

where $c_a(t) = \Delta\mathbf{r}(t) \cdot \mathbf{x}_a$. Suppose that the system is perfectly harmonic, i.e., the system follows linearized equations of motion:

$$\zeta \frac{d\Delta\mathbf{r}}{dt} = -H \cdot \Delta\mathbf{r}. \quad (B2)$$

Then each mode decays exponentially with $c_a(t) = c_a(0)e^{-\lambda_a t}$ and the equipartition law $\langle c_a(0)^2 \rangle = (T/\lambda_a)$ holds (The brackets in Appendix B represent an average over initial configurations with the same inherent structure).

In this harmonic approximation, the potential energy decreases with time as

$$\begin{aligned} \langle E(t) - E(t \to \infty) \rangle &= \frac{1}{2} \sum_a \langle c_a^2(0) \rangle \lambda_a \exp(-2\lambda_a t) \\ &= \frac{NT}{2} \int d\lambda \rho(\lambda) \exp(-2\lambda t), \end{aligned} \quad (B3)$$

where $\rho(\lambda) = 1/N \sum_a \delta(\lambda - \lambda_a)$ is the density of eigenvalues.

Let us assume that the density of state has the contributions from the phononic modes following the Debye law and quasilocalized modes following the non-Debye quartic law, i.e., $g(\omega) = A_0\omega^{D-1} + A_4\omega^4$ [52,80,81]. Then the density of eigenvalues reads $\rho(\lambda) = g(\omega)[(d\omega)/(d\lambda)] = A_0\lambda^{D/2-1}/2 + A_4\lambda^{3/2}/2$. Thus,

$$\begin{aligned} \langle E(t) - E(t \to \infty) \rangle &\sim \frac{NT}{4} \int d\lambda (A_0\lambda^{D/2-1} + A_4\lambda^{3/2})e^{-2\lambda t} \\ &\sim t^{-D/2} + O(t^{-5/2}). \end{aligned} \quad (B4)$$

Therefore, the energy relaxation is dominated by $t^{-D/2}$ when $D \le 5$. Since, for the steepest descent dynamics with the equations of motion given by Eq. (1), the energy decay can be related to the velocity decay, we finally obtain $\langle|\mathbf{v}(t)|\rangle \sim t^{-\beta_{\text{harm}}}$ with $\beta_{\text{harm}} = D/4 + 1/2$ for $D \le 5$.

[1] *Complex Systems*, Lecture Notes of the Les Houches Summer School, 2006, Session LXXXV. edited by J.-P.
Bouchaud, M. Mézard, and J. Dalibard (Elsevier, New York, 2007).

[2] D. Wales, *Energy Landscapes: Applications to Clusters, Biomolecules and Glasses*, Cambridge Molecular Science (Cambridge University Press, Cambridge, England, 2004).

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

[4] F. H. Stillinger, *A Topographic View of Supercooled Liquids and Glass Formation*, Science **267**, 1935 (1995).

[5] A. Auffinger, G. B. Arous, and J. Černý, *Random Matrices and Complexity of Spin Glasses*, Commun. Pure Appl. Math. **66**, 165 (2013).

[6] V. Ros, G. B. Arous, G. Biroli, and C. Cammarota, *Complex Energy Landscapes in Spiked-Tensor and Simple Glassy Models: Ruggedness, Arrangements of Local Minima, and Phase Transitions*, Phys. Rev. X **9**, 011003 (2019).

[7] F. Sciortino, *Potential Energy Landscape Description of Supercooled Liquids and Glasses*, J. Stat. Mech. (2005) P05015.

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

[9] F. Krzakała, A. Montanari, F. Ricci-Tersenghi, G. Semerjian, and L. Zdeborová, *Gibbs States and the Set of Solutions of Random Constraint Satisfaction Problems*, Proc. Natl. Acad. Sci. U.S.A. **104**, 10318 (2007).

[10] Y. LeCun, Y. Bengio, and G. Hinton, *Deep Learning*, Nature (London) **521**, 436 (2015).

[11] M. Baity-Jesi, L. Sagun, M. Geiger, S. Spigler, G. B. Arous, C. Cammarota, Y. LeCun, M. Wyart, and G. Biroli, *Comparing Dynamics: Deep Neural Networks versus Glassy Systems*, J. Stat. Mech. (2019) 124013.

[12] J. Kent-Dobias and J. Kurchan, *Complex Landscapes*, Phys. Rev. Research **3**, 023064 (2021).

[13] M. Goldstein, *Viscous Liquids and the Glass Transition: A Potential Energy Barrier Picture*, J. Chem. Phys. **51**, 3728 (1969).

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

[15] A. Cavagna, I. Giardina, and G. Parisi, *Stationary Points of the Thouless-Anderson-Palmer Free Energy*, Phys. Rev. B **57**, 11251 (1998).

[16] V. Ros, G. Biroli, and C. Cammarota, *Complexity of Energy Barriers in Mean-Field Glassy Systems*, Europhys. Lett. **126**, 20003 (2019).

[17] T. Rizzo, *Path Integral Approach Unveils Role of Complex Energy Landscape for Activated Dynamics of Glassy Systems*, Phys. Rev. B **104**, 094203 (2021).

[18] F. H. Stillinger and T. A. Weber, *Hidden Structure in Liquids*, Phys. Rev. A **25**, 978 (1982).

[19] L. Berthier and J. P. Garrahan, *Nontopographic Description of Inherent Structure Dynamics in Glassformers*, J. Chem. Phys. **119**, 4367 (2003).

[20] J. C. Dyre, *Colloquium: The Glass Transition and Elastic Models of Glass-Forming Liquids*, Rev. Mod. Phys. **78**, 953 (2006).

[21] G. Biroli and R. Monasson, *From Inherent Structures to Pure States: Some Simple Remarks and Examples*, Europhys. Lett. **50**, 155 (2000).

[22] A. Cavagna, *Fragile vs. Strong Liquids: A Saddles-Ruled Scenario*, Europhys. Lett. **53**, 490 (2001).

[23] L. Angelani, R. Di Leonardo, G. Ruocco, A. Scala, and F. Sciortino, *Saddles in the Energy Landscape Probed by Supercooled Liquids*, Phys. Rev. Lett. **85**, 5356 (2000).

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

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

[26] L. Berthier and D. Coslovich, *Novel Approach to Numerical Measurements of the Configurational Entropy in Supercooled Liquids*, Proc. Natl. Acad. Sci. U.S.A. **111**, 11668 (2014).

[27] M. Ozawa, A. Ikeda, K. Miyazaki, and W. Kob, *Ideal Glass States Are Not Purely Vibrational: Insight from Randomly Pinned Glasses*, Phys. Rev. Lett. **121**, 205501 (2018).

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

[29] S. Gelin, H. Tanaka, and A. Lemaître, *Anomalous Phonon Scattering and Elastic Correlations in Amorphous Solids*, Nat. Mater. **15**, 1177 (2016).

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

[31] E. Lerner, G. Düring, and E. Bouchbinder, *Statistics and Properties of Low-Frequency Vibrational Modes in Structural Glasses*, Phys. Rev. Lett. **117**, 035501 (2016).

[32] D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J.-L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner, and M. L. Manning, *Predicting Plasticity in Disordered Solids from Structural Indicators*, Phys. Rev. Mater. **4**, 113609 (2020).

[33] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in Fortran 90* (Cambridge University Press, Cambridge, England, 1996).

[34] R. N. Chacko, P. Sollich, and S. M. Fielding, *Slow Coarsening in Jammed Athermal Soft Particle Suspensions*, Phys. Rev. Lett. **123**, 108001 (2019).

[35] G. Folena, S. Franz, and F. Ricci-Tersenghi, *Rethinking Mean-Field Glassy Dynamics and Its Relation with the Energy Landscape: The Surprising Case of the Spherical Mixed p-Spin Model*, Phys. Rev. X **10**, 031045 (2020).

[36] K. González-López and E. Lerner, *An Energy-Landscape-Based Crossover Temperature in Glass-Forming Liquids*, J. Chem. Phys. **153**, 241101 (2020).

[37] P. Charbonneau and P. K. Morse, *Memory Formation in Jammed Hard Spheres*, Phys. Rev. Lett. **126**, 088001 (2021).

[38] G. Folena, S. Franz, and F. Ricci-Tersenghi, *Gradient Descent Dynamics in the Mixed p-Spin Spherical Model:*
*Finite-Size Simulations and Comparison with Mean-Field Integration*, J. Stat. Mech. (2021) 033302.

[39] E. Stanifer and M. L. Manning, *Avalanche Dynamics in Sheared Athermal Particle Packings Occurs via Localized Bursts Predicted by Unstable Linear Response*, arXiv: 2110.02803.

[40] A. Manacorda and F. Zamponi, *Gradient Descent Dynamics and the Jamming Transition in Infinite Dimensions*, arXiv:2201.01161.

[41] J. Kurchan and L. Laloux, *Phase Space Geometry and Slow Dynamics*, J. Phys. A **29**, 1929 (1996).

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

[43] P. Olsson, *Relaxation Times and Rheology in Dense Athermal Suspensions*, Phys. Rev. E **91**, 062209 (2015).

[44] R. Mari and J. Kurchan, *Dynamical Transition of Glasses: From Exact to Approximate*, J. Chem. Phys. **135**, 124504 (2011).

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

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

[47] L. Berthier and T. A. Witten, *Glass Transition of Dense Fluids of Hard and Compressible Spheres*, Phys. Rev. E **80**, 021502 (2009).

[48] W. Kob and H. C. Andersen, *Testing Mode-Coupling Theory for a Supercooled Binary Lennard-Jones Mixture. II. Intermediate Scattering Function and Dynamic Susceptibility*, Phys. Rev. E **52**, 4134 (1995).

[49] F. Giacco, L. de Arcangelis, M. P. Ciamarra, and E. Lippiello, *Rattler-Induced Aging Dynamics in Jammed Granular Systems*, Soft Matter **13**, 9132 (2017).

[50] F. Krzakala and L. Zdeborová, *Hiding Quiet Solutions in Random Constraint Satisfaction Problems*, Phys. Rev. Lett. **102**, 238701 (2009).

[51] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevX.12.021001 for further discussions on the replica liquid theory of the MK model, velocity decay, and defects in real space, which includes Refs. [34,52–59].

[52] H. Mizuno, H. Shiba, and A. Ikeda, *Continuum Limit of the Vibrational Properties of Amorphous Solids*, Proc. Natl. Acad. Sci. U.S.A. **114**, E9767 (2017).

[53] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, *Structural Relaxation Made Simple*, Phys. Rev. Lett. **97**, 170201 (2006).

[54] P. Charbonneau, Y. Jin, G. Parisi, and F. Zamponi, *Hopping and the Stokes–Einstein Relation Breakdown in Simple Glass Formers*, Proc. Natl. Acad. Sci. U.S.A. **111**, 15025 (2014).

[55] H. Ikeda, F. Zamponi, and A. Ikeda, *Mean Field Theory of the Swap Monte Carlo Algorithm*, J. Chem. Phys. **147**, 234506 (2017).

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

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

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

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

[60] The exponent $\beta$ for the spin glass model is obtained from the exponent for the energy decay reported in Ref. [35] which is given by $2\beta - 1$.

[61] The term “harmonic” here means that the system energy is expanded up to the quadratic term. This harmonic approximation is applicable to any smooth interaction potential, and not only to the system of harmonic spheres.

[62] M. L. Falk and J. S. Langer, *Dynamics of Viscoplastic Deformation in Amorphous Solids*, Phys. Rev. E **57**, 7192 (1998).

[63] J. T. Parley, S. M. Fielding, and P. Sollich, *Aging in a Mean Field Elastoplastic Model of Amorphous Solids*, Phys. Fluids **32**, 127104 (2020).

[64] N. Yanagisawa and R. Kurita, *Size Distribution Dependence of Collective Relaxation Dynamics in a Two-Dimensional Wet Foam*, Sci. Rep. **11**, 2786 (2021).

[65] J. C. Qiao, Q. Wang, J. M. Pelletier, H. Kato, R. Casalini, D. Crespo, E. Pineda, Y. Yao, and Y. Yang, *Structural Heterogeneities and Mechanical Behavior of Amorphous Alloys*, Prog. Mater. Sci. **104**, 250 (2019).

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

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

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

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

[70] L. Angelani, G. Ruocco, M. Sampoli, and F. Sciortino, *General Features of the Energy Landscape in Lennard-Jones-like Model Liquids*, J. Chem. Phys. **119**, 2120 (2003).

[71] J. Nocedal and S. J. Wright, *Numerical Optimization*, Springer Series in Operations Research and Financial Engineering (Springer, New York, 2006).

[72] F. H. Stillinger, *Supercooled Liquids, Glass Transitions, and the Kauzmann Paradox*, J. Chem. Phys. **88**, 7818 (1988).

[73] L. Berthier, M. Ozawa, and C. Scalliet, *Configurational Entropy of Glass-Forming Liquids*, J. Chem. Phys. **150**, 160902 (2019).

[74] H.-B. Yu, R. Richert, and K. Samwer, *Structural Rearrangements Governing Johari-Goldstein Relaxations in Metallic Glasses*, Sci. Adv. **3**, e1701577 (2017).
[75] B. Guiselin, C. Scalliet, and L. Berthier, *Microscopic Origin of Excess Wings in Relaxation Spectra of Deeply Supercooled Liquids*, arXiv:2103.01569.

[76] R. H. Kraichnan, *J. Math. Phys. (N.Y.)* **3**, 475 (1962).

[77] C. S. O'Hern, L. E. Silbert, J. Liu, and S. R. Nagel, *Jamming at Zero Temperature and Zero Applied Stress: The Epitome of Disorder*, Phys. Rev. E **68**, 011306 (2003).

[78] Y. Nishikawa, A. Ikeda, and L. Berthier, *Relaxation Dynamics of Non-Brownian Spheres Below Jamming*, J. Stat. Phys. **182**, 37 (2021).

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

[80] G. Kapteijns, E. Bouchbinder, and E. Lerner, *Universal Nonphononic Density of States in 2D, 3D, and 4D Glasses*, Phys. Rev. Lett. **121**, 055501 (2018).

[81] M. Shimada, H. Mizuno, L. Berthier, and A. Ikeda, *Low-Frequency Vibrations of Jammed Packings in Large Spatial Dimensions*, Phys. Rev. E **101**, 052906 (2020).