Rethinking mean-field glassy dynamics and its relation with the energy landscape: the awkward case of the spherical mixed p-spin model

TL;DR

This paper investigates the out-of-equilibrium dynamics of spherical mixed p-spin models, revealing a novel dynamical phase transition and complex aging behavior relevant to glass physics and machine learning problems.

Contribution

It uncovers a new dynamical phase transition in mixed p-spin models and analyzes the aging dynamics and energy landscape connection, extending understanding beyond pure models.

Findings

Discovery of a dynamical phase transition above the dynamical transition temperature.

Asymptotic energy drops below the threshold energy of marginal minima.

Identification of complex aging behavior and deep marginal minima relaxation.

Abstract

The spherical p-spin model is not only a fundamental model in statistical mechanics of disordered system, but has recently gained popularity since many hard problems in machine learning can be mapped on it. Thus the study of the out of equilibrium dynamics in this model is interesting both for the glass physics and for its implications on algorithms solving NP-hard problems. We revisit the long-time limit of the out of equilibrium dynamics of mean-field spherical mixed p-spin models. We consider quenches (gradient descent dynamics) starting from initial conditions thermalized at some temperature in the ergodic phase. We perform numerical integration of the dynamical mean-field equations of the model and we find an unexpected dynamical phase transition. Below an onset temperature, higher than the dynamical transition temperature, the asymptotic energy goes below the "threshold energy" of the dominant marginal minima of the energy function and memory of the initial condition is kept. This behavior, not present in the pure spherical p-spin model, resembles closely the one observed in simulations of glass-forming liquids. We then investigate the nature of the asymptotic dynamics, finding an aging solution that relaxes towards deep marginal minima, evolving on a restricted marginal manifold. Careful analysis, however, rules out simple aging solutions. We compute the constrained complexity in the aim of connecting the asymptotic solution to the energy landscape.

Figures

Equations

$$H_J[\sigma ]=-\sum _p\sqrt{a_p}\sum _{i_1<...<i_p}^N{J}_{i_1...i_p}^{(p)}{\sigma }_{i_1}...{\sigma }_{i_p}$$

$$\overline{H[\sigma ]H[\tau ]}=Nf(q_{\sigma \tau })\text{where}$$

$$f(q)=\frac{1}{2}\sum _p{a}_p{q}^p\text{and}{q}_{\sigma \tau }=\frac{1}{N}\sum _{i=1}^N{\sigma }_i{\tau }_i.$$

$$\frac{\partial H[\sigma ]}{\partial {\sigma }_i}+\mu {\sigma }_i=H_i^{\prime }+\mu {\sigma }_i=0.$$

$$\mu =-\frac{1}{N}\sum _i{\sigma }_i{H}_i^{\prime }.$$

$$M_{ij}=\frac{\partial H[\sigma ]}{\partial {\sigma }_i\partial {\sigma }_j}+\mu {\delta }_{ij}=H_{ij}^{\prime \prime }+\mu {\delta }_{ij}$$

$$\Sigma (E,\mu )=D(\mu )-\frac{E^2{f}^{\prime \prime }(1)+(E+\mu {)}^2{f}^{\prime }(1)}{2\{f(1)[f^{\prime \prime }(1)+f^{\prime }(1)]-f^{\prime }(1{)}^2\}}.$$

$${\partial }_t{\sigma }_i(t)=-\mu (t){\sigma }_i(t)-\frac{\partial H}{\partial {\sigma }_i}(t)+{\xi }_i(t)$$

$$⟨{\xi }_i(t){\xi }_j(t^{\prime })⟩=2T_f{\delta }_{ij}\delta (t-t^{\prime })$$

$$P[\sigma (0)]=\frac{e^{-\beta H(\sigma (0))}}{Z(\beta )}$$

$$\mu (t)=T_f-\frac{1}{N}\sum _i{\sigma }_i(t)\frac{\partial H}{\partial {\sigma }_i}(t).$$

$${\partial }_{t}C(t,t^{\prime })=-\mu (t)C(t,t^{\prime })+{\int }_0^{t}ds{f}^{\prime \prime }(C(t,s))R(t,s)C(s,t^{\prime })+{\int }_0^{t^{\prime }}ds{f}^{\prime }(C(t,s))R(t^{\prime },s)+\beta f^{\prime }(C(t,0))C(t^{\prime },0)$$

$${\partial }_{t}R(t,t^{\prime })=-\mu (t)R(t,t^{\prime })+{\int }_{t^{\prime }}^{t}ds{f}^{\prime \prime }(C(t,s))R(t,s)R(s,t^{\prime })$$

$$\mu (t)=T_f+{\int }_0^{t}ds{f}^{\prime \prime }(C(t,s))R(t,s)C(t,s)+{\int }_0^{t}ds{f}^{\prime }(C(t,s))R(t,s)+\beta f^{\prime }(C(t,0))C(t,0)$$

$$\text{with energy given by}E(t)=-{\int }_0^t{f}^{\prime }(C(t,s))R(t,s)ds-\beta f(C(t,0)).$$

$${\chi }_{\text{EA}}=\int d\lambda \frac{\rho (\lambda )}{\lambda }=\frac{\mu +\sqrt{{\mu }^2-4f^{\prime \prime }(1)}}{2f^{\prime \prime }(1)},$$

$$F(t)$$

$$\equiv E(t)-T_{f}S(t)$$

$$y_0{f}^{\prime }(q_0)q_0=\beta f^{\prime }(q_{12})q_{12}.$$

$$\displaystyle\frac{dC(t)}{dt}=-TC(t)+\beta\int_{0}^{t}ds\;f^{\prime}\big{(}C(t-s)\big{)}\frac{dC(s)}{ds}\,.\qquad$$

$${\beta }^2{f}^{\prime }(q_1)(1-q_1)=q_1.$$

$${\Delta }_{r}x(t)=\frac{x(t;\Delta t)-x(t;\Delta t=0.1)}{x(t;\Delta t)}.$$

$$E(\mu )=E_{\infty }+A({\mu }_{mg}-\mu {)}^{\alpha },$$

$$-2V_{\text{1RSB}}(q_{12},\chi ,q_0,y)=\chi f^{\prime }(1)+y(f(1)-f(q_0))+\frac{1}{y}\log \left(\frac{\chi +y(1-q_0)}{\chi }\right)+\frac{q_0-q_{12}^2}{\chi +y(1-q_0)}+2\beta f(q_{12}).$$

$$\{\begin{array}{ll}\text{hfill}{\partial }_{\chi }{V}_{\text{1RSB}}=0& ⟹\chi (1-q_{12}^2)+y(1-q_0{)}^2-\chi (\chi +y(1-q_0){)}^2{f}^{\prime }(1)=0\text{vspace}2mm\\ \text{hfill}{\partial }_{q_0}{V}_{\text{1RSB}}=0& ⟹q_0-q_{12}^2-(\chi +y(1-q_0){)}^2{f}^{\prime }(q_0)=0\text{vspace}2mm\\ {\partial }_{q_{12}}{V}_{\text{1RSB}}=0& ⟹q_{12}-\beta (\chi +y(1-q_0))f^{\prime }(q_{12})=0\text{vspace}2mm\\ \text{hfill}\text{marginality}& ⟹{\chi }^2{f}^{\prime \prime }(1)=1\end{array}$$

$$E$$

$$\mu$$

$$\chi$$

$$y$$

$$E=E_{th}\equiv -{\chi }_{mg}{f}^{\prime }(1)-y_{0}f(1),\mu ={\mu }_{mg}\equiv 2\sqrt{f^{\prime \prime }(1)}.$$

$$y(T_{\text{SF}})=\frac{f^{\prime \prime \prime }(1)}{2[f^{\prime \prime }(1){]}^{3/2}}$$

$$-2V_{\text{RS}}(q_{12},\chi )=\chi f^{\prime }(1)+\frac{1-q_{12}^2}{\chi }+2\beta f(q_{12}),$$

$$\{\begin{array}{l}1-q_{12}^2-{\chi }^2{f}^{\prime }(1)=0\\ q_{12}-\chi \beta f^{\prime }(q_{12})=0\end{array}$$

$$q_{12,\text{SF}}=\sqrt{1-\frac{f^{\prime }(1)}{f^{\prime \prime }(1)}}.$$

$$\mu =\chi f^{\prime \prime }(1)+(\chi +y)f^{\prime }(1)-y{f}^{\prime }(q_0)q_0+\beta f^{\prime }(q_{12})q_{12}$$

Full text

We revisit the long-time limit of the out of equilibrium dynamics of mean-field spherical mixed $p$-spin models. We consider quenches (gradient descent dynamics) starting from initial conditions thermalized at some temperature in the ergodic phase. We perform numerical integration of the dynamical mean-field equations of the model and we find an unexpected dynamical phase transition. Below an onset temperature $T_{\text{onset}}$, higher than the dynamical transition temperature $T_{\text{\tiny MCT}}$, the asymptotic energy goes below the “threshold energy” of the dominant marginal minima of the energy function and memory of the initial condition is kept. This behavior, not present in the pure spherical $p$-spin model, resembles closely the one observed in simulations of glass forming liquids. We then investigate the nature of the asymptotic dynamics, finding an aging solution that relaxes towards deep marginal minima, evolving on a restricted marginal manifold. Careful analysis however rules out simple aging solutions. We compute the constrained complexity in the aim of connecting the asymptotic solution to the energy landscape.

I Introduction

Understanding the relation between the dynamical behavior and the underlying energy landscape is a fundamental question in the physics of glassy systems. The same question arises in a broader interdisciplinary perspective in optimization, information theory and computer science. In these fields, the comprehension of the performances of classes of iterative algorithms optimizing complicated functions represents both a practical and a fundamental theoretical problem. One notable example, among many, is provided by supervised learning, where one formalizes data learning as the problem of minimizing a loss function that depends parametrically on the data Friedman et al. (2001).

In the years, concepts originated from glassy physics have been proven useful in the description of neural networks Amit and Amit (1992), protein folding Frauenfelder et al. (1991); Sfatos and Shakhnovich (1997), supervised learning Seung et al. (1992); Engel and Van den Broeck (2001), information and computational complexity theory Monasson et al. (1999); Mezard et al. (2009), and much more.

Fundamental unifying concepts have been provided by theory of Spin Glasses. Its interaction with computer science has produced and refined powerful optimization algorithms, from the ones that aim to mimic physical dynamics, such as Simulated Annealing Kirkpatrick et al. (1983), Parallel Tempering Earl and Deem (2005) and Quantum Annealing Finnila et al. (1994), to a family of message passing algorithms, such as Survey Propagation Mézard et al. (2002) and Approximate Message Passing Donoho et al. (2010), which are based on the the mean field equations for spin glasses. Finally, Spin Glass theory has provided a collection of models that, disregarding the finite-dimensional embedding of real materials, allow to study glassy phenomena in a simplified setting.

The model that more than any other has shaped our ideas about the relation between landscapes and dynamics has been the spherical $p$-spin model Crisanti and Sommers (1992), a simple system of spherical spins that interact through a disordered and long-range $p$-body Hamiltonian. The model exact solution offers a detailed description for the long time dynamics Crisanti et al. (1993); Cugliandolo and Kurchan (1993), the thermodynamics Crisanti and Sommers (1992) and the structure of metastable states Crisanti and Sommers (1995). The ensuing physical picture is at the basis of the Random First Order Transition (RFOT) theory Castellani and Cavagna (2005).

The spherical $p$-spin model has been used as starting point to understand complex dynamics in several interdisciplinary contexts, for example, it has been proposed as a simplified model to understand the properties of the loss function and the dynamics of learning in deep networks Choromanska et al. (2015). More generally in Machine Learning many fundamental problem involving tensors are NP-hard Hillar and Lim (2013), e.g. tensor factorization, tensor completion and tensor Principal Component Analysis Montanari and Richard (2014); Ben Arous et al. (2018). In this class of tensor problems the hardness in solving the problem has been directly connected to the properties of the underlying energy landscape Arous et al. (2017); Ros et al. (2019); Biroli et al. (2019). Thus better understanding relaxation in the spherical $p$-spin model has direct consequences on fundamental NP-hard problems in Machine Learning.

A widely used statistical characterization of the energy landscape sampled by the glassy dynamics is the numerical study of the inherent structures (ISs), defined as the local minima of the energy potential reached by a steepest descent procedure Sastry et al. (1998); Sciortino et al. (1999). In simulations of glass forming liquids, it is observed that starting from thermalized configurations at temperature $T$, the energy of the corresponding ISs concentrates around a value $E_{\text{\tiny IS}}(T)$ that depends on $T$ only below a characteristic onset temperature $T_{\text{onset}}$. Remarkably, this temperature is higher than the estimated ‘mode coupling temperature’ $T_{\text{\tiny MCT}}$ where Mode Coupling Theory would predict structural arrest. Not limited to liquids, a similar effect is also observed in 3D Heisenberg spin glasses Baity-Jesi and Martin-Mayor (2019). From the algorithmic point of view, this dependence shows the effectiveness of Simulated Annealing: that is, in order to minimize the energy, equilibrating at lower and lower temperatures is a better strategy than a sudden quench from a high energy state. A similar effect is known in optimization, where a ‘smart initialization’ from carefully chosen initial conditions can greatly improve the quality of the solutions. A strategy of this kind, e.g. a spectral initialization, is often adopted to allow simple algorithms to solve non-convex problems Lu and Li (2017). However a general theory based on the energy landscape explaining the effect of a ‘smart initialization’ is lacking. We provide in this work a first comprehensive study of the effect of starting from better than random initial configurations.

According to the common belief, the temperature dependence of the IS energy is a consequence of activation and should not appear in models with long range interactions: initial conditions thermalized above $T_{\text{\tiny MCT}}$ would all have the same value of the IS energy, any memory of the initial condition would be lost after a long enough time, and Simulated Annealing would not give any benefit respect to a sudden quench. This would also mean that in this class of NP-hard problems any smart initialization would be useless. However this belief is based on the solution of the spherical pure $p$-spin model Crisanti and Sommers (1992), which, in a nutshell, gives the following dynamical picture: (1) From the point of view of equilibrium dynamics, the model provides an exact realization of Mode Coupling Theory scenario. One finds dynamical freezing at a temperature $T_{\text{\tiny MCT}}$ where the free-energy does not present singularities Crisanti et al. (1993). (2) Randomly chosen initial conditions evolved according to the Langevin dynamics at temperature $T_{f}<T_{\text{\tiny MCT}}$ fail to equilibrate. They fall in an aging state which fails to attain a time translation invariant long time regime and present non trivial response to perturbations Cugliandolo and Kurchan (1993). (3) Any memory of the initial condition is lost in the aging state, the choice of an initial condition thermalized at temperatures $T_{in}>T_{\text{\tiny MCT}}$ would not change the asymptotic of a dynamics in presence of a bath at temperature $T_{f}<T_{\text{\tiny MCT}}$. There is an universal threshold energy $E_{th}$ that attracts the zero temperature dynamics of any such initial conditions, i.e. $E_{IS}(T)=E_{th}$ for all $T>T_{\text{\tiny MCT}}$. (3bis) Initial conditions thermalized at $T<T_{\text{\tiny MCT}}$ live in metastable states that can be followed down in temperature Franz and Parisi (1995); Barrat et al. (1996, 1997). In that case one reaches an energy $E_{0}(T)$ which is monotonous in $T$ and smaller than then threshold value $E_{th}$. Unfortunately, such initial conditions cannot be generated in a finite time.

Such a detailed dynamical picture is deeply rooted in the structure and organization of the stationary points of the energy and the free-energy landscapes Cavagna et al. (1998). (A) For $E<E_{th}$ almost all stationary points of the energy are isolated and stable minima. The stability-gap of the minima (i.e. the minimum eigenvalue of the Hessian) is uniquely a function of the energy. Saddle points are exponentially suppressed and have only a finite number of unstable directions. (B) For $E>E_{th}$ almost all stationary points are saddles. (C) Minima with $E=E_{th}$ are marginal. Dynamics is attracted by these threshold states, no matter the effort of thermalization at $T>T_{\text{\tiny MCT}}$. It is clear that Simulated Annealing is of no use for such a system. (D) Below-threshold energy minima are separated by extensive barriers. There is a one-to-one correspondence between energy minima and finite temperature metastable states (TAP free-energy minima). Metastable states can be followed up and down in temperature.

The above general picture, both for the dynamics and for the structure of the landscape, is not restricted to the spherical $p$-spin model, it is generic to mean-field models with a discontinuous glass transition (one step replica symmetry breaking or 1RSB in spin glass jargon) and has been confirmed by the exact description of glasses of particles in the limit of infinite dimension Charbonneau et al. (2014). However there are important aspects of the physics of the spherical $p$-spin that are deeply model-specific. Its Hamiltonian is homogeneous and does not allow for bifurcations or merging of metastable states while changing the temperature. The energy landscape is simpler than more generic mean-field models. It is well known that in other mean-field modes like the Ising $p$-spin model, the Potts glass or even spherical models with mixtures of multi body interactions Barrat et al. (1997) when metastable states are followed down in temperature bifurcations are encountered. This leads to the glass-to-glass Gardner phase transition that has received a lot of attention in the last years Berthier et al. (2019).

The study of metastable states below such a transition reveals puzzling aspects. The mean-field solutions describing the evolution of states in temperature may unexpectedly disappear Sun et al. (2012) as the temperature is lowered. This anomalous behavior has not received a satisfactory explanation, and its implications remain little studied. In this paper we re-examine dynamics in models where these phenomena occur and challenge crucial aspects of the $p$-spin dynamical picture: the emergence of a universal dynamical threshold and the loss of memory of high temperature initial conditions.

The paper is organized as follows. In Sec. II we define of the model. In Sec. III we discuss the properties of the landscape and the number of stationary points of the energy surface of the model. Sec. IV presents the equations describing the off-equilibrium dynamics and the effect of a non-random initial condition. Sec. V constitutes the core of the manuscript and discusses the energy of inherent structures, the puzzling properties of the zero temperature dynamics and its relation with the landscape. The results are extended to relaxation in presence of a small thermal noise in Sec. VI. We discuss an approximated dynamical solution and the emergence of an onset temperature in Sec. VII. We comment on the results and their consequences in Sec. VIII. The paper comprises a large number of appendixes devoted to the discussion of more technical aspects of the work that support and complement the main text. In App. A we discuss finite temperature dynamics in equilibrium. In App. B we describe the numerical solution of the dynamical equations. App. C deals with the properties of the response function at short times in the aging regimes. In App. D we present some formal solutions of the long time asymptotic aging regime and the connections with the effective potential. An approximate solution to the asymptotic dynamics is analysed in App. E. Finally we show the computation of the constrained complexity of energy minima in App. F.

II Model definition

We concentrate on the mixed $p$-spin spherical model. This model, while presenting a more generic phenomenology, keeps much of the analytical simplicity of the pure spherical $p$-spin. In thi model states do bifurcate, and the above mentioned anomalous behavior occurs. The Hamiltonian of the mixed $p$-spin model is just a sum of $p$-spin interacting terms with different $p$ values

[TABLE]

Choosing the couplings $J^{(p)}_{i_{1}...i_{p}}$ as independent Gaussian variables of zero mean and variance $\mathbb{E}(J_{p}^{i_{1}...i_{p}})^{2}=\frac{p!}{2N^{p-1}}$, and the coefficients $a_{p}\geq 0$ with $\sum_{p}a_{p}<\infty$, we have that $H_{J}$ is a random Gaussian function on the $N$-dimensional sphere ($\sum_{i}\sigma_{i}^{2}=N$) with covariance

[TABLE]

We choose $f(q)$ as a polynomial in a way that the model has a random first order (RFOT) transition in thermodynamics, associated to an ideal mode coupling transition in equilibrium dynamics. The pure $p$-spin model corresponds to the case of a single monomial $f(q)=\frac{1}{2}q^{p}$; we study a more generic class of polynomials $f$ where at least two coefficients are positive (mixed $p$-spin model). While we keep the function $f$ generic in formulas, all our numerical results are presented for $f(q)=\frac{1}{2}(q^{3}+q^{4})$, and compared with the pure case with $p=3$. We checked in other examples, that all the results we get in this particular case are typical of general mixed models.

To investigate how much of the above picture based on the solution of the pure $p$-spin model is generic, we perform a high precision integration of the dynamical equations describing the evolution of the system in the thermodynamic limit, starting from a thermalized initial condition at temperature $T$ and subsequently undergoing a gradient descent dynamics. Our main finding is that the above picture is not generic and inhomogeneous models show several unexpected features. The value of the asymptotic energy predicted under the assumption of aging with loss of memory is only found starting from high enough temperatures. We find strong evidence for the existence of an onset temperature $T_{\text{onset}}>T_{\text{\tiny MCT}}$ below which the energy of the IS does depend on the initial temperature $T$ and goes below the threshold value.

III The energy landscape

In order to put in the right frame our dynamical results, we would like first to discuss the structure and organization of the stationary point of the mixed Hamiltonian. The stationary points of the Hamiltonian $H[\sigma]$ on the sphere $\sum_{i}\sigma_{i}^{2}=N$ verify

[TABLE]

The parameter $\mu$, that we call radial reaction force, is a Lagrange multiplier that insures the spherical constraint. In any stationary point, it takes the value

[TABLE]

The radial reaction is directly related to the nature and the stability of the stationary points, in fact, the Hessian matrix (intended to be restricted to fluctuations on the sphere) reads

[TABLE]

It is well known Cavagna et al. (1998) and rigorously proven Auffinger et al. (2013) that $H^{\prime\prime}_{ij}$ is a GOE random matrix with variance $\text{Var}[{H^{\prime\prime}_{ij}}]=\frac{1}{N}f^{\prime\prime}(1)$. Thus the matrix $M$ has a semicircular spectral distribution $\rho(\lambda)=\sqrt{(\lambda-\mu)^{2}-4f^{\prime\prime}(1)}/(\pi\sqrt{f^{\prime\prime}(1)})$ which has a lower band edge at $\lambda_{min}=\mu-2\sqrt{f^{\prime\prime}(1)}$. We notice that the value of the radial reaction determines the nature of the stationary points. 111As well know in random matrix theory, while macroscopic fluctuations of the spectral bulk are suppressed by factors scaling exponentially in $N^{2}$, the probability of having a finite number of eigenvalues at finite distance from the spectral support is exponentially small in $N$. As in the case of the pure model, given value of $\mu$ such that the spectral support is positive, besides the dominating minima we have saddle points with a finite number of negative directions. It is natural to define a marginal value of the radial reaction $\mu_{mg}=2\sqrt{f^{\prime\prime}(1)}$ that separates minima from saddles with an extensive number of negative directions. Stationary points of $H$ are overwhelmingly minima if $\mu>\mu_{mg}$ and they are saddles or maxima if $\mu<\mu_{mg}$. Notice that all marginal minima with vanishing spectral gap lie on the manifold specified by $\mu(\sigma)=\mu_{mg}$. In the pure model, in any stationary point the radial reaction and the energy verify $\mu=-pE$. In the mixed model, there is not such a one-to-one relation between $\mu$ and $E$ in stationary points. This can be seen computing the number of stationary points in (2), for which both $E$ and $\mu$ are fixed ${\cal N}(E,\mu)=e^{N\Sigma(E,\mu)}$.

The computation of the complexity $\Sigma(E,\mu)$ was performed in Ref. Arous et al. (2018) and we just quote the result here (for completeness, a physicist’s derivation and the explicit formula for $D(\mu)$ is provided in Appendix F)

[TABLE]

The above expression holds when returns a positive value.

Fixing $\mu$, the complexity $\Sigma(E,\mu)$ is a parabola in $E$ with a curvature that is finite for mixed models (see Fig. 1 where we show the complexity for the (3+4)-spin model) and that tends to infinity in the pure limit. For $\mu<\mu_{max}$, with $\mu_{max}$ implicitly defined by equations $0=\Sigma(E_{max},\mu_{max})=\partial_{E}\Sigma(E_{max},\mu_{max})$, there is a finite interval of energies $[E_{-}(\mu),E_{+}(\mu)]$ such that within the interval $\Sigma(E,\mu)>0$. Most importantly, this is true for the threshold value of the radial reaction $\mu_{mg}$. Seen from the point of view of $\mu$ the organization of the stationary points resembles the one of the pure model, however from the point of view of the energy there are deep differences. In the pure model, for a given level of energy, the index of any stationary point 222The index counts the number of negative eigenvalues of the Hessian at a stationary point. differs at most by a finite number from the the index of dominating ones. In the mixed models, for a given energy, stationary points coexist with macroscopically different indices: i.e. their spectra are shifted by a finite amount. A unique relation between $E$ and $\mu$ only holds for the exponentially dominating stationary points.

IV Dynamical equations in the out-of-equilibrium regime

Given the Hamiltonian in Eq. (1), we consider the following Langevin dynamics

[TABLE]

where $T_{f}$ is the temperature of the thermal bath (later set to zero) and $\beta=1/T$ is the inverse temperature at which the initial condition is equilibrated. $\mu(t)$, the time dependent radial reaction force, is the Lagrange multiplier that constrains the dynamics on the sphere, and at time $t$ it takes the value

[TABLE]

Taking the thermodynamic limit, $N\rightarrow\infty$, at finite times, Eq. (6) imply closed integro-differential equations Sompolinsky and Zippelius (1982); Crisanti et al. (1993); Mezard et al. (1987) for the correlation $C(t,t^{\prime})\equiv\langle\sigma_{i}(t)\sigma_{i}(t^{\prime})\rangle$ and for the response $R(t,t^{\prime})\equiv\frac{\partial\langle\sigma_{i}(t)\rangle}{\partial\xi_{i}(t^{\prime})}$. For an initial temperatures $T$ greater than the Kauzmann temperature of thermodynamic phase transition $T_{\text{\tiny K}}$ the dynamical equations read Barrat et al. (1996)

[TABLE]

The initial condition enters the above equations in a rather simple way. For $\beta=\beta_{f}$ one gets equilibrium dynamics, whose MCT transition is illustrated in Appendix A. For $\beta=0$ the equations reduce to the usual form valid for uncorrelated initial condition. In this case, for $T_{f}<T_{\text{\tiny MCT}}$ the equations admit a Cugliandolo-Kurchan weak ergodicity breaking aging solution completely analogous to the one of the pure model. It is possible to show that in this solution, for $T_{f}\to 0$ the radial reaction $\mu(t)$ tends to the marginal value $\mu_{mg}$ for $t\to\infty$, while the energy tends to the threshold value $E_{th}=\frac{-f^{\prime}(1)^{2}+f(1)[f^{\prime}(1)+f^{\prime\prime}(1)]}{f^{\prime}(1)\sqrt{f^{\prime\prime}(1)}}$. In correspondence of this value of the energy, the complexity $\Sigma(E_{th},\mu)$, as a function on $\mu$, is maximum for $\mu=\mu_{mg}$. In other words, $E_{th}$ is not the energy of the most numerous marginal states, but the energy where marginal states exponentially dominate the complexity.

It is possible to see that for $T<\infty$ the same aging solution solves the dynamical equations in the long time limit, if one assumes that memory of the initial condition is lost and $C(t,0)\to 0$ at a large time. If this would be the correct description for all $T>T_{\text{\tiny MCT}}$, then the asymptotic energy would be independent of $T$ in this domain.

We study the dynamical equations in (8) by numerical integration. The simplest possibility, which we adopt here, is a fixed time step $\Delta t$ first-order Euler discretization algorithm, that gives very reliable results at short times Franz and Mezard (1994) and allows for reliable extrapolations in the $\Delta t\to 0$ limit (see Appendix B). Other algorithms with variable time step have been proposed and used in the literature Kim and Latz (2001); Berthier et al. (2007). Unfortunately in the case of mixed $p$-spin model with correlated initial condition ($\beta>0$), these algorithms appear to be unstable at short times and do not allow any improvement over the simplest one. With the Euler integration scheme and a maximum step of $\Delta t=0.1$, we reach times of order $10^{3}$. These time scales are often not large enough to allow for a simple naive extrapolation to the infinite time limit. Nonetheless, under some conditions, will allow us to make clear claims on the large time dynamics. We concentrate in this paper on the case of zero temperature dynamics, where the energy monotonically decreases and once below threshold it cannot rise up again. We show however in Section VI that this $T_{f}=0$ result extends to a range of positive temperatures, $T_{f}>0$.

V Affinities and divergences between mixed and pure models

We discuss here the core question of this paper, and compare the energy reached by the dynamics from different initial conditions, to the threshold energy that can be computed from the asymptotic solution with $\beta=0$ Cugliandolo and Kurchan (1993) or by computation of the complexity of minima Crisanti and Sommers (1995). For the 3+4 model the value of the threshold energy is $E_{th}=-71/42\simeq-1.6905$. In Fig. 2 we show the curves of the energy as a function of time for different $T$, ranging from $T=\infty$ to $T=0.796$, the dynamical transition temperature of the model being $T_{\text{\tiny MCT}}=0.805166$ (see Appendix A). Data are obtained from an integration with step $\Delta t=0.1$, reaching the maximum time of $2500$. This step size gives a relative integration error smaller than $10^{-5}$ (see detailed discussion in the Appendix B).

On the time scale we can reach and for $T\gtrsim 1$ the energy seems to have reached an asymptotic behavior well fitted by a power law $E(t)=E_{th}+a/t^{\gamma}$ with $a$ slightly dependent on temperature and $\gamma=0.66\pm 0.01$, making us confident that $E_{th}$ is the asymptotic value of the energy in this range of temperatures. On the other hand, if $T$ is small enough but still larger than $T_{\text{\tiny MCT}}$, the curves of the energy go below $E_{th}$. For comparison in the upper inset we present the same curves for the pure model with $p=3$. In that case manifestly the energy tends to the threshold value for all $T\geq T_{\text{\tiny MCT}}$. These results suggest a scenario with a new dynamical transition in the mixed model, with an onset temperature $T_{\text{onset}}$ separating a memoryless phase for $T>T_{\text{onset}}$ from a memorious phase at $T<T_{\text{onset}}$. It is also interesting to look at $T<T_{\text{\tiny MCT}}$. In our data we do not observe anything special happening at $T_{\text{\tiny MCT}}$, it is only below a state following temperature $T_{\text{\tiny SF}}\approx 0.798$ that we observe relaxation becoming fast and the energy decaying exponentially to its asymptote. We will see that we can predict the asymptotic energy in this region through a quasi-static solution of the dynamical equations (‘state following’ with memory of the initial condition).

We will see in the next section that a small temperature $T_{f}$ in the dynamical equations does not change this behavior.

In order to understand the dynamical mechanisms that allow to beat the threshold, we look at the behavior of response and correlations. A first quantity we would like to understand is the correlation with the initial state, $C(t,0)$. According to the usual weak ergodicity breaking scenario this quantity should vanish at large $t$. Both in the mixed and in the pure model we observe that the relaxation of $C(t,0)$ is slower and slower as $T$ decreases. As for the energy we can identify three temperature regimes for $T>T_{\text{\tiny K}}$: (I) a high temperature regime $T>T_{\text{onset}}$ where memory of the initial condition is lost and $C(t,0)\to 0$; (II) an intermediate ‘hic sunt leones’ regime $T_{\text{\tiny SF}}<T<T_{\text{onset}}$ where the relaxation is slow and any extrapolation is difficult; (III) a low temperature regime $T<T_{\text{\tiny SF}}$ where $C(t,0)$ relaxes exponentially fast to a non zero value $q_{12}$.

While no aging is found in regime III, in regimes I and II we find aging, qualitatively similar in the two cases. In Fig. 3 we show $C(t,s)$ as a function of $t-s$ for several times $t$ and initial temperature $T=1.01T_{\text{\tiny MCT}}$. The curves show a clear aging behavior with the dynamics decorrelating slower and slower as time passes, followed by a fast drop of the correlation for small values of $s$ to the value $C(t,0)$.

Despite the energy goes below threshold, we find compelling evidence that the dynamics remain critical, and approaches asymptotically marginally stable minima. We reach this conclusion, studying the asymptotic behavior of the radial reaction $\mu(t)$, that for $t\to\infty$ determines the spectral gap of the asymptotically visited states. In Fig. 4 we can clearly see that both in region I and II, the radial reaction tends to the marginal value $\mu_{mg}\equiv 2\sqrt{f^{\prime\prime}(1)}$ of marginal minima. In the region III by contrast, the system reaches stable minima and $\mu(\infty)>\mu_{mg}$.

Aging behavior is often qualified studying the relation between response and correlation Cugliandolo et al. (1997). Fig. 5 shows the parametric plot of the integrated response $\chi(t,s)\equiv\int_{s}^{t}R(t,u)\;du$ versus $C(t,s)$ as a function of $s$ for various values of $t$. In zero temperature dynamics, in the large time limit, the integrated response has a jump in $C=1$ to the intrastate susceptibility $\chi_{\text{\tiny EA}}$. In a given minimum

[TABLE]

and in marginal minima $\mu=\mu_{mg}$ implies $\chi_{\text{\tiny EA}}=\chi_{mg}\equiv 1/\sqrt{f^{\prime\prime}(1)}$. The curves in Fig. 5 show the formation of a jump in $C=1$ very well compatible with this value.

The memoryless solution Cugliandolo and Kurchan (1993) predicts a linear behavior for $C<1$ with slope given by the ‘fluctuation-dissipation ratio’ $y_{0}\equiv\sqrt{f^{\prime\prime}(1)}/f^{\prime}(1)-\chi_{mg}$.

On the time scales we can access, both the pure and the mixed models show, for all the values of $T>T_{\text{\tiny SF}}$, after the jump to $\chi\simeq\chi_{mg}$, an approximately linear part. The linear behaviour continues till a rather well defined time dependent value of the correlation $q_{0}(t)>C(t,0)$ where the response essentially stops to increase. The value $q_{0}(t)$ identified in Fig. 3 marks the value that separates slow dynamics, where an effective temperature emerges, from fast partial decorrelation from the initial condition, during which the system does not respond (see Appendix C).

Having established that the system relaxes towards marginal minima, it comes the question of which marginal minima are selected. In an attempt to address this question, we study the constrained complexity Capone et al. (2006) $\Sigma(E,\mu,q_{12};T)$ of energy minima of fixed radial reaction $\mu$, energy $E$ and correlation $q_{12}$ from a reference configuration thermalized at temperature $T$. This quantity for $q_{12}=0$ reduces to the one in Eq. (5) studied in section III. The detailed computation of this constrained complexity is presented in appendix F. The computation reveals that for $q_{12}>0$ we have qualitatively the same picture illustrated in Fig. 1, but with $\Sigma(E,\mu,q_{12};T)<\Sigma(E,\mu)$, being $\Sigma(E,\mu,q_{12};T)$ monotonically decreasing in $q_{12}$. One can define a generalized threshold energy that depends on $T$ and $q_{12}$, where marginal states dominate the constrained complexity. It turns out that $E_{th}(q_{12},T)$ is a decreasing function of $q_{12}$ and $T$. Unfortunately, as explained in detail in appendix F, while the resulting $E_{th}(q_{12},T)$ has the right qualitative behavior, the corresponding minima have features that are incompatible with the actual attractors of the dynamics. The complete static characterization of the attractors of the dynamics remains an important open problem.

VI Extending results to finite temperature relaxation dynamics

As we said in section IV, for $T_{f}=0$ the energy monotonically decreases, while if the bath temperature $T_{f}$ is positive this is not necessarily the case. If we integrate the dynamical equations for a small temperature $T_{f}$, we observe that the finite time energy is continuous, $E(t,T_{f})=E(t,T_{f}=0)+o(T_{f})$. We would like to show that, for small enough $T_{f}$, this implies that $\lim_{t\to\infty}E(t,T_{f})<E_{th}(T_{f})$.

At $T_{f}>0$ the Lyapunov function of the Langevin dynamics, that decreases monotonically in time, is the free-energy functional Kubo et al. (2012):

[TABLE]

where $P_{t}[\sigma]$ is the probability distribution of configurations at time $t$ induced by the initial condition and the dynamics.

In principle the energy function $E(t)$ could be non-monotonous in time for $T_{f}>0$, because an entropy increase $\Delta S(t)=S_{\infty}-S(t)>0$ can induce an energy increase $\Delta E>0$. However, being the free-energy a decreasing function of time, then we have $\Delta F=\Delta E-T_{f}\Delta S<0$, and the energy increase is upperbounded by $\Delta E<T_{f}\Delta S$. This is a quantity that tends to zero for $T_{f}\to 0$. Consequently, if the temperature is small enough, the energy $E(t)$, which is below threshold by a $O(1)$ amount, cannot be pushed back to the threshold value.

In Fig. 6 we show the energy reached at several times $t\in\{125,250,500,1000\}$ during the integration of the dynamical equations as a function of the temperature $T_{f}$. The topmost full line is the analytic prediction for the threshold energy as a function of $T_{f}$. We observe that the dynamical energy lies below the threshold energy and shows a continuous and nicely linear behavior in $T_{f}$. The weak time dependence of the data shown in Fig. 6 suggests the energy is close to its asymptotic value. Even willing to insist that an eventual entropy increase may induce an energy increase at later times, this increase would be tending to zero for small $T_{f}$ and thus, by continuity, there would be a finite $T_{f}$ range where the energy will certainly remain below the threshold value even at infinite time.

VII An approximate asymptotic solution to the dynamics

The dependence of the final energy on $T$ in regime II testifies a form of memory of the initial condition. Differently thermalized configurations lie in basins of attraction of different marginal states. The asymptotic ansatz of Cugliandolo and Kurchan has been generalized in Barrat et al. (1997) supposing ‘aging within a metabasin’ and a non vanishing correlation with the initial state $\lim_{t\to\infty}C(t,0)=q_{12}>0$. Unfortunately, the resulting equations for $\chi,y,q_{12},q_{0}$ only have an aging solution different from the ‘amnesic’ one ($q_{12}=q_{0}=0$) in the interval $T_{\text{\tiny SF}}\leq T\leq T_{0}=0.803<T_{\text{\tiny MCT}}$. Even when the solution exists, the values of the various parameters do not match the observations, in particular one finds $y\approx 0.3$ while in dynamics $y\approx 0.52$. This aging solution terminates at the temperature $T_{\text{\tiny SF}}$ where $q_{0}$ tends to one. Below that temperature, there is a stable ‘state following’ solution with no aging, $q_{12}>0$, $\mu>\mu_{mg}$ and $\chi<\chi_{mg}$ (see Appendix D for a detailed derivation). This solution correctly reproduce the asymptotic values of the energy, correlation, radial reaction of the large time dynamics for temperatures $T\lesssim T_{\text{\tiny SF}}$. Despite all efforts, we have been unable to find a different asymptotic dynamical ansatz describing aging in the ‘hic sunt leones’ regime II and/or to find in the numerics an indication on where the memory of the initial condition could affect the asymptotic solution. The lack of simple aging solutions is certainly related, and equally paradoxical, to the phenomenon of lost of solution in the ‘state following’ procedure of Sun et al. (2012).

Insisting in using a 1RSB dynamical ansatz with $\chi_{\text{\tiny EA}}=\chi_{mg}$, $\mu=\mu_{mg}$ and $y\simeq y_{0}$, which is supported by the numerical solution, one can easily derive the following relation

[TABLE]

A second relation can be obtained from the observation that $\chi_{\text{tot}}(t)\equiv\int_{0}^{t}ds\;R(t,s)\approx\chi_{mg}+y_{0}(1-q_{0}(t))$ depends approximately linearly on $C(t,0)$ even at finite times (see dotted lines in Fig. 5). Moreover this relation does not present an appreciable temperature dependence in the range $T_{\text{\tiny SF}}\leq T\lesssim 1$ (data shown in Appendix E) and implies a linear relation between $q_{12}$ and $q_{0}$ that can be estimated from two analytically known limits: for $T$ large enough $q_{12}\!=\!q_{0}\!=\!0$, while at $T\!=\!T_{\text{\tiny SF}}$ we have $q_{0}\!=\!1$ and $q_{12,\text{\tiny SF}}=\sqrt{1-f^{\prime}(1)/f^{\prime\prime}(1)}$. Assuming the relation $q_{12}=q_{12,\text{\tiny SF}}\,q_{0}$, Eq. (10) admits a solution with $q_{0}>0$ if $T<T_{\text{onset}}\equiv{q_{12,\text{\tiny SF}}^{k}}/{y_{0}}$ where $k$ is defined by $f(q)\propto q^{k}$ for $q\to 0$ (in the 3+4-model $T_{\text{onset}}=0.91$). This approximated solution gives an asymptotic energy close to the one extrapolated from the numerical integration (see Fig. 7). Despite this is not an exact solution of the asymptotic equations, it gives us a strong indication that the passage from memoryless to memorious aging could be marked by a phase transition. In addition, it shows that the exact solution cannot be too far from a 1RSB weak ergodicity breaking solution inside a basin.

VIII Discussion and perspectives

Solving the out of equilibrium relaxation dynamics in the spherical mixed $p$-spin model, starting from a thermalized configuration at temperature $T$, we have uncovered several interesting and unexpected features. While the known solution of the relaxation dynamics in the pure $p$-spin model suggested the existence of a unique dynamical phase transition at $T_{\text{\tiny MCT}}$, where ergodicity breaks down, and of a unique threshold energy where the dynamics relaxes below $T_{\text{\tiny MCT}}$, our results about the mixed $p$-spin model reveal a much richer scenario.

Our main result is that there is a temperature range above $T_{\text{\tiny MCT}}$ where the final energy depends on $T$ and memory of the initial condition is kept. A temperature $T_{\text{onset}}$ marks the onset of this effect. This resembles realistic glass-formers, but was unexpected in mean field models, because such a phenomenon was ascribed to activated processes, which are absent in the out of equilibrium relaxation of mean-field models. Our results show, contrary to common wisdom, that at $T_{\text{onset}}$ activated processes do not play any crucial role, and the dependence on the initial configuration can arise even in purely relaxation dynamics.

From the numerical solution of the dynamical mean-field equations, this dynamical onset of memory may be a cross-over or a new genuine phase transition. We have found an approximate solution of the relaxation dynamics at large times that describes the onset as a phase transition and provides, an analytic prediction for the onset temperature. This provides finally a solvable (although approximate) model whose predictions can be compared with numerical experiments.

Another piece of common wisdom that the present work demystifies is the mantra that ‘relaxation dynamics in mean-field models relaxes to the threshold energy where the most numerous metastable states are found’. Calling threshold energy the energy where the relaxation dynamics converge if started from a random configuration ($T=\infty$), we have shown that (i) metastable states of infinite lifetime (i.e. locally stable states) do exist also above the threshold energy and (ii) the most numerous are above the threshold energy. The threshold energy can be thus understood as the energy value where the measure is dominated by marginal states. And this can be computed via the unconstrained complexity. Above the threshold energy, an exponential number of marginal states already exist, but the measure over stationary points is dominated by saddles. The same argument trivializes in the pure $p$-spin model since in the latter marginal states do exist only at the threshold energy, and so assuming convergence to marginal states determines uniquely the corresponding threshold energy.

We have show evidences that the correct recipe for understanding the large time limit of the relaxation dynamics is to assume that such a dynamical out of equilibrium process always tends to marginal states. For $T<T_{\text{onset}}$ the persistent memory of the initial configuration suggests to compute a constrained complexity of these marginal states that might be correlated with the initial condition. We have performed this task with different levels of accuracy and success.

For an initial temperature below $T_{\text{\tiny SF}}$, which is strictly smaller than $T_{\text{\tiny MCT}}$ in the mixed $p$-spin model, the starting configuration lies in a well defined metastable state that can be followed down to zero temperature. In this situation the dynamics never enters an aging regime and its asymptotic behavior can be computed from a standard state following computation.

The most interesting regime is the one reached with an initial temperature $T_{\text{\tiny SF}}<T<T_{\text{\tiny MCT}}$. The relaxation process falls asymptotically in an aging regime, but remains correlated with the initial configuration. In this regime, the constrained complexity, while confirming the existence of marginal states below threshold and correlated with the initial condition, fails to identify the ones chosen by the dynamics.

All our results remain true for dynamics performed at a small temperature $T_{f}$.

We believe that, contrary to the pure $p$-spin model, the mixed $p$-spin model we studied is generic enough that the properties we uncovered should typically hold in models with a mean-field RFOT, such as Potts glasses Gross et al. (1985) Ising p-spins Gardner (1985) and, importantly, glasses of spherical particles in the infinite dimensional limit Charbonneau et al. (2014).

In conclusion, on the basis on the above results, we have to severely revisit our common beliefs about the relaxation dynamics in $p$-spin models and more generically in models presenting a mean-field RFOT. The physical picture which has been spread and given credit in our whole scientific community, where a simple connection between the relaxation dynamics and the complexity of metastable states exists, is unfortunately valid only in the case of the pure $p$-spin model and is false in more general models.

In pure $p$-spin models strong symmetries lead to the equalities $T_{\text{\tiny SF}}=T_{\text{onset}}=T_{\text{\tiny MCT}}$ and these make disappear many of the interesting phenomena that we described in this work and have misled us for two decades with a too simplistic connection between dynamics and energy landscape. Mixed models appear to be well suited to describe the dynamical onset of glassy phenomenology, they have much more complex energy landscapes which deserve to be better studied in the future. The exact solution to the asymptotic dynamics remains unknown and the connection between asymptotic dynamics and energy landscape is still very open in spherical mixed $p$-spin models. Research is needed in two directions: one should find algorithms capable of reaching large times in solving the dynamical equations and new theoretical ideas on possible structures of the memorious aging are necessary to understand the asymptotic aging regime.

The richer scenario we have uncovered in this study is likely to have impact also on the many applications in machine learning we discussed in the introduction. Given away the simplest connection between relaxation dynamics and energy landscape, it is clear that more in depth studies of the energy landscape will be needed in order to predict the performance of algorithms (a first example can be found in Ref. Mannelli et al. (2019)), especially of those that do not start from random configuration, but use some smart initialization to improve their performances.

Acknowledgements.

We thank P. Charbonneau and G. Parisi for illuminating discussions. The research has been supported by the Simons Foundation (grants No. 454941, S. Franz; No. 454949, G. Parisi), by the European Research Council under the European Unions Horizon 2020 research and innovation programme (grant No. 694925, G. Parisi) and by the project “Meccanica statistica e complessità”, a research grant funded by PRIN 2015 (Agreement no. 2015K7KK8L). SF is a member of the Institut Universitaire de France.

Appendix A Equilibrium dynamics

The onset temperature $T_{\text{onset}}$ marks the point where landscape influences gradient descent dynamics. In realistic model glasses it has been found that the same temperature marks the onset of non-exponential relaxation, thus making a direct relation between equilibrium dynamics and energy landscape Sastry et al. (1998). It is natural therefore to search for a signature of $T_{\text{onset}}$ in equilibrium dynamics.

The aim of this Appendix is to briefly review the equilibrium dynamics at finite temperature and show that, while the usual MCT-like transition takes place at $T_{\text{\tiny MCT}}$, no anomalies around $T_{\text{onset}}$ are found. The equilibrium dynamical equation for the correlation function reads

[TABLE]

This is readily obtained imposing $T_{f}=T$ in Eq. (8), together with time translation invariance in the correlation function $C(t,s)=C(t-s)$ and the fluctuation-dissipation relation $R(t)=-\beta\frac{dC(t)}{dt}$. At high temperature, in the ergodic phase, $C(t)$ vanishes at large times. The dynamical MCT-like ergodicity breaking transition is signaled by the appearance of a non zero limit $q_{1}=\lim_{t\to\infty}C(t)$, with $q_{1}$ discontinuous at the transition. Evaluating (11) at large time, one finds that $q_{1}$ is solution of

[TABLE]

This equation coincides with the equation defining the non-trivial minimum of the equal-temperature effective potential discussed in detail in Appendix D.

The transition temperature, the highest temperature for which Eq. (12) has a non zero solution, satisfies $\beta^{2}f^{\prime\prime}(q_{1})(1-q_{1})^{2}=1$. Standard mode coupling analysis of the dynamics Götze (2008) describes the approach to the transition in terms of the formation of a plateau in the correlation function as a function of log-time. This analysis in particular shows that the approach to $q_{1}$ at $T_{\text{\tiny MCT}}$ follows a power law with an exponent $a$ which is determined by $\Gamma(1-a)^{2}/\Gamma(1-2a)=\frac{T}{2}f^{\prime\prime\prime}[q]/f^{\prime\prime}[q]^{3/2}$. For the mixed model $f(q)=\frac{1}{2}(q^{3}+q^{4})$, the transition temperature is T_{\text{\tiny MCT}}=21/\big{(}2\sqrt{37\sqrt{37}-55}\big{)}\approx 0.805166, corresponding to q_{1}=\big{(}1+\sqrt{37}\big{)}/12\approx 0.59023 and $a=-0.38504$. In Fig. 8 we display data from the numerical integration of (11) with a simple first order Euler algorithm similar to the one we use in the off-equilibrium case. The curves show that theoretical expectations are met: a transition occurs at $T_{\text{\tiny MCT}}$ according to the expected pattern and below $T_{\text{\tiny MCT}}$ the system is in a non-ergodic state.

We then investigate the temperature range around $T_{\text{onset}}$ to see if this temperature corresponds to a qualitative change in the equilibrium dynamics. Following Sastry et al. (1998), we fit the correlation function at intermediate times with a stretched exponential f(t)=a\exp\big{(}-(t/\tau)^{\hat{\beta}}\big{)} and study the behavior of the stretching exponent $\hat{\beta}$ as a function of the temperature 333We use $\hat{\beta}$ instead of the standard symbol $\beta$ to avoid confusion with the inverse temperature.. In Sastry et al. (1998), for Lennard-Jones liquids, $\hat{\beta}$ was shown to become sensibly different from 1 around $T_{\text{onset}}$. While the long time behavior of the correlation function is a pure exponential, we can fit the decay at intermediate times. In Fig. 9 we see that the stretched exponential form provides a decent fit if we do not approach too much close $T_{\text{\tiny MCT}}$.

We plot the stretching exponent $\hat{\beta}$ for the mixed model as a function of the temperature in Fig. 10, together with the same exponent in the pure model. In both cases we observe that while $\hat{\beta}\approx 1$ at high enough temperature and it starts to decrease from that value as temperature is lowered. We could not find any qualitative difference here in the behavior of the pure and of the mixed model, and we observe that in the mixed model the stretching starts for $T\simeq 1.2$, well above the estimated value of $T_{\text{onset}}\simeq 0.91$. Concluding, we do not see any signature of $T_{\text{onset}}$ in the equilibrium dynamics of the mixed model.

Appendix B Numerical extrapolations

The dynamical equations for the correlation and the response functions, written in Eq. (4), were integrated via a simple Euler algorithm with a fixed integration step $\Delta t$. In this way we get extremely precise results at finite time, but the times we can reach are limited. We tried more sophisticated integration schemes Kim and Latz (2001), where the integration steps is increased during the evolution but unfortunately, as soon as a mixture is used we have met numerical instabilities at very short times. In practice only the pure $p$-spin model seems to allow those integration schemes to work.

The results presented in the main text have been obtained integrating with an optimal integration step $\Delta t=0.1$, which allows us to reach the largest times without facing any numerical instability (consider that at short times the differential equations we are solving have a natural times scale of order 1, and so $\Delta t$ cannot be much larger than the value we used).

The integration error in the Euler algorithm is linear in the integration step $\Delta t$, and we are interested in understanding how much physical quantities computed with $\Delta t=0.1$ do differ from the corresponding $\Delta t\to 0$ limit. To this purpose we study the relative errors in one-time quantities defined as follows

[TABLE]

By definition this relative error is zero for $\Delta t=0.1$, should be linear in $\Delta t$ for $\Delta t$ small enough and in the limit $\Delta t\to 0$ provides the relative error committed in using an integration step $\Delta t=0.1$ instead of the exact integration ($\Delta t\to 0$).

In Fig. 11 we show the relative errors in some one-time quantities, $C(t,0)$, $\chi(t,0)$, $E(t)$ and $\mu(t)$, obtained in the integration of the dynamical equations for the (3+4)-spin with $T=T_{\text{\tiny MCT}}$. We clearly see that at any time $t$ the relative error scales roughly linear with $\Delta t$. For any time $t$ the relative errors are extremely small: of order $10^{-5}$ for correlation and response, and of order $10^{-7}$ for energy and radial reaction. Moreover the coefficient of the linear relation does not seem to grow indefinitely with time, but actually show a maximum around $t=O(100)$.

In order to better appreciate the non monotonic behavior of the relative errors, we plot in Fig. 12 the integration error in the energy obtained starting from $T=T_{\text{\tiny MCT}}$ as a function of time. The fact that relative errors decrease at very large times is probably related to another important observation: trajectories integrated with different $\Delta t$ may differ sensibly at short times, but at large times they seem to converge to the same asymptotic solution, which thus are very attractive and stable (this may be the reason why integration errors at very large times are smaller than those at intermediate times). The results of the integration shown in the present work are very reliable and stable.

Despite the high precision at finite times, the extrapolation to large $t$ of physical quantities, notably $\mu(t)$, $E(t)$ or $C(t,0)$ is delicate. While $C(t,0)$ decays very slow and we do not attempt any extrapolation in the $t\to\infty$ limit, both $\mu(t)$ and $E(t)$ converge fast enough and their large time limits can be estimated.

We already noticed in Fig. 4 that for $T\geq T_{\text{\tiny SF}}$ the radial reaction $\mu(t)$ is perfectly compatible with an extrapolation to the marginal value $\mu_{mg}=6$. At very high temperatures the decay is a simple power law, while at lower temperatures we observe pre-asymptotic effects, and an eventual crossover to the asymptotic decay. The latter seems to be characterized by an exponent which is roughly temperature independent.

Given the critical character of the dynamics in the aging regimes I and II, we have estimated the asymptotic energy $E_{\infty}(T)$ via power law fits to $E(t)$ data and we show in Fig. 13 the results. While at high temperatures the data follow a nice power law in time for about a couple of decades, at lower temperatures the presence of the pre-asymptotic regime and the crossover to the asymptotic one, makes very hard to assess the reliablity of the value of $E_{\infty}(T)$ that we extrapolate in this way.

However, comparing the data shown in Fig. 4 for $\mu(t)$ and in Fig. 13 for $E(t)$ we notice a very similar behavior (including the crossover) and thus it could be useful to plot parametrically $E(t)$ versus $\mu(t)$ varying $t$, and see whether a better estimation of $E_{\infty}(T)$ could be obtained. This is done in Fig. 14. The behavior of the energy as a function of the radial reaction is very smooth, practically linear at high temperatures and close to quadratic approaching $T_{\text{\tiny SF}}$. We can then fit the relation $E(\mu)$ and obtain a very good estimate of $E_{\infty}(T)$ just assuming that $\mu(t)$ converges to the marginal value $\mu_{mg}$ in the large time limit. We have interpolated the data via the function

[TABLE]

reporting with dotted curves in Fig. 14 the best fits.

In Fig. 15 we show the two estimates of $E_{\infty}(T)$ obtained from the procedures just described: the power law fit of energy versus time and the power law fit of energy versus radial reaction $\mu_{mg}-\mu(t)$. We see they are perfectly compatible.

Although we do not want to estimate the large time limit of $C(t,0)$ directly extrapolating in time (in Section E an estimate will be provided based on a different approach), we believe it is useful to show the $C(t,0)$ curves such that the reader can make her own idea. In Fig. 16 the data for $C(t,0)$ is shown, measured both in the pure 3-spin model and in the mixed (3+4)-spin model. We have decided to show these data as a function of the scaling variable $t^{-1/3}$ that describes pretty well the decay starting from a random configuration ($T=\infty$). En passant, we notice that the exact value for this decay is not known. We see in Fig. 16 that, while data for the 3-spin model can be easily extrapolated to zero in the large times limit, the data for the (3+4)-spin model seems much more compatible with a non-zero limit $\lim_{t\to\infty}C(t,0)$ when the temperature gets close to $T_{\text{\tiny MCT}}$.

The reason why we have not been able to integrate the dynamical equations for a longer time, resides in the difficulties in using integration algorithms with adaptive integration steps. This kind of adaptive algorithms have been used in the past to integrate the dynamical equations in the case of the pure p-spin model starting from random initial conditions Kim and Latz (2001); Berthier et al. (2007). Unfortunately, we uncovered that this adaptive algorithm is not robust to changes in the details of the dynamics, and its stability strictly depends on the regime under consideration. At each contraction of the grid, some information on the correlation and response functions computed at the previous steps is lost. This does not seem to be relevant if one starts from a random state and memory of the initial conditions is lost in the asymptotic state. However, in the mixed p-spin model, where memory of the initial configuration crucially affects the large time behavior, the error induced by the time-step adaptation give rise to strong numerical instabilities. Increasing the grid size pushes the instability to later times, but unfortunately, even the largest grid size we could used ($8192\times 8192$) did not allow substantial improvements with respect to the simple Euler algorithm.

Appendix C Response at short times

A natural hypothesis about the memory effect in region II is that it could be associated with an increased response $R(t,s)$ to perturbations applied at short times $s$. Here we show that short time dynamics provides hints against this hypothesis. Let us then compare the response $R(t,s)$ starting from infinite temperature and the one for an initial temperature in the region II. In Fig. 17 we plot the response as a function of the time difference $t-s$ starting from infinite temperature and from a temperature $T=1.01\;T_{\text{\tiny MCT}}$. We clearly see that starting from finite temperature the response at short times $s$ is much smaller than the one starting from infinite temperature. It is also interesting to study the behavior in time of $R(t,0)$, that we plot in Fig. 18. It is well known that starting from a uniformily distributed initial condition one has that the response to the initial condition coincide with the remanent magnetization, namely, $R(t,0)=C(t,0)$ Franz and Mezard (1994). For finite $T$ this relation does not hold, in fact it is possible to show that $R(t,0)=C(t,0)-\beta\int_{0}^{t}ds\;f^{\prime}(C(0,s))R(t,s)$. We see that a finite limiting value for $C(t,0)$ does not necessarily imply a finite response to the initial condition. Our data show that the two terms largely compensate and the larger $\beta$, the larger the compensation effect. Both for the mixed and the pure model, the lower the starting temperature, the lower $R(t,0)$. The decay of $R(t,0)$ becomes exponential only for $T=T_{\text{\tiny SF}}$. In the pure model this coincide with $T_{\text{\tiny MCT}}$, while in the mixed model $T_{\text{\tiny SF}}<T_{\text{\tiny MCT}}$.

Appendix D Asymptotic Solutions

In Ref. Barrat et al. (1997) it was shown that the asymptotic solution to the dynamical equations assuming a simple memorious aging within a state with a unique effective temperature has parameters $\chi$, $q_{0}$, $q_{12}$ and $y$ that satisfy the same equations that can the obtained extremizing the Franz-Parisi (FP) potential computed in the one step replica symmetry breaking (1RSB) scenario, provided that the parameter $y$, is fixed by a marginality condition.

In our case the 1RSB FP potential has to be computed at zero temperature ($T_{f}=0$), while the reference configuration is in equilibrium at temperature $T=1/\beta$

[TABLE]

The saddle point equations and the marginality condition thus read

[TABLE]

We notice en passant that equations in (16) do not always select minima of the FP potential, because for that one needs to set to zero the total derivative with respect to $q_{12}$ and not the partial derivative. The energy and the radial reaction are given respectively by

[TABLE]

The equations in (16) always admit the solution with $q_{12}=q_{0}=0$, representing the memoryless or ‘amnesic’ aging solution with parameters

[TABLE]

and consequently energy and radial reaction are given by

[TABLE]

For the (3+4)-spin model the numerical values for the above parameters are $\chi_{mg}=1/3$, $y_{0}=11/21\simeq 0.52381$, $E_{th}=-71/42\simeq-1.69048$ and $\mu_{mg}=6$.

A non trivial ‘memorious’ aging solution with strictly positive values for $q_{12}$ and $q_{0}$, disconnected from the amnesic aging solution, exists only for temperatures below $T_{0}$. Moreover at an even lower temperature $T_{\text{\tiny SF}}$ we have that $q_{0}\to 1$ and this solution becomes replica symmetric. In Fig. 19 we report the values of $q_{0}$, $q_{12}$ and $y$ in the memorious aging solution. In the (3+4)-spin model the existence domain for this solution is upper bounded by $T_{0}=0.8031557$ where the solution disappears by a square root singularity and lower bounded by $T_{\text{\tiny SF}}=0.7982754$, where $q_{0}\to 1$ and the solution becomes replica symmetric (RS). At $T=T_{\text{\tiny SF}}$ the parameters extremizing the potential can be computed analytically. The value of $q_{12}$ is more easily obtained from the RS solution (see below), while the value of $y$ can be obtained from the third order expansion in $\varepsilon=1-q_{0}$ of the difference between the two equations in (16), thus getting

[TABLE]

The numerical integration of the off-equilibrium dynamics does not give any indication in favour of this solution, in particular for the (3+4)-model, this solution has $y\approx 0.3$ in the whole range of validity, which is incompatible with the value $y\approx 0.52$ we obtain from the numerical solution of the dynamics. Even if we were disposed to believe that the dynamics would slowly cross-over to this solution on time scales we cannot observe numerically, still this would not solve the puzzle of the behavior of the asymptotic energy going below threshold for $T>T_{0}$. In this solution, the value of $q_{0}$ tends to 1 at temperature $T_{\text{\tiny SF}}$, well above the temperature of the static (Kauzmann) phase transition of the model. The one above is the only exact aging solution that we found; we searched, without success, more complex solutions with more than one effective temperatures.

When $q_{0}\to 1$ the 1RSB potential reduces to the RS potential

[TABLE]

and saddle point equations read

[TABLE]

with the energy taking the value $E=-\chi f^{\prime}(1)-\beta f(q_{12})$. In the 3+4 model, the energy of this solution at $T_{\text{\tiny SF}}$ is equal to $E_{\text{\tiny SF}}=-1.69981$ and decreases below. The susceptibility $\chi$ takes the marginal value $\chi_{mg}$ at $T_{\text{\tiny SF}}$ and becomes smaller than this value at lower $T$, as shown in the upper panel of Fig. 20. Our estimate $T_{\text{\tiny SF}}=0.7982756$ comes exactly from imposing $\chi=\chi_{mg}$ in this RS solution. In the lower panel of Fig. 20 we plot the value of $q_{0}$ and $q_{12}$ obtained in the RS ‘state following’ solution on the left of the vertical axis, which is located exactly at $T=T_{\text{\tiny SF}}$. On the right of the vertical axis the report the corresponding values computed in the aging RSB solution (the same plotted in Fig. 19). We observe that $q_{12}$ is a smooth function at $T_{\text{\tiny SF}}$. Its value right at $T_{\text{\tiny SF}}$ can be computed analytically from the first equation in (20) by imposing $\chi=\chi_{mg}$, leading to

[TABLE]

For $T\leq T_{\text{\tiny SF}}$ the FP potential has a secondary minimum described by a replica symmetric ansatz. So we are effectively describing the quenching process from temperature $T$ to zero temperature as a ‘state following’ process Sun et al. (2012): the observation that both at the starting and ending temperatures the state we are following is well described by a replica symmetric ansatz is an evidence that static-dynamics equivalence hold in this case. And this is what we indeed observe by comparing the asymptotic dynamics obtained integrating numerically the dynamical equations to the values derived here from the thermodynamic FP potential. For $T\leq T_{\text{\tiny SF}}$ aging disappears and one finds a simple agingless relaxation within a state described by the parameters computed in the secondary minimum of the FP potential.

Appendix E Semi-empirical relations to describe the large times dynamics

We have seen that there is no simple aging solution to the asymptotic equations with $q_{12}>0$ that is consistent with the finite time numerical integration of the dynamical equations. Nonetheless the FD plot in the mixed (3+4)-spin model do not differ in any appreciable way from the usual ones of the pure $p$-spin model with a single effective temperature and a value of $y$ equal to the one that holds for the $\beta=0$ initial condition. Insisting in representing the best that we can the data in the hic sunt leones region II with a single effective temperature dynamical ansatz, we ask what kind of relation can be derived between the overlaps $q_{0}$ and $q_{12}$. We suppose that $\mu$, $\chi$ take the marginal values $\mu_{mg}\equiv 2\sqrt{f^{\prime\prime}(1)}$ and $\chi_{mg}\equiv 1/\sqrt{f^{\prime\prime}(1)}$. Moreover we assume the FD slope in the range $C(t,t^{\prime})\in[q_{0}(t),1]$ is independent of the initial temperature, thus it is equal to value obtained with $\beta=0$, i.e. $y_{0}\equiv\sqrt{f^{\prime\prime}(1)}/f^{\prime}(1)-1/\sqrt{f^{\prime\prime}(1)}$. Finally we assume the response is null in the range $C(t,t^{\prime})\in[C(t,0),q_{0}(t)]$, where the correlation decays extremely fast.

Within this 1RSB ansatz the asymptotic value for the radial reaction is given by

[TABLE]

and by imposing $\mu=\mu_{mg}$, $\chi=\chi_{mg}$ and $y=y_{0}$ we get the relation

[TABLE]

that must hold on the asymptotic solution. We can notice at this point, that using (17) the energy can be written as $E=E_{th}+\Delta E$ with $\Delta E=yf(q_{0})-\beta f(q_{12})$. Only in the pure model Eq. (22) implies $\Delta E=0$, in all the other cases it gives a non vanishing $\Delta E$.

A second relation between $q_{0}$ and $q_{12}$ can be derived from the observation of the FD plots in Fig. 21, where each curve is a parametric plot of $\chi(t,t^{\prime})$ versus $C(t,t^{\prime})$ at fixed $t$ varying $t^{\prime}$. We notice that, below the onset temperature, the points \big{(}C(t,0),\chi(t,0)\big{)} follows closely the dash-dotted curve that we are going to derive analytically. It is important to remark that the dash-dotted curve is the same in all the panels of Fig. 21, so it is a temperature independent relation. Also the dashed line is temperature independent and corresponds to the FD relation in the aging solution with $\beta=0$, that is a line of slope $-y_{0}$ connecting the points $(0,\chi_{mg}+y_{0})$ and $(1,\chi_{mg})$.

In order to obtain the dash-dotted line we assume that in the aging regime (i.e. for $\chi(t,0)\geq\chi_{mg}$) the relation between $C(t,0)$ and $\chi(t,0)$ is linear, while in the regime $\chi(t,0)<\chi_{mg}$ the dynamics is asymptotically exploring a state and thus we assume the relation that holds within the RS ‘state following’ solution. The latter can be easily derived from the first equation in (20)

[TABLE]

that thus holds for $q_{12}\in[q_{12,\text{\tiny SF}},1]$, while the linear part has slope $-y_{0}/q_{12,\text{\tiny SF}}$ and passes through the points $(0,\chi_{mg}+y_{0})$ and $(q_{12,\text{\tiny SF}},\chi_{mg})$. We notice that the dash-dotted curve has a continuous first derivative at the point $(q_{12,\text{\tiny SF}},\chi_{mg})$, as can be easily checked by taking derivatives.

Our asymptotic ansatz thus implies a very simple relation between the overlaps describing the asymptotic aging regime, namely

[TABLE]

and plugging this relation inside Eq. (22) it is easily to find that a solution with $q_{12}>0$ can exists only if

[TABLE]

where $k$ is defined by $f(q)\propto q^{k}$ for $q\to 0$ (in the 3+4-model $T_{\text{onset}}=0.91$). Despite this is not an exact solution of the asymptotic equations, it is a strong indication that there is a phase transition between a memoryless phase where dynamics decorrelates from the initial condition and falls over the ‘usual’ threshold states with $E=E_{th}$ and a phase where aging takes place in a confined space, with an asymptotic energy below threshold and depending on $T$.

In Fig. 22 we report an even stronger evidence that this approximate solution provides a very good description of the asymptotic dynamics obtained by numerically integrating the dynamical equations. We plot numerically extrapolated values as data points with errors, while full lines are prediction from the approximate solution presented in this Appendix and dashed lines correspond to the standard 1RSB aging solution discussed in the previous Appendix. Needless to comment on which one is better in describing the asymptotic dynamics.

Appendix F Counting the minima

In this Appendix we try to answer to the question whether the attractors of the dynamics can be well described in terms of typical marginal saddles and minima that lie close to the initial configuration. Let us consider the stationary points of the Hamiltonian $H[\sigma]$ on the sphere $\sum_{i}\sigma_{i}^{2}=N$:

[TABLE]

As in dynamics, the radial reaction $\mu$, takes, in any stationary point, the value

[TABLE]

We wish to classify the stationary points according to their energy $E=\frac{1}{N}H[\sigma]$ and the value of the radial reaction $\mu$.

Differently from the pure models where $\mu=pE$, the relation between $E$ and $\mu$ here is not univocal and stationary points are found in a whole region of the $(E,\mu)$ plane.

We have seen that for $T>T_{\text{\tiny SF}}$ dynamics is attracted by some family of marginal minima. In order to characterize these minima, we count the number of stationary points of the Hamiltonian $H[\bm{\sigma}]$ with fixed energy $E$ and radial reaction $\mu$ that lie at a fixed overlap $\bm{\sigma}\cdot\bm{\sigma}_{0}=Nq_{12}$ from a reference configuration $\bm{\sigma}_{0}$, weighted with a Gibbs measure $e^{-\beta H[\bm{\sigma}_{0}]}$. Since the complexity, i.e. the logarithm of their number, is self-averaging, we write

[TABLE]

The computation of $\Sigma$ is standard, and can be performed in several steps. First of all, since the matrix $\bm{H}^{\prime\prime}$ is a GOE random matrix, the distribution of eigenvalues of $\mu\bm{I}+\bm{H}^{\prime\prime}$ is self-averaging and is a shifted semicircular $\rho(\lambda)=\frac{2}{\pi 4f^{\prime\prime}(1)}\sqrt{4f^{\prime\prime}(1)-(\lambda-\mu)^{2}}$. The modulus of the determinant $|\det(\mu\bm{I}+\bm{H}^{\prime\prime})|$ is self-averaging and its logarithm reads

[TABLE]

which only depends on $\mu$. The imaginary part of this expression is the proportion of negative eigenvalues.

To evaluate the remaining terms we use replicas and write $\Sigma\equiv\overline{\log(\mathcal{N})}=\lim_{n\rightarrow 0}\frac{\overline{\mathcal{N}^{n}}-1}{n}$. We concentrate to the case of temperatures greater that the static transition temperature ($T>T_{K}$) of the model where the partition function appearing in the denominator of (28) is self-averaging and takes its annealed value $\overline{Z_{\beta}}=e^{\frac{N}{2}\beta^{2}f(1)}$. One can then average over the disorder and the configuration $\sigma_{0}$ at the same time. Opening the delta function in the Fourier basis,

[TABLE]

where $\int\mathscr{D}s=\int_{S}\mathscr{D}\bm{\sigma}_{0}\prod_{a}\left(\int_{S}\mathscr{D}\bm{\sigma}_{a}\int\mathscr{D}\hat{\bm{\sigma}}_{a}\int\hat{\beta}_{a}\right)$.

And since the disorder is Gaussian :

[TABLE]

Now we define overlap variables $NQ_{ab}=\bm{\sigma}_{a}\cdot\bm{\sigma}_{b}$, $N\chi_{ab}=i\bm{\sigma}_{a}\cdot\hat{\bm{\sigma}}_{b}$ and $NR_{ab}=-\hat{\bm{\sigma}}_{a}\cdot\hat{\bm{\sigma}}_{b}$, and the overlaps with the reference configuration $Nq_{12}=\bm{\sigma}_{a}\cdot\bm{\sigma}_{0}$, $N\chi_{p}=i\hat{\bm{\sigma}}_{a}\cdot\bm{\sigma}_{0}$. This change of variables defines a matrix:

[TABLE]

where $Q_{aa}=1$ due to the spherical constraint. And, from the equivalence between replicas, we fix $i\hat{\beta}_{a}=y$ and $\chi_{aa}=\chi$ $\forall a$. With this change of variables Eq. (30) becomes

[TABLE]

where $\frac{1}{2}\log(\det\mathscr{Q})$ is the volume factor that comes from the change of variables from spins to overlaps. To get the leading $N$ contribution we must extremize with respect to all the overlap parameters $\mathscr{Q}$ and $y$. We notice that a further simplification to the expression (31) can be obtained by firstly extremizing with respect to $R_{ab}$. Assuming a replica symmetric ansatz for the overlap matrices $Q$ and $\chi$, i.e. $Q=\delta_{ab}+(1-\delta_{ab})q_{0}$ and $\chi_{ab}=\delta_{ab}\chi+(1-\delta_{ab})\chi_{1}$ we get, in the limit $n\rightarrow 0$:

[TABLE]

where

[TABLE]

This can be extremized explicitly with respect to $y,\chi,\chi_{1},\chi_{p}$, while $q_{0}$-extremization has to be done numerically. For $q_{12}=0$ the solution is $q_{0}=0$ and we recover the expression found in Arous et al. (2018)

[TABLE]

This unconstrained complexity is plotted in Fig. 1. For any value of $\mu$, $\Sigma(E,\mu)$ has a concave parabolic shape as a function of the energy and many values of the energy are possible for the same $\mu$ value. In particular this is true for $\mu=\mu_{mg}$: there is not a unique threshold energy (at variance to pure $p$-spin models), but a whole interval for which $\Sigma(E,\mu_{mg})>0$. We define the dominating stationary points at a given value of the energy as the ones that maximize $\Sigma$ as a function of $\mu$. The threshold energy as defined from dynamics, corresponds to the values of the energy that separates minima from saddles on the dominating line. Notice that this value does not correspond to the most numerous marginal minima, which occurs for $\mu=\mu_{mg}$ and $E>E_{th}$, but to the point where the most numerous stationary saddles become minima. This observation sheds some light on the dynamics from a random initial condition; minima are not accessible at levels of the energy $E>E_{th}$ where saddles dominate the landscape, even if many minima are present.

Extending the concept of threshold energy to the case $q_{12}>0$ is straightforward. As shown in Fig. 23 for any fixed value of $q_{12}$ we can plot the curve corresponding to dominating minima. This curve ends at an energy value that we call $E_{th}(q_{12})$ and this is the best candidate for the asymptotic energy in the relaxation dynamics. Indeed for $E>E_{th}(q_{12})$ dominating stationary points are saddle and with high probability the dynamics will be able to relax further until the energy level $E_{th}(q_{12})$ is reached. However, even assuming the dynamics converges to $E_{th}(q_{12})$, we are still left with an unknown: the value of $q_{12}$, i.e. the asymptotic value for $C(t,0)$. This is the same problem the authors of Ref. Capone et al. (2006) faced after computing the constrained complexity for the free-energy at a non-zero temperature.

In Fig. 24 we plot the quenched constrained complexity of dominant states in the (3+4)-spin model, computed with $T=T_{\text{\tiny MCT}}$ (left panel) and $T=1.02\,T_{\text{\tiny MCT}}$ (right panel). Please notice that the choice of the temperature is not crucial: as long as $T\in(T_{\text{\tiny SF}},T_{\text{onset}})$ the plot would be very similar. The almost-horizontal red curve in the plot represents the threshold energy $E_{th}(q_{12})$ defined above. We notice that the range of energies with a positive complexity is very large compared to variations in $E_{th}$. The green ellipses represents our best estimate for the large time limit of the actual dynamics solved numerically: while the energy can be very well estimated, the limit of $C(t,0)$ is plagued by a large uncertainty due to its slow convergence (until now we have not attempted such an extrapolation, but here we want to be more speculative and we take the risk). From the plot one is tempted to conjecture that the dynamics always converges to marginal states with a threshold energy $E_{th}(q_{12})$. However, we have not found any principle to fix the value of $q_{12}$ solely from the complexity curve and further studies are needed to better match the large time limit of the dynamics to the energy landscape.

Although Fig. 24 may suggest a relation between dynamics and the generalized threshold energies, a more careful analysis reveals its limitations. Assuming that at large times the relaxation dynamics converges to the manifold of marginal states belonging to the curve $E_{th}(q_{12})$ one could estimate the point reached by the dynamics by extrapolating the asymptotic energy $E_{\infty}(T)$ and estimating $q_{12}$ from the equality $E_{th}(q_{12})=E_{\infty}$. Having thus fixed the values of $E$ and $q_{12}$ one can proceed estimating the remaining parameters of the asymptotic aging dynamics, $q_{0}$ and $y=\partial_{E}\Sigma$, from the saddle point equations used to compute the quenched complexity. The result of this computation is shown in Fig. 25 with full lines, and compared to the (very uncertain) extrapolation of $C(\infty,0)$, shown by red points with errors. We clearly see that, while the estimate of $q_{12}$ is compatible with the actual dynamics, the other two parameters are far from the values measured in the numerical solution to the dynamics. Indeed $q_{0}$ becomes smaller than $q_{12}$ while in the actual dynamics the inequality $q_{12}<q_{0}$ is always satisfied, and $y$ becomes much smaller than $y_{0}$ (marked by a dotted horizontal line in Fig. 25) which is good descriptor of the actual aging in the whole temperature range studied. So, the present computation of the quenched complexity of marginal minima of the energy function in mixed $p$-spin models does not allow to identify the attractors of the dynamics. More work will be required to connected the large time aging dynamics to the properties of the energy landscape.