Path Integral Approach Unveils the Role of Complex Energy Landscape for Activated Dynamics of Glassy Systems
PPath Integrals for Activated Dynamics in Glassy Systems
Tommaso Rizzo
ISC-CNR, UOS Rome, Universit`a “Sapienza”, Piazzale A. Moro 2, I-00185, Rome, Italy andDip. Fisica, Universit`a “Sapienza”, Piazzale A. Moro 2, I–00185, Rome, Italy
The Random-First-Order-Transition theory of the glass transition stems from the fact that mean-field models of spin-glasses and supercooled liquids display an exponential number of metastablestates that trap the dynamics. In order to obtain quantitative dynamical predictions to assesthe validity of the theory I discuss how to compute the exponentially small probability that thesystem jumps from one metastable state to another in a finite time. This is expressed as a pathintegral that can be evaluated by saddle-point methods in mean-field models, leading to a boundaryvalue problem. The resulting dynamical equations are solved numerically by means of a Newton-Krylov algorithm in the paradigmatic spherical p -spin glass model. I discuss the solutions in theasymptotic regime of large times and the physical implications on the nature of the ergodicity-restoring processes. More than thirty years after its formulation theRandom-First-Order-Transition (RFOT) theory [1, 2] isstill one of the major competing theories in the ongoingdebate on the nature of the Glass transition. In a nutshellthe theory posits that the physics of supercooled liquids isthe same of Spin-Glass (SG) models displaying one-stepof Parisi’s Replica-Symmetry-Breaking (1RSB) [3]. Themean-field versions of these models display an ergodicity-breaking transition at a dynamical temperature T d wherethe phase space splits into many metastable states thattrap the dynamics; at lower temperatures, the configura-tional entropy, i.e. the log of the number of metastablestates, decreases eventually vanishing at a static tem-perature T s . Ergodicity breaking between T d and T s isa mean-field artifact and one expects that in real sys-tems ergodicity is restored through droplet-like excita-tions. Furthermore the size of these excitations mustdiverge as the configurational entropy vanishes leadingeventually to a genuine ergodicity-breaking transition at T s . RFOT originated from the realization [4] that simi-lar features had been discussed in various unrelated (andthemselves controversial) theories of supercooled liquids,most notably: 1) dynamics at the ergodicity-breakingtransition is the same of the (avoided) Mode-Coupling-Theory (MCT) of supercooled liquids [5], 2) within theAdams-Gibbs-Di Marzio theory [6, 7], the glass transi-tion is driven by a correlation length that diverges atthe Kauzmann temperature where the configurational en-tropy vanishes.Efforts to validate the theory have been driving the-oretical, numerical and experimental research for years.At the theoretical level RFOT has been substantiatedby a number of results arguing that mean-field mod-els of supercooled liquids exhibit 1RSB [8, 9] includingthe solution of the limit of infinite physical dimensions[10, 11]. Numerically, the MCT phenomenology is welldocumented [12, 13] as well as the increase of dynamic[14, 15] and static correlation lengths [16, 17]. Experi-mentally, the observation of a decreasing configurationalentropy dates back to the 40’s while more recently RFOThas inspired measurements of non-linear susceptibilities[18]. Furthermore the analogy with supercooled liquids FIG. 1: Supercooled liquids display caging and hopping inreal space (left) that correspond to evolution in the ruggedphase space of mean-field models (right). has also led to the discovery that off-equilibrium relax-ational dynamics of mean-field Spin-glass models displayaging [19] and it is an active line of research [20, 21].In spite of this huge body of work, consensus on thevalidity of the theory is still lacking. One of the problemis that RFOT-inspired literature often focuses on quanti-ties, namely point-to-set correlation lengths and configu-rational entropy, whose actual relevance for the problemof the glass transition, that instead is dynamical in itsessence, can also be questioned. Besides, predictions areoften merely qualitatively, which is a problem given that e.g. the observed static length increases are too mod-est to convince the community that they actually drivethe slowing down of the dynamics. To make progressit would be important to obtain RFOT predictions thatare both quantitative and dynamical . The main challengeis that current theoretical knowledge is mostly limitedto mean-field models that display ergodicity-breaking at T d : in order to recover ergodicity between T d and T s and obtain realistic predictions we have to go beyondmean-field and this line of research is currently very ac-tive [22–28]. In recent years progress has been made forthe region close to the dynamical temperature T d [29–31].It is now possible to describe qualitatively and quanti-tatively how the ergodicity-breaking MCT transition isturned into a dynamical crossover in mean-field SG [32](due to finite-size effects) and most importantly in some a r X i v : . [ c ond - m a t . d i s - nn ] J a n FIG. 2: left to right, top to bottom: density plot of the sixfunctions C ( t, t ), C dt ( t, t ), T ˆ R ( t, t ), T ˆ R dt ( t, t ), T ˆ X ( t, t ), T ˆ X dt ( t, t ) for −T ≤ t, t ≤ T and T = 40 for the Spherical3-SG model at T = 1 / . < T d , (numerical solution with∆ t = T / finite-dimensional models [33]. In this letter I considerinstead the region between T d and T s where metastablestates are present and discuss how to compute the expo-nentially small probability of a jump from an equilibriumstate to another occurring in a finite time. To makecontact with the phenomenology of supercooled liquidswe have to remember that below the experimental MCTtransition temperature a particle is trapped most of thetime into a cage formed by the surrounding particles anddiffusion occurs through hopping i.e. sudden rare jumpsfrom a cage to another (see fig. 1). Mean-field mod-els capture caging through the appearance of metastablestates and the study of jumps in the free energy land-scape initiated in this paper is essential to a quantitativedescription of hopping in real space.The main results to be discussed are: i) a pathintegral method to compute the transition rate ii) aNewton-Krylov algorithm yielding the numerical solu-tion of the corresponding dynamical equations iii) anasymptotic analysis of the solutions in the regime oflarge times and iv) some non-trivial implications on theergodicity-restoring processes. At the methodological level the problem is formulated as a path integral overLangevin dynamic trajectories that can be computedthough saddle-point methods in mean-field models. Theproblem displays some important differences with re-spect to the standard relaxational dynamics [19, 34, 35],namely the use of replicas and the need to explicitly han-dle the divergent path integral. I have focused on theparadigmatic spherical p -spin SG model [35, 36] but themethod is fairly general and the equations can be de-rived with some effort for other mean-field systems e.g. supercooled liquids in large dimension [37–39]. Anotherimportant difference with ordinary relaxational dynam-ics is the fact that one ends up with a boundary valueproblem instead of an initial value problem [46]. Thusthere is no trivial iterative algorithm to solve the corre-sponding dynamical equations numerically. A successfulsolution strategy has been developed based on three ele-ments: 1) Newton’s method on discretized equations, 2)Krylov methods with physical preconditioning to invertthe Jacobian, 3) Richardson extrapolation to reach thecontinuum limit. The whole technology can be again ex-ported to other problems, with possible computationalcomplexity issues due to the specific dynamical order pa-rameter.The main object considered is the transition rate T T ( σ | τ ) defined as the probability that the system isin configuration σ at time t fin = T given that it was inconfiguration τ at time t in = −T . It is convenient to ac-tually consider the following object that, due to detailedbalance, is symmetric with respect to the exchange of σ and τ : ˆ T T ( σ, τ ) ≡ T T ( σ | τ ) e β H ( σ ) − β H ( τ ) . (1)An integral representation of Langevin dynamics is usedand σ , τ are chosen as generic equilibrium configurations.At low temperatures ˆ T ( σ, τ ) is exponentially small inthe system size N in mean-field models, therefore an an-nealed average P eq ( σ ) P eq ( τ ) ˆ T ( σ, τ ) would interfere withthe equilibrium measure of σ and τ and the correct pro-cedure is to consider the quenched average[ln ˆ T ] ≡ X σ,τ P eq ( σ ) P eq ( τ ) ln ˆ T ( σ, τ ) . (2)One can resort to the replica method to eliminate the log-arithm. If quenched disorder is present the correspondingaverages (represented by an overline in the following) re-quire the introduction of additional replicas of the initialand final configurations. One can argue that the rate isself-averaging, meaning that most couples ( σ, τ ) displaythe same rate [47], thus, given an initial configuration σ ,the total number of configurations with rate equal to theaverage [ln ˆ T T ] is equal to the total number of equilib-rium configurations τ , i.e. the exponential of the entropy e S . Therefore the total probability of jumping to one ofthe configurations with typical rate is e [ln ˆ T ]+ S and mustbe smaller than one leading to the bound:[ln ˆ T T ] + S ≤ correlated with theinitial one, therefore the probability to be in a generic equilibrium configuration (that is uncorrelated to the ini-tial one) must be smaller than one meaning that at anyfinite T the above bound should not be saturated. Onthe other hand ergodicity implies that when T goes toinfinity the probability measure will be flat over the e S equilibrium configurations and the average rate shouldbecome equal to e − S saturating the bound.I have obtained an expression for the average transitionrate of the spherical p -spin model for T > T s (see the sup-plemental material (SI)) that depends on six real func-tions C ( t, t ), ˆ R ( t, t ), ˆ X ( t, t ), C dt ( t, t ), ˆ R dt ( t, t ) andˆ X dt ( t, t ) defined for −T ≤ t, t ≤ T , see fig. 2. Two ad-ditional functions µ ( t ) and ˆ µ ( t ) enforce the spherical con-straint leading to C ( t, t ) = 1 and ˆ R ( t, t ) = 1 / t .Extremization of the expression leads to eight non-linearintegro-differential equations (albeit two are redundantdue to symmetries). The physical meaning of C ( t, t )and C dt ( t, t ) is straightforward: C ( t, t ) = [ h s i ( t ) s i ( t ) i ] (4) C dt ( t, t ) = [ h s i ( t ) ih s i ( t ) i ] (5)where the square brackets mean averages with respectto the dynamical trajectories at fixed initial and finalconfigurations ( σ, τ ) and fixed disorder. Thus C ( t, t )is the average correlation between configurations vis-ited by the same trajectory at times t and t while C dt ( t, t ) is the correlation between configurations visitedby different trajectories (hence the suffix dt ). It followsthat the equations must satisfy the boundary conditions C ( ±T , ±T ) = C dt ( ±T , ±T ) = 1 and C ( ±T , ∓T ) = C dt ( ±T , ∓T ) = 0.While ergodicity implies that the rate must got to − S in the T → ∞ limit, the use of the saddle-point methodimplies that the thermodynamic limit N → ∞ is alwaystaken first. Quite naturally the two limits cease to com-mute for T < T d : above T d one can show that in the T → ∞ limit the saddle-point equations admit a solutionin which the six functions are expressed in terms of theequilibrium correlation C eq (0 , t ) and that the rate tendsto − S . This is possible because C eq (0 , ∞ ) = 0 and theboundary condition C ( ±T , ∓T ) = 0 can be satisfied for T → ∞ by the equilibrium solution, while it is no longertrue for
T < T d because equilibrium dynamics is trappedand C eq (0 , ∞ ) = 0, as a consequence for T < T d the aver-age log-rate is smaller than − S also in T → ∞ limit. Forany finite T the solutions cannot be expressed in terms ofthe equilibrium correlation. The equations can be solvedanalytically in the free case ( T = ∞ ) in which the systemperforms a Brownian motion on the N − - -
20 0 20 400.00.20.40.60.81.0 - - FIG. 3: Spherical 3-SG at T = 1 / . < T d . Top, left: C ( −T , t ) vs. t + T for T = 8 , , , ,
40, the dashed lineis the equilibrium dynamics that is trapped inside a state andhas a plateau at 0 . t/ T . Bottom, left: C (0 , t ) vs. t for T = 8 , , , ,
40, right: same vs. t/ T .Each curve was obtained from finite N t solutions throughRichardson extrapolation (see text). the harmonic oscillator formulas to derive the expressionof the rate as a function of T . Knowledge of the rate inthe non-interacting case provides an alternative way tocompute the rate at finite temperature without having todeal with the infinite normalization factor. Differentiat-ing the saddle-point expression of the rate with respectto the inverse temperature β one gets indeed a finite ex-pression that yields the correct derivative when evaluatedon the solutions of the saddle-point equations. The rateat finite T and β can then be obtained by integration in β starting from the exact β = 0 result.To solve the equations numerically one discretizes timewith ∆ t = T /N t for some integer N t . The six functionscan then be jointly represented as two real matrices ofsize 2(2 N t + 1) so that the total number of variables is O (32 N t ) (symmetries reduces this number by four). Us-ing standard discrete formulas for integrals and second-order derivatives one obtains expressions with an O (∆ t )error. Newton’s method turned out to be able to solvethe equations. A starting function not too far away fromthe correct solution can be obtained from the analyticsolutions in the free case or the T → β = 0 at finite T and switch on the temper-ature, then at fixed the temperature, T can be changedchanging the discretization parameter ∆ t at fixed N t bysmall amounts. At fixed T and β , ∆ t can be reducedextrapolating the result of a coarser grid to a finer grid(larger N t ) and using it as a starting point for Newton’smethod at the new N t . The main technical problem isthat every iteration of Newton’s method requires to in-vert the Jacobian of the equations, a (32 N t ) × (32 N t )matrix. Even exploiting the symmetries of the problem,with current technology exact inversion of the Jacobianbecomes unfeasible for N t of the order 70 −
80 limitingthe values of T that can be studied. N t could be in-creased up to 2000 using an approximate method for thewell-studied problem [42] of solving a very large linearsystem A x = b . Specifically I used the Generalized Min-imal Residue (GMRES) algorithm [42, 43] that requiresthe computation of the Krylov subspace of order k , de-fined by the vectors { b, A b, A b, A b, . . . , A k b } . The ad-vantage of Krylov methods is that one always work withvectors v and has to perform matrix multiplications A v without the need to store the full (32 N t ) × (32 N t ) matrix A . In GMRES one searches for an approximate solutionin the Krylov sub-space of order k , an orthonormal basisof the Krylov space is obtained using a numerically stableGram-Schmidt orthogonalisation called Arnoldi iterationand the solution is then found by least squares minimiza-tion of the linear equations in this space. The advantageof the method is that the error decreases systematicallyincreasing k , the drawback is that it requires to gener-ate and store all the k vectors of the basis and k cannotbe too small to obtain accurate solutions. This effectscan be reduced using the aforementioned symmetries ofthe problem that allows to work with smaller values of k and give a four-fold decrease of the memory requiredto store the vectors. Efficient Krylov methods often re-quire preconditioning which amount to reduce that thespan of the eivengalues of A . In our case the equationsare badly ill-conditioned due to the presence of second-derivatives leading to an unbounded spectrum in the con-tinuum limit. This can be solved multiplying them bythe inverse of the free part prior to iterative solution tohave a Jacobian with a discrete and bounded spectrum.To compute the equations and the products J c involv-ing the Jacobian we have to store huge matrices andperform multiplications and element-wise operations, forwhich many highly-optimized and parallelized librariesexist. I wrote a code using Mathematica (see the ancil-lary files section) being able to reach values of N t = 2000with k = 400 using up to 90% memory on cluster with256 Giga of RAM. Note that for more complex problemone could use a Jacobian-free method (involving only theequations) in which the product J v is approximated by E ( c + (cid:15) v ) − E ( c )) /(cid:15) for small (cid:15) . Other algorithms existthat do not require to store the whole Krylov Space, e.g. Biconjugated Gradient [42]; in general they are less safeand controlled than GMRES, but an efficient algorithmrequiring less memory would be welcomed. FurthermoreGMRES allows to have an approximation for the spec-trum of the Jacobian, which is also very useful informa-tion to devise alternative algorithms. Once the the nu-merical solutions at fixed β and T is obtained for various N t , a polynomial (Richardson) extrapolation is essentialto reach the ∆ t = 0 limit and remove a few pathologiesof the finite ∆ t solutions (SI). Extrapolations have somedegree of arbitrariness but, with some care, yield reliableresults.The solutions at finite values of T allow to guess the asymptotic behavior for T → ∞ in the activated re-gion, see fig. 3. For finite time differences | t − t | = O (1) (cid:28) T the functions C ( t, t ), ˆ R ( t, t ) and ˆ X ( t, t )verify the equilibrium relationships corresponding to fluctuation-dissipation theorem and time-translationalinvariance meaning that on finite time-scales the trajec-tories are essentially equilibrium trajectories in a self-induced slowly-varying field. On the large time scales | t − t | = O ( T ) the solutions are expressed in terms of uni-versal functions as C ( t, t ) = C u ( t/ T , t / T ), C dt ( t, t ) = C dtu ( t/ T , t / T ), ˆ R ( t, t ) = T − ˆ R u ( t/ T , t / T ), ˆ R dt ( t, t ) = T − ˆ R dtu ( t/ T , t / T ), ˆ X ( t, t ) = T − ˆ X u ( t/ T , t / T ),ˆ X dt ( t, t ) = T − ˆ X dtu ( t/ T , t / T ), µ ( t ) = µ u ( t/ T ), ˆ µ ( t ) = T − ˆ µ u ( t/ T ). Note that ˆ R ( t, t ) and ˆ X ( t, t ) are smallwhile they are finite close to the diagonal, see fig. 2.Plugging the above expressions in the equations one seesthat the universal functions solve the dynamical equa-tions with the second derivates dropped, similarly towhat happens for equilibrium dynamics at T d and foroff-equilibrium dynamics [19, 35] . Closing those equa-tions would allow to work directly at T → ∞ and is anopen problem that is left for future work.One should note that the asymptotic structure of thesolutions is utterly different from that of metastabilityin ferromagnetism. The corresponding computation de-scribes the transition rate from the metastable minimumto the stable one and leads to an instantonic equation inwhich the second order derivatives is not dropped in theasymptotic limit. As a consequence even if
T → ∞ thejump effectively occurs in time window centered around t = 0 that remains finite and does not scale with T .Another more striking difference with metastability inferromagnetism occurs when we consider the ergodicity-restoring processes. In this work I focused on the ex-ponentially small probability that the system jumps toanother equilibrium state in a finite time, this is com-plementary to the problem of determining the exponen-tially large time-scale τ erg over which the system jumps toanother equilibrium state with finite probability. Giventhat the exponentially small probability to jump to anequilibrium state with typical rate is p = e [ln T ]+ S it isnatural to expect that such a probability becomes finiteon a time-scale of order 1 /p . This is indeed what hap-pens for ferromagnetism, the following discussion showsthat instead it is not true in presence of an exponentialnumber of metastable states, i.e. a finite configurationalentropy. The computations presented here can indeedbe generalized to give the total probability of jumpingto states with higher energy, taking into account theirdifferent entropy and rate. In particular the derivativeof p ( E ) with respect to E at the equilibrium value canbe easily expressed in terms of the solutions and turnsout to be positive. This means that there is a range ofvalues of the energy E > for which p ( E > ) (cid:29) p ( E ) andthe system has a finite probability to jump to one ofthese states on a time-scale 1 /p ( E > ) (cid:28) /p . On theother hand detailed balance implies that the probabilityto jump from a higher energy state back to an equilib-rium state is p ( E > ) e ∆ S − β ∆ E and thus it is exponentiallylarger than p ( E > ) given that, by definition, S − βE ismaximal on the equilibrium states. Therefore, after thesystem has jumped to a higher energy metastable state itwill jump back to a generic equilibrium state on a scale ex-ponentially smaller than 1 /p ( E > ), i.e. instantaneously.This implies that an intermediate jump to one of theexponentially many metastable states provides a moreefficient path for restoring ergodicity than a direct jumpto another equilibrium state and thus the ergodic scaleis smaller than 1 /e [ln T ]+ S , at variance with magnetism.Therefore to determine τ erg one should (at least) maxi-mize p ( E ), a detailed analysis is left for future work. Acknowledgments
I acknowledge the financial support of the SimonsFoundation (Grant No. 454949, Giorgio Parisi). I thankE. Zaccarelli for substantial help with computing re-sources. [1] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes,Physical Review A , 1045 (1989).[2] P. G. Wolynes and V. Lubchenko, Structural glasses andsupercooled liquids: Theory, experiment, and applications (John Wiley & Sons, 2012).[3] M. M´ezard, G. Parisi, and M. Virasoro,
Spin glass the-ory and beyond: An Introduction to the Replica Methodand Its Applications , vol. 9 (World Scientific PublishingCompany, 1987).[4] T. Kirkpatrick and D. Thirumalai, Phys. Rev. B , 5388(1987).[5] W. G¨otze, Complex dynamics of glass-forming liquids: Amode-coupling theory , vol. 143 (OUP Oxford, 2008).[6] J. H. Gibbs and E. A. DiMarzio, The Journal of ChemicalPhysics , 373 (1958).[7] G. Adam and J. H. Gibbs, The journal of chemicalphysics , 139 (1965).[8] M. M´ezard and G. Parisi, Structural Glasses and Super-cooled Liquids: Theory, Experiment, and Applicationspp. 151–191 (2012).[9] R. Monasson, Physical review letters , 2847 (1995).[10] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, andF. Zamponi, Nature communications , 1 (2014).[11] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, andF. Zamponi, Annu. Rev. Condens. Matter Phys. , 265(2017).[12] W. Kob and H. C. Andersen, Physical Review E , 4626(1995).[13] W. Kob, J. Condens. Matter Phys. , R85 (1999).[14] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, andS. C. Glotzer, Physical review letters , 2827 (1997).[15] E. Flenner, H. Staley, and G. Szamel, Physical reviewletters , 097801 (2014).[16] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti,D. El Masri, D. L’Hˆote, F. Ladieu, and M. Pierno, Sci-ence , 1797 (2005).[17] G. Biroli, J.-P. Bouchaud, A. Cavagna, T. S. Grigera,and P. Verrocchio, Nature Physics , 771 (2008).[18] S. Albert, T. Bauer, M. Michl, G. Biroli, J.-P. Bouchaud,A. Loidl, P. Lunkenheimer, R. Tourbot, C. Wiertel-Gasquet, and F. Ladieu, Science , 1308 (2016).[19] L. F. Cugliandolo and J. Kurchan, Physical Review Let-ters , 173 (1993).[20] G. Folena, S. Franz, and F. Ricci-Tersenghi, Physical Re-view X , 031045 (2020).[21] A. Altieri, G. Biroli, and C. Cammarota, Journal ofPhysics A: Mathematical and Theoretical , 375006(2020).[22] M. Baity-Jesi, A. Achard-de Lustrac, and G. Biroli, Physical Review E , 012133 (2018).[23] M. Baity-Jesi, G. Biroli, and C. Cammarota, Journalof Statistical Mechanics: Theory and Experiment ,013301 (2018).[24] M. R. Carbone, V. Astuti, and M. Baity-Jesi, PhysicalReview E , 052304 (2020).[25] I. Hartarsky, M. Baity-Jesi, R. Ravasio, A. Billoire, andG. Biroli, Journal of Statistical Mechanics: Theory andExperiment , 093302 (2019).[26] V. Ros, G. Biroli, and C. Cammarota, EPL (EurophysicsLetters) , 20003 (2019).[27] D. A. Stariolo and L. F. Cugliandolo, EPL (EurophysicsLetters) , 16002 (2019).[28] D. A. Stariolo and L. F. Cugliandolo, Phys. Rev. E ,022126 (2020).[29] T. Rizzo, EPL (Europhysics Letters) , 56003 (2014).[30] T. Rizzo and T. Voigtmann, EPL (Europhysics Letters) , 56008 (2015).[31] T. Rizzo, Phys. Rev. B , 014202 (2016).[32] T. Rizzo, Philosophical Magazine , 636 (2016).[33] T. Rizzo and T. Voigtmann, Physical Review Letters , 195501 (2020).[34] H. Sompolinsky and A. Zippelius, Physical Review B ,6860 (1982).[35] A. Crisanti, H. Horner, and H.-J. Sommers, Zeitschriftf¨ur Physik B Condensed Matter , 257 (1993).[36] A. Crisanti and H.-J. Sommers, Zeitschrift f¨ur Physik BCondensed Matter , 341 (1992).[37] J. Kurchan, T. Maimbourg, and F. Zamponi, Journalof Statistical Mechanics: Theory and Experiment ,033210 (2016).[38] T. Maimbourg, J. Kurchan, and F. Zamponi, Physicalreview letters , 015902 (2016).[39] A. Manacorda, G. Schehr, and F. Zamponi, The Journalof Chemical Physics , 164506 (2020).[40] J. Zinn-Justin, Quantum Field Theory and Critical Phe-nomena (Oxford Science Publications, 2002).[41] G. Parisi,
Statistical field theory (Addison-Wesley, 1988).[42] Y. Saad,
Iterative methods for sparse linear systems (SIAM, 2003).[43] Y. Saad and M. H. Schultz, SIAM Journal on scientificand statistical computing , 856 (1986).[44] A. Lopatin and L. Ioffe, Physical Review B , 6412(1999).[45] A. Lopatin and L. Ioffe, Physical review letters , 4208(2000).[46] This is also the main difference with the earlier methodof [44, 45] that allows to study some special activatedprocesses through relaxational dynamics equations but cannot be used to study transition rates between genericequilibrium states.[47] This can be shown computing O ( n ) corrections to the rate where n is the replica number. Supplemental Materials: Path Integrals for Activated Dynamics in Glassy Systems
Contents
I. Path Integral Expression of the Rate p -spin Spherical Model: The Action and Saddle-Point Equations 51. Integrals 8 II. Analysis of The Saddle-Point equations of the Spherical p -Spin-Glass Model III. Numerical Solution of the Equations References I. PATH INTEGRAL EXPRESSION OF THE RATEA. Replicas
We want to compute ˆ T ( σ, τ ) when σ and τ are equilibrium configurations. Furthermore ˆ T ( σ, τ ) is exponentiallysmall in the system size N in the activated region, therefore an annealed average P eq ( σ ) P eq ( τ ) ˆ T ( σ, τ ) would interferewith the measure of σ and τ . Thus, as explained in the main text, the correct procedure is to consider the quenchedaverage X σ,τ P eq ( σ ) P eq ( τ ) ln ˆ T ( σ, τ ) . (1)The logarithm can be eliminated with the replica method. We will consider systems with quenched disorders whoseequilibrium properties can also be studied by the replica method, thus we introduce the following object: Z ˆ T ( m, m , n ) ≡ X σ ...σ m X τ ...τ m exp − β m X i =1 H J ( σ i ) − β m X i =1 H J ( τ i ) ˆ T ( σ , τ ) n (2)and we have: X σ,τ P eq ( σ ) P eq ( τ ) ln ˆ T ( σ, τ ) = lim m → lim m → lim n → ddn ln Z ˆ T ( m, m , n ) (3)As usual in the context of mean-field model we will perform the thermodynamic limit before the above limits. Notethat, as usual, the study of lim m → ,m → Z ˆ T ( m, m , n ) at finite n allows to study the large deviations of ln ˆ T ( σ, τ ) a r X i v : . [ c ond - m a t . d i s - nn ] J a n and determine if it is self-averaging with respect to the equilibrium configurations σ and τ . As we will recall in thefollowing it is possible to obtain an integral representation of the dynamics (see [1], chapter 17):ˆ T ( σ , τ ) = Z s ( −T )= σ ,s ( T )= τ ds exp (cid:20) Z d d s (1)Γ(1 , s (2) − β Z d H J ( s (1)) (cid:21) . (4)That leads to: Z ˆ T ( m, m , n ) = X σ ...σ m X τ ...τ m n Y a =1 "Z s a ( −T )= σ ,s a ( T )= τ [ ds a ] exp " n X a =1 Z d d s a (1)Γ(1 , s a (2) − β m X i =1 H J ( σ i ) − β m X i =1 H J ( τ i ) − β n X a =1 Z d H J ( s a (1)) , (5)this in turn can be written in a compact form as: Z ˆ T ( m, m , n ) = Z [ d s ] exp (cid:20) Z d1d2s ( ) Γ ( , ) s ( ) − β Z d H J ( s ( )) (cid:21) (6)where the bold index runs over the replicas and the dynamical indexes s ( ) = { σ , . . . , σ m , s (1) , . . . , s n (1) , τ , . . . , τ m } (7)and we have: Z f ( s ( )) d = m X i =1 f ( σ i ) + m X i =1 f ( τ i ) + n X a =1 Z d f ( s a (1)) . (8)The dynamical operator Γ ( , ) is diagonal withe respect to the replica indexes and associates each replicas with theboundary conditions at σ and τ . As usual the great advantage of having the above compact representation is thatone can perform the average of the disorder and then perform standard manipulations yielding an expression formallyidentical to the one obtained in the case of a static replica computation. B. Path Integral Representation of Langevin Dynamics
We start from the Langevin equation1Γ ˙ q = − β dHdq + ξ , h ξ ( t ) ξ ( t ) i = 2Γ δ ( t − t ) . (9)and we discretize it as: 1Γ q i + i − q i ∆ t = − β (cid:18) c dHdq i + (1 − c ) dHdq i +1 (cid:19) + ξ i (10)where 0 ≤ c ≤ c = 1, in Stratonovich we have c = 1 / T ( σ | τ ) = Z N − Y i =1 dq i ! N − Y i =0 d ˆ q i π ! exp[ − N − X i =0 (cid:18) ˆ q i E i ∆ t − ln (cid:20) − (1 − c ) β d Hdq i +1 ∆ t (cid:21)(cid:19) ] (11)where the time is divided in N intervals of size ∆ t . The logarithm comes from the determinant that gets contributiononly from the diagonal since the Jacobian is a triangular matrix. Expanding the Jacobian at first order in ∆ t weobtain: T ( σ | τ ) = Z N − Y i =1 dq i ! N − Y i =0 d ˆ q i π Γ ! exp[ −L ] (12)The Lagrangian reads: L = N − X i =0 (cid:18) ˆ q i E i ∆ t − (1 − c ) β d Hdq i +1 Γ ∆ t (cid:19) (13)and in the continuum limit we have L = Z dt (cid:20) (cid:0) ˆ q ˙ q − ˆ q (cid:1) + ˆ qβ dHdq − (1 − c )Γ β d Hdq (cid:21) . (14)Note that unexpectedly the continuum limit expression depends on the microscopic parameter c of the discretizationwhile one would expect it to be irrelevant. This is a well-known ambiguity of path integral representation of stochasticequations. One can choose to use the Ito discretization corresponding to c = 1 and neglect it but it will appear laterin the computation. In the following we prefer to keep it also to remind us that the continuum limit of stochasticequations is delicate, the ordinary rules of calculus (integration by parts, differentiation, chain rules) are modified forstochastic processes. Beside we will use Hamiltonians where the interaction part is just linear (the p -spin interactions)and thus in the end we will go back to special Langevin equation for single variable. Let us consider the symmetricrate defined as ˆ T ( σ, τ ) ≡ T ( σ | τ ) e β ( H ( σ ) − H ( τ )) (15)In the continuum limit we would expect the following to be an equality β H ( σ ) − H ( τ )) = β Z T − T dt ˙ q ∂Hdq . (16)Instead in order to get the correct expression we should go back to the discretized expression. We have: H ( q i +1 ) − H ( q i ) = dHdq i ∆ q i + 12 d Hdq i ∆ q i (17)then we have to use the fact that in the Lagrangian we use the following discretized definition of dH/dqc dHdq i + (1 − c ) dHdq i +1 . (18)By rewriting the differential as dHdq i = (cid:18) c dHdq i + (1 − c ) dHdq i +1 (cid:19) − (1 − c ) d Hdq i ∆ q i (19)we obtain dH = dHdq dq + (cid:18) c − (cid:19) d Hdq dq . (20)We can see that for c = 1 we recover Ito’s lemma while for the Stratonovich prescription c = 1 / O ( dt ) contribution but we canmake the replacement dq = 2 Γ dt (21)and obtain β H ( σ ) − H ( τ )) = β Z TT dt (cid:18) dHdq ˙ q + (cid:18) c − (cid:19) d Hdq (cid:19) . (22)The same result can also be obtained reabsorbing the term ∆ q in the term ˙ q below (see e.g. [1], section 4.6). Makingthe following change of variable ˆ q = ˆ x + ˙ q L ( t ) = 1Γ (cid:18) ˙ q − ˆ x (cid:19) + ˆ xβ dHdq −
12 Γ β d Hdq (24)we can also integrate out the ˆ x and obtain:ˆ L ( t ) = 1Γ ˙ q (cid:18) β dHdq (cid:19) −
12 Γ β d Hdq (25)See also [1], pag. 70. Note that the integration over ˆ x leads to a divergent prefactor to the path integral:12 π Γ Z d ˆ xe ˆ x dt = 12 π Γ (cid:18) π Γ dt (cid:19) / . (26)The divergent prefactor is usually buried into the expression [ dq ] defined as:[ dq ] ≡ dq (cid:18) / (2Γ )2 πdt (cid:19) N/ (27)see eq. 2.20 in [1]. C. Path Integral Representation for Models with Multi-Linear Interactions
Typical mean-field SG models have multi-linear p -spin interactions, while the non-linear part of the Hamiltonian islocal and does not need to be decoupled in order to obtain a saddle point expression. This allows to use simplifiedintegral representations of the dynamics in which fewer variables are introduced with respect to the general casediscussed in [1], chapter 17. In the following we introduce a bosonic variable η that behaves as the product of twoGrassmann variables, that is we have: η = 0 , Z dη = 0 , Z η dη = 1 (28)with it we define a new coordinate 1 ≡ ( t, η ) and field s (1): s i (1) ≡ s i ( t ) + ˆ x i ( t ) η . (29)With the above definitions the interacting part of the multi-linear Hamiltonian can be written as: β Z dt X i ˆ x i ( t ) dHds i = β Z d H ( s (1)) . (30)A similar formulation is also useful in the spherical model. In this case however the non-linear part of the Hamiltonianis due to the spherical constraint which is not local and must be treated appropriately. Let us consider the problemof the integral representation of Langevin dynamics of N real spins s i constrained on a N − G ( s ) = 0. The statics of the problem can be written as Z d N s δ ( G ) |∇ G | e − βH ( s ) (31)A convenient way to define Langevin dynamics on the surface is to relax the delta function replacing it with a Gaussianof infinitesimal variance (cid:15) . Then one have to compute the standard dynamical integral in presence of a Hamiltonian H ( x ) = G ( x ) / (2 N (cid:15) ). For the spherical constraint on N continuous spins s i we have G ≡ N X i =1 s i − N (32)the term |∇ G | is exactly equal to N and can be ignored, For H ( x ) = G ( x ) / (2 N (cid:15) ) we have: X i ˆ x i β dHdq i −
12 Γ β X i d Hdq i = 2 N (cid:15) Gβ X i ˆ x i s i − Γ N (cid:15) Gβ = βN (cid:15) G X i ˆ x i s i − Γ N ! (33)plus o ( N ) terms. The expression can be decoupled through a Hubbard-Stratonovich transformation in terms of twoadditional fields µ ( t ) and ˜ µ ( t ) and taking the limit (cid:15) → − Z dt ˆ µ ( t ) X i s i − N ! − Z dtµ ( t ) X i ˆ x i s i − Γ N ! (34)Introducing the variable µ (1) = µ ( t ) + ˆ µ ( t ) η we can write the above term as − Z d µ (1) X i s (1) − N ! + 12 Γ Z dtµ ( t ) (35)The second term is essential to get the correct saddle-point equations and the correct value of the rate. D. The p -spin Spherical Model: The Action and Saddle-Point Equations The Hamiltonian of the spherical model is given by: H = ∞ X p =1 X i <...
12 Tr ln( Γ + µ + Λ ) + const. + 12 Z d a µ ( a ) + n X a =1 Γ Z µ a ( t ) dt (38)where f ( x ) ≡ ∞ X p =1 µ p x p (39)The expression const. above collects a number of terms coming from the explicit Gaussian integration, it is divergentand cancels the divergences associated to the expression Tr ln( Γ + µ + Λ ). These are pathologies of the path integralrepresentation that are fixed going back to the discrete times, we will further discuss them later. Note that theapplication of the standard manipulations to the form (6) would lead to the above expression without the last term,the last term appears instead if we want to use the simplified formulation which is suitable for multi-linear interactions.The above expression has to be extremized with respect to Q ( ab ), Λ ( ab ) and with respect to µ ( a ) . Extremizationwith respect to Q ( ab ), Λ( ab ) leads to the saddle point equations: Λ ( ab ) = − β f ( Q ( ab )) (40) Q = 1 Γ + Λ + µ (41)The last equation can be rewritten as: Γ( ac ) Q ( cb ) + µ ( a ) δ ( ac ) Q ( cb ) + Λ( ac ) Q ( cb ) = δ ( ab ) (42)where there integration over the variable c is implicit. The above expression is extremely compact, in the followingwe wil see that it encodes eight integro-differential equations.From now on we specialize to the case of a 1RSB transition. We start noticing that the initial and final configurationare weigthed with the equilibrium Gibbs measure and their properties must not depend on the dynamics. This isgranted by the fact that the terms depending on the dynamics in the equations for the replicas of the initial and finalconfigurations are O ( n ) and disappear in the n → Q ( ab ) with a and b corresponding to the a, b = 1 , . . . , m + m replicas associated to the equilibrium boundary conditions is theordinary equilibrium replica matrix. For simplicity we will work in zero field and zero random field: f (0) = 0 , (43)in this way the overlap between different equilibrium configurations is zero. Therefore above T d the solution is Q ( ab ) = δ ( ab ) that plugged into the above equations leads to: µ = 1 + β f (1) . (44)Below T d the solution is actually 1RSB. The replicas determining the boundary conditions should naturally belongto different RSB blocks ensuring that we are studying the transition between different states. On the other hand theequations for the dynamical part will have a non-vanishing correlation with the remaining x − x − x − T s < T < T d we haveexactly x = 1 and thus this contribution vanishes. Therefore we will safely use the high temperature RS solution Q ( ab ) = δ ( ab ) also in the region T s < T < T d . This is consistent with the known result that annealed and quenchedaverages are equivalent above T s .The global order parameter Q ( ab ) can be divided into a static part corresponding to the m + m replicas, a dynamicpart describing the the n replicas of the dynamics and a mixed part. We have seen before that, as it should, the staticpart does not depend on the dynamics part due to the n → n dynamical replicas, therefore the dynamicalcomponent of Q ( ab ) will be characterized by two matrices Q ( ab ) = δ αβ Q ( ab ) + (1 − δ αβ ) Q dt ( ab ) (45)where we have moved from the full coordinates ( ab ) to purely dynamical coordinates ( ab ) and replica coordinates α, β = 1 , . . . , n . The RS ansatz also implies µ α ( t ) = µ ( t ) , ˆ µ α ( t ) = ˆ µ ( t ) . (46)We also have: Γ( ab ) = δ αβ Γ( ab ) (47) δ ( ab ) = δ αβ δ ( ab ) . (48)The dynamical components of equations (42) can then be rewritten in a form that can be analitically continued toreal values of the replica number n Γ( ac ) Q ( cb ) + βµ ( a ) δ ( ac ) Q ( cb ) + (Λ Q ) st = δ ( ab ) (49)Γ( ac ) Q dt ( cb ) + βµ ( a ) δ ( ac ) Q dt ( cb ) + (Λ Q ) dt = 0 (50)where again the integration over the variable c is implicit and(Λ Q ) st = Λ( ac ) Q ( cb ) + ( n − dt ( ac ) Q dt ( cb ) + Λ( a Q (1 b ) + Λ( a Q (2 b ) (51)(Λ Q ) dt = Λ dt ( ac ) Q ( cb ) + Λ( ac ) Q dt ( cb ) + ( n − dt ( ac ) Q dt ( cb ) + Λ( a Q (1 b ) + Λ( a Q (2 b ) (52)where 1 and 2 label the final and initial conditions, note that the corresponding terms appear when we integrate overthe full coordinate c in eq. (42). In the above equations we have naturally:Λ( ab ) = − β f ( Q ( ab )) (53)Λ dt ( ab ) = − β f ( Q dt ( ab )) . (54) Q ( ab ) and Q dt ( ab ) can be expressed in terms of four two-time functions as: Q ( ab ) ≡ C ( t a , t b ) + ˆ R ( t a , t b ) η a + ˆ R ( t a , t b ) η b + ˆ X ( t a , t b ) η a η b (55) Q dt ( ab ) ≡ C dt ( t a , t b ) + ˆ R dt ( t a , t b ) η a + ˆ R dt ( t a , t b ) η b + ˆ X dt ( t a , t b ) η a η b (56)The same representation can be obtained from any function A ( ab ) as A ( ab ) ≡ C A ( t a , t b ) + ˆ R ,A ( t a , t b ) η a + ˆ R ,A ( t a , t b ) η b + ˆ X A ( t a , t b ) η a η b . (57)from which we obtain C Λ ( t a , t b ) = − β f ( C ( t a , t b )) (58)ˆ R , Λ ( t a , t b ) = − β f ( C ( t a , t b )) ˆ R ( t a , t b ) (59)ˆ R , Λ ( t a , t b ) = − β f ( C ( t a , t b )) ˆ R ( t a , t b ) (60)ˆ X Λ ( t a , t b ) = − β f ( C ( t a , t b )) ˆ X ( t a , t b ) − β f ( C ( t a , t b )) ˆ R ( t a , t b ) ˆ R ( t a , t b ) (61) C Λ dt ( t a , t b ) = − β f ( C dt ( t a , t b )) (62)ˆ R , Λ dt ( t a , t b ) = − β f ( C dt ( t a , t b )) ˆ R dt ( t a , t b ) (63)ˆ R , Λ dt ( t a , t b ) = − β f ( C dt ( t a , t b )) ˆ R dt ( t a , t b ) (64)ˆ X Λ dt ( t a , t b ) = − β f ( C dt ( t a , t b )) ˆ X dt ( t a , t b ) − β f ( C dt ( t a , t b )) ˆ R dt ( t a , t b ) ˆ R dt ( t a , t b ) (65)The operator Γ( ab ) is defined from12 Z dadb Γ( ab ) q ( a ) q ( b ) ≡ Z (cid:18) ˙ q − ˆ x (cid:19) dt , (66)that leads to Γ( ab ) ≡ (cid:20) − δ ( t a − t b ) η a η b − δ ( t a − t b ) (cid:21) . (67)For Γ = 1 we have:Γ( ac ) Q ( cb ) = − R ( t a , t b ) − d dt a C ( t a , t b ) η a − X ( t a , t b ) η b − d dt a ˆ R ( t a , t b ) η a η b . (68)The corresponding expression in eq. (50) is obtained replacing Q ( ab ) with Q dt ( ab ). The term depending on µ leadsto: µ ( a ) Q ( ab ) = µ ( t a ) C ( t a , t b ) ++ η a (cid:16) µ ( t a ) ˆ R ( t a , t b ) + ˆ µ ( t a ) C ( t a , t b ) (cid:17) ++ η b µ ( t a ) ˆ R ( t a , t b ) ++ η a η b (cid:16) ˆ µ ( t a ) ˆ R ( t a , t b ) + µ ( t a ) ˆ X ( t a , t b ) (cid:17) . (69)As above, the corresponding expression in eq. (50) is obtained replacing Q ( ab ) with Q dt ( ab ). In order to completethe derivation of the saddle-point equations we need the expression for a product of the form A ( ac ) B ( cb ): A ( ac ) B ( cb ) = C A ( t a , t c ) ˆ R ,B ( t c , t b ) + ˆ R ,A ( t a , t c ) C B ( t c , t b ) ++ η a (cid:16) ˆ R ,A ( t a , t c ) ˆ R ,B ( t c , t b ) + ˆ X A ( t a , t c ) C B ( t c , t b ) (cid:17) ++ η b (cid:16) ˆ R ,A ( t a , t c ) ˆ R ,B ( t c , t b ) + C A ( t a , t c ) ˆ X B ( t c , t b ) (cid:17) ++ η a η b (cid:16) ˆ R ,A ( t a , t c ) ˆ X B ( t c , t b ) + ˆ X A ( t a , t c ) ˆ R ,B ( t c , t b ) (cid:17) . (70)For the contributions of the initial and final configurations in the interaction term we have:Λ( a Q (1 b ) + Λ( a Q (2 b ) = C Λ ( t a , −T ) C ( −T , t b ) + C Λ ( t a , t f ) C ( T , t b )+ η a (cid:16) ˆ R , Λ ( t a , −T ) C ( −T , t b ) + ˆ R , Λ ( t a , T ) C ( T , t b ) (cid:17) ++ η b (cid:16) C Λ ( t a , −T ) ˆ R ( −T , t b ) + C Λ ( t a , T ) ˆ R ( T , t b ) (cid:17) ++ η a η b (cid:16) ˆ R , Λ ( t a , −T ) ˆ R ( −T , t b ) + ˆ R , Λ ( t a , T ) ˆ R ( T , t b ) (cid:17) . (71)We also have: δ ( ab ) = ( η a + η b ) δ ( t a − t b ) . (72)Collecting the various components in eq. (49) and (50) we obtain eight integro-differential equations that will bestudied in the next section. The extremization of expression (38) with respect to µ ( t ) and ˆ µ ( t ) gives the followingconditions: C ( t, t, ) = 1 , ˆ R ( t, t ) = ˆ R ( t, t ) = 12 . (73)note that the last term in (38) is essential to obtain the correct expression for ˆ R ( t, t ) and ˆ R ( t, t )
1. Integrals
It is useful to express the sum of the elements of a generic object B ( ab ) in terms of its components. We considerthe case in which it has essentially the structure of the order parameter Q ( ab ) . In particular we have for the m + m static replicas: B ( α, β ) = δ αβ C B ( T , T ) (74)similarly the n dynamical replicas are correlated only with the static replicas associated to the initial and finalconditions and their correlations have a RS form. We first consider A ( a ) ≡ Z B ( ba ) db (75)For a corresponding to any of the m + m static replicas other than those fixing the initial and final condition wesimply have A ( a ) = C B ( T , T ) (76)For a corresponding to either one of the two static replicas controlling the initial and final conditions we have: A ( T ) = Z B ( T , b ) db = C B ( T , T ) + C B ( T , −T ) + n Z ˆ R ,B ( T , t ) dt (77)For a given by the η component of one of the dynamical replicas we have:ˆ a ( t ) = ˆ R ,B ( t, T ) + ˆ R ,B ( t, −T ) + Z ˆ X B ( t, t ) dt + ( n − Z ˆ X dtB ( t, t ) dt . (78)Putting everything together we have: Z B ( ab ) dadb = Z A ( a ) da = ( m + m ) C B ( T , T ) ++ n (cid:18)Z ˆ R ,B ( t, −T ) dt + Z ˆ R ,B ( t, T ) dt + Z ˆ R ,B ( T , t ) dt + Z ˆ R ,B ( −T , t ) dt ++ Z ˆ X B ( t, t ) dt dt + ( n − Z ˆ X B dt ( t, t ) dt dt (cid:19) (79) II. ANALYSIS OF THE SADDLE-POINT EQUATIONS OF THE SPHERICAL p -SPIN-GLASS MODELA. The order parameter and its meaning In the previous section we have derived saddle-point equations for the computation of the log-rate. The orderparameter is a couple of 2 × t and t on the square −T ≤ t, t ≤ T : C ≡ (cid:18) C ( t, t ) ˆ R ( t, t )ˆ R ( t, t ) ˆ X ( t, t ) (cid:19) (80) C dt ≡ (cid:18) C dt ( t, t ) ˆ R dt ( t, t )ˆ R dt ( t, t ) ˆ X dt ( t, t ) (cid:19) (81)As mentioned in the main text he physical meaning of the order parameter is straightforward. We have C ( t, t ) = [ h s i ( t ) s i ( t ) i ] (82) C dt ( t, t ) = [ h s i ( t ) ih s i ( t ) i ] (83)Where the square brackets represent the averages with respect to dynamical trajectories connecting configurations σ and τ , the square bracket represent average with respect to σ and τ and the overline represents the average over thequenched disorder (if present). Thus C ( t, t ) is the correlation on the same trajectory and therefore C ( t, t ) = 1 atall times in the spherical and Ising model. C dt ( t, t ) measures the correlations between the configurations visited bydifferent trajectories (hence the suffix dt ). Since by definition all trajectories have the same initial and final conditionwe have C dt ( −T , t ) = C ( −T , t ) (84) C dt ( T , t ) = C ( T , t ) (85)Furthermore the correlations are symmetric with respect to ( t, t ) → ( t , t ): C ( t, t ) = C ( t , t ) (86) C dt ( t, t ) = C dt ( t , t ) (87)Given that the measure over the trajectories is invariant under time reversal we have an additional symmetry withrespect to the exchange of the initial and final configuration. Given that t in = − t fin this symmetry translate into: C ( t, t ) = C ( − t, − t ) (88) C dt ( t, t ) = C dt ( − t, − t ) (89)For the ˆ R components of the order parameter we have:ˆ R ( t, t ) = [ h ˆ x i ( t ) s i ( t ) i ] (90)ˆ R dt ( t, t ) = [ h ˆ x i ( t ) ih s i ( t ) i ] (91)where ˆ x is the auxiliary variable of the dynamics. They translate into:ˆ R ( t, t ) = (cid:20) δ h s i ( t ) i β δh ( t ) (cid:21) − ddt C ( t, t ) (92)ˆ R dt ( t, t ) = (cid:20) δ ln D ( σ, τ ) β δh ( t ) h s i ( t ) i (cid:21) − ddt C dt ( t, t ) (93)Thus ˆ R ( t, t ) is connected to the response of the measure over trajectories to a field h ( t ). The functions ˆ R ( t, t )and ˆ R dt ( t, t ) are equal to the l.h.s.’s of the above equations with the exchange t ↔ t . Thus while neither function issymmetric with respect to t ↔ t we have: ˆ R ( t, t ) = ˆ R ( t , t ) (94)ˆ R dt ( t, t ) = ˆ R dt ( t , t ) (95)0which implies that the matrices (80) and (81) are symmetric. On the other hand time-reversal invariance impliesthat ˆ R ( t, t ), ˆ R ( t, t ), ˆ R dt ( t, t ), ˆ R dt ( t, t ), are symmetric with respect to ( t, t ) → ( − t, − t ) (because t in = − t fin asabove). The symmetries implies that we effectively have two functions ˆ R ( t, t ) and ˆ R dt ( t, t ) such thatˆ R ( t, t ) = ˆ R ( t , t ) = ˆ R ( t , t ) (96)ˆ R dt ( t, t ) = ˆ R dt ( t , t ) = ˆ R dt ( t , t ) (97)For the ˆ X components we have: ˆ X ( t, t ) = [ h ˆ x i ( t )ˆ x i ( t ) i ] (98)ˆ X dt ( t, t ) = [ h ˆ x i ( t ) ih ˆ x i ( t ) i ] (99)The physical meaning can also be associated to particular responses. The above formulas imply that both ˆ X ( t, t )and ˆ X dt ( t, t ) are symmetric with respect to t ↔ t and to time-reversal ( t, t ) → ( − t, − t ). B. The equations in compact form
In order to write the saddle-point equations in a form suitable for numerical integration we consider the space of2 × t and t on the square −T ≤ t, t ≤ T . The genericelement of this space can be written as A ≡ (cid:18) C A ( t, t ) ˆ R ,A ( t, t )ˆ R ,A ( t, t ) ˆ X A ( t, t ) (cid:19) (100)Given two elements A and B in the above space we have a natural definition of the product that generalizes thematrix product (it corresponds to exactly to ordinary matrix products if times are discretized). For a real function B ( x ) we also define the element-wise function B [ A ] according to: B [ A ] ≡ (cid:18) B ( C A ( t, t )) B ( C A ( t, t )) ˆ R ,A ( t, t ) B ( C A ( t, t )) ˆ R ,A ( t, t ) B ( C A ( t, t )) ˆ X A ( t, t ) + B ( C A ( t, t )) ˆ R ,A ( t, t ) ˆ R ,A ( t, t ) (cid:19) . (101)We define also: M ≡ (cid:18) − δ ( t, t ) + ˆ µ ( t ) δ ( t, t ) µ ( t ) δ ( t, t ) µ ( t ) δ ( t, t ) − δ ( t, t ) (cid:19) (102)and T ≡ (cid:18) δ ( t, t ) δ ( t, t ) o (cid:19) . (103)In order to write down the saddle-point equations it is useful to introduce two additional objects Λ and Λ dt that arealso 2 × δ ∓ ≡ (cid:18) ˆ R Λ , ( t, ∓T ) C ( ∓T , t ) ˆ R Λ , ( t, ∓T ) ˆ R ( ∓T , t ) C Λ ( t, ∓T ) C ( ∓T , t ) C Λ ( t, ∓T ) ˆ R ( ∓T , t ) (cid:19) . (104)With the above definitions the saddle-point equations of the spherical model derived in the previous section read:Λ = − β f [ C ] , Λ dt = − β f [ C dt ] (105) M C + T Λ T C + ( n − T Λ dt T C dt + δ − + δ + = I (106) M C dt + T Λ dt T C + T Λ T C dt + ( n − T Λ dt T C dt + δ − + δ + = 0 (107)1The quantity µ ( t ) and ˆ µ ( t ) must be determined self-consistently in order to satisfy the conditions: C ( t, t ) = 1 , ˆ R ( t, t ) = 12 ∀ t (108)Due to the presence of the operator M that contains second-order derivatives the first equation must be also supple-mented with the following boundary conditions on the components of C : C ( −T , t ) = C ( t , −T ) (109) C ( T , t ) = C ( t , T ) (110)ˆ R ( −T , t ) = ˆ R ( t , −T ) (111)ˆ R ( T , t ) = ˆ R ( t , T ) (112)Identical boundary conditions must be enforced for the corresponding components of C dt to complete the definitionof the second equation. C dt ( −T , t ) = C dt ( t , −T ) (113) C dt ( T , t ) = C dt ( t , T ) (114)ˆ R dt ( −T , t ) = ˆ R dt ( t , −T ) (115)ˆ R dt ( T , t ) = ˆ R dt ( t , T ) (116)Additional boundary conditions that are a consequence of the properties of the initial and final configurations are: C ( −T , −T ) = C ( T , T ) = 1 , C ( −T , T ) = C ( T , −T ) = 0 (117) C dt ( −T , −T ) = C dt ( T , T ) = 1 , C dt ( −T , T ) = C dt ( T , −T ) = 0 (118)The above boundary conditions are necessary to compute the r.h.s. of the saddle point equations (106,107) for ageneric C and C dt . But given the structure of the equation and the physical meaning of the order parameter thefollowing symmetries must also be obeyed by the solution: C ( t, t ) = C ( t , t ) , ˆ X ( t, t ) = ˆ X ( t , t ) , ˆ R ( t, t ) = ˆ R ( t , t ) (119) C dt ( t, t ) = C dt ( t , t ) , ˆ X dt ( t, t ) = ˆ X dt ( t , t ) , ˆ R dt ( t, t ) = ˆ R dt ( t , t ) (120)The symmetry of the problem under the exchange between the initial and final configuration leads to the additionalsymmetries that in the case t in = − t fin make all components of C and C dt symmetric under the transformation( t, t ) → ( − t, − t ). Similarly we have: µ ( t ) = µ ( − t ) ˆ µ ( t ) = ˆ µ ( − t ) (121)Finally since all replicas of the dynamics have the same initial conditions the corresponding components of C and C dt are identical on the sides of the square −T ≤ t ≤ T , −T ≤ t ≤ T . We have fore instance C ( −T , t ) = C dt ( −T , t ).In a numerical treatment times are discretized and C and C dt become actual matrices, the above symmetries allow afour-fold reduction of the memory required to store them. C. The equations in expanded form
Expressions (106) and (107) are compact and useful for a numerical treatment. They correspond to eight integro-differential equations that we write in the following in explicit form for completeness. They have to be supplementedwith the definitions (58-65) and the boundary conditions of the previous subsection.0 = − R ( t, t ) + µ ( t ) C ( t, t ) ++ Z + T−T (cid:16) C Λ ( t, t ) ˆ R ( t , t ) + ˆ R , Λ ( t, t ) C ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) C Λ dt ( t, t ) ˆ R dt ( t , t ) + ˆ R , Λ dt ( t, t ) C dt ( t , t ) (cid:17) dt ++ C Λ ( t, −T ) C ( −T , t ) + C Λ ( t, T ) C ( T , t ) . (122)2 δ ( t − t ) = − d dt C ( t, t ) + µ ( t ) ˆ R ( t, t ) + ˆ µ ( t ) C ( t, t ) ++ Z + T−T (cid:16) ˆ R , Λ ( t, t ) ˆ R ( t , t ) + ˆ X Λ ( t, t ) C ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ R dt ( t , t ) + ˆ X Λ dt ( t, t ) C dt ( t , t ) (cid:17) dt ++ ˆ R , Λ ( t, −T ) C ( −T , t ) + ˆ R , Λ ( t, T ) C ( T , t ) . (123) δ ( t − t ) = − X ( t, t ) + µ ( t ) ˆ R ( t, t ) ++ Z + T−T (cid:16) ˆ R , Λ ( t, t ) ˆ R ( t , t ) + C Λ ( t, t ) ˆ X ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ R dt ( t , t ) + C Λ dt ( t, t ) ˆ X dt ( t , t ) (cid:17) dt ++ C Λ ( t, −T ) ˆ R ( −T , t ) + C Λ ( t, T ) ˆ R ( T , t ) . (124)0 = − d dt ˆ R ( t, t ) + µ ( t ) ˆ X ( t, t ) + ˆ µ ( t ) ˆ R ( t, t ) ++ Z + T−T (cid:16) ˆ R , Λ ( t, t ) ˆ X ( t , t ) + ˆ X Λ ( t, t ) ˆ R ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ X dt ( t , t ) + ˆ X Λ dt ( t, t ) ˆ R dt ( t , t ) (cid:17) dt ++ ˆ R , Λ ( t, −T ) ˆ R ( −T , t ) + ˆ R , Λ ( t, T ) ˆ R ( T , t ) . (125)0 = − R dt ( t, t ) + µ ( t ) C dt ( t, t ) ++ Z + T−T (cid:16) C Λ dt ( t, t ) ˆ R ( t , t ) + ˆ R , Λ dt ( t, t ) C ( t , t ) (cid:17) dt ++ Z + T−T (cid:16) C Λ ( t, t ) ˆ R dt ( t , t ) + ˆ R , Λ ( t, t ) C dt ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) C Λ dt ( t, t ) ˆ R dt ( t , t ) + ˆ R , Λ dt ( t, t ) C dt ( t , t ) (cid:17) dt ++ C Λ ( t, −T ) C ( −T , t ) + C Λ ( t, T ) C ( T , t ) . (126)0 = − d dt C dt ( t, t ) + µ ( t ) ˆ R dt ( t, t ) + ˆ µ ( t ) C dt ( t, t ) ++ Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ R ( t , t ) + ˆ X Λ dt ( t, t ) C ( t , t ) (cid:17) dt ++ Z + T−T (cid:16) ˆ R , Λ ( t, t ) ˆ R dt ( t , t ) + ˆ X Λ ( t, t ) C dt ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ R dt ( t , t ) + ˆ X Λ dt ( t, t ) C dt ( t , t ) (cid:17) dt ++ ˆ R , Λ ( t, −T ) C ( −T , t ) + ˆ R , Λ ( t, T ) C ( T , t ) . (127)30 = − X dt ( t, t ) + µ ( t ) ˆ R dt ( t, t ) ++ Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ R ( t , t ) + C Λ dt ( t, t ) ˆ X ( t , t ) (cid:17) dt ++ Z + T−T (cid:16) ˆ R , Λ ( t, t ) ˆ R dt ( t , t ) + C Λ ( t, t ) ˆ X dt ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ R dt ( t , t ) + C Λ dt ( t, t ) ˆ X dt ( t , t ) (cid:17) dt ++ C Λ ( t, −T ) ˆ R ( −T , t ) + C Λ ( t, T ) ˆ R ( T , t ) . (128)0 = − d dt ˆ R dt ( t, t ) + µ ( t ) ˆ X dt ( t, t ) + ˆ µ ( t ) ˆ R dt ( t, t ) ++ Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ X ( t , t ) + ˆ X Λ dt ( t, t ) ˆ R ( t , t ) (cid:17) dt ++ Z + T−T (cid:16) ˆ R , Λ ( t, t ) ˆ X dt ( t , t ) + ˆ X Λ ( t, t ) ˆ R dt ( t , t ) (cid:17) dt ++ ( n − Z + T−T (cid:16) ˆ R , Λ dt ( t, t ) ˆ X dt ( t , t ) + ˆ X Λ dt ( t, t ) ˆ R dt ( t , t ) (cid:17) dt ++ ˆ R , Λ ( t, −T ) ˆ R ( −T , t ) + ˆ R , Λ ( t, T ) ˆ R ( T , t ) . (129) D. The solutions on the corners and on the diagonal
In the following we discuss a number of properties that follows from the equations. The equation for the η b component implies that ˆ X ( t, t ) = − δ ( t − t ) + δ ˆ X ( t, t ) (130)where δ ˆ X ( t, t ) is a bounded function. In the limit t, t → ±T the interaction part in the equations greatly simplifiesbecause we have C ( −T , t ) = C dt ( −T , t ) , ˆ R ( t, −T ) = ˆ R dt ( t, −T ) , δ ˆ X ( −T , t ) = ˆ X dt ( −T , t ) . (131)This allows to characterize the order parameters on the corners of the t, t domain. In particular one easily sees thatfor n = 0 the interactions terms for the dynamical replicas cancel in the limit t, t → −T and only δ in remains. Theequations for the scalar component then reads: − µ ( −T ) − β f (1) = 0 → lim t →±T µ ( t ) = µ eq (132)which implies that µ ( t ) goes continuously to the equilibrium value as t → ±T . The equation for the η b componentleads to: δ ˆ X ( ±T , ±T ) = ˆ X dt ( ±T , ±T ) = µ eq β f (1) (133)Similarly for n = 0 all interaction terms cancels in the limit t → −T , t → T and the equations imply:0 = C ( ±T , ∓T ) = C dt ( ±T , ∓T ) (134)0 = ˆ R ( ±T , ∓T ) = ˆ R dt ( ±T , ∓T ) (135)0 = ˆ R ( ±T , ∓T ) = ˆ R dt ( ±T , ∓T ) (136)0 = δ ˆ X ( ±T , ∓T ) = ˆ X dt ( ±T , ∓T ) (137)4 E. Derivatives of the log-rate
The expression of the logarithm rate in terms of the order parameter obtained in the path integral formulation isdivergent as will be discussed in the free case. Its computation requires to go back to the discretized case. Howevera convenient alternative is to compute the derivative of the log rate with respect to the temperature and integrateusing the knowledge of the infinite temperature limit that can be obtained explicitly. The partial derivative withrespect to the temperature of expression (38) is simply given by β ddn R d a d b f ( Q ( a , b )) and it coincides with thetotal derivative when computed on the solution of the saddle-point equations. Using the formulas for the integralscomputed in subsection I D 1 we obtain:1 N ∂ [ln ˆ T ] ∂β = β Z + T−T ˆ R ,f [ C ] ( t, −T ) dt + Z + T−T Z + T−T ˆ X f [ C ] ( t, t ) dtdt + ( n − Z + T−T Z + T−T ˆ X f [ C dt ] ( t, t ) dtdt ! . (138)As we mentioned in the main text it is interesting to consider a generalized [ln ˆ T ] in which the final configuration attime T is selected with the Gibbs weight corresponding to a different temperature β . Following the same argumentabove we have that the total derivative of the log-rate with respect β coincide with the partial derivative. One easilyfinds: 1 N ∂ [ln ˆ T ] ∂β = β Z + T−T ˆ R ,f [ C ] ( t, T ) dt (139)The probability of jumping to the β configurations with typical rate isexp h βE ( β ) / T ] − βE ( β ) / S ( β ) i (140)and the derivative of the logarithm of the above expression with respect to β for β = β is: β Z + T−T ˆ R ,f [ C ] ( t, T ) dt − β f (1) . (141) F. The energy
An important observable is the instantaneous energy on the trajectory: E ( t ) ≡ ∞ X p =1 X i <...
1. The Saddle-Point equations and their solutions
In the infinite temperature limit the interaction vanishes and the system perform a free Brownian motion on the N − T ] by integrating with respect to the temperature.The saddle point equations (106),107) become in this case: M C = I , M C dt = 0 (146)Explicitly the first equation read: − R ( t, t ) + µ ( t ) C ( t, t ) = 0 (147) − d dt C ( t, t ) + µ ( t ) ˆ R ( t, t ) + ˆ µ ( t ) C ( t, t ) = δ ( t − t ) (148) − X ( t, t ) + µ ( t ) ˆ R ( t, t ) = δ ( t − t ) (149) − d dt ˆ R ( t, t ) + µ ( t ) ˆ X ( t, t ) + ˆ µ ( t ) ˆ R ( t, t ) = 0 (150)Note that the equations for C does not depend on C dt . We start noticing the the conditions: C ( t, t ) = 1 , ˆ R ( t, t ) = 12 ∀ t , (151)plugged into (147) lead to µ ( t ) = 1 ∀ t , (152)then eq. (147), the symmetries ˆ R ( t, t ) = ˆ R ( t , t ), C ( t, t ) = C ( t , t ) and eq. (149) lead toˆ R ( t, t ) = 12 C ( t, t ) (153)ˆ R ( t, t ) = 12 C ( t, t ) (154)ˆ X ( t, t ) = − δ ( t − t ) + 14 C ( t, t ) . (155)As a consequence the correlation C ( t, t ) is determined by the following equation: − d dt C ( t, t ) + (cid:18)
12 + ˆ µ ( t ) (cid:19) C ( t, t ) = δ ( t − t ) (156)To be solved with boundary conditions C ( ±T , t ) = C ( t, ±T ) and C ( ±T , ±T ) = 1, C ( ±T , ∓T ) = 0. The equationhas the solution C ( t, t ) = C ( | t − t | ) , ˆ µ ( t ) = ˆ µ . (157)To determine the function C ( x ) it is convenient to define a quantity −∞ < a < T = 12 √ a arctanh √ a (158)where T is half the time difference between the initial and final times. We then have, see fig (1):ˆ µ = a −
12 (159) C ( x ) = cosh √ ax − √ a sinh √ ax . (160)6 FIG. 1: Left: The function C ( t, t ) in the free case for a = 0 .
999 corresponding to T = 2 . C dt ( t, t ). Note that we have; C (0) = 1 , ˙ C (0) = − . (161)This is a general result that remains true also at finite temperature because it is a consequence of the presence of thedelta function and of the symmetry C ( t, t ) = C ( t , t ). It is interesting to discuss the limits of large and small T . Forlarge values of T we have: a ≈ − e − T , ˆ µ ≈ − e − T for T → ∞ (162)Note that in the large T limit ˆ µ tends to zero and the function C ( x ) tends to the equilibrium solution C eq ( x ) = e − x (163)On the other hand a becomes negative for T < / a ≈ (cid:18) T − (cid:19) for T ≈
12 (164)Note that a / is imaginary and the hyperbolic functions in (160) become ordinary trigonometric functions: C ( x ) = cos p | a | x − p | a | sin p | a | x (165)In the small T limit a tends to minus infinity as a ≈ − π T , ˆ µ ≈ − π T for T → C ( x ) ≈ cos πx T for T → x = 2 T but it does not display the correct behavior C ( x ) ≈ − x at small x . Indeed it is only valid for x = O ( T ) and the correct linear behavior is recovered for x (cid:28) T . We now turn to theequations for C dt that read: − R ( t, t ) + µ ( t ) C dt ( t, t ) = 0 (168) − d dt C dt ( t, t ) + µ ( t ) ˆ R dt ( t, t ) + ˆ µ ( t ) C dt ( t, t ) = 0 (169) − X dt ( t, t ) + µ ( t ) ˆ R dt ( t, t ) = 0 (170) − d dt ˆ R dt ( t, t ) + µ ( t ) ˆ X dt ( t, t ) + ˆ µ ( t ) ˆ R dt ( t, t ) = 0 (171)7The condition µ ( t ) = 1 derived earlier plus eq. (168), the symmetries ˆ R dt ( t, t ) = ˆ R dt ( t , t ), C dt ( t, t ) = C dt ( t , t ) andeq. (170) lead to: ˆ R dt ( t, t ) = 12 C dt ( t, t ) (172)ˆ R dt ( t, t ) = 12 C dt ( t, t ) (173)ˆ X dt ( t, t ) = 14 C dt ( t, t ) (174)and C dt ( t, t ) obeys the equation − d dt C dt ( t, t ) + (cid:18)
12 + ˆ µ (cid:19) C dt ( t, t ) = 0 (175)To be solved with boundary conditions C dt ( ±T , t ) = C dt ( t, ±T ) and C dt ( ±T , ±T ) = 1, C dt ( ±T , ∓T ) = 0.The solution reads: C dt ( t, t ) = √ − aa cosh( √ a ( t + t )) − − aa cosh( √ a ( t − t )) T > / C dt ( t, t ) = √ − aa cos( p | a | ( t + t )) − − aa cos( p | a | ( t − t )) 0 < T < / . (177)In the limit T → ∞ we have C dt ( t, t ) = 0 for any finite t, t . Close to the initial and finite time we have instead C dt ( ±T + δt, ±T + δt ) = C eq ( | δt + δt | ) for T → ∞ . (178)In the limit T → C dt ( t, t ) ≈ cos π ( t − t )4 T for T → C dt ( t, t ) ≈ C ( t, t ) meaning that all trajectories tend to follow the same path.
2. The Transition Rate
In the free case the expression for the Replicated logarithm transition rate (38) simplifies considerably due to β = 0and Λ( ab ) = 0. In the RS ansatz we have: n (cid:18) Z dt ˆ µ ( t ) + Γ Z µ ( t ) dt + Z dσ √ π e − σ Z dτ √ π e − τ ln Z ( σ, τ ) (cid:19) (180)where Z ( σ, τ ) = Z q ( T )= τq ( −T )= σ [ dq ] exp (cid:20)Z dt (cid:18) − ˙ q − Γ µ q −
12 ˆ µq (cid:19)(cid:21) (181)We recognize the path integral representation of the Harmonic oscillator that is usually written as Z ( σ, τ ) = Z q ( T )= τq ( −T )= σ [ dq ] exp (cid:20) − Z dt (cid:18) m ˙ q + 12 mω q (cid:19)(cid:21) (182)with the identification m = 1 / (2Γ ) that leads to the same ˙ q factor and the same [ dq ]. As discussed in classictextbooks the above path integral is ill defined. This is easily seen switching to a frequency representation where itis ultraviolet divergent as ln Z ∝ R ∞ dk ln( ω + k ). The actual quantity Z ( σ, τ ) is finite because the differential [ dq ]includes a prefactor diverging as 1 / ∆ t as we have seen in section (I B). A careful computation leads to the followingexpression (see eq. 2.23 in Zinn-Justin’s [1], eq. 13.45 in Parisi’s [2] ):ln Z ( σ, τ ) = 12 ln mω π sinh ωτ − mω ωτ (cid:2) cosh ωτ ( σ + τ ) − στ (cid:3) (183)8with m = 12 , ω = √ a , τ ≡ = t fin − t in = 2 T (184)where we have fixed Γ = 1 and used the previous results µ = 1 and a ≡ µ . Performing the averages over σ and τ we obtain: 1 N [ln ˆ T ] = −
12 ln 2 π −
12 + (1 + ˆ µ ) T + 14 ln − ˆ µ T goes to zero:1 N [ln ˆ T ] ≈ − π T , T ≈ . (186)For T going to infinity we have ˆ µ ≈ − e − T but the O ( T ) divergences in the last two terms cancel and the rate hasa finite limit. This limit is exactly equal to minus the entropy S ( β ) of the spherical model at infinite temperature( β = 0), ı.e. the logarithm of the surface of the N dimensional sphere of radius √ N :lim T →∞ N [ln ˆ T ] = −
12 ln 2 π −
12 = − S (0) . (187)This is expected in the T → ∞ limit at any finite N and implies that the two limits commute. One can show thatthe approach to the T → ∞ limit is exponential:1 N [ln ˆ T ] ≈ − S (0) − e − T , T (cid:29) . (188) H. The Ergodic Phase
As we mentioned in the main text in the ergodic phase defined by
T > T d the T → ∞ limit commutes with the N → ∞ limit. In order to characterize the solutions in this regime it is convenient to analyze first the ergodic limitat finite N . While the present formalism is fully invariant under time reversal to discuss this limit it is convenient toreintroduce the arrow of time. In the ergodic limit the dynamics looses any dependence on the initial configurationand the transition rate obeys lim T →∞ T T ( σ | τ ) = e − βH ( σ ) Z (189)that, according to the definition given in the main text, leads to:lim T →∞ ˆ T T ( σ, τ ) = 1 Z e − β H ( σ ) − β H ( τ ) (190)and the following exact result: lim T →∞ [ln ˆ T T ( σ, τ )] = − βE − ln Z = − S (191)In the ergodic limit we expect that trajectories close to the initial condition at time −T are not influenced by thefact that we are fixing the final configuration at time T . This implies that they are typical trajectories and sincethe initial configuration is weighted with the equilibrium weight we expect that correlations and response are thosevalid at equilibrium and in particular satisfy time-translational-invariance (TTI) and fluctuation-dissipation theorem(FDT). To see the implications on the functions it is convenient to start from the expressions of the six functionsas averages of s i ( t ) and ˆ x i ( t ) obtained in subsection II A. It is then convenient to transform back to the variable ˆ q introduced in section I B by writing ˆ x i = ˆ s i ( t ) − ˙ s i ( t ) . (192)It is well known the equilibrium averages of ˆ s i ( t ) are associated to responses to a field at time t , since according to(189) the distribution at σ is completely independent of what happens at any finite time t (cid:28) T we have for thesetimes h ˆ s ( t ) i = 0 , h ˆ s ( t )ˆ s ( t ) i = 0 (193)9The above relationships plus FDT allow to derive the following relationships expected in the ergodic limit: C ( t, t ) = C eq ( | t − t | ) (194)ˆ R ( t, t ) = ˆ R ( t , t ) = sign( t − t ) 12 ddt C ( t, t ) (195)ˆ X ( t, t ) = − d dt dt C ( t, t ) (196)Note that ˆ R ( t, t ) is symmetric (which is not true in general) and it is equal to 1 / X ( t, t ) is valid also for t = t and thus it may be integrated, thisis consistent with the fact that ˆ X ( t, t ) has a term proportional − δ ( t − t ) / t, t (cid:28) T , in order to study the region of times close to T it is better to reverse the arrow oftime, following the same arguments we obtain the above equations with ( t, t ) → ( − t, − t ) and since the functions aresymmetric with respect to this transformation we conclude that they are valid at all times.The correlation between different trajectories are different from zero only for times t, t that are both close to either −T or T . To determined them we can use the following relationship that follows from detailed balance: X τ e − βH ( τ ) Z T ∆ t ( τ | τ ) T ∆ t ( τ | τ ) = e − βH ( τ ) Z T ∆ t +∆ t ( τ | τ ) (197)Furthermore responses between different trajectories vanish because of h ˆ s ( t ) i = 0. In particular if we write times as t = ∓T ± ∆ t (∆ t ≥
0) we have: C dt ( t, t ) = C eq (∆ t + ∆ t ) (198)ˆ R dt ( t, t ) = ˆ R dt ( t , t ) = ∓ ddt C dt ( t, t ) (199)ˆ X dt ( t, t ) = 14 d dt dt C dt ( t, t ) . (200)Note that the expression for ˆ R dt ( t, t ) changes sign depending on weather we are close to ∓T , indeed the second case(+ T ) is obtained from the transformation ( t, t ) → ( − t, − t ) applied to the first case ( −T ). Note that a solution of thesaddle-point equations with the above structure can be found only if C eq ( ∞ ) = 0, thus the ergodic solution does notexist below T d where C eq ( ∞ ) = q in the thermodynamic limit and the condition C ( ±T , ∓T ) = 0 cannot be fulfilled.The above expressions are valid for both Ising and spherical models, in addition for the spherical model we have: µ ( t ) = µ eq , ˆ µ ( t ) = 0 . (201)Note that all components are parameterized by a single function C eq ( t ). In general one can introduce a genericequilibrium object A ( ab ) whose components A ( ab ) and A dt ( ab ) are obtained from the above formulas from a function C A ( t ). One can then show that the equilibrium structure is preserved by the application of a function, i.e. f ( A ( ab ) )has also the equilibrium structure with: f ( A ( ab ) ) → C f [ A ] ( t ) = f ( C A ( t )) (202)Similarly given another equilibrium object B ( ab ) parameterized by a function C B ( t ) one can show through a tediouscomputation that the product has also the equilibrium structure with: Z A ( ac ) B ( cb ) dc → C AB ( t ) = C A (0) C B ( t ) − Z t C A ( t − s ) dC B ( s ) ds ds (203)The above equation holds for C A ( ∞ ) = C B ( ∞ ) = 0 which is granted by the fact that we are working in zero fieldand f (0) = 0 so that the overlap between different equilibrium states is zero. The last two equations allow to deriveeasily the equilibrium equation: ˙ C ( t ) = − C ( t ) − β Z t f ( C ( t − s )) ˙ C ( s ) ds (204)where we have used µ = 1 + β f (1) that follows from the condition C (0) = − ˙ C (0) = 1. The above equations isusually written as [3–5]:˙ C ( t ) = − C ( t ) + β f ( C ( t ))(1 − C ( t )) − β Z t ( f ( C ( t − s )) − f ( C ( t ))) dCds ds (205)0Where C ( ∞ ) is the solution of the equation C = β f ( C )(1 − C ) (206)that admits the solution C = 0 at all temperature develops a non-zero solution at T d specified by the condition:1 = β d f ( C )(1 − C ) − f ( C )) (207)In the pure p -spin models f ( x ) = x p we have T d = s p ( p − p − p − p − , C d ( ∞ ) = p − p − . (208)Using the above equilibrium formulas one can show that the formula (79) for the integrals takes a very simple form: Z B ( ab ) dadb = ( m + m ) C B (0) + n (cid:18) Z R ,B ( −T , t ) dt + Z ˆ X B ( t, t ) dt dt + ( n − Z ˆ X B dt ( t, t ) dt dt (cid:19) == ( m + m + n ) C B (0) (209)We are now in position to show that for T > T d the limits commute and we havelim N →∞ S ( β ) N + lim T →∞ lim N →∞ N [ln ˆ T ] = 0 for 0 < β < β d (210)In the previous section we have shown explicitly that the relationship is satisfied for β = 0 and thus we can just showthat its derivative with respect to the β is also zero. As we discussed in sec. II E, the derivative of the rate withrespect to β is β dd n R d a d b f ( Q ( a , b )) that combined with (209) yields: ddβ lim T →∞ lim N →∞ N [ln ˆ T ] = β f (1) , (211)to be compared with 1 N dS ( β ) dβ = − β f (1) . (212)Thus we have lim T →∞ lim N →∞ N [ln ˆ T ] = − S ( β ) = − S (0) + β f (1) . (213)The above relationship is valid also for Ising systems. I. The Activated Phase
In the activated phase T s < T < T d the saddle-point equations have been solved numerically for the classic pure p -spin with f ( x ) = x where β d = 1 . β s = 1 . β = 1 .
695 for values of T = 8 , , , ,
40. In fig. (2) we clearly see the features discussed in the main text, inparticular: i) the overlap between the initial (final) configuration and the t = 0 intermediate configuration C ( ±T , T → ∞ limit, ii) the functions ˆ R ( t, t ), ˆ R dt ( t, t ), ˆ X ( t, t ), ˆ X dt ( t, t ) decrease withincreasing T for | t − t | = O ( T ).In fig. (2) the same data are rescaled to demonstrate the asymptotic limit T → ∞ discussed in the main text. Thesame-trajectory functions ˆ R ( t, t ) and ˆ X ( t, t ) deviate from this scaling in the region | t − t | = O (1) that goes to zero onthe scale of the plots for T → ∞ . In that region ˆ R ( t, t ) and ˆ X ( t, t ) converge to a finite limit asymptotically as shownin fig. (4). One can also check that in that region C ( t, t ), ˆ R ( t, t ) and ˆ X ( t, t ) satisfy the equilibrium relationships(195) and (196). The auxiliary functions µ ( t ) and ˆ µ ( t ) are shown in fig. (5). Note that the T − scaling of ˆ µ ( t ) is1consistent with the fact that C ( t, t ), ˆ R ( t, t ) and ˆ X ( t, t ) satisfy equilibrium relationship on the diagonal according to(201). However µ ( t ) changes with time implying equilibrium on finite time-scales but not on the global scale O ( T ).The asymptotic scaling discussed in the main text suggests that the leading corrections are O (1 / T ). In fig. (6) weplot the the data for µ ( t ) together with a 1 / T fit whose quality is excellent for T ≥
24. In fig. (7) we plot thederivative of h ln ˆ T i + S with respect to the inverse temperature computed according to expression (138). Accordingto the result of the previous section this quantity should go to zero in the T → ∞ limit for β < β d = 1 . T is larger than the equilibrium relaxation time that diverges as | T − T d | − γ where γ = 1 .
765 is a MCT exponent [4, 6, 7]. Thus at any T , no matter how large, there is always a rangeof temperatures ∆ T ∝ T − /γ such that for T − T d < ∆ T , h ln ˆ T i + S is smaller than zero. The effect decreases withincreasing T as the figure shows but is still significant at the values of T that we could study. On the other hand thecurves seems to have converged to − .
47 at β = 1 .
695 (as data from T = 40 also suggest) and this result, supplementedwith the information that it must be zero asymptotically at β d , allows for a rough estimate of the integral leading to h ln ˆ T i + S = − .
014 at β = 1 . III. NUMERICAL SOLUTION OF THE EQUATIONS
In the following we will add a few more details on the the algorithm developed to solve numerically the equationsas discussed in the main text. The equations in the form presented in section II B can be easily written in discretizedform, in particular the generic object (100) becomes a 2(2 N t + 1) × N t + 1) matrix and so does the operator M (102). For the second derivatives appearing in the equations we have used the formula f ( t ) ≈ f ( t + ∆ t ) + f ( t − ∆ t ) − f ( t )∆ t (214)where ∆ t = T /N t . (215)The above formula has a O (∆ t ) error if f ( t ) has continuous derivatives up to the third order. The integrals havebeen written as Z + T−T f ( t ) dt ≈ ∆ t f ( −T ) + f (+ T )) + N t − X i = − N t +1 f ( i ∆ t )∆ t . (216)The above trapezoidal rule also has a O (∆ t ) error for a continuous function. As mentioned in the main text, thefull algorithm has a O (∆ t ) error, however this is not trivial because the functions C ( t, t ), ˆ R ( t, t ) and ˆ X ( t, t ) havediscontinuous odd derivatives for t = t , while the discontinuity of the first derivative is canceled by the delta functionthe discontinuity of the third derivative leads to a O (∆ t ) error in expression (214) on the diagonal t = t . However onthe diagonal the order of the equation is not O (1) but O (1 / ∆ t ) due to the presence of the delta function in eq. (123)and of ˆ X ( t, t ) (that can be written as a delta function on the diagonal plus a regular part) in eq. (125), thus even ifthe absolute error on the diagonal is O (∆ t ), the relative error is O (∆ t ). In Fig. (8) we demonstrate that the erroris indeed O (∆ t ) for µ ( t ) and show how that the ∆ t → N t setavailable (up to N t = 2000 for T = 32 ,
40) using in most cases a fourth order form c + c ∆ t + c ∆ t + c ∆ t wherethe vanishing of the linear term ∆ t was imposed. One should note that the finite N t curves often display pathologiesdue to the discretisation that tend to be less severe increasing N t . For instance in the right of fig. (8) we see that µ ( t )at finite N t displays cusps close to t = ±T that are absent for N t → ∞ . It is impressing how a simple polynomialextrapolation over N t leads to the disappearance of these spurious features and allows to obtain accurate predictionswith relatively small values of N t that are indistinguishable from extrapolations obtained from considerably largervalues of N t . The same features are also seen in fig. (9), in this case corrections higher that O (∆ t ) are so small thattoo high-order interpolation functions overfit that data and it is convenient to use the form c + c ∆ t . Besides wesee that, even if the finite N t results for ˆ µ ( t ) have pronounced spurious cusps close to t ± T and may have the wrongsign, the extrapolations are cuspless and negative for all t .As we mentioned in the main text a moderate span of the eigenvalues of the Jacobian is essential for the successof Krylov methods. The equations in the form (49) and (50) are ill-conditioned because the operator M containssecond-order derivatives leading to a unbounded continuous spectrum. To overcome this problem we have consideredthe equations that one obtains multiplying (49) and (50) by M − whose Jacobian turns out indeed to have a discrete2and bounded spectrum as can be seen numerically using the fact that the Arnoldi diagonalization allows to obtainan approximate set of eigenvalues and eigenvectors. More details on the procedure can be found in the commentedcodes provided in the ancillary files section.From the computational point of view the procedure devised requires the computation of matrix products andelement-wise operations. One can solve eq. (49) and (50) for generic C and C dt and the algorithm converge tosolutions with the required symmetries discussed in sec. II A . However it is useful to work in the subspace ofsolutions with the required symmetries obtaining a four-fold reductions of the memory required to store C and C dt and thus the generic element of the Krylov subspace. Besides it turns out that this allows to consider a smallerKrylov subspace to obtain the same accuracy. Explicit use of the symmetry however requires an efficient procedure tocompress and decompress C and C dt and the elements of the Krylov subspace. The choice depends on the particularprogramming tool used, see the commented codes provided in the ancillary files section. [1] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena (Oxford Science Publications, 2002).[2] G. Parisi,
Statistical field theory (Addison-Wesley, 1988).[3] A. Crisanti and H.-J. Sommers, Zeitschrift f¨ur Physik B Condensed Matter , 341 (1992).[4] A. Crisanti, H. Horner, and H.-J. Sommers, Zeitschrift f¨ur Physik B Condensed Matter , 257 (1993).[5] T. Castellani and A. Cavagna, Journal of Statistical Mechanics: Theory and Experiment , P05012 (2005).[6] F. Caltagirone, U. Ferrari, L. Leuzzi, G. Parisi, F. Ricci-Tersenghi, and T. Rizzo, Physical review letters , 085702 (2012).[7] U. Ferrari, L. Leuzzi, G. Parisi, and T. Rizzo, Physical Review B , 014204 (2012). C ( - , - + Δ t ) vs. Δ t -
20 0 20 40 600.00.20.40.60.81.0 C ( - / ) vs. t + / - -
20 0 20 400.00.20.40.60.81.0 C ( ) vs. t C dt ( - ,t ) vs. t + -
20 0 20 40 600.00.10.20.30.40.5 C dt ( - / ) vs. t + / - -
20 0 20 400.000.050.100.150.200.250.300.35 C dt ( ) vs. t - R ( - ,t ) vs. t + -
20 0 20 40 600.0000.0050.0100.0150.0200.025 R ( - / ) vs. t + / - -
20 0 20 400.0000.0050.0100.0150.0200.025 R ( ) vs. t - R dt ( - ,t ) vs. t + -
20 0 20 40 60 - R dt ( - / ) vs. t + / - -
20 0 20 400.0000.0010.0020.0030.0040.0050.006 R dt ( ) vs. t X ( - ,t ) vs. t + -
20 0 20 40 600.00000.00050.00100.00150.00200.0025 X ( - / ) vs. t + / - -
20 0 20 400.00000.00050.00100.00150.00200.0025 X ( ) vs. t X dt ( - ,t ) vs. t + -
20 0 20 40 60 - X dt ( - / ) vs. t + / - -
20 0 20 400.000000.000020.000040.000060.00008 X dt ( ) vs. t FIG. 2: Spherical p -Spin Glass for β = 1 .
695 and p = 3. Solutions for T = 8 ( blue ), 16 ( yellow ) ,
24 ( green ) ,
32 ( red ) ,
40 ( purple )on various strips of the −T ≤ t, t ≤ T plane. The vertical scale of the ˆ R ( t, t ) and ˆ X ( t, t ) functions does not allows to see thediagonal region t − t = O (1). C ( - , - + Δ y ) vs. Δ y - C ( - / ) vs. t / + / - - C ( ) vs. t / C dt ( - ,t ) vs. t / + - C dt ( - / ) vs. t / + / - - C dt ( ) vs. t / - - R ( - ,t ) vs. t / + - R ( - / ) vs. t / + / - - R ( ) vs. t / - - R dt ( - ,t ) vs. t / + - - R dt ( - / ) vs. t / + / - - R dt ( ) vs. t / X ( - ,t ) vs. t / + - X ( - / ) vs. t / + / - - X ( ) vs. t / X dt ( - ,t ) vs. t / + - - X dt ( - / ) vs. t / + / - - X dt ( ) vs. t / - - C dt ( t,t ) vs. t / - - - R dt ( t,t ) vs. t / - - X dt ( t,t ) vs. t / FIG. 3: Spherical p -Spin Glass for β = 1 .
695 and p = 3. Scaled solutions for T = 8 ( blue ),16 ( yellow ) ,
24 ( green ) ,
32 ( red ) ,
40 ( purple ) on various strips of the −T ≤ t, t ≤ T plane. R ( - ,t ) vs. t + - - - R ( - / ) vs. t + / - - - R ( ) vs. t δ X ( - ,t ) vs. t + - - - δ X ( - / ) vs. t + / - - - δ X ( ) vs. t FIG. 4: Spherical p -Spin Glass for β = 1 .
695 and p = 3. The functions ˆ R ( t, t ) and δ ˆ X ( t, t ) = ˆ X ( t, t ) + δ ( t − t ) / T = 8 , , , ,
40 on various slices of the −T ≤ t, t ≤ T plane. The vertical and horizontal ranges are appropriate tovisualize the diagonal region t − t = O (1). On this scale the functions for different T ’s are almost indistinguishable. - -
20 0 20 404.95.05.15.25.3 μ ( t ) vs. t - -
20 0 20 40 - - - - μ ( t ) vs. t - - μ ( t ) vs. t / - - - - - - - - μ ( t ) vs. t / FIG. 5: Spherical p -Spin Glass for β = 1 .
695 and p = 3. The auxiliary functions µ ( t ) and ˆ µ ( t ) for T =8 ( blue ) ,
16 ( yellow ) ,
24 ( green ) ,
32 ( red ) ,
40 ( purple ). The T = 40 value in the rescaled plot for ˆ µ ( t ) has not been plot be-cause the extrapolation from N t = 2000 seems not accurate enough on that scale. - - μ ( t ) vs. t / μ ( ) vs. 1 / FIG. 6: Spherical p -Spin Glass for β = 1 .
695 and p = 3. Left: µ ( t ) for T = 8 ( yellow ),16 ( green ) ,
24 ( red ) ,
32 ( purple ) ,
40 ( brown ), the dashed line is an extrapolation to T = ∞ obtained from a linear fit in 1 / T .Right: The linear fit and the data for µ (0). = = - - - - N dd β ( Ln [ T ] + S ) vs. β FIG. 7: Spherical p -Spin Glass with p = 3. Derivative of the total transition probability with respect to the inverse temperaturefor T = 24 ,
32 and β = 1 .
55, 1 .
6, 1 .
65, 1 . fit2fit1 μ ( ) vs. Δ t - -
10 0 10 205.055.105.155.205.255.30 μ ( t ) vs. t FIG. 8: Spherical p -Spin Glass with p = 3, numerical solution for β = 1 . T = 24 and N t = 600, 800, 1000, 1200, 1400, 1600,1800. Left: numerical values of µ (0), the lines are two polynomial interpolation of the form c + c ∆ t + c ∆ t + c ∆ t overthe data for { } and { } . The polynomial extrapolations give respectively c = 5 . c = 5 . N t (individual points are not distinguishable at thescale of the plot). At the bottom there are two (indistinguishable) lines obtained from polynomial interpolations performedseparately for each time t , one of the form c + c ∆ t + c ∆ t + c ∆ t over the data for { } , and the otherof the form c + c ∆ t + c ∆ t over the data for { } . fit2fit1 μ ( ) vs. Δ t - -
10 0 10 20 - - μ ( t ) vs. t FIG. 9: Spherical p -Spin Glass with p = 3, numerical solution for β = 1 . T = 24 and various N t . Left: numerical values ofˆ µ (0) for N t = 600, 800, 1000, 1200, 1400, 1600, 1800, the lines are two polynomial interpolation of the form c + c ∆ t overthe data for { } and { } . Right: from top to bottom data for N t = 1000 (blue), 1400 (yellow), 1800 (green)(individual points are not distinguishable at the scale of the plot). At the bottom there are two (indistinguishable) lines (red,purple) obtained from polynomial interpolations of the form c + c ∆ t performed separately for each time t , using respectively N t = { } and N t = { }}