Mode-Dependent nonequilibrium temperature in aging systems

We introduce an exactly solvable model for glassy dynamics with many relaxational modes, each one characterized by a different relaxational time-scale. Analytical solution of the aging dynamics at low temperatures shows that a nonequilibrium or effective temperature can be associated to each time-scale or mode. The spectrum of effective temperatures shows two regions that are separated by an age dependent boundary threshold. Region I is characterized by partially equilibrated modes that relax faster than the modes at the threshold boundary. Thermal fluctuations and time-correlations for modes in region I show that those modes are in mutual thermal equilibrium at a unique age-dependent effective temperature $\Theta (s)$. In contrast, modes with relaxational timescales longer than that of modes at the threshold (region II) show diffusive properties and do not share the common temperature $\Theta (s)$. The shift of the threshold toward lower energy modes as the system ages, and the progressive shrinking of region II, determines how the full spectrum of modes equilibrates. As is usually done in experiments, we have defined a frequency-dependent effective temperature and we have found that all modes in region I are mutually equilibrated at the temperature $\Theta (s)$ independently of the probing frequency. The present model aims to explain transport anomalies observed in supercooled liquids in terms of a collection of structurally disordered and cooperative rearranging mesoscopic regions.


I. INTRODUCTION
Finding a general theory for nonequilibrium systems constitutes one of the major challenges in modern physics. Despite of the many different existing approaches to nonequilibrium phenomena 1,2 , at present there is not yet a general theory applicable to systems arbitrarily far from equilibrium. During many decades, physicists have studied different sorts of nonequilibrium systems. A common strategy has been to adapt well-known thermodynamic concepts to the nonequilibrium regime with the hope that the most fundamental assumptions still hold in a varied range of nonequilibrium conditions. One of the concepts that physicists have tried to rescue repeatedly is that of local equilibrium 3 and the possibility of defining a nonequilibrium temperature 4 . Different definitions of a nonequilibrium temperature have been explored for nonequilibrium liquids 5 , steady-state driven systems 6,7 , granular media 8 and glassy systems 9,10 . Admittedly, the concept of a nonequilibrium temperature in realistic systems still remains elusive, mainly because of the scarce number and difficulty of the experiments that have to address fundamental predictions by the theory 11,12 .
Compelling results have been obtained over the past two decades in the context of glassy systems. Several concepts (e.g. mutual equilibration) have been rationalized by borrowing mathematical tools used in the context of spin glasses (such as reparametrization-time invariance or supersymmetry 13 ). Glassy systems are non-stationary systems which relax very slowly toward equilibrium 14 . Relaxational processes can be characterized by the age of the glass, which is the time elapsed since the system was prepared in the non-equilibrium state. A key signature of aging is the fact that correlations and responses decay in a time-scale which is roughly proportional to the age of the system, thereby showing strong and complicated hysteretic effects 15 . Theoretical studies in mean-field spin-glasses have shown 9 that a nonequilibrium temperature can be defined in the aging regime. This entails an extension of the fluctuation-dissipation theorem (FDT) to the aging regime. This nonequilibrium temperature (usually called effective temperature in the glass jargon) is defined by, where C(t, s) is a generic two-times correlation function and G(t, s) is the corresponding response of the system to an external perturbation applied at a given previous time s. At first glance, (1) is just a useful way to quantify deviations from thermal equilibrium (the equilibrium regime is characterized by T eff (t, s) = T ). In experiments 11 the effective temperature is usually defined in terms of quantities measured in frequency space 16 : where Ω is the probing frequency, S(Ω, s) is the power spectrum of fluctuations at time s and χ ′′ (Ω, s) is the dissipative part of the response. In equilibrium T eff = T independently of the probed frequency Ω. Over the past years the validity of this parameter from a statistical and thermometric point of view has been widely studied 17,18 . Yet, there are still several debated issues. Of the utmost importance is to understand under which conditions the effective temperature governs the heat flow between different aging systems, whether a zero-th law is applicable between different systems sharing a common effective temperature 19 and whether such temperature is experimentally measurable 12 . More speculative but not less important is whether we could use such measurements to trace back the age of a man-made glass, i.e. to approximately determine the time elapsed since it was quenched. Our capability to answer such a question might help us to understand the limitations and difficulties in the experimental measurement of effective temperatures.
No doubt that many aspects related to the nonequilibrium properties of glasses are tightly dependent on the microscopic details of the system under study (a polymer blend, a structural glass former, a disordered magnet or a dirty superconductor). However, generic and universal aspects of the nonequilibrium thermal properties in glassy systems are expected beyond the specific properties of the system under study. This has motivated the appearance over the past decades of several phenomenological approaches to the glassy state, from the old free volume theories to the Adam-Gibbs-DiMarzio entropic theories 20 , the free-energy landscape approaches 21 until the more recent mosaic pictures 22,23 .
A major question that has recently attracted much attention in experimental studies of supercooled liquids is the microscopic origin of the non-exponential behavior observed in the relaxation functions. According to some views, a glass is spatially homogeneous and the relaxation of local regions in the system is intrinsically non-exponential. Another view supports the idea of a glass as an intrinsically heterogeneous system made out of a collection of cooperative mesoscopic regions in the system, each one showing exponential relaxation with its own characteristic relaxation time. Ergodicity also implies that regions must have a spectrum of finite lifetimes and that these regions intermittently switch from mobile to frozen. The superposition of the relaxation of the many local regions gives rise to the non-exponential behavior observed in the relaxation of bulk quantities. Concomitant to the experiments there has been an upsurge in the theoretical study of heterogeneities in simple models using analytical methods [24][25][26] and numerical simulations 27 . In general these models support the heterogeneous scenario in agreement with what has been observed in experiments [28][29][30] . If the heterogeneous scenario is valid then one can ask whether it is possible to define a mesoscopic effective temperature locally defined for each region and intermittently fluctuating at time intervals roughly equal to the lifetime of the region. In such case it seems more reasonable to talk about a fluctuating effective temperature field rather than a single well defined value. Recent AFM 29 , fluorescence spectroscopy 30 and Nyquist noise measurements 11,31 seem to endorse the heterogeneous scenario. Moreover, Nyquist noise experiments purport evidence of frequency dependent effective temperatures (2) which can be interpreted in terms of degrees of freedom that have different effective temperatures and not being in mutual thermal equilibrium.
Inspired by these facts, we introduce an exactly solvable model with glassy dynamics, the Disordered Oscillator Model (DOM). The model is an extension of a previous one introduced in Ref. 32 as a simple model for a glass with slow dynamics in the presence of entropy (rather than energy) barriers, the so-called Oscillator Model (OM). The OM represents a simple example of cooperative slow dynamics due to entropy barriers where kinetics is constrained by the dynamical rules rather than by the Hamiltonian. In the present paper we show how a mode dependent effectivetemperature can be defined in the presence of many different time-scales. This is naturally accomplished by adding structural disorder in the OM (a similar strategy has been followed for the Backgammon model 33 ). In the presence of disorder the system preserves aging dynamics at low temperatures and a mode-dependent effective temperature which is higher than the bath temperature drives the relaxation of the system toward equilibrium. The scenario of relaxation in this model is as follows: relaxational modes can be split in two regions (I and II) separated by a threshold boundary. Region I contains mutually thermalized modes at an age-dependent (denoted by s) but frequency independent effective temperature Θ(s) that lies above the bath temperature T . In contrast, modes in region II are non-thermalized and diffusive. Relaxation takes place by the progressive mutual thermalization of degrees of freedom in region I to the single temperature Θ(s). As the system ages, the threshold region shifts, Θ(s) decreases further and more modes enter in region I. Full thermalization occurs when Θ has relaxed to T and region II has vanished.
The paper is organized as follows. In Section II we introduce the model. In Section III we present the formal solution of the model. In Section IV we derive the dynamical evolution for the energy, correlations and responses.
In Section IV we also analyze the behavior of the non-equilibrium effective temperatures for the different modes in the system. In Section V we derive the effective temperature using fluctuation relations. Section VI summarizes the results and Section VII presents a discussion of the implications of this work to our understanding of the glassy state. Technical aspects are described in the Appendices. The original OM model was introduced in 32 as a simple entropic model for a glass. It consists of N identical and uncoupled harmonic oscillators with Hamiltonian: where K is the stiffness or spring constant and the x i are the positions of the oscillators. The dynamics of the model is stochastic of Monte Carlo type. A Monte Carlo trial move is defined as follows. The position of all oscillators x i are simultaneously shifted by x i + r i / √ N where r i are random variables extracted from a Gaussian distribution with zero mean and variance ∆ 2 . The global move is accepted according to a transition probability W (∆E) which satisfies detailed balance: W (∆E) = W (−∆E) exp(−β∆E), where ∆E is the change in the energy (3) and β = 1/(k B T ) (from now on we set the Boltzmann constant equal to 1, k B = 1). We use the Metropolis algorithm for the dynamics.
Note that whenever a move is accepted, all oscillators are updated in parallel which means that the motion of all oscillators is cooperative and coupled by the dynamical rule. Therefore, although they are non-interacting in the Hamiltonian (3), the dynamical rule couples all the oscillators in a non-trivial way.
At low temperatures the OM shows typical non-equilibrium features of glassy systems such as aging in the correlation and response functions. Indeed, at very low temperatures the acceptance of the Monte Carlo updates decreases very fast (like 1/t). In this way relaxation generates entropic barriers which make dynamics extremely slow (the energy decreases as 1/ log t). The simplicity of this model, however, makes it exactly solvable. For a compilation of other results see Refs. ( 16,34,35 ). The OM has also been extended to include non-harmonic generalized potentials 36 .
In the present work we add different time-scales to the system by introducing disorder in the different degrees of freedom that enter into the Hamiltonian. The disordered oscillator model (DOM) is defined by the following Hamiltonian: where Kǫ i is the spring constant of each oscillator. As in the original model, the system consists of a set of noninteracting harmonic oscillators evolving in parallel according to a Monte Carlo dynamics (4). Again, the dynamics is cooperative and couples the evolution of all oscillators. In this disordered version, the spring constants of the oscillators have a probability density given by, We assume that ǫ i > 0, ∀i. The OM (3) is a particular case of the DOM with g(ǫ) = δ(ǫ − 1). In Fig. 1 we show a schematic picture of the model where a set of different oscillators move coupled according to the dynamical rule. The inclusion of disorder in the DOM allows us to explicitly consider nonequivalent degrees of freedom in the system having different relaxational properties. Thermodynamics of the DOM is trivial. The free energy per oscillator f is given by, where g(ǫ) is defined in (6). This expression is valid for any finite N if ǫ i > 0, ∀i. In the large N limit, g(ǫ) may become a continuous function of ǫ. The last integral term in (7) can diverge depending on the behavior of g(ǫ) when ǫ → 0. Whatever value this last term takes, it only represents an additive constant to the total entropy. It may be taken as defining the entropy of a reference state without influencing the rest of thermodynamic quantities. As we will see later, the behavior of the oscillators with coupling constants ǫ → 0 determines the global relaxational properties of the system.

III. FORMAL SOLUTION OF THE MODEL
In this section we derive the dynamical equations for the energy and the acceptance rate of the system. To this end we must compute the probability distribution of having a change δE in the energy E = H 0 /N (5) in one Monte Carlo step. In a Monte Carlo step the N oscillators are moved according to the rule: where r i are random variables Gaussian distributed with zero mean and variance ∆ 2 . The change in energy is given by, In what follows we identify N i=1 by i without risk of confusion. The probability of having a change δE can be written as, Using the integral representation of the Dirac's delta function and performing the Gaussian integrals we obtain: where we have defined the following quantities: where the brackets < .. > denote an average over initial conditions and dynamical histories of the system. Note that (11) is analogous to that obtained in the OM 32 (where ǫ i = 1 , ∀i, and Q k = KE, ∀k). In the DOM the energy is given by E = Q 1 . Comparing to the OM (11) now depends on Q 2 rather than on the energy E. As in the case without disorder, the evolution for the energy is given by the expression: where we have used (4). Straightforward calculations give, where we have defined the following quantities It is easy to see that for ǫ i = 1, ∀i, the results of Ref. 32 are recovered. The acceptance or average fraction of accepted moves is Inserting (11) in (17) we get: To solve (14) and (18) we note that the equations for the energy and the acceptance are not closed because they depend on the quantity Q 2 (t). The equation for Q 2 turns out to depend on Q 3 and Q 2 . In general, the equation for Q k depends on Q k+1 and Q 2 . In Appendix A we present the dynamical equations for Q k . The final result (A7) is: where we have defined the moments of the coupling distribution: In what follows we will identify ǫ = ǫ 1 . We must notice that the hierarchy (19) can be solved in the same way as in the case of the spherical Sherrington-Kirkpatrick model 37 or the Backgammon model 38 by introducing a generating function and solving the resulting partial differential equation by the method of the characteristics. We will not pursue this strategy here, rather we will analyze the simpler adiabatic approximation.
A. Adiabatic approximation: the zero-temperature limit Equation (14) can be solved at zero temperature by means of an adiabatic approximation. It states that some quantities evolve much slower than the energy so we can consider them as partially equilibrated over the constant energy surface. At equilibrium the derivative of the energy (14) is zero and Q eq 2 is given by: The adiabatic approximation becomes asymptotically valid for long enough times. It relates Q 2 (t) with the energy E(t) by assuming that the system has partially equilibrated over the hypersurface of energy E(t): Within the adiabatic approximation (14) is closed and can be solved analytically. In fact, this equation is identical to the equation of the energy in the OM model. As was done in Ref. 32 the time-evolution equation for α(t) (15) at zero temperature is: where we have expanded the error function in the limit α → ∞. For large times, the parameter α grows logarithmically as a function of time: therefore, the acceptance is independent of the disorder distribution g(ǫ) and decays as, plus subdominant logarithmic corrections. The energy decays logarithmically in time and depends only on the mean value of the distribution, ǫ: In Fig (26) while the continuous line is the numerical solution of (14). In the inset we show the acceptance. The dashed line with circles correspond to the adiabatic solution (25) and the continuous line is the exact solution of the numerical equations. To solve (14,18) we have integrated up to 500 terms in the hierarchy (19). (22) we can define a time-dependent effective temperature (that we denote as Θ(t)) given by:

From (21) and the adiabatic approximation
At this stage of the calculations this is just a definition which denotes the temperature that we could assign to the system if it were partially in equilibrium over the hypersurface of constant energy E(t) (meaning that it reaches the maximum entropy state compatible with that energy). In the following sections we will show how this effective temperature drives the relaxation of some degrees of freedom in the system. In addition, we will show that Θ(t) can be exactly derived by assuming that energy fluctuations are described by a non-equilibrium measure reminiscent of the Edwards measure introduced in the context of granular media 39 .

IV. ANALYTICAL SOLUTION OF MODE-DEPENDENT QUANTITIES
In this section we analyze the relaxation of the different degrees of freedom in the system. We define the modeenergy density which is the energy density of those oscillators with a given spring constant value Kǫ, and compute its time evolution. Then we compute the correlations and responses for the different modes. Finally, we study the effective temperature for the different degrees of freedom in the system.

A. MODE-ENERGY DENSITY RELAXATION
A physical quantity which provides deeper insight into the behavior of the system is the mode-energy density w(ǫ, t) defined as follows: This quantity corresponds to the mode-energy density of those oscillators with the same spring constant Kǫ. The total energy of the system can be written in terms of the mode-energy density as: The characteristic relaxation time τ (ǫ) of a given mode is inversely proportional to Kǫ (see the computations in Appendix B). Therefore, w(ǫ, t) is the energy density of those oscillators with a given characteristic relaxation time τ (ǫ) (B8). In what follows we consider the time evolution of w(ǫ, t). The computation is outlined in Appendix A and the result (A13) is: In equilibrium one-time quantities are stationary. So the derivatives of both the energy and w(ǫ, t) are zero. From (14) and imposing stationarity for the energy, we obtain the relation (21) between the temperature and Q eq 2 . From (30) and imposing that the derivative of the w(ǫ, t) is zero, we obtain the following relation: valid for any ǫ in equilibrium. Eq. (31) is the equipartition theorem for each energy mode ǫ. It states that the average energy per degree of freedom is equal to T 2 . Next we consider the mode-density relaxation at T = 0. Eq. (30) can be solved in the asymptotic long-time limit using the results from the adiabatic approximation (25,26) by noting that ∂E ∂t << A(t) and ∂w(ǫ,t) ∂t → 0 yielding: which can also be written using the adiabatic approximation (22) as: This last result, which is valid for all values of ǫ in the asymptotic (adiabatic) long time limit, is based on the fact that the energy is a monotonic decreasing function of time. On the other hand, for low enough values of ǫ the first term in the rhs of (30) can be neglected and we get: We conclude that, for long-enough times and for continuous distributions of g(ǫ) with support on ǫ-values extending down to ǫ = 0, there must exist a value ǫ * such that ∂w(ǫ * ,t) ∂t = 0, i.e. ∂w(ǫ,t) ∂t changes sign at ǫ = ǫ * . From (30,33,34) we can distinguish three different regions of ǫ values depending on the sign of ∂w(ǫ,t) ∂t : Threshold Boundary : Equation (37) defines the time-dependent threshold ǫ * . Modes above the threshold are described by (35) and have a negative time-derivative of the mode-energy density w(ǫ, t) while modes below the threshold are described by (36) and have a positive time-derivative of the corresponding mode-energy density. For distributions g(ǫ) whose support of ǫ values does not extend down to zero it does not necessarily exist a value ǫ * which satisfies (37). In such cases the two regions (35,36) may still exist but not the threshold boundary. At zero temperature the system is always out-of-equilibrium: the acceptance decays as A(t) ≈ 1 t log t (25) and the energy decays as E(t) ≈ 1 log t (26), (these results being valid at leading order in 1/ log t). Therefore in the long-time limit: The asymptotic time-dependence of the threshold (37) is then given by the simple expression: From the dynamical equations (30) we can check that modes with ǫ > ǫ * verify the equipartition theorem at all times (32), so these modes are partially equilibrated at the temperature Θ(t) (27), i.e they have equilibrated over the hyper surface (in phase space) of constant energy E(t) = Θ(t) 2 . From (22), the value of Q 2 (t) decreases as the energy does. For distributions g(ǫ) with support on ǫ-values extending down to ǫ = 0, the threshold value ǫ * decays asymptotically to zero as a function of time. The time decay depends in a non trivial way on the form of the distribution g(ǫ) for ǫ → 0. For a distribution with finite weight at zero ǫ the decay seems to be a power law (see below). Unfortunately we have been unable to find a general analytical expression for such decay.
By numerically solving (39) we can compute the time-evolution (in Monte Carlo steps) for the energy threshold. In Fig. 3 we show the time evolution of ǫ * for a system that consists of 100 types of oscillators with g(ǫ) = 1 and spring constants uniformly distributed in the range ǫ = 0.01 − 1. It is clear that, as time goes on, the threshold energy diminishes and more and more oscillators partially equilibrate. As the system under consideration has 100 different types of oscillators then the solutions of (39) are discretized also. This is the reason why the time evolution for ǫ * shows a stepwise shape. As the energy of the mode gives the characteristic time-scale of the corresponding oscillators (see Sec.IV B below and Appendix B, in particular (B8)), those oscillators which have a higher value of the spring constant Kǫ also have a shorter relaxational time-scale. As time goes on, the modes with ǫ > ǫ * are the first to partially equilibrate, while those modes with lower ǫ equilibrate over longer time scales as they are cached up by the decreasing threshold value ǫ * (39).  ǫ2)) with ǫ1 = 1 and ǫ2 = 0.1. We also plot w(ǫ, t) coming from the corresponding equipartition relations (40) for the two modes of the system. The two continuous curves in the top correspond to ǫ1 while the two curves at the bottom correspond to ǫ2. The symbols correspond to the Monte Carlo simulation for a system with 5000 oscillators in each mode. The continuous lines correspond to the numerical solution of (30) and (40) corresponding to the exact and the adiabatic solution respectively. Different regions I (35) and II (36) according to the sign of the time derivative of w(ǫ, t) are also shown.
A summary of this behavior is shown in Fig. 4 for the case of a binary distribution of the disorder g(ǫ) = (1/2)(δ(ǫ − ǫ 1 ) + δ(ǫ − ǫ 2 )). For this particular example we note that equation (39) is only defined for ǫ = ǫ 1 , ǫ 2 so it does not have a well defined threshold solution ǫ * . However, as we have already mentioned earlier, we can distinguish the two regions (35) and (36) corresponding to different signs of ∂w(ǫ,t) ∂t . This is shown in Fig.4 where we show the time evolution of w(ǫ, t) at zero temperature (30). We also plot the solution for w(ǫ, t) using the equipartition theorem (40) for each of the two modes. We can clearly see that the oscillators with the largest value of ǫ partially equilibrate very fast to the temperature Θ(t) (27) (the two curves at the bottom), while the two upper curves deviate from each other even at long times. In Fig.4 we also plot the results obtained from a Monte Carlo simulation (symbols) and compare them with the numerical solution of the exact equation (30)  To sum up, in the non-equilibrium regime we have found two different classes of oscillators in the system: on the one hand there are partially equilibrated oscillators (ǫ > ǫ * ) for which the equipartition theorem (40) holds out of equilibrium as in the model without disorder 32 ; on the other hand, there are non-equilibrated modes (ǫ < ǫ * ) for which the equipartition theorem does not hold. This fact has been illustrated with the binary case of a system with two different values of ǫ, i.e two different time-scales.
In the next section we will define an effective temperature for each energy mode (each time-scale) in the system. Again we will distinguish between the partially equilibrated modes and the rest.

B. MODE CORRELATIONS AND RESPONSES
We can define the correlation and the response functions for each energy mode. Let us define the mode-correlation C(ǫ, t, s) as: which is the correlation between the positions of the oscillators which have the same spring constant Kǫ. In this case, similar calculations as in 32 , show that the dynamical evolution of the mode-correlation function is given by: We can compute the dynamic evolution for the mode-response after perturbing the Hamiltonian: where H 0 is defined in (5) and the perturbation is coupled to the magnetization of that particular mode of energy ǫ.
We define the mode-magnetization as: The mode-response function G(ǫ, t, s) is defined as usual, Computing the dynamical equation for the mode-magnetization and deriving it respect to the external field h, we obtain the dynamical equation for the mode-response function: Note that the evolution of the mode-correlation and the mode-response are identical and differ only from a delta function term. Equations (42,46) can be solved exactly and allow us to make a deeper analysis of the non-equilibrium behavior of the system. The equilibrium behavior of these functions is analyzed in the Appendix B.
In the general nonequilibrium regime the correlation function and the response function can be written (B2,B11): where F (ǫ, s) is given in (B3) and we have defined a time-dependent relaxation time τ (ǫ, t, s) (B4): The relaxation time of a given oscillator is inversely proportional to the value of its spring constant ǫ. This result indicates that each type of oscillator (each ǫ) can be considered as a single variable (evolving according Langevin dynamics) with a time-dependent relaxation time τ (ǫ, t, s) given by the expression (49). Note that the distribution of disorder fully determines correlations and responses through the functions f (x) and Q 2 (x). The full correlation function is then given by, (an analogous expression is valid for the response function). The relaxation of the full correlation function depends strongly on the disorder distribution g(ǫ). In equilibrium, the full correlation function can be written using (B7,B8,B9), Defining the parameter τ = (t−s)K γ this correlation (51) is given by the following Laplace transform, where we have defined h(ǫ) as: The equilibrium correlation function C eq (τ ) is well defined whenever the integral (52) converges, i.e. for g(ǫ → 0) ∼ ǫ a with a > 0. In this case, also the entropy constant in (7) is finite. When the integral (52) does not converge it means that the largest contribution to the full correlation function is determined by diffusive behavior of the ǫ → 0 modes. A relevant feature of the DOM is that it can describe the heterogeneous behavior observed in glass formers (see the final discussion in Section VII). The decorrelation time in the DOM can be defined in terms of the integrated relaxation time 16 : where we have used (20) (with ǫ = ǫ 1 ) and the expression for the viscosity γ (B10). This expression shows that the relaxation time is thermally activated and that the main effect of the disorder on the relaxation time is through the activation barrier which is proportional to ǫ. The other moments of the disorder distribution ǫ −1 , ǫ −2 only affect the prefactor of the relaxation time.
In general, the distribution of disorder g(ǫ) can be chosen in order to reproduce any given relaxation function by computing the inverse Laplace transform of the full correlation C eq (τ ). As an example we consider a few cases: • Binary distribution. Consider a binary system with one half of the oscillators having ǫ = ǫ 1 and the other half having ǫ = ǫ 2 where ǫ 2 > ǫ 1 . In this case, the full correlation function in equilibrium (51) is simply the sum of two exponentials, For long times, the full correlation (55) is dominated by the first exponential. For a general uniform discrete distribution with n finite different values of ǫ, the relaxation of the full correlation will be the sum of n exponentials of the type (B7).
• Continuous distribution. Consider now the continuous distribution g(ǫ) = Aǫ α exp(−ǫ) where A is the normalization factor. In this case, the full correlation function in equilibrium can be computed from (50,51) giving a power law decay: • Stretched exponential relaxation. Let us consider now the inverse problem of finding the disorder distribution for a correlation of the type: Taking the inverse Laplace transform, we get for the disorder distribution: where A is a positive constant. Note that this distribution is non-normalizable since it decays too slowly in the ǫ → ∞ limit. This is not a problem since the large ǫ limit only influences the short (t − s) behavior of the correlation function. Therefore by introducing a cutoff for g(ǫ) in the large ǫ limit we do not modify the stretching of the relaxation function. For the large-t generic stretching exponential behavior C eq (τ ) ∼ exp −Aτ b we have, in the ǫ → 0 limit, g(ǫ) ∼ exp(−B/ǫ b 1−b ) with B a positive constant. Taking b = 1/2 we recover the result (58) in the ǫ → 0 limit. The exponential relaxation corresponds to the case b → 1 and the g(ǫ) has a gap. In general, the large t behavior of C eq (τ ) is determined by the ǫ → 0 dependence of g(ǫ).
From this analysis we can see that the DOM is a very general model. The distribution of the disorder can be chosen in order to fit any experimental data just computing the inverse Laplace transform of the correlation.

C. MODE EFFECTIVE TEMPERATURES
Building up on the ideas by Cugliandolo, Kurchan and Peliti 9 we can define an effective temperature for each mode and analyze how the different degrees of freedom in the system relax. From (1) we can define an effective temperature for each mode: The mode effective temperature T eff (ǫ, s) only depends on the lowest or waiting time s but not on the largest time t. This is a typical feature of mean-field structural glass models 18 . As mentioned in Sec I, we can also define an effective temperature in the Fourier space (2) for each mode: where S(ǫ, Ω, s) and χ(ǫ, Ω, s) are defined as 16 : In equilibrium the power spectrum can be written as: where C eq (ǫ, t, s) is given in (B7). The result is a Lorentzian function: The full power spectrum is expressed in terms of the full correlation function (50) and is given by the integral of the Lorentzian weighed by g(ǫ): As for the full correlation function the disorder distribution g(ǫ) can be chosen to fit the experimental power spectrum data by using the relation (65). In (62) we have used the susceptibility, which in terms of the response is given by, Note that there is a factor 2 missing in (60) as compared to the standard definition of the FDT in frequency space 16 . This is due to the fact that we are assuming that the integral (62) extends from s (when the measure of the relaxation function starts) to ∞ rather than from −∞ to ∞ as in the equilibrium case. The calculations are shown in Appendix D. With these definitions we obtain for the effective temperature in Fourier space (D6): Interestingly, the effective temperature in Fourier space does not depend on the frequency and only depends on the waiting time s and the mode ǫ considered. We can observe that the result in Fourier space (67) is equal to the first term of (59) in the time domain so the subleading corrections introduced by the second term in (59) are absent in T F eff (ǫ, s). Let us note that it is also possible to define a "full" effective temperature T eff (s) (either in the time or frequency domains) from the full correlation function (50) and the analogous expression for the response function. From the analysis carried out in Section IV A we have that ǫ * → 0 asymptotically. The full effective temperature is then determined by the contribution to the full correlation and full response from the modes in region I, ǫ > ǫ * . In Fig.5 we can see the time-evolution (in Monte Carlo sweeps) of the normalized full effective temperature in time domain for the system considered in Fig.3, with 100 types of oscillators with spring constants uniformly distributed in the range ǫ = 0.01 − 1. As time goes on, more and more oscillators partially equilibrate and the full effective temperature approaches asymptotically to Θ(s) (27).  (27)) for the same system considered in Fig.3.
In the glass literature it is sometimes useful to construct the so-called fluctuation-dissipation (FD) plots 40 . Generally, a FD plot is a plot of the integrated-response versus the correlation function. In the DOM, the full correlation (50) at equal times C(s, s) depends on s. Therefore, for a meaningful FD plot we should consider the normalized full correlation C(t, s) = C(t, s)/C(s, s) and the full integrated-response given by: The FD plot consists on representing χ as a function of C for a fixed waiting time s and varying t. At equilibrium this plot gives a straight line with slope −1/T where T is the temperature of the bath. In the DOM, for asymptotic waiting times (s → ∞) we find, for a fixed s, Therefore, for a given asymptotic waiting time s the FD plot gives a straight line with negative slope equal to the inverse of the asymptotic effective temperature (27). In Fig.6 we represent the resulting asymptotic FD plot which has also been obtained before in the OM model 36 as a particular case. Straight FD plots are characteristic of the one-step behavior of structural glasses. In Fig.6 we can also see the absence of an equilibrium relaxation (usually called β-relaxation) for small responses. This unrealistic feature is due to the mean-field nature of the model. Next we consider in detail the effective temperature for the different types of modes. As we will see all modes in region I give the same mode effective temperature value which is equal to the value of the full effective temperature.
1. Region I of partially equilibrated modes: ǫ > ǫ * For the partially equilibrated modes, i.e ǫ > ǫ * , we just keep the first term in (59) (the derivative of the w(ǫ, s) goes to zero faster than w(ǫ, s) itself). Therefore, and the effective temperature reads: which coincides with the expression derived within the adiabatic approximation (27) (up to subleading logarithmic corrections of order 1/ log(t)). For the effective temperature in Fourier space we get the same result but without logarithmic corrections: All modes which are partially equilibrated share the same effective temperature that we define as Θ(s). As the time s goes on, more and more oscillators achieve this value Θ(s) and, in the asymptotic limit, all oscillators will thermalize to this effective temperature. In Fig.7 we can see the dependence of the effective temperature (59) on the energy mode at a given time. We plot the effective temperature normalized by Θ(s) for different values of s at zero temperature. As we increase the waiting time s those oscillators with higher value of ǫ (shorter relaxation time according to (B8)) tend to achieve an effective temperature T eff (ǫ, s) equal to Θ(s). Whereas those oscillators with lower values of ǫ have an effective temperature greater than Θ(s). From the figure we note that the ratio T eff (ǫ, s)/Θ(s) for ǫ > ǫ * is smaller than 1 due to the logarithmic correction coming from the second term in the r.h.s of (59) whereas such corrections are not present in T F eff (ǫ, s). Therefore, as the waiting time increases, lower energy modes progressively reach partial equilibrium. In the inset of Fig.7 we plot the normalized effective temperature as a function of ǫ/ǫ * showing good collapse. The shape of Fig.7 is expected to be quite general. The same qualitative picture was found for a Gaussian distribution for the disorder or an exponential one (data not shown).
Modes around the threshold, ǫ ≃ ǫ * , are not yet partially equilibrated. By plotting ∂w(ǫ,s) ∂s as a function of ǫ we can identify the width ∆ǫ of the region separating the partially equilibrated modes from the rest. This is shown in Fig.8 for the same system considered in Fig.7. There, we have multiplied ∂w(ǫ,s) ∂s by a factor s log(s) which is the inverse of the leading order of the decay of ∂w(ǫ,s) ∂s computed in the adiabatic approximation. In the inset we can clearly see that this quantity again scales as a function of ǫ ǫ * . From such scaling we note that the width ∆ǫ around the threshold scales as ǫ * and therefore, as ǫ * decays to 0, the region around the threshold becomes progressively thinner relative to the support of g(ǫ). This allows us to define the threshold boundary as the asymptotically thin sector of ǫ-values that separates region I from region II. The derivative of the mode-energy density multiplied by s log(s) as a function of ǫ for the same system as in Fig.7 In the inset we show the same plot as a function of ǫ/ǫ * .
2. Region II of diffusive modes: ǫ < ǫ * Using (30), the expression for the effective temperature (59) can also be written as: As a consequence those degrees of freedom with larger time-scales are hotter than partially equilibrated ones. This result can be checked in Fig.7 where the oscillators in the region ǫ < ǫ * have an effective temperature higher that Θ(s). For finite times, in the limit ǫ → 0, the second term in the rhs of (73) vanishes and the effective temperature is just twice the asymptotic effective temperature (the effective temperature for the partially equilibrated modes). Note that this value corresponds to the effective temperature of a diffusive process 41 , whose value is twice the temperature as given by the Einstein relation D = 2T where D is the diffusion constant associated to the diffusive variable x, (∆x 2 (t)) ∼ Dt. This result can be intuitively understood. Modes with ǫ small do not contribute much to the total energy (5), so their motion is completely uncorrelated from the rest of modes and become diffusive. In order to illustrate this result we show in Fig.9 the effective temperature for a binary system. In this case, 20% of the oscillators have ǫ = 0.1 while 80% have ǫ = 1. We can see how the oscillators with larger ǫ are partially equilibrated to the effective temperature Θ(s) (71) while the oscillators with smaller ǫ have an effective temperature approximately equal to 2Θ(s) as given by (73).  9. We plot the effective temperature for a binary system: 20% of oscillators have ǫ = 0.1 (triangles) while the rest have ǫ = 1 (circles). We also plot the time-evolution, measured in Monte Carlo steps, of Θ(s) and 2Θ(s) (71) (full lines) for the two modes.

V. FLUCTUATION RELATIONS AND THE EFFECTIVE TEMPERATURE
The effective temperatures for the partially equilibrated modes can be obtained by using a recently proposed fluctuation relation 42 . In 43 it has been suggested that effective temperatures can be derived by assuming microscopic reversibility between surfaces of constant energy. This relation has been also proposed in the framework of mean-field spin glass models 44 and applied to generalized oscillator models 36 . The effective temperature Θ(s) can be computed from the ratio between the probabilities of having a change in energy equal to δE and (−δE). From (11), where we have used (11). This relation between forward and reverse transitions only holds at the level of energies rather than configurations. It should not be confused with the detailed balance or reversibility property (4) of microscopic dynamics which obviously holds at the level of configurations in and out of equilibrium (4). Equation (75) defines a nonequilibrium temperature that is equal to the effective temperature obtained from the fluctuation-dissipation relation for the partially equilibrated modes (67,71). We can also extend this argument and compute the probability P s (Q ǫ ) of having an instantaneous change in the energy Q ǫ of a given mode ǫ at a given time s, The probability P s (δQ ǫ ) is computed in (A9). It is given by: and the fluctuation relation for each mode is given by: Only for partially equilibrated modes, we can use (71) and get, This relation emphasizes the fact that only for the partially equilibrated modes a unique effective temperature can be properly defined. Otherwise, for non-equilibrated modes 2ǫw(ǫ,s)

g(ǫ)
= Θ(s) Recent work 36,43 has shown that this effective temperature can also be obtained from microcanonical relations relating observable changes. In the present case, we consider the mode-magnetization (44). The probability of having a change δM ǫ in the mode-magnetization is computed in the Appendix C (C6). The rate of having a positive and a negative change gives: At difference with (78) now the effective temperature Θ(s) describes mode-magnetization fluctuations for all values of ǫ. This is not unexpected due to the linear dependence of the observable M ǫ on the configurational variables (44) and the Gaussian character of the distribution (C1) which determine that (C6) depends on Q 2 (i.e. on Θ(s)) for all values of ǫ. Although not proven we suspect that this general dependence might not hold in general, e.g. for generic non harmonic potentials 36 .

VI. SUMMARY OF RESULTS
In this paper we have introduced an exactly solvable disordered model for glassy dynamics, the disordered oscillator model (DOM). The model is the disordered version of the oscillator model (OM) previously introduced by one of us 32 as a simple model with entropic barriers for aging dynamics. The DOM consists of a collection of non-interacting oscillators or modes with different spring constants. The model has trivial statics but is kinetically constrained and has cooperative dynamics displaying glassy behavior at low temperatures. In the DOM the relaxation time associated to each mode is inversely proportional to the effective spring constant. This allows us to study the relaxational properties (both in equilibrium and in the aging regime) associated to each mode and the value of the effective temperature (1) for the different degrees of freedom.
Dynamics in the DOM can be solved writing down a hierarchy of coupled first-order differential equations (19). These set of equations can be closed by introducing a generating function and using the method of characteristics or the adiabatic approximation at low temperatures. The latter approach is much simpler and allows us to define an effective temperature Θ(s) (27) in the aging regime (s being the waiting time). The relaxation of one-time quantities appears similar to the OM, just replacing the spring constant in the OM, K, by the average spring constant K ǫ in the DOM.
To investigate how different energy modes relax to equilibrium we introduce the mode-energy density (28), the mode-magnetization (44) and the mode-correlations and mode-responses (41,45). From the dynamical equations equations (30) we can prove the existence of time-dependent threshold located at ǫ * (s) separating two regions (I and II) characterized by the sign of the time derivative of the mode-energy density (35,36,37). Only the modes in region I satisfy the equipartition relation at the effective temperature Θ(s) (40). The threshold value ǫ * (s) decays with the age of the system with more modes entering region I as the system approaches equilibrium. In Fig.10 we show an schematic picture of the different mode regions in the model. Furthermore, for each mode we can define an equilibrium relaxation time (B8) in terms of a temperature dependent viscosity that diverges when T → 0 (B10). Full correlation functions diverge if too many diffusive modes are present, i.e. if g(ǫ) vanishes too slowly in the limit ǫ → 0. For any equilibrium correlation function C eq (τ ) there exists a distribution g(ǫ) that in the limit ǫ → 0 reproduces the long-time decay of C eq (τ ).
Mode-dependent effective temperatures can be also defined (59) from mode-correlations and mode-responses. They agree with the previous scenario. In region I modes are mutually equilibrated at the temperature Θ(s) whereas in region II modes are out-of-equilibrium. The latter show diffusive behavior and are characterized by an effective temperature larger than Θ(s) and equal to 2Θ(s) for ǫ → 0. This last result is characteristic of purely diffusive behavior (X ∞ = 1/2, see 16,18,41 ). The effective temperature can also be computed in frequency space (2) corroborating the results obtained in the time-domain. It is noteworthy to mention that no frequency dependence is found (67). This result is in contradiction to what is observed in experiments and might be consequence of the simplicity of the model (see the discussion in the next section).
Finally, we have corroborated the scenario of mode-dependent effective temperatures using the fluctuation relations recently proposed in the framework of aging systems 36,42,43 . Interestingly, for both partially equilibrated modes (region I) and diffusive modes (region II) the fluctuation relation for the mode-magnetization gives the unique value Θ(s). It is unclear to us whether this result is of general validity or restricted to the DOM. It would be interesting to put under scrutiny other types of models. Region I (ǫ > ǫ * ) consists of partially equilibrated modes at the temperature Θ(s). In region II (ǫ < ǫ * ) there are the non-equilibrated or diffusive modes. In between the two regions there is a boundary region of width ∆ǫ around the threshold ǫ * where modes start to equilibrate. As time goes on, the threshold value ǫ * decreases with the width ∆ǫ being proportional to ǫ * . For sake of clarity the size of the regions in the picture is not in scale (region II and ∆ǫ should be tiny as compared to region I).

VII. CONCLUDING REMARKS
In which ways can the DOM improve our current understanding of the glassy state? A lacking piece in the glass transition puzzle is to understand whether it is possible to identify effective parameters (such as local temperatures or stresses) that reflect the nonequilibrium relaxation toward equilibrium of macroscopic or bulk properties (e.g. enthalpy, specific heat, volume, pressure, viscosity, thermal conductivity,...). However, if dynamics and transport phenomena are strongly heterogeneous then it is hard to understand how a description of the macroscopic time evolution of bulk quantities can be achieved from a few number of parameters as proposed by several phenomonological approaches.
It would be very interesting to get experimental evidence that a unique nonequilibrium or effective temperature can be identified in glass formers in the aging state. Such temperature would constitute a candidate parameter to describe many physical nonequilibrium properties in the aging state. This point of view has been advocated a few years ago by Nieuwenhuizen 45 . The existence of a unique effective temperature appears to be an inherent property of mean-field systems. The DOM is yet another example of a mean-field system with a well defined nonequilibrium or effective temperature (that we have denoted by Θ(s)). Nevertheless, the current experimental evidence points toward the fact that dynamics is strongly heterogeneous and a unique nonequilibrium temperature might not exist. Rather, the effective temperature appears more as a local quantity that changes from region to region inside the glass and which takes a value that depends directly on the degree of mobility of the region.
The DOM is a mean-field model with cooperative glassy dynamics. Alone it cannot explain the heterogenous behavior observed in glasses because it does not incorporate the fact that different modes intermittently switch for mobile to frozen behavior. The DOM describes the dynamics of an isolated cooperatively rearranging region (CRR) maybe of a few nanometers of spatial extension that includes a few hundreds of atoms or molecules. The different energies ǫ of the vibrational modes in each CRR are determined by its volume (the number of atoms it contains) and its surface through possible interactions with nearby CRRs. The glass could then be viewed as a mosaic made out of many CRRs each one described by a DOM with a given energy levels distribution g(ǫ). Such distribution strongly fluctuates among CRRs. In the scenario where CRRs intermittently appear and dissapear at a given rate, each time a new CRR is produced it carries a new energy level distribution g(ǫ). Consequently, at a given time each region shows a different ǫ. According to (54), fluctuations in ǫ will lead to a great disparity of relaxation times among CRRs and to the presence of slow and fast regions. Also, the effective temperature will be inversely proportional to ǫ (71) leading to an effective-temperature age-dependent field Θ( ǫ, s) as well as a Ω-dependent effective temperature, T F eff (Ω, ǫ, s). The latter dependence of the effective temperature on the probing frequency would agree with the experimental results 11,31 . A schematic picture of this scenario is shown in Fig. 11. This microscopic scenario is similar to that proposed by Xia and Wolyness 22 from the random first order transition theory for mean-field spin glasses 15,46 applied to the glass transition. A related scenario in terms of aggregating activated regions has been also described by Crisanti and Ritort 47 . A more elaborated and compelling description of the mosaic picture has been recently proposed by Bouchaud and Biroli 23 .
Merging of random energy levels and cooperative dynamics in spatially localized regions (CRRs) appears a promising route to provide a phenomenological description of the glassy state. It remains to be seen whether such view is capable to provide new predictions that can be tested in the laboratory. This maybe difficult until we get a more clear picture and better understanding about the nature of the heterogeneities uncovered by the experimentalists. In turn, this might provide a definitive experimental evidence in favour of the intuitive notion of CRRs put forward by Adam-Gibbs-Di Marzio many years ago. 11. Schematic picture of a piece of glass made out of different CRRs, each one characterized by a value of ǫ (we show 9 labelled regions in the picture). In such a mesoscopic description regions are maybe of a few nanometers of spatial extension and typical activation energies K∆ 2 ǫ are of the order of a few kBT . ǫ values fluctuate from region to region. According to (54) dynamics is highly heterogenous with relaxation times in local regions differing by many orders of magnitude. The effective temperature (71) is also inversely proportional to ǫ leading to an effective-temperature field.

APPENDIX A: DYNAMICAL EQUATIONS FOR QK
To solve the hierarchy of equations for Q k we start by computing the evolution for the quantity Q 2 (t) (12) which has appeared in the dynamic evolution for the energy (14). In order to do that we compute the joint probability distribution of having a change of δE in the energy (5) and δQ 2 in Q 2 (12) for an elementary move (8), This probability factorizes and can be written as follows: where x = δE, y = δQ 2 and P (x) is given in (11). The other quantities can be defined through the general expressions: From the joint probability we can derive the equation for Q 2 (t) using (4): The integral over y is just the quantity J 1 and does not depend on J 2 . Expressed in terms of the dynamical evolution of the energy ∂E ∂t and the acceptance A we find: As we expect, the evolution of Q 2 depends also on the new quantity Q 3 . Therefore we have to generalize these computations in order to obtain the whole set of coupled equations. The result is: These equations are, for a given probability distribution of the disorder, the complete solution for the dynamical one-time quantities.
At zero temperature only negative changes in energy are accepted. The equations turn out to be a bit simpler: We now outline the computations for the mode-energy density. As for the quantities Q k we must first compute the joint probability of having a change of δE in the energy (5) and δw in w(ǫ) (28), the result is: where now x = δE, y = δw and P (x) is given in (11). The quantities S 1 (ǫ) and S 2 (ǫ) are given by, From the equation, we get: ∂w(ǫ, t) ∂t = ǫw(ǫ, t) Q 2 (t) ∂E ∂t + K∆ 2 g(ǫ) 2 − a c ǫw(ǫ, t) Q 2 (t) A(t) . (A13)

APPENDIX B: EQUILIBRIUM BEHAVIOR OF THE CORRELATIONS AND RESPONSES
In this appendix we will solve the dynamical equations for the mode-correlation and the mode-response functions and analyze their equilibrium behavior. First of all, let us define: In terms of this function the solution for the mode-correlation (42) is given by: Where the initial condition for the mode-correlation is: C(ǫ, s, s) = 2 K w(ǫ, s). In a more compact way we can write: In general we can define a time-dependent relaxation time, and the correlation function can be written as: In equilibrium Q eq 2 is given by (21) and f eq (B1) is given by: f eq = a c erfc(α eq ) , where α eq is given in (15). Therefore, the mode-correlation in equilibrium has an exponential behavior: C eq (ǫ, t, s) = 2w eq K exp [−(t − s)/τ eq (ǫ)] , where w eq is given in (31) and the characteristic equilibration time of the energy mode ǫ is: τ eq (ǫ) = 1 Kǫ 2T ∆ 2 erfc(α eq ) .
As we have mentioned before the characteristic relaxation time of all the oscillators belonging with energy mode ǫ is inversely proportional to the value of ǫ. Therefore, the higher is the value of ǫ the shorter is the characteristic time-scale of the oscillator. Therefore we can define the friction coefficient γ using the relation τ eq (ǫ) = γ Kǫ for an overdamped harmonic oscillator in thermal equilibrium: which, for low enough temperatures behaves, Therefore, the divergence of the relaxation time is of an activated type just as in the model without disorder. We can also derive the equilibrium relaxation time from the response function: with F (ǫ, x) given in (B3). In equilibrium the mode-response becomes: G eq (ǫ, t, s) = g(ǫ) ǫτ eq (ǫ)K exp [−(t − s)/τ eq (ǫ)] θ(t − s) .

APPENDIX C: MODE-MAGNETIZATION RATES
In this section we derive the main results for the mode-magnetization rates. The joint probability of having a change in the mode magnetization (44) δM ǫ and in the energy δE can be computed as in the original model 34,43 . The result is: where we have defined the quantities: S 3 (ǫ) = ǫM K∆ 2 − h∆ 2 g(ǫ) , (C3) At T = 0 the probability that an attempted change δM ǫ is accepted is given by: Assuming that the value of the magnetization is zero (condition of neutrality 43,48 ) we find: where and therefore,