Low-frequency vibrational spectrum of mean-field disordered systems
Eran Bouchbinder, Edan Lerner, Corrado Rainone, Pierfrancesco Urbani, Francesco Zamponi
LLow-frequency vibrational spectrum of mean-field disordered systems
Eran Bouchbinder, Edan Lerner, Corrado Rainone, Pierfrancesco Urbani, and Francesco Zamponi Chemical and Biological Physics Department, Weizmann Institute of Science, Rehovot 7610001, Israel Institute of Theoretical Physics, University of Amsterdam,Science Park 904, 1098 XH Amsterdam, the Netherlands Universit´e Paris-Saclay, CNRS, CEA, Institut de physique th´eorique, 91191, Gif-sur-Yvette, France Laboratoire de Physique de l’Ecole Normale Sup´erieure, ENS, Universit´e PSL,CNRS, Sorbonne Universit´e, Universit´e de Paris, F-75005 Paris, France
We study a recently introduced and exactly solvable mean-field model for the density of vibra-tional states D ( ω ) of a structurally disordered system. The model is formulated as a collection ofdisordered anharmonic oscillators, with random stiffness κ drawn from a distribution p ( κ ), subjectedto a constant field h and interacting bilinearly with a coupling of strength J . We investigate thevibrational properties of its ground state at zero temperature. When p ( κ ) is gapped, the emergent D ( ω ) is also gapped, for small J . Upon increasing J , the gap vanishes on a critical line in the ( h, J )phase diagram, whereupon replica symmetry is broken. At small h , the form of this pseudogapis quadratic, D ( ω ) ∼ ω , and its modes are delocalized, as expected from previously investigatedmean-field spin glass models. However, we determine that for large enough h , a quartic pseudogap D ( ω ) ∼ ω , populated by localized modes, emerges, the two regimes being separated by a specialpoint on the critical line. We thus uncover that mean-field disordered systems can generically displayboth a quadratic-delocalized and a quartic-localized spectrum at the glass transition. Introduction —
The vibrational spectrum of struc-tural glasses displays a series of universal features in dif-ferent frequency ranges, which are responsible for impor-tant material properties, such as wave attenuation, heattransport, and plasticity [1–3]. Motivated by these ob-servations, several authors have constructed simple mod-els of the non-phononic vibrational density of states ofstructurally disordered systems, D ( ω ) [4–20]. Mean-fieldmodels typically display a quadratic spectrum, D ( ω ) ∼ ω [10, 12], of delocalized and featureless modes [21]; thisdelocalization is inherently different from the one associ-ated with phononic excitations in solids (which are ab-sent in the mean-field limit), and is a manifestation of themarginal stability associated to replica symmetry break-ing [22, 23].On the contrary, numerical simulations of model glassformers in finite dimension have revealed that non-phononic excitations in those systems are quasi-localized in nature, that they emerge from self-organized glassyfrustration, and that they follow a seemingly universal quartic law D ( ω ) ∼ ω [24–29]. Given these discrep-ancies with the mean-field scenario detailed above, andthe localized nature of these excitations, the naive ex-pectation is that the modes that populate the quarticlaw would disappear in the mean-field limit, and thatmean-field models are therefore unable to tell much aboutthe physics responsible for the ω spectrum of structuralglasses [30–32].Nearly two decades ago, Gurevich, Parshin andSchober (GPS) proposed a three-dimensional latticemodel [6] for this glassy density of states, formulated interms of interacting anharmonic oscillators, with a cou-pling strength that decays with distance as ∼ r − , where r is the distance between the oscillators. GPS showednumerically that D ( ω ) ∼ ω emerges in that model, andproposed a phenomenological theory [6, 8], which has later been investigated by other authors [20]. A similarmean-field model was studied by K¨uhn and Horstmann(KH) [4], who however did not investigate the vibrationalspectrum; as we will show below, their model’s spectrumfollows D ( ω ) ∼ ω .In this Letter, we study a recently introducedmodel [33], which corresponds to both the infinite-dimensional, mean-field version of the GPS model, andto a generalization of the KH model. Following Ref. [33],we hereafter refer to it as the KHGPS model . The modelis formulated as a collection of N interacting anharmonicoscillators, each represented by a generalized coordinate x i and stiffness κ i [34]; the model’s Hamiltonian reads H ≡ (cid:88) i 0, in such away that all the oscillators have a single minimum at J = 0. An external constant “magnetic” field h is addedin order to break the spurious x i → − x i symmetry thathas no counterpart in amorphous solids [35, 36]. Themodel is related to a soft-spin version of the Sherrington-Kirkpatrick model [37].In what follows we describe the exact solution of theKHGPS model using the replica method [4, 22]. We con-struct the phase diagram of the model, in the plane ofthe applied magnetic field h and the coupling strength J .We rigorously show that, for κ m > 0, the model’s spec-trum is gapped at small enough coupling, in the replicasymmetric phase where the energy landscape is convex. a r X i v : . [ c ond - m a t . d i s - nn ] D ec Upon increasing the coupling strength, a phase transi-tion is encountered, whereupon replica symmetry is bro-ken, the energy landscape becomes rough, and the gapin the spectrum closes. On this critical line, the spec-trum behaves as D ( ω ) ∼ ω at small h , a typical mean-field scenario. Conversely, for large h , the spectrum be-haves as D ( ω ) ∼ ω and its modes are partially local-ized. The two regimes are separated by a special pointon the critical line, whose location is determined. Allin all, our work demonstrates that disordered mean-fieldmodels can display a quartic density of states of localizedmodes in certain regions of their phase diagram, includ-ing critical lines whereupon replica symmetry is broken.Related results have been reported in Ref. [38] for theXY model defined on a random graph, which is how-ever much more difficult to analyze. This result opensnew perspectives for the microscopic understanding ofthe universal D ( ω ) ∼ ω law in finite-dimensional glassysystems. Furthermore, it shows that replica symmetrybreaking (RSB) phase transitions can present profoundlydifferent characteristics from the marginal stability sce-nario usually associated to it, even at a mean-field level. Vibrational spectrum — The Hessian M ij ≡ ∂ H/∂x i ∂x j corresponding to H takes the form M ij = J ij + δ ij (cid:18) κ i + 12 x i (cid:19) ≡ J ij + δ ij a i , (2)which is the sum of a member ( J ij ) of the Gaussian Or-thogonal ensemble (GOE) of random matrices [21] andof a diagonal matrix A ij = a i δ ij , with diagonal elements a i ≡ κ i + x i / 2. Assuming that there is no statisticalcorrelation between these two matrices, calculating thespectrum of their sum becomes a standard problem inrandom matrix theory [21], which only requires knowl-edge of the statistics of the diagonal part. Past effortsin calculating typical ground-state spectra of mean-fielddisordered systems [10] indicate that this assumption isvalid, hence we adopt it here and proceed as follows.Assuming that the statistics p ( a ) of the diagonal el-ements is known and has support in [ a m , a M ] (we willcompute it below), one can compute the density of eigen-values ρ M ( λ ) of M by defining the resolvent [39], g M ( z ) ≡ (cid:90) d λ ρ M ( λ ) z − λ , (3)which implies ρ M ( λ ) = 1 π lim η → + Im [ g M ( λ − iη )] . (4)The resolvent of M is then implicitly expressed in termsof the spectrum of the diagonal part, ρ A ( a ) = p ( a ), as [39] g M ( z ) = (cid:90) a M a m d a p ( a ) 1 z − a − J g M ( z ) . (5)This equation requires a numerical solution, but theposition and shape of the lower edge of the spec-trum can be worked out analytically. Let us define g ( z ) ≡ g M ( z ) − z/J ; one can then recast Eq. (5) as z = − J (cid:90) a M a m d a p ( a ) (cid:20) g + 1 a + J g (cid:21) ≡ F ( g ) . (6)We next argue as follows: for λ outside of the supportof the spectrum ρ M ( λ ) and in the limit η → 0, Eq. (4)implies that g ( z ) cannot have an imaginary part. Wetherefore expect a band of values of the function F ( g ) tobe forbidden for real g , and to correspond to the supportof the spectrum. Let us then consider which values F ( g )can attain for real g . This function is obviously not de-fined to the left of − a M /J or to the right of − a m /J ,and intuitively, we expect the branch for g > − a m /J to be the one controlling the lower edge; therefore, thisbranch needs to be bounded from above. There are thenonly two possibilities: (i) The function F ( g ) has a maxi-mum g m for g > − a m /J , meaning that F (cid:48) ( g m ) = 0. Thecorresponding value of F , λ m = F ( g m ), is then the loweredge of the spectrum. In this case, the support of the di-agonal elements has no influence on the lower edge, andthe GOE part of M dominates: close to the edge theeigenvectors are delocalized and ρ M ( λ ) ∝ ( λ − λ m ) / [40].We dub this a GOE-like spectrum. (ii) The function F ( g )has no maximum for g > − a m /J . In this case, the valueof g that corresponds to the edge must be g m = − a m /J ,and the lower edge itself is λ m = F (cid:16) − a m J (cid:17) = a m − J (cid:90) a M a m d a p ( a ) 1 a − a m . (7)The edge is then determined by the support of p ( a ), anddominated by the diagonal part of M , and the eigenvec-tors near the edge are partially localized [40]. We dub thisa DIAG-like spectrum. Furthermore, if p ( a ) ∼ ( a − a m ) ν near its lower edge, one can show by an expansion near g m that ρ M ( λ ) ∼ ( λ − λ m ) ν .The value of coupling that separates the two regimes issuch that F (cid:48) ( g = g m ) = 0, which gives the self-consistentequation:Λ = 1 − J (cid:90) a M a m d a p ( a ) 1( a − a m ) = 0 , (8)such that Λ > < p ( a ) and its edges depend on J , this equation defines thecritical value of coupling at which the spectrum changesshape only implicitly. Replica method — We next aim at determining thestatistics p ( a ) of the diagonal elements a i appearingin Eq. (2), which depend on the oscillator positions x i in the ground state of the model, and can be de-termined by solving its thermodynamics in the zero-temperature ( T → 0) limit. We do so by employing thereplica method [22]; we assume that the ground state isunique, which corresponds to a replica-symmetric (RS)ansatz. This picture is expected to be justified as longas the coupling strength J is below a critical threshold J c ( h ), which we self-consistently determine below.We briefly delineate the steps in obtaining the RS so-lution of the model, leaving the details to Appendix B.The mean-field nature of the model allows one to writeits solution as a problem of decoupled oscillators in an ef-fective, self-consistent random external potential, whichin the T → v eff ( x ) ≡ x 4! + m x − ( f + h ) x , m = κ − J χ , (9)where f is a Gaussian random force of zero mean andvariance J ˜ q , and the new parameters χ and ˜ q emergefrom the correlations between different replicas gener-ated by the disorder, and have to be determined self-consistently [22].Depending on the value of the coefficients, the effectivepotential can be either an asymmetric single well (SW)or double well (DW), with two minima separated by anenergy barrier. In particular, if the effective stiffness m is negative, there is always some value of the field f forwhich the potential is a DW. We thus conclude that if m m = κ m − J χ < 0, DWs appear with finite probabil-ity, and we show in Appendix B 4 that in this case the RSsolution is always unstable towards RSB. Consequently,we now restrict ourselves to the case m m ≥ 0, which isrealized at small enough J if κ m > 0, and we discuss theRS phase of the model.Under the assumption m m ≥ 0, we show in the Ap-pendix that the parameters ˜ q and χ are self-consistentlydetermined through the equations χ = (cid:28) v (cid:48)(cid:48) eff ( x ∗ ( f, m )) (cid:29) m,f , ˜ q = (cid:104) ( x ∗ ( f, m )) (cid:105) m,f , (10)where x ∗ ( f, m ) denotes the point of absolute min-imum of the effective potential, and the aver-age is taken over the random effective stiffnesses m ∼ U ( m m = κ m − J χ, m M = κ M − J χ ) and random fields f ∼ N (0 , J ˜ q ). The self-consistency of this picture istested by verifying the positivity of the replicon eigen-value λ R of the Hessian matrix of the replica action [22,23]. The definition of the replica action and the compu-tation of the replicon can be found in Appendix B 3. Thefinal result reads λ R = 1 − J (cid:28) v (cid:48)(cid:48) eff ( x ∗ ( f, m )) (cid:29) m,f . (11)Recalling that a = κ + x / v (cid:48)(cid:48) eff ( x ) + J χ , cf. Eq. (2),where x has to be evaluated in x ∗ ( f, m ), we can express χ and λ R as χ = (cid:90) a M a m d a p ( a ) 1 a − J χ ,λ R = 1 − J (cid:90) a M a m d a p ( a ) 1( a − J χ ) , (12)where λ R differs from Λ in Eq. (8) by the replacementof a m by J χ in the denominator of the integrand. Notethat a = m + J χ + x / 2, hence under the assumptionthat m m ≥ 0, we have a m ≥ J χ , which imples λ R ≥ Λ. FIG. 1. Phase diagram of the model for κ m = 0 . κ M =1 in the ( h, J ) plane. The line J c ( h ) separates the convex-landscape RS phase from the rough-landscape RSB phase.Along the dotted line, D ( ω ) ∝ ω at the transition, whilst D ( ω ) ∝ ω on the solid line. Phase diagram — We are now in a position to deter-mine the phase diagram of the model in the ( h, J ) plane.We fix here κ m = 0 . κ M = 1, but the qualitativepicture is independent of this choice as long as κ m > J = 0, the curvatures v (cid:48)(cid:48) eff ( x ∗ ) are finite and conse-quently χ is finite. Hence, at J = 0 we have m = κ andconsequently for small enough J , the condition m m > m m > 0, we have a m > J χ andthe integral that appears in Eq. (12) is finite, leading to λ R ≈ J ≈ 0. We thus conclude that the RS phase isstable at small J . While it is easy to show that the spec-trum is always gapped in this phase, the sign of Λ, andthus the shape of the spectrum near its edge, depends onthe behavior of p ( a ) near a m and on the values of a m and a M . For example, the KH model studied in Ref. [4] has κ m = κ M = 1 and in that case a m = a M leading formallyto Λ = −∞ , in such a way that the spectrum is alwaysGOE-like. With our choice of p ( κ ), instead, the integralin Eq. (8) is finite and the spectrum is always DIAG-likeat low enough J .The RS phase can then become unstable in two ways:(i) The replicon can vanish, while m m remains strictlypositive. In this case, at the transition point we haveΛ ≤ λ R = 0, hence the spectrum is GOE-like. Fora GOE-like spectrum to be gapless, the two equations λ m = F ( g m )=0 and F (cid:48) ( g m )=0 must hold, which is equiv-alent to Eqs. (12) with λ R =0 and g m = − χ . Hence, thespectrum is gapless at the critical point and ρ M ( λ ) ∼ λ / , which is equivalent to D ( ω ) ∼ ω . Just above thecritical point, the replicon becomes negative. This is astandard RSB transition, observed in several spin glassmodels. (ii) The lower bound of the effective stiffnesscan vanish, m m = 0, while the replicon is still positive, λ R > 0. When m m < 0, there is a finite probabilityof having DWs in the ensemble of effective potentials,and we show in Appendix B 4 that this formally implies λ R = −∞ . Hence, the replicon jumps discontinuously tominus infinity beyond this transition. At the transitionpoint, m m = 0 implies (see Appendix C 2 for details)that a m = J χ , which implies that Λ = λ R > p ( a ) ∼ ( a − J χ ) / = ⇒ p (˜ a ) ∼ ˜ a / , (13)where ˜ a = v (cid:48)(cid:48) eff ( x ∗ ) = a − J χ is the curvature of theeffective potential at its minimum. Note that Eq. (10)then gives χ = (cid:104) / ˜ a (cid:105) and from Eq. (7) it follows that λ m = J [ χ − (cid:104) / ˜ a (cid:105) ] = 0. We conclude that the spec-trum is gapless and DIAG-like, i.e. ρ M ( λ ) ∼ λ / , orequivalently D ( ω ) ∼ ω .The phase diagram obtained by solving numericallyEqs. (10) is reported in Fig. 1. We observe that the glasstransition line J c ( h ) falls into case (i) for small h , andinto case (ii) for large h . The two lines are separatedby a special point at which m m = 0 and λ R = 0 simul-taneously. We also verify these predictions numerically,by directly calculating the spectrum D ( ω ) of the Hessianin the minima of the Hamiltonian in Eq. (1), obtainedby means of a gradient descent algorithm [41]. Thesenumerical results, which confirm our theoretical predic-tions, are reported in Fig. 2. We note that when κ m isreduced and approaches zero, the line J c ( h ) moves to-wards the left, i.e. towards smaller values of J , and the ω region increases; when κ m = 0, the model is in theRSB phase at all J . This regime was studied numericallyand through a scaling theory in Ref. [33]. Discussion — We studied a mean-field model of in-teracting disordered anharmonic oscillators [33] having,in absence of coupling, a gapped spectrum. We showedthat at small coupling the spectrum remains gapped [42],and that at the glass transition point it can display ei-ther the universal D ( ω ) ∼ ω localized spectra observedin finite-dimensional computer glass models [11, 25–29]and in the random graph XY model [38], or the stan-dard D ( ω ) ∼ ω observed in most mean-field spin glassmodels and jammed sphere packings [9, 10, 12, 30]. Theimmediate implication of our results is that systems at aRSB transition, and possibly even deep within the RSBphase, can in fact exhibit localized excitations, even atthe mean-field level. The class of models to which theKHGPS model studied here belongs is expected to berather broad — according to existing evidence [5, 6, 20]— and largely robust to changes in these models’ input.We note that an effective potential in the form of aquartic polynomial, Eq. (9), naturally emerges from ourtheory. This effective potential, which resembles the SoftPotential Model framework [3, 43, 44] that also predictsa ω nonphononic spectrum under some nontrivial as-sumptions (spelled out, e.g., in [5]). In light of our re-sults, the Soft Potential Model can be viewed as an effec-tive description of the collective, many-body statistical-mechanics of the KHGPS model.Moreover, we note that the zero temperature limit of -2 -1 -6 -5 -4 -3 -2 -1 -2 -4 -3 -2 -1 -6 -5 -4 -3 -2 -1 -1 -4 -2 FIG. 2. Numerical results for the vibrational spectrumof the model, at two selected values of h (top: h = 0, J c ( h ) = 0 . h = 0 . J c ( h ) = 0 . the spin glass susceptibility behaves very differently onthe two parts of the critical line, being divergent whenthe spectrum at the transition is ω , and finite when thespectrum is ω [45]. Finally, we stress that our resultsapply upon approaching the transition at strictly zerotemperature, and therefore it is important to investigatethe model’s behavior at finite temperature.In this work, we limited ourselves to the investigationof the RS phase of the model with κ m > 0, up to thecritical line whereupon replica symmetry is broken and aglassy phase appears. A natural direction for future re-search is to investigate the vibrational spectrum deep inthe glass phase. One might expect, by continuity argu-ments, that the quartic spectrum extends into the glassphase, hence being valid in a finite region of the phase di-agram. This point of view seems to be supported by thenumerical results of Ref. [33], but whether this intuitionis correct can only be confirmed by an investigation ofthe RSB equations of the model. The gradient descentdynamics might also display interesting features in theglass phase [46, 47] and, if minima reached by quenchingdynamics retain the properties of the model at the tran-sition, one could expect different (or even the absence of)aging dynamics. Acknowledgements .—We benefited from discussionswith Giulio Biroli, Jean-Philippe Bouchaud, GustavoD¨uring, Eric De Giuli, and Guilhem Semerjian. Thisproject has received funding from the European ResearchCouncil (ERC) under the European Union’s Horizon2020 research and innovation programme (grant agree- ment n ° [1] G. Ruocco and F. Sette, Journal of Physics: CondensedMatter , 9141 (2001).[2] T. Nakayama, Reports on Progress in Physics , 1195(2002).[3] U. Buchenau, Y. M. Galperin, V. L. Gurevich, and H. R.Schober, Phys. Rev. B , 5039 (1991).[4] R. K¨uhn and U. Horstmann, Phys. Rev. Lett. , 4067(1997).[5] V. Gurarie and J. T. Chalker, Phys. Rev. B , 134207(2003).[6] V. L. Gurevich, D. A. Parshin, and H. R. Schober, Phys.Rev. B , 094203 (2003).[7] T. S. Grigera, V. Mart´ın-Mayor, G. Parisi, and P. Ver-rocchio, Nature , 289 (2003).[8] D. A. Parshin, H. R. Schober, and V. L. Gurevich, Phys.Rev. B , 064206 (2007).[9] E. DeGiuli, A. Laversanne-Finot, G. During, E. Lerner,and M. Wyart, Soft Matter , 5628 (2014).[10] S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Proc.Natl. Acad. Sci. U.S.A. , 14539 (2015).[11] M. Baity-Jesi, V. Mart´ın-Mayor, G. Parisi, and S. Perez-Gaviro, Phys. Rev. Lett. , 267205 (2015).[12] A. Sharma, J. Yeo, and M. A. Moore, Phys. Rev. E ,052143 (2016).[13] F. P. C. Benetti, G. Parisi, F. Pietracaprina, and G. Si-curo, Phys. Rev. E , 062157 (2018).[14] Y. V. Fyodorov and P. Le Doussal, Journal of Physics A:Mathematical and Theoretical , 474002 (2018).[15] E. Stanifer, P. K. Morse, A. A. Middleton, and M. L.Manning, Phys. Rev. E , 042908 (2018).[16] H. Ikeda, Phys. Rev. E , 050901 (2019).[17] M. Baggioli and A. Zaccone, Phys. Rev. Lett. ,145501 (2019).[18] M. Shimada, H. Mizuno, and A. Ikeda, Soft Matter ,7279 (2020).[19] Y. V. Fyodorov and P. Le Doussal, Journal of StatisticalPhysics , 176 (2020).[20] P. Das, H. G. E. Hentschel, E. Lerner, and I. Procaccia,Phys. Rev. B , 014202 (2020).[21] G. Livan, M. Novaes, and P. Vivo, Introduction to ran-dom matrices: theory and practice (Springer, 2018).[22] M. Mezard, G. Parisi, and M. A. Virasoro, Spin glasstheory and beyond (World Scientific, Singapore, 1987).[23] G. Parisi, P. Urbani, and F. Zamponi, Theory of sim-ple glasses: exact solutions in infinite dimensions (Cam-bridge University Press, 2020).[24] B. B. Laird and H. R. Schober, Phys. Rev. Lett. , 636(1991).[25] E. Lerner, G. D¨uring, and E. Bouchbinder, Phys. Rev.Lett. , 035501 (2016). [26] H. Mizuno, H. Shiba, and A. Ikeda, Proc. Natl. Acad.Sci. U.S.A. , E9767 (2017).[27] G. Kapteijns, E. Bouchbinder, and E. Lerner, Phys. Rev.Lett. , 055501 (2018).[28] L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel,and E. Flenner, Nat. Comm. , 26 (2019).[29] D. Richard, K. Gonz´alez-L´opez, G. Kapteijns, R. Pater,T. Vaknin, E. Bouchbinder, and E. Lerner, Phys. Rev.Lett. , 085502 (2020).[30] P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet, andF. Zamponi, Phys. Rev. Lett. , 045503 (2016).[31] H. Ikeda and M. Shimada, arXiv:2009.12060 (2019).[32] M. Shimada, H. Mizuno, L. Berthier, and A. Ikeda, Phys.Rev. E , 052906 (2020).[33] C. Rainone, P. Urbani, F. Zamponi, E. Lerner, andE. Bouchbinder, arXiv:2010.11180 (2020).[34] For our abstract mathematical model, all quantities areassumed to be dimensionless.[35] P. Urbani and G. Biroli, Phys. Rev. B , 100202 (2015).[36] S. Albert, G. Biroli, F. Ladieu, R. Tourbot, and P. Ur-bani, arXiv:2010.03294 (2020).[37] H. Sompolinsky and A. Zippelius, Phys. Rev. B , 6860(1982).[38] C. Lupo, arXiv:1706.08899 (2017).[39] J. Bun, J.-P. Bouchaud, and M. Potters, Phys. Rep. ,1 (2017).[40] J. O. Lee and K. Schnelli, Probab. Theory Relat. Fields , 165 (2016).[41] We use the BFGS algorithm implemented in the GSLlibrary [48].[42] W. Ji, T. W. J. de Geus, M. Popovi´c, E. Agoritsas, andM. Wyart, arXiv:1912.10537 (2020).[43] U. Buchenau, Y. M. Galperin, V. L. Gurevich, D. A.Parshin, M. A. Ramos, and H. R. Schober, Phys. Rev.B , 2798 (1992).[44] V. Gurarie and J. T. Chalker, Physical Review B ,134207 (2003).[45] The spin glass susceptibility is defined by χ SG = (cid:80) ij ( (cid:104) x i x j (cid:105) − (cid:104) x i (cid:105)(cid:104) x j (cid:105) ) /N , and in the zero temperaturelimit it becomes χ SG = N Tr M − .[46] L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. ,173 (1993).[47] G. Folena, S. Franz, and F. Ricci-Tersenghi, PhysicalReview X , 031045 (2020).[48] M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman,P. Alken, M. Booth, F. Rossi, and R. Ulerich, GNUscientific library (Citeseer, 2002).[49] https://en.wikipedia.org/wiki/Cubic_equation . Appendix A: The spectrum The Hessian matrix, evaluated in a minimum x ∗ i of the Hamiltonian, is the sum of a GOE matrix J ij and of adiagonal matrix A ij = a i δ ij with diagonal elements a i = κ i + ( x ∗ i ) . The random variable a i is distributed accordingto p ( a ) in the interval [ a m , a M ]. In the following, we assume for simplicity that a i and J ij are uncorrelated; it can beproven both analytically and numerically that this assumption is correct [10]. 1. Resolvent equation We want to calculate the density of eigenvalues ρ M ( λ ) of the matrix M . This can be defined in terms of theresolvent (or rather, the trace of the resolvent in the thermodynamic limit) g M ( z ) ≡ lim N →∞ N Tr( z − M ) − = (cid:90) d λ ρ M ( λ ) z − λ , (A1)which implies ρ M ( λ ) = 1 π lim η → + Im g M ( λ − iη ) . (A2)The resolvent of M can be obtained in terms of the resolvent of the diagonal matrix A via the fixed-point equation g M ( z ) = g A ( z − J g M ( z )) , (A3)where we used the fact that J ij is a GOE matrix [39]. The resolvent of a diagonal matrix is trivial because ρ A ( λ ) = p ( a ),hence the fixed-point equation is explicitly written as g A ( z ) = (cid:90) d a p ( a ) 1 z − a ⇒ g M ( z ) = (cid:90) d a p ( a ) 1 z − a − J g M ( z ) . (A4) 2. Location of the spectrum edge To investigate the low-frequency tail of the spectrum we start from Eq. (A4), and we define g ( z ) = g M ( z ) − z/J so we get z = − J (cid:90) d a p ( a ) (cid:20) g + 1 a + J g (cid:21) ≡ F ( g ) . (A5)Outside the support of the spectrum, for z = λ − iη with η → g ( z ) needs to be real. We expect that there is aband of values of F ( g ) that are forbidden for real g , which correspond to the support of the spectrum. So, we studythe function F ( g ) for real g . We have F (cid:48) ( g ) = − J (cid:90) d a p ( a ) (cid:20) − J ( a + J g ) (cid:21) , F (cid:48)(cid:48) ( g ) = − J (cid:90) d a p ( a ) (cid:20) J ( a + J g ) (cid:21) . (A6)Note that if p ( a ) has support in [ a m , a M ], then F ( g ) is only defined for g / ∈ [ − a M /J , − a m /J ] on the real axis.There are two possibilities for the spectrum:GOE-like– Suppose that the function F ( g ) has a minimum for g M < − a M /J and a maximum for g m > − a m /J . In thiscase, if g m , g M are the solutions of F (cid:48) ( g ) = 0, then λ m = F ( g m ) and λ M = F ( g M ) are the edges of the spectrum.In the vicinity of the edges we can expand, e.g. for z = λ M − ε and g = g M + δg , and we get λ M − ε ∼ F ( g M + δg ) ∼ F ( g M ) + 12 F (cid:48)(cid:48) ( g M ) δg + · · · ⇒ ε = − F (cid:48)(cid:48) ( g M ) δg + · · · ⇒ δg = (cid:115) − λ M − z F (cid:48)(cid:48) ( g M ) . (A7)Clearly if F (cid:48)(cid:48) ( g M ) (cid:54) = 0 we get ρ ( λ ) ∼ (cid:112) λ M − λ in the vicinity of λ M . The same happens for the lower edge.DIAG-like– It can happen however that F ( g ) has no maximum for any g > − a m /J . In this case, the value of g thatcorresponds to the edge is g m = − a m /J , and the location of the edge is λ m = F ( − a m /J ) = a m − J (cid:90) a M a m d a p ( a ) (cid:20) a − a m (cid:21) . (A8)Critical J – The value of coupling that separates the GOE and DIAG regimes is such that F (cid:48) ( g = g m ) = 0, which gives theself-consistent equation: Λ = 1 − J (cid:90) a M a m d a p ( a ) 1( a − a m ) = 0 . (A9)Note that because p ( a ) and its edges depend on J , this equation defines the critical value at which the spectrumchanges shape only implicitly. The case Λ > F (cid:48) ( g = g m ) < 0, hence to a DIAG-like spectrum,while the case Λ < F (cid:48) ( g = g m ) > 0, hence to a GOE-like spectrum. 3. Shape of the edge and prefactors We now focus on the DIAG-like spectrum and we study in more details the behavior close to the edge. The analysisdepends on the details of p ( a ), so we will assume the power-law form p ( a ) ∼ A d ( a − a m ) / (A10)close to the lower edge. The analysis is performed similarly for other values of the exponent ν (cid:54) = 3 / F ( g ) is in g m = − a m /J and we want to expand F ( g ) around it. From Eq. (A6)we observe that F (cid:48) ( g m ) = − J Λ is finite, while F (cid:48)(cid:48) ( g m ) is divergent, which suggests a non-analytic behavior for F ( g )with exponent 3 / g m , as we now show. We define δg = g − g m and δ F ( δg ) = F ( g m + δg ) − F ( g m ) − F (cid:48) ( g m ) δg = − J δg (cid:90) d a p ( a ) 1( a − a m ) ( a − a m + J δg ) . (A11)For small δg we have, defining q ( a ) = p ( a ) / ( a − a m ) / → A d for a → a m , and changing variable to x = ( a − a m ) /δg , B = lim δg → (cid:112) δg (cid:90) d a q ( a ) √ a − a m ( a − a m + J δg ) = lim δg → (cid:90) d x q ( a m + xδg ) √ x ( x + J ) = A d πJ . (A12)Collecting all together these results, we have for small δg : δz = z − λ m = F ( g ) − F ( g m ) ∼ − J Λ δg − J π A d δg / + · · · . (A13)Inverting this relation we obtain δg ( z ) = − J Λ δz − π A d Λ / ( − δz ) / + · · · . (A14)If we choose δz = δλ to be real and positive, we getIm g M ( z ) = Im g ( z ) = π A d Λ / δλ / , ⇒ A g = A d Λ / . (A15)Similar results are obtained for other values of ν . 4. Summary So far, we have obtained the following results, for a yet unknown p ( a ): • There exist a critical value of coupling (or of other parameters) defined by the conditionΛ = 1 − J E [ a − a m ) ] = 0, which separates a DIAG-like spectrum from a GOE-like spectrum. • When Λ < g m of F (cid:48) ( g ) = 0 and λ m = F ( g m ).The spectrum is ρ ( λ ) ∼ √ λ − λ m close to the edge. • When Λ > p ( a ).Assuming p ( a ) ∼ A d ( a − a m ) / , we find that the lower edge is λ m = a m − J E (cid:104) a − a m (cid:105) and ρ ( λ ) ∼ A g ( λ − λ m ) / with A g = A d / Λ / .We now need to obtain information on p ( a ), i.e. on the statistics of x ∗ i in the minima of the Hamiltonian. We do soby solving the thermodynamics of the model in the T → Appendix B: Replica-symmetric solution of the model1. The partition function and the free energy The replicated partition function at finite temperature T = 1 /β (the Boltzmann constant is set to k B = 1), afterhaving averaged over the disorder in the couplings J ij and stiffnesses κ i (whose distribution p ( κ ) we leave unspecifiedfor now), and introduced the overlap matrix Q ab , is [22] Z n = (cid:90) d Q ab e − ( βJ )24 N (cid:80) nab Q ab (cid:34)(cid:90) d p ( κ ) (cid:90) d n x exp (cid:34) − βκ n (cid:88) a =1 x a + βh n (cid:88) a =1 x a − β n (cid:88) a =1 x a + ( βJ ) n (cid:88) ab x a Q ab x b (cid:35)(cid:35) N . (B1)We now assume a RS form for the Q ab matrix, Q RS ab ≡ (˜ q − q ) δ ab + q , which gives for Z n Z n = (cid:90) d Q RS ab e − ( βJ )24 Nn [˜ q +( n − q ] × (cid:90) d p ( κ ) (cid:90) d n x exp − βκ n (cid:88) a =1 x a + βh n (cid:88) a =1 x a − β n (cid:88) a =1 x a + ( βJ ) q − q ) n (cid:88) a =1 x a + ( βJ ) q (cid:32) n (cid:88) a =1 x a (cid:33) N . (B2)We rewrite the last term using an Hubbard-Stratonovich transformationexp ( βJ ) q (cid:32) n (cid:88) a =1 x a (cid:33) = (cid:90) d z √ π exp (cid:34) − z z ( βJ ) √ q n (cid:88) a =1 x a (cid:35) ≡ (cid:90) D z exp (cid:34) z ( βJ ) √ q n (cid:88) a =1 x a (cid:35) , (B3)where z ∼ N (0 , 1) is a random variable distributed according to a standard normal distribution. This relation allowsus to write Z n = (cid:90) d q d˜ qe − ( βJ )24 Nn [˜ q +( n − q ] (cid:20)(cid:90) d p ( κ ) (cid:90) D z (cid:18)(cid:90) d x e − βv eff ( x ) (cid:19) n (cid:21) N , (B4)with the definition of the effective potential : v eff ( x ) ≡ κ x + 14! x − βJ q − q ) x − ( J √ qz + h ) x . (B5)Note that the random force f = J √ qz is Gaussian distributed with zero mean and variance J q . The replicatedpartition function can finally be written as Z n = (cid:90) d q d˜ qe − NG (˜ q,q ) , (B6)with the replica-symmetric action G defined as G (˜ q, q ) = ( βJ ) n [˜ q + ( n − q ] − log (cid:20)(cid:90) d p ( κ ) (cid:90) D z (cid:18)(cid:90) d x e − βv eff ( x ) (cid:19) n (cid:21) . (B7)Assuming that q and ˜ q have been already selected using the saddle point method, we can then write the replica-symmetric free energy using the replica trick [22] log Z = lim n → n Z n , (B8)which, once the n → f RS (˜ q, q ) = ( βJ ) q − q ) − (cid:90) d p ( κ ) (cid:90) D z log (cid:18)(cid:90) d x e − βv eff ( x ) (cid:19) . (B9) 2. Saddle-point equations The saddle point equations for ˜ q and q can be found by differentiating f RS , Eq. (B9), with respect to ˜ q and q .The term to the left is trivial, whilst the second requires one to keep in mind the definition Eq. (B5) of the effectivepotential v eff ( x ) and its dependence on q and ˜ q . One gets ∂f∂ ˜ q = 0 = ⇒ ˜ q = (cid:90) d p ( κ ) (cid:90) D z (cid:10) x (cid:11) , (B10) ∂f∂q = 0 = ⇒ q = (cid:90) d p ( κ ) (cid:90) D z (cid:28) x − zx √ qβJ (cid:29) , (B11)where the bracket (cid:104)•(cid:105) denote a Gibbs average over the effective potential v eff ( x ), (cid:104)O ( x ) (cid:105) ≡ (cid:82) d x O ( x ) e − βv eff ( x ) (cid:82) d x e − βv eff ( x ) . (B12) 3. The replicon We also need to determine the transition line to the RSB phase. This is done by calculating the replicon eigenvalueof the matrix of second derivatives of the replica action [22]. The replica action is S ( Q ab ) = ( βJ ) n (cid:88) ab Q ab − log (cid:90) d n x exp (cid:32) − βκ n (cid:88) a =1 x a − β n (cid:88) a =1 x a + βh n (cid:88) a =1 x a + ( βJ ) n (cid:88) ab x a Q ab x b (cid:33) , (B13)where the overline denotes an average over p ( κ ). We wish to calculate the tensor of second derivatives of this actionwith respect to Q ab , M ab ; cd ≡ ∂S∂Q ab ∂Q cd = M (cid:18) δ ac δ bd + δ ad δ bc (cid:19) + M (cid:18) δ ac + δ bd + δ ad + δ bc (cid:19) + M , (B14)where δ ij is simply a Kronecker delta, and the last expression is the most general form that can be taken by a replica-symmetric tensor with four indices (here grouped as ab ; cd to emphasize that the first two indices are related to thefirst derivative with respect to Q ab , and the other two to the second derivative) [23]. We recall that the replicon modeis simply given by [23] λ R = M . (B15)The derivatives of the first (kinetic) term are easy to take, and one easily gets ∂S kin ∂Q ab ∂Q cd = ( βJ ) (cid:18) δ ac δ bd + δ ad δ bc (cid:19) . (B16)The derivatives of the second (interaction) term are more cumbersome. But using a compact notation, one can writethem down as ∂S pot ∂Q ab ∂Q cd = ( βJ ) (cid:104) x a x b x c x d (cid:105) n − (cid:104) x a x b (cid:105) n (cid:104) x c x d (cid:105) n ] , (B17)0with the (cid:104)•(cid:105) n averages defined, for a replica-symmetric Q ab , as (cid:104)O ( x , . . . , x n ) (cid:105) n ≡ Z eff ) n (cid:90) d n x O ( x , . . . , x n ) n (cid:89) i =1 e − βv eff ( x a ) , (B18)and v eff ( x ) has been defined in Eq. (B5). Now the overline also indicates an average over the Gaussian measure D z .In order to calculate M (and therefore the replicon), we can just observe that M = 2 M − M + 2 M , (B19)which comes from direct inspection of Eq. (B14). The kinetic part is trivially obtained from Eq. (B16) M kin1 = ( βJ ) , (B20)while, for the interaction part, we need to compute the averages M int1 = ( βJ ) (cid:104) x x x x (cid:105) n − (cid:104) x x (cid:105) n (cid:104) x x (cid:105) n ) − (cid:104) x x x x (cid:105) n − (cid:104) x x (cid:105) n (cid:104) x x (cid:105) n )+2( (cid:104) x x x x (cid:105) n − (cid:104) x x (cid:105) n (cid:104) x x (cid:105) n )] , (B21)at RS level, and for n → 0. Thanks to replica symmetry, one has (cid:104) x a x b (cid:105) n = (cid:104) x c x d (cid:105) n , ∀ a, b, c, d : a (cid:54) = b, c (cid:54) = d, (B22)so the expression for the replicon reduces to M int1 = ( βJ ) (cid:104) x x x x (cid:105) n − (cid:104) x x x x (cid:105) n + 2 (cid:104) x x x x (cid:105) n ] . (B23)We now need to calculate these averages for n → 0. The first one islim n → (cid:104) x x x x (cid:105) n = lim n → Z eff ) n (cid:18)(cid:90) d x x e − βv eff ( x ) (cid:19) (cid:18)(cid:90) d x x e − βv eff ( x ) (cid:19) (cid:18)(cid:90) d x e − βv eff ( x ) (cid:19) n − = (cid:104) x (cid:105) , (B24)where the (cid:104)•(cid:105) average is defined as in Eq. (B12). For the second term, we havelim n → (cid:104) x x x x (cid:105) n = lim n → Z eff ) n (cid:18)(cid:90) d x x e − βv eff ( x ) (cid:19) (cid:18)(cid:90) d x x e − βv eff ( x ) (cid:19) × (cid:18)(cid:90) d x x e − βv eff ( x ) (cid:19) (cid:18)(cid:90) d x e − βv eff ( x ) (cid:19) n − = (cid:104) x (cid:105) (cid:104) x (cid:105) , (B25)and for the third, one can easily get by the same logiclim n → (cid:104) x x x x (cid:105) n = (cid:104) x (cid:105) . (B26)In summary, one has for M int1 M int1 = ( βJ ) (cid:104) x (cid:105) − (cid:104) x (cid:105) (cid:104) x (cid:105) + 2 (cid:104) x (cid:105) ] = ( βJ ) ( (cid:104) x (cid:105) − (cid:104) x (cid:105) ) , (B27)and the final expression of the replicon eigenvalue, factoring out a positive constant, reads: λ R ∝ − ( βJ ) ( (cid:104) x (cid:105) − (cid:104) x (cid:105) ) . (B28)1 4. The T → limit In order to compute the spectrum, we need to find the ground state of the system in the athermal limit. In thatlimit, one can easily see that q → ˜ q linearly in T , so we define the following “athermal” overlaps and their associatedsaddle point equations, χ = β (˜ q − q ) = (cid:90) d p ( κ ) (cid:90) D z z (cid:104) x (cid:105)√ ˜ qJ , ˜ q = (cid:90) d p ( κ ) (cid:90) D z (cid:10) x (cid:11) . (B29)In the zero temperature limit, the equilibrium averages (cid:104) x (cid:105) on the effective potential v eff ( x ) become dominated by itsground state. One has then lim T → (cid:104) x (cid:105) = x ∗ ( z, κ ) , (B30)where x ∗ is the absolute minimum of the potential in Eq. (B5), which in this limit reads v eff ( x ) ≡ m x + 14! x − ( f + h ) x , (B31)with the definitions m ≡ κ − J χ and f ≡ Jz √ ˜ q , as given in the main text. The Eqs. (B29) can then be written asfollows χ = (cid:90) d p ( m ) (cid:90) d p ( f ) f x ∗ ( f, m ) J ˜ q , ˜ q = (cid:90) d p ( m ) (cid:90) d p ( f ) ( x ∗ ( f, m )) , (B32)where p ( f ) ≡ N (0 , J ˜ q ) , p ( m ) = U ( κ m − J χ, κ M − J χ ) . (B33)To solve them, one can proceed as follows. Starting from a guess for χ and ˜ q , one first generates the two randomparameters ( z, κ ), and for each realization, one finds the minimum of the effective potential, by solving the cubicequation v (cid:48) eff ( x ) = 0 . (B34)Because the potential is quartic, an analytical solution of the cubic equation can be obtained and is given explicitly inAppendix D. One then averages over the random variables ( z, κ ) to compute the r.h.s. of Eqs. (B32) and obtain newestimates of χ and ˜ q . The procedure is iterated until convergence. In Appendix D we provide the detailed algorithmswe used to obtain the phase diagram reported in the main text.We note that an alternative equation for χ , which we report in the main text and is more useful when it comes tounderstating the location λ m of the spectrum’s lower edge, can be obtained. We start from the first of Eqs. (B29), atfinite temperature, and we rewrite it as χ = (cid:90) d p ( κ ) (cid:90) D z zβJ ˜ q dd z log (cid:90) d xe − βv eff ( x ) = (cid:90) d p ( κ ) (cid:90) D z βJ ˜ q d d z log (cid:90) d xe − βv eff ( x ) = (cid:90) d p ( m ) (cid:90) d p ( f ) 1 β d d f log (cid:90) d xe − βv eff ( x ) = (cid:90) d p ( m ) (cid:90) d p ( f ) d (cid:104) x (cid:105) d f , (B35)where we used the following relation, easily obtained by integration by parts and valid for any function g ( z ): (cid:90) D zzg (cid:48) ( z ) = (cid:90) D zg (cid:48)(cid:48) ( z ) . (B36)The T → (cid:104) x (cid:105) → x ∗ (the absolute minimum of the effectivepotential) in that limit, and x ∗ is not guaranteed to be a smooth function of f . In fact, if the effective potential v eff ( x ; f, m ) has multiple minima (i.e. it is a double well), then x ∗ will jump discontinuously when the sign of the2linear term f + h changes, as the absolute minimum switches from one side of the origin to the other. This willhappen as soon as m m < x ∗ ( f ) is the solution of v (cid:48) eff ( x ∗ , f ) = 0, one has0 = dd f v (cid:48) eff [ x ∗ ( f ) , f )] = v (cid:48)(cid:48) eff [ x ∗ ( f ) , f )] d x ∗ d f − ⇒ d x ∗ d f = 1 v (cid:48)(cid:48) eff ( x ∗ ) . (B37)Adding the singular term, the proper limit of d (cid:104) x (cid:105) d f therefore isd x ∗ d f = 1 v (cid:48)(cid:48) eff ( x ∗ ) + [ x ∗ ( − h + ) − x ∗ ( − h − )] δ ( f + h ) . (B38)As long as m m ≥ 0, no DW are present, the second term vanishes and one has the equation for χχ = (cid:90) d p ( m ) (cid:90) d p ( f ) 1 v (cid:48)(cid:48) eff ( x ∗ ( f, m )) = (cid:20) v (cid:48)(cid:48) eff ( x ∗ ( f, m )) (cid:21) , (B39)i.e., in absence of DWs, χ is the average of the inverse of the curvature of the effective potential in its minimum. Thisis the equation that we report in the main text and we use to prove that λ m = 0 on the RSB transition line. Notethat the overline denotes the average over f, m which is indicated as (cid:104)•(cid:105) f,m in the main text.The last ingredient we miss is the T → βJ ) ( (cid:104) x (cid:105) − (cid:104) x (cid:105) ) = J (cid:18) d (cid:104) x (cid:105) d f (cid:19) . (B40)Therefore one has, for T → λ R = 1 − J (cid:18) d x ∗ d f (cid:19) . (B41)Notice that then, when m m < −∞ , hence the replica symmetry is automatically broken. As westate in the main text, the presence of DWs in the ensemble of effective potentials is a sufficient condition for a RSBglass transition to take place in our model, with the replicon jumping to minus infinity rather than vanishing. On thecontrary, when m m > 0, the second term in Eq. (B38) vanishes and one can simply write for the replicon λ R = 1 − J (cid:20) v (cid:48)(cid:48) eff ( x ∗ ( f, m )) (cid:21) . (B42)which is the expression given and used in the main text, valid in the RS phase and in absence of double wells. Appendix C: Analysis of the effective potential We now focus on the effective potential. In particular, we want to compute the statistics of the diagonal elements a ≡ m + J χ + x / 1. Effective potential and ground state The effective potential has the form v eff ( x ) = 14! x + m x − Hx , v (cid:48) eff ( x ) = 13! x + mx − H , v (cid:48)(cid:48) eff ( x ) = 12 x + m , (C1)with the two effective parameters m = κ − J χ , p ( m ) = U ( κ m − J χ, κ M − J χ ) ,H = h + f , p ( H ) = N ( h, J ˜ q ) . (C2)3The equation for the stationary points of the effective potential is a depressed cubic [49] of the form v (cid:48) eff ( x ) = 16 (cid:0) x + P x + Q (cid:1) = 0 , P = 6 m , Q = − H , (C3)whose discriminant is ∆ = 4 P + 27 Q ∝ m + 98 H . (C4)Hence the solutions are organised as follows: • For ∆ < x k = 2 (cid:114) − P (cid:34) 13 arccos (cid:32) Q P (cid:114) − P (cid:33) − πk (cid:35) , k = 0 , , , (C5)Note that ∆ < P < x ∗ = − Q ) (cid:114) | P | (cid:34) 13 arccos (cid:32) | Q | | P | (cid:115) | P | (cid:33)(cid:35) . (C6) • For ∆ > x ∗ = − Q ) (cid:113) − P cosh (cid:104) arccosh (cid:16) − | Q | P (cid:113) − P (cid:17)(cid:105) for P < , − (cid:113) P sinh (cid:104) arcsinh (cid:16) Q P (cid:113) P (cid:17)(cid:105) for P > . (C7) • For ∆ = 0 there are two possibilities: – P = Q = 0 and x = 0 is a triple root, i.e. v eff ( x ) = x / – P (cid:54) = 0 and then x = 3 Q/P is a single root and x = − Q/ (2 P ) is a double root.To summarize, the ground state can be written as follows: x ∗ ( m, H ) = 2sgn( H ) (cid:112) | m |F sgn( m ) (cid:18) | H | (2 | m | ) / (cid:19) , H = h + f , m = κ − J χ , (C8)with F + ( ξ ) = sinh (cid:20) 13 arcsinh ξ (cid:21) , F − ( ξ ) = (cid:40) cos (cid:2) arccos ξ (cid:3) ξ < , cosh (cid:2) arccosh ξ (cid:3) ξ > . (C9) 2. Double wells and distribution of curvatures Using the auxiliary formula (cid:90) (cid:15) − (cid:15) d p ( H ) = 12 (cid:20) erf (cid:18) h + (cid:15)J √ q (cid:19) − erf (cid:18) h − (cid:15)J √ q (cid:19)(cid:21) ≈ (cid:15) e − h J q (cid:112) πJ ˜ q , (cid:15) (cid:28) J (cid:112) ˜ q , (C10)we get as a first result the fraction of double well potentials in the ensemble, which is given by: p dw = p (∆ < 0) = (cid:90) min { ,κ M − J χ } κ m − J χ d mκ M − κ m (cid:90) √ − m / − √ − m / d p ( H ) . (C11)In particular, when κ M > h > 0, and k m − J χ → − , we have p dw ≈ √ 215 ( J χ − κ m ) / κ M − κ m p (0) , p (0) = e − h J q (cid:112) πJ ˜ q . (C12)4 FIG. 3. Sketch of the integration region for p (˜ a ). The second result is the distribution of the curvatures ˜ a in the ground state, related by a simple shift to thedistribution of diagonal elements a , ˜ a = v (cid:48)(cid:48) eff ( x ∗ ) = m + ( x ∗ ) a − J χ . (C13)We are interested in the small ˜ a behavior, which is obtained following similar steps as in the soft potential modelanalysis [3, 43, 44]. First of all, we note that the distribution of curvatures can be either gapped or gapless (thecurvature cannot be negative). The only possibility to have ˜ a = 0 (i.e., a quartic potential) is to have P = Q = 0, orequivalently m = H = 0. Hence, if m m > m M < 0, the distribution of ˜ a is gapped.We then assume that m m = κ m − J χ ≤ m M = κ M − J χ ≥ 0; in this case, the distribution of ˜ a is gapless,which implies a m = J χ and λ m = 0, as we state in the main text. We shall focus on this particular case, which isthe one relevant for the ω -transition line. In this case, the integration domain over the random variables m and H can be decomposed as sketched in Fig. 3. For m M > 0, the contribution of positive m to the cumulative distributioncan be written, when ˜ a → 0, as G + (˜ a ) = (cid:90) ˜ a d m (cid:90) d p ( H ) θ [ m + 12 x < ˜ a ] ∼ p (0)˜ a / (cid:90) d η (cid:90) d ξξ / θ [ ξ (1 + 12 y ) < p (0) 25 ˜ a / (cid:90) d η (cid:18) y (cid:19) − / = 1 . . . . × p (0)˜ a / , (C14)where we introduced ξ = m/ ˜ a , η = H/ ( ξ ˜ a ) / and y = x √ ξ ˜ a = − √ (cid:20) 13 arcsinh (cid:18) − η √ (cid:19)(cid:21) . (C15)5The contribution of negative m with ∆ > G − (˜ a ) = (cid:90) m< d m (cid:90) d p ( H ) θ [ m + 12 x < ˜ a ] θ (∆ > ∼ p (0)˜ a / (cid:90) η > / d η (cid:90) ξ> d ξξ / θ [ ξ ( − y ) < p (0) 25 ˜ a / (cid:90) η > / d η (cid:18) − y (cid:19) − / = 0 . . . . × p (0)˜ a / ,y = 2 √ η ) cosh (cid:20) 13 arccosh (cid:18) | η | √ (cid:19)(cid:21) . (C16)Finally, the contribution of negative m with ∆ < G − (˜ a ) = (cid:90) m< d m (cid:90) d p ( H ) θ [ m + 12 x < ˜ a ] θ (∆ < ∼ p (0)˜ a / (cid:90) η < / d η (cid:90) ξ> d ξξ / θ [ ξ ( − y ) < p (0) 25 ˜ a / (cid:90) η < / d η (cid:18) − y (cid:19) − / = 0 . . . . × p (0)˜ a / ,y = 2 √ η ) cos (cid:20) 13 arccos (cid:18) | η | √ (cid:19)(cid:21) . (C17)Collecting these results, we obtain G (˜ a ) = 1 . . . . × e − h J q (cid:112) πJ ˜ q ˜ a / , for κ m < J χ < κ M ,p (˜ a ) = A d ˜ a / , A d ( J ) = 3 . . . . × e − h J q (cid:112) πJ ˜ q , for κ m < J χ < κ M . (C18)This applies whenever m m < m M > 0, and implies p ( a ) ∼ A d ( a − J χ ) / , which leads to the results of themain text in terms of location and shape of the spectrum edge. Note that if ( m m = 0 , m M > 0) or ( m m < , m M = 0),one also obtains the ˜ a / law, but with a different prefactor because the contribution of positive (or negative) m isabsent. Appendix D: Drawing the phase diagram In this section we report the algorithm used to determine the ω - and ω -transition lines in the phase diagram,building up from the equations derived in the previous section. We place ourselves in the RS region of the phasediagram, with the aim of determining its boundaries. We start from the form of the RS Eqs. (B32) for ˜ q and χ ,with x ∗ ( z, m ) being the unique ground state of the effective potential (having double wells would automatically implyRSB, as detailed in Appendix B 4 and the main text), given by Eq. (C8). Because we also want to explore the limits J → h → q → J ˜ q → ˜ q, (D1)so that the equations take the form˜ q = J (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / ( x ∗ ( z, m )) ,χ = (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / zx ∗ ( z, m ) √ ˜ q , (D2)and x ∗ is now the unique solution of the equation x mx − h + (cid:112) ˜ qz = 0 , (D3)given by an expression similar to Eq. (C8). Furthermore, the equation for χ can be rewritten in the following way χ = (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / m + ( x ∗ ( z, m )) , (D4)6which is the rescaled form of Eq. (B39) and completely equivalent to its form in Eqs. (D2) everywhere in the RSphase. However, we found that this form is better behaved under numerical resolution.There are two ways to break the replica symmetry at T = 0, as discussed in the main text:(i) The replicon vanishes continuously, λ R = 0, but the minimal effective stiffness stays positive, m m > 0. Thiscase corresponds to having a GOE-like spectrum, with an ω low-frequency tail populated by delocalized modes.This is a standard RSB transition.(ii) The minimal effective stiffness vanishes, m m = κ m − J χ = 0 with λ R > 0. At this point, DW effective potentialsappear. The replica symmetry is then broken via a discontinuity in replicon eigenvalue, which jumps to −∞ beyond the transition. This case corresponds to having a DIAG-like spectrum, with a ω low-frequency tail andpartially localized modes near the edge. 1. Finding the transition lines a. ω -transition On the ω -transition line one has m m = 0. Therefore one can take Eqs. (D2), set m m = 0,˜ q = J (cid:90) m M d p ( m ) (cid:90) ∞−∞ d z √ π e − z / ( x ∗ ( z, m )) ,κ m = J (cid:90) m M d p ( m ) (cid:90) ∞−∞ d z √ π e − z / m + ( x ∗ ( z, m )) , (D5)and solve them to find h c ( J ) and ˜ q at fixed J . h c ( J ) is the critical line in this case, and we re-mind that κ m is a fixed model parameter. This can be achieved via the following numerical scheme: Algorithm 1: ω -transition linefix J ;initialize ˜ q and h ; while ˜ q and h not converging do ˜ q ← (damped) J (cid:82) m M d p ( m ) (cid:82) ∞−∞ d z √ π e − z / ( x ∗ ( z, m )) ; h ← (damped)( h − κ m + J (cid:82) m M d p ( m ) (cid:82) ∞−∞ d z √ π e − z / m + ( x ∗ ( z,m )) ) ; endResult: (˜ q, h = h c ( J )) are the values of the corresponding parameters at the transition point, for each J .In order to ensure that the transition is ω -like, one needs to prove that λ R > λ R however contains integrable singularities that could make its numerical computation unstable. We derive belowan expression for λ R that does not suffer from these problems, and furthermore proves that λ R is indeed positive atthe transition. b. ω -transition In this case, one has λ R = 0 at the transition, but differently from the previous case one still needsto determine both ˜ q and χ trough Eqs. (D2), and only then get the transition point from the λ R = 0condition. Therefore we need to find also J c ( h ) at fixed h . A slightly more complicated numeri-cal scheme, which we report below, is needed (note the update for J c which avoids bisection methods):7 Algorithm 2: ω -transition linefix h ;initialize ˜ q , χ and J c ; while ˜ q , χ and J c not converging do ˜ q ← (damped) J c (cid:82) m M m m d p ( m ) (cid:82) ∞−∞ d z √ π e − z / ( x ∗ ( z, m )) ; χ ← (damped) (cid:16) J c (cid:82) m M m m d p ( m ) (cid:82) ∞−∞ d z √ π e − z / m + ( x ∗ ( z,m )) (cid:17) ; J c ← (damped) (cid:18) J c + (cid:20) − J c (cid:82) m M m m d p ( m ) (cid:82) ∞−∞ d z √ π e − z / [ m + ( x ∗ ( z,m )) ] (cid:21)(cid:19) endResult: (˜ q, χ, J c ( h )) are the values of the corresponding parameters at the transition point. 2. Numerically stable expression for the replicon We recall the expression in Eq. (B42) of the replicon eigenvalue in the RS phase, in explicit form: λ R = 1 − J (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / (cid:104) m + ( x ∗ ( z, m )) (cid:105) . (D6)We remark that the denominator appearing in the integral is essentially ˜ a . The expression could then be equivalentlyrewritten as λ R = 1 − J (cid:90) d˜ a (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / (cid:104) m + ( x ∗ ( z, m )) (cid:105) δ (cid:18) ˜ a − m − 12 ( x ∗ ( z, m )) (cid:19) = 1 − J (cid:90) d˜ a p (˜ a )˜ a . (D7)At the transition line, the distribution of ˜ a is gapless and follows p (˜ a ) ∼ ˜ a / near its edge, as discussed in section C 2.The above expression above highlights the singularity of the integrand at ˜ a = 0; this singularity is integrable, whichproves that λ R is finite at the transition, and that the jump to −∞ is due to the singular δ ( x + h ) term in Eq. (B38),while the v (cid:48)(cid:48) eff ( x ∗ ) term always stays finite and positive. Still, the singularity could cause problems when evaluatingthe replicon numerically. Is is possible to manipulate this expression to obtain an alternative one that, while beingmore cumbersome, contains no singularities. Using the fact that in the RS phase m m ≥ x ∗ ( z, m ) is the unique solution of Eq. (D3), we can rewrite I R = (cid:82) d˜ ap (˜ a ) / ˜ a as I R = (cid:90) ∞−∞ d x (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / m + x δ (cid:18) x mx − h + (cid:112) ˜ qz (cid:19) = (cid:90) ∞−∞ d x (cid:90) ∞−∞ dˆ x π (cid:90) m M m m d p ( m ) (cid:90) ∞−∞ d z √ π e − z / m + x exp (cid:20) i ˆ x (cid:18) x mx − h + (cid:112) ˜ qz (cid:19)(cid:21) = (cid:90) ∞−∞ d x √ π ˜ q (cid:90) m M + x / m m + x / d p ( m ) 1 m exp (cid:34) − q (cid:18) x (cid:18) m − x (cid:19) x − h (cid:19) (cid:35) . (D8)If m m = 0, this expression has an integrable singularity at m = 0, which we can eliminate with an integration byparts. In doing so we obtain I R = I + I + I , I = 1∆ κ (cid:90) ∞−∞ d x √ π ˜ q ln (cid:18) m M + x (cid:19) exp (cid:34) − q (cid:18) x m M x − h (cid:19) (cid:35) , I = − κ (cid:90) ∞−∞ d x √ π ˜ q ln (cid:18) m m + x (cid:19) exp (cid:34) − q (cid:18) x m m x − h (cid:19) (cid:35) , I = − (cid:90) ∞−∞ d x √ π ˜ q (cid:90) m M + x / m m + x / d m ∆ κ (ln m ) dd m exp (cid:34) − q (cid:18) x (cid:18) m − x (cid:19) x − h (cid:19) (cid:35) , (D9)8where ∆ κ ≡ κ M − κ m . The first integral I is perfectly convergent assuming m M > 0. The second integral I , however,has some integrable singularity if m m = 0. By assuming m m = 0, we can rewrite the integral as I = 1∆ κ (cid:90) ∞−∞ d x √ π ˜ q (cid:20) ln 2 + 2 ( x ln | x | − x ) dd x (cid:21) exp (cid:34) − q (cid:18) x − h (cid:19) (cid:35) , (D10)which now contains no singularities. Finally we need to consider I . Integrating again the logarithmic singularity byparts, one obtains I = − κ (cid:90) ∞−∞ d x √ π ˜ q (cid:20)(cid:18) m M + x (cid:19) ln (cid:18) m M + x (cid:19) − (cid:18) m M + x (cid:19)(cid:21) dd m exp (cid:34) − q (cid:18) x mx − h (cid:19) (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m = m M + 1∆ κ (cid:90) ∞−∞ d x √ π ˜ q (cid:20)(cid:18) m m + x (cid:19) ln (cid:18) m m + x (cid:19) − (cid:18) m m + x (cid:19)(cid:21) (cid:34) dd m exp (cid:34) − q (cid:18) x mx − h (cid:19) (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m = m m (cid:35) + (cid:90) ∞−∞ d x √ π ˜ q (cid:90) m M m m d m ∆ κ (cid:18)(cid:18) m + x (cid:19) ln (cid:18) m + x (cid:19) − (cid:18) m + x (cid:19)(cid:19) d d m exp (cid:34) − q (cid:18) x mx − h (cid:19) (cid:35) . (D11)In summary, we have reduced the computation of the replicon integral I R to the sum of perfectly convergent,singularity-free integrals that can be easily evaluated numerically. For the purpose of numerical integration (suchas in the case of the algorithms reported above), it is convenient to rescale the integration variable x by √ ˜ q and alsodo the same on h : x → x (cid:112) ˜ q h = Γ (cid:112) ˜ qq