Stationarization and Multithermalization in spin glasses

We develop further the study of a system in contact with a multibath having different temperatures at widely separated timescales. We consider those systems that do not thermalize in finite times when in contact with an ordinary bath but may do so in contact with a multibath. Thermodynamic integration is possible, thus allowing one to recover the stationary distribution on the basis of measurements performed in a `multi-reversible' transformation. We show that following such a protocol the system is at each step described by a generalization of the Boltzmann-Gibbs distribution, that has been studied in the past. Guerra's bound interpolation scheme for spin-glasses is closely related to this: by translating it into a dynamical setting, we show how it may actually be implemented in practice. The phase diagram plane of temperature vs"number of replicas", long studied in spin-glasses, in our approach becomes simply that of the two temperatures the system is in contact with. We suggest that this representation may be used to directly compare phenomenological and mean-field inspired models.Finally, we show how an approximate out of equilibrium probability distribution may be inferred experimentally on the basis of measurements along an almost reversible transformation.


I. INTRODUCTION
Glasses are dynamical objects, the properties that are relevant for them such as large viscosity and aging, are essentially dynamical in nature. Somewhat surprisingly, a fruitful theoretical approach to them has been to study the proprieties of the energy landscape through the equilibrium properties: this led to the Parisi scheme, and the discovery of hierarchical organization of states-within-states that it implies. On the other side, the mean-field solution of the dynamics starting from a random configuration is largely well-understood: it is characterized by the emergence of widely separated timescales, each with a different characteristic temperature: a situation we shall denote as 'multithermalization'. The Parisi and dynamic multi-thermalization constructions are such that, even if we do not know the actual solution of any finite-dimensional glass model, we do know what both would imply for it.
The purpose of this paper is to establish a stronger connection between static solution and dynamic multithermalization. In section II we review the properties of a system in contact with a bath having different temperatures in the limit of widely-separated timescales: in short, a 'multibath'. We establish the notion of 'multi-thermalization', as the condition in which the system is stationary and has, for all its observables, the same fluctuation-dissipation temperature as the bath, this equality holding at each timescale. A system which would equilibrate with the fastest of these baths (a liquid, a paramagnet) will also multithermalize, but the multibath will in itself generate some slow tails of correlation function syncronized with it. These clearly disappear when the coupling to the bath is weak. A more dramatic situation is known to arise in mean-field glass models: the system, starting form a high temperature configuration, never becomes stationary after being placed in contact with a low-temperaure bath: it 'ages'. Instead, a weak multibath may make it stationary: how weak it may be to do so depends on its timescale -the longer the timescale the weaker the bath needed. We discuss this situation in detail in Sections II A and II B. Beyond mean field models, one may prove that the multi-thermalization situation is still valid, provided the system satisfies in equilibrium a Parisi scheme (although we do not have at present any model for which we may prove that this happens): we know this by extending trivially the result of Franz et al. [19]. This point will be discussed in Sec. IV B. In Section V we show how to extend the classical thermodynamical notion of reversible transformation into a multi-reversible one. This allows us, both in theory and in practice, to transform multi-reversibly a system and then infer the probability distribution it follows at each step. We come back to this in Sec VIII, where we show how this procedure may be implemented (at least numerically) in a simulation of a realistic structural glass. Thus, from the dynamic Fluctuation-Dissipation data one may reconstruct the distribution, which is a generalization of the Boltzmann-Gibbs one. Note that, for this to be the case, it is necessary that the system admits multi-thermalization at each step of the multi-reversible transformation. At this stage, one recognizes that the constructions we are using are closely related to the construction that Guerra [14] used to prove a bound on the free energy of the Sherrington-Kirkpatrick and other models. In Section VII, we uncover the dynamical content of Guerra's scheme: this gives a physically appealing -and numerically realizable -implementation of the procedure.

II. MULTIBATH AND THERMALIZED DISORDER
In a spin-glass system such as the Sherrington-Kirkpartick(SK) model of spins σ i interacting through a random interaction J ij one needs to compute the averaged logarithm of the partition function E ln σ e −βH(σ,J) (2) where E denotes expectation with respect to both the two-body interactions J ij , and the external magnetic fields J i . The former expression provides in fact the correct generating functional for the quenched moments of the Hamiltonian where the quenched measure is defined, for an observable A, as E σ A(σ, J) e −βH(σ,J) σ e −βH(σ,J) .
In other words at fixed disorder J a Boltzmann-Gibbs computation is made on σ and later averaged on the disorder. The (easier to compute) annealed measure has instead a different physical significance, and corresponds on a standard Boltzmann-Gibbs computation on (J, σ): .
inferring (and guessing with a suitable ansatz) the expression (2) by means of continuation for n → ζ, a real positive number, and its limit when ζ → 0. This was accomplished by Parisi in a remarkable series of papers [27,32]. Integrals like (5) are ubiquitous in Parisi's construction, as intermediate steps, with generic n. One may ask if there is a more physical way to interpret them. Indeed this is so. As we shall see in detail below, if evolving the spins at temperature T , and, at a much slower rate, evolving the J's at temperature T /n, we precisely obtain the averages generated by (5). This relation between a replica computation and a 'two-bath' computation is just one instance of a deeper and more general one [5,6,10,13,21,36].
Consider a system with two sets of variables x 1 and x 2 and Hamiltonian H(x 1 , x 2 ). We assume that the variable x 2 reaches equilibrium with relaxation time τ 2 while in contact with a thermal bath at temperature T 2 . Similarly the variable x 1 has relaxation time τ 1 at temperature T 1 .
In the limit τ 2 τ 1 we may formalize an equilibrium theory for the system described by the following thermodynamic functions. Setting β 1 = 1/KT 1 , β 2 = 1/KT 2 (where K is the Boltzmann constant) the free energy is obtained in two steps : The previous thermodynamic expressions, leading to a nested Gibbs-Boltzmann structure, assume that one evolution is adiabatic with respect to the other and that both, the slow and the fast one, have enough time to reach equilibrium.
Remarks: For T 1 = T 2 this whole construction reduces to the standard Gibbs-Boltzmann measure since ζ 1 = 1. Identifying x 2 with σ and x 1 with J and choosing ζ 1 = n we find (5). Moreover for real ζ 1 we recover the quenched free energy in the limit ζ 1 = 0 while ζ 1 = 1 corresponds to the annealed case.
The construction may be generalized to an arbitrary number of timescales r > 2, with their corresponding variables and temperatures, such that each evolution is adiabatic with respect to the previous. Namely we consider an Hamiltonian H(x 1 , x 2 , ..., x r ) where the degree of freedom x a has relaxation time τ a and is contact with a bath β a . Assuming widely separated timescales τ r τ r−1 . . . τ 1 one obtain the full measure recursively. In terms of free energy starting from F r = H we define for any 1 ≤ a ≤ r. Defining for any a = 1, . . . , r the parameter ζ a = β a /β r then the free energy at the final step can be written as where x r e −βrH(x1,x2,x3,...,xr) It is straightforward to show that the equivalent iterative expression for the generating functional (also called pressure in the mathematical literature): turns out to satisfy We will call multibath measure the measure induced by the generating functional P 0 .
Generating functionals of this kind were introduced by Parisi and Virasoro [31], without connection to any dynamics, as a concrete way to construct order parameters conjugate to replica symmetry breaking. It was later discussed in a dynamic context in [9,13], and more recently in [5]. Moreover the recursion (13) is the core of Guerra's interpolation scheme [14] (see section VII). We shall derive dynamically the above measure in detail in two examples below.
A. Thermalized external fields Referring to the above notation we will consider in this section x 1 , ..., x r as magnetic fields and x r+1 as the spin variables. To this purpose we consider a system made of N interacting spins σ = (σ i ) i≤N , σ i ∈ R , with quenched twobody coupling (J ij ) i,j≤N and also coupled with a family of dynamic external fields (J a i ) a≤r i≤N through the Hamiltonian where γ, γ a are non negative real parameters, A is a large positive constant that forces σ i ∼ ±1 .
We use the notation J = (J 1 , . . . , J r ) and J a = (J a i ) i≤N for any a = 1, . . . , r. In the Hamiltonian we kept the explicit dependence of only the dynamical variables (σ, J) for fixed realization of (J ij ) i,j≤N .
The main assumption on the dynamics is that the degrees of freedom (σ, J) have widely separated timescales: denoting by τ σ , τ a the relaxation times of σ and J a respectively, we assume that τ 1 >> τ 2 >> . . . >> τ r >> τ σ . Clearly these scales may depend on the size of the system and become infinite as N → ∞. The dynamic is described by a system of r + 1 Langevin equations: From now on we will set the Boltzmann constant equal to 1. Here T = 1 β and T a = 1 βa , for a = 1, . . . , r, are the temperatures of the thermal bath of the spins and the fields respectively. For all practical purposes one can also think σ i evolving instead following a Glauber or Monte Carlo dynamic with energy (14).
Our aim is to show that the stationary measure of such as dynamical system coincides with the multibath previously introduced. Let us start with the case r = 1. The spins have temperature T 2 = T and timescale τ σ , in addition there is a single family of external fields J ≡ J 1 at temperature T 1 = T with timescales τ 1 τ σ . On a short time-scale compared to τ 1 the spins evolve while the J's are nearly constant. Hence on timescales τ 1 τ τ σ the solution of (15a) is the usual Gibbs measure given J: On the other hand, (15b) is linear, and its solution is Using the assumption of adiabaticity, we may substitute σ i by its fast-time average given by the measure (16): which depends on time through J i . We get: We now use the identity to transform equation (19) into: This is a Langevin equation with temperature T and potential −T ln Z. In equilibrium, it leads to the distribution: where N is the normalization factor and ζ = T /T . Thus one obtains the multibath measure generated by (8) Notice that the term 1 2 (J 1 i ) 2 in (14) carries in (23) as a centered Gaussian measure with variance T . The general case with several J a i with nested timescales is obtained by iteration, i.e. keeping at each step some variables as constants, and identifying the conditional distribution µ(J a |J a−1 , ..., J 1 ). For a given a the free energy F a = − 1 β ln Z a acts as a potential for J a i , in the sense that in (15b) one can make the substitution Therefore one obtains again a Langevin equation, which leads to the conditioned equilibrium µ(J a |J a−1 , ..., J 1 ).

B.
Correlation and Response for the multibath.
Another way to express the dynamics is to transform the problem with extra fields into a problem with no fields, in contact with a multibath. We start by writing (15b) where we have defined where C a (z) = T a e − |z| τa . Response and correlation satisfy the identities We notice that multibath acts on the spin dynamics as a memory term [10]. Indeed for (15a) one has: where R(z) = a γ 2 a R a (z), the memory kernel, is the response function of the bath. The combined noiseρ i (τ ) = a γ aρ a i (τ ) is correlated as which defines the correlation of the multibath. Inserting (27) in the definition of R and comparing with (31), we obtain for the correlation and response of the multibath the relation where T (z) is constant within every one of the nested scales: T (z) = T a for z τa away from zero and of order one to stay within the time-scale. The function T (z) defines the effective temperature of the multibath [11][12][13]21] Since β a = ζ a β then one can write the inverse effective temperature β(z) = 1/T (z) as Notice that the function x can be viewed as the dynamical analogous of the Parisi order parameter in Guerra's bound as explained in detail in section VII. The physical reason for the fact that effective temperatures increases with increasing timescales was discussed by [6].
Here it is important to remark that although we used a particular set of couplings, leading to exponential decays in time, any bath with this fluctuation-dissipation relation and nested timescales will do for the purposes of this paper. is plotted against the autocorrelation function C(τ, τ ) for a system of N = 2 10 spins with random bimodal coupling constants Jij = ±1/ √ N . The external fields are bimodal variables ± γ 2 1 . The temperature of the spins is T = 1/2 and that of the magnetic field is is T = T /ζ, with ζ = 1/2. Curves averaged over 8 − 19 · 10 3 realisations, for the various cases (in the main figure and in the insets). The perturbing field for the computation of χ is set to H = 0.1. The timescale of the slow bath is τ1 = 10 for all the curves. τ is set to τ = 1500 in order to be in the stationary state. The dotted line is the equilibrium slope −1 while the dashed red one are the expected slope −1/2. Lower inset: comparison, for H = 1, between the case with ζ = 1/2 (green curve, same as in the main figure) and ζ = 2/3 (blue curve). For the latter it is τ = 600. Upper inset: comparison, for H = 1 and ζ = 1/2, between the case with τ1 = 10 (green curve, same as in the main figure) and with τ1 = 100. For the latter it is τ = 6000.
C. The SK model with a multithermalized random field In this section we briefly discuss the equilibrium properties of the static analogous of the Hamiltonian (14), i.e. a SK model with gaussian external fields coupled with different thermal baths. The Ising spins case has already been introduced in [31] where a Parisi-like formula has been obtained within the replica framework. More recently the solution has been extended and rigorously proved for generic spin distribution showing an interesting link with Hamilton-Jacobi PDE framework [23]. Here, we briefly review the simpler case of Ising spins and gaussian external field coupled with an additional thermal bath (r=1). The Hamiltonian of the system can be written as where σ i = ±1 and the J's are independent standard Gaussian random variables. The Hamiltonian (34) can be seen as the static analogous of (14) where we take Ising spins and the Gaussian weight is absorbed in the distribution of the J's variables. The factor 1 √ N in (34) ensures a well defined thermodynamic limit. Notice that H N (σ) is a Gaussian process with covariance where The generating functional of the multibath measure (23) (also called pressure density) is where E 0,1 denotes the average w.r.t the quenched coupling and the fields respectively and If ζ 0 → 0 we recover the SK model in a quenched random field with variance γ 2 1 . For real 0 < ζ < 1 one can prove that, following the same procedure presented in [4], the limiting value may be represented as a Parisi-like variational problem: where and f (q, y; x) satisfy a suitable Parisi's PDE. The infimum is taken over x ∈ X ζ which is the space of distribution functions containing the point ζ in the image. This constraint is due to the fact that external field is not quenched and it is averaged out according to (36). The optimal x solution of (38) represents the limiting distribution of the overlap w.r.t. the multibath measure induced by (36).
In order to explain the relation between the multibath measure and the dynamics described in the previous section we recall the definition of effective temperature. Given the two-time correlations C(τ ) of the spins and the associated response R(τ ) and integrated response χ(τ ) =´τ 0 R(τ − τ )dτ then the effective temperature can be defined by the relation [8,12,17] 1 Now we consider the integrated response as function of the correlation in the large N, τ limit It has been proved in [19] that for a class of spin glass models in finite dimensions and under suitable assumptions (see the discussion in section IV B) the quantity χ(C) provides a direct link between static and dynamics, more precisely wherex is the distribution of the overlap of the system w.r.t. the equilibrium measure. We have studied this dynamical quantities by means of numerical simulations, which we now detail (such description applies to the simulations of Secs. II D,II E,III A as well). We considered a system with bimodal couplings J ij = ±1/ √ N and similarly for the fields J i = ± γ 2 1 . Such bimodal forms are chosen for numerical convenience, we don't expect major differences with respect to the Gaussian case considered insofar. We set N = 10 10 (except for the data in Fig. 4 where N = 10 12 ). We have checked that larger values of N do not yield substantial differences. Starting with a random initial distribution of the spins and of the J's the system is evolved in time with Montecarlo rules using Glauber transition rates: a spin σ i is chosen at random and flipped with probability w(σ i → −σ i ) = (1/2)(1 − tanh ∆E/T ). The same transition rate is used for updating the fields J i , with a different temperature T = T /ζ and an extra factor τ −1 1 realising the slower evolution. The dynamical average . . . is taken over a large number of initial configurations and thermal histories, namely Montecarlo trajectories. Due to this large averaging procedure statistical errors on the data are rather small, typically of the order of the symbols used to draw the data. The most significative source of errors is represented by systematic effects due, e.g., to finite times, specifically τ 1 . The response function is computed routinely by running in parallel a copy of the original system perturbed by a small (in principle H → 0) external magnetic field H (constant is space and time). Since data get noisier for smaller values of H we fix the perturbation to the largest value above which we start seeing a dependence of the results on H. Stationarization is achieved by waiting a sufficiently long time such that both C and χ are observed to depend only on the time difference τ − τ . Specific values of the various parameters are given in the captions of the figures.
The result for T χ vs C are shown in Fig. 1. In main part of the figure a comparison is shown between two cases with a relatively strong and weak field intensity, γ 1 = 1 and γ 1 = 0.1, respectively, for a given value ratio ζ = 1/2 between the spin and field temperatures. In the case with strong field intensity γ 2 1 = 1, after a short FDT regime for C 1 a constant slope roughly of order −ζ is approximately observed in a wide range of C. For smaller values of γ 1 such slope is present only in a narrow range of C after the FDT regime. In the lower inset we compare two cases with strong field but with two different values of ζ. We observe that in both cases the slope of the approximately linear part to the left of the FDT regime roughly agrees with the values of ζ. Fitting in the range x ∈ [0.6, 0.8] we find slopes −0.6 and −0.5 for ζ = 2/3 and 1/2, respectively. In the upper inset we compare two cases with strong field and the same value of the temperatures (ζ = 1/2) but with two different timescales τ 1 of the slow bath. The two curves behave similarly, but the transition from the FDT slope to the non-trivial one on the left part of the plot is more sharp for larger τ 1 .
Assuming that the static-dynamic correspondence discussed above holds also for SK, these features can be related to the overlap distribution properties. For γ 1 γ, i.e. weak external field, we expect that the system and its overlap distribution behaves like a standard SK model at inverse temperature βγ. Similarly for γ γ 1 , or equivalently weak two-body coupling, the system is driven mostly by the multibath measure induced by the external field. Here we expect that the overlap distributions develops a plateau at height ζ, or equivalently T χ vs C has slope −ζ at the origin. Besides that, Figure 1 shows that, for finite γ 1 the extension of the plateau is an increasing function of γ 1 .

D. Thermalized interaction parameters
It is also possible that the set of interaction parameters themselves are slow, multithermalized variables [5,15,26]. Referring to our general notation we choose x 1 = (J ij ) i,j≤N and x 2 = (σ i ) i≤N and the Hamiltonian function for some γ, k > 0. The interacting part of (43) is the same of an SK model but here the J ij are not quenched variables but evolve in time as the spin variables . The σ i follows a Langevin dynamic at temperature T and energy (43) and the J ij evolve with a slower timescale at a different temperature T : The noise ρ ij is centered with variance 2T , where T is the temperature of the second equilibrium bath. We write, as before: If τ 1 τ τ σ we may replace σ i (τ )σ j (τ ) with its average respect to the stationary µ(σ|J) defined as in (16): Therefore the solution of (45) gives the stationary measure for the J: where N is the normalization, matching the definition of multibath measure (23) with The generating functional of the measure is where E 0 denotes the average w.r.t. J ij that are independent Gaussian random variables with variance T /k. Taking γ = 1 √ N and k = T one gets the standard setting of the SK model. If ζ → 0 then (48) gives the quenched pressure in the same spirit of the replica trick (5). The thermodynamic limit of (48) for ζ > 0 has been rigoursly studied in [4,26]. It turns out that (48) is represented as a Parisi-like problem of the form (38) where the infimum is taken over the space of distribution functions that have a jump discontinuity with gap ζ at the origin. Notice that this constraint is harder then than the one obtained in the multithermalized external field case. The reason is that here there are no quenched variables coupled with the spins. Keeping in mind the connection between statics and dynamics provided by (40) and (42), we can conclude that the integrated response function reaches the origin with slope ζ . In simple words, the fact that interactions are at temperature T has the effect of 'killing' all effective temperatures T ef f > T in the original problem. This means that in the fluctuation-dissipation diagram T χ vs C the branch to the left of the point with slope ζ is straight with slope ζ. We have checked this fact with numerical simulations, which have been previously detailed in Sec. II C. The only difference is that what is coupled to the slow bath are the coupling constants J ij , and the fields J i are set to zero. The result of the simulations are showed in Figure 2. In the main figure, curves for τ 1 = 10 2 and three choices of ζ are shown. We have checked that different values of the spin temperature T and of ζ yield similar results. The overall behavior is similar to the one observed in Fig. 1, with an FDT part of slope 1 on the right sector of the plot and a different slope on the left. However in this case one observes much nicer straight behaviors in the latter sector and, in addition, the agreement between the observed slope and the expected one (i.e. −ζ) is much better, except for ζ = 1/4. However, for this value of ζ, the comparison shown in the inset between the cases τ 1 = 10 2 and τ 1 = 10 3 indicates that the discrepancy is due to an insufficient value of τ 1 and that there is a convergence to the expected slope increasing τ 1 . Fitting the data (using τ 1 = 10 3 for ζ = 1/4) in the range x ∈ [0, 0.1] we find slopes −0.68, −0.54, and −0.38, for ζ = 2/3, 1/2, and 1/4, respectively. The curves are straight lines with good approximation: indeed the fitted slope is rather stable upon changing the fit interval in the range x ∈ [0, 0.4], the difference being on the third significant figure.

E. Timescales of a system
Let us take advantage of this construction to discuss the question of timescale within a system. Let us consider the behavior of the autocorrelation function C which, for the SK model thermalized couplings J ij , is shown in Fig. 3 in the case with T = ∞. Suppose the decay up to a plateau q EA is independent of τ 1 (for large enough τ 1 ): q EA is then defined as the Edwards-Anderson parameter for the multithermalized system. For correlations below this value, the decay scales non trivially with τ 1 . One possibility is: for a(τ 1 ) a suitable growing function of τ 1 . We have put to the test the scaling of C by means of numerical simulations.
The results for the SK model with spin temperature T = 1/2 and coupling constants J ij coupled to a slow bath at temperature T = ∞ are reported in Fig. 3. In the upper left panel we see that, for the SK model, the scaling (49) does not work at all. A scaling that indeed seems to work is the following: for b(τ 1 ) another suitable growing function of τ 1 , and g is a decreasing function [22]. This can be observed in the upper right panel of Fig. 3. It turns out that b(τ 1 ) increases as the logarithm of τ 1 , see inset and fit described in the caption.
To understand the meaning of this, following [8] we consider triangles of correlations at three large times t 1 < t 2 < t 3 . Scaling (49) implies: an isomorphism of the sum. The range of values where the 'triangle relation' takes this form is usually called "a timescale', because all times involved are commensurate. Instead, scaling (50) is: which is ultrametricity in time: there are infinitely many timescales. In reference [8] a complete classification of all possibilities for the triangle relations of the form C(t 3 , t 1 ) = f (C(t 3 , t 2 ), C(t 2 , t 1 )) is made, based on the fact that in general f must be an associative function.
Notice that the linear scale of the y axis in the two upper panels of Fig. 3 is only suited for the inspection of the scaling properties for relatively large values of C. For the smallest values of correlation, using a logarithmic y scale ( lower panel) allows one to appreciate that, for with an exponential scaling function f (x) for large x. We have checked that all the results discussed in this section are independent of the number N of spins, provided it is sufficiently large.

A. Stationarization
Some systems never thermalize in the thermodynamic limit: their equilibration time diverges with the size. The simplest case is phase-separation starting from a mixed situation: the domains of the phases grow with time, and take a time that depends on the size of the sample to achieve their stationary situation. Another example is the case of spin-glasses, where the divergence with system size is much more rapid. There are some systems like ordinary structural glasses, of which we do not know if the equilibration time is infinite or just longer than we can measure. In all these cases, an autocorrelation function of any observable has a form that is not time-translational invariant, such as for example in domain growth where C(τ, τ ) = C 1 (τ − τ ) + C 2 (τ /τ ) for τ > τ . Such a situation (with C 2 = 0) is called 'aging'.
Systems of this kind may achieve stationarity once subjected to a random, time-dependent interaction [16], even in the thermodynamic limit. By this we mean that all time correlations and response functions, averaged over randomness, become time-translational invariant (TTI). This does not mean that the system is in thermal equilibrium, but rather in a non-equilibrium stationary state. A very intuitive example of this is the case of thermalized couplings discussed in section II D. The J ij evolve with a slow timescale τ 1 in contact with a bath of a higher temperature than that of the spins, and follows a Monte Carlo with renewals every τ σ . If the fast bath has very low temperature, the spins will be 'trying to optimize" the free energy, but they will be 'chasing' a continuously changing optimum, and will not be able to improve beyond some (τ σ -dependent) level. The dynamics then becomes stationary with a timescale of order τ 1 .
A more subtle situation is that of a system subjected to a linear field , itself slowly evolving as in section II A. Because the low temperature configurations depend strongly on the fields J i ('chaos in field'), again one would expect that the slightest field would stationarize the evolution. However, this cannot be true for arbitrary τ 1 , since for τ 1 small the field merely contributes to the fast bath, and amounts only to a rise in its temperature: aging does not disappear if this change is not strong enough. We thus have to expect stationarity to happen in a region of the plane 1/γ 1 , 1/τ 1 around the corner (0, 0). We have studied numerically this plane Checking for stationarity can be an hard and ambiguous task, particularly for large τ 1 . In order to do that we used the following criterion: we stipulate that the system stationarizes if the energy H becomes time independent and/or the autocorrelation becomes stationary and/or it shows an exponential decay (indeed we noticed that the decay of C as a function of τ − τ for fixed τ is approximately exponential or much slower if the system is stationary or ages). The result of our studies is shown in Fig. 4 which confirms what we expected.
Finally, let us note that more complicated situations are possible. A multibath may partially stationarize a system which originally had autocorrelations decaying from q d to q in a time-translational manner, and from q to zero in an aging one, by enlarging the range of stationarity q d − q [9]. Let us also mention here that some systems, specifically ferromagnets, typically do not stationarize when detailed balance is broken by coupling to different baths [33] or when mechanically driven [34]. The reason seems to be that their effective temperature is infinite.

B. Multithermalization with a multibath
Let us consider a multibath with two-time correlations C, response R and integrated responseχ(τ ) =´τ 0 R(τ )dτ and a stationary system with two-time correlations C and response R and integrated response χ(τ ) =´τ 0 R(τ )dτ . For example one may think of the multibath as realized with external fields and the stationary system being a pure SK model as in section II C. At each time scale we consider the effective temperatures (40) where z = τ − τ is the time difference. We say that the multibath and the stationary system are multithermalized if at each z the temperatures are the same T (z) =T (z), for all the pair of observables of the system used to define C, R. In other words, we have a multi-fluctuation-dissipation relation consistent with that of the bath (32). We shall see in section V that this has strong implications for then equilibrium measure. Let us anticipate when we expect multithermalization to happen. i) Any system with short timescales in contact with a multibath (with suitably separated timescales) develops the scales that are thermalized with those of the multibath. ii) Systems that do not become stationary (they age) in contact with an ordinary bath, may become stationary in contact with a multibath. iii) However, only if the multibath's temperatures coincide with the natural aging temperatures of the system multithermalization may be achieved with minimal energy transport, as we shall see in the dynamical version of Guerra's scheme in Sec. VII. The possibility that the system synchronizes its timescales with those of the bath so as to make temperatures match, exists if the system has reparametrization invariances [7].

C. Work and power of a multibath
Let us compute the work per unit time (power) after the classical definition W = force × velocity. Using the quantities introduced in section II B and (26) we obtain The average in (54) is computed over the dynamics (i.e. over the noises). Since the noise is Gaussian one can rewrites the terms ρ a iσ i dyn using integration by parts. Let us start recalling that the evolution of σ follows (15a): Rewriting the solutions of (55) using the Fourier representation of the Dirac deltâ where the other terms in the exponent does not depend on the Gaussian noise ρ a i . The covariance of ρ a i is given by (28), hence integration by parts leads to Now collecting all terms and keeping in mind the definitions of R, C given in (30) and (31), we get where are the correlation and response functions of the system. We may obtain a more explicit expression by integrating by parts in time, and introducing the (time-dependent) effective temperatures R(z) =β(z)∂ z C(z) and R(z) = β(z)∂ z C(z) so that (60) becomes This indeed looks like a conduction term. Since C = a γ 2 a C a , we may discriminate the work done at each timescale writing W = a W a with where we used the fact the [β −β](τ − τ ) is constant within each timescale. This is an energy per unit time. If we adimensionalize separately every timescale, we get This is energy transferred by the a-bath during a time τ a , and they all vanish upon the multithermalization condition β ≡β . At this point we need to be more precise about what we mean by 'large times' and their structure. If the system is finite, the answer is straightforward: we need that at each timescale the system has had enough time to reach the final distribution of all the 'faster' variables, while the 'slower' ones are still substantially unchanged. All this fixes a hierarchy. What about systems in the thermodynamic limit N → ∞ ?

A. Mean-field systems
Consider the paradigmatic case of the Langevin dynamics on a p-spin glass with energy i1...ip J i1,...,ip σ i1 · · · σ ip for p > 2. This example is interesting because it shows us where things can go differently in the thermodynamic limit. The complexity landscape [20] is constituted as in Fig. 5, the states have a density e N Σ(f ) and stop abruptly at a 'threshold level'. In times of order one, the dynamics age -without becoming stationary -just over the threshold energy of the highest, and overwhelmingly more numerous, states. At times t ∼ e KN the dynamics penetrate down to a K-dependent free energy density below the threshold level. Note that each 'step down' in free energy density takes an exponentially (in N ) time longer than the previous one. At each free-energy, we define the effective temperature 1 T ef f = ∂Σ(f ) ∂f as in Fig. 5. Let us now couple the system to a slow bath of intensity γ 1 , temperature T 1 and timescale τ 1 . This is the dynamic counterpart of a well-studied procedure, see [30] (see also [29] and [31]). For very long τ 1 , longer than any timescale of a finite system, and small γ 1 (again, but not vanishing with N ), the system multithermalizes and we find that it eventually sticks at a level with T ef f = T 1 . Starting from a high energy situation, this takes a long time: the system has to age its way down to the appropriate multi-equilibrium level, and this takes exponentially long in N .
If we consider times of order one with the thermodynamic limit taken first, the system becomes stationary for an arbitrarily weak bath provided its temperature is T 1 ≥ T threshold ef f . It may furthermore multithermalize in times that are large but still do not diverge with N , but if and only if T 1 = T threshold ef f . We have hence discovered that if the thermodynamic limit is taken first, in this case a stationary situation exists with energy densities higher than the equilibrium one, even with a field of low amplitude and timescale of order one. This, we shall claim, cannot happen in finite dimensions.

B. Finite-dimensional systems
In a remarkable paper, Franz et al [19] argued that by using the property of stochastic stability introduced in [1] the quantities calculated out of equilibrium in an aging finite dimensional system with short-range interactions at long times coincide with the ones at equilibrium. For a rigorous study of stochastic stability and its consequences in finite dimensional systems see [2,3]. In particular, this is expected (see [35] for a discussion of this hypothesis) for the susceptibilities [8] associated with a perturbation H → H + i1,...,ir J i1,...,ir σ i1 ...σ ir : The essence of their argument is the following: if we are guaranteed that the dynamics lead in times of order one (finite as N → ∞) to an energy density that coincides with the one of the statics, and this for any field J i1,...,ir , then, by simple derivatives of the target values we obtain the same susceptibilities in an aging system as in an equilibrated one. Hence, static and dynamic susceptibilities coincide at long (but finite in the thermodynamic limit) times. To show the convergence, they use a nucleation argument that shows that metastable states with higher free-energy densities are not possible in finite dimensions. What is important to us here is that under basically the same assumptions, we may apply their arguments to a system under the action of a multibath, to show that the timescale-separations required for the multibath, though large, need not diverge in the thermodynamic limit, if the system is finite-dimensional. For example a term like will have the time to relax in finite times to its asymptotic value (τ a → ∞), provided that γ a is of order one (also an assumption related to stochastic stability [19]). In other words, the situation we met above for the mean-field p-spin model (p > 2), where one needs τ σ that are exponentially long in N to reach true stationarity, cannot happen in finite dimensions.

V. MULTI-REVERSIBLE TRANSFORMATIONS, THERMODYNAMIC INTEGRATION AND MEASURE
The appearance of effective temperatures in aging glassy systems, even in the absence of a multibath, has long been known [12,17]. The real power of the multi-bath appears when we use the fact of multi-thermalization with a slowly evolving bath. This allows us to infer the underlying probability distribution of a system, even in a numerical simulation of a realistic system, as we shall see in Section VIII.

A. Ordinary reversible transformations, thermodynamic integration and Fluctuation-Dissipation relation
Consider a system which depends upon variables that we shall denote collectively by σ. Given two observables A, B the correlation function, denoting by · dyn the dynamic average, is: Given a perturbation of the type H → H − δh(τ )B the response function is where · dyn δh is the average under perturbation. We shall consider transformations of the energy function (via its parameters) that are reversible, by which we mean that: • they are quasi-static: if at any step of the transformation we were to stop and wait until the average values of all observables (over time-windows, or over several copies of the system following the same protocol) does not evolve, and, furthermore, two-time correlations depend exclusively on time-differences.
• the Fluctuation-Dissipation relation holds at each step: In practice, this means that duplicating all times involved does not change the result. Let us show that this process leads to the Boltzmann-Gibbs distribution. The idea is to consider a perturbation δβ β H, so the expectations will be those associated with the measure e −β[H+ δβ β H] or equivalently β → β + δβ. Consider an arbitrary observable A. Then by definition of the response function (69) and by FDT (70) we have that where we have used the clustering property at widely separate times. Because this holds for every A, one can determine µ(σ), the distribution of σ. Indeed we can choose A = δ(σ − σ ) and take the average over σ obtaining Imposing the normalization, this gives the Gibbs-Boltzmann measure. What we have done is to reconstruct this measure by means of a 'thermodynamic integration' in β.

B. Multi-reversible transformations, thermodynamic integration and multi-Fluctuation-Dissipation relation
Suppose now we have a system in contact with a multibath with sufficiently separated timescales. We consider transformations that change the Hamiltonian slowly enough, so that the transformation is multireversible, namely • it is quasi-static: if at any step of the transformation we were to stop and wait, with the multibath still on, the average values of all observables (over time-windows, or over several copies of the system following the same protocol) would not evolve, and, furthermore, two-time correlations and response functions depend exclusively on time-differences.
• the Fluctuation-Dissipation Ratio of all the observables of the system and of the multibath coincide at each timescale with a single β(τ ): This is just the 'multi' version of ordinary reversibility. We shall show, with a procedure that is a direct generalization of the one above, that we obtain the multibath measure introduced in section II.
Let us assume that we evolve T = 1/β, the fast temperature, from T = 0 to any finite temperature. Consider a system with fast σ and slow J variables at inverse temperature β 1 and relaxation time τ 1 . We shall assume multi (bi) thermalization namely with effective temperature β(τ − τ ) = βx(τ − τ ) where for some τ * such that 1 τ * τ 1 . Hence we are assuming that the system spontaneously respects (72), for all observables A, B, namely FDT with temperatures β in the fast timescales, and β 1 in the slow timescales. This corresponds to µ(σ|J) and µ(J), the former being the distribution reached by σ before J had the time to move. These are the distributions we wish to compute. We proceed as above, treating energy as a perturbation, but this time the r.h.s. of (71) can be split into timescales. We choose the time τ * such that σ has performed all its fast relaxation, but J has not had the time to change. Linear response (69) reads as Now, we make the crucial assumption that the time-difference τ − τ * is large enough that σ is able to thermalize at given J, yet small enough that J hasn't changed. Then, the expectation A(τ )H(τ * ) dyn partially clusters, and (74) corresponds to δ A δβ = −ˆdσdJ AH µ(σ|J) µ(J) where r(β) ≡´dσdJ H µ(σ/J) µ(J). Integrating once over σ, and rearranging, we get the two equations: and subtracting them we find The solution of the first equation is µ(σ|J) = g(J, β)e −βH(σ,J) for some g(J, β) that may be fixed by normalization: Plugging this into the second equation, it becomes: which implies: which matches (22). The generalization to r nested timescales is straightforward, one must be able to choose times τ * 1 , . . . , τ * r such that at τ − τ * a the variables J 1 , ..., J a did not have the time to move, while J a+1 , . . . , J r , σ have reached their equilibrium distribution. The effective temperature is a staircase function taking values βζ a for a = 1 . . . , r. As mentioned above, the average (5) was originally only motivated by the value around ζ = 0. At ζ = 1 we have the annealed average, and it is easy to check that there is no transition at any temperature for a spin-glass. Clearly, at low enough temperatures there is a spin-glass phase as ζ → 0. The question then arose of what is the transition line in the ζ − T plane: this was computed for the random energy model by Gardner and Derrida [18] and later by Guerra and Talagrand by rigorous methods [24,26], for the SK model by Kondor using replica approach (see [4,26] for a rigorous treatment). In both cases, the phase diagram looks like the sketch in Fig 6. This may be obtained by a two-temperature multibath, where the spins are at temperature T and the couplings at temperature T = T /ζ, and evolve at a much slower rate τ 1 as explained in section II D. The interpretation of the phase diagram is physically appealing, and may be seen in the sketch of Figure 7 and 8. The values of effective temperatures βx of the system are cut off at the level of T = T /ζ as follows: the slope of the curve χ vs C is approximately independent of ζ from q EA down to the point q min where the tangent dχ dC = − ζ T ; and continues as a straight line down to the minimal value of C (see Fig. 7). Starting from zero and increasing ζ, the transition takes place at the value of ζ = x crit (at which T matches the lowest available effective temperature for the system): at this point the relaxation takes place with only two possible values for x(q) (see Fig. 8). At the transition point one should observe in Fig. 7 two straight lines with slope −ζ and −1, that intersect with each other at the point C = q EA . Within the two-temperature interpretation of this diagram, it seems possible to make a phenomenological description in terms of droplets within the glass phase of this diagram, and this is a framework that would render the different approaches directly comparable.

VII. PHYSICALLY IMPLEMENTING GUERRA'S INTERPOLATION AS A MULTIREVERSIBLE TRANSFORMATION
The multibath measure generated by (13) is the core of the Guerra's interpolation scheme for the SK model [14]. In section II A we showed that a multibath measure can be viewed as a stationary measure for a dynamical system in contact with different thermal baths and widely separated timescales, hence it is natural to look for a dynamical analogous of Guerra interpolation. In this section we investigate this analogy by addressing in particular the following question: is there a dynamical counterpart of the positivity property in Guerra's scheme? We will show that if the system multithermalizes (see section III B for the precise meaning) along the interpolating path then the answer to the previous question is positive thanks to the property of a multi-reversible transformation (described in section V B).
Let us start by briefly sketching Guerra's construction, we refer to the original work for the details [14]. Consider a system of N spins and two independent random Hamiltonian H(σ) andH(σ) with centered gaussian disorder and covariances where q 12 = 1 N i≤N σ 1 i σ 2 i is the overlap. In other words H is the Hamiltonian of the SK model whileH is a gaussian external field. Let t ∈ (0, 1) be an interpolation parameter and r be an integer. Consider a non decreasing sequence q = (q a ) a≤r with q 0 = 0, q r = 1. Let (H a ) 1≤a≤r be a family of i.i.d. copies ofH and define and for t ∈ (0, 1) the interpolating Hamiltonian Then H 1 is the Hamiltonian of the SK model while H 0 contains just one-body interactions. Following Guerra, for H t we assume a multibath measure associated to a given non decreasing sequence ζ = (ζ a ) a≤r . The generating functional or pressure density for this measure is where E averages the quenched variables and P 0 (t) is obtained trough the recursion (13) starting with We denote by t the average w.r.t. the multibath measure induced by (84) and x(q) denotes the discrete distribution associated to the sequences q and ζ . Then one can prove a crucial inequality for any choice of the sequences q and ζ. Integrating both sides of (86) from t = 0 to t = 1 on gets the celebrated Guerra's Replica Symmetry Broken bound [14]. Notice that this procedure can be viewed as a thermodynamic integration. Next we will show how one can obtain an inequality analogous to (86) in the dynamical setting.
A. Dynamic realization Consider a dynamical system with Hamiltonian (14) and set γ = t N and γ a = We assume as before that the quenched couplings J ij are i.i.d. standard gaussian. We write the Hamiltonian as H t = H 0 + H 1 + H 2 : where A is a large constant forcing σ 2 i ∼ 1, so that H 1 is for large A essentially a constant. We recall that H t is Guerra's Hamiltonian (83) apart from the term H 1 that allows us to use Langevin for σ i .
The evolution of σ follows, see eq. (15a), where we have chosen τ σ = 1 and T is the temperature of the main bath. The J evolve according to (15b): Consider now the quantity In complete analogy with Guerra's method, we wish to compute the average of (90) over the dynamics (i.e. over the noises). We shall need the measure for this. For the fields J a i it is simply the Wiener measure while for the σ we complete the measure with a Fourier representation of Dirac-δ for equation (88): dσ e (´dτ ia γaσi(τ )J a i (τ )+γ ij Jijσi(τ )σj (τ )...+ other terms) where we have only specified the terms containing J a i and J ij . We denote by dyn the average with respect to the measure Eˆdσ dσ dJ e −´dτ S(σ,σ,J) + other terms where E is the average on the quenched variables and

B. Integrations by parts
We start noticing that − ∂Ht ∂t dyn is the average of a sum of terms containing the variables J ij and J a i . We will rewrite those terms using integration by parts.

The Jij
The variables J ij are quenched then can use integration by parts for Gaussian vectors. Hence where we have used (88) and (92). Integration by parts in this case means to replace in the average J ij → ∂ ∂Jij , and then we get: There are two interpretations of this term. The most general one is to notice that this is a response to fields acting on products σ i σ j in the links that are coupled. That is, we add a term in the energy ij h ij σ i σ j and compute ij δhij δhij σ i σ j h=0 . Using the fact that correlations among sites vanishes, we rewrite the last term in (96) as where are the correlation and response functions of the system to the perturbation h ij .

The J a i
The core of Guerra's construction is the evaluation of the terms J a i σ i dyn . In order to do this, one should keep in minds formulas (26)- (29). By definition 93 we have that J a σ dyn =ˆdσ dσ dJ e − 1 2Ta τa´d τ (τaJ a +J a −γaσ) 2 +γa´dτ aσ (τ )J a (τ ) J a σ where for simplicity we omit here the subindex i. The integral is Gaussian, we wish again to integrate it by parts. The variation of the exponent is: (100) so that we may replace in the average note the sign in the bracket in the last term, which may be adjusted by integrating by parts. Now we use (29) in equation (101) to write the inverses obtaining Now, reinstating the indices i , keeping in mind that γ 2 a → (1 − t)(q a − q a−1 ) and using (97) we obtain: Going back to (90) and putting all terms together: This is closely analogous to Guerra's expression for the remainder [14]. To see this bear in mind the connection between statics and dynamics discussed in the previous sections. We consider a process, where the parameters are adiabatically varying. We assume the this process is stationary step by step, and satisfies at each step: for two different effective temperatures β(τ ),β(τ ). One can substitute the above relations in (105). As example one can writeˆτ When the multithermalization condition β(τ ) =β(τ ) = βx(τ ) holds for all τ one gets and since dx dτ has by construction a negative sign, the negativity of the Guerra's remainder for dynamical average is obtained in dynamical setting with the assumption of multithermalization. We may now perform thermodynamic integration of the l.h.s., and because of multithermalization we obtain a dynamical version of the Guerra's bound (86).
The relevance of such positivity in a dynamical setting is still to be understood. On one hand, in the equilibrium picture the positivity of the remainder has provided an excellent guide to search for the rigorous proof of the Parisi solution in the mean field case. On the other hand the dynamical setting described here provides a bridge with experimentally accessible computations and thus makes possible to test the robustness of the positivity property also beyond the assumption of multi-reversible thermalization.

VIII. BEING REALISTIC: PRACTICAL MEASURES FOR GLASSES
A realistic glass may be modelled as a system of particles, of different sizes to avoid crystallisation. One may subject such a system to a multibath, by applying uncorrelated fields to each particle, themselves in contact with a slow thermal baths.
On the other hand, several developments in the 90's [10,29,30] based on the Random First Order scenario, pointed to the fact that the structure of an aging glass could be reproduced by the measure (5) with a Hamiltonian with ζ = T T ef f and T ef f a free parameter, adjusted to represent the out of equilibrium system at its age, of the order of what the temperature was at the moment it fell out of equilibrium.
As one can see, ζ is not the only parameter because there is also the intensity γ, and herein lies the entire problem. The auxiliary slow bath selects the states with the appropriate effective temperature, but in order to do so, it needs to have an intensity γ that scales with the rate of escape from those states, their inverse lifetime. In a mean-field situation, as explained in Section IV A, in which states with higher free-energy density have an exponentially large lifetime ∼ e aN , one can let γ → 0 at the end of the calculation, because states have zero escape rate in the thermodynamic limit. In other words, the thermodynamic limit and the γ → 0 limit do not commute. In a finite-dimensional case, where the escape rate is finite, the limit γ → 0 sends us back to the usual Gibbs-Boltzmann measure, and we get nothing. In other words, the thermodynamic limit and the γ → 0 limit do not commute.We need a value of γ that is as small as possible, but large enough to compensate for the escape rates (i.e. the finite lifetime) of the states. Clearly, the construction is not without ambiguities.
When we work with a multibath, we make the same construction, and of course we have the same problem with γ. However, here we have a direct experimental test of our assumptions. Consider the following protocol: we let our system age. Assume that at time t w it is still out of equilibrium. Regardless of whether it will eventually equilibrate or not, we wish to characterize the measure that describes it such as it is. We check by measuring correlation and response that the system has a two-temperature fluctuation-dissipation behavior, a fact that is well-attested numerically, at least as a good approximation. Now we apply a weak multibath such that: i) timescale is of the order of the α correlation decay (from q EA to zero) scale of the glass, ii) it has the same temperature as the fluctuation-dissipation one of the system, and crucially, iii) its amplitude is just sufficient so that we verify that the system becomes stationary by virtue of its interaction with the bath (the α timescale ceases to grow, as it does in an aging system). Next, we slowly change parameters of multibath, temperature and intensity γ, always verifying that the timescale of the slow bath is of the order of the α timescale, and that the effective slow temperatures of bath and system are the same. If such a procedure is possible, and we can take the system to the liquid situation in which it is in ordinary equilibrium, then a (multi)thermodynamic integration is legitimate, and we have in effect experimentally proven, following the results in the previous sections, that the system as it was at the 'age' at which we started, may indeed be described by the multibath measure, with the amplitude we needed to ascribe to it so that the system remained stationary from the moment we connected the multibath, and was multithermalized by it.
Note that it is both the parameters T ef f and γ that play a role. Indeed, we may describe this as a (Gedanken) experiment to measure these two parameters.

IX. CONCLUSIONS
The Parisi construction, although often described as a solution, is in fact something wider: a symmetry-breaking scheme [32]. In this sense it is more akin to the general ideas of ferromagnetism from Curie to Landau, than to the Onsager solution for the Ising model. This is why we may ask if it applies to other problems, such as finite-dimensional spin glasses. Indeed, the fact that it is a scheme where a symmetry group is broken into symmetry subgroups means that it is possible to propose a solution of this form in any model. The problem is, however, that the symmetry in question (replica symmetry) is a very bizarre one, and is too closely dependent on one particular formalism.
On the other hand, the dynamics with many timescales -and one temperature scale at each -can be thought of as as a symmetry-breaking situation as well. Consider first ordinary thermalization of Hamiltonian dynamics. When a system is thermalized at a given temperature and then isolated, it has a dynamic time-reversal symmetry with temperature as its parameter, that implies the fluctuation-dissipation and Onsager reciprocity relations. The Hamiltonian dynamics itself, before choosing a temperature, had a larger group of symmetry, and thermalization can be seen as the act of breaking down this symmetry to a subgroup labeled by the temperature. Multithermalization of widely separated timescales corresponds to a more complicated breaking of the large symmetry group, with a parameter for each timescale. This may seem an unnecessarily pompous and abstract way of putting things, but, again, it allows us to see that any system with slow dynamics and widely separated timescales may possess a solution with multithermalization.