Relaxation dynamics in the energy landscape of glass-forming liquids

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 non-trivial 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 localised glassy excitations. Our work unifies several recent lines of studies aiming at a detailed characterisation of the complex potential energy landscape of glass-formers, and challenges both mean-field and real space descriptions of glasses.


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 characterisation and dynamic exploration of complex potential energy landscapes are important problems for amorphous materials [4,7,8], optimisation 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 fullyconnected disordered spin models, such as p-spin models.In this limit, the phase space can be divided into longlived metastable states (or, pure states), and both freeenergy 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 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][23][24][25].However, these approaches do not have the same level of rigour 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], 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 few studies have been dedicated to the physical processes at play during the energy minimization itself [34][35][36][37][38][39][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, 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 non-interacting 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 onset , and it decreases with decreasing T when T ≤ T onset .This sharp onset temperature does not affect the relaxation dynamics itself which obeys a nontrivial power-law time dependence as long as T ≥ T SF .By contrast, the decay is exponentially fast below T SF (initials stand for 'State Following').This implies that the system is always close to an energy minimum for T ≤ T SF in which it converges very quickly by steepest descent.The critical temperature T SF < T onset also reflects a change in the structure of the potential energy landscape, as inherent states have a marginal density of states above T SF , which becomes gapped below.Within p-spin models, these two characteristic temperatures are unrelated to the equilibrium dynamics, which becomes non-ergodic at the mode-coupling temperature T MCT , distinct from both T SF and T onset , 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 non-trivial, 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-from low-temperature behaviours.These interesting studies do not provide a complete physical picture of the relaxation dynamics towards energy minima, neither do they assess the existence of the critical temperature T SF found in mean-field approaches.The universality of the powerlaw 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 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 SF 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 non-trivial limits that we analyse in detail.We show that this crossover is controlled by a finite population of localised defects where particle rearrangements take place during the minimisation, with the overall concentration of these defects decreasing at lower temperatures.Therefore, finite dimensional glassforming systems at finite temperatures cannot be seen as inherent structures excited by small thermal fluctuations, since they are neither described by mean-field energy landscapes nor by a simple picture of non-interacting localised defects.Our results provide a complete picture of the relaxation dynamics in glassy landscapes, and illuminate the role, nature and interactions of localised defects in finite-dimensional structural glasses.

A. Steepest descent dynamics
We numerically solve the equations of motion of steepest descent dynamics, starting at time t = 0 from an equilibrium configuration prepared at initial temperature T , where ζ is the damping coefficient and E is the potential energy.The time unit is where 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 E(t) and the root mean squared velocity, 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 β for the time decay [34] as For the dynamics in Eq. ( 1), the energy decay is exactly related to the velocity decay as As a result, the energy decay can be expressed using the same exponent: 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 → ∞ 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 Appendix A and Supplemental Material [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 SI, and obtain the dynamical mode-coupling transition temperature: T MCT 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 0.015.
We study the steepest descent starting from equilibrium configurations in the range T ∈ [0.0001, 0.015].Figure 1 shows the velocity decay |v(t)| for various temperatures and system sizes.Figure 1(a) shows that the relaxation dynamics strongly depends on initial equilibrium temperature.For high temperatures, T 0.006, |v(t)| follows a clear power-law decay with an exponent that we estimate as β 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 analyse finite-size effects.In Fig. 1(b), we show the velocity |v(t)| for several system sizes at two selected temperatures.At high initial temperature T = 0.015, |v| has a strong system-size dependence.Larger systems take longer times to reach energy minima and follow the power-law decay with β 0.75 over a broader time window.The system-size dependence suggests that in the thermodynamic limit, N → ∞, the velocity decay has a genuine power-law behaviour with a diverging time scale.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 characterising the nature of the relaxation dynamics akin to the behaviour reported at the temperature T SF discussed in p-spin models.Interestingly, for the MK model, the T SF is noticeably smaller than the estimated T MCT ≈ 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 β 0.75 at high temperatures observed for the Mari-Kurchan model differs slightly from the spin glass model where β 0.83 [52], or the random Lorentz gas in d → ∞ where β ≈ 1 [40].A broader class of meanfield 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 → ∞, 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 β depends on D. We find β 0.92 for D = 2 and β 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 β 0.75.This suggests that a meanfield value β = 0.75 in the high temperature regime is reached for D ≥ 4. Interestingly, this exponent is close to the one observed in the Mari-Kurchan model at high temperatures.In SI, we discuss the system-size dependence of |v(t)| .Larger systems always take a longer time to reach energy minima, consistently with a pure power law decay in the thermodynamic limit N → ∞ at large times.
To investigate universality, we look at the results at high temperatures for various models.Figures 2(b-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 β 0.92 and 0.85 in D = 2 and D = 3, respectively.Therefore, we conclude that the value of β 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 ∈ [0.062, 1.0], which includes T MCT 0.1 and T onset 0.18 determined by standard methods [46].The model reproduces the universal exponent in D = 3, β = 0.85, at finite but high temperatures, see Fig. 2(b).With decreasing temperature, the velocity relaxation |v(t)| 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 |v(t)| displays a power-law decay, but with a larger apparent exponent, β 1.25.The same exponent at low T is found in the D = 3 Kob-Andersen model [see Fig. 2(c)], suggesting that β 1.25 is also universal.We vary the initial stability for D = 2 harmonic spheres, soft spheres, and the Kob-Andersen model in Fig. 2(d), (e), and (f).The data demonstrate the same trend as in D = 3 models, yet the low-temperature velocity decay exponent is now β 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 finitedimensional 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 characterised 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 [53] 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 with ∆r(t) = r(t → ∞) − r(t) and H the Hessian matrix in the energy minimum.Let us assume that the phononic modes following the Debye law and quasi-localised modes following the non-Debye quartic law coexist in the lowfrequency region of the vibrational density of states [54].
By linearising 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 β harm = D/4 + 1/2, yielding β harm = 1 and 1.25 for D = 2 and D = 3, respectively.(See Appendix A 1 d for a more detailed discussion of the harmonic approximation, which also shows that quasilocalised 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. Localised 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 non-harmonic effects becoming predomi- nant 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.Figure 3(a-c) show the displacement and non-affine displacement fields [55] between initial and final configurations for the D = 2 soft sphere system.For T = 0.035, where the harmonic description works well (and β = β harm = 1.0 is measured), particle displacements are indeed very small, implying that most particles interact with the same neighbours 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 localised particle rearrangements occurring during the steepest descent, which we call 'defects'.As the temperature is increased further, more particles have large non-affine displacements, see Fig. 3(c), and the initial and final configurations become substantially different.
To quantify particle rearrangements during minimisation, we introduce a variable φ i for each particle defined such that φ i = 0 if particle i neither loses nor gains any neighbour during steepest descent, and φ i = 1 otherwise (see SI for precise definitions), and denote φ = 1 N i φ i the concentration of such defects in a given configuration with N particles.The field φ i thus identifies the location of particle rearrangements, as shown in Fig. 3(d-f).Red particles with φ i = 1 are found in high D 2 min regions, which validates the proposed identification of defects.The defects are also observed in the displacement field |r i (t) − r i (∞)|, see SI for further discussion.In Ref. [34], similar defects were visualised using the nonaffine velocity field.In Fig. 4, we show the average concentration of defects, φ , for D = 2 and D = 3 soft-sphere models, and the collective susceptibility of defects, χ = N ( φ 2 − φ 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 φ ) 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 → 0. It displays an inflection point at a temperature T def that also corresponds to the peak of the susceptibility.Physically, T def represents the temperature where φ varies more strongly with T and has the largest fluctuations, thus separating the high-T regime where φ approaches unity, from low-T where it is very small.The gradual disappearance of localised 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 quasi-localised 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.

III. DISCUSSION
We studied the physical dynamics during steepest descent energy minimisation 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 β characterising 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 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 β from its high-temperature universal value towards a larger harmonic value β harm = D/4+1/2, observed at low-temperatures, reflects in fact the gradual suppression of a population of localised defects with decreasing T .The relative importance of defects and plane waves explains the observed evolution of β.Since β 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 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 elasto-plastic models [56], to account better for our numerical observations, in particular the value of the exponent β.
Our results have broad physical consequences.First, they imply that the defect dynamics leading to the coarsening of the non-affine velocity field described in Ref. [34] (see also SI) is actually relevant for generic finite D glassforming liquids, and is unrelated to the athermal jamming transition.The observed universal exponent β implies similarly universal geometrical features of the potential energy landscapes of generic structural glasses.Interestingly, an experimental realisation of the steepest descent dynamics has recently been proposed [57].By perturbing a stable foam configuration in two dimensions, localised defects during the relaxation were also observed.Such experiments could validate our numerical findings, especially the universal exponent β found at high initial temperatures.More generally, our observations about defects at lower initial temperatures is another supporting evidence of the existence of localised excitations in stable glasses relevant for metallic and molecular glasses [58].
Second, together with recent analytic and numerical works [35,59], 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 localised defects in finite dimensions found here is consistent with the recent conclusion [59] that the critical transition at T SF is replaced in finite D by a smooth crossover.A similar scenario controlled by noninteracting localised 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 [60].However, the power law decay revealed by our study cannot result from the de-excitation of a non-interacting gas of isolated defects and steepest descent dynamics in kinetically constrained models would instead be unremarkable.It is also unclear whether elasto-plastic models where relaxation events are coupled by elasticity can account for our findings.
Third, our finding that a finite concentration of defects controls the non-harmonic 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 [61,62].The strong chaoticity of the minimisation dynamics implies that the energy minimum reached from a given equilibrium configuration in fact strongly depends on the minimisation algorithm itself [63].The steepest descent (SD) dynamics we used is just the simplest algorithm for numerical optimisation, but there are several other (usually more efficient) ways to reach the bottom of the potential energy landscape, such as conjugate-gradient (CG) [64] and fast inertial relaxation engine (FIRE) [65].Indeed, we find that starting from the same initial configuration, the SD, CG and FIRE dynamics typically converge to different inherent structures, as quantified by their mutual distances (see SI).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 localised defects 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 minimisation algorithm is made.Localised defects, which had been used by Stillinger [66] to construct an argument against the existence of a Kauzmann transition (see the discussion in [67]), instead weaken the thermodynamic significance of a tiling of configuration space directly based on inherent states.
In recent years, localised glassy defects have been reported from the study of harmonic [31] and nonharmonic excitations, in the fields of plasticity [32] and low-temperature transport properties [30], and in connection with secondary relaxations in deeply supercooled liquids [68,69] and dynamic facilitation [60].Steepest descent dynamics thus corresponds to another situation where localised excitations control structural rearrangements at the particle scale and reveal that they interact in a non-trivial manner.Future work should establish the similarities and differences between these disparate observations.Ultimately, we expect that a unifying picture of localised defects with specific interactions will soon become available and applicable to a host of different physical situations.monic 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,70].

a. Soft spheres
The two-and three-dimensional soft sphere models [45,46] where v 0 sets the unit of energy (and of temperature with the Boltzmann constant k B ≡ 1) and = 0.2 quantifies the degree of nonadditivity of particle diameters.We introduce > 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 cut-off 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 ρ = 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 δ max = 0.12d.

b. Harmonic spheres
We study the harmonic sphere model [47,71] in two, three, four, and eight dimensions.The harmonic sphere model has an interaction potential where v 0 is again the unit of the energy scale.For the two dimensional model, to avoid crystallisation 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 = 0.2 in two dimensions.The unit length scale for the twodimensional model is 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, crystallisation is highly suppressed, and the simple additive ( = 0) monodisperse model is enough to study the relaxation dynamics to disordered states.Due to 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 towards zero energy states [72].
Since in other models we study the relaxation dynamics towards 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 φ = 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: for r < r cutoff = 2.5σ 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: AA = 1.0,AB = 1.5 AA , BB = 0.5 AA ; σ AA = 1.0, σ AB = 0.8σ AA , σ BB = 0.88σ AA .The mixture has 80:20 composition in 3D and 65:35 composition in 2D, to optimise glass-forming ability.In D = 2, we study a system consisting 125000 particles and for D = 3, we study a system of 76800 particles.Additionally, in D = 3, we study a system of 27135 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. [73].

d. Mari-Kurchan model
We study the mean-field Mari-Kurchan (MK) model [44,70] in three dimensions with the simple monodisperse soft-sphere interaction in Eq. (A1) with = 0 and the cutoff length r cutoff = 4d, where d is the diameter of particles.The volume fraction is φ = 0.5.The MK model has quenched randomness in the particle distance, and the interaction potential is thus U (|r i − r j + A ij |), where 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,74].For systems with general isotropic interactions, the cubic shape of the box complicates the direct sampling of the random shifts from the Boltzmann distribution .
(A5) We thus use the Markov chain Monte Carlo method to sample the random shifts {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 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 {λ a } and corresponding eigenvectors {x a }.Since the Hessian matrix is real symmetric, eigenvectors are orthogonal; x a • x b = δ ab , where δ ab is the Kronecker's delta.Using the eigenvectors, we have the particle dis-placement written as where c a (t) = ∆r(t)•x a .Suppose that the system is perfectly harmonic, i.e. the system follows linearised equations of motion, Then each mode decays exponentially with c a (t) = c a (0)e −λat and the equipartition law c a (0) 2 = T λa holds.
In this harmonic approximation, the potential energy decreases with time as 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 quasi-localised modes following the non-Debye quartic law i.e. g(ω) = A 0 ω D−1 + A 4 ω 4 [54,75,76].Then the density of eigenvalues reads ρ Therefore, the energy relaxation is dominated by t −D/2 when D ≤ 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 |v(t)| ∼ t −β harm with β harm = D/4 + 1/2 for D ≤ 5.
configurations.Let us first imagine a perfect crystal quenched from a small, finite temperature to zero temperature.The two configurations before and after the quench should be essentially the same except for small thermal fluctuations.In this ideal situation, a given particle neither loses nor gains any neighbour during the quench dynamics.In that case, the defect order parameter should produce zero for all particles.Instead, for a disordered liquid at finite temperature, the steepest descent dynamics involves collective rearrangements and a reshuffling of neighboring particles.When a given particle loses some of the neighboring particles, this process is referred to as a "bond-breaking", as studied extensively [S3, S4].We then assign a positive value to the defect order parameter for a particle involved in a bond-breaking process.We also expect that particles can gain new neighbours during the dynamics.Thus, a defect order parameter should also take into account this process which we call a "bond-insertion".These two processes are sketched in Fig. S2.
To define neighboring particles in practice, we monitor the radial distribution functions in Fig. S3.We compute a radial distribution function as a function of a normalized distance, x ij = r ij /d ij , which is a suitable representation for continuously polydisperse systems [S5].In this 2D polydisperse model, g(x) for the equilibrium state changes quite a lot because the swap Monte Carlo allows us to sample an extremely wide range of temperatures.Consequently, the first local minimum of g(x), x EQ min , varies x EQ min ≈ 1.34 − 1.48 with varying temperature, T = 0.035 − 0.800.On the other hand, g(x) for the inherent structure barely changes except for the amplitude of the first peak.In particular, the first local minimum is located in a nearly flat region (x ≈ 1.23 − 1.50).
Based on the observation of g(x), we consider the definition of neighbors for the bond-breaking process.For each particle, we define the neighboring particles in the initial equilibrium configurations as particles located within a cutoff, x EQ cut .We set x EQ cut = x EQ min .Similarly, we define neighboring particles in the inherent structure by introducing a cutoff x IS cut .Ideally, one would like to use again x IS cut = x EQ min , but then one finds that some neighbors have a displacement which is just above that threshold, and using the same x EQ min leads to false neighbor changes.Therefore, one needs to use x IS cut > x EQ cut to remove the false positives (see Fig. S2(a)), as already discussed in several previous studies [S3, S4].Thus, we set x IS cut = 1.5.Next, we consider the bond-insertion process.This can be viewed as a bond-breaking process for the time-reversal of the minimisation dynamics.Applying the same reasoning as above, we impose the condition x IS cut < x EQ cut to remove the false positives.Thus we set x IS cut = 1.23.Once neighbors are defined, we determine the list of neighbors for each particle in both equilibrium (L eq ) and inherent structure (L IS ) configurations: L eq = {l eq } and L IS = {l IS }, where l eq and l IS denote the identities of the neighbours.For each particle, say i, we define the number of bond breakings, B i by counting particles with l eq / ∈ L IS .We also count the number of bond insertions, I i , which is the number of particles with l IS / ∈ L EQ .We note that our choice of cutoff produces nearly identical mean values for B = 1 N i B i and I = 1 N i I i at all temperatures, meaning that the bond-breaking and bond-insertion processes take place with the same frequency, as expected.
Finally we define the defect order parameter φ i as follows: φ i = 1 if (B i + I i ) > 0 and φ i = 0 if (B i + I i ) = 0.It is useful to make the variable φ i binary, as it allows us to determine the absolute concentration of defects, defined as φ = 1 N i φ i .In the main text, we report the behaviour of the average defect density, φ , and of the corresponding collective susceptibility: χ = N ( φ 2 − φ 2 ).

IV. DISPLACEMENT FIELD
In Fig. S4, we show the displacement field |r i (t) − r i (t → ∞)| at temperatures T = 0.2 and 0.035 and times t = 10 −1 , 10 0 , 10 1 , and 10 2 .These two temperatures belong to the high-and low-temperature regimes, respectively (see Fig. 2(e)).At high temperature T = 0.2, the displacement field reveals many defects at short times.At longer timescales, the number of localised defects decreases, but some of them survive for a very long time.When a defect disappears, particles around the defect rearrange and the velocity |v(t)| has a large sudden increase.The time-evolution of the displacement field at high temperature is consistent with the coarsening picture studied in the T → ∞ limit for harmonic spheres in Ref. [S6].
At lower temperature T = 0.0035, on the other hand, the displacement field is much smoother even at short times, and we hardly see localised defects.This smooth displacement field is consistent with a harmonic description, as shown in Sec.II C.

A. Inherent structures from different algorithms
We show that localised defects affect the mapping between an equilibrium liquid microstate and a 'corresponding' minimum in the energy landscape, i.e. an inherent structure.Starting from the same initial equilibrium state, we minimize the energy using different algorithms available in LAMMPS [S8], namely conjugate gradient, FIRE [S7], and steepest descent.This exercise is done for the 2D system of polydisperse soft spheres.In all cases, we use the same tolerance thresholds (energy tolerance of 10 −16 and force tolerance of 10 −20 x cut EQ  x cut IS FIG.S2.Bond-breaking (a) and bond-insertion (b) processes.Schematic pictures show a particle (particle 0) loses a bond (with particle 4) during the minimisation (a), or gains a new bond (with particle 8) (b).Neigbors are defined as a pair of particles closer than a cutoff x EQ cut (which depends on T ) and x IS cut (which depends on the dynamic process under study).down to the same minimum, independent of the algorithm used and the mapping between an equilibrium liquid state and a local minimum in the energy landscape becomes unique.However, when increasing the system size and/or the temperature, these minima are always distinct and we suspect that this is always the case in the thermodynamic limit at any finite temperature.The difference between pairs of minima obtained from a given initial configuration is visualized in the main text.

FIG. 1 .
FIG. 1.(a) Velocity |v(t)| 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 β = 0.75.

FIG. 3 .
FIG. 3. 2D soft spheres.(a-c) Non-affine displacements D 2 min and (d-f) defects for three different configurations.In (df), defects are shown in red, other particles are in blue.In all panels, the displacement vectors ri(t = 0) − ri(t → ∞) (arrows) is amplified by a factor of 3. The temperature is T = 0.035 in (a, d), 0.2 in (b, e), and 0.8 in (c, f).

FIG. 4 .
FIG. 4. Soft spheres.(a) Average concentration φ and (b) collective susceptibility χ of defects, as a function of initial equilibrium temperature.The temperature evolution of the defect concentration is smooth, with a maximum variation near T def ≈ 0.25 and T def ≈ 0.15 for D = 2 and 3, respectively.

FIG. 5 .
FIG. 5. Displacement fields, r CG i − r FIRE i , between pairs of minima obtained via conjugate gradient ({r CG i }) and FIRE algorithms ({r FIRE i }), 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.
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 ∈ [d min , d max ], where A is a normalization constant.The size polydispersity is quantified by δ = (d 2 − d 2 ) 1/2 /d, where the overline denotes an average over the distribution f (d).Here we choose δ = 0.23 by imposing d min /d max = 0.449.The average diameter, d, sets the unit of length.The softsphere interactions are pairwise and described by an inverse power-law potential