Nonlinear response and fluctuation dissipation relations
Eugenio Lippiello, Federico Corberi, Alessandro Sarracino, Marco Zannetti
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p Nonlinear response and fluctuation dissipation relations
Eugenio Lippiello † Dipartimento di Scienze Fisiche, Universit´a di Napoli “Federico II”, 80125 Napoli, Italy.
Federico Corberi ‡ , Alessandro Sarracino ¶ and Marco Zannetti § Dipartimento di Matematica ed Informatica via Ponte don Melillo,Universit`a di Salerno, 84084 Fisciano (SA), Italy
A unified derivation of the off equilibrium fluctuation dissipation relations (FDR) is given forIsing and continous spins to arbitrary order, within the framework of Markovian stochastic dynamics.Knowledge of the FDR allows to develop zero field algorithms for the efficient numerical computationof the response functions. Two applications are presented. In the first one, the problem of probingfor the existence of a growing cooperative length scale is considered in those cases, like in glassysystems, where the linear FDR is of no use. The effectiveness of an appropriate second order FDRis illustrated in the test case of the Edwards-Anderson spin glass in one and two dimensions. Inthe second one, the important problem of the definition of an off equilibrium effective temperaturethrough the nonlinear FDR is considered. It is shown that, in the case of coarsening systems, theeffective temperature derived from the second order FDR is consistent with the one obtained fromthe linear FDR. † [email protected] ‡ [email protected] ¶ [email protected] § [email protected]: 05.70.Ln, 75.40.Gb, 05.40.-a I. INTRODUCTION
The statistical mechanics of systems out of equilibrium is a rapidly evolving subject, due to the intensive researchin the slow relaxation phenomena arising in several different contexts, such as coarsening systems, glassy and granularmaterials, colloidal systems etc. Understanding the basic mechanism underlying the slow relaxation is an issue ofmajor importance. In particular, a key question is whether the large time scales are due to cooperative effects onlarge length scales. For non disordered coarsening systems this is certainly the case, since the observed power lawrelaxation can be directly related to the growth of the time dependent correlation length, or domain size [1]. In thecase of disordered or glassy systems the establishment of such a connection is much more problematic, due to thedifficulty of pinpointing the observables fit to the task.The use of the nonlinear susceptibilities has been advocated [2, 3, 4] as experimental or numerical probe apt todetect the heterogeneous character of the glassy relaxation and, possibly, to uncover the existence of the growing lengthscale responsible of the slowing down of the relaxation. However, this requires to establish clearly the relationshipbetween non linear susceptibilities and correlation functions, in order to make sure what actually do the nonlinearsusceptibilities probe. In other words, the problem of the derivation of the nonlinear fluctuation dissipation relations(FDR) in the out of equilibrium regime needs to be addressed. As a matter of fact, this has been one of the most fruitfullines of investigation in the field of slow relaxation, although mostly limited to the domain of linear response [5, 6].In this paper we approach the problem on fairly general grounds. Within the framework of Markovian stochasticevolution, we bring to the fore the structural elements which are common to discrete and continous spins. We developthe formal apparatus necessary for a unified derivation of the FDR in the two cases and to arbitrary order. We alsoshow, expanding on the work of Semerjian et al. [7], how the nonlinear FDR of arbitrary order can be derived from afluctuation principle [8] also in the off equilibrium regime. This allows to regard the FDR as a manifestation of theconstraint imposed on the dynamics by the requirement of microscopic reversibility.The immediate application of the FDR is in the development of algorithms for the computation of the responsefunctions without the imposistion of an external perturbation, the so called zero field algorithms. The numericaladvantages of a zero field algorithm are remarkable. These have been illustrated and discussed in detail, in the linearcase, in a recent paper [9]. In the present paper we apply the zero field algorithm to the computation of the nonlinearresponse functions. We consider two cases, where knowledge of a nonlinear response function is required. The firstone arises in the search for a growing length scale in the context of glassy systems, as mentioned above. The presenceof quenched or self-induced disorder makes the linear response function short ranged, compelling to resort to thenonlinear ones. The second one deals with the extension of the effective temperature concept [10] to nonlinear order.This is an important and difficult problem. Here, we consider it in the context of non disordered coarsening systems,showing the consistency of the effective temperatures derived from the linear and the second order FDR.The paper is organised as follows: the formal developments concerning the time evolution, the response functions,the FDR and the fluctuation principle are presented in sections 2, 3, 4 and 5, respectively. In section 6 the problemof the detection of a cooperative length trough a nonlinear susceptibility is addressed, while the effective temperatureto nonlinear order is treated in section 7. Concluding remarks are presented in section 8.
II. FORMALISM AND TIME EVOLUTION
Let us consider a system in contact with a thermal reservoir, whose microscopic states are the configurations σ = [ σ i ] of the N degrees of freedom, discrete or continous, placed on the discrete set of sites i = 1 , .., N . Assuminga Markovian stochastic dynamics, the time evolution of the system is fully specified once the probability distribution P ( σ, t ) at some initial time t is given, together with the transition probability P ( σ, t | σ ′ , t ′ ) for any pair of times t ≤ t ′ ≤ t . Observables, or random variables, are functions A ( σ ) defined over the phase space of the microscopicstates, whose expectations are given by h A ( t ) i = X σ A ( σ ) P ( σ, t ) (1)where P ( σ, t ) = P σ ′ P ( σ, t | σ ′ , t ) P ( σ ′ , t ).In order to keep the notation compact, it is convenient to switch to the operator formalism [11], whereby themicroscopic states introduced above are the set of labels of the basis vectors | σ i = O | σ i i , i = 1 , ...N (2)of a vector space. The single site states obey the orthonormality relation h σ i | σ ′ j i = δ σ i ,σ ′ j δ i,j . With this notation, theprobability distributions P ( σ, t ) become the time dependent vectors | P ( t ) i such that P ( σ, t ) = h σ | P ( t ) i (3)and the transition probability is associated to the propagator of the process ˆ P ( t | t ′ ), which is the operator whosematrix elements are given by P ( σ, t | σ ′ , t ′ ) = h σ | ˆ P ( t | t ′ ) | σ ′ i . (4)For stochastically continous processes, through the first order expansionˆ P ( t + ∆ t | t ) = ˆ I + ˆ W ( t )∆ t + O (∆ t ) (5)where ˆ I is the identity operator, there remains defined the generator of the process ˆ W ( t ). Then, the pair | P ( t ) i , ˆ W ( t )contains all the information on the process, since the propagator can be written asˆ P ( t | t ′ ) = T (cid:18) e R tt ′ ds ˆ W ( s ) (cid:19) (6)where T is the time ordering operator. In differential form this is equivalent to the equations ∂∂t ˆ P ( t | t ′ ) = ˆ W ( t ) ˆ P ( t | t ′ ) (7) ∂∂t ′ ˆ P ( t | t ′ ) = − ˆ P ( t | t ′ ) ˆ W ( t ′ ) (8)from the first one of which follows the equation of motion of the state vector ∂∂t | P ( t ) i = ˆ W ( t ) | P ( t ) i . (9)Notice that the normalization of probabilities imposes on the generator the condition h−| ˆ W ( t ) = 0 (10)where h−| = P σ h σ | is called the flat vector. We shall assume that at each instant of time the generator satisfies thedetailed balance condition e β ˆ H ( t ) ˆ W ( t ) e − β ˆ H ( t ) = ˆ W † ( t ) (11)where ˆ H ( t ) denotes the Hamiltonian of the system, with a possible time dependence. This implies that the istantaneousGibbs state | P ( t ) i β = 1 Z ( t ) e − β ˆ H ( t ) |−i (12)is an invariant state in the sense that ˆ W ( t ) | P ( t ) i β = 0 . (13)The time dependent partition function is given by Z ( t ) = h−| e − β ˆ H ( t ) |−i , as it follows from the normalization of thestate h−| P ( t ) i β = 1. From now on we shall adopt the notation |·i β for the Gibbs states.Notice that Eq. (1) requires that observables correspond to diagonal operators. Each random function A ( σ ) ismapped into the operator ˆ A = A (ˆ σ ), where the operators ˆ σ i are defined byˆ σ i | σ i = σ i | σ i (14)in such a way that the expectation (1) can be written as h ˆ A ( t ) i = h−| ˆ A | P ( t ) i . (15)If ˆ W and | P i β are time independent, one can show that the multi-time expectations of observables in the stationarystate obey the Onsager relation [12]. Namely, if ˆ A , ˆ A , ..., ˆ A n is a set of observables and t n ≥ t n − ... ≥ t an orderedsequence of instants of time, then h ˆ A n ( t n ) ˆ A n − ( t n − ) ... ˆ A ( t ) ˆ A ( t ) i β = h ˆ A ( t n ) ˆ A ( t n − ( t − t )) ... ˆ A n − ( t + ( t n − t n − )) ˆ A n ( t ) i β (16)where h ˆ A n ( t n ) ˆ A n − ( t n − ) ... ˆ A ( t ) ˆ A ( t ) i β = h−| ˆ A n ˆ P ( t n | t n − ) ˆ A n − ... ˆ A ˆ P ( t | t ) ˆ A | P i β . (17)Finally, using Eqs. (7) and (8), it is straightforward to show that the time derivative of a multi-time expectation isgiven by ∂∂t k h ˆ A n ( t n ) .. ˆ A k ( t k ) .. ˆ A ( t ) i = h ˆ A n ( t n ) .. [ ˆ A k , ˆ W ( t k )]( t k ) .. ˆ A ( t ) i (18)namely, the time derivative in front of the expectation amounts to the insertion, at the time t k , of the commutator[ ˆ A k , ˆ W ( t k )] inside the expectation. This applies for all k = 1 , .., n . In particular, in the case of a single observable ∂∂t h ˆ A ( t ) i = h−| ˆ A ˆ W ( t ) | P ( t ) i = h−| [ ˆ A, ˆ W ( t )] | P ( t ) i (19)where the second equality is a consequence of Eq. (10). It should be noticed that, in general, [ ˆ A, ˆ W ( t )] is not anobservable. III. RESPONSE FUNCTIONS
Let us assume that the generator ˆ W ( t ) depends on time through an external field h i ( t ). Then, upon varying h i ( t ),there remains defined a family of stochastic processes and one may ask whether these stochastic processes are relatedone to the other. In particular, one would like to know whether the process with the generic h i ( t ) can be reconstructedfrom the properties of the unperturbed process, the one with the particular choice h i ( t ) ≡
0. In order to answer thequestion, let us start from the statement that all the information in the process, with a given h i ( t ), is contained inthe full hierarchy of the time dependent moments M i ,...,i n ( t , ..., t n , [ h i ( t ′ )]) = h ˆ σ i ( t ) ... ˆ σ i n ( t n ) i (20)each of which is a functional of h i ( t ). Assuming analiticity, one can write the formal expansion M i ..i n ( t , .., t n , [ h i ( t ′ )]) = M i ..i n ( t , .., t n ) + (21) ∞ X m =1 m ! X j ..j m Z tt w dt ′ ... Z tt w dt ′ m R ( n,m ) i ..i n ; j ..j m ( t , .., t n ; t ′ , .., t ′ m ) h j ( t ′ ) ...h j m ( t ′ m )where M i ..i n ( t , .., t n ) is the unperturbed moment, ( t w , t ) is the time interval over which the action of the externalfield is considered, and R ( n,m ) i ..i n ; j ..j m ( t , .., t n ; t ′ , .., t ′ m ) = δ m M i ..i n ( t , .., t n , [ h i ( t ′ )]) δh j ( t ′ ) ...δh j m ( t ′ m ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 (22)is the m -th order response function of the n -th moment. Therefore, the question asked above can be positivelyanswered if the response functions can be expressed in terms of quantities computable in the unperturbed process,namely if the FDR of arbitrary order can be obtained.Without loss of generality, we limit ourselves to work out the response functions for the first moment. For theresponses of the higher moments there is nothing conceptually different, just the formalism gets more involved. As amatter of fact, in section VI we shall deal with the second order response of the second moment. Coming back to thefirst moment, from M i ( t, [ h j ( t ′ )]) = h−| ˆ σ i ˆ P h ( t | t w ) | P ( t w ) i (23)where ˆ P h ( t | t w ) is the propagator in the presence of the field, follows R (1 ,m ) i ; j ..j m ( t, t , .., t m ) = h−| ˆ σ i δ m ˆ P h ( t | t w ) δh j ( t ) ...δh j m ( t m ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 | P ( t w ) i (24)which obviously vanishes if t j ( t w , t ) for any j = 1 , .., m . Let us write explicitely the first two derivatives δ ˆ P h ( t | t w ) δh j ( t ) = T " e R ttw ds ˆ W ( s ) ∂ ˆ W ( t ) ∂h j ( t ) = ˆ P h ( t | t ) ∂ ˆ W ( t ) ∂h j ( t ) ˆ P h ( t | t w ) (25)and δ ˆ P h ( t | t w ) δh j ( t ) δh j ( t ) = T ( e R ttw ds ˆ W ( s ) " ∂ ˆ W ( t ) ∂h j ( t ) ∂ ˆ W ( t ) ∂h j ( t ) + ∂ ˆ W ( t ) ∂h j ( t ) δ (12) = ˆ P h ( t | t M ) ∂ ˆ W ( t M ) ∂h j M ( t M ) ˆ P h ( t M | t m ) ∂ ˆ W ( t m ) ∂h j m ( t m ) ˆ P h ( t m | t w )+ ˆ P h ( t | t ) ∂ ˆ W ( t ) ∂h j ( t ) ˆ P h ( t | t w ) δ (12) (26)where t M = max( t j ), t m = min( t j ), j M , j m are the sites where the field acts at the times t M or t m , respectively, and δ ( np ) = δ j n ,j p δ ( t n − t p ). The third order derivative is written down in Appendix I.In order to go further, we must specify how ˆ W ( t ) depends on the field h i ( t ) and that is where the distinctionbetween discrete spins and continous spins comes in. In the following, in order to keep the derivation as simple aspossible we shall specialize to the case of Ising spins with single spin-flip dynamics. The generalization to q -state spins(e.g. clock models) and to dinamical rules with conservation laws, such as Kawasaki spin-exchange, turns out to bestraightforward. A. Ising spins
For Ising spins the state vectors | σ i = ± i are represented by the column vectors | σ i = 1 i = (cid:18) (cid:19) , | σ i = − i = (cid:18) (cid:19) . A generator of single spin flip dynamics is of the typeˆ W = 1 N N X i =1 ˆ W i (27)where ˆ W i has non vanishing matrix elements only between states which differ for the value of σ i . The form of ˆ W i isobtained imposing the detailed balance condition (11). Singling out the site i , the Hamiltonian can be written in theform ˆ H ( t ) = H ([ˆ σ z ] i , t ) + { h W ([ˆ σ z ] i ) − h i ( t ) } ˆ σ zi , (28)where ˆ σ zi is the z Pauli matrix ˆ σ zi = (cid:18) − (cid:19) , [ˆ σ z ] i stands for the set of all spins except for the one on the i -th site and ˆ h Wi = h W ([ˆ σ z ] i ) is the Weiss field on thesite i . Inserting into the detailed balance condition (11), one finds the generator of the Glauber type [13]ˆ W i ( t ) = (ˆ σ xi − ˆ I ) e β [ˆ h Wi − h i ( t )]ˆ σ zi (29)where ˆ σ xi is the x Pauli matrix ˆ σ xi = (cid:18) (cid:19) . Hence, the derivatives are given by ∂ n ˆ W ( t ) ∂h nj ( t ) = ( − β ) n ˆ W j ( t )(ˆ σ zj ) n (30)with (ˆ σ zj ) n = (cid:26) ˆ I, for n evenˆ σ zj , for n odd . (31)For later reference, let us write here the identityˆ W i ˆ σ zi = 12 [ ˆ W i , ˆ σ zi ] + 12 { ˆ W i , ˆ σ zi } (32)obtained by adding the commutator and the anticommutator. It is straightforward to verify that the anticommutatoris a diagonal operator and, therefore, is an observable which we will denote byˆ B i ( t ) = { ˆ σ zi , ˆ W ( t ) } . (33)Due to Eq. (10), the average of the left hand side of Eq. (32) vanishes and Eq. (19) for ˆ σ zi can be rewritten as ∂∂t h ˆ σ zi ( t ) i = h ˆ B i ( t ) i . (34) B. Continuous spins
In order to avoid confusion, here we shall denote by ϕ = [ ϕ i ] the set of the N continuous variables. Let us assumean equation of motion of the Langevin type ∂∂t ϕ i ( t ) = B i ( t ) + η i ( t ) (35)where the drift B i is related to the Hamiltonian, or free energy functional H [ ϕ ], by B i ( t ) = − ∂ H [ ϕ ( t )] ∂ϕ i (36) η i ( t ) is a gaussian white noise with expectations h η i ( t ) i = 0 , h η i ( t ) η j ( t ′ ) i = 2 T δ i,j δ ( t − t ′ ) (37)and T is the temperature of the thermal reservoir. The generator of the corresponding Markov process is the Fokker-Planck operator ˆ W F P = X i ˆ W F Pi (38)where ˆ W F Pi = −{ T ˆ p i + iB i [ ˆ ϕ ]ˆ p i + D i [ ˆ ϕ ] } , (39) B i [ ˆ ϕ ] is defined through Eq. (36) and D i [ ˆ ϕ ] = − ∂ H [ ϕ ] ∂ϕ i (cid:12)(cid:12)(cid:12)(cid:12) ϕ = ˆ ϕ . (40)The conjugated operators ˆ ϕ i and ˆ p i , defined by h ϕ ′ | ˆ ϕ i | ϕ i = ϕ i h ϕ ′ | ϕ i (41) h ϕ ′ | ˆ p i | ϕ i = − i ∂∂ϕ i h ϕ ′ | ϕ i (42)obey the commutation relation [ ˆ ϕ i , ˆ p j ] = iδ i,j (43)and satisfy the equalities h−| ˆ p i = 0 (44) h−| ˆ W F Pi = 0 . (45)The external field enters the free energy functional linearly H h [ ϕ, t ] = H [ ϕ ] − X i h i ( t ) ϕ i (46)yielding B h,i [ ϕ, t ] = B i [ ϕ ] + h i ( t ) , (47)while D i [ ϕ ] remains unaltered. Therefore, the generator is changed intoˆ W F Ph,i = ˆ W F Pi − ih i ( t )ˆ p i (48)and the derivatives are given by ∂ n ˆ W F Ph,i ∂h ni ( t ) = (cid:26) − i ˆ p i , for n = 10 , for n > . (49)Obviously, the identity (32) holds also in this caseˆ W F Pi ˆ ϕ i = 12 [ ˆ W F Pi , ˆ ϕ i ] + 12 { ˆ W F Pi , ˆ ϕ i } (50)and, using the definiton (39) together with the commutation relation (43), it is not difficul to rewrite it in the moreconvenient form i ˆ p i = β n [ ˆ W F Pi , ˆ ϕ i ] + B i [ ˆ ϕ ] o . (51)Again, since the average of the left hand side vanishes, we get the analogous of Eq. (34) ∂∂t h ˆ ϕ i ( t ) i = h B i [ ˆ ϕ ]( t ) i (52)showing that the anticommutator (33) for Ising spins and the drift (36) for continous spins play the same role in theevolution, thus justifying the choice of the same notation. IV. FLUCTUATION DISSIPATION RELATIONS
Let us now return to the general treatment, valid for discrete and continous spins alike, with the operator ˆ σ standingeither for ˆ σ z or for ˆ ϕ , depending on the context. Which is the case, will be specified whenever necessary.Comparing the left hand sides of Eqs. (32)and (51) with the first derivatives (30) and (49), we can write the basicequation in this paper ∂ ˆ W∂h i = − β (cid:16) [ ˆ W i , ˆ σ i )] + B i [ˆ σ ] (cid:17) (53)where, in the case of continous spins, the superscript FP on the Fokker-Planck operator has been dropped and, asspecified above, ˆ B i stands for ˆ B i = (cid:26) { ˆ σ i , ˆ W i } , for Ising spins B i [ ˆ ϕ ] , for continous spins . (54)We are now in the position to derive the FDR, by going through the following steps:1. according to Eq. (24), the response function is obtained by appropriate insertions of derivatives of the generatorin between propagators, as exemplified in Eqs. (25) and (26)2. according to Eq. (53), each first derivative of ˆ W can be repalced by the sum of the commutar and the driftoperator ˆ B i
3. according to Eq. (18), the insertion of the commutator amounts to a time derivative acting in front of theaverage.The above steps exaust all there is to do in the continous spin case, since there are no derivatives of the generatorof order higher than the first. Higher derivatives do, instead, appear in the Ising case producing singular terms. Inorder to see how the procedure works in practice, let us carry out the computation up to second order. The explicitresult for the third order response function is presented in Appendix I.
A. Linear FDR
Let us begin with the simplest case of the linear response function. From Eqs. (24) and (25) follows R (1 , i ; j ( t, t ) = h−| ˆ σ i ˆ P ( t | t ) ∂ ˆ W ( t ) ∂h j ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 ˆ P ( t | t w ) | P ( t w ) i (55)and, using Eqs. (53) and (18), this becomes the linear FDR R (1 , i ; j ( t, t ) = β (cid:20) ∂∂t M ij ( t, t ) − h ˆ σ i ( t ) ˆ B j ( t ) i (cid:21) (56)due to the appearence of unperturbed correlation functions of observables in the right hand side. This result has beenknown for some time for continous spins, see e.g. Ref. [6], while for Ising spins has been derived for the first time inRef. [9], where it has been exploited to develop the zero field algorithm for the computation of R (1 , i ; j ( t, t ), mentionedin the Introduction.If the system is in equilibrium, the averages in the right hand side are equilibrium averages and, invoking theOnsager relation (16), one gets h ˆ σ zi ( t ) ˆ B j ( t ) i β = h−| [ˆ σ j , ˆ W ] ˆ P ( t | t )ˆ σ i | P i β = − ∂∂t M ij ( t, t ) (57)after using space and time translation invariance in the last equality, Eqs. (32) or (51) to eliminate ˆ B j , together withthe normalization conditions (10) or (44). Inserting this into Eq. (56), the equilibrium fluctuation dissipation theorem(FDT) is recovered R (1 , ij ( t, t ) = β ∂∂t C ij ( t, t ) (58)where we have introduced the pair correlation function C ij ( t, t ) = M ij ( t, t ) − M i ( t ) M j ( t ) (59)since it is clear that in equilibrium the one time quantities are time independent and the time derivative of C ij ( t, t )coincides with that of M ij ( t, t ). B. Second order FDR
From Eqs. (24) and (26), the second order, or two kicks nonlinear response function, is given by R (1 , i ; j j ( t, t , t ) = h−| ˆ σ i ˆ P ( t | t M ) ∂ ˆ W ( t M ) ∂h j M ( t M ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 ˆ P ( t M | t m ) ∂ ˆ W ( t m ) ∂h j m ( t m ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 ˆ P h ( t m | t w ) | P ( t w ) i + h−| ˆ σ i ˆ P ( t | t ) ∂ ˆ W ( t ) ∂h j ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 ˆ P ( t | t w ) | P ( t w ) i δ (12) . (60)Substituting, next, Eq. (53) for the first derivatives and the identityˆ W i = 12 ˆ σ zi ˆ B i + 12 ˆ σ zi [ˆ σ zi , ˆ W i ] (61)for the manipulation of the singular term, eventually one finds the second order FDR R (1 , i ; j j ( t, t , t ) = ( β/ n ∂∂t M ∂∂t m M ij M j m ( t, t M , t m ) − ∂∂t M h ˆ σ i ( t )ˆ σ j M ( t M ) ˆ B j m ( t m ) i − ∂∂t m h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j m ( t m ) i + h ˆ σ i ( t ) ˆ B j M ( t M ) ˆ B j m ( t m ) i o + ( β / n h ˆ σ zi ( t )ˆ σ zj m ( t M ) ˆ B j m ( t m ) i + ∂∂t m M ij m j m ( t, t M , t m ) o δ j M ,j m δ ( t M − t m ) (62)where, it is worth recalling, the last singular contribution in the braces is present only for Ising spins.Again, at stationarity the Onsager relation can be used to eliminate the ˆ B ’s entering with the shortest time in favorof time derivatives. In so doing, the second and the fourth contribution in the right hand side become identical to thefirst and the third one, respectively. Furthermore, in the stationary state the last contribution in the right hand sidevanishes, eventually yielding what we may call the second order FDT R (1 , i ; j j ( t, t , t ) = ( β / (cid:26) ∂∂t M ∂∂t m C ij M j m ( t, t M , t m ) − ∂∂t m h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j m ( t m ) i β (cid:27) (63)where we have substituted the time derivative of the moment with that of the correlation function. Notice that thereremains an ineliminable presence of ˆ B in the second term in the right hand side, which carries the information on thespecific rule governing the time evolution. This is a distinctive feature of all FDR of order higher than linear, makingthem less universal than the linear one, which is ˆ B dependent only off equilibrium. V. FLUCTUATION PRINCIPLE
In the following we show that the FDR obtained in the previous section arise as a consequence of a fluctuationprinciple [8]. For convenience, the derivation is presented for the continous spin case, but it can be worked outalong the same lines also for Ising spins. The idea of a derivation from the fluctuation principle was implemented bySemerjian et al. [7] in the stationary case. Here we extend the derivation to the more general off equilibrium case.Let us call experimental protocol the assigned time dependence of the external field h ( t ) in some time interval( t , t F ). Then, from the detailed balance condition follows that the probability P β ([ ϕ ( t )] | ϕ , [ h ( t )]) of a path [ ϕ ( t )],taking place under the protocol [ h ( t )] and conditioned to the initial value ϕ ( t ) = ϕ is related to the probability ofthe reverse path e ϕ ( t ) = ϕ ( e t ) , e t = t F − t + t (64)under the reverse protocol [ e h ( t )] and conditioned to e ϕ ( t ) = ϕ F by P β ([ ϕ ( t )] | ϕ , [ h ( t )]) exp (cid:26) − β Z t F t dt h ( t ) ˙ ϕ ( t ) (cid:27) = P β ([ e ϕ ( t )] | ϕ F , [ e h ( t )]) exp { β [ H ( ϕ ) − H ( ϕ F )] } (65)where the subscript β is there to remind that the evolution takes place while the system is in contact with a singlethermal reservoir at the inverse temperature β . Multiplying both sides by ϕ i,F P I ( ϕ ), where P I ( ϕ ) is an arbitraryprobability distribution, and summing over the set C ( t , t F ) of all paths in the interval ( t , t F ) one finds Z C ( t ,t F ) d [ ϕ ( t )] ϕ i,F P β ([ ϕ ( t )] | ϕ , [ h ( t )]) exp (cid:26) − β Z t F t dt h ( t ) ˙ ϕ ( t ) (cid:27) P I ( ϕ )= Z β, Z C ( t ,t F ) d [ ϕ ( t )] P I ( ϕ ) e β H ( ϕ ) P β ([ e ϕ ( t )] | ϕ F , [ e h ( t )]) ϕ i,F P β, ( ϕ F ) (66)0where P β, and Z β, denote the equilibrium distribution and the corresponding partition function, in absence of theexternal field. Hence, the above result can be rewritten more compactly as (cid:28) ϕ i ( t F ) exp (cid:26) − β Z t F t dt h ( t ) ˙ ϕ ( t ) (cid:27)(cid:29) I → β, [ h ( t )] = Z β, D P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) E β, → β, [ e h ( t )] (67)where h·i I → β, [ h ( t )] stands for the average in the process starting with the initial condition P I , thereafter in contact withthe thermal resorvoir β and evolving with the protocol [ h ( t )], while h·i β, → β, [ e h ( t )] stands for the process in contact withthe thermal resorvoir β , starting with the unperturbed equilibrium distribution P β, and evolving with the reverseprotocol [ e h ( t )].The next step is to expand both sides in powers of [ h ( t )] about h ( t ) ≡ Z β, D P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) E β, = h ϕ i ( t F ) i I → β, (68)where h·i I → β, stands for the average in the off-equilibrium process starting with P I and evolving in contact with thethermal reservoir in absence of the external field. At higher orders one gets Z β, δ n (cid:10) P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) ...δ e h j n ( e t n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 = n X p =0 ( − β ) n − p X P ( j ,..,j n − p | j n − p +1 ,..,j n ) δ p h ϕ i ( t F ) ˙ ϕ j ( t ) ... ˙ ϕ j n − p ( t n − p ) i I → β, [ h ( t )] δh j n − p +1 ( t n − p +1 ) ...δh j n ( t n ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 (69)where the second summation is over all the distinct permutations P ( j , .., j n − p | j n − p +1 , .., j n ) between the two sets of indeces ( j , ..., j n − p ) and ( j n − p +1 , ..., j n ). For instance, at the firsttwo orders the above formula reads Z β, δ (cid:10) P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 = − β ∂∂t h ϕ i ( t F ) ϕ j ( t ) i I → β, + δ h ϕ i ( t F ) i I → β, [ h ( t )] δh j ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 (70)and Z β, δ (cid:10) P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) δ e h j ( e t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 = β ∂∂t ∂∂t h ϕ i ( t F ) ϕ j ( t ) ϕ j ( t ) i I → β, − β ∂∂t δ h ϕ i ( t F ) ϕ j ( t ) i I → β, [ h ( t )] δh j ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 − β ∂∂t δ h ϕ i ( t F ) ϕ j ( t ) i I → β, [ h ( t )] δh j ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 + δ h ϕ i ( t F ) i I → β, [ h ( t )] δh j ( t ) δh j ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 . (71)From these two equations the off equilibrium FDR (56) and (62) can be recovered. Here we show this explicitely inthe linear case, the extension to higher orders being straightforward.Recalling the definition (22) of the response functions and making the replacements t F → t and t → t w , Eq. (70)can be rewritten as R (1 , i ; j ( t, t ) = β ∂∂t h ϕ i ( t ) ϕ j ( t ) i I → β, + Z β, δ (cid:10) P I ( ϕ ( t )) e β H ( ϕ ( t )) ϕ i ( t w ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 . (72)1According to Eq. (49), the derivative with respect to e h j ( e t ) can be replaced by the insertion of − i ˆ p j ( e t ) and, usingEq. (51), the second term in the right hand side can be rewritten as Z β, δ (cid:10) P I ( ϕ ( t )) e β H ( ϕ ( t )) ϕ i ( t w ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 = Z β, β ∂∂ e t D P I ( ϕ ( t )) e β H ( ϕ ( t )) ϕ j ( e t ) ϕ i ( t w ) E β, − Z β, β D P I ( ϕ ( t )) e β H ( ϕ ( t )) B j ( e t ) ϕ i ( t w ) E β, . (73)Furthermore, since after setting to zero the external field the averages become equilibrium averages, using the Onsagerrelation we get Z β, δ (cid:10) P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 = − Z β, β ∂∂t D ϕ i ( t ) ϕ j ( t ) P I ( ϕ ( t w )) e β H ( ϕ ( t w )) E β, − Z β, β D ϕ i ( t ) B j ( t ) P I ( ϕ ( t w )) e β H ( ϕ ( t w )) E β, . (74)The next step consists in the recognition that for an arbitrary observable AZ β, D A ( t ) P I ( ϕ ( t w )) e β H ( ϕ ( t w )) E β, = h A ( t ) i I → β, (75)since the factor Z β, P I ( ϕ ( t w )) e β H ( ϕ ( t w )) in the left hand side has the effect of undoing the equilibrium initial conditionand replacing it with P I . Hence, Eq. (74) can be rewritten as Z β, δ (cid:10) P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) (cid:11) β, → β, [ e h ( t )] δ e h j ( e t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)e h =0 = − β ∂∂t h ϕ i ( t ) ϕ j ( t ) i I → β, − β h ϕ i ( t ) B j ( t ) i I → β, (76)and inserting it into Eq. (72), eventually one finds R (1 , i ; j ( t, t ) = β ∂∂t h ϕ i ( t ) ϕ j ( t ) i I → β, − β h ϕ i ( t ) B j ( t ) i I → β, (77)thus recovering Eq. (56). VI. NONLINEAR SUSCEPTIBILITY AND GROWING LENGTH SCALE
In this section we consider the use of the response functions in the diagnostics of cooperative effects taking placeover large length scales during the relaxation, referring to the case of Ising spin systems. A partial and preliminaryaccount of the material in this section has been presented in Ref. [14].In non disordered coarsening systems, such as a ferromagnet quenched to or to below the critical point, the existenceof a growing dynamical correlation length is well captured through the scaling properties of the equal time correlationfunction C ij ( t ) = M ij ( t, t ) [1], recalling that in this kind of processes M i ( t ) ≡
0. In glassy systems, instead, quenchedor self-induced disorder renders the two body correlation function short ranged, making it necessary to resort to higherorder correlation functions. Attention has been particularly focussed on the four-point correlation function [15] C (4) ij ( t, t w ) = M iijj ( t, t w , t, t w ) − M ii ( t, t w ) M jj ( t, t w ) (78)2which describes the so called heterogeneities [16], namely the space fluctuations associated to the local time decor-relation. Although quite convenient in the numerical simulations, C (4) ij ( t, t w ) has the shortcoming of being hardlyaccessible in the experiments. Conversely, susceptibilities are more easily measurable and this has prompted theinvestigation of the nonlinear FDR.Here, we expand on the proposals [2, 3, 4, 14] to investigate dynamic scaling through higher order susceptibilities.Let us begin by considering the second order response of the second moment at equal times R (2 , ij ; j j ( t, t ; t , t ) = δ M ij ( t, t, [ h i ( t ′ ]) δh j ( t ) δh j ( t ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 . (79)Proceeding exactly as in the derivation of Eq. (62), the corresponding FDR is given by R (2 , ij ; j j ( t, t ; t , t ) = ( β/ n ∂∂t M ∂∂t m M ijj M j m ( t, t, t M , t m ) − ∂∂t M h ˆ σ zi ( t )ˆ σ zj ( t )ˆ σ zj M ( t M ) ˆ B j m ( t m ) i− ∂∂t m h ˆ σ zi ( t )ˆ σ zj ( t ) ˆ B j M ( t M )ˆ σ zj m ( t m ) i + h ˆ σ zi ( t )ˆ σ zj ( t ) ˆ B j M ( t M ) ˆ B j m ( t m ) i o + ( β / n h ˆ σ zi ( t )ˆ σ zj ( t )ˆ σ zj m ( t M ) ˆ B j m ( t m ) i + ∂∂t m M ijj m j m ( t, t, t M , t m ) o δ j M ,j m δ ( t M − t m ) . (80)Looking at the above equation, it may seem farfetched to relate the properties of C (4) to those of R (2 , , since the fourthorder moment appears explicitly only in the first term in the right hand side and in the singular term. Nonetheless,information on the existence of a growing correlation length can be obtained through a scaling argument, as it willbe shown below. Let us consider the more manageable integrated response function − χ (2 , ij ( t, t w ) = Z tt w dt Z tt w dt R (2 , ij ; ij ( t, t ; t , t ) − χ (1 , i ( t, t w ) χ (1 , j ( t, t w ) (81)where the χ (1 , i in the subtraction are the time integrals of the linear response function (56), i.e. χ (1 , i ( t, t w ) = Z tt w dt R (1 , i ; i ( t, t ) . (82)The reason for this subtraction and for the minus sign on the left hand side will be clear shortly. Assuming thateventually an equilibrium state is reached, from equilibrium statistical mechanics follows T lim t →∞ χ (2 , ij ( t, t w ) = C ij,eq (83)where, for simplicity, we have taken h ˆ σ zi i eq = 0. The scaling relation C ij,eq = ξ − d − η F C,eq ( | i − j | /ξ ), where ξ ( T ) isthe equilibrium correlation length, suggests a finite time scaling behavior of the form T χ (2 , ij ( t, t w ) = ξ x F (cid:18) | i − j | ξ , t /z ξ , t w t (cid:19) (84)where x = 4 − d − η and z is the dynamical exponent.Another quantity, which has been recently [3] considered in relation to C (4) , is the third order integrated responseof the first moment χ (1 , ij ( t, t w ) = − Z tt w dt dt dt R (1 , i ; ijj ( t ; t , t , t ) (85)where R (1 , i ; ijj is given by Eq.(111) of Appendix I. Again, from the large time result T lim t →∞ χ (1 , ij ( t, t w ) = C ij,eq (86)3we may infer the scaling behavior T χ (1 , ij ( t, t w ) = ξ x G (cid:18) | i − j | ξ , t /z ξ , t w t (cid:19) . (87)Notice that the prefactor in the definition (85) has been arranged in such a way that the large time limits of T χ (2 , and T χ (1 , are the same.These patterns of behavior have been checked numerically in the one-dimensional Ising model. The simulationhas been carried out through standard Montecarlo techniques with Glauber transition rates, where ˆ B i = ˆ σ zi − tanh( β P
0) and χ (1 , ij ( t, h → t w = 0 in order to reduce the variablesfrom three to two in the right hand sides. Then, recalling that for the 1 d Ising model η = 1 implies x = 0 and varyingthe temperature and the distance in such a way to keep | i − j | /ξ fixed, it is matter of showing that T χ (2 , ij ( t,
0) and T χ (1 , ij ( t,
0) are functions only of t /z /ξ . This is shown, with good accuracy, in Fig. (1) for the quenches to the threefinal temperatures T = 0 . , T = 0 . T = 1 . | i − j | /ξ . Both T χ (2 , ij ( t, T χ (1 , ij ( t,
0) grow from zero to the same asymptotic value C ij,eq = exp {− | i − j | /ξ } on the same timescale, asexpected. Note that the data for χ (1 , are much more noisy than those for χ (2 , . Since the same numerical resourceshave been allocated in the computation of each of these two quantities, we conclude that investigations based on χ (2 , are more efficient, at least numerically.We now discuss how χ (2 , can be effectively used for the measurement of a cooperative length in disorderedsystems. In order to improve the statistics further, it is convenient to consider the ~k = 0 component of the spaceFourier transform χ (2 , ~k =0 ( t, t w ) = (1 /N ) X i,j χ (2 , ij ( t, t w ) (88)which, using Eq. (84), scales as χ (2 , ~k =0 ( t, t w ) = ξ − d − η F (cid:18) t /z ξ , t w t (cid:19) . (89)In Ref. [14] we have computed numerically this quantity, with t w = 0, in the Edwards-Anderson (EA) model withHamiltonian ˆ H = P ij J ij ˆ σ zi ˆ σ zj and d = 1 ,
2. The data have been analysed as follows: for large t the curves saturateto the equilibrium value χ (2 , ~k =0 ,eq ∼ ξ − d − η . Using the known values of η , one can extract ξ . In the off-equilibriumregime t /z ≪ ξ the growing correlation length L ( t ) ∼ t /z is expected not to depend on ξ . Enforcing this conditionfrom Eq. (89) one must have F ( t /z /ξ, ∼ ( ξ/L ( t )) η + d − , which yields χ (2 , ~k =0 ( t, ∼ L ( t ) − d − η . (90)This allows to determine L ( t ). After doing this, we checked for the data collapse by plotting ξ − d +2 η χ (2 , k =0 ( t,
0) vs L ( t ) /ξ for all the temperatures considered (see Figs. 2,3).Let us consider first the d = 1 EA model with bimodal distribution of the coupling constants J ij = ±
1. This systemis considered in order to test the method, since it can be mapped onto the ferromagnetic system, just considered above.Moreover, in this simple case, in addition to χ (2 , k =0 ( t, C k =0 ( t ), obtaining independent determinations of L ( t ) and ξ to compare with. This shows that the sets of data for L ( t )and ξ , obtained in both ways, are in agreement with each other and with the analytical behaviors up to the numericaluncertainty. The data collapse of χ (2 , k =0 ( t,
0) and of C k =0 ( t ) are shown in Fig. 2. Here, one clearly observes the offequilibrium regime, characterized by the power-law behavior of χ (2 , k =0 ( t,
0) and C k =0 ( t ) with exponents 4 − d − η and2 − η , respectivly. In the large time regime equilibration takes place with the convergence of χ (2 , k =0 ( t,
0) and C k =0 ( t )to χ (2 , k =0 ,eq and to C k =0 ,eq ( t ).4 FIG. 1: (Color online). Upper panel: The quantity T χ (2 , ij ( t,
0) is plotted against t /z /ξ ( z = 2) for six fixed values of r = | i − j | /ξ ( r = 0 , . , . , . , . , . T = 0 . , T = 0 . , T =1 . ξ ( T ) /ξ ( T ) = 2 and ξ ( T ) /ξ ( T ) = 4. Lower panel: Same for T χ (1 , ij ( t, t ξ -1 -2 -1 ξ η + d - χ k = ( t , ) T=1.2T=1T=0.8T=0.6T=0.4T=0.3 -2 -1 t ξ -1 -2 -1 ξ η - C k = ( t ) x η x η ( , ) FIG. 2: (Color online). Data collapse of χ (2 , k =0 ( C k =0 in the inset) for several temperatures in the d = 1 EA model ( η = 1 , z = 2).The dashed lines are the expected power-laws in the non equilibrium regime. After this explicit verification, we have turned to the d = 2 case, where the independent information on the structurefactor is not available. Both with bimodal and Gaussian distributions of J ij , the behavior of ξ extracted from χ (2 , ~k =0 ,eq ,using η = 0 [17], has been found consistent with previous results [17, 18]. The non-equilibrium behavior is compatiblewith a power law L ( t ) ∼ t /z ( T ) with a temperature dependent exponent in agreement with z ( T ) ≃ /T , as reportedin Ref. [19]. The data collapse of χ (2 , k =0 ( t,
0) is shown in Fig. 3. Notice also the additional collapse of the curves withbimodal and Gaussian bond distribution, further suggesting that the two models may share the same universalityclass at finite temperatures [17].
VII. EFFECTIVE TEMPERATURE
One of the most interesting developments in the study of the linear FDR out of equilibrium has been the introductionof the concept of effective temperature T eff [10]. The idea is that the off equilibrium behavior observed during slowrelaxation can be accounted for by the separation of the time scales for different subsets of degrees of freedom. Eachof these is regarded as in equilibrium with a different virtual thermostat at some appropriate effective temperature,which depends on the time scale and is different from the physical temperature of the real thermostat driving therelaxation. The value of T eff can be inferred by forcing the off equilibrium linear FDR in the form of the equilibriumFDT. Although appealing, this idea has turned out not to be applicable tout court , since T eff might turn out to beobservable dependent [20]. Nonetheless, with the proper caveats, the concept remains quite useful and suggestive. Inthis section we make a preliminary exploration of another open end in the important question of how general T eff canbe, investigating whether it is possible to extend to the nonlinear FDR the effective temperature concept, consistentlywith what it is done in the linear case.Let us first recall how T eff is defined from the linear FDR. For definiteness, the Ising spin case will be considered.Writing explicitely the time integral in Eq.(82), one has χ (1 , i ( t, t w ) = β Z tt w dt (cid:20) ∂∂t M ii ( t, t ) − h σ zi ( t ) ˆ B i ( t ) i (cid:21) . (91)Assuming M i ( t ) ≡ M ii ( t, t ) with the autocorrelation function6 t ξ −1 -2 -1 ξ η + d - χ k = ( t , ) T=1.8T=1.7T=1.6T=1.5T=1.4T=1.3T=1.2T=1.1T=1T=0.9T=0.8T=0.7 x η ( , ) FIG. 3: (Color online). Data collapse of χ (2 , k =0 for several T in the d = 2 EA model with bimodal (open symbols) or Gaussian(filled symbols) bond distribution, with z ( T ) = 4 /T and η = 0. The dashed line is the expected power-law in the non equilibriumregime. C ( t, t ), the quantity ψ (1) ( t, t w ) = Z tt w dt ∂∂t C ( t, t ) = 1 − C ( t, t w ) (92)for fixed t w is a monotonously increasing function of time, which allows to reparametrize t in terms of ψ (1) and towrite χ (1 , i ( t, t w ) in the form χ (1 , i ( t, t w ) = χ (1 , i ( ψ (1) , t w ) . (93)In equilibrium, where time translation invariance holds, the dependence on t w disappears and the parametric repre-sentation becomes linear χ (1 , i ( ψ (1) ) = βψ (1) (94)with the obvious consequence β = dχ (1 , i ( ψ (1) ) dψ (1) . (95)Off equilibrium, the parametric representation won’t be linear and an effective temperature can be defined by thegeneralization of the above relation β eff ( ψ (1) , t w ) = ∂χ (1 , i ( ψ (1) , t w ) ∂ψ (1) (96)with β eff = 1 /T eff .In order to see how T eff = T arises in a simple context, let us consider the relaxation to a low temperature phasecharacterized by ergodicity breaking and, therefore, by a non vanishing Edwards-Anderson order parameter q EA .In particular, let us think of the already mentioned coarsening process, like in a ferromagnet quenched to below thecritical point and relaxing via domain growth. In that case q EA coincides with the spontaneous magnetization squared M eq .7As t w → ∞ , the separation of time scales takes place. The short, or quasiequilibrium, time regime holds for C > M eq , that is ψ (1) < − M eq , while the large time scale sets in when C < M eq , or ψ (1) > − M eq . The existenceof this latter regime makes it evident the failure of equilibration, even in the t w → ∞ limit, since the autocorelationfunction falls below the Edwards-Anderson plateau. The behavior of ψ (1) , obtained from numerical simulations ofthe Ising model in d = 2, is shown in the inset of Fig. 4. Starting from zero, there is a fast growth in the short timeregime, followed by a plateau, more evident for large t w , and finally there is convergence, with the power law behavior1 − ψ (1) ( t, t w ) ∼ t − λ/z , toward the asymptotic value ψ (1) = 1. Notice that λ is the Fisher-Huse exponent [21] andthat the plateau flattens over the asymptotic value 1 − q EA as t w → ∞ .Correspondingly, the integrated response function can be written as the sum of two pieces [22] χ (1 , i ( ψ (1) , t w ) = χ (1 , st ( ψ (1) ) + χ (1 , ag ( ψ (1) , t w ) (97)where χ (1 , st is the stationary contribution arising from the equilibrated bulk of domains, while χ (1 , ag is the agingcontribution due to the off equilibrium domain walls. The stationary contribution obeys Eq. (94) in the short timeregime, saturates to its equilibrium value and remains constant in the large time regime, while the aging contributionvanishes as t w → ∞ according to χ (1 , ag ( ψ (1) , t w ) = t − aw F ( ψ (1) ) (98)where a > t w →∞ χ (1 , i ( ψ (1) , t w ) = (cid:26) βψ (1) , for 0 ≤ ψ (1) ≤ − M eq β (1 − M eq ) , for 1 − M eq < ψ (1) ≤ t w →∞ β eff ( ψ (1) , t w ) = (cid:26) β, for 0 ≤ ψ (1) ≤ − M eq , for 1 − M eq < ψ (1) ≤ . (100)Namely, the effective temperature coincides with the physical temperature in the short time regime, where the sytemappears equilibrated, while it is drastically different from it in the off equilibrium large time regime.Let us now carry out the parallel analysis on the ~k = 0 component of the second order integrated response function χ (2 , ~k =0 ( t, t w ) = 1 N X i,j Z tt w dt Z tt w dt R (2 , ij ; ij ( t, t ; t , t ) . (101)Notice that, although we use the same notation, this quantity differs from the one in Eq. (88), because of the overallsign and of the absence of the subtraction. Let us introduce the quantity ψ (2) ( t, t w ) = 12 1 N X i,j Z tt w dt Z tt w dt h ∂ ∂t M ∂t m M ijij ( t, t, t M , t m )+ h σ zi ( t ) σ zj ( t ) ˆ B i ( t M ) ˆ B j ( t m ) i i (102)with properties similar to those of ψ (1) ( t, t w ), as shown in Fig. 4. The main feature is the monotonous increasefrom zero to an asymptotic value well above the limiting value (1 − q EA ) that one would get from the equilibriumcalculation. This is the value at which a plateau develops as t w gets large, signaling the separation of time scales.Using ψ (2) to reparametrize the time t , from Eq. (63) follows that at equilibrium the FDR becomes linear χ (2 , ~k =0 ( ψ (2) ) = β ψ (2) . (103)Hence, by following the same reasoning as in the linear case, in the off equilibrium regime an effective temperaturecan be introduced by the analogue of Eq. (96) β eff ( ψ (2) , t w ) = ∂χ (2 , ~k =0 ( ψ (2) , t w ) ∂ψ (2) . (104)The question, now, is whether the two β eff defined by Eqs. (96) and (104) are consistent or not, that is whether theequality lim t w →∞ β eff ( ψ (2) , t w ) = lim t w →∞ β eff ( ψ (1) , t w ) (105)8 FIG. 4: (Color online). The quantity ψ (2) is plotted against t − t w for t w = 100 , , , , , , d Ising system of size 1800 quenched to T = 2 ( T c ≃ . t w height of the plateau ( y = (1 − q EA ) ). In the inset the quantity ψ (1) is shown (same t w of the main figure). holds or not. This is a difficult question to answer in general, we shall limit to the consideration of the particularcoarsening process analysed above in the linear case.The short and the large time scales, in terms of ψ (2) , correspond to ψ (2) smaller or larger than (1 − M eq ) ,respectively. Writing χ (2 , ~k =0 as the sum of two pieces, as in Eq. (97), χ (2 , ~k =0 ( ψ (2) , t w ) = χ (2 , st ( ψ (2) ) + χ (2 , ag ( ψ (2) , t w ) (106)the same considerations made on χ (1 , st apply exactly to χ (2 , st , since this is an equilibrium contribution. Namely,after obeying Eq. (103) in the short time regime, saturates to the equilibrium value (1 − M eq ) and then remainsconstant in the large time regime. For χ (2 , ag ( ψ (2) , t w ) there are no previous results to rely on. We have, then,measured χ (2 , ag ( ψ (2) , t w ) numerically in the quench of a two dimensional Ising model below T C evolving with Glauberdynamics. The aging contribution χ (2 , ag has been isolated using the method based on the no-bulk-flip dynamicsdiscussed in [23, 25, 26, 27].In order to check if a scaling form of the type (98) is obeyed χ (2 , ag ( ψ (2) , t w ) = t − a w F ( t/t w ) (107)we have plotted χ (2 , ag ( t, t w ) for a fixed value of t/t w against t w , as shown in the inset of the upper panel of Fig. 5.From the observed power law behavior we have extracted the exponent a , finding values in the range [0 . − . t a w χ (2) ag ( t, t w ) against t/t w . For large t w the collapse (Fig. 5)is quite good, confirming that the scaling form (107) is obeyed with an exponent a ≃ . − .
62. For small values of t w and t/t w deviations are observed due to preasymptotic effects, similarly to what was already observed in the linearcase [23, 27]. Notice, also, that for large t/t w one has a power law decay of the scaling function F ( t/t w ) ∼ ( t/t w ) − a with the same exponent a entering Eq. (107), exactly as it was observed in the linear case [27].9 FIG. 5: (Color online). t a w χ (2 , ag ( t, t w ) is plotted against t/t w for t w = 100 , , , , , , , T = 1 , . , a is extracted as discussed in thetext, finding a = 0 .
61 for T = 1, and a = 0 .
62 for T = 1 . T = 2. The dashed lines represent the asymptotic behavior χ (2 , ag ( t, t w ) ∼ ( t/t w ) − a for large t/t w . In the inset the quantity χ (2 , ag ( t, t w ) is plotted against t w with t/t w fixed. FIG. 6: (Color online). The parametric plot of χ (2 , k =0 ( t, t w ) against β ψ (2) ( χ (1 , i ( t, t w ) against βψ (1) in the inset) is shownfor T = 2 and t w = 100 , , , , , , , In conclusion, like in the linear case, the existence of the scaling behavior (107) with a > t w →∞ χ (2 , ~k =0 ( ψ (2) , t w ) = (cid:26) β ψ (2) , for ψ (2) ≤ (1 − M eq ) β (1 − M eq ) , for (1 − M eq ) < ψ (2) . (108)The approach to this asymptotic behavior is shown in Fig. 6. Hence, using the definition (104), we findlim t w →∞ β eff ( ψ (2) , t w ) = (cid:26) β, for ψ (2) ≤ (1 − M eq ) , for (1 − M eq ) < ψ (2) . (109)The comparison with Eq. (100) suggests that the consistency condition (105) is satisfied. VIII. CONCLUSIONS
In this paper we have derived the off equilibrium FDR of arbitrary order for systems evolving with Markovianstochastic dynamics. The main effort has been to put the FDR in the same form for both continous and discretespins. In order to stress the generality of the result, we have also shown how the whole hierarchy of FDR can bemade to descend from the fluctuation principle. Once the FDR are available, response functions of arbitrary orderare expressed in terms of unperturbed correlation functions of observables. The payoff is in the development of simpleand efficient zero field algorithms for the numerical simulations.As an application, we have considered the problem of detecting the existence of a growing length in those cases,like in glassy systems, where standard methods based on two-point correlation functions and the corresponding linearresponse functions are of no use. In these cases the simplest object carrying useful information, in principle, would1be a four-point correlation function which, however, is not directly accessible to experiment. Instead, experimentallyaccessible are the nonlinear response functions involving the four-point correlation function through the nonlinearFDR. The choice of which response function and, therefore, of which FDR to use is not univocal, once the realm ofthe nonlinear response functions is entered. Then the choice is matter of convenience. We have made the proposalto use the second order response of a two-point correlation function, rather than the third order response of themagnetization, as advocated elsewhere in the literature. We have, then, demonstrated the numerical advantage of ourchoice through the implementation of the zero field algorithm.Finally, we have made a first step into the important but difficult problem of the definition of the effective tem-perature through the nonlinear FDR. We have considered the domain coarsening process ensuing the quench of aferromagnet below its critical point. Indeed, in that case we have found that it is possible to extract from the nonlin-ear FDR an effective temperature which is consistent with the effective temperature obtained from the much studiedlinear FDR.
IX. APPENDIX I
Proceeding like in the derivation of Eq. (26), the third order derivative of the propagator is given by δ ˆ P h ( t | t w ) δh j ( t ) δh j ( t ) δh j ( t ) =ˆ P h ( t | t M ) ∂ ˆ W ( t M ) ∂h j M ( t M ) ˆ P h ( t M | t I ) ∂ ˆ W ( t I ) ∂h j I ( t I ) ˆ P h ( t I | t m ) ∂ ˆ W ( t m ) ∂h j m ( t m ) ˆ P h ( t m | t w )+ ˆ P h ( t | t M ) ∂ ˆ W ( t I ) ∂h j I ( t I ) ˆ P h ( t I | t m ) ∂ ˆ W ( t m ) ∂h j m ( t m ) ˆ P h ( t m | t w ) δ j M ,j I δ ( t M − t I )+ ˆ P h ( t | t M ) ∂ ˆ W ( t M ) ∂h j M ( t M ) ˆ P h ( t M | t I ) ∂ ˆ W ( t I ) ∂h j I ( t I ) ˆ P h ( t I | t w ) δ j I ,j m δ ( t I − t m )+ ˆ P h ( t | t ) ∂ ˆ W ( t ) ∂h j ( t ) ˆ P h ( t | t w ) δ (12) δ (23) (110)where t M = max( t j ), t m = min( t j ), t m ≤ t I ≤ t M , j M , j I , j m are the sites where the field acts at the times t M , t I or t m , respectively, and δ ( np ) = δ j n ,j p δ ( t n − t p ). Inserting this into Eq. (24) and using Eq. (53), we get the third orderresponse of the first moment R (1 , i ; j j j ( t, t , t , t ) = (cid:18) β (cid:19) n ∂ ∂t M ∂t I ∂t m M ij M j I j m ( t, t M , t I , t m ) − ∂ ∂t M ∂t I h ˆ σ i ( t )ˆ σ j M ( t M )ˆ σ j I ( t I ) ˆ B j m ( t m ) i − ∂ ∂t M ∂t m h ˆ σ i ( t )ˆ σ j M ( t M ) ˆ B j I ( t I )ˆ σ j m ( t m ) i− ∂ ∂t I ∂t m h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j I ( t I )ˆ σ j m ( t m ) i + ∂∂t m h ˆ σ i ( t ) ˆ B j M ( t M ) ˆ B j I ( t I )ˆ σ j m ( t m ) i + ∂∂t I h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j I ( t I ) ˆ B j m ( t m ) i + ∂∂t M h ˆ σ i ( t )ˆ σ j M ( t M ) ˆ B j I ( t I ) ˆ B j m ( t m ) i− h σ i ( t ) B j M ( t M ) B j I ( t I ) B j m ( t m ) i o + β n(cid:16) ∂∂t m h ˆ σ i ( t )ˆ σ j I ( t I ) ˆ B j I ( t I )ˆ σ j m ( t m ) i + ∂∂t m ∂∂t I M ij M j I j m ( t, t M , t I , t m ) − h ˆ σ i ( t )ˆ σ j I ( t I ) ˆ B j I ( t I ) ˆ B j m ( t m ) i− ∂∂t I h ˆ σ i ( t )ˆ σ j M ( t M )ˆ σ j I ( t I ) ˆ B j m ( t m ) i (cid:17) δ j M j I δ ( t M − t I ) + (cid:16) ∂∂t M h ˆ σ i ( t )ˆ σ j M ( t M )ˆ σ j I ( t I ) ˆ B j I ( t I ) i + ∂∂t M ∂∂t m M ij M j I j m ( t, t M , t I , t m ) − h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j I ( t I ) ˆ B j I ( t I ) i− ∂∂t m h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j I ( t I )ˆ σ j m ( t m ) i (cid:17) δ j I j m δ ( t I − t m ) o + β (cid:16) ∂∂t M ij ( t, t ) − h ˆ σ i ( t ) ˆ B j ( t ) i (cid:17) δ j j δ j j δ ( t − t ) δ ( t − t ) . (111)2At stationarity this becomes R (1 , ij j j ( t, t , t , t ) = β n ∂ ∂t M ∂t I ∂t m M ij M j I j m ( t, t M , t I , t m ) − ∂ ∂t M ∂t m h ˆ σ i ( t )ˆ σ j M ( t M ) ˆ B j I ( t I )ˆ σ j m ( t m ) i− ∂ ∂t I ∂t m h ˆ σ i ( t ) ˆ B j M ( t M )ˆ σ j I ( t I )ˆ σ j m ( t m ) i + ∂∂t m h ˆ σ i ( t ) ˆ B j M ( t M ) ˆ B j I ( t I )ˆ σ j m ( t m ) i o + β n ∂∂t m h ˆ σ i ( t )ˆ σ j I ( t I ) ˆ B j I ( t I )ˆ σ j m ( t m ) i + ∂∂t m ∂∂t I M ij M j I j m ( t, t M , t I , t m ) o δ j M j I δ ( t M − t I )+ β ∂∂t M ij ( t, t ) δ j j δ j j δ ( t − t ) δ ( t − t ) . (112)It should be recalled that the singular terms in the last two equations are present only in the Ising spin case. X. APPENDIX II
The expansion of the left hand side of Eq. (67) can be done in two steps. Expanding first the exponential we get (cid:28) ϕ i ( t F ) exp (cid:26) − β Z t F t dt h ( t ) ˙ ϕ ( t ) (cid:27)(cid:29) I → β, [ h ( t )] = ∞ X m =0 ( − β ) m m ! X j ...j m Z t F t dt ...dt m h ϕ i ( t F ) ˙ ϕ j ( t ) ... ˙ ϕ j m ( t m ) i I → β, [ h ( t )] h j ( t ) ...h j m ( t m ) (113)and expanding the individual averages h ϕ i ( t F ) ˙ ϕ j ( t ) ... ˙ ϕ j m ( t m ) i I → β, [ h ( t )] = ∞ X p =0 p ! X q ...q p Z t F t dt ′ ...dt ′ p δ p h ϕ i ( t F ) ˙ ϕ j ( t ) ... ˙ ϕ j m ( t m ) i I → β, [ h ( t )] δh q ( t ′ ) ...δh q p ( t ′ p ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 h q ( t ′ ) ...h q p ( t ′ p ) (114)all together these two contributions give (cid:28) ϕ i ( t F ) exp (cid:26) − β Z t F t dt h ( t ) ˙ ϕ ( t ) (cid:27)(cid:29) I → β, [ h ( t )] = ∞ X m =0 ∞ X p =0 ( − β ) m m ! p ! X j ...j m X q ...q p Z t F t dt ...dt m dt ′ ...dt ′ p δ p h ϕ i ( t F ) ˙ ϕ j ( t ) ... ˙ ϕ j m ( t m ) i I → β, [ h ( t )] δh q ( t ′ ) ...δh q p ( t ′ p ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 h j ( t ) ...h j m ( t m ) h q ( t ′ ) ...h q p ( t ′ p ) . (115)Reorganizing the double sum ∞ X m =0 ∞ X p =0 ( − β ) m m ! p ! X j ...j m X q ...q p = ∞ X n =0 n X p =0 ( − β ) n − p ( n − p )! p ! X j ...j n − p X q ...q p (116)3the above result can be rewritten as (cid:28) ϕ i ( t F ) exp (cid:26) − β Z t F t dt h ( t ) ˙ ϕ ( t ) (cid:27)(cid:29) I → β, [ h ( t )] = h ϕ i ( t F ) i I → β, ∞ X n =1 n ! n X p =0 ( − β ) n − p (cid:20) n !( n − p )! p ! (cid:21) X j ...j n − p X q ...q p Z t F t dt ...dt n − p dt ′ ...dt ′ p δ p h ϕ i ( t F ) ˙ ϕ j ( t ) ... ˙ ϕ j n − p ( t n − p ) i I → β, [ h ( t )] δh q ( t ′ ) ...δh q p ( t ′ p ) (cid:12)(cid:12)(cid:12)(cid:12) h =0 h j ( t ) ...h j n − p ( t n − p ) h q ( t ′ ) ...h q p ( t ′ p ) . (117)Notice that the combinatorial factor in the square bracket gives the number of the distinct permutations among thetwo sets of indeces ( j , ..., j n − p ) and ( q , ..., q p ).Going over to the right hand side of Eq. (67) and introducing the shorthand D P I ( ϕ ( t F )) e β H ( ϕ ( t F )) ϕ i ( t ) E β, → β, [ e h ( t )] = h RHS i β, → β, [ e h ( t )] (118)one gets h RHS i β, → β, [ e h ( t )] = h RHS i β, + ∞ X n =1 n ! X j ...j n Z t F t dt ...dt n δ n h RHS i β, → β, [ e h ( t )] δh j ( t ) ...δh j n ( t n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 h j ( t ) ...h j n ( t n ) (119)and, since e h j ( e t ) = h j ( t ), this can be rewritten as h RHS i β, → β, [ e h ( t )] = h RHS i β, + ∞ X n =1 n ! X j ...j n Z t F t dt ...dt n δ n h RHS i β, → β, [ e h ( t )] δ e h j ( e t ) ...δ e h j n ( e t n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h =0 h j ( t ) ...h j n ( t n ) (120)where h·i β, stands for the equilibrium average at the temperature β and without external field. Therefore, comparingwith Eq. (117), one arrives at Eqs. (68) and (69). [1] For a review see A.J.Bray, Adv.Phys. , 357 (1994).[2] D.A. Huse, J.Appl.Phys. , 5776 (1988).[3] J.P.Bouchaud and G.Biroli, Phys.Rev.B , 064204 (2005).[4] L. Berthier, G. Biroli, J.P. Bouchaud, L.Cipelletti, D. El Masri, D. L Hˆote,F. Ladieu and M. Pierno, Science , 1797(2005).[5] L.F.Cugliandolo and J.Kurchan, Phys.Rev.Lett. , 173 (1993); L.F.Cugliandolo and J.Kurchan, J.Phys.A: Math.Gen. , 5749 (1994); Philos.Mag. , 501 (1995); S.Franz, M. M´ezard, G.Parisi and L.Peliti, Phys.Rev.Lett. , 1758 (1998);J.Stat.Phys. , 459 (1999); G.Parisi, F.Ricci-Tersenghi and J.J.Ruiz-Lorenzo, Eur. Phys. J. B , 317 (1999).[6] L.F.Cugliandolo, J.Kurchan and G.Parisi, J.Phys. I France, , 1641 (1994).[7] G.Semerjian, L.F.Cugliandolo and A.Montanari, J.Stat.Phys. , 493 (2004).[8] J.Kurchan, J.Phys. A , 3719 (1998). For e reviews see F.Ritort, Nonequilibrium fluctuations in small systems: Fromphysics to biology , arXiv: 0705.0455v1; U.Marini Bettolo Marconi, A.Puglisi, L.Rondoni and A.Vulpiani, Phys.Reports , 111 (2008).[9] The zero field algorithm based on the FDR (56) was presented in E.Lippiello,F.Corberi and M.Zannetti, Phys.Rev.E ,036104 (2005). Zero field algorithms, different from ours, were previously proposed by C.Chatelain, J.Phys.A , 10739(2003) and F.Ricci-Tersenghi, Phys.Rev.E , 065104(R) (2003).[10] L.F.Cugliandolo, J.Kurchan and L.Peliti, Phys. Rev. E , 3898 (1997).[11] L.P.Kadanoff and J.Swift, Phys. Rev. , 310 (1968); K.Kawasaki, Phase Transitions and Critical Phenomena vol. 2, ed.C.Domb and M.S.Green, p.443 (New York, Academic, 1972). [12] See for instance Ref. [7], where the derivation of the Onsager relation is carried out for continous variables, but it worksexactly in the same way also for Ising spins.[13] R.J. Glauber, J.Math.Phys. , 294 (1963).[14] E. Lippiello, F. Corberi, A. Sarracino, and M. Zannetti, Phys.Rev. B , 212201 (2008).[15] C. Donati, S.C. Glotzer and P. Poole, Phys.Rev.Lett. , 5064 (1999); S. Franz, C. Donati, G. Parisi and S.C. Glotzer,Phil.Mag.B , 1827 (1999); S. Franz and G. Parisi, J.Phys.:Condens.Mat. , 6335 (2000). See also Ref. [3] for a discussion.[16] See L.F.Cugliandolo, Heterogeneities and local fluctuations in glassy systems , cond-mat/0401506v1 and references quotedtherein.[17] T. Jorg, J. Lukic, E. Marinari and O.C. Martin, Phys.Rev.Lett. , 237205 (2006).[18] H.G. Katzgraber, L.W. Lee and A.P. Young, Phys. Rev. B , 014417 (2004); H.G. Katzgraber, L.W. Lee and I.A. Camp-bell, Phys.Rev.B 75, 014412 (2007).[19] H. Rieger, B. Steckemetz and M. Schreckenberg, Europhys. Lett. , 485 (1994); H.G. Katzgraber, and I.A. Campbell,Phys.Rev.B , 014462 (2005).[20] S.Fielding and P.Sollich, Phys.Rev.Lett. , 050603 (2002); P. Calabrese and A. Gambassi, J.Stat.Mech., P07013 (2004).[21] D.S. Fisher and D.A. Huse, Phys.Rev. B , 373 (1994).[22] J.P.Bouchaud, L.F.Cugliandolo, J.Kurchan and M.Mezard, in Spin glasses and random fields , A.P.Young ed. (WorldScientific, 1998); L.F.Cugliandolo, in
Slow relaxation and non-equilibrium dynamics in condensed matter, Les HouchesSession LXXVII , J.L.Barrat, J.Dalibard, J.Kurchan and M.V.Feigel’man eds., (Springer-Verlag, Heidelberg, 2002).[23] The actual value of this exponent is controversial. Our results can be found in Ref. [24, 25]. What is uncontroversial, andwhat matters here, is that for d > a > , 061506 (2001); F. Corberi, E. Lippiello, and M. Zannetti,Phys. Rev. E , 041113 (2006).[25] F. Corberi, E. Lippiello and M. Zannetti, Phys. Rev. E , 056103 (2005);[26] F. Corberi, E. Lippiello, M. Zannetti, Eur. Phys. J. B (2001), 359; F. Corberi, E. Lippiello, and M. Zannetti, Phys. Rev.E , 041106 (2006); R. Burioni, D. Cassi, F. Corberi, and A. Vezzani, Phys. Rev. Lett. , 235701 (2006); R. Burioni,D. Cassi, F. Corberi, and A. Vezzani, Phys. Rev. E , 011113 (2007).[27] F. Corberi, E. Lippiello and M. Zannetti, Phys.Rev. E68