Probing chaos in the spherical p-spin glass model

We study the dynamics of a quantum $p$-spin glass model starting from initial states defined in microcanonical shells, in a classical regime. We compute different chaos estimators, such as the Lyapunov exponent and the Kolmogorov-Sinai entropy, and find a marked maximum as a function of the energy of the initial state. By studying the relaxation dynamics and the properties of the energy landscape we show that the maximal chaos emerges in correspondence with the fastest spin relaxation and the maximum complexity, thus suggesting a qualitative picture where chaos emerges as the trajectories are scattered over the exponentially many saddles of the underlying landscape. We also observe hints of ergodicity breaking at low energies, indicated by the correlation function and a maximum of the fidelity susceptibility.


Introduction
Characterizing the dynamics of a system with multiple equilibrium configurations is a challenging problem encompassing several branches of natural sciences, from biology [1] to economy [2] and theoretical ecology [3].The nontrivial interplay between stable and unstable configurations makes the dynamical evolution of these systems extremely sensitive to the initial conditions and thus unpredictable: since the pioneering works of the 19th [4,5] and 20th centuries [6,7], this phenomenon has been known as chaos.
Among the physical systems displaying multiple equilibria, a prominent role is played by spin-glasses, a class of models originally introduced to describe magnetic alloys [8,9]: at low temperature, the dynamics is trapped into one of the exponentially many possible metastable states [10][11][12][13] for a long-time and explores the phase space through rare, activated jumps between different states [14][15][16], leading to ergodicity-breaking phenomena [17][18][19][20].The inclusion of quantum fluctuations usually opens up a new route for thermalization in the equilibrium dynamics due to tunnelling between metastable states [21][22][23][24], even though counter-intuitive effects of quantum fluctuations inducing glassiness has been found in quantum quench protocols [25] or more recently in presence of more complex energy landscapes [26].Within such a rich variety of phenomena, the recent theoretical [27][28][29][30][31] and experimental [32] observation of many-body chaos in quantum spin-glasses has drawn particular attention.Specifically, in Refs.[29,30] chaos was detected in the dynamics of a quantum spherical p-spin glass model [23], from the exponential growth of an out-of-time-order correlator (OTOC) [33][34][35][36][37][38][39]: the corresponding Lyapunov exponent λ L exhibits a single broad peak at a temperature scale where the dynamics is still ergodic, but at the onset of slow relaxation.
The behaviour of the Lyapunov exponent λ L extracted from the OTOCs was connected in Ref. [29] to a dynamical crossover between two step and one step time relaxation in the spin-spin correlation functions.This behaviour suggests a deeper connection between this dynamical feature and the structure of the underlying stationary configurations of the model.In particular, if such relation exists other features associated to the peak in other chaos indicators would be observed.Among these chaos indicators, one example is the Kolmogorov-Sinai (KS) entropy [40][41][42], that is the sum of all the positive Lyapunov exponents describing growth of the sub-leading terms appearing either in the OTOC for quantum systems [43,44] or in the distance between nearby trajectories for classical chaotic systems [45].The KS entropy would indeed provide a deeper insight on the emergence on the strength of chaos and to the entanglement production in a wide range of quantum systems [43,44,[46][47][48].A second, example of indicator focusing on the transition from integrability to chaos in genuine quantum spin systems is the fidelity susceptibility [49,50].Such quantity, despite being closely connected to the magnetic response obtained in quantum spin-glass experiments [21], has never been investigated in the context of spin-glasses.
The purpose of this work is to develop a qualitative understanding of chaotic behavior in the p-spin model by using phase space methods to compute λ L , at fixed energy and in the classical limit.We find that λ L exhibits a broad peak around E = 0, at the onset of slow dynamics, and vanishes both at low and high energies, compatibly with previous studies.Working at fixed E we observe that the maximum in λ L occurs when the energy landscape has the maximum possible number of stationary points, as characterized by the complexity (see Ref. [51]).We find that the profile of the Kolmogorov-Sinai entropy is qualitatively identical to that of λ L .To further investigate the chaotic features of our model, we also introduce a quantity classically equivalent to the fidelity susceptibility χ in the ergodic phase finding that, as a function of E, it has a single peak, aligning with the onset of non-ergodic behavior in the dynamics.Our results can be easily interpreted in terms of the typical behaviour of the underlying trajectories at different energy scales, which in turn determine the profile of the correlation function: While the low and high-energy limits respectively correspond to trajectories performing small oscillations around a local minimum or a uniform circular motion on the N -sphere, for intermediate energies the nontrivial interplay between the saddles of the landscape enhance dynamical chaos.Our method is straightforwardly generalized to other spin-glass models, for example describing spins 1/2 in a transverse field [28], where a similar behaviour for the trajectories is expected both for the low and high-energy limit.

Classical dynamics of the p-spin glass spherical model
Throughout this work, we will focus on the isolated dynamics of the p-spin glass Spherical Model (PSM), whose Hamiltonian with describes a system of N spins interacting through random, all-to-all couplings J i 1 ,...,i p , independently sampled from a Gaussian distribution with zero mean and variance J 2 = 2p!/N p−1 .
Here the spins σi are treated as continuous variables, obeying the spherical constraint i 〈 σ2 i 〉 = N [52], and quantum fluctuations are implemented by the canonical quantization relations [ σi , Πj ] = iħ hδ i j .The term "spin-glass" is attributed to the quantum PSM model, despite its use of position and momentum operators, due to historical reasons.The introduction of the quantum PSM can be traced back to Ref. [23], where it was revealed to manifest a thermodynamic phase transition from a paramagnetic to a spin-glass state, driven by the temperature T and by a dimensionless parameter Γ = ħ h 2 /M J, which quantifies the strength of quantum fluctuations.The transition line Γ c (T ) that separates these two phases shares qualitative features with the glass transition line experimentally observed in disordered spin 1/2 systems [21,53], coupled to quantum fluctuations through an homogeneous transverse field.The nomenclature "spin-glass" was then associated to the quantum PSM due to these observed similarities.The thermodynamic phase transition in the PSM is either of the first or second order depending on the strength of Γ [23,24] and its corresponding transition line can be parametrized also in terms of by a critical temperature T c (ħ h).At temperature T d (ħ h) ≳ T c (ħ h), a dynamical, ergodicity-breaking transition is also observed [22], whose properties are sensitive to quenches in ħ h in a non-trivial way [25].In refs.[29,30], many-body chaos in the quantum PSM was studied using an OTOC in the large-N limit and for a unitary evolution of the system from a thermal initial state.The OTOC grows exponentially at any temperature T , with a quantum Lyapunov exponent λ L (T ), which displays qualitatively the same profile for a wide range of fixed values of ħ h and in particular in the classical limit ħ h → 0, having a single maximum at T m (ħ h) > T d (ħ h) and vanishing in the low and high-temperature limits.In contrast to its fermionic counterpart, the Sachdev-Ye-Kitaev (SYK) model [34,54,55], in the PSM the bound on chaos λ L ≤ 2πT /ħ h, proven for a general quantum many-body system in Ref. [37], is never saturated.A similarity between the Lyapunov exponent of the SYK model and that of the PSM only arises for ħ h → 0, where the bound becomes trivial.In this case, the λ L (T ) grows linearly with temperature T at low T for both cases [29,56].
The aim of this work is to gain physical insight on the behaviour of λ L investigating the dynamics of the PSM in the classical limit ħ h → 0 and at fixed energy-density E, a framework usually used for classical dynamical systems [42].This goal can be achieved in the framework of Truncated Wigner Approximation (TWA) [57,58], briefly reviewed in the following (see Refs. [59,60] for further details).We begin by observing that, for a generic system with Ndimensional position and momentum-like degrees of freedom, like σ and Π, the Heisenberg dynamics from an initial state ρ of any operator Ô = O( σ, Π) can be represented as where is called Weyl symbol of the operator Ô, while the Wigner function W ρ (σ, Π) is the Weyl symbol of the density matrix ρ.The TWA is then a classical approximation for the dynamics of O W (σ, Π, t), summarized as where (σ(t), Π(t)) is the classical trajectory evolving from the initial condition (σ, Π) and determined by the following Hamilton equations: The Eqs. ( 6) are essentially derived by adding to the Hamiltonian a Lagrange multiplier term z(t)(σ 2 − N ), where z(t) is determined in a self-consistent way so that the dynamics satisfies the spherical constraint at any time t (see Appendix A for more details).The TWA is believed to be valid at least up to an Ehrenfest time t Ehr ∼ log ħ h −1 [33] diverging for small-ħ h.
To investigate the dynamics at fixed energy density E, the ideal setup would be to fix ρ as a micro-canonical state and compute the corresponding Wigner function, which is in general a formidable task.Instead, we fix an initial condition as an ensemble of coherent states , defined as eigenstates of the operators âj = σj /( 2l) + il Πj /( 2ħ h), for j = 1, . . ., N and a free parameter l, with eigenvalues N }, respectively.For each state |α〉, the Wigner function is the Gaussian wave-packet In the following, we fix l = ħ h, to have the same uncertainty in the variables σ and Π.The centers (σ α , Π α ) of the various wave-packets are chosen in a way that 〈α| ĤJ |α〉 /N = E.For the small values of ħ h we are interested in, it can be easily shown that such choice is equivalent to fix the classical energy density so that we can determine the centers of the wave-packets in Eq. ( 8) using the following "classical annealing" algorithm: 1. First extract a random configuration (σ 0 , Π 0 ) in phase space, with σ 0 uniformly sampled on the N -sphere and Π 0 sampled from the normal distribution with zero mean and unit variance.
2. To bring the system in a configuration at the desired energy density E, we integrate the dynamics starting from (σ 0 , Π 0 ) and defined by the equations where a dissipative term of strength γ > 0 has been added.Notice that γ > 0 if H cl (σ 0 , Π 0 ) > N E and γ < 0 otherwise.
3. Finally stop the integration as soon as the system reaches a configuration (σ 1 , Π 1 ) such that H cl (σ 1 , Π 1 ) = N E. Afterwards we may set γ = 0 and integrate the Hamilton dynamics (Eq.( 6)) from (σ 1 , Π 1 ) for a time t eq , to let the system reach a typical configuration on the corresponding classical microcanonical manifold, which we finally take as the center (σ α , Π α ) of the wave-packet in Eq. ( 4).Throughout the rest of this work, we fix γ = 0.5 and t eq = 5.
In practice, the Wigner function for the initial state in Eq. ( 7) is obtained taking the average over N s different wave-packets, sampled from the classical annealing algorithm for the same fixed configuration of the disorder {J i 1 ,...,i p }, as We repeat the algorithm for each of the N s states: the resulting set of points {(σ (l) c , Π (l) c )} l=1...N s is then a non-uniform sampling of the classical microcanonical manifold at energy density E. As long as ħ h is very small, TWA enables us to sample orbits evolving from a neighborhood of each of each configuration (σ α , Π α ), a feature which we will use in Section 3 to investigate classical chaos in the PSM from the average growth of the width of the wave-packets.We compute all the observables we are interested in by averaging over an ensemble of trajectories, evolving according to Eq. ( 6) from an initial condition (σ, Π) sampled from the distribution in Eq. (11).In the rest of this work, we will denote by 〈 • 〉 the average over the sampled trajectories and  12), reported for few paradigmatic energy densities .The average over the initial condition is performed extracting 25 random configurations from each of the N s = 20 wave-packets, defined in Eq. ( 8) and obtained through the classical annealing algorithm.The dynamics is integrated up to a time t max = 80.The logaverage over the disorder is taken over N d = 30 configurations.(Right) Lyapunov exponent λ L , defined in Eq. ( 13), obtained through a linear fit of log d J (t) over a time window [t I , t F ], defined in such a way that, for each E, t I is beyond the oscillations at early times displayed by log d J (t) and t F is smaller than the saturation time-scale.by ( • ) the one over the disorder configurations {J i 1 ...i p }.We fix ħ h = 10 −8 , N = 100 and we focus on the paradigmatic case of p = 3.
We conclude this section with two technical remarks.First, it is worth noting that the classical annealing algorithm described above allows us to sample energies that are arbitrarily high, but not arbitrarily low.Specifically, in the PSM, the minima of the potential V J (σ) are situated below a characteristic energy scale of E th = − 2(p − 1)/p [10] (also discussed in Sec. 4 and detailed in Appendix E).Due to this limitation, our classical annealing approach cannot explore phase space configurations with energies E < E th , as the dissipative dynamics of Eq. (10) becomes trapped in the vicinity of the first encountered local minimum.As a result of this constraint on sampling low-energy configurations, we are unable to investigate the relationship between λ L and E in the vicinity of the ground state.Consequently, we are also unable to make a comparison against the linear dependence of λ L on T , for small T , observed in Ref. [29].The second remark is that our definition of the Wigner function is not rigorous for the degree of freedom σ lying on a compact configuration space (see also Ref. [62]).Specifically, the configurations σ sampled from the distribution in Eq. (11) are not confined to the N -sphere.Although our approximation fails in capturing precise quantum dynamics at finite ħ h, it is expected to be reliable in the limit of ħ h → 0. For our investigative purposes, it serves merely as a tool to examine the classical trajectories evolving from nearby initial configurations sampled at low ħ h (thus testing classical chaos).Furthermore, as ħ h tends to 0, the configurations we extract are expected to asymptotically lie with the N -sphere.This is true as long as the center σ c of each sampled wave-packet W α (l) (σ, Π) also resides on the N -sphere, a condition consistently met within the classical annealing algorithm.

Results: Chaos estimators
We present here our results for the chaos estimators of the PSM, starting from the Lyapunov exponent λ L .In principle, λ L can be obtained as a function of E by computing the Weyl symbol, defined in Eq. ( 4), of the OTOC [28].However, in the ħ h → 0 limit we are interested in, λ L becomes the exponential rate of divergence of pairs of nearby orbits of the emerging classical dynamics 1 and can be computed from where ŷ = ( σ1 , . . ., σn , Π1 , . . ., ΠN ) is the set operators corresponding to a classical phase space configuration.The quantity d J (t) is easier to compute than the OTOC in the TWA framework and, in Appendix B, we show explicitly that both d J (t) and the OTOC grow exponentially with the same rate in the classical limit.To get rid of the small time scale fluctuations appearing in the dynamics of log d J (t), we compute its average log d J (t) over the disorder, and retrieve a smooth linear growth log on intermediate time scales, as shown Fig. 1 (left).The corresponding Lyapunov exponent λ L , shown in Fig. 1 (right) against the energy density E, has a clear peak close to E = 0, while it asymptotically vanishes at low and high energies.In Refs.[29,30], it was shown that λ L has the same qualitative profile as a function of T , displaying for small ħ h a single maximum around T m (ħ h) ≃ 1: this maximum is consistent with the one we find at E = 0, since in the paramagnetic phase the classical energy density E and the temperature T are related by the equation E = T /2 − 1/(2T ) (as shown in Ref. [63]).Our computed results are also in close numerical agreement with those obtained in Refs.[29,30], where the estimated maximum for the Lyapunov exponent was around λ L ≃ 0.6, like in our findings.However, it is important to note that the order of operations in Eq. ( 13), involving a logarithm and an average over disorder, is reversed compared to previous studies (see Ref. [64] for a more general discussion).Consequently, we do not expect a perfect match between the λ L we compute here and results from Refs.[29,30].In summary, the exponent λ L we calculate is essentially classical, and the introduction of small fluctuations in the initial conditions is merely a tool we use to sample nearby trajectories starting from the same wave-packet.It is worth mentioning that, while we use quantum fluctuations to sample nearby configurations in phase space, classical chaos can also be probed using different kind of fluctuations (see Ref. [65] for an example).
In the classical limit, the strength of chaos at different energy densities can be also connected to entropy generation by looking at the Kolmogorov-Sinai (KS) entropy [40,41].This is defined from the observation that, in general, N exponential terms in the form of exp(λ i t) contribute to the growth of the distance in Eq. ( 12), with non-negative Lyapunov exponents hierarchically ordered as 66,67].The KS entropy (density) is then just the sum Λ KS = N i=1 λ i /N : classically, Λ KS quantifies the rate of spreading of coarse-grained volumes in phase space [68], while its quantum counterpart is believed to describe the entanglement growth at early times for a wide range of systems [46][47][48] and to be a stronger indicator of scrambling dynamics than the single λ L [44].As shown in Refs.[46,48], the eigenvalues of the symmetric fluctuation matrix 1 It is important to notice that the classical limit of λ L we are investigating corresponds to twice the classical Lyapunov exponent λ cl [38], as the classical limit of the OTOC is actually the square of a typical distance between the underlying trajectories (see also Appendix B.1).  14) and ( 15), reported for few paradigmatic energy densities.The average over the initial condition is performed extracting 30 random configurations from each of the N s = 30 wave-packets, defined in Eq. ( 8) and obtained through the classical annealing algorithm.The dynamics is integrated up to a time t max = 80.The log-average over the disorder is taken over N d = 96 configurations.(Right) Kolmogorov-Sinai entropy Λ KS , defined in Eq. ( 15), obtained through a linear fit of S KS over a time window [t I , 10], where t I is beyond the scale of oscillations at early time displayed by S KS , for each E. diverge as exp(λ i t) in the limit ħ h → 0.Then, the KS entropy is straightforwardly extracted from the growth rate of the quantity on intermediate time-scales.We plot S KS (t) inf Fig. 2 (left) and, in Fig. 2 (right), its corresponding slope Λ KS , showing that the maximal chaos located at E = 0 can be detected also by the KS entropy.Notably, a similar result was recently derived also for a classical spin system without disorder [69].

Chaos, ergodicity and energy landscape
In this section we elaborate further on the results of the previous section and try to provide a qualitative interpretation for the observed maximal chaos around energy E = 0 in the PSM.
In particular, a natural question is whether this result can be understood from the relaxation dynamics of the PSM at fixed energy density, as probed for example from the spin correlation function, or related to properties of the energy landscape.As we are going to show, the maximal chaos around E = 0 in the PSM occurs when the spin relaxation is the fastest and when the complexity of the underlying energy landscape is maximal.
The relaxation dynamics for both classical [17,18] and quantum [22] spin glasses is most often studied in presence of a finite temperature bath.Relaxation is usually defined in terms of the (symmetric) correlation function At high temperatures, the function C(t w , t w + τ) becomes approximately time-translation invariant for moderately high values of t w , and decays to zero for large τ, indicating that the underlying dynamics of the system is ergodic.However, at sufficiently low temperatures, the system may exhibit the so-called 'weak ergodicity breaking scenario' [17,22,70] (see Ref. [18] for a review), meaning that lim with a finite dynamical overlap q 1 > 0 and where again C st (τ) vanishes for τ → ∞, determining a non-ergodic dynamics.In spin glasses, ergodicity breaking is usually accompanied by a breaking of time-translational invariance in C(t w , t w + τ) a phenomenon usually referred to as aging [18,71]: C(t w , t w +τ) has a plateau around q 1 , whose finite length increases as t w grows (and diverging for t w → ∞), before eventually decaying to zero for longer time-scales.Here we are interested instead in the Hamiltonian relaxation dynamics, starting from fixed energy initial conditions.We note that the Hamiltonian dynamics of both classical and quantum PSMs starting from a finite temperature state has been studied recently [25,63].To this extent we compute C(t w , t w + τ) in the TWA formalism at finite energy density E, making use of Eq. ( 3) and of the identity The results in Fig. 3-(a) show that the correlation function undergoes a temporal crossover from high energies, where it displays wide oscillations, to low energies, where the dynamics slows down.Quite interestingly, at the maximally chaotic point E = 0 we observe the fastest relaxation of the correlation function, decaying to zero with few oscillations, again compatibly with Refs.[29,30].Upon decreasing further the energy below E = 0 the dynamics slows down and we expect a finite plateau around q 1 to to appear.Detecting this intermediate plateau within the given simulation time-window is challenging.Consequently, we compute q 1 as the value of the correlation function at t w = τ = t max /2: our results are expected to converge to the ones predicted by Eq. ( 17) in the limit of t max → ∞.In Figure 3-(b) we show that q 1 becomes nonzero below a certain energy threshold, estimated to be around E d ≃ −0.38, and that it increases as E further decreases below the threshold.At the same time, the profiles of the correlation function shown in Fig. 3-(c) lose time-translation invariance again below E ≃ −0.38.Notably, these findings are also compatible with the ones obtained from the same simulations performed for a larger t max , which are reported in Appendix C and highlight that both a finite q 1 and a loss of time-traslational invariance are retrieved at the same energy scale.The results discussed above indicate signs of ergodicity breaking below the energy E d ≲ −0.38.However, we recognize that our findings, including our estimated value for E d , might undergo quantitative changes with a more extensive analysis using a significantly larger t ma x than the values considered in this manuscript.Therefore, the ergodicity breaking we observe here should be interpreted as a regime where the dynamics is sufficiently slow for the thermalization time to extend beyond t max .We now argue that both maximal chaoticity and fastest spin relaxation emerge alongside maximal complexity in the topology of the potential V J (σ), from Eq. (2), at the E = 0 level surface.To understand this connection, we first observe that the profile of the correlation function in Fig. 3-(a) can be associated to the typical behaviour of the underlying trajectories, as sketched in Fig. 3-(d).While the regular oscillations at high energies are due to an underlying uniform circular motion on the N -sphere, in the limit of low energies the trajectories oscillate in a well around a local minimum of V J (σ), whose amplitude can be roughly estimated as the typical distance between two configurations of the same trajectory, observed at large times  Asymptotic value q 1 of the correlation functions C(t w , t w + τ) reported in panel (a), obtained through the time-average of the latter over τ ∈ [39,40].(c) Comparison between the profiles of C(t w , t w = τ) obtained fixing different values of t w .Each panel correspond to a different energy density E. (d) Time-evolution of three typical orbits, whose initial condition are obtained extracting one configuration from the distribution in Eq. (11), at different energy densities E and for the same configuration of the disorder.The orbits evolve in a 200-dimensional phase space and are projected over the two axes defined by the initial spin configuration σ(0) and the initial momentum Π(0).
Chaos and relaxation emerge in between these two trivial limits, where the trajectories are scattered in neighbors of the stationary configurations of the dynamics in Eq. ( 6) and explore the whole configuration space.The stationary configurations can be defined as solutions of the equations (see also Ref. [72]) where in the first equation we used the identity j σ j ∂ V J /∂ σ j = pV J (σ), holding for the potential V J (σ) defined in Eq. ( 2).As discussed in Appendix E the average number of solutions of Eq. ( 20), lying on the microcanonical manifold at energy density , is in a one-to-one correspondence with the stationary points of the potential V J (σ) on the Nsphere.The average number of such stationary points can then be expressed as (see again Appendix E): For the classical potential V (σ) in Eq. ( 2) and in the large-N limit, the number of stationary points scales exponentially as N (E) ≃ exp{N Σ(E)} [10,51], where Σ(E) is usually referred to as complexity.At the same time, the stability of such stationary points is characterized by the index k(E), where N k(E) is the average number of unstable directions around every stationary points.In Appendix E, we also derive the analytical expressions for both Σ(E) and k(E) as functions of E, finding that: where and Intuitively, the value E th = − 2(p − 1)/p appearing in the previous formulas represents the threshold energy density below which the stationary points are typically local minima of V (σ), so that k(E) = 0 (similarly, −E th is the energy density above which all stationary points of V (σ) are typically local maxima).Both Σ(E) and k(E) are plotted in Fig. 4. The first observation that we make is that the complexity Σ(E) has a maximum on the E = 0 surface, where stationary points are predominantly saddles of V J (σ), surrounded on average by half stable and half unstable directions, as k(E = 0) = 1/2.Second, we notice that the complexity Σ(E) vanishes at two points E = ±E 0 , with E 0 > |E th |.Beyond E 0 , we have Σ(E) < 0, which implies a vanishing number of stationary configurations, so we interrupt the plot at E 0 .Intuitively, −E 0 (E 0 ) can be interpreted as the typical value of the absolute minimum (maximum) of V J (σ)/N [72], which is always finite for σ lying on the N -sphere.We observe that our results for the complexity at fixed energy are compatible with the ones in the literature obtained using different methods [51,73].The maximum of Σ(E) at E = 0 unveils a correlation between the number of saddles in the energy landscape and the maximal chaos, detected by λ L .This observation suggests that the scattering of trajectories against a maximal number of saddles could offer a potential explanation for the emergence of maximal chaos at E = 0. We conclude our analysis by discussing the connection between chaos and ergodicity in the classical PSM.In the TWA framework we employed, chaos is probed by the Lyapunov exponents while ergodicity is determined by the long-time behavior of the correlation function.However, recent works [50,74] have shown a possibly universal connection between ergodicity in quantum systems, defined there as the emergence of a random matrix behaviour [75], and the fidelity susceptibility [49], which characterizes the sensitivity of the eigenstates to an external perturbation.Specifically they shown that an ergodic Hamiltonian exhibits a scaling behavior of χ ∼ ω −1 L against the mean level spacing ω L , a behaviour which is consistent with random matrix theory predictions, while for an integrable [76][77][78][79][80] or disorder-localized [80][81][82][83] system has a value of order one.Then, a stronger divergence of the fidelity as χ ∼ ω −α L , with α > 1, indicates the approach to a non-ergodic (either integrable or disorder-localized) point.Such scaling is expected to be retrieved in a region around the non-ergodic point which become exponentially small in the system size, when taking the thermodynamic limit.Despite such a complete picture, the fidelity susceptibility has never been tested in classical systems.Here we use it as a tool to identify ergodicity in our p-spin model, where ergodicity and its breaking are controlled by the energy density E. In particular, here we perturb the Hamiltonian in Eq.( 1) with local magnetic fields B i , summarized in the extra term Then, by perturbation theory, the sensitivity χ (i) n of the n-th eigenstate to the magnetic field B i (posing B j = 0 on every other site j), is defined by and we define the fidelity susceptibility as its average χ over the initial condition in Eq. ( 7) and the disorder configurations: In principle, the fidelity defined in Eq.( 27) can be computed classical framework because, as already observed in Ref. [50], χ can be expressed in terms of the spectral function Cav (ω), that is the Fourier transform of the time-averaged correlation function The two are related by the equation (proven in detail in Appendix D): where ω L is again the average quantum level spacing.
In practice, two obstacles prevent us to use directly the Eq.( 29).The first one is that we do not have directly access to the spacing ω L .We then study regularized fidelity (used also in Ref. [74]) where we introduced the cutoff µ to suppress the contribution to the integral from frequencies |ω| ≲ µ and plays a role equivalent to the one played by ω L in genuinely quantum systems.We determine the asymptotic profile of χ µ by analyzing its behaviour as µ → 0. The second, more practical obstacle is that in our simulations we do not have access to infinite-time average, so that the integral in Eq. ( 28) is performed over a finite time-window [0, T ].While in the ergodic phase this is a good approximation of the long-time average, due to time-translation invariance, in the non-ergodic one χ always depends on choice of the time-window and converges to the definition in Eq. ( 27) only in the limit T → ∞.However, as shown in Appendix D, the qualitative profile we retrieve for χ µ is the same for a wide range of T between 0 and the maximum integration time t ma x , so that here we can focus on the specific case of T = t max /2.We plot χ µ , as a function of E and for several values of µ, in Fig. 5-(a): its profile has a peak close to the estimated ergodicity breaking energy scale E ≃ −0.38 and the maximum point has a little drifting to the left approaching small values of µ.We also observe that, in general, a natural low-frequency cutoff ∆ω ∼ 2π/t max emerges when discretizing the integral in Eq. (30) in our finite-time simulations, so that the asymptotic behaviour of χ µ can be studied only up to µ ≳ ∆ω.Thus, to refine our analysis, we compute χ µ for a dynamics integrated for a larger t ma x so that µ > 3∆ω for all the values of µ we investigate: the new results, shown in Fig. 5-(b), are qualitatively the same of the one shown in Fig. 5-(a), further validating our analysis.To complete the comparison between our classical analysis and the one performed in Ref. [50], we also analyze the scaling behavior of χ µ against µ and find that the fidelity susceptibility scales as χ µ ∼ 1/µ α in the range of values of µ explored, as shown in Fig. 5-(c).We compute the corresponding exponent α as a function of the energy density E and plot it in Fig. 5-(d): in the ergodic phase, we find that α is slightly greater than 1 (see inset in Fig. 5-(d)), resulting in a scaling approximately consistent with random matrix theory,2 while α exhibits a maximum of α ≃ 1.8 at E ≃ −0.4,close to the point where the thermalization time exceeds the simulation time window.Further insight is gained when investigating the profiles of the rescaled fidelities µχ µ against E, as done in the inset Fig. 5-(d): at high energy densities, the profiles collapse in a region which expands toward lower energies as we decrease µ.This collapse will in general break down at low energies, in particular where the maximum of χ µ is expected to occur.It is also worth mentioning that, from a deeper inspection collapse of the various profiles, we could in principle extract the Thouless time [84],defined as the typical relaxation time-scale τ th (E) of the correlation function at fixed energy density E, using the following prescription: first we define the cut-off µ E as the largest number such that the profiles of χ µ collapse for all energy densities greater than E and for all µ < µ E ; then the Thouless time can be easily obtained as τ th (E) ∼ µ −1 E [74].We expect that τ th (E) diverges as we the dynamics approaches the ergodicity breaking point from above.
As the cutoff µ plays the same role of the lowest level spacing for our analysis, these findings are consistent with those in Refs.[50,74], which show that the fidelity exhibits the strongest divergence with µ when the corresponding spin relaxation dynamics slows down.We therefore conclude that the fidelity susceptibility could be a good indicator of ergodicity even in classical systems, which warrants further investigation.

Conclusions
In this work we have investigated the unitary dynamics of the quantum p-spin glass spherical model in the limit of small ħ h, where the dynamics is effectively classical, and for the paradigmatic case of p = 3. Fixing the initial condition as an ensemble of narrow wavepackets centered on a fixed classical energy shell, at energy density E, we have investigated the chaotic dynamics of the model by the exponential divergence of the classical trajectories evolving from the same wave-packets.We have found that the corresponding exponent λ L is maximised when E is close to 0. We have found a similar behavior in a different chaos estimator, the Kolmogorov-Sinai entropy which also shows a pronounced maximum around E = 0.
To gain further insights into this result we have investigated the relaxation dynamics of the classical PSM at fixed energy density.We have shown that the spin correlation function displays a crossover in energy, from wide oscillations at high energies to aging low energies, consistently with the known results obtained at finite temperature.Interestingly we have found that around E = 0, where the chaos is maximized, the spin relaxation dynamics is the fastest.We give a physical interpretation of all our results in terms of the typical behaviour of the underlying trajectories, which either perform a uniform circular motion at asymptotically high energies or oscillate, at low energies, around a local minimum of the energy landscape.In between these two limits, we suggest that chaos emerges as the trajectories are scattered over the exponentially many saddles of the underlying landscape.Indeed a calculation of the number of stationary configurations shows that the complexity is also maximal at the same energy.Finally, we also gave a classical definition of the fidelity susceptibility.We found that the fidelity has a single maximum, as function of E, corresponding to the observed slowing down of the dynamics, a result reminiscent of the ones found in Refs.[50,74].
The results presented in this study hold true in the ħ h → 0 limit, where the TWA faithfully reproduces the dynamics for the initial condition defined in Eq. (11).Our findings can be potentially extended beyond the realm of small ħ h, by utilizing a Wigner state that reproduces the same fluctuations of a realistic quantum micro-canonical or canonical state.In this context, the quantum Lyapunov exponent λ L can be computed even at finite ħ h, derived from the exponential growth of a classical analog of the OTOC [28,33].A similar rationale applies to the correlation function and subsequently to the fidelity susceptibility, where the latter can be computed using Eq. ( 30) for both classical and quantum dynamics.Our analysis can also be extended to the transverse-field counterpart of the p-spin glass model, where chaos has been recently observed experimentally [32] and where the energy minima exhibit a more complicate structure [85].

A The truncated Wigner approximation for the p-spin spherical model
In this Appendix we show that, for the quantum spin glass defined in Eq. ( 1) of the main text in the Truncated Wigner Approximation (TWA) the dynamics is ruled by Eqs.( 6) and that such approximation becomes exact in the classical limit.We begin by rewriting the Weyl symbol for a generic operator evolving in the Heisenberg picture.As shown in Ref. [60], the Eq.(A.1) can be represented in a path integral form, suitable to study both the ħ h → 0 and the thermodynamic limit.Without reproducing the details of the calculation, here we just quote the final result: with initial conditions σ(0) = σ, ξ(0) = ξ and Π(0) = Π.The Weyl symbol of the Hamiltonian is defined as H W σ, Π = Π 2 /2M +V J (σ), and we inserted a term proportional to the Lagrange multiplier z(τ) in the action, to enforce the spherical constraint (see also Ref. [22]).
It is straightforward to see that TWA is equivalent to a (classical) expansion of the action in the integrand of Eq. (A.3): indeed we obtain, at leading order, that Integrating out the variables η(τ) and ξ(τ), we obtain a δ-function constraint on the trajectories σ(τ) and Π(τ): The Lagrange multiplier is obtained by imposing the constraint σ 2 = N , that by imposing that the total radial force appearing in the second of Eqs.(A.5) is the correct centripetal one: Replacing the expression of z(t) obtained in this way into Eqs.(A.5), we obtain exactly the Eqs.( 6) of the main text.Such expansion at leading order is exact in the limit ħ h → 0 where the dynamics is determined by the saddle point of the action appearing in the path-integral.

B Comparison with exponential growth of the OTOC
In this appendix, we show that, in the classical limit of ħ h → 0, the distance d J (t) defined in Eq. ( 12) of the main text diverges with the same Lyapunov exponent λ L of the corresponding out-of-time-order correlator (OTOC).

B.1 Classical chaos in dynamical systems
For a classical dynamical system, chaos is defined as the exponential divergence in time of the distance between orbits starting at nearby initial conditions.For example, we consider the system of Eqs. ( 6) of the main text and consider two solutions of such equations, a reference trajectory ỹ(y, t) evolving respectively form a generic initial condition y = (σ, Π), and a perturbed one ỹ(y + δy, t), evolving from y + δy for |δy| ≪ 1.Then the reference trajectory is said to be chaotic if the square distance ∆(t) = |ỹ(y + δy, t) − ỹ(y, t)| 2 grows exponentially in time or, in a more mathematical language, if the corresponding Lyapunov exponent [42] is positive.In practice, for ∆(0) → 0 we have that so that λ L is also detected by the exponential growth at long times of the matrix elements of the derivative matrix M(y, t) = ∂ ỹ(y, t)/∂ y.
We remark that, for classical systems, the Lyapunov exponent λ cl is defined from the exponential growth of the distance ∆(t), instead of ∆(t), so that we have λ L = 2λ cl .However, as explained in the next section, the definition in Eq. (B.1) is the classical limit of the quantum Lyapunov exponent detected by the out-of-time-order correlator (OTOC) and thus more suitable to make a comparison with the chaotic dynamics in a quantum system.

B.2 The exponential growth of the OTOC
For quantum systems, a generalization of the Lyapunov exponent can be defined from the exponential growth of the average square commutators [33] retrieved also in the corresponding OTOC [37] Such a generalization comes from the observation that, in the limit of small ħ h, the commutator [ σi (t), Πj (0)]/iħ h can be replaced by the corresponding Poisson parenthesis: where the quantum average 〈•〉 is replaced by a suitable average over the classical trajectories.Then, for ħ h → 0, F(t) is expected to grow with the same Lyapunov exponent λ L characterizing the underlying classical dynamical system.

B.3 The exponential growth of the distance
In this subsection we show that, for ħ h → 0, the distance d J (t) defined in Eq. ( 12) of the main text also grows exponentially in time, with exponent given by λ L .We begin by rewriting each of the average over a single wave-packet |α〉, appearing on the right-hand side of Eq. ( 12) as where we remind that ỹ(y, t) is a trajectory evolving from y = (σ, Π).The average 〈•〉 α on the right-hand side of Eq. (B.6) is performed over the coherent wave-packet For ħ h → 0, W α (y) becomes a δ-function of 0 width.To investigate the behaviour of Eq. (B.6) in this limit, it is useful to perform the change of variable y = y α + ħ hx.Then we have: Then, plugging Eqs.(B.8) into Eq.(B.6) and performing a little algebra, we obtain As explained in the first subsection of this Appendix, terms in the sum on the right-hand side of Eq. (B.9) growth exponentially in time, for a chaotic system.Equation (B.9) reproduces the dynamics of the derivative matrix only for a finite time-window: the noise produced by the fluctuations of order ħ h eventually lead to a saturation of the typical distance between trajectories evolving from a neighbourhood of y α .To obtain the Lyapunov exponent within such finite time-scale, we also average over the possible configurations of y α on the same energy shell, as explained in the main text and obtain which for small ħ h is expected to grow exponentially with exponent λ L , as confirmed also by Fig. 1 of the main text.

C Correlation function on a longer time-scale
In this Appendix, we discuss some results obtained from the dynamics of the p-spin model on larger time-scales.In particular, we integrate the dynamics according to the protocol described in Section 2 of the main text, and compute the corresponding correlation function C(t, t ′ ) in Eq. ( 16).Here the averages are performed on fewer realizations, with respect to the results presented in the main text.The results in Fig. 6 show that the correlation function seems to break time-translation invariance on energy scales between −0.25 and −0.42, even though the profiles of C(t w + τ, t w ), for various t w , are less smooth due to the noise induced by the fewer realization we took.The results shown here are anyway compatible with the ones presented in the main text.The results shown in Fig. 6 are anyway compatible with the ones presented in the main text and are also qualitatively similar with the plots describing the correlation function of a classical spin-glass in the non-ergodic phase and in the large-N limit, see for example Fig. (10) of Ref. [18].

D Details on the fidelity susceptibility
In this appendix, we prove the expression in Eq. ( 27) of the main text for the average fidelity susceptibility χ.We also discuss how the qualitative profile of the fidelity χ µ , defined in Eq. (30) and shown in Figg.5-(a) and (b) of the main text, is independent both of the time window [0, T ], over which we average the average correlation function in Eq. ( 28), and of the value of the cut-off µ, provided that the latter is sufficiently small.
We begin by recalling our definition for the fidelity susceptibility: We perturb the Hamiltonian ĤJ = 1 2M N i=1

Π2
i + V J ( σ) of the PSM with local magnetic fields, as and define the local susceptibilities of the as where |n(B)〉 is the n-th eigenstate of the perturbed Hamiltonian Ĥ(B) and |n〉 = |n(0)〉.By standard calculations of non-degenerate perturbation theory [86], we have where E n is the n-th unperturbed energy level of the Hamiltonian ĤJ .For the initial state ρ defined in Eq. ( 7) of the main text, we define the fidelity susceptibility χ as the weighted average performed over the sites, the eigenstates and the disorder configurations.Then, we observe that χ is connected to the Fourier transform of the time-averaged correlation function In particular, C(t w , t w + τ) and C av (τ) can be represented as follows: (D.6) The second line leads immediately to the Lehmann representation of the average correlation function, which reads as3 The latter is immediately related to the typical susceptibility χ in Eq.(D.4) via the expression where ω L is the average spacing of the unperturbed energy levels [50].
As stated in the main text, within the TWA framework we can not compute the fidelity directly from Eq. (D.8) for two reasons: first, our simulations are performed up to a finite maximum time t ma x , so that we can perform the average in Eq. ( 28) only on a finite time window [0, T ], with T < t ma x ; second, in the limit ħ h → 0 we do not have access to the spacing ω L and have a frequency cut-off set by ∆ω = 2π/t max .The second issue was already solved by using the regularized fidelity χ µ , from Eq. ( 30) of the main text, in place of χ.In Fig. 7 we also show that the profile of χ µ is qualitatively the same for a wide range of T between 0 and t ma x , so that the first issue is actually irrelevant for our results.

E Calculation of the complexity
In this appendix, we compute in detail the average number of stationary configurations of the dynamics induced by Eq. ( 6) of the main text.We start by rewriting the Eqs.(20), defining such stationary points: More precisely, we look for a solution of Eq. (E.1) lying on a microcanonical manifold defined by the equation Thus, as the kinetic energy vanishes (as Π i = 0 for every i), we can rewrite the equations for the stationary configurations in the following equivalent form: The system of Eqs.(E.3) can be further simplified by the following two observations.The first one is that, as we are interested in counting the number of the solutions of Eqs.(E.3), we can just focus on the number of solutions of the equations involving σ, as the trivial equations Π i = 0 do not bring any degeneracy.The second one is that the equations involving σ, written in the first three lines of Eqs.(E.3), are equivalent to the following reduced system of equations The equivalence can be seen by multiplying the first line of the system of Eqs.(E.4) by σ i and summing over i: by making use of the spherical constraint, one recovers the equation E = V J (σ)/N ; then, by substituting back E = V J (σ)/N in the first line of Eqs.(E.4), we obtain the first line of Eqs.(E.3).
To summarize, the number of stationary points of Eqs.(20) lying on a manifold at energy density E coincides with the number of solution of the equations for i = 1, . . ., N , lying on the N -sphere.In the spirit of Ref. [72], in what follows we will often write the indices for the p = 3 case, such that J i1...ip becomes J i jk .However, to give formulas that are valid even in the general case, we will write all the factors containing a term p for the generic p.Then, the average number of solutions of Eq. (E.5) on the N -sphere then reads: the overbar denoting the average over the disorder and the spherical constraint will be from now on hidden in the integration measure Dσ = δ( i σ 2 i − N ) N i=1 dσ i .The absolute value of the determinant appearing on the right-hand side is just a Jacobian factor appearing because the Eqs.(E.5) are written in an implicit form.
To compute N (E) in the large-N limit, we need two approximations.The first one consists in assuming that there is no correlation between the last two terms in the right-hand side of Eq. (E.6) [11,72], that is so that each the average of the δ-function and the one of the determinant can be computed independently from each other.We compute the average δ-function by using the exponential representation and averaging out the disorder after a proper symmetrization of the exponent.Taking into account the spherical constraint i σ 2 i = N and posing J = 1 for simplicity, we get To get rid of the term ( j µ j σ j ) 2 , we perform and extra Hubbard-Stratonovich transformation and perform straightforwardly all the remaining Gaussian integrals, posing again i σ 2 i = N in all our expressions.The final result is: up to a multiplicative constant which becomes irrelevant in the thermodynamic limit.The integration of the Jacobian factor is a bit more tricky.To perform it, we make our second approximation by assuming that the sign of the determinant, for any fixed configuration of the disorder, is given by the average number of negative eigenvalues of the corresponding Hessian matrix at energy density E, that we write as N k(E).In formulas, this is equivalent to: ) , (E.11) k(E) being the average fraction of negative eigenvalues.Once we got rid of the modulus, we rewrite the average of the determinant using a fermionic representation [87]: (E.12) Integrating out the disorder and using again an Hubbard-Stratonovich transformation to get rid of quartic fermionic terms (see Ref. [72] for more details about this step), we arrive at det where G(z) = z 2 p(p−1) + log(z − pE).Plugging everything together, we have that (up to an irrelevant prefactor) Thus, we are left with the computation of the integral along the imaginary axis Γ in the complex plane.In the thermodynamic limit N → ∞, this goal can be achieved by using the saddle-point method [88], which we briefly review in the following.First, we observe that, for generic z = x + i y (for x, y real numbers), the function G(z) = u(x, y) + i v(x, y) can be decomposed in its real and imaginary parts as x y + arctan y x − pE + πΘ(pE − x)sign( y) . (E.16) In summary, the saddle point method states that if we find a deformation γ of Γ in the context plane such that: 1. v(x, y) is constant over γ, 2. u(x, y) has a global maximum along γ at some point z = z 0 , 3. G(z) is analytic in the closed domain encompassed by the curves Γ and γ, we have that where the last asymptotic relation holds in the N → ∞ limit and is known as Laplace method [88].It is easy to see that the first condition is equivalent to state that γ is parallel to ∇u(x, y), as the relation ∇u•∇v = 0 holds for the holomorphic function G(z).Then if ∇u(x, y) vanishes along γ, it vanishes in the whole 2 plane, so that any maximum of u(x, y) along γ is a stationary point of u(x, y) in 2 .However, the choice of a suitable γ depends on the value of the energy density E, because G(z) has a branch-cut on the half-line {z = x|x < pE} (see the expression of v(x, y) in eq (E.16)) and the position in the complex plane of the stationary points of u(x, y), given by also depends on E. In particular, we can identify four relevant energy windows, which we treat separately: For E < E th , we can take the curve γ as the level curve v(x, y) = v(z + (E)) = 04 (shown in Fig. 8-( where we posed the phase k(E) = 0 in Eq. (E.14),as N (E) has to be a positive real number.The physical meaning of having k(E) = 0 is that, in this energy range, the integral in Eq. (E.6) is dominated by local minima of the potential V J (σ), where the Hessian is positive definite.For E th < E < 0, the only suitable deformation of Γ is γ = γ + ∪ γ − , where γ + and γ − are respectively the two level curves v(z) = v(z + (E)) and v(z) = v(z − (E)) of v(z).The two curves intersect respectively the points z + (E) and z − (E) in the complex plane (see Fig. 8-(b)), which in turn are maxima of u(x, y) along each of the two curves.Then, as u(z + ) = u(z − ) and v(z + ) = −v(z − ), the N → ∞ asymptotic value of I Γ is the sum of two contributions: e iN v(z + ) + e −iN v(z + ) 2 e N u(z + ) .(E.21) To give a physical interpretation to our result, let us first note that the function N k(E), defined in Eq. (E.11), is a non-negative integer.This is because N k(E) is the average number of Since N (E) is a positive real number, any phase obtained from the integral I Γ must compensate for the one coming from k(E).Therefore, with a bit of lack of rigor, we can conclude that for all physically meaningful values of k(E), the following equalities hold: ) is dominated by saddles having an extensive number of unstable direction: this is why in literature E th is referred to as the energy threshold where the local minima cease to dominate [10,72].For 0 < E < |E th |, the branch-cut crosses the integration path Γ and we can not use the saddle-point method straightforwardly.Instead, we first notice that it exist two curve in the complex plane, γ + and γ − , passing respectively through the saddle points z + (E) and z − (E) (both being maxima of u(z) along such curves) and ending up in a point z = x 1 (E) on the branch-cut.However, we observe that if we split Γ in two curves, Γ + and Γ − , defined in such a way that Finally, for E > |E th | both the z + (E) and z − (E) lie on the branch-cut, like in Fig. 8-(d).By some algebraic manipulations, one can show that u(x, y) has an absolute maximum at the point z − (E) along the level curves γ + and γ − of v(x, y) that intersect z − (E), while u(x, y) has a minimum at the point z + (E) along the level curves of v(x, y) that intersect z + (E).As done for the case 0 < E < |E th | we divide once again Γ in the curves Γ + and Γ − , both ending in z − (E) and deform each of them respectively along γ + and γ − (see Fig. 8-(d)).By observing that v(z) = ±π respectively along γ + and γ − , we conclude that k(E) = 1 in the energy range taken into exam, meaning that for E > |E th | the majority of stationary points of V J (σ) are local maxima.At the same time, the application of the Laplace method gives us: which is the result plotted in Figure 4 of the main text.

Figure 1 :
Figure1: (Left) Time-evolution of the log-average over the disorder of the distance d J (t) in Eq. (12), reported for few paradigmatic energy densities .The average over the initial condition is performed extracting 25 random configurations from each of the N s = 20 wave-packets, defined in Eq. (8) and obtained through the classical annealing algorithm.The dynamics is integrated up to a time t max = 80.The logaverage over the disorder is taken over N d = 30 configurations.(Right) Lyapunov exponent λ L , defined in Eq. (13), obtained through a linear fit of log d J (t) over a time window [t I , t F ], defined in such a way that, for each E, t I is beyond the oscillations at early times displayed by log d J (t) and t F is smaller than the saturation time-scale.

Figure 2 :
Figure 2: (Color online) (Left) Time-evolution of the quantity S KS defined in Eqs.(14) and (15), reported for few paradigmatic energy densities.The average over the initial condition is performed extracting 30 random configurations from each of the N s = 30 wave-packets, defined in Eq. (8) and obtained through the classical annealing algorithm.The dynamics is integrated up to a time t max = 80.The log-average over the disorder is taken over N d = 96 configurations.(Right) Kolmogorov-Sinai entropy Λ KS , defined in Eq. (15), obtained through a linear fit of S KS over a time window [t I , 10], where t I is beyond the scale of oscillations at early time displayed by S KS , for each E.

Figure 3 :
Figure 3: (a) Time-dependence of the correlation function, at fixed t w = 40.(b)Asymptotic value q 1 of the correlation functions C(t w , t w + τ) reported in panel (a), obtained through the time-average of the latter over τ ∈[39,40].(c) Comparison between the profiles of C(t w , t w = τ) obtained fixing different values of t w .Each panel correspond to a different energy density E. (d) Time-evolution of three typical orbits, whose initial condition are obtained extracting one configuration from the distribution in Eq.(11), at different energy densities E and for the same configuration of the disorder.The orbits evolve in a 200-dimensional phase space and are projected over the two axes defined by the initial spin configuration σ(0) and the initial momentum Π(0).

Figure 5 :
Figure 5: (Color online) Fidelity susceptibility χ µ from Eq. (30).(a) χ µ is shown as a function of E and fixing several values of the cut-off µ.The data are obtained from a dynamics up to time t ma x = 80, with the same parameters described in Fig. 3.The average time window for the correlation function in Eq. (28) is set to [0, T ], with T = t ma x = 40.(b) Same plot of panel (a), for a dynamics integrated up to time t ma x = 320.Here the average over the initial condition is performed over 5 random configurations extracted from each the N s = 10 wave-packets constructed by the classical annealing algorithm.We average over N d = 42 disorder configurations.(c) Same data of panel (b).χ µ is shown as a function of the cutoff µ, on a log-log scale, for some fixed values of the energy density E. (d) Exponent α(E) obtained by a linear fit of log χ µ against − log µ, at several fixed values of the energy density E. For each E, the data used for the fit are the ones from panel (b).

Figure 6 :
Figure 6: Several plots of the correlation function C(t w + τ, t w ), for different fixed values of the energy density E. For all the panels, the data are obtained from a dynamics up to time t ma x = 320, with the same simulation parameters described in Fig. 5-(b) of the main text.The data are presented on a log-log scale.

Figure 7 :
Figure 7: Fidelity susceptibility χ µ from Eq. (30) of the main text, shown as a function of E for a fixed cutoff µ = 0.07.The average time window for the correlation function in Eq. (28) is set to [0, T ], for several values of T .The data are obtained from a dynamics up to time t ma x = 320, with the same parameters described in Fig.5-(b) of the main text.

Figure 8 :
Figure 8: (Color online) Color map in the complex plane representing the values assumed by the function v(x, y) defined in the appendix.Each panel corresponds to a different, fixed value of the energy density E, each representative of one of the four, qualitatively different cases studied in the appendix.Here we plot the results for p = 3.The white thick line correponds to the branch-cut of v(z) = v(x, y) in the complex plane, where z = x + i y, while the black lines represent the integration contour Γ or its deformation Γ + ∪ Γ − .The blue and the red lines correspond to the level curve of v(z), passing respectively through z − (E) and z + (E).(a) E < E th .(b)E th < E < 0.(c) 0 < E < |E th |.(d) E > |E th |.