VOLUME 88, NUMBER 23 P H Y S I C A L R E V I E W L E T T E R S 10 JUNE 2002

237201-1

Heterogeneous Aging in Spin Glasses

Horacio E. Castillo,1 Claudio Chamon,1 Leticia F. Cugliandolo,2,3 and Malcolm P. Kennett41Physics Department, Boston University, Boston, Massachusetts 02215

2 Laboratoire de Physique Théorique de l’École Normale Supérieure, Paris, France3 Laboratoire de Physique Théorique et Hautes Energies, Jussieu, Paris, France

4Department of Physics, Princeton University, Princeton, New Jersey 08544(Received 13 December 2001; published 20 May 2002)

We introduce a set of theoretical ideas that form the basis for an analytical framework capable ofdescribing nonequilibrium dynamics in glassy systems. We test the resulting scenario by comparing itspredictions with numerical simulations of short-range spin glasses. Local fluctuations and responses areshown to be connected by a generalized local out-of-equilibrium fluctuation-dissipation relation. Scalingrelationships are uncovered for the slow evolution of heterogeneities at all time scales.

DOI: 10.1103/PhysRevLett.88.237201 PACS numbers: 75.10.Nr, 75.10.Hk, 75.10.Jm, 75.50.Lk

Very slow equilibration and sluggish dynamics are char-acteristics shared by disordered spin systems and by otherglassy systems such as structural and polymeric glasses.The origin of this dynamic arrest near and below the glasstransition is currently poorly understood. Studies of thetime evolution of many quantities, such as the remanentmagnetization, the dielectric constant, or the incoherentcorrelation function, have shown that below the glass tran-sition the system falls out of equilibrium [1,2]. This isevidenced by the presence of aging, i.e., the dependenceof physical properties on the time since the quench intothe glassy state, and also by the breakdown of the equilib-rium relations dictated by the fluctuation-dissipation theo-rem (FDT) [3,4].

Most analytical progress in understanding nonequilib-rium glassy dynamics has been achieved in mean-fieldfully connected spin models [3], while numerical simula-tions have addressed both structural glasses [5] and short-range spin glass models [6]. Until recently, however,experimental, numerical, and analytical studies havemainly focused on global quantities, such as globalcorrelations and responses, which do not directly probelocal relaxation mechanisms. Local regions that behavedifferently from the bulk, or dynamic heterogeneities,could be crucial to understand the full temporal evolu-tion, and have received considerable experimental [7–9]and numerical [10] attention lately. However, no cleartheoretical picture has yet emerged to describe the localnonequilibrium dynamics of the glassy phase.

Here we introduce such a theoretical framework, andtest its predictions via numerical simulations of a short-range spin glass model. We show that local correlationsand responses are linked, and we find scaling propertiesfor the heterogeneities that connect the evolution of thesystem at different times. This universality may providea general basis for a realistic physical understanding ofglassy dynamics in a wide range of systems.

The framework that we propose is motivated by ananalogy [11] between aging dynamics and the well-known

0031-9007�02�88(23)�237201(4)$20.00

statics of Heisenberg magnets. For concreteness, wetest its predictions against Monte Carlo simulations onthe prototypical spin glass model, the three-dimensionalEdwards-Anderson (3DEA) model, H �

P�ij� Jijsisj ,

where si � 61 and the nearest-neighbor couplings areJij � 61 with equal probability. We argue that twodynamical local quantities, the coarse-grained localcorrelation Cr�t, tw � � 1

V

Pi[Vr

si�t�si�tw� and integrated

response xr�t, tw � � 1Nf

PNf

k�11V

Pi[Vr

si�t�jh�k� 2si�t�h

�k�i

, areessential to understand the mechanisms controlling thedynamics of glassy systems. The spins are representedby si in the absence of an applied field and by sijh�k� inthe presence of an applied field. si�t� � 1

t

Pt0�t21t 0�t2t si�t0�

is the result of coarse graining the spin over a small timewindow [typically, t � 1000 Monte Carlo steps (MCS)].Vr is a cubic box with volume V centered at the point r.By taking V to be the volume of the whole system, thebulk or global correlation C�t, tw� and response x�t, tw �are recovered. Two generic times after preparation arerepresented by tw and t, with tw # t. When the systemis not in equilibrium, time dependences do not reduce toa dependence on the time difference t 2 tw. We measurea staggered local integrated linear response by applying a

bimodal random field on each site h�k�i � 6h during the

time interval �tw, t�. Linear response holds for the valuesof h that we use. The index k � 1, . . . , Nf labels differentrealizations of the perturbing field. We use random initialconditions. The thermal histories, i.e., the sequences ofspins and random numbers used in the MC test, are thesame with and without a perturbing field.

In a disordered spin model, the coarse-grained localmagnetization typically vanishes, but the local correlationis nontrivial. Averaged over disorder and the thermal his-tory, this correlator defines the Edwards-Anderson parame-ter qEA when tw ! `, and t 2 tw ! ` subsequently.Can we detect the growth of local order [12] by analyz-ing the evolution of the local correlator, as one easily canfor a system undergoing ferromagnetic domain growth? In

© 2002 The American Physical Society 237201-1

VOLUME 88, NUMBER 23 P H Y S I C A L R E V I E W L E T T E R S 10 JUNE 2002

rtw=410000 MCs

0.75 0.5 0.25 0

-0.25 -0.5

0 16 32 48 64x 016

3248

64

y

−1−0.5

00.5

1

C

FIG. 1. The local correlations on a 2D cut of the cubic cell.The linear size is L � 64. T � 0.72Tc , V � 33, tw � 4.1 3105 MCS, and t � 2.8 3 106 MCS.

Fig. 1 we show the local correlation for fixed tw and t on a2D cut of the 3DEA model. Regions with large values ofCr are intertwined with regions with a small value of Cr asshown by the contour levels. This behavior persists for alltw and t that we can reach with the simulation, and a moresophisticated analysis is necessary to identify a growingorder in this system.

It is clear from Fig. 1 that different sites have distinctdynamics. Analysis of the local correlation for fixed tw

as a function of t shows that in general the relaxation isnonexponential; this is often ascribed to the presence ofheterogeneous dynamics. How can one characterize theheterogeneous dynamics and determine its origin? We ar-gue that relevant quantities are the probability distributionfunction (PDF) of the local correlation, r���Cr �t, tw����, thePDF of the local integrated response, r���xr�t, tw����, andthe joint PDF r���Cr �t, tw�, xr �t, tw ����, and we start by dis-cussing the latter.

The FDT relates the correlation of spontaneous fluctua-tions to the integrated linear response of a chosen observ-able [e.g., Cr�t, tw� and xr �t, tw� averaged over thermalhistories] at equilibrium. Glassy systems modify theFDT in a particular way first obtained analytically formean-field models [3], later verified numerically in anumber of realistic models [13,14], and more recentlytested experimentally [15]. A parametric plot of the bulkintegrated response, x�t, tw�, against the bulk correlation,C�t, tw�, for fixed and long waiting time tw and using t as aparameter, approaches a nontrivial limit, x�C�, repre-sented by the crosses in Fig. 2(b). The curve has a straightsection, for which t 2 tw ø tw and the correlation decaysfrom 1 to qEA with a slope 21�T as found with the equilib-rium FDT. Beyond this point, as t increases towards infin-ity, the curve separates from the FDT line. Now, considereach lattice site for fixed times tw and t. If we plot pointsfor the pairs �Cr �t, tw�, xr �t, tw ��, where will they lie?

When tw ! ` and t 2 tw ø tw , all local correlatorssatisfy the FDT strictly once averaged over thermal histo-ries, since the magnitude of local deviations from the FDThas an upper bound [16]. We have checked that Cr�t, tw�

237201-2

and xr �t, tw� obey the FDT for an individual thermal his-tory apart from small fluctuations [see Fig. 2(b)].

For the regime of widely separated times we proposean analysis similar in spirit to the one that applies to thelow energy excitations of the Heisenberg model. There,the free energy for the coarse-grained magnetization �m� �r�is F �

Rddr ���� �=�r �m� �r����2 1 V ���j �m��r�j����. A spontaneous

symmetry breaking signals the transition into the orderedphase � �m� � �m0 fi 0, in which the order parameter hasboth a uniform length [the radius of the bottom of theeffective potential V�j �mj�], and a uniform direction. F isinvariant under uniform rotations �m��r� ! R �m��r�. Thelowest energy excitations (spin waves) are obtained fromthe ground state by leaving the length of the vector invari-ant and applying a slowly varying rotation to it: �m��r� �R��r� �m0. These are massless transverse fluctuations(Goldstone modes). In contrast, longitudinal fluctuations,which change the magnitude of the magnetization vector,are massive and energetically costly.

Let us now apply the same kind of analysis to the dy-namics of the spin glass. Here, the relevant fluctuatingquantities are the coarse-grained local correlations Cr andtheir associated local integrated responses xr . In Ref. [11]we derived an effective action for these functions thatbecomes invariant under a global time reparametrizationt ! h�t� in the aging regime. This symmetry leaves thebulk relation, x�C�, invariant. A uniform reparametriza-tion is analogous to a global rotation in the Heisenbergmagnet, and the curve x�C� is analogous to the surfacewhere V �j �mj� is minimized. Hence, we expect that, forfixed long times tw and t in the aging regime, the lo-cal fluctuations in Cr and xr should be given by smoothspatial variations in the time reparametrization, hr �t�, i.e.,Cr�t, tw� � CSP���hr �t�, hr �tw���� C���hr �t�, hr �tw����, whereCSP is the global correlation at the saddle-point level whichin the numerical studies we approximate by the actualglobal correlation C, and similarly for xr . These trans-verse fluctuations are soft Goldstone modes. Longitudinalfluctuations, which move away from the x�C� curve, aremassive and penalized. This implies the first testable pre-diction of our theoretical framework: the pairs �Cr , xr �should follow the curve x�C� for the bulk integrated re-sponse against the bulk correlation.

In Fig. 2, we test this prediction by plotting the distri-bution of pairs �Cr , xr �. We find, as expected, that forlong times the dispersion in the longitudinal direction [i.e.,away from the bulk x�C� curve] is much weaker than inthe transverse direction [i.e., along the bulk x�C� curve].In the coarse-grained aging limit we expect the former todisappear while the latter should remain. (This limit corre-sponds to the way actual measurements are performed: thethermodynamic limit is taken first to eliminate finite sizeeffects and undesired equilibration; then the large tw limitis taken to reach the asymptotic regime, and, finally, thelimit V ! ` serves to eliminate fluctuations through thecoarse-graining process; in the figure we used a large vol-ume V � 133 to approach the latter limit though we found

237201-2

VOLUME 88, NUMBER 23 P H Y S I C A L R E V I E W L E T T E R S 10 JUNE 2002

0

0.5Cr 0

0.5

1

χr

5

15

25

ρ

(a)

0

0.5Cr 0

0.5

1

χr

5

15

25

ρ

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

++

++++

+ Bulk

FDT(b)

0 0.5 1Cr

0

0.5

1

χr

FIG. 2. L � 64, T � 0.72Tc , V � 133, tw � 4 3 104. (a) Surface plot of the joint PDF for t�tw � 16. Points within the showncontour account for 2�3 of the total probability. (b) Time evolution of r���Cr�t, tw�, xr�t, tw����, for t�tw � 1.00005, 1.001, 1.06, 2, 8,32 (from right to left). Each contour contains 2�3 of the total probability for the given pair of times. The line marked FDT is theequilibrium relation between xr and Cr . The crosses are the bulk x�C� for the same pairs of times.

a similar qualitative behavior for smaller V .) Figure 2(a)displays the joint PDF r�Cr , xr � for a pair of times �tw , t�that are far away from each other. Figure 2(b) showsthe projection of a set of contour levels for tw fixed andsix values of t. Even though the data for each contourcorrespond to a single pair of times �tw , t�, the fluctua-tions span a range of values that, for the bulk quantities,would require a whole family of pairs �tw, t�. This re-veals that the aging process is nonuniform across a finite-range model.

We now turn our attention to a more detailed exami-nation of the time dependences in the local dynamics. Agood t�tw scaling that breaks down only for very largevalues of the subsequent time t 2 tw has been obtainedfor the bulk thermoremanent magnetization experimentally[2] and for the bulk correlation numerically [6] once thestationary (t 2 tw ø tw) part of the relaxation is sub-tracted, as suggested by the solution to mean-field models[3]. For systems that display this particular dependenceon t�tw for the bulk correlator, a second prediction canbe extracted from our theoretical framework: the distri-bution r���Cr �t, tw���� should depend only on the ratio t�tw .Even further, if the bulk correlator has a simple power lawform CSP�t, tw� qEA�t�tw�2r , an approximate treatmentof fluctuations leads to a rescaling and collapse of r���Cr �t,tw���� even for pairs of times with different ratios t�tw .

Since we are dealing with ratios of times, it is convenientto define hr �t� � ewr�t�, so that Cr �t, tw � � CSP���hr �t��hr �tw���� � CSP�ewr �t�2wr�tw��. Therefore the statistics oflocal correlations are determined from the statistical dis-tribution of distances between two “surfaces,” wr �t� 2

wr�tw�. In this form, a dynamic theory of short-range spinglasses is not different from a theory of fluctuating geome-tries or elasticity. We propose a simple reparametrizationinvariant effective action for wr�t� � lnt 1 dwr�t�, ex-panding around dwr�t� � 0, with no zeroth or first or-der term in dwr�t�. We ensure that the effective action isreparametrization invariant by taking one time derivativefor each time variable. Thus [17]

237201-3

S �qEAr

2

Zddr

Z `

0dt

Z `

0dt 0 = �wr�t�= �wr�t0�

3 e2rjwr �t�2wr�t0�j,

where the last factor penalizes fast time variations ofwr and the = ensures that spatial variations are smooth.Expanding to lowest order in dwr�t� yields wr�t� 2wr�tw� � ln�t�tw� 1 dwr�t� 2 dwr�tw� � ln�t�tw� 1

�a 1 b ln�t�tw��aXr �t, tw�, where a and b are determinedby the magnitude of the fluctuations, and Xr �t, tw� is arandom variable drawn from a time-independent PDF thatgoverns the fluctuations of the surfaces. In our approxi-mation, which describes uncorrelated drift between twosurfaces (i.e., a random walk), a � 1�2.

Figure 3 displays r���Cr�t, tw ���� for several choices of theratio t�tw. Interestingly enough, all the curves have anoticeable peak at a value of Cr that is independent oft and tw, with a height that decreases significantly with

0

0.5

1

1.5

2

2.5

-1 -0.5 0 0.5 1

ρ(C

r)

Cr

t/tw = 2 →

t/tw = 8 →

t/tw = 32 →

nw = 3nw = 7nw = 3nw = 5nw = 3

Response

FIG. 3. PDF of the local correlations for the ratios t�tw � 2, 8,and 32, averaged over 36 realizations of disorder. L � 32, T �0.72Tc , and V � 33. The waiting-times are tw � 104�2nw21�MCS, with nw given in the key. We include the PDF of thelocal integrated responses, xr �t, tw�, for nw � 3 and t�tw � 2.[r�xr� is divided by 4 to fit on the figure.]

237201-3

VOLUME 88, NUMBER 23 P H Y S I C A L R E V I E W L E T T E R S 10 JUNE 2002

0

0.4

0.8

1.2

-3 -2 -1 0 1

ρ(s-1

ln C

r/C

typ)

s-1 ln Cr/Ctyp

t/tw = 2, 4, 8, ..., 128

tw = 104, ..., 1.28 106

FIG. 4. Scaling of r���Cr �t, tw����. L � 32, T � 0.72Tc ,V � 133, and 28 pairs �t, tw�, with t � 104�2n21� MCS, n �1, . . . , 8, and tw � 104�2nw21� MCS, nw � 1, . . . , n 2 1.Ctyp � exp�lnC�, where �· · ·� denotes an average over thedistribution and s rescales the maximum value of each curve.With this choice of V we eliminate the Cr , 0; cf. Fig. 3.

increasing ratio t�tw . The form derived above for wr �t� 2

wr�tw� explains the approximate collapse of r���Cr �t, tw����for a fixed ratio t�tw, as shown in Fig. 3 for a small valueof V . (Because of mixing with the stationary part, the t�tw

scaling worsens when V increases.) Barely noticeable inFig. 3 is a slow drift of the curves for increasing valuesof tw at fixed ratio t�tw such that the height of the peakdecreases while the area below the tail at lower values ofCr increases. This trend leads to the “subaging” scalingobserved for bulk quantities [2,6].

Furthermore, the above expression for wr�t� 2 wr�tw�implies that the PDFs for all of the 28 pairs of times �t, tw�should collapse by rescaling with two parameters: lnCtypand s, corresponding, respectively, to the nonrandom partin wr�t� 2 wr�tw� and to the width of the random part(see Fig. 4). The scaling curve itself gives the PDF forXr �t, tw �. The rather good collapse of the curves shouldbe improved by further knowledge of CSP.

A local separation of time scales leads to a re-parametrization invariant action [11] with a soft modethat controls the aging dynamics. Our framework, basedon the analogy with Heisenberg magnets, predicts theexistence of a local relationship between Cr and xr ,expressed by r���Cr �t, tw �, xr �t, tw ���� being sharply concen-trated along the global x�C� in the �C, x� plane. Underadditional assumptions, we obtained the scaling behaviorof r���Cr �t, tw ���� for all �t, tw�. Our simulations bothconfirm these predictions and uncover striking regularitiesin the geometry of local fluctuations [17]. These resultsopen the way to a systematic study of local dynamicfluctuations in glassy systems and they suggest a numberof exciting avenues for future research. On the theoreticalside, this framework could be applied to a large variety ofglassy models, including those without explicit disorder.More interesting still are the experimental tests suggested

237201-4

by this work. For instance, the local correlations ofcolloidal glasses are accessible experimentally with theconfocal microscopy technique [8]. Similarly, cantilevermeasurements of noise spectra [9] allow probing of localfluctuations in the glassy phase of polymer melts. Theseare just two examples: any experiment that measures localfluctuations in glassy systems is a potential candidate fortesting our ideas.

We thank D. Huse and J. Kurchan for useful discus-sions. This work was supported in part by the NSF GrantNo. DMR-98-76208, the NSF-CNRS Grant No. 12931,and the Alfred P. Sloan Foundation. Supercomputing timewas allocated by the Boston University SCF.

[1] L. C. E. Struick, Physical Aging In Amorphous Polymersand Other Materials (Elsevier, New York, 1978).

[2] E. Vincent et al., in Proceedings of the Sitges Conference,edited by E. Rubi (Springer-Verlag, Berlin, 1997).

[3] L. F. Cugliandolo and J. Kurchan, Phys. Rev. Lett. 71, 173(1993); L. F. Cugliandolo and J. Kurchan, J. Phys. A 27,5749 (1994).

[4] J-P Bouchaud et al., in Spin Glasses and Random Fields,edited by A. P. Young (World Scientific, Singapore, 1998).

[5] J-L Barrat and W. Kob, Eur. Phys. J. B 13, 319 (2000).[6] M. Picco, F. Ricci-Tersenghi, and F. Ritort, Eur. Phys. J. B

21, 211 (2001).[7] M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000).[8] A. van Blaaderen and P. Wiltzius, Science 270, 1177

(1995); E. R. Weeks et al., Science 287, 627 (2000);W. K. Kegel and A. van Blaaderen, Science 287, 290(2000).

[9] E. Vidal-Russell and N. E. Israeloff, Nature (London) 408,695 (2000).

[10] P. H. Poole et al., Phys. Rev. Lett. 78, 3394 (1997);A. Barrat and R. Zecchina, Phys. Rev. E 59, R1299(1999); F. Ricci-Tersenghi and R. Zecchina, Phys. Rev. E62, R7567 (2000); C. Bennemann et al., Nature (London)399, 246 (1999); W. Kob et al., Phys. Rev. Lett. 79, 2827(1997).

[11] C. Chamon et al., cond-mat/0109150.[12] D. S. Fisher and D. A. Huse, Phys. Rev. Lett. 56, 1601

(1986).[13] S. Franz and H. Rieger, J. Stat. Phys. 79, 749 (1995);

E. Marinari et al., J. Phys. A 33, 2373 (2000); W. Koband J-L. Barrat, Eur. Phys. J. B 13, 319 (2000); A. Barratet al., Phys. Rev. Lett. 85, 5034 (2000); H. Makse andJ. Kurchan, Nature (London) 415, 614 (2002); J-L. Barratand L. Berthier, cond-mat/0110257; J-L. Barrat, Phys. Rev.E 57, 3629 (1998).

[14] A. Barrat and L. Berthier, Phys. Rev. Lett. 87, 087204(2001).

[15] T. S. Grigera and N. E. Israeloff, Phys. Rev. Lett. 83, 5038(2000); L. Bellon, S. Ciliberto, and C. Laroche, Europhys.Lett. 53, 511 (2001); D. Herisson and M. Ocio, cond-mat/0112378.

[16] L. F. Cugliandolo, D. S. Dean, and J. Kurchan, Phys. Rev.Lett. 79, 2168 (1997).

[17] H. E. Castillo, C. Chamon, L. F. Cugliandolo, and M. P.Kennett (to be published).

237201-4