Dynamical transition in the D = 3 Edwards-Anderson spin glass in an external magnetic field
Janus Collaboration, M. Baity-Jesi, R. Alvarez Baños, A. Cruz, L.A. Fernandez, J. M. Gil-Narvion, Gordillo-Guerrero, D. Iñiguez, A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S.F. Schifano, B. Seoane, A. Tarancon, R. Tripiccione, D. Yllanes
aa r X i v : . [ c ond - m a t . d i s - nn ] M a r Dynamical transition in the D = M. Baity-Jesi,
1, 2, 3
R. A. Baños,
3, 4
A. Cruz,
4, 3
L.A. Fernandez,
J. M. Gil-Narvion, A. Gordillo-Guerrero,
D. Iñiguez,
A. Maiorano,
F. Mantovani, E. Marinari, V. Martin-Mayor,
J. Monforte-Garcia,
A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro,
M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo,
S.F. Schifano, B. Seoane,
2, 3
A. Tarancon,
R. Tripiccione, and D. Yllanes
2, 3 Departamento de Física Teórica I, Universidad Complutense, 28040 Madrid, Spain. Dipartimento di Fisica, La Sapienza Università di Roma, 00185 Roma, Italy. Instituto de Biocomputación y Física de Sistemas Complejos (BIFI), 50009 Zaragoza, Spain. Departamento de Física Teórica, Universidad de Zaragoza, 50009 Zaragoza, Spain. D. de Ingeniería Eléctrica, Electrónica y Automática, U. de Extremadura, 10071, Cáceres, Spain. Fundación ARAID, Diputación General de Aragón, Zaragoza, Spain. Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, and INFN, Ferrara, Italy. Dipartimento di Fisica, IPCF-CNR, UOS Roma Kerberos and INFN, La Sapienza Università di Roma, 00185 Roma, Italy. D. de Ingeniería, Electrónica y Comunicaciones and I3A, U. de Zaragoza, 50018 Zaragoza, Spain. Departamento de Física, Universidad de Extremadura, 06071 Badajoz, Spain. Dipartimento di Matematica e Informatica, Università di Ferrara and INFN, Ferrara, Italy
We study the o ff -equilibrium dynamics of the three-dimensional Ising spin glass in the presence of an externalmagnetic field. We have performed simulations both at fixed temperature and with an annealing protocol.Thanks to the Janus special-purpose computer, based on FPGAs, we have been able to reach times equivalentto 0 .
01 seconds in experiments. We have studied the system relaxation both for high and for low temperatures,clearly identifying a dynamical transition point. This dynamical temperature is strictly positive and dependson the external applied magnetic field. We discuss di ff erent possibilities for the underlying physics, whichinclude a thermodynamical spin-glass transition, a mode-coupling crossover or an interpretation reminiscent ofthe random first-order picture of structural glasses. PACS numbers: 75.50.Lk, 75.40.Mg
I. INTRODUCTION
The glass transition is a ubiquitous but still mysterious phe-nomenon in condensed matter physics [1–3]. Indeed, manymaterials display a dramatic increase in their relaxation timeswhen cooled down to their glass temperature T g . Exam-ples include fragile molecular glasses, polymers, colloids or,more relevant here, disordered magnetic alloys named spinglasses [4]. The dynamic slowing down is not accompaniedby dramatic changes in structural or thermodynamic proper-ties. In spite of this, quite general arguments suggest thatthe sluggish dynamics must be correlated with an increasinglength scale [5]. This putative length scale can be fairly dif-ficult to identify (a significant amount of research has beendevoted to this problem, see, e.g., [6–8]).In this context, spin glasses are unique for a number of rea-sons. To start with, they provide the only example of a mate-rial for which it is widely accepted that the sluggish dynamicsis due to a thermodynamic phase transition at a critical tem-perature T c = T g [9–11]. This phase transition is continu-ous, and the time-reversal symmetry of the system is sponta-neously broken at T c .Spin glasses are also remarkable because special tools areavailable for their investigation. On the experimental side,time-dependent magnetic fields are a very flexible tool toprobe their dynamic response, which can be very accuratelymeasured with a SQUID (for instance, see [12]). For verylow applied magnetic fields, one can measure glassy mag-netic domains of a diameter up to ξ ∼
100 lattice spac-ings [13, 14], much larger than any length scale identified for structural glasses [8]. On the theoretical side, spin glassesare simple to model, which greatly eases numerical simu-lation. In fact, special-purpose computers have been builtfor the simulation of spin glasses [15–19]. In particular, theJanus computer [17, 18] has allowed us to simulate the non-equilibrium dynamics from picoseconds to a tenth of a sec-ond [19, 20], which has resulted in detailed connections be-tween non-equilibrium dynamics at finite times and equilib-rium physics in finite systems [21, 22] (see also [23]).Now, Mean Field provides compelling motivation to inves-tigate spin glasses in an externally applied magnetic field. Al-though a magnetic field explicitly breaks time-reversal sym-metry, it has been shown that a phase transition should stilloccur when cooling the mean-field spin glass in an externalfield [24–26], which leads us to expect a sophisticated andstill largely unexplored behavior for short-range spin glasses.In this respect, we remark as well an intriguing suggestionby Moore, Drossel and Fullerton. These authors speculatethat the spin glass in a field sets the universality class for thestructural-glass transition [27, 28] (the statement is supposedto hold in d = is crucial).However, two major problems hamper further progress: (i)the mean-field theory, which is expected to be accurate abovethe upper critical dimension d > d u =
6, predicts a ratherdi ff erent behavior for spin and structural glasses; (ii) the be-havior of spin glasses in a field in d = p -spin glassmodel, with p -body interactions [29, 30] (for odd p , the time-reversal symmetry is explicitly broken). In the mean-field ap-proximation, the odd- p models display a dynamic phase tran-sition in their paramagnetic phase. Reaching thermal equi-librium becomes impossible in the temperature range T c < T < T g . The dynamic transition at T g is identical to the idealMode Coupling transition of supercooled liquids [31]. Thethermodynamic phase transition at T c is analogous to the idealKauzmann thermodynamic glass transition [3]. The thermo-dynamic transition is very peculiar: although it is of the sec-ond order (in the Ehrenfest sense), the spin-glass order param-eter jumps discontinuously at T c from zero to a non-vanishingvalue. On the other hand, for spin glasses in a field, mean-fieldtheory predicts a single transition (i.e., T c = T g ) and an orderparameter that behaves continuously as a function of temper-ature.However, mean-field theory is accurate only in high-enoughspatial dimensions, d > d u =
6, hence it is legitimate to won-der about our three-dimensional world. In fact, the ideal ModeCoupling transition is known to be only a crossover for su-percooled liquids. The power-law divergences predicted byMode Coupling theory hold when the equilibration time liesin the range 10 − s < τ < − s. Fitting to those power-laws, one obtains a Mode Coupling temperature T MC . How-ever, τ is finite at T MC (typically T MC is 10% larger than theglass temperature T g where τ ∼ seconds). A theory fora thermodynamic glass transition at T c < T g has been putforward [32–35], but it has still not been validated (however,see [36]).Regarding now our second problem, we recall that whetherspin glasses in a magnetic field undergo a phase transition hasbeen a long-debated and still open question. There are mainlytwo conflicting theories. The above mentioned mean-field pic-ture is the replica symmetry breaking (RSB) theory [25]. Onthe other hand, the droplet theory [37–40] predicts that thesystem’s behavior is akin to that of a disguised ferromagnet(i.e., no phase transition in a field). Recent exponents of thiscontroversy are [41–43]. We elaborate on concrete predictionsby both theories in Section II C. For now, we note that recentnumerical simulations in d = < d u =
6) [50].A di ff erent numerical approach has consisted in the study ofLevy graphs: one-dimensional systems where the interactiondecays with distance as a power-law J ( r ) ∼ / r ˆ a [51]. It hasbeen suggested that the critical behavior in d spatial dimen-sions can be matched with that of a Levy graph with an ap-propriate choice of the decay exponent ˆ a . The spin-glass tran-sition has been investigated in this way, both with and withoutan externally applied magnetic field. The latest result withinthis approach suggests that the de Almeida-Thouless line ispresent in four spatial dimensions, but not in three dimen-sions [52]. However, it has been pointed out that the match-ing between the decay exponent ˆ a and the spatial dimension d changes when a magnetic field is applied [53], a possibilitynot considered in [52]. Our scope here is to explore the dynamical behavior of d = L =
80, where we expect finite-sizee ff ects to be negligible [19]. Our time scales will range from1 picosecond (i.e., one Monte Carlo full lattice sweep [4]) to0 .
01 seconds. Hence, if the analogy with structural glassesholds, we should be able to identify the Mode Couplingcrossover. Our study will be eased by the rather deep theo-retical knowledge of the relevant correlation functions [54].Hence, we shall be able to correlate the equilibration time τ with the correlation length ξ .The layout of the remaining part of this work is as follows.In Section II we describe the model and our observables andwe give an overview of the di ff erent theoretical pictures thathave been put forward to explain the dynamics of the spinglass in a field. We pay particular attention to the specificpredictions for the observables that we are going to study(Sections II C and Section II D). In Section III we describethe di ff erent protocols that we have considered (simulating adirect quench and annealings with di ff erent temperature varia-tion rates). Next, in Section III B, we turn our attention to thecrucial and delicate issue of identifying intrinsic time scalesin the dynamics.Once this foundation has been laid, we delve into thephysical analysis and interpretation of our results. First, weconsider the high-temperature regime, where our simulationsreach equilibrium, in Section IV. We thus try to approach thecritical region from above. Afterwards, in Section V, we studythe low-temperature regime, where the relaxation times aremuch longer than our simulations (perhaps infinite). We con-sider two complementary approaches: in Section V A we trythe simplest ansatz for a low-temperature behavior compati-ble with a RSB spin-glass transition, while in Section V B weassume from the outset that no transition occurs. The spatialcorrelation of the system is studied in Section VI, where wecarry out a direct analysis of the correlation length. Finally, inSection VII we present our conclusions, evaluating the di ff er-ent physical scenarios on the light of our data.We also include two appendices: one describing some de-tails of our implementation and the other presenting a moredetailed look at the correlation functions. II. MODEL, OBSERVABLES AND THE DROPLET-RSBCONTROVERSY
In Section II A we describe the model that we have simu-lated. The observables that we consider are defined in Sec-tion II B. At that point, we will be ready to describe the di ff er-ent predictions of the droplet and RSB theories in Section II C.However, it has been recently suggested that spin-glass on afield behave just as supercooled liquids [27, 28]. The currenttheoretical predictions for supercooled liquids do not quitematch neither the droplet nor the RSB theories for spin glasseson a field. Hence, we briefly recall those predictions in Sec-tion II D. A. Model
We have studied a three-dimensional cubic lattice systemwith volume V = L ( L is the linear size) and periodic bound-ary conditions. On every node of the lattice there is an Isingspin σ x and nearest neighbors are joined by couplings J xy .The spins are dynamic variables, while the couplings J xy arefixed during the simulation (this is the so-called quenched dis-order [55]). The couplings are independent random variables: J xy = ± h x on every node. The magnetic field is Gaussianlydistributed with zero mean and variance H . The Hamiltonianof the model is H = − X h x , y i J xy σ x σ y − X x h x σ x , (1)where h x , y i means sum over nearest neighbors. A given re-alization of couplings, J xy , and an external field, h x , definesa sample . All of our results will be averaged over many sam-ples.It is important to realize that both the Hamiltonian (1) andthe probability distribution functions (pdf) for the couplingsand the magnetic field are invariant under a gauge symmetry.Let be ǫ x = ± σ x → ǫ x σ x , h x → ǫ x h x , J xy → ǫ x ǫ y J xy . (2)The observables defined in Section II B must be invariant un-der this symmetry. As explained in the literature (see, e.g.,[20, 21]), one forms gauge-invariant observables by consider-ing real replicas, namely copies of the system that evolve in-dependently under the same coupling constants and magneticfields. In particular, we consider four replicas per sample.The model (1) has been simulated on the Janus special-purpose computer [17, 18, 56, 57]. Janus’ limited memorydoes not allow us to study one independent, real-valued mag-netic field per site h x . We circumvented the problem usingGauss-Hermite integration [58]. Details of our implemen-tation, as well as comparison with PC simulations for real-valued h x , can be found in Appendix A.From the theoretical point of view, two main and contradic-tory theoretical predictions for the finite-dimensional modelin an external magnetic field can be found in the literature:the droplet theory (DT) and the Replica Symmetry Breakingtheory (RSB).In the droplet framework, the H = n , of the e ff ective field theory be zero, which im-plies the above cited degeneration between the anomalous andlongitudinal sectors of the theory.However, [61] started using the most general Hamiltoniancompatible with the symmetries of a replica symmetric phaseand relaxed the n = d . B. Observables
We are going to consider the temporal evolution of dif-ferent physical quantities in simulation protocols with tem-perature variation. To this end, we define t tot as the totaltime elapsed during the simulation (measured in Monte CarloSweeps, MCS, i.e., updates of the whole lattice) and the wait-ing time t w as the total time elapsed since the last change oftemperature. When this distinction is not important, we use t as a generic time variable.First, a couple of useful definitions of local quantities. Onevery node x of the lattice we have the local overlap: q x ( t ) = σ (1) x ( t ) σ (2) x ( t ) , (3)where the superscripts are the replica indices. The total over-lap is written as q ( t ) = V X x q x ( t ) , (4)where ( · · · ) means sample average (over the J ’s and h ’s), and V = L is the total number of spins in the lattice.In addition, we have focused in this work on the magneticenergy defined as E mag ( t ) = V X x h x σ x ( t ) . (5)and W ( t ) = − T E mag ( t ) / H . (6)We note that if the magnetic field is Gaussian distributed, thefollowing identity is fulfilled in thermal equilibrium : W = h q i . (7)To derive Eq. (7) integrate by parts the term h x e − h x / (2 H ) that appears on the disorder average for E mag , and recallthe equilibrium fluctuation-dissipation relation d h s x i / d h x = (1 − h s x i ) / T .Both q ( t ) and W ( t ) disregard the spatial dependence. Toaddress this issue, we should consider correlation functions.In a magnetic field, one may consider three di ff erent correla-tors G , G and G at equal time. The three of them can becomputed with four replicas: G ( r , t ) = V X x h σ x ( t ) σ x + r ( t ) i = V X x h σ (1) x ( t ) σ (1) x + r ( t ) σ (2) x ( t ) σ (2) x + r ( t ) i , (8) G ( r , t ) = V X x h σ x ( t ) σ x + r ( t ) ih σ x ( t ) ih σ x + r ( t ) i = V X x h σ (1) x ( t ) σ (1) x + r ( t ) σ (2) x ( t ) σ (3) x + r ( t ) i , (9) G ( r , t ) = V X x h σ x ( t ) i h σ x + r ( t ) i (10) = V X x h σ (1) x ( t ) σ (2) x + r ( t ) σ (3) x ( t ) σ (4) x + r ( t ) i . To gain statistics, we average over all the replica-index per-mutations that yield the same expectation values in the aboveequations.Now, because of the magnetic field, none of the correlators G , G and G tend to zero for large r , hence one needs to de-fine connected correlators. This problem was faced long ago,and the answer is in the three basic propagator of the repli-cated field theory namely the longitudinal ( G L ), anomalous( G A ) and replicon ( G R ) (see [24, 59]): G R = G − G + G , (11) G L = G + n − G +
12 ( n − n − G , (12) G A = G + ( n − G − ( n − G , (13)where n is the number of replicas of the e ff ective replica fieldtheory ( n should not be confused with the number of realreplicas that we hold fixed to four). Quenched disorder is re-covered from replica field theory in the limit n →
0. In thislimit, Eqs. (11, 12, 13) yield the replicon G R ( r , t ) = V X x ( h σ x ( t ) σ x + r ( t ) i − h σ x ( t ) ih σ x + r ( t ) i ) , (14)while G A and G L coalesce in the limit to a single propagator: G L = G A = G − G + G , (15)(strictly speaking, replica field theory applies only to equi-librium; the reader will forgive us for borrowing the natural equilibrium definitions for our dynamic computation at finite t ). Having in our hands two correlation functions (replicon andlongitudinal / anomalous), it is natural to compute the associ-ated susceptibilities χ ( t ) = Z d r C ( r , t ) , (16)where C stands for any of the correlators introduced in thissection.At finite t w , correlations surely are sizeable only withinsome characteristic correlation length ξ ( t w ). If one expectsa typical scaling form for such correlators C ( r ) ∼ f ( r /ξ ( t w )) r a , (17)where f is some long-distance cuto ff function, it is natural toextract the correlation length from integral estimators [19, 20] ξ k , k + ( t w ) ≡ I k + ( t w ) I k ( t w ) ∝ ξ ( t w ) . (18)where I k ( t w ) ≡ Z L / d r r k C ( r , t w ) . (19)The r in Eq. (19) is a shorthand for r = ( r , ,
0) and permu-tations. The signal-to-noise ratio is increased if one imposesa long-distance cuto ff in Eq. (19). In particular, we stop in-tegrating when the relative error in C ( r , t w ) grows larger than1 / ff ect of the tail with an exponentialextrapolation, although this is a very minor e ff ect, see [20] andAppendix B 1).Notice that a similar method can be used to obtain estima-tors for the susceptibility. Assuming isotropy (see [20]) wecan write Eq. (16) as χ ( t ) = π Z L / d r r C ( r , t ) = π I ( t ) , (20)which can be computed with the long-distance cuto ff .Finally, let us mention that quantities such as C ( r , t w ), due tothe slow temporal evolution and the strong sample-to-samplefluctuations, turn out to be very rugged functions of t w . This isparticularly bad for the computation of ξ , since it can havean unpredictable e ff ect on the cuto ff . In order to avoid thisproblem, we perform a temporal binning, averaging C ( r , t w )in blocks of 4 consecutive times and considering it as a func-tion of the geometrical mean of the times in each block. Thissmoothing procedure yields a significant error reduction. C. The droplet versus RSB controversy in terms of ourobservables
We summarize here the major di ff erences among the pre-dictions from the droplet and RSB theories with an emphasison the quantities defined in Section II B. We discuss as wellthe di ff erent predictions for the characteristic time scale forequilibration, τ (see Section III B). As we shall see, the pre-dictions of the droplet theory are more detailed, hence let usstart there.
1. Predictions from the droplet theory
According to the droplet theory, there is no phase transitionin a field. This means that the long-time and large-system lim-its, t w → ∞ and L → ∞ , can be taken in any order. In thiswork, we will always take first the limit of large L (because L ≫ ξ ( t w ) in our simulations, see [19]). Thus, it is a conse-quence of the droplet picture that Eq. (7) should always holdfor large lattices and then large times.The physics at low temperatures is predicted to be ruledby a fixed point at T = H = y T = − θ , (21) y H = d − θ , (22)where θ is the droplet sti ff ness exponent and d is the spacedimension ( d = t w → ∞ . A scaling law follows fromEqs. (21,22), which should hold at least for small temperaturesand magnetic fields: ξ ( T , H ) = H / ( d − θ ) F ( T H θ/ ( d − θ ) ) . (23) F ( x ) is a scaling function, which is supposed to remainbounded when x → ξ Ψ [ Ψ is asecond droplet exponent which should satisfy θ ≤ Ψ ≤ ( d −
1) [40]]. Hence, since ξ is expected to remain bounded atall temperatures, the dynamics is predicted to be of Arrheniustype τ ∝ exp (cid:20) ξ Ψ ( T , H ) T (cid:21) . (24)Of course, if ξ grows significantly with T for some tem-perature interval, Eq. (24) would predict an apparent super-Arrhenius behavior.Now, in order to make some use of Eqs. (23,24), it ismandatory to have an estimate for the exponents θ and Ψ . Forthe barrier exponent Ψ , we may quote an experimental deter-mination on the Ising system Fe . Mn . TiO , Ψ ∼ .
03 [14].As for the θ exponent, we are aware of two types of com-putations. On the one hand, the properties of ground stateshave been analyzed. The comparison of periodic and antiperi-odic boundary conditions yields θ ∼ . θ = . θ = θ ina rather more direct way, by considering the behavior of thespatial correlation functions. This approach yields θ = . θ ,we shall use both θ = .
24 and θ = .
61 when trying to assessEq. (23).
2. Predictions from the RSB theory
The RSB theory is based on the solution of a mean-fieldmodel (the Sherrington-Kirkpatrick model). Extending thetheory below its upper critical dimension d u = d = T c ( H ), the locus of the deAlmeida-Thouless line, in the form of a power law: ξ ( T , H ) ∝ | T − T c ( H ) | ν , (25)Consequently, the relaxation time should diverge at T c ( H ). Itmay do so in the form of critical slowing down: τ ∝ ξ z , (26)which defines the dynamic critical exponent, z . This is consis-tent with the mean-field analysis of the dynamics [69]. How-ever, it is possible that the relevant fixed point be at T =
0, yet,at variance with the droplet theory, at a finite magnetic field H c (see Figure 1 in [42]). This is precisely the situation encoun-tered in the random-field Ising model (see, e.g., [70]), whichsuggests that an activated dynamics, as in Eq. (24), might ap-ply instead of (26).At any rate, a defining feature of the RSB picture is the non-triviality for the probability distribution of the spin overlap q ,Eq. (4). If, for T < T c ( H ), one takes first the limit of largetimes t w → ∞ and only afterwards lets L → ∞ , all valuesof the overlap in an interval q min ≤ q ≤ q max can be found.Now, we shall be studying the problem with the reverse orderof limits (i.e., first L → ∞ , and only afterwards t w → ∞ ).Since our simulations will start from a disordered state, with q =
0, for temperatures below the dAT line one expectslim t w →∞ q ( t w ) = q min . (27)From Eqs. (27) and (7) one can conclude:lim t w →∞ (cid:0) W ( t w ) − q ( t w ) (cid:1) = h q i − q min . (28)In the droplet model, the rhs of the previous equation is justzero, while it is non-zero in a RSB spin-glass phase. D. A new scenario: supercooled liquids
Probably, the most successful theory for supercooled liquiddynamics is the Mode-Coupling theory (MCT) [31, 71]. Thisis a very rich theory, with many predictions. We shall contentourselves by recalling the results most directly relevant to ourdiscussion.It is now well understood that MCT is a Landau-like or clas-sical theory [72–74]. As such, it is exact above eight spatialdimensions. In our d = T d τ ∝ T − T d ) γ . (29)The correlation length diverges as well at T d as ξ ∝ T − T d ) / . (30)So τ ∝ ξ γ as in Eq. (26). The critical exponent γ is given interms of two other critical exponents a and b , whose precisedefinition is not relevant to us now (for instance, see [73]): γ = a + b . (31)The two exponents a and b are related to a crucial exponent λ through the equation ( Γ ( x ) is Euler’s Γ function [58]): λ = Γ (1 + b ) Γ (1 + b ) = Γ (1 − a ) Γ (1 − a ) . (32)Although λ is often treated as an adjustable parameter, we re-mark that it is actually a static renormalized coupling con-stant [75]. As such, it may be systematically computed (forinstance, in an hypernetted-chain approximation [73, 74] or ina numerical simulation). A typical value for supercooled liq-uids is λ ∼ . a and b only if λ ≤
1. Hence, if one finds λ > β -relaxation inliquids. Roughly speaking this regime covers the temperaturerange corresponding to 10 < τ < , as measured in MonteCarlo steps (see, e.g., [76] whose numerical findings are con-sistent with λ = . T d neither τ nor ξ diverge. In-stead, activated processes enter the stage, playing the role of anon-perturbative phenomenon that erases the mode-couplingtransition [77]. The dynamics is expected to become of acti-vated type, as in Eq. (24). The behavior of ξ upon lowering the temperature is still unclear. Some expect that ξ will diverge atthe Kautzman temperature [32–35], but the issue is still underactive investigation [36]. III. DYNAMIC PROTOCOLS AND THE IDENTIFICATIONOF THE RELEVANT TIME SCALEA. The dynamic protocols
We have performed several independent sets of simula-tions, both at a fixed temperature ( direct quench ) and withtemperature changes ( annealing ) . Our aim was to identifytemperature-dependent, intrinsic properties (by intrinsic wemean independent of the dynamic protocol that we followed).In all cases we have simulated four replicas for each sam-ple. We have consider external fields H = .
1, 0 . . L =
80. We storethe full configuration of the system for all times of the form t w = [2 i / ], where [ · · · ] denotes the integer part. From thesestored configuration we can compute any physical quantity of-fline.In the simulations at a fixed temperature we ran 462 sam-ples for each external field at T = . H = . MCS.
L T H
MCS N
80 0 . . . . . . . . . . . . N is the number of samplessimulated. L [ T init , T end ] H t base ∆ T MCS N
80 [2 . , .
4] 0 . –10 . × . , .
4] 0 . –10 . × . , .
85] 0 . × . × . , .
4] 0 . –10 . × T ini and T end mark the initialand final temperatures of the annealing procedure. We decrease thetemperature in ∆ T increments and we run for a time of t base × (cid:18) T init − T ∆ T (cid:19) MCS at each temperature T . For H = . cold and the hot annealings. For H = . , . The second set of simulations was performed with an an-nealing algorithm . We started the simulation at a high tem- W , q T q ( t tot ) W ( t tot ) 0.1 0.3 0.5 0.710 W , q T t tot FIG. 1. (Color online) Overview of our two annealing simulationsfor H = .
2. We show the evolution of q ( t tot ) [the red (lower) curve]and W ( t tot ) [the blue (uppe) curve] during the whole simulation. Thecontinuous black line indicates the temperature throughout the sim-ulation (see right-hand vertical axis for the scale). The top panelcorresponds to the cold annealing and the bottom panel to the hot annealing (see Table II). Notice how, in the former, the system is notable to reach equilibrium for the lowest temperatures (as signaled bydi ff erent values of W and q at the end of each temperature step). perature T init = T . After t base MCS, we change the tempera-ture to a new one ∆ T cooler, i.e., T = T − ∆ T , and we take2 t base steps. We iterate this procedure, decreasing the temper-ature by a fixed step ∆ T and increasing the number of stepsat fixed temperature in a geometric progression until we reachour lowest temperature T end . That is, for a given temperature T k the total elapsed time is in the range t k ≤ t tot ≤ t k + , where t k = (2 k − t base . (33)We performed in every case the annealing from T = T init = . T end = . t base = i , i = , . . . , ∆ T = . . × × t base MCS in each sample and replica.See Table II for more details. Throughout the paper we shallrefer to these runs as the cold annealings . If we do not say oth-erwise, a mention in the paper to a cold annealing will alwaysrefer to the slowest one, with t base = .Finally we have run a yet slower annealing constrained tothe high-temperature region in order to obtain equilibrium re-sults. In this case, only for H = .
2, we go from T init = . T end = .
85, with t base = × and ∆ T = . hot annealing .We show in Figure 1 an overview of our two sets of anneal-ing runs for H = .
2. We represent both W ( t tot ) and q ( t tot ),which, for large t w , should converge to the same value (in the W , q H = 0.1, T =1.3 q ( t w ) W ( t w ) 0.47 0.51 0.5510 W , q t w H = 0.1, T =0.8 FIG. 2. (Color online) Detail of our annealing simulation for H = . ff erent temperatures, T = . T = .
8. Notice how inthe second case, despite the longer simulation time, the system doesnot reach equilibrium. paramagnetic phase). As we can see, in the cold annealing thiscondition is not satisfied for several of the lower temperatures,signaling that the system has fallen out of equilibrium. Thehot annealing reaches the equilibrium regime for the wholetemperature range.A more detailed picture of these two regimes can be seenin Figure 2, which represents W ( t w ) and q ( t w ) for two tem-peratures and H = .
1. As we can see, for the higher tem-perature both observables converge to the same long- t w limit.For the lower temperature, however, the two quantities are farapart during the whole simulation and seem to have a di ff er-ent asymptote. This indicates that either the equilibration timeis much larger than our simulation or we are in a spin-glassphase. Deciding between these two possibilities is the maingoal of this paper. B. Identification of intrinsic time scales
Our strategy will rely on the study of the di ff erence between W ( t w ) and q ( t w ) as a function of time, which should decay tozero in the high-temperature phase. Some inspiration comesfrom the classic paper by Ogielski [16], although we shall useother approaches as well.We show in Figure 3 the behavior of W ( t w ) − q ( t w ) for T = . H = . W − q di ff ers greatly from one protocol to the other. It is very large W – q t w QuenchCold annealingHot annealing
FIG. 3. (Color online) Our three di ff erent simulations at T = . H = .
2. We show the di ff erence W − q as a function of thewaiting time t w , whose decay to zero determines the time scale ofthe simulations. Despite having very di ff erent initial conditions, thethree simulation protocols need roughly the same time to reach equi-librium. We shall quantify this assertion in this Section III B. Noticehow, for the two annealing protocols, the behavior of W − q is linearin log( t w ) for a long time range. (well out of the graph’s scale) for the direct quench, whichstarts with a random configuration (representing a very hightemperature) and it is smallest for the hot annealing, whichhas a small temperature step.Nevertheless, all protocols seem to need roughly the samenumber of MCS to reach equilibrium, as evinced by the merg-ing of the curves at t w ∼ . This observation gives us somehope of determining an intrinsic time scale τ , depending onlyon the system’s temperature and not on its history.In principle, the robust way to compute τ would be to per-form a calculation analogous to Eq. (18), replacing C by W − q and r by t w . Unfortunately, in the interesting temperaturerange the maximum of the integrand t k w [ W − q ]( t w ) is alwaysin a region with a dismal signal-to-noise ratio (cf. Figure 1 in[20]). Therefore, we have to resort to more phenomenologicaldeterminations.A traditional way to identify this time scale, which wasfound adequate in the absence of a field [16], is fitting thedi ff erence W ( t w ) − q ( t w ) to a stretched exponential decay: W ( t w ) − q ( t w ) = A t x w exp " − (cid:18) t w τ ′ (cid:19) β . (34)From this fit one gets a characteristic time τ ′ , which we canuse as our time scale.Computing a fit to Eq. (34) is di ffi cult due to the extremecorrelation of our data, which prevents us from inverting itsfull covariance matrix (necessary to define the χ goodness-of-fit indicator). Therefore, we consider only the diagonal partof the matrix in order to minimize χ and take correlations intoaccount by repeating this procedure for each jackknife blockin order to estimate the errors in the parameters. This is, ofcourse, only an empirical procedure, but one that has beenshown to work well under these circumstances (see, e.g., [20], β TH = 0.1 H = 0.2 H = 0.3 FIG. 4. (Color online) Behavior of the stretching exponent β (seeEq.(34)) as a function of T for our three simulated magnetic fields. especially sections 2.4 and 3.2).The results of these fits are gathered in Table III. We do notreport the value of the (diagonal) χ / d.o.f. because it is, inall cases, χ / d.o.f. < τ ′ longer than the simulation time).A possible source of uncertainty in our determination of τ ′ is the dependence of the fit on the value of β . Indeed, for each T we are fitting simultaneously for x , A , τ ′ and β in (34).However, a small variation in β can have a large e ff ect on τ ′ ,which may lead us to think that the fit is unstable and unre-liable. Fortunately (see Fig. 4), β is actually a very smoothmonotonic function of T , which leads us to believe that thefitting procedure is sound.There is a final di ffi culty with this functional form: τ ′ onlyhas a straight interpretation as a correlation time (i.e., as anestimator for τ ) if β ≈
1. However, in the interesting tempera-ture range, 0 . . β . .
3. This means that the actual value of τ ′ cannot be interpreted directly as an estimator for τ , but stillwe expect its divergence at the dynamical transition point tobe intrinsic, as discussed in Section IV.Notice, finally, that the values of τ ′ and β computed at thesame temperature for the hot and cold annealing protocols for H = . τ , completely phe-nomenological, comes to mind from a visual inspection ofFigure 3. Indeed, we can see that for both annealing proto-cols, the di ff erence W − q is linear in log( t w ) for a very widetemperature range. This behavior is more clearly shown inFigure 5, which represents this quantity for two di ff erent tem-peratures.Moreover, if we represent this linear behavior as W ( t w ) − q ( t w ) ≃ B " − log( t w )log( τ ′′ ) , (35) H T
Annealing Stretched exponential, Eq. (34) Linear fit, Eq. (35) τ ′ / β log( τ ′′ ) τ ′′ / t w ) (35) (log being the natural logarithm). For each value of H we report the last few temperatures before the system fallsout of equilibrium. All the reported fits have (diagonal) χ / d.o.f. < W – q T = 1.1, Cold annealing T = 1.1, Hot annealing 0 0.01 0.02 0.03 0.04 5 10 15 20 W – q log( t w ) T = 0.9, Cold annealing T = 0.9, Hot annealing FIG. 5. (Color online) Detail of our two annealing simulations for H = . ff erent temperatures, T = . T = .
9. Werepresent the di ff erence W − q as a function of log( t w ) (we have usedlog as the natural logarithm throughout the paper). As pointed out inFigure. 3, in this representation W − q is linear for a very long range.More interestingly, if we compute linear fits to (35) (continuous blacklines), the intercepts with the horizontal axis are independent of thesimulation protocol. we can see from Figure 5 that the value of τ ′′ does not de-pend on the annealing rate and is therefore dependent only onthe temperature. Furthermore, since there is no stretching ex-ponent, we can take τ ′′ directly as an estimate of the actual intrinsic time scale of the system: τ ′′ ≈ τ .Of course, Eq. (35) is only empirical and cannot be correctfor very long times (it would predict an unphysical negativevalue of W − q for t w > τ ′′ ), but its simplicity and robustnesscompensate for this problem. We give the values of τ ′′ for sev-eral temperatures in Table III (in accordance with the previousdiscussion on the meaning of β , notice that τ ′′ is more than anorder of magnitude larger than τ ′ ). In the following, we shalluse both τ ′ and τ ′′ to study the possible critical behavior ofthe system. IV. THE EQUILIBRIUM REGIME
In this section we consider the temperature dependence ofthe characteristic time scales τ identified in Sect. III B. Fol-lowing Ogielski [16] as well as experimental studies (e.g.,[78]), we shall fit our equilibrium data to power-law diver-gence: τ = τ (cid:2) T − T c ( H ) (cid:3) ν z , (36)where τ is a microscopical time. We use the traditional nota-tion, where T c ( H ) is a critical temperature, ν , is the correlationlength critical exponent and z is the dynamical critical expo-nent. However, as we shall see in Section VI, this divergencemay or may not correspond to an actual phase transition at T c ( H ).In Figure 6 we show the relaxation time τ ′ computed withthe stretched exponential (34) in Table III as a function oftemperature for the three simulated magnetic fields. We alsoshow fits to the power-law divergence at finite T highc ( H ) (thesuperscript “high” refers to the fact that we are using high-temperature data, cf. Section V). We have obtained the fol-lowing values:0 τ ' T H = 0.1 H = 0.2 H = 0.3 FIG. 6. (Color online) Behavior of the characteristic time τ ′ of thestretched exponential (34) as a function of the temperature for thethree magnetic fields simulated. We also plot our fits to a power-lawdivergence of τ ′ for a finite temperature, Eq. (36). • H = . T highc = . z ν = . . . ≤ T ≤ . χ / d . o . f . = . / • H = . T highc = . z ν = . . . ≤ T ≤ . χ / d . o . f . = . / • H = . T highc = . z ν = . . ≤ T ≤ . χ / d . o . f . = . / τ ′ reliably. The choice of the fitting rangeis not critical, several temperatures can be added or elimi-nated without altering the results significantly (especially inthe high-temperature end). For H = .
2, we use the τ ′ for thecold annealing, since they have smaller error bars than thosefor the hot annealing. We have checked that the τ ′ at T = . τ ′ doesnot have a straightforward interpretation as a relaxation time,since β ,
1. Therefore, the above values of T highc might be anartifact of our way of estimating τ . In order to dispel this pos-sibility, we have recomputed the fits to (36), this time using therelaxation time τ ′′ computed with the linear fit in log t w (35).Now the fit parameters are: • H = . T highc = . z ν = . . ≤ T ≤ . χ / d . o . f . = . / • H = . T highc = . z ν = . . ≤ T ≤ . χ / d . o . f . = . / • H = . T highc = . z ν = . . ≤ T ≤ . χ / d . o . f . = . / χ for all fits. The values of T c ( H ) are consistent for both sets of fits, while the values of z ν are a little higherfor the fit with τ ′′ (2 . H = .
3, 1 . H = . . H = . W − q . A. The relaxation time in the supercooled liquids approach
As discussed in Section II D, the MCT predicts a divergenceof the autocorrelation time at a temperature T d , as in Eq. (29).In principle, Eq. (29) is exactly the same as Eq. (36), whichwe have just used to characterize the growth of the relaxationtimes. The crucial di ff erence is that we have used our lowestthermalized temperatures and have assumed that the growthof τ was related to an actual critical divergence (as evincedby our notation of z ν for the exponent). On the other hand, inthe supercooled liquids literature, Eq. (29) is used in a highertemperature range corresponding to 10 < τ < (notice, forinstance, that our values of z ν are very high compared to thevalues of γ that can be found in the MCT literature). For lowertemperatures, the behavior of τ deviates from (29), because ofthe emergence of activated processes.Therefore, if we wanted to follow a supercooled liquids ap-proach, we should first fit τ ′ to (29) in the high temperaturerange and then move on to an exponential growth: τ ′ = C ( T − T d ) γ , τ ′ . , (37) τ ′ = exp (cid:2) D / ( T − T c ) γ ′ (cid:3) , τ ′ & . (38)Unfortunately, our simulations are not suited to the determi-nation of small τ , so the fit to (37) will probably be plaguedby strong systematic e ff ects.We consider only H = . H = . H = . τ ′ are ratherunstable for T > . τ ′ to (37) in the range1 . ≤ T ≤ . < τ ′ < ). Theresult is T d = . γ = . χ / d.o.f. = . / τ ′ in the range T ≤ . z ν ≈
8) and attempt a fit to (38) instead. Unfortunately,fitting for T c and γ ′ simultaneously is simply not possible withour data (the resulting error in T c is greater than 100%). Inshort, we can only say that a temperature dependence of τ according to (37) and (38) cannot be excluded, but we cannotmake this statement more quantitative. Nevertheless, we shallreturn to the possibility of activated dynamics in Sections V Band VI. V. NON-EQUILIBRIUM REGIME
As already explained in Section III B, our annealing rateeventually becomes too fast, compared to the growth of τ ( T )1upon cooling. At that point, the simulation falls out of equi-librium. We enter here the reign of extrapolation, which isalways rather risky.We shall extrapolate our data to long times following twovery di ff erent strategies. In Section V A we extrapolate usingpower laws. The outcome will be consistent with the RSBtheory. On the other hand, in Section V B we use the linear-log extrapolation (see Sect. III B), which assumes from theoutset that no phase transition occurs. A. Power-law extrapolations to long times
So far, we have been working in the high-temperaturephase, where W − q goes to zero for long times. We saw that,as we lower the temperature, the associated relaxation timegrows very quickly and eventually becomes much larger thanour simulations. This rapid growth was actually consistentwith a power-law divergence of τ at finite T .In this section we shall take a complementary approach. Wenow work in the low-temperature regime, where τ is eitherinfinite or, at the very least, much larger than our simulationtimes. In this regime, rather than assuming that W − q goes tozero for long times, we can try to extrapolate for a (possibly)non-zero asymptote.Following the literature [79, 80], we shall first attempt astudy in the total annealing time t tot , considering di ff erent an-nealing rates (Section V A 1). Then we shall repeat the analy-sis using only t w (as in the rest of the paper) in Section V A 2.
1. Study in t tot for di ff erent annealing rates In the limit of a very slow annealing, the simplest ansatz forthe low-temperature behavior of W − q is a power-law decay(cf. [16]): W ( t tot ) − q ( t tot ) = a ( T , H ) + bt x tot . (39)Eventually b and x could also depend on the temperature [16]and on the external magnetic field. Notice that, in contrastwith the rest of the paper, here we are considering the totaltime t tot since the simulation started [79, 80], not just the time t w since the last temperature change.Should the system experience a spin-glass transition, wewould expect a ( T , H ) > T [recall Eq. (28)].This asymptote would decrease as we increase the tem-perature until eventually, at some temperature T lowc ( H ), a (cid:0) T lowc ( H ) , H (cid:1) =
0. This description is consistent with a qual-itative look at W − q (recall, for instance, Figure 2).If the RSB picture is correct, we would expect T lowc ( H ) tocoincide with the divergence of τ and signal a thermodynamicphase transition, that is, T lowc ( H ) = T highc ( H ) = T c .Equation (39), with a >
0, should hold only deep in thespin-glass phase. If we approach the transition from below (intemperature) we would start to see the critical e ff ects of the(thermodynamical) critical point, and the exponent x would x H = 0.3 H = 0.2 H = 0.1 0 0.05 0.1 0.15 0.4 0.6 0.8 1 a T FIG. 7. (Color online)
Top:
Exponent x of the extrapolation of thedi ff erence between W ( t tot ) and q ( t tot ), Eq. (39), as a function of tem-perature, for our three external magnetic fields. Value computed fromthree-parameter fits. Bottom:
As above, for the asymptote a . begin to be controlled by this critical point and not by the“critical” spin-glass phase (Goldstone phase). So, in the criti-cal region we should expect: W ( t tot ) − q ( t tot ) = ft x c tot , (40)where in general x c (driven by the critical point) [81] shouldbe di ff erent from x (driven by the spin-glass phase).From the previous discussion, and assuming the onset of aphase transition, it is clear that the x exponent should take aconstant value at lower temperatures (here we are assumingthat the phase transition is universal in the magnetic field),then change as we reach the critical region.Now, in order to study the decay of W − q we obviouslyneed to follow the time evolution along several orders of mag-nitude. However, as described in Section III A, for a fixedtemperature T k = T init − k ∆ T the total time elapsed in ourannealing simulations varies in the range t k ( t base ) ≤ t tot ≤ t k + ( t base ) ∼ t k ( t base ), which is too narrow in a logarithmicscale. Therefore, instead of analyzing the data for all the t tot in a given annealing simulation, we shall combine all our coldannealing simulations for di ff erent values of t base . That is, foreach temperature T k we take the value of W − q at t k + ( t base ),for t base = i , i = , . . . ,
5. Thus, we get for each temper-ature T k a series of six values of [ W − q ]( t tot ) in the range t k + ( t base = ≤ t tot ≤ t k + ( t base = ), which we fit to (39)[recall that t k + ( t base ) = (2 k + − t base ].The resulting values of a ( T ) and x ( T ) are plotted in Fig-ure 7. As we can see, the qualitative picture is very much2what we painted above. In particular, we obtain a positivevalue of the asymptote a for low temperatures, which goesto zero at a temperature T lowc not very di ff erent from T highc —we can estimate T lowc ( H = . ≈ . , T lowc ( H = . ≈ . , T lowc ( H = . ≈ . x is roughly constant (and independent of H ) at lowtemperatures, while it grows noticeably as we approach T lowc .However, we must caution the reader that the fits we havejust discussed are rather delicate. In particular, even afterdiscarding the two smallest t tot (so we are left with a three-parameter fit to four points), we find that values of the χ goodness-of-fit estimator are sometimes very high. In par-ticular, the fits for H = . , . χ / d.o.f. ≤ . / H = . χ / d.o.f. that can be in excess of 8 /
1, clearlyunacceptable. More worryingly, if we shift the fitting windowwe find that the fitted values for x and a decrease noticeablywith increasing t tot (the change in x can be as high as 50% justby shifting the fitting window so that we discard the longesttime but include an extra point in the lower end of the range).Still, for H = . , . t tot gives reasonable fits (for H = . t w .Of course, from the above arguments one could think that,for long enough t ttot , the exponent could decrease so much thatthe asymptote a would become zero. In order to check againstthat possibility, and to obtain a sort of lower bound for a , wehave also attempted fits to the following function W ( t tot ) − q ( t tot ) = a ′ + b ′ log( t tot ) x ′ . (41)Using this fitting function we get good fits with values of x ′ ≈ a ′ for T < T highc , so the qualitative picture isthe same. The determination of T lowc is a little lower, but stillcompatible with our T highc .The next step in this study is the investigation of the scalingin W and q separately. To this end, we are going to considerfits of the form W ( t tot ) = h q i + b ′′ t x tot , q ( t tot ) = q min + b ′′′ t x tot , (42)where for both observables we take the same value of x thatwe computed in Figure 7. The resulting plots of h q i ( T )and q min ( T ) can be seen in Figure 8 for H = .
2. Ex-cept for H = .
2, we obtain excellent fits both for q andfor W ( χ / d.o.f. < T ≥ . H = .
2) we do not need to extrapolate, since we reach equi-librium in our annealing simulations. There is a small gapwith no extrapolations, corresponding to temperatures above T c , where (39) does not work but we cannot reach equilib-rium. The qualitative picture is what one would expect in theRSB scenario. T q min h q i FIG. 8. (Color online) Extrapolations to infinite time of W ( t tot ) → h q i and q ( t tot ) → q min for H = .
2, according to (42). We use the sameexponent x computed in Figure 7. q m i n H B T = 0.4 T = 0.5 T = 0.6 FIG. 9. (Color online) Scaling of the minimum overlap with themagnetic field. A simple power-law behavior, q min ∼ H B , workswell for all our subcritical temperatures, even though we are far fromthe small- q min limit. Finally, we can consider the scaling of q min with the mag-netic field at fixed temperature. As the magnetic field goes tozero, so must q min and we could expect a behavior of the kind q min ( H ) ∼ H B . Indeed, a rough dimensional analysis [80] tellsus that B = θ (0) / [ D − θ (0) / D = θ (0) = . B = . q ( H ) = C H B (43)for T = . , . . T c for our three magnetic fields). The results are B ( T = . = . B ( T = . = . B ( T = . = . B = . q min that we are seeing). In all cases we obtain excellent3 x H = 0.3 H = 0.2 H = 0.1 0 0.1 0.2 0.4 0.6 0.8 a T FIG. 10. (Color online)
Top:
Exponent x of the extrapolation ofthe di ff erence between W ( t w ) and q ( t w ), Eq. (44), as a function oftemperature, for our three external magnetic fields. Values computedin a three-parameter fit. Bottom:
As above, for the asymptote a .Lines are linear fits to the points where a ( T , H ) > H = . T = . values of the χ estimator. We have plotted q B min in Figure 9,using for all fields an intermediate value of B = .
2. Study in t w As discussed above, the total time t tot since the simula-tion started is the more physical variable to conduct the low-temperature study. However, we have seen in Sections III Aand IV that we can also study the relaxation of the system in t w in a consistent way. Therefore, it is interesting to repeat thestudy of Section V A 1 taking only our slowest cold anneal-ing (with t base = ) and studying for each temperature therelaxation in t w .To this end we consider an equation analogous to (39): W ( t w ) − q ( t w ) = a + bt x w . (44)We have computed fits to (44) for all our lower temperatures,finding that the power-law decay describes the behavior of W − q with great precision. Indeed we find for the standardfigure of merit χ / d . o . f . < . χ only with the diagonal part of the covariance matrix). Weshow in Figure 10 the fit parameters of x and a as a functionof temperature (upper and lower panels respectively). Clearly,this study leads to much lower values of x than the analysis W – q log( t w ) T = 0.8 T = 0.7 T = 0.6 FIG. 11. (Color online) W − q as a function of log( t w ) for several lowtemperatures and H = .
2. Even in this temperature range we canidentify a wide linear regime. in t tot (but recall that in the previous section the value of x decreased if we shifted the fitting window to longer times, aproblem we do not have here). Indeed, the values of a aremore similar to those computed with (41). Finally, let us re-call that relaxation exponents of O (10 − ) have already beenseen in the H = x are much lower.This is because, since we are using the run with the longest t base , the e ff ective time at t w = W − q is already very low at t w =
0, since it has already evolved fora considerable time at higher temperatures).Again, we can use the temperature at which a becomes zeroas our estimate of T lowc ( H ). We note in Figure 10–bottom thatthe statistical errors for a are small only when it is positive.Hence we have located the zeros by performing first a linearfit to these points, and then finding the root of the linear func-tion. In order to take care of the extreme data correlation,we use a jackknife procedure: fit with the diagonal part ofthe covariance matrix, but then perform separate fits for eachjack-knife block [20]. We obtain T lowc ( H = . = . T lowc ( H = . = . T lowc ( H = . = . q ( t w ) and W ( t w ) to infinitetime separately. Again, we obtain χ / d.o.f. < t w & q min at low T ).In short, we can say that assuming a power-law decay of W − q at low temperatures leads to a picture consistent with aRSB spin-glass transition at T lowc ( H ) ≈ T highc ( H ). B. Assuming no phase transition: linear-log extrapolations tolarge times
In the previous section we assumed that the decay of W − q followed a power law at low temperatures and tried to deter-4 H T log( τ ′′ ) T log( τ ′′ )0.3 0.7 27.9(4) 19.55(26)0.6 43.0(7) 25.8(4)0.5 80.2(1.3) 40.1(6)0.4 179(4) 71.5(1.5)0.2 0.8 28.4(6) 22.7(4)0.7 42.6(9) 29.8(6)0.6 78.8(1.8) 47.3(1.1)0.5 178(6) 89.1(2.8)0.4 409(13) 163(6)0.1 1.0 26.7(6) 26.7(6)0.9 39.9(9) 35.9(8)0.8 67.6(1.6) 54.1(1.3)0.7 115(3) 80.6(2.3)0.6 226(6) 135(4)0.5 450(14) 225(8)0.4 — —TABLE IV. Computation of a lower bound for the relaxation timeat low temperature using the linear fit in log( t w ) of Eq. (35). Foreach H we include the values of log( τ ′′ ) for temperatures in the non-equilibrium regime (i.e., lower than those included in Table III). For H = . T = . W − q is flat within errors and wecannot determine any log( τ ′′ ). Data from the cold annealings. mine the point where the asymptote became positive. In thissection we take the opposite approach and will assume thatthere is no phase transition, that is, that W − q goes to zero forall T > W − q was linear in log( t w ). In Sec-tion III B we used this functional form to estimate a character-istic time scale τ ′′ . Naturally, the real curve W − q must deviatefrom (35), otherwise it would become negative for t > τ ′′ , butat high temperature we found that the curvature was notice-able only at the very end of the simulation, where W − q wasalready very small (even compatible with zero).In this section, we are going to use (35) in order to obtaina lower bound for the relaxation time of the system. Indeed,if we look at Figure 11, we can see that even for very low T the di ff erence W − q behaves linearly in log( t w ) for a longtime scale, before slowing down its decay. Therefore, in thevery reasonable assumption that there is no convexity change,we can use the parameter log( τ ′′ ) of the linear fit in order toobtain a lower bound for the actual relaxation time τ of thesystem.With this procedure, we obtain a finite lower bound for τ even for very low temperatures (see Table IV). The resultingvalues of log( τ ′′ ) are enormous (as an amusing comparison,the age T of the universe measured in MCS is log( T ) ∼ τ ′′ ). Therefore, evena rather loose lower bound gives us a wildly growing timescale.In order to make this statement more quantitative, let us
10 100 1000 0.3 0.5 0.7 1 1.3 T l og ( τ '' ) T H = 0.1 H = 0.2 H = 0.3 FIG. 12. (Color online) Plot of a lower bound for T log( τ ), computedusing (35), against T for our three magnetic fields. At very low tem-peratures this quantity is compatible with a Vogel-Fulcher-Tammanndivergence at T =
0, Eq. (47). Notice the change of regime at around T highc ( H ) (recall Figure 6). recall that, in the droplet picture, ξ ( H , T ) is finite even at T = H >
0. Therefore, for low temperatures we would expectlog τ = ξ (0 , H ) Ψ / T , (45)or, in other words, we would expect T log( τ ) to be constant.However, see Table IV, we find that even our lower boundlog( τ ′′ ) grows much faster than predicted by the droplet the-ory.We can take this one step further. In Figure 12 we haveplotted T log( τ ′′ ) against T in a log-log scale. As we can see,the data for low temperature are well described by a Vogel-Fulcher-Tammann divergence at T = T log( τ ′′ ) = AT c . (46)A fit to (46) for our three magnetic fields gives: • H = . c ( H = . = . T ≤ .
7, with χ / d.o.f. = . / • H = . c ( H = . = . T ≤ .
7, with χ / d.o.f. = . / • H = . c ( H = . = . T ≤ .
6, with χ / d.o.f. = . / H = . H = . τ ′′ ) grows a little more slowly for H = . H = . T log( τ ′′ ) = A ( T − T ∗ ) c , (47)5 χ G R , H = 0.1 G L , H = 0.1 0 30 60 90 0.5 0.8 1.1 1.4 1.7 2 χ T G R , H = 0.3 G L , H = 0.3 FIG. 13. (Color online) Replicon and longitudinal susceptibilities inour cold annealings for H = . , . t w for each T ).The plot contains both equilibrium and non-equilibrium data. where | T ∗ | is very small but T ∗ could even be negative. Un-fortunately, we do not have enough degrees of freedom to fitsimultaneously for c and T ∗ . We shall discuss the possibleimplications of this T ∗ in the following section.Notice, finally, that in Figure 12 we can appreciate a sharpchange of regime precisely around the temperature where weidentified a power-law divergence of τ , fitting from the high-temperature phase.In the following section we shall introduce a more directstudy of ξ and try to combine the results of Sections IV and Vin a consistent physical picture. VI. DYNAMICS AND THE CORRELATION LENGTH
Up to now, we have focused on the determination of char-acteristic times and their temperature dependence. However, aproper discussion of any phase transition requires as well theconsideration of spatial correlation and the correlation length(this is, in fact, a longstanding obstacle in the investigation ofstructural glasses [5, 7, 8]). In the framework of spin glasseswe are advantaged, because the structure of correlators hasbeen investigated in detail (see Section II B and referencestherein). In particular, there are two types of correlation func-tions to deal with, the replicon and the longitudinal / anomalouscorrelator. We shall first decide which of the two correlatorsis worth studying and check that equilibrium results can in-deed be obtained in some temperature range (Section VI A).At that point, we shall revisit the the dynamics on the view ofthe correlation length (Section VI B). ξ t w T = 1.0, H = 0.2 T = 0.9, H = 0.2 FIG. 14. (Color online) Equilibration of the correlation length: for H = .
2, we show the ξ correlation length as computed from thereplicon correlator G R as a function of t w . We display data for twoconstant-temperature steps in the cold annealing (recall that t w is thetime elapsed since the last temperature drop). The constant tempera-ture time step is longer than τ only for T = . A. Which correlation function?
We compare in Figure 13 the replicon and longitudi-nal / anomalous susceptibilities [recall Eq. (16)], for H = . , . .
3. The susceptibilities are shown as a func-tion of temperature, as computed for the latest time on eachtemperature step along the cold annealing. This means thatFigure 13 contains both equilibrium and non-equilibrium data(depending on whether the constant-temperature step is muchlarger than τ , or not). In either case, it is rather obvious thatsignificant correlations appear only on the replicon correla-tor (in agreement with equilibrium, mean-field computations[24]). Therefore, we focus on the replicon correlator fromnow on. The anomalous sector is studied in more detail inAppendix B, Section B 3.We note as well that the failure of the longitudi-nal / anomalous correlator to display the correlations relevantto the problem might be related to analogous failures in ex-perimental investigations of structural glasses [1, 2].Specifically, we consider the ξ correlation length as com-puted from the replicon correlator G R using Eq. (18). Itstime evolution for two constant-temperature steps in the coldannealing run is displayed in Figure 14. For both tempera-tures we identify three regimes. For short t w the correlationlength basically remains constant (the fact that ξ does notdecrease at t w ∼ ff ects are weak). Then, the time evolutionstarts to be noticeable and ξ starts to increase. Finally, when t w ≫ τ , the correlation length becomes time-independent,which is consistent with our physical interpretation in Sec-tion IV that thermal equilibrium has been reached. We noteas well that the equilibrium regime is barely reachable for T = . τ ′′ ( T = . = . × , while τ ′′ ( T = . = . × . In the next paragraph, we shall6 τ '' ξ H = 0.1 H = 0.2 H = 0.3 H = 0, Critical dynamics FIG. 15. (Color online) Logarithmic plot of τ ′′ , computed in Sec-tion III, versus the correlation length ( ξ ) for the three magneticfields simulated. Data are in equilibrium. For comparison, we alsoshow the critical dynamics for H = τ ∼ ξ z , with z = .
86 [20]. discuss ξ as a function of temperature, but only for thosetemperatures where thermal equilibrium can be reached. B. Dynamics from the point of view of the correlation length
An important question is: is the dynamics activated [ τ ∼ e ξ Ψ / T , Eq. (24)], or critical [ τ ∼ ξ z , Eq. (26)]?As explained in Section II C, the droplet theory supports ac-tivated dynamics. On the other hand, the RSB theory is some-what ambiguous on this point [at mean-field level the dynam-ics is critical, this is the rationale for using z ν as the criticalexponent in Eq. (36)]. Furthermore, current theories for su-percooled liquid relaxations predict both types of behaviors(Section II D). These theories expect critical dynamics at hightemperatures, with an e ff ective exponent z MCT = γ [recallEqs. (29,30)]. However, at lower temperatures the dynamicsshould crossover to an activated behavior.At this point, we have in our hands equilibrium determina-tions for both τ ( T , H ) and ξ ( T , H ). Therefore, we can try toassess Eqs. (24,26) directly. This is attempted in Figure 15,where we used τ ′′ from fits to Eq. (35). Although this choiceis arbitrary to some extent, we recall that the critical diver-gence studied in Section IV turned out to be independent onthe choice of τ .We note two di ff erent regimes in Figure 15. For high-temperatures, data follow a critical dynamics. However, atlower temperatures (i.e., larger ξ ) τ starts to grow much fasterwith ξ . In fact, the e ff ective exponent z e ff = d log ξ/ d log τ be-comes as large as z e ff ≈
14, which clearly suggest that the dy-namics is becoming activated. Overall, this crossover exem-plifies the behavior expected for a supercooled liquid. In fact,critical dynamics is found in the range 10 < τ < , whichis also the range where MCT applies for simple supercooledliquids [76]. However, an alternative interpretation is possi-ble. First, one may note that we identified in Section IV A ξ T H = 0.1 H = 0.2 H = 0.3 FIG. 16. (Color online) Logarithmic plot of the equilibrium corre-lation length ξ versus the temperature for the three magnetic fieldssimulated. Lines are fits to ξ ( T ) = a H / T b H , see Eq. (48). When as-sessing the high-temperature data for H = .
1, recall that the criticaltemperature without a field is T H = = . an exponent γ ≈
2. Considering the large uncertainty in thisdetermination, this is consistent, via Eqs. (29) and (30), withthe z e ff ≈ z H = = .
86, the value forthe critical dynamics in the absence of a magnetic field [20].Hence, the crossover can be also due to the proximity of theRenormalization-Group fixed point at ( T c , H = H is the smaller the ξ needed to find activated dynamics(see Figure 15).Let us consider the temperature evolution of ξ in equilib-rium, see Figure 16. We do not find good fits to critical di-vergences, ξ ∝ / | T − T c ( H ) | ν , with T c ( H ) compatible withthe characteristic temperatures identified in Section IV. Forinstance, a fit for H = . χ value only for T c . .
5, while at H = . T c . . T c = ξ ( T ) = a H T b H , (48) b H = . = . . ≤ T ≤ . , χ / dof = . / . (49) b H = . = . . ≤ T ≤ . , χ / dof = . / . (50) b H = . = . . ≤ T ≤ . , χ / dof = . / . (51)Hence, at least within the temperature range that can be equi-librated, our data are compatible with a divergence of ξ onlyfor very low (perhaps even vanishing) T c ( H ). Notice, how-ever, that we are always working with ξ <
6, so this kind of fitis rather dangerous.At this point, is is only natural to ask whether the droplettheory, see Eq. (23), describes our data. The answer is nega-tive, see Figure 17.Therefore, none of the available theories provide a satis-factory description of our simulation. Of course, this mightbe due to the fact that we have not reached the regime where7 ξ H / ( d - θ ) T H θ /( d -2 θ ) θ = 0.61 0 0.4 0.8 1.2 1.6 0.7 1 1.3 1.6 ξ H / ( d - θ ) T H θ /( d -2 θ ) θ = 0.24 H = 0.1 H = 0.2 H = 0.3 FIG. 17. (Color online) Test of the droplet scaling law, Eq. (23).For each of our magnetic fields, we plot ξ ( T , H ) H / ( d − θ ) , as a func-tion of T H θ/ ( d − θ ) . For a proper choice of the sti ff ness exponent θ ,data should collapse in a single curve. We try (and fail within thereachable temperature window) to collapse the data with two di ff er-ent estimates of θ , see Section II C these theories apply (low enough temperatures, or low enoughmagnetic fields). However, we should stress that our data spana rather significant range of time scales (from one picosec-ond to a hundredth of a second). Hence, we dare say that oursimulations are of direct experimental relevance. The issue isdiscussed at length in the Conclusions.A final remark. One can be tempted to compare Eq. (48),which works for our equilibrium data, with the analysis inSection V B giving T log τ ∼ / T c (which is based on an ex-trapolation to times beyond our simulated timescales). Thiscomparison tells us that the droplet exponent Ψ ≈ c / b . There-fore, our data for H = . Ψ ≈ .
5, while our resultsfor H = . , . Ψ ≈
2. These values are ratherlarge, as compared to the value Ψ ∼ .
03 found in [14].
VII. CONCLUSIONS
In this work, we have investigated the approach to equilib-rium and the building-up of spin-glass order for the d = L =
80) than the correlation length, ξ ( t w ). Hence, we thinkthat our results are representative of the thermodynamic limit.Our time scales range from the picosecond to one hundredthof a second. We are thus approaching the experimental scale.However, when the temperature was low enough, we havebeen unable to reach the thermalization time scale, τ . We havemonitored this e ff ect carefully. Therefore, in this work weare presenting in a controlled way both equilibrium and non-equilibrium data. Our results have been analyzed on the lightof the two major theories on the market, the replica-symmetrybreaking and the droplet theory. On the view of recent claims[27, 28], we have also analyzed our data as suggested by cur-rent theories for relaxation in supercooled liquids. None ofthese three approaches was fully satisfying.The problem with the droplet / RSB theories was in the cor-relation length: the growth of ξ upon lowering the tempera-ture is too fast to fit the droplet theory and too slow to fit RSB.We summarize now the strengths and weaknesses of each ap-proach. We start with the droplet theory, then consider RSBand, finally, the supercooled liquids point of view.As we show in Section VI, the dynamics really seems tobe of activated type, as predicted by the droplet theory. How-ever, the scaling law predicted by the theory is not fulfilled byour data. In fact, see Section V B, the dynamics is of super-Arrhenius type at least down to temperature T = . T c = . H = equilibrium estimate of ξ does not seem todiverge at the de Almeida-Thouless line (also, the Fisher-Sompolinsky scaling[60] is not verified, as the reader mayeasily check). There are a number of possible explanationsfor our failure to find the divergence: • For all three magnetic fields, we have been able to equi-librate the system only down to T ≈ . T c ( H ). Perhapsthe critical growth of ξ ( T , H ) starts only closer to the deAlmeida-Thouless line. • It is by no means guaranteed that we are looking atthe right correlation function. We have shown in Sec-tion VI A that some correlators might display sizeablecorrelations while others do not. In fact, the quest forsensible correlators is a longstanding problem in the in-vestigation of supercooled liquids [5, 7, 8, 83–86]. Alsoin the field of spin glasses it has been suggested that en-ergy and link-overlap correlators deserve more attention[21, 87–89]. • As explained in Section II C, it is very possible thatthe physics in d = T =
0. If this is the case, activated dynamics is tobe expected also in the RSB theory. Under these cir-cumstances, the de Almeida-Thouless line identified in8Section IV and V A might well represent a dynamicglass transition. In fact, our data for ξ ( T , H ) are con-sistent with a critical divergence below the de Almeida-Thouless line. The divergence could take place in therange 0 ≤ T c ≤ . H = .
1, 0 ≤ T c ≤ . H = . ≤ T c ≤ . H = . • Finally, as explained above, it is possible that thedroplet theory could be correct and no transition takesplace for H > d =
3. Notice as well that the lower critical dimensionin a field is well below d =
4, since very clear evidenceof a second-order phase transition has been obtained in d = T d (Section IV A) and a cross-over to activateddynamics (Sects. IV A and VI B). However, the descriptionremains qualitative since our numerical accuracy does not al-low a strict test of basic relations among critical exponents.We note as well that the would-be Mode Coupling tempera-ture for H = . T d ≈ .
22 is rather large as compared tothe de Almeida-Thouless line T c ( H = . ≈ . T d ≈ . T g ( T g is the dynamic glasstemperature where τ ∼ ACKNOWLEDGMENTS
The Janus project has been partially supported by the EU(FEDER funds, No. UNZA05-33-003, MEC-DGA, Spain);by the European Research Council under the EuropeanUnion’s Seventh Framework Programme (FP7 / E , W , q t w E (G.) E ( n =2) E ( n =5) W (G.) W (n=2) W (n=5) q (G.) q ( n =2) q ( n =5) FIG. 18. (Color online) Thermal energy E ( t w ), magnetic energy W ( t w ) and overlap q ( t w ) as a function of time for L = T = . H = .
3. We have plotted the results from a full Gaussian distri-bution (G.), as well as discretizations with n = n =
5. Noticethat all three simulations have produced the same equilibrium valuesfor these observables.
Appendix A: Discretization of the Gaussian Magnetic Field
In this appendix we describe the procedure we have usedto discretize (using the Gauss-Hermite quadrature [58]) theGaussian magnetic field in order to implement our simula-tions on the Janus dedicated computer (which cannot handlenon-integer arithmetic e ffi ciently). This implementation wasintroduced (but not explained in detail) in [90].The Gauss-Hermite quadrature formula can be interpretedas an approximation formula for probability distributions . Infact, if we multiply either of the two distributions related inEq. (A1), below, by an arbitrary polynomial of order 4 n −
1, and integrate this product through x ∈ ( −∞ , ∞ ), identicalresults are obtained [58]: e − x ≈ n X k = w k δ ( x − x k ) + δ ( x + x k )] . (A1)where x k are the positive zeros of the 2 n -th Hermite polyno-mial H n ( x ) and the weights, w k , are given by w k = n − (2 n )! √ π (2 n ) H n − ( x k ) . (A2)In our numerical simulations we need to compute integralslike O ≡ H √ π Z ∞−∞ dh O ( h ) e − h H (A3)Furthermore, the above equation can be further simplified us-ing the gauge symmetry (2), which allows us to consider onlypositive magnetic fields, so O = H √ π Z ∞ dh O ( h ) e − h H (A4)9and using Eq. (A1), one can finally write O = √ π n X k = w k O ( √ Hx k ) . (A5)We shall limit ourselves to n = w = w / √ π the field is h = √ Hx (itis set to h = √ Hx otherwise). Note that two bits per lat-tice site are enough to code this n = n = x = . x = . w = . w = . n = L = n = H ≤ .
3. Inparticular, we have compared the results of simulations with n = n = W (see Figure 18). We have checked that the di ff er-ences between the observables computed at finite n and thosecomputed with full Gaussian magnetic field are statisticallycompatible with zero.Another strong test of our implementation is the agree-ment in the asymptotic values of q ( t w ) and W ( t w ) in the high-temperature phase (see Figure 2).In conclusion: we have checked that, for the main quantitiesconsidered in this work (computed with n = Appendix B: Spatial correlation functions
In this appendix we discuss three di ff erent features of thespatial correlation functions of the d = direct quenches (this anomaly seems to be absent fromannealing protocols), and (iii) the comparison of the repliconand the longitudinal / anomalous correlators.
1. Long-distance behavior of the replicon propagator
As explained in Section II B, when computing the integrals I k , see Eq. (19), it is crucial to impose a long distance cuto ff .Otherwise, the I k integrals become non-self-averaging objectsthat can be accurately computed only with a huge numberof samples. However, one may try to correct the systematice ff ects induced by the cuto ff by studying the long-distancebehavior of the propagator G R . One fits the curve to a suit-able, simple functional form and then computes by hand theremaining part of the integral. The contribution to I k from r > r cuto ff is usually tiny, but we prefer to monitor it. This -5 -4 -3 -2 -1
0 5 10 15 20 25 G R ( r , t w ) r T = 0.90 A exp(- r / ξ ) B exp(- r / ξ ’) / r FIG. 19. (Color online) Semilogarithmic plot of the equilibriumreplicon correlator G R as a function of distance r , for temperature T = . H = . r , , G R ( r ) ≈ A e − r /ξ ( A and ξ arefitting parameters), obtained for r ≥
15 and to G R ≈ B e − r /ξ ′ / r for r ≥
10. These two functional forms are indistinguishable for large r . issue has been discussed at length in [19, 20], where the spinglass without a field was studied.Here, we show in Figure 19 that, in the temperature regimewhere we manage to equilibrate the system in a field, G R ( r ≫ ξ ) decays exponentially. Therefore, estimating the contribu-tion from r > r cuto ff to the integrals I k is fairly easy (the result-ing correction is smaller than the error bar). We can includean algebraic prefactor in the fitting function the better to fitthe small- r sector, but this is irrelevant for the tails (see Fig-ure 19).
2. Overshooting in the direct quench
The direct quench is an idealized temperature-variation pro-tocol: one takes a fully disordered spin glass (i.e.. T = ∞ )and places it instantaneously at the working temperature. Itis clear that, in the laboratory, temperature should vary gradu-ally. Hence, the annealing protocols described in Section III Aare closer to the temperature variations that one can realizeexperimentally. On the other hand, the direct quench is thesimplest protocol in a computer simulation.In fact, we have found with some surprise that the non-equilibrium behavior of the replicon correlator is rather dif-ferent in a direct quench and in an annealing protocol. In Fig-ure 20 we show the time evolution of the correlation length ξ ( t w ) for two temperatures in our hot annealing, the initialone T = .
2, and the second temperature T = .
15. Note thatthe time evolution at the very first temperature in the annealingcan be aptly described as a direct quench. Indeed, in Figure 20we notice an overshooting of ξ ( t w ) in the direct quench: wellbefore equilibrium is reached, a maximum is found which liesabove the equilibrium correlation length. No such maximumarises in the lower temperatures of the annealing. We have0 ξ ( t w ) t w T = 1.20 T = 1.15 FIG. 20. (Color online) Time evolution of the correlation length ξ as computed for the replicon correlator for temperatures T = . T = .
15 in our hot annealing, see Section III A. Note that thetemperature step T = . T = .
15 already belongs to the annealing part of the run. The max-imum for T = .
2, where ξ ( t w ) is larger than its equilibrium value,seems to be a generic feature of any direct quench for a spin glass ina field in three dimensions. On the other hand, during the annealing,the time evolution of ξ ( t w ) is monotonic. -5 -4 -3 -2 -1
0 2 4 6 8 10 12 14 16 18 G R ( r , t w ) rt w = 339958, T = 1.20 t w = 42469097, T = 1.20 FIG. 21. (Color online) Comparison of the non-equilibrium correla-tor G R ( r ; t w ) with the corresponding equilibrium value for the initialtemperature, T = .
2, in our hot annealing (i.e., in a direct quenchrun). The time t w corresponds to the maximum correlation length inFigure 20. checked that this overshooting is characteristic of the directquench, as it happens basically for all temperatures and mag-netic fields.We can look at this overshooting in greater detail in Fig-ure 21, where we compare the equilibrium G R ( r ) with the non-equilibrium G R ( r ; t w ) at the t w corresponding to the maximumin Figure 20. The two correlators are remarkably featurelessas a function of r , but the overshooting e ff ect is also clear from G R ( r ).
3. The anomalous / longitudinal sector The longitudinal / anomalous correlator defined in Eq. (15)appears naturally in the analysis of the mean-field approx-imation [24]. To the best of our knowledge, the longitudi-nal / anomalous correlator has not been studied in three spatialdimensions, in the presence of a field. We recall that fromthese correlators, one may obtain associated susceptibilities,see Eq. (16).We showed in Figure 13 that the replicon susceptibilitygrows significantly upon lowering the temperature, while thelongitudinal / anomalous susceptibility does not. However,when looking at the plot of χ , which is a spatial integral of G , we are losing information on the shape of the correlationfunction.Here we perform a more detailed comparison of both corre-lators by studying their ratio as a function of r in Fig. 22, left. G L / G R was computed at the longest time available, i.e., asclose as possible to thermal equilibrium. We identify two dif-ferent regimes, at high and low temperatures. At high temper-atures, G L / G R decreases exponentially in r , see the right panelin Fig. 22. In fact, barring unavoidable di ff erences on the al-gebraic prefactors, G L ( r ) / G R ( r ) ∝ exp (cid:2) − r (cid:0) ξ L − ξ R (cid:1)(cid:3) . Hencethe exponential decrease in Fig. 22–right implies ξ R > ξ L . Onthe other hand, at the lowest temperatures that we reach, i.e., T = . . G L / G R becomes essentially constant at large r , suggesting that the correlation length is the same for bothcorrelators. This is quite surprising: if a de Almeida-Thoulessline exists one expects 0 < ξ L /ξ R < T = G L ( r , t tot ) G R ( r , t tot ) = G L ( r , t tot = ∞ ) G R ( r , t tot = ∞ ) + A ( r ) t x tot . (B1)In the above equation, the exponent x was allowed to dependon r and T (we found that it barely depended on r for a giventemperature). We were able to carry out this extrapolationsafely down to temperature T = .
0, see right panel in Fig. 23.In the limit of long times, the ratio of propagators does de-crease exponentially with r , which confirms that ξ R > ξ L (atleast down to temperature T = .
0, for H = . H = . G L / G R are not relevant).1 G L ( r ) / G R ( r ) r T = 1.8 T = 1.6 T = 1.4 T = 1.2 T = 1.0 T = 0.8 T = 0.6 T = 0.4 0.1 1 0 5 10 15 G L ( r ) / G R ( r ) r T = 1.8 T = 1.6 T = 1.4 T = 1.2 T = 1.0 FIG. 22. (Color online) Comparison of the replicon and the longitudinal propagator for H = .
1. Data from the longest t w at each temperaturestep, in our slowest annealing for each temperature. Notational conventions are as in Figure 19. The left panel shows the quotient G L / G R for the whole temperature range in a linear scale, while the right panel shows the same quantity in a logarithmic scale (removing the lowesttemperatures). In all cases we cut each graph at the point where the statistical error becomes greater than 50%. For low T , the quotient seemsto reach a non-zero asymptotic value. For high temperatures, it seems to decay exponentially, although the decay slows down for large r ,suggesting a possible non-zero plateau for large r . G L ( r ) / G R ( r ) t tot T = 1.0, r = 2 T = 1.0, r = 4 0.01 0.1 1 0 2 4 6 8 10 12 14 G L ( r ) / G R ( r ) r T = 1.4 T = 1.2 T = 1.0 FIG. 23. (Color online)
Left:
We plot the quotient G L / G R for T = . H = . t tot for two di ff erent values of r . It is clear that even our longest times have small thermalisation e ff ects. The long- t tot limit may be estimated with a power-law fit (continuouslines). Right: G L / G R , extrapolated to infinite t tot for several temperatures. The curvature observed in Figure 22 has disappeared and a pureexponential decay is now apparent.[1] P. G. Debenedetti, Metastable Liquids (Princeton UniversityPress, Princeton, 1997).[2] P. G. Debenedetti and F. H. Stillinger, Nature , 259 (2001).[3] A. Cavagna, Physics Reports , 51 (2009), arXiv:0903.4264.[4] J. A. Mydosh,
Spin Glasses: an Experimental Introduction (Taylor and Francis, London, 1993).[5] A. Montanari and G. Semerjian, J. Stat. Phys. , 23 (2006),arXiv:cond-mat / , 139 (1965).[7] E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, andD. A. Weitz, Science , 627 (2000).[8] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D. El Masri,D. L’Hôte, F. Ladieu, and M. Pierno, Science , 1797 (2005).[9] K. Gunnarsson, P. Svedlindh, P. Nordblad, L. Lundgren,H. Aruga, and A. Ito, Phys. Rev. B , 8199 (1991). [10] M. Palassini and S. Caracciolo,Phys. Rev. Lett. , 5128 (1999), arXiv:cond-mat / , 14237 (2000),arXiv:cond-mat / , 257202 (2002),arXiv:cond-mat / , 438 (1999).[14] F. Bert, V. Dupuis, E. Vincent, J. Hammann, and J.-P.Bouchaud, Phys. Rev. Lett. , 167203 (2004).[15] A. Cruz, J. Pech, A. Tarancon, P. Tellez, C. L. Ul-lod, and C. Ungil, Comp. Phys. Comm , 165 (2001),arXiv:cond-mat / [16] A. T. Ogielski, Phys. Rev. B , 7384 (1985).[17] F. Belletti, F. Mantovani, G. Poli, S. F. Schifano, R. Tripic-cione, I. Campos, A. Cruz, D. Navarro, S. Perez-Gaviro,D. Sciretti, A. Tarancon, J. L. Velasco, P. Tellez, L. A. Fernan-dez, V. Martin-Mayor, A. Muñoz Sudupe, S. Jimenez, A. Maio-rano, E. Marinari, and J. J. Ruiz-Lorenzo (Janus Collabora-tion), Computing in Science and Engineering , 41 (2006).[18] F. Belletti, M. Cotallo, A. Cruz, L. A. Fernandez, A. Gordillo,A. Maiorano, F. Mantovani, E. Marinari, V. Martin-Mayor,A. Muñoz Sudupe, D. Navarro, S. Perez-Gaviro, J. J.Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Tarancon,R. Tripiccione, and J. L. Velasco (Janus Collaboration),Comp. Phys. Comm. , 208 (2008), arXiv:0704.3573.[19] F. Belletti, M. Cotallo, A. Cruz, L. A. Fernandez,A. Gordillo-Guerrero, M. Guidetti, A. Maiorano, F. Man-tovani, E. Marinari, V. Martin-Mayor, A. M. Sudupe,D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo,S. F. Schifano, D. Sciretti, A. Tarancon, R. Tripiccione,J. L. Velasco, and D. Yllanes (Janus Collaboration),Phys. Rev. Lett. , 157201 (2008), arXiv:0804.1471.[20] F. Belletti, A. Cruz, L. A. Fernandez, A. Gordillo-Guerrero,M. Guidetti, A. Maiorano, F. Mantovani, E. Marinari,V. Martin-Mayor, J. Monforte, A. Muñoz Sudupe, D. Navarro,G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schi-fano, D. Sciretti, A. Tarancon, R. Tripiccione, and D. Yl-lanes (Janus Collaboration), J. Stat. Phys. , 1121 (2009),arXiv:0811.2864.[21] R. Alvarez Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano,F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Taran-con, R. Tripiccione, and D. Yllanes (Janus Collaboration),J. Stat. Mech. , P06026 (2010), arXiv:1003.2569.[22] R. Alvarez Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maiorano,F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Taran-con, R. Tripiccione, and D. Yllanes (Janus Collaboration),Phys. Rev. Lett. , 177202 (2010), arXiv:1003.2943.[23] S. Franz, M. Mézard, G. Parisi, and L. Peliti,Phys. Rev. Lett. , 1758 (1998).[24] J. R. L. de Almeida and D. J. Thouless,J. Phys. A , 983 (1978).[25] M. Mézard, G. Parisi, and M. Virasoro, Spin-Glass Theory andBeyond (World Scientific, Singapore, 1987).[26] E. Marinari, G. Parisi, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, and F. Zuliani, J. Stat. Phys. , 973 (2000),arXiv:cond-mat / , 217202 (2002), arXiv:cond-mat / , 5388 (1987).[30] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Phys. Rev.A , 1045 (1989).[31] W. Götze and T. F. Sjögren, Rep. Prog. Phys. , 241 (1992).[32] M. Mézard and G. Parisi, Phys. Rev. Lett. , 747 (1999).[33] M. Mézard and G. Parisi, J. Chem. Phys. , 1076 (1999).[34] B. Coluzzi, G. Parisi, and P. Verrocchio,Phys. Rev. Lett. , 306 (2000).[35] G. Parisi and T. Rizzo, J. Phys. A , 235003 (2010).[36] S. Singh, M. Ediger, and J. de Pablo, Nature Mater. , 139 (2013).[37] W. L. McMillan, J. Phys. C: Solid State Phys. , 3179 (1984).[38] A. J. Bray and M. A. Moore, in Heidelberg Colloquium onGlassy Dynamics , Lecture Notes in Physics No. 275, edited byJ. L. van Hemmen and I. Morgenstern (Springer, Berlin, 1987).[39] D. S. Fisher and D. A. Huse, Phys. Rev. Lett. , 1601 (1986).[40] D. S. Fisher and D. A. Huse, Phys. Rev. B , 386 (1988).[41] M. A. Moore and A. J. Bray, Phys. Rev. B , 224408 (2011),arXiv:1102.1675.[42] G. Parisi and T. Temesvári, Nucl. Phys. B , 293 (2012),arXiv:1111.3313.[43] J. Yeo and M. A. Moore, Phys. Rev. E , 052501 (2012),arXiv:1208.3044.[44] A. P. Young and H. G. Katzgraber,Phys. Rev. Lett. , 207203 (2004), arXiv:cond-mat / , 197202 (2008), arXiv:0712.2009.[46] P. E. Jönsson, H. Takayama, H. A. Katori, and A. Ito,Phys. Rev. B , 180412(R) (2005), arXiv:cond-mat / , 5130 (1999), arXiv:cond-mat / , 207206 (2002), arXiv:cond-mat / , 123704 (2010), arXiv:1009.6115.[50] R. A. Baños, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, D. Iniguez,A. Maiorano, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane,A. Tarancon, P. Tellez, R. Tripiccione, and D. Yllanes,Proc. Natl. Acad. Sci. USA , 6452 (2012).[51] G. Kotliar, P. W. Anderson, and D. L. Stein,Phys. Rev. B , 602 (1983).[52] D. Larson, H. G. Katzgraber, M. A. Moore, and A. P. Young,Phys. Rev. B , 024414 (2013), arXiv:1211.7297.[53] L. Leuzzi and G. Parisi, Phys. Rev. B , 224204 (2013),arXiv:1303.6333.[54] C. de Dominicis and I. Giardina, Random Fields and SpinGlasses (Cambridge University Press, Cambridge, England,2006).[55] G. Parisi,
Field Theory, Disorder and Simulations (World Sci-entific, 1994).[56] F. Belletti, M. Guidetti, A. Maiorano, F. Mantovani, S. F. Schi-fano, R. Tripiccione, M. Cotallo, S. Perez-Gaviro, D. Sciretti,J. L. Velasco, A. Cruz, D. Navarro, A. Tarancon, L. A.Fernandez, V. Martin-Mayor, A. Muñoz-Sudupe, D. Yl-lanes, A. Gordillo-Guerrero, J. J. Ruiz-Lorenzo, E. Marinari,G. Parisi, M. Rossi, and G. Zanier (Janus Collaboration),Computing in Science and Engineering , 48 (2009).[57] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernan-dez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti,D. Iniguez, A. Maiorano, F. Mantovani, E. Marinari,V. Martin-Mayor, J. Monforte-Garcia, A. Munoz Sudupe,D. Navarro, G. Parisi, M. Pivanti, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane,A. Tarancon, P. Tellez, R. Tripiccione, and D. Yllanes,Eur. Phys. J. Special Topics , 33 (2012), arXiv:1204.4134.[58] M. Abramowitz and I. A. Stegun, Handbook of MathematicalFunctions: with Formulas, Graphs, and Mathematical Tables. ,ninth ed. (Dover Publications, New York, 1972).[59] A. J. Bray and S. A. Roberts, J. Phys. C: Solid St.Phys. , 5405(1980). [60] D. S. Fisher and H. Sompolinsky,Phys. Rev. Lett. , 1063 (1985).[61] T. Temesvári and C. De Dominicis,Phys. Rev. Lett. , 097204 (2002), arXiv:cond-mat / , 220401 (2008).[63] A. K. Hartmann, Phys. Rev. E. , 84 (1999),arXiv:cond-mat / , 5126 (1999).[65] S. Boettcher, Eur. Phys. J. B , 83 (2004),arXiv:cond-mat / , 3017 (2000).[67] F. Krzakala and O. C. Martin, Phys. Rev. Lett. , 3013 (2000).[68] D. Yllanes, Rugged Free-Energy Landscapes in DisorderedSpin Systems (Ph.D. thesis, UCM, 2011) arXiv:1111.0266.[69] H. Sompolinsky and A. Zippelius, Phys. Rev. B , 6860(1982).[70] T. Nattermann, in Spin glasses and random fields , edited byA. P. Young (World Scientific, Singapore, 1998).[71] W. Götze,
Complex dynamics of glass-forming liquids: A mode-coupling theory (Oxford University Press, USA, 2009).[72] A. Andreanov, G. Biroli, and J.-P. Bouchaud, Europhysics Let-ters , 16001 (2009).[73] S. Franz, H. Jacquin, G. Parisi, P. Urbani, and F. Zamponi, J.Chem. Phys. , 12A540 (2012), arXiv:1210.7304.[74] S. Franz, H. Jacquin, G. Parisi, P. Urbani, and F. Zam-poni, Proc. Natl. Acad. Sci. USA , 18725 (2012),arXiv:1206.2482.[75] F. Caltagirone, U. Ferrari, L. Leuzzi, G. Parisi, F. Ricci-Tersenghi, and T. Rizzo, Phys. Rev. Lett. , 085702 (2012),arXiv:1111.6420.[76] W. Kob and H. C. Andersen, Phys. Rev. Lett. , 1376 (1994).[77] V. Lubchenko and P. G. Wolynes, Annu. Rev. Phys. Chem. ,235 (2007). [78] T. Jonsson, K. Jonason, P. E. Jönsson, and P. Nordblad,Phys. Rev. B , 8770 (1999).[79] G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo,Phys. Rev. B , 13617 (1998), arXiv:cond-mat / , 1181 (1998),arXiv:cond-mat / x c can be expressed as: x c = ( d − + η ) / (2 z ), where d is the di-mensionality of the system, z is the dynamical critical exponentand η is the anomalous dimension.[82] M. Baity-Jesi, R. A. Baños, A. Cruz, L. A. Fernandez, J. M.Gil-Narvion, A. Gordillo-Guerrero, D. Iniguez, A. Maiorano,F. Mantovani, E. Marinari, V. Martin-Mayor, J. Monforte-Garcia, A. Muñoz Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, M. Pivanti, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F.Schifano, B. Seoane, A. Tarancon, R. Tripiccione, and D. Yl-lanes (Janus Collaboration), Phys. Rev. B , 224416 (2013),arXiv:1310.2910.[83] S. Franz and G. Parisi, J. of Physics: Condensed Matter ,6335 (2000).[84] T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. A , 4439(1988).[85] A. Cavagna, T. S. Grigera, and P. Verrocchio, Phy. Rev. Lett. , 187801 (2007), arXiv:cond-mat / , 771 (2008), arXiv:cond-mat / , 14852 (1998).[88] P. Contucci and C. Giardinà, Phys. Rev. B , 014456 (2005).[89] P. Contucci, C. Giardinà, C. Giberti, and C. Vernia,Phys. Rev. Lett. , 217204 (2006).[90] L. Leuzzi, G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo,Phys. Rev. Lett.103