The disordered Backgammon model

In this paper we consider an exactly solvable model which displays glassy behavior at zero temperature due to entropic barriers. The new ingredient of the model is the existence of different energy scales or modes associated to different relaxational time-scales. Low-temperature relaxation takes place by partial equilibration of successive lower energy modes. An adiabatic scaling solution, defined in terms of a threshold energy scale $\eps^*$, is proposed. For such a solution, modes with energy $\eps\gg\eps^*$ are equilibrated at the bath temperature, modes with $\eps\ll\eps^*$ remain out of equilibrium and relaxation occurs in the neighborhood of the threshold $\eps\sim \eps^*$. The model is presented as a toy example to investigate conditions related to the existence of an effective temperature in glassy systems and its possible dependence on the energy sector probed by the corresponding observable.


I. INTRODUCTION
The study of exactly solvable models has been always an active area of research in the field of statistical physics. They help us to grasp general principles governing the physical behavior of realistic systems which, due to the complicated interactions among the different constituents, cannot be predicted using standard perturbative techniques. Glasses in general are systems falling into this category. The slow relaxation of glasses observed in the laboratory is a consequence of the simultaneous interplay of its constituents that yields a very complex and rich phenomenology.
It is well known that glasses fall out of equilibrium when the characteristic observation time is larger than their relaxation time. Because the relaxation time is strongly dependent on temperature, it turns out that glasses are immediately out of equilibrium as soon as the temperature is few degrees below the glass transition. Well below the glass transition temperature T g no time evolution is apparently observed in the glass and one is tempted to conclude that the glass is in a stationary state. Nothing more far from the truth. Glasses still relax but slow enough for any change to be observable in laboratory time-scales. Old experiments on polymers reveal that the slowly relaxing state corresponds to an aging state 1 . That is, if the system is perturbed while being in its aging state, then the characteristic time associated to the response of the system scales with the age of the system (i.e. the time elapsed since it was quenched). Another way to look at this aging phenomena is to evaluate the time autocorrelation function. It is observed that the typical decorrelation time scales with the age of the system 2 .
A simple scenario to explain these results is the following. Consider a liquid well above T g where correlations decay exponentially with time. One may consider the resultant behavior of the liquid as the superposition of different and independent harmonic modes. Each of this energy modes corresponds to a normal mode of a system and describes a collective oscillation of N atoms around their local minimum. This is the harmonic approximation which is known to work quite well in liquids. Nevertheless, already as the temperature goes below a critical temperature T c (the transition temperature predicted by the mode-coupling theory 3 ) other collective modes different from the standard vibrational ones become important. The nature of these modes is quite different from the usual harmonic normal modes because they do not represent oscillations around a given configuration within a metastable well but transitions among different wells. These modes are reminiscent of some type of instanton solutions recently computed in the framework of some spin glass models 4 . Below T c relaxational dynamics proceeds by activation over the barriers characterizing these collective modes. Now, the main difference between these collective modes and the usual harmonic modes relies on how they relax to equilibrium when put in contact with a thermal bath at temperature T . Relaxation to equilibrium is determined by the height of the energy barriers separating different modes. Suppose a given normal mode has frequency ω k and energy E k ∝ ω 2 k . The relaxation time for each of these modes is typically of order τ k ∼ exp(E k /K B T ). Therefore, as the energy E k is lower than the thermal bath temperature this mode rapidly equilibrates. On the contrary, if E k ≫ K B T , this mode remains frozen. Collective modes are different. As the reference energy of the collective modes deplete, the typical barrier separating these modes increases leading to the contrary behaviour and to super-activation effects. While high energy collective modes are separated by low barriers, low energy collective modes are separated by high barriers. A simple schematic representation of this scenario in a one dimensional configurational space is shown in figure 1. This behavior is common to the majority of exactly solvable glassy models [5][6][7] , phenomenological trap models 8 and kinetically constrained models 9 to cite a few. In what follows we will use the generic word mode to refer to this kind of collective excitation. Let us label the modes with the integer variable r and let us denote their energy by ǫ r . Let us suppose that the energy levels are ordered from lower to higher energies according to the label r. It is natural to assume that there is a characteristic mode r * with associated energy ǫ * such that, all modes with ǫ r ≫ ǫ * have already relaxed, while in the other limit, ǫ r ≪ ǫ * , all modes are frozen. If the system is quenched well below T c then equilibration cannot be achieved in laboratory time scales, this means that all modes below ǫ * remain frozen while modes above ǫ * remain equilibrated at the bath temperature. The energy threshold ǫ * decays with time because, as time goes by, higher barriers are accessible to the system. The resulting scenario is that of a liquid where collective modes above ǫ * are in some sort of local equilibrium at the temperature of the bath T while modes below ǫ * are frozen. This scenario, as it stands, is too naive because it is based on the assumption that there are no dynamical correlations between the different modes, i.e. high energy modes do not influence low energy modes. While this is true in the equilibrium state it may not be valid (an indeed it is not) when any type of local dynamics induces correlations between different modes. Still it is interesting to investigate in which conditions such a description turns out to be correct. Examples where these type of description holds are mean-field models 5,10-13 . By mean-field we mean those models where there is no spatial dimensionality associated with the set of interactions. In these cases, different modes are not related to different length scales. Therefore local dynamical rules do not necessarily induce correlations between the different modes. Generally speaking, the identification of this energy scale remains an open problem for which we do not have yet a complete understanding.
This energy scale ǫ * is related to what has received the name of effective or fictive temperature T eff in the most recent literature about glassy dynamics [14][15][16][17]19,20 . T eff is an effective time-dependent parameter describing equilibrium fluctuations for those thermalized modes with ǫ ≫ ǫ * = k B T eff . While a thermometer coupled to the fast modes is expected to measure the bath temperature, a thermometer coupled to the slower ones should yield the effective temperature. Hence, the effective temperature corresponds to an energy threshold which separates collective modes which are frozen from those which have relaxed to the bath temperature. Relaxation, then, takes place at energy scales in the neighborhood of the threshold value ǫ * . For real systems this energy threshold is related to the typical volume of the drop which is able to release strain energy during its relaxation to equilibrium. A scenario on this type of physical mechanism for glassy relaxation has been recently introduced in 21 where this threshold energy ǫ * has been related to the size of the cooperative region as an explanation of the fragility and super-activation anomalies in real glasses.
The existence on an effective temperature is tightly related to the validity of some approximations used in the context of slowly relaxing systems such as the adiabatic approximation. Within the adiabatic approximation one obtains a Markovian description for the dynamics. Again, this Markovian description encodes within a single parameter (the effective temperature) all the complicated previous past of the system. The adiabatic approximation has been shown to give the correct asymptotic dynamical behavior in some simple models of glasses such as the Backgammon model 22,23 while for the majority of the other most famous models in the literature (for instance the p-spin model 5 ) the implementation of such approximation is generally not yet known. Quite probably a Markovian description for glassy dynamics is unrealistic and the original idea of experimentalists 24 to encode the dynamical behavior into effective parameters such as the effective pressure or the effective (else said fictive) temperature could be only labels without a deep physical meaning.
These ideas have been contrasted in several exactly solvable models. All these models have the advantage of being mean-field, hence dynamical equations can be closed in one or another way. Still, the mean-field character of these models does not generally allow to investigate in a simple way the relaxational properties of the different modes of the spectrum. In general, the dynamics of all these models has been studied by addressing the closure of the dynamics associated to global quantities (such as the energy or correlation-response functions) which do not discern the contribution to the global relaxation of the different energy modes. Examples of these models are the p-spin model 5,10,13,16 , the Backgammon 20,22,23,25 the harmonic oscillator and other spherical spin models 7,16,17 . The question of the existence of an effective temperature is also tightly related to the particular way in which the fluctuation-dissipation theorem (FDT) is violated 13,26,27 . The question whether there exists a single effective temperature describing the violation of FDT is still controversial (see, for instance, 28 ). We believe that this and other related questions can be better addressed, even at the level of mean-field models, by analyzing the contribution of all different modes. The analysis of the harmonic vibrational modes in glasses is already an interesting step in this direction 29 although here we have in mind other type of collective modes.
The purpose of this paper is the study of a new model where the dynamical relaxation of the different energy modes can be made explicitly clear. This model is what we refer to as the disordered Backgammon model (DBG model) and consists in a generalization of the Backgammon model to allow different energies for different boxes. Again, just like its predecessor, the slow relaxation of this model is due to entropic barriers. For the DBG model we show the existence of an energy threshold ǫ * (t) which separates equilibrated form non-equilibrated modes. The DBG model (with its associated threshold ǫ * (t)) provides a microscopic realization that is reminiscent of some phenomenological models proposed in the past such as the trap model in a tree considered by Bouchaud and Dean 30 . The advantage in the DBG model is that now one can exhaustively investigate the distinct relaxation of each of the different energy modes verifying whether the scenario of the effective temperature presented before holds. This will help to better understand the meaning of the effective temperature and its relation with the violation fluctuation-dissipation theorem.
In section II we introduce the disordered Backgammon (DBG) model. Its thermodynamic properties are reported in section III, where the special case of very low temperatures is explicitly worked out. In section IV we present the dynamical equations whose solutions are found within a new type of adiabatic approximation in section V. There the features of such an approximation are carefully analyzed. In section VI numerical results are presented for two specific models belonging to the family of DBG models . Finally, in section VII, a novel method is introduced in order to estimate the threshold energy scale ǫ * directly from the dynamics. In section VIII we present our conclusions. Some technical issues are presented in four different appendices.

II. THE DISORDERED BACKGAMMON MODEL
A. Definition of the model Let us take N particles which can occupy N boxes, each one labeled by an index r which runs from 1 to N . Suppose now that all particles are distributed among the boxes. A given box r contributes to the Hamiltonian with an energy −ǫ r only when it is empty. In this case the total Hamiltonian of the system reads where δ is the Kronecker delta and n r denotes the occupancy or number of particles in box r. The ǫ r are quenched random variables extracted from a distribution g(ǫ) which we assume to be defined only for ǫ ≥ 0. The interest of this definition will be discussed below. Like in the original Backgammon (BG) model 6 we consider Monte Carlo mean-field dynamics where a particle is randomly chosen in a departure box and a move to an arrival box a is proposed. If d denotes the departure box, the proposed change is accepted according to the Metropolis rule with probability W (∆E) = Min(1, exp(−β∆E)) where ∆E = ǫ a δ na,0 − ǫ d δ n d ,1 . Note that the departure box satisfies n d ≥ 1 and departure boxes are chosen with probability n d /N . This dynamics corresponds to Maxwell statistics where particle are distinguishable and differs from the one corresponding to Bose statistics 31,23,32,20 where particles are indistinguishable and arrival boxes are chosen with uniform probability 1/N . In the dynamics the total number of particles is conserved so that the occupancies satisfy the closure condition N r=1 n r = N . (2) Now we want to show that the interesting case corresponds to the situation where g(ǫ) is only defined for ǫ ≥ 0. In this case, the dynamics turns out to be extremely slow at low temperatures, similarly to what happens for the original BG model. The difference lies in the type of ground state. The ground state of (1) corresponds to the case where all particles occupy a single box, the one with the smallest value of ǫ. Let us denote by ǫ 0 this smallest value. Then the ground state energy is given by the relation E GS = − N r=1 ǫ r + ǫ 0 . Since all the ǫ are positive no other configuration can have a lower energy. If g(ǫ) is a continuous distribution the ground state is also unique. Now it is easy to understand that, during the dynamical evolution at zero temperature, all boxes with high values of ǫ become empty quite soon and the dynamics involves boxes with progressively lower values of ǫ. The asymptotic dynamics is then determined by the behavior of the distribution g(ǫ) in the limit ǫ → 0. If g(ǫ) ∼ ǫ α , for ǫ → 0, we will show that the asymptotic long-time properties only depend on α. Note that the normalization of the g(ǫ) imposes α > −1. This classification includes also the original BG model where there is no disorder at all. In that case g(ǫ) = δ(ǫ − 1) so the distribution has a finite gap at ǫ = 0. The behavior corresponding to this singular energy distribution can be obtained from the previous one in the limiting case α → ∞.
One important aspect of the model (1) is that, in the presence of disorder, it is not invariant under an arbitrary constant shift of the energy levels. Actually, by changing ǫ r → ǫ ′ r = ǫ r + c with c ≥ 0, the model turns out to be a combination of the original model plus the classical BG model. After shifting, the new distribution g(ǫ ′ − c) has a finite gap (equal to c plus the gap of the original distribution). The new model corresponds again to the α → ∞ case and the asymptotic dynamical behavior coincides with that of the standard BG model. As we will see later the present model is characterized by an energy threshold ǫ * which drives relaxation to the stationary state. Only when the energy threshold can go to zero we have a different asymptotic behavior. For all models with a finite gap, ǫ * cannot be smaller than the gap, hence asymptotically sticks to the gap and the relaxational behavior of the DBG model with a finite gap corresponds to that of the standard BG model.
One of the outstanding features of this model is that a description of the dynamics in the framework of an adiabatic approximation turns out to be totally independent from the type of distribution g(ǫ) (and hence on α) despite of the fact that the asymptotic long-time behavior of the effective temperature and of the internal energy depend on the value of α.

B. Observables
Like in the original BG model we define the occupation probabilities, P k , that a box contains k particle, and the corresponding densities that a box of energy ǫ contains k particles The P k and the g k are related by and the conservation of particles reads The energy can be expressed in terms of the density g 0 (ǫ) as This set of observables depend on time through the time evolution of the occupancies n r of all boxes. In the next section we analyze the main equilibrium properties of the model.

III. EQUILIBRIUM BEHAVIOR
The solution of the thermodynamics proceeds similarly as for the case of the original BG model. The partition function can be computed in the grand partition ensemble. It reads where z = exp(βµ) is the fugacity, µ is the chemical potential and Z C (N ) stands for the canonical partition function of a system with N particles. The canonical partition function can be written as where δ i,j is the Kronecker delta. Introducing this expression in (9) we can write down Z GC as an unrestricted sum for all the occupancies n r The factor N !/ r n r ! in the partition function Z C is introduced to account for the distinguishability of particles. This factor leads to an over-extensive entropy (i.e., the Gibbs paradox) which can be cured eliminating from Z GC the over-counting term N ! in the numerator. The final result is yielding the grand-canonical potential energy per box where F is the Helmholtz free-energy per box. The fugacity z is determined by the conservation condition (2) which reads yielding the closure condition This equation gives the fugacity z as function of β and, from (13) and its derivatives, the whole thermodynamics. In particular, the equilibrium expressions for g k (ǫ) are the corresponding P k being given by eq.(6) which, together with the closure relation (15), leads to the expression Starting from eq.(8) the equilibrium energy density is obtained as All these expressions can be evaluated at finite temperature. Note that, although the values of P eq k in (17) are independent of the disorder distribution g(ǫ), they directly depend on that distribution through the equilibrium value of z (which obviously depends on the g(ǫ)). Of particular interest for the dynamical behavior of the model are the low-temperature properties that we analyze below.

A. Thermodynamics at low temperatures
A perturbative expansion can be carried out close to T → 0 to find the leading behavior of different thermodynamic quantities. Let us start analyzing the closure condition (15). Doing the transformation ǫ ′ = βǫ, eq.(15) can be rewritten as In the limit T → 0 the fugacity z depends on the behavior of g(T ǫ ′ ) in the limit T = 0, i.e. on the behavior of g(ǫ) The integral can be expanded around T = 0 by taking successive derivatives of the function h: Using the asymptotic result z → ∞, when T → 0, everything reduces to estimate the following integral in the large z limit: The term k = 0 in the series yields the leading behavior for z, which turns out to be In a similar way the energy can be computed to leading order in T , giving a finite specific heat at low temperatures. In section VI we will show explicitly such a behaviour for two specific DBG models, one with α = 1 (defined in eq. (48), sec. VI) and the other with α = 0 (eq. (49), sec. VI). Solving eq. (15) numerically for z(T ) in each specific model and inserting z(T ) in the expression (18) for the equilibrium energy density, we get the energy dependence on the temperature. In all cases, at equilibrium, the energy is linear for low T , as predicted in eq. (23). It yields, therefore, a finite specific heat as corresponds to a classical model with Maxwell-Boltzmann statistics (see fig. (2)). Left: energy dependence on temperature at equilibrium for two DBG models, with α = 1 and α = 0. Also the energy for the standard BG model (α → ∞) is plotted. Right: specific heat for the same models. For all cases (disordered or not) it turns out to be finite at T = 0.

IV. DYNAMICAL EQUATIONS
Here we consider the dynamical equations for the occupation probabilities P k and their associated densities g k (ǫ). The dynamical equations in this model are derived in a similar way as for the standard BG model. The main difference is that in the DBG the equations for the occupancies probabilities P k do not generate a closed hierarchy of equations. Only for T = 0 such a closed hierarchy is obtained. As we will see later, this has important consequences when we become interested in the zero-temperature relaxation.
A hierarchy of equations can only be obtained at the level of the occupation probability densities g k (ǫ). A detailed derivation of these equations is reported in the Appendix A. Here we show the final result, The equations for the P k are directly obtained by integrating the g k (ǫ) according to (6). The result is (27) with P −1 = 0. It is easy to check that the equilibrium solutions (16) are indeed stationary solutions. As previously said, for general β, the equations for the P k do not generate a hierarchy by themselves but depend on the g k (ǫ) through the distribution g 0 (ǫ) in (27). Nevertheless, a remarkable aspect is that they generate a well defined hierarchy at T = 0 which coincides with the equations of the original BG model 25 . These are It is easy to understand why at T = 0 the dynamical equations are independent on the density of states g(ǫ). The reasoning is as follows. For T = 0, all moves of particles between departure and arrival boxes with different energies ǫ d and ǫ a depend on the precise values of these energies only when the departure box contains a single particle and the arrival box is empty. But such a move does not lead to any change in any of the P k , hence dynamical equations for the P k remain independent of g(ǫ). Obviously this does not hold for other observables such as the energy E (i.e. the mean value of ǫ over the density g 0 (ǫ)) and higher moments of g k (ǫ).
This observation is crucial, since we look at the glassy behaviour at T = 0. It turns out that the analysis of the dynamical equations for the DBG at T = 0 decomposes into two parts. On the one hand, the equations for the P k coincide with those of the original non-disordered BG model. On the other hand, in order to analyze the behavior of the energy one must analyze the behavior of the hierarchy of equations for the g k (ǫ) which is quite complicated. In the next section we will see how the analysis of these equations can be done within the framework of a generalized adiabatic approximation.

V. THE ADIABATIC APPROXIMATION
A. Standard adiabatic solution for the P k In this section we are interested in the solution of the dynamical equations at T = 0. We already saw that the equations (28) for the P k coincide with those of the original non-disordered BG model. Consequently, the same adiabatic approximation used for the P k in the original BG model is still valid for the DBG. Let us remind the main results 22,23 . The key idea behind that approximation is that, while P 0 constitutes a slow mode, the other P k with k > 0 are fast modes. Hence they can be considered as if they were in equilibrium at the hyper-surface in phase space P 0 = constant, this constant being given by the actual value of P 0 at time t. In the original BG model P 0 = −E, hence thermalization of the fast modes P k (k > 0) occurs on the hyper-surface of constant energy. For the DBG model this is not true, the hyper-surface where equilibration of fast modes occurs does not coincide with the constant energy hyper-surface simply because the energy and P 0 are different quantities. Indeed, we will see later that their leading time behavior is different.
At T = 0 the equation for P 0 (eq. (28) for k = 0) reads If local equilibrium is reached on the hyper-surface of constant P 0 we can relate P 1 to P 0 using eqs. (17). The simplest way of dealing with (29) is to relate both P 1 and P 0 to the time-dependent fugacity z * writing down a dynamical equation for z * . Using eq.(17) this yields .
In the large time limit z * diverges and the leading asymptotic behavior is given by z * ≃ log(t) + log(log(t)). The occupation probabilities are, then, given by the relations (17) replacing z by the time-dependent fugacity z * , as long as z * ≫ 1, Hence, according to (22), the inverse effective temperature is given by the relation α+2 α+1 (32) and the effective temperature depends on the properties of the disorder distribution g(ǫ) in the limit ǫ → 0 through the value of the exponent α. Clearly, when the density of levels decreases as we approach ǫ = 0, the relaxation turns out to be slower; the limiting case being the original BG model for which α → ∞ and β eff ∼ log(t). In the other limit α → −1, when disorder becomes unnormalized, the inverse effective temperature diverges very fast. Already from (23) one can anticipate that the same asymptotic behavior holds for the energy (see eq. (44), hence α interpolates between fast relaxation (α = −1) and very slow relaxation (α = ∞). A relaxation slower than logarithmic is not possible in the present model. The equations for the g k (ǫ) at T = 0 are To solve the dynamical equations for the g k (ǫ) in the adiabatic approximation we note that, contrarily to the global quantities P k , they cannot be equilibrated among all different modes. The reason is that, due to the entropic character of the relaxation, very low energy modes are rarely involved because the time needed to empty one further box increases progressively as time goes by, hence they cannot be thought as effectively thermalized. Note that in the original BG model all boxes have the same energy, hence there is a unique class of modes. For the general disordered model we expect the existence of a time dependent energy scale ǫ * separating equilibrated and non-equilibrated modes. The mechanism of relaxation is the one we reminded in the introduction when speaking about the behavior of collective modes. At zero temperature there is no thermal activation and the equilibrated modes are in the sector ǫ ≫ ǫ * while the non-equilibrated modes are in the other sector ǫ ≪ ǫ * . The value of ǫ * can be easily guessed. After quenching to zero temperature the system starts to relax to its ground state. Because there is no thermal activation, relaxation is driven by entropic barriers, i.e. flat directions in configurational space through which the system diffuses. Entropic relaxation is energy costless so its rate is determined by the number of available configurations with energy smaller or equal to the actual energy. A simple microcanonical argument giving the relaxation rate goes as follows. Let us denote M the number of occupied boxes and Ω(M ) the number of configurations with M occupied boxes. The typical time to increase by one unity the number of empty boxes is given by, .
For distinguishable particles we have Ω(M ) ≃ M N . In the large N, M limit we get Using the relation P 0 = 1 − M/N and (36) we find where we have used the fact that at zero temperature the number of occupied boxes can only decrease by one unity, thus ∆M = −1. This yields the result P 0 ≃ 1 − 1 log(t)+log(log(t)) in agreement with the adiabatic approximation. Actually, from the solution in (30) for z * and using the adiabatic relation (31), we obtain P 0 ≃ 1 − 1/z * , which yields the same result. The typical relaxation time (37) behaves, then, like τ ≃ exp( 1 1−P0 ) ≃ exp(z * ). If the threshold ǫ * plays the role of an energy barrier and β eff accounts for the effective thermal activation due to entropic effects, we obtain, for the typical relaxation time, τ ≃ exp(β eff ǫ * ). This expression is only valid to leading order. As we will see below there are sub-leading corrections to this expression arising for the fact that the relaxation time is better described by the expression τ ≃ exp(β eff ǫ * )/(β eff ǫ * ) (see eq. (45)). Hence, at a given time-scale t (i.e., the time elapsed since the system was quenched) all modes where τ ≪ t are equilibrated at zero temperature (which in this case is the temperature of the thermal bath) and therefore frozen. Modes with τ ≫ t, although dynamically evolving, are also frozen because the barriers (in this case entropic barriers) are too high to allow for relaxation within the time-scale t.
Only those modes whose characteristic time is τ ∼ t are relaxing at a given time-scale t. We get for the time dependent energy scale ǫ * and the effective temperature the relation ǫ * ∼ log(t)/β eff and this yields the leading behavior According to what has been said one can impose the following Ansatz solution for the g k (ǫ). If g eq k stands for the equilibrium density at T = 0 (i.e., according to eq. (16), g eq k = g(ǫ)δ k,0 ) then we have, where ∆P k ≡ P k − P eq k = P k − δ k,0 andĝ k (x) decays pretty fast to zero for x > 1. This expression tells us the following. Above ǫ * the g k (ǫ) have relaxed to their corresponding equilibrium distributions at the temperature of the bath (in this case the bath is at zero temperature). On the other hand, in the sector of the energy spectrum where ǫ < ǫ * , the densities g k are still relaxing (specially in the region ǫ/ǫ * ∼ O(1)). Since the relaxation is driven by the shift in time of the threshold energy ǫ * the proposed scaling solution Ansatz seems quite reasonable. The prefactor ∆P k /ǫ * is introduced to fulfill condition (6). Furthermore the condition ∞ 0 dxĝ k (x) = 1 is imposed on the scaling functionĝ k .
In the Appendix B we show how this Ansatz closes the set of equations (33) reproducing also the leading asymptotic behavior for ǫ * and z * which turns out to be, For later use we define the following function which scales as function of ǫ/ǫ * . The scaling relation (40) yields the leading asymptotic behavior of all observables different from the occupation probabilities P k . For instance, the energy is given by E = − ∞ 0 dǫ ǫ g 0 (ǫ); using the scaling relation (40) and the asymptotic expression (41) we get for the leading term Note that the asymptotic scaling behavior of the energy is the same as for the effective temperature T eff = 1/β eff in agreement with the quasi-equilibrium hypothesis (see (23)). An important result is that the threshold ǫ * decays slower to zero than the effective temperature. A case where this difference can be clearly appreciated corresponds to the case where the density of states vanishes exponentially fast g(ǫ) ∼ exp(−A/ǫ). In this case ǫ * decays slower than logarithmically, namely like 1/ log(log(t)) (see appendix B for details).

C. Relaxational spectrum in equilibrium
One of the crucial features behind the applicability of the adiabatic approximation is that the long-time behavior at zero temperature finds its correspondence with the low-temperature relaxational properties of the equilibrium state.
To analyze the spectrum of relaxation times τ eq (ǫ) in equilibrium we expand up to first order in perturbation theory the dynamical equations for the g k (ǫ) around their equilibrium solutions g eq k (ǫ). Using the expansion g k (ǫ) = g eq k (ǫ) + δg k (ǫ) we get a set of equations for the variations δg k (ǫ). These are shown in the Appendix C. A complete derivation of the relaxation time τ (ǫ) in equilibrium is complicated. But it is easy to convince oneself that the relaxation time is asymptotically (in the limit T → 0) strongly peaked around the threshold energy ǫ * . For ǫ ≫ ǫ * the relaxation time is small because the population of high energy boxes in equilibrium is rather small. On the other hand, for ǫ/ǫ * ≪ 1 the relaxation is estimated to be finite and independent of T . This result is derived in the aforementioned Appendix C where we show that the maximum relaxation time occurs for ǫ around ǫ * . Starting from eqs. (69) for δg 0 (ǫ) and δg 1 (ǫ) and making use of the adiabatic Ansatz (40), we find, for ǫ ≃ ǫ * , where ǫ * (T ) ∼ T 1/(2+α) , is the asymptotic temperature dependence of the threshold energy at low temperature. This yields for the temperature dependence of the relaxation time, showing there is activated behavior as a function of the temperature but with a relaxation time which increases slower than Arrhenius as T → 0. Note that for the standard BG model corresponding to α → ∞ we obtain an Arrhenius behavior and in the opposite limit α → −1 the relaxation time does not diverge anymore.

VI. NUMERICAL RESULTS
In this section we numerically check the main results obtained in the previous sections. In particular, we want to show the existence of the threshold energy ǫ * separating equilibrated from non-equilibrated energy modes. We have compared three different models characterized by three different types of distributions ( figure (3)). All three distributions were chosen to satisfy the conditions in such a way that the ground state has energy E GS = −1 in the limit N → ∞ for all three cases. The models are the following ones: • Case A: Non disordered model with a gap ( figure (3a)). This is the original BG model where g(ǫ) = δ(ǫ − 1). This case corresponds to α → ∞, therefore ǫ * = 1 and the threshold energy is time independent. The energy is expected to decay like E + 1 ∼ T eff ∼ 1/ log(t). As previously said in section II, the same behavior is expected for any disorder distribution g(ǫ) with a finite gap.
• Case C: Disordered Model without gap and g(0) finite (figure (3c)). We have considered the distribution This case corresponds to α = 0. The energy threshold ǫ * scales like 1/ log(t) and the effective temperature and the energy scale like E + 1 ∼ T eff ∼ 1/(log(t)) 2 . In figure (4) we plot the decay of the energy for all three models. Simulations were done for N = 10 4 , 10 5 , 10 6 boxes (the number of particles is identical to the number of boxes) showing that finite-size effects are not big in the asymptotic regime. We show data for one sample and N = 10 6 . We plot the energy as function of time starting from a random initial condition (particles randomly distributed among boxes, E(t = 0) = −1/e). As clearly seen from the figure relaxation is faster for the Case C and slower for the standard BG model (Case A). The different asymptotic behaviors are shown in figure (5). There we plot (E(t) − E GS )(log(t)) λ with λ ≡ α+2 α+1 . To avoid finite-size corrections when the energy is close to its ground state we computed exactly E GS = 1 N (− N r=1 ǫ r + ǫ min ) where ǫ min is the minimum value among the ǫ ′ s. The different curves saturate at a finite quantity corresponding to the asymptotic leading constant. Note that convergence is slow, showing the presence of subleading logarithmic corrections to the leading behavior. Let us now analyze the shape of the probability densities g k (ǫ). For these distribution we also only show results for N = 10 6 because a smaller number of boxes results in higher noise in the curves. The distribution probabilities were numerically computed by binning the ǫ axis from ǫ = 0 up to ǫ = ǫ max where ǫ max is the maximum value of ǫ r among all the N boxes. One hundred bins are enough to see the behavior of the time evolution of the different distributions. In figures (6) and (7) we show the g 0 (ǫ) for cases B and C respectively. Note that the g 0 (ǫ) converge to the asymptotic result g(ǫ) for ǫ > ǫ * in agreement with the adiabatic solution (40) while they are clearly different for ǫ < ǫ * . The value of ǫ * where g 0 (ǫ) deviates from the asymptotic curve g(ǫ) shifts slowly to zero (like 1/(log(t)) 1 2 or 1/ log(t) for cases B and C respectively), as can be seen in figures (6,7)).  We do not show results for the other g k (for instance g 1 ) because they decay very fast to zero (already for t = 2 17 there are no occupied boxes with more than one particle). Instead, in figures (8) and (9) we verify the adiabatic Ansatz, eqs.(40,43), for the densities g 0 and g 1 in the two models B and C. Figure (8) plots G 0 (ǫ) for both models. Figure (9) plots G 1 (ǫ) for both models. We have used the relation (39) together with z * = log(t) + log(log(t)) yielding G 0 (ǫ) = ∆g 0 (ǫ) log(t) + log(log(t)) Note that the scaling is pretty well satisfied and that theĝ k (x) indeed vanishes for x ≃ 1 yielding an estimate for ǫ * is both cases. We find, ǫ * ≃ 6/ log(t) for case B and ǫ * ≃ 12/ log(t) for case C. Note also that the quality of the collapse of the G 0 is slightly worse for case B than for case C. We think this is due to the stronger subleading corrections to the shift of ǫ * which decays slower to zero for case B. Hence the asymptotic regime is reached only for later times. Indeed, as figure (6) shows, the value of ǫ * obtained within our time-scales has not yet reached the maximum of the distribution g(ǫ), so that we are still far from the asymptotic behavior g(ǫ * ) ∼ ǫ * . Yet, it is remarkable how well the scaling Ansatz eqs.(40,43) fits the numerical data.

VII. A METHOD TO DETERMINE THE THRESHOLD ENERGY SCALE ǫ *
In this section we are interested in the following question. Is there a general method to determine the energy scale ǫ * without having any precise information about the adiabatic modes present in the system ? In the previous sections we addressed this question by proposing an adiabatic scaling Ansatz to the dynamical equations. Here we propose a general method to determine the energy scale ǫ * from first principles without the necessity of knowing the nature of the slow modes present in the system. Obviously for models such as the standard BG model this energy scale has no role since we know from the beginning that relaxation takes place on a single energy scale.
Consider the following quantity P (∆E) defined as the normalized probability density of having a first accepted energy change ∆E at time t. Let us consider the case of zero temperature where this probability density is defined only for ∆E ≤ 0. If Q(∆E) denotes the probability of proposing an energy change at time t (the move is not necessarily accepted), it is easy to show that P and Q are proportional to each other: where A = 0 −∞ Q(∆E)d∆E is the acceptance rate. The expression for Q(∆E) (and therefore P (∆E)) can be exactly computed. Note that computing Q(∆E) yields all information about the statistics of energy changes, in particular the evolution equation for the energy 1 . On the contrary, given the time evolution for the energy this does not necessarily yield the distribution Q(∆E). For the DBG this function can be exactly derived (its derivation is shown in the Appendix D). Here we quote the result, Using the scaling Ansatz eq.(40) we obtain the simple scaling scaling relation, A collapse of different P (∆E) for different times can be used to determine the time evolution of ǫ * . In figure (10) we show the scaling of P (∆E) for the model B for N = 10 4 and different times t = 10 2 , 10 3 , 10 4 , 10 5 . Starting from a random initial configuration, statistics has been collected over approximately 30000 jumps for every time. In figure (11) we check the scaling relation (55) plotting P (∆E)ǫ * as function of ∆E/ǫ * where we have taken ǫ * ∼ 1/ log(t). Note also that the range where P (∆E) is finite corresponds to the region where ǫ ∼ ǫ * . In figure (11) this corresponds to ǫ * ≃ 6/ log(t) in agreement with what was observed in figures (8,9). The scaling works pretty well showing how this method could be used to guess the time evolution of the energy threshold ǫ * in general glassy models in those cases where there exists an energy threshold.

VIII. CONCLUSIONS
In this paper we have considered a new solvable glass model, the disordered Backgammon model (DBG). The new ingredient of this model is that each box has associated a positive random energy which is obtained from a distribution g(ǫ). Again, similarly to its predecessor (the Backgammon model (BG)), the model displays slow relaxation due to the presence of entropic barriers. Actually, it turns out that the relaxation at T = 0 of the number of empty boxes and all other occupation probabilities P k is exactly the same as the original BG model, and, in particular, independent of the disorder distribution g(ǫ). In general, the relaxation of other quantities such as the energy and other disorder dependent observables, displays an asymptotic relaxation which depends on the statistical properties of g(ǫ) in the limit ǫ → 0. In the asymptotic long-time regime, relaxation takes place by diffusing particles among boxes with the smallest values of ǫ. Therefore the asymptotic decay of the energy as well as that of other observables only depends on the exponent α which describes the limiting behavior g(ǫ) → ǫ α . The original BG model is recovered as a special case in the limit α → ∞.
We have written integral equations for the densities g k (ǫ). These equations form a hierarchy of dynamical equations which can be closed by introducing a suitable generating function. We focused on the solution of this hierarchy in the particular case of zero temperature. In this case the analytical solution of these equations proceeds in two steps. Firstly, the equations for the occupancies P k are exactly the same as in the original BG model and they can be solved by using known analytical methods. Secondly, this information is used to guess an adiabatic solution for the g k (ǫ) in terms of a time-dependent energy threshold ǫ * .
All densities are out of equilibrium but admit a scaling solution of the typeĝ k (ǫ/ǫ * ), with the condition that, for ǫ > ǫ * ,ĝ k decays to zero fast enough to guarantee that ∞ 0 dxĝ k (x) = 1. That means that for energies above the threshold the modes are almost completely thermalized.
This Ansatz solution yields two types of leading behaviors: on the one hand, they yield the asymptotic long-time behavior of ǫ * , z * at zero temperature; on the other hand they produce the low-temperature behavior of the values of ǫ * , z * in the stationary equilibrium limit t → ∞. The adiabatic approximation is nothing else but stating the validity of the complementary description of these two very different regimes. On the one hand, the equilibrium regime, where first the limit t → ∞ is taken and later T → 0. On the other hand, the far from equilibrium regime, where the limit T → 0 is taken first and later t → ∞. The commutation of these two limits allows for the interchange between different variables such as energy, temperature and times when expressed in terms of their asymptotic leading behavior. Knowing the leading behavior of the quantities z * and ǫ * , dimensional reasoning as presented in Appendix B yields the leading behavior of β eff which turns out to be proportional to z * /ǫ * . The DBG model offers a scenario where there are two energy sectors separated by the energy scale ǫ * which have very different physical properties. These two sectors manifest in the behavior of observables like the probability densities g k (ǫ) where the time dependent threshold ǫ * separates the equilibrated modes (ǫ > ǫ * ) from the modes which stay off-equilibrium (ǫ < ǫ * ). In the off-equilibrium regime entropic barriers are typically higher than the time dependent barrier at the threshold level ǫ * . For ǫ > ǫ * barriers are lower and equilibrium is fastly achieved.

IX. DISCUSSION
What is the interest of the existence of this energy scale in the context of real glasses? One of the most interesting properties of the present model is that it introduces, in a very simple way, the concept of a threshold energy scale. To our knowledge such concept has never been discussed in any one among the plethora of mean-field glassy models studied during the last years. In those cases one studies relaxation of global quantities which get contribution of all possible energy scales involved into the problem. In principle, nothing is wrong with that since macroscopic observables are those quantities which are always measured in the laboratory. The problem arises when tackling issues related to the violation of the fluctuation-dissipation theorem and concepts such as the effective temperature and partial equilibration. Usually an effective temperature is defined in terms of the measured dissipation (response) and fluctuations (correlations) in the aging regime. This effective temperature is supposed to quantify the amount of energy transfer when the system is put in contact with a thermal bath and behaves in several aspects as a real temperature 14,16,34 . There is a fundamental problem with this definition which is the following. Suppose one takes a piece of silica well below T g (for instance, at room temperature). The piece of silica is not in equilibrium (actually, it is always relaxing even if the relaxation rate can be extremely small and unobservable at room temperature) so one would be tempted to claim that its effective temperature (that describing equipartition among the set of nonthermalized degrees of freedom) is around T g well above the room temperature. Obviously, if we touch a piece of silica, then the hand plays the role of a thermal bath at the room temperature. Therefore, why we do not feel the effective temperature which may be hundreds of degrees above the room temperature? Note that the content of energy of the glass is still very large. Indeed, if the glass suddenly crystallized it would liberate all its latent heat 33 . There are two possible explanations to explain this discrepancy. The first one was analyzed in the context of the oscillator model and assumes that the thermal conductivity is so small that heat transfer is negligible such a short time-scale 34 . The other explanation is that, when touching the glass, we are not touching the slow collective degrees of freedom which still contain a lot of energy but the thermalized degrees of freedom. These two explanations are not totally exclusive. Assuming the existence of an energy threshold ǫ * such that, above ǫ * , all collective modes are thermalized at the bath temperature and below that threshold they are off-equilibrium, this offers an explanation about why when touching a piece of glass we feel it at the room temperature. Our hand only couples to the higher energy degrees of freedom (phonons) and not to the (much hotter) collective excitations. Yet, we cannot exclude that even if we coupled the bath to the hottest collective degrees of freedom then the conductivity would be extremely small and no heat transfer would be measured. From a different point of view, the two different explanations for the small amount of heat transfer established between a bath and the "hot" glassy system reduce into a single one: the existence of a threshold scale ǫ * is consequence of the highly different orders of magnitude for the conductivities in the two energy sectors. Future studies of other glassy models will better clarify this issue 35 .
Finally, we have proposed a method to determine the threshold energy scale ǫ * by computing the general probability distribution Q(∆E). Preliminary investigations in other glassy models show that this distribution provides a general way to determine the threshold scale ǫ * . Moreover it gives interesting information about fluctuations in the aging state although future work is still needed to understand better its full implications in our understanding of the aging regime.
Further investigations in this model will address other issues such as the measurement of effective temperatures. For instance, it would be interesting to understand how the effective temperature, defined as the temperature of the thermal bath which does not produce a net thermal current when put in contact with the system, depends on the energy sector ǫ probed by the bath. By coupling the bath with a selected set of modes of energy ǫ we can understand whether there is a single effective temperature T eff for all modes or rather, there is an ǫ dependent temperature. Note that local equilibrium in this model is only valid in the energy sector ǫ > ǫ * and it could well be that there does not exist a well defined effective temperature in the other sector ǫ < ǫ * . Nevertheless, the most natural possibility is that thermal fluctuations in the off-equilibrium sector ǫ < ǫ * are determined by the effective temperature (see eq.(32)) which determines the relaxation rate of the slow collective mode P 0 . Future studies should enlighten this and other related questions.
Using eqs. (3)(4)(5)(6) and the following identities, we get the equation of motion for g 0 (ǫ) (namely eq. (24)): We then consider the evolution of the probability density for boxes containing one particle. In the following table we list the processes contributing to the evolution of such occupation probability density.
Departure boxes are chosen with probability n d /N . Arrival box are chosen with uniform probability 1/N . Using again eqs. (3)(4)(5)(6) and eqs. (58-60) we are able to derive the equation of motion for the probability density of boxes with one particles and energy equal to ǫ: For densities of boxes with k > 1 particle the scheme of the contributions is presented is the following table: Using the standard adiabatic results (31), P 0 = 1 − 1/z * , P 1 = 1/(exp(z * )) we obtain the results (41,42). Note that the set of equations for h k are still impossible to solve. Only in certain regimes such as ǫ ≪ ǫ * it may be possible to obtain results. There is a set of equations which couples the different h k . But this set of equations is time independent and should yield all the scaling functionsĝ k (x) once appropriate treatment is taking of the amplitude constant which fixes the leading behavior of ǫ * .
We also consider, as an example, the case in which the probability distribution of the quenched disorder becomes exponentially high at high values of ǫ and zero for low values, namely we choose For this choice ǫ 0 dǫ ′ g(ǫ ′ ) ∼ −ǫ exp − A ǫ − A Γ 0, A ǫ , where the generalized Euler function Γ(0, x) goes to zero as x → ∞. In order to estimate ǫ * from eq. (63) we notice now that for P 1 eq. (65) is still valid, while for ∆P 0 we obtain eventually yielding .
In order to estimate the relaxation characteristic time to equilibrium at low temperature we can expand eqs. (69). First we introduce the asymptotic threshold energy ǫ * (T ) as the energy discriminating between thermalized and not thermalized collective modes at temperature T . If we define it through the relation ǫ * (T ) = T z(T ) and we use the relation (22) obtained by doing a low T expansion then we get, where z 0 is the coefficient of the leading term of z(T ) at low T (see eq. (22): z(T ) = z 0 T 1+α 2+α . Then we expand eqs. (69), take ǫ ≃ ǫ * and introduce the following adiabatic Ansatz, Note that this solution is equivalent to the Ansatz eq. (40) introduced for the asymptotic dynamics at zero temperature but with a static ǫ * (T ) now replacing the dynamical threshold. Now consider eq.(69) for δg 0 (ǫ). Because δP k = dǫδg k (ǫ) it can be shown that the slowest mode corresponds to k = 0, i.e. δg 0 (ǫ) ≫ δg k (ǫ) for k > 0. Therefore the second term in the right hand side of eq.(69) dominates the first and the second terms. Introducing (76) into eq.(69) we get that the relaxation time behaves like, τ eq (ǫ * ) ∝ e βǫ * βǫ * For ǫ >> ǫ * the relaxation time is much smaller, since those are the modes with lower energy barriers.
Appendix D: probability distribution of proposed energy updates.
occupation contribution to E ′ − E probability n d = 1 n a = 0 −ǫ d + ǫ a g 1 (ǫ d ) g 0 (ǫ a ) n a > 0 −ǫ d g 1 (ǫ d ) [g(ǫ a ) − g 0 (ǫ a )] n d > 1 n a = 0 ǫ a g 0 (ǫ a ) 1 N p n p [g(ǫ d ) − g 1 (ǫ d )] n a > 0 0 [g(ǫ a ) − g 0 (ǫ a )] 1 N p n p [g(ǫ d ) − g 1 (ǫ d )] The probability distribution Q(∆E) of proposed energy updates is the average of all possible changes, each computed with its probability: where ∆E is the proposed update, E is the energy of the system before the updating and E ′ the energy afterwards. This means The term with δ(∆E) is the term responsible for diffusive motion of particles. Such a contribution does not actually give any contribution to the relaxation of the system and therefore we will not consider it from now on. The probability distribution of accepted energy changes is given by where W (β∆E) is the Metropolis function The normalization factor is As T → 0 the distribution P becomes with A = ∞ 0 dǫ ′ ∞ ǫ ′ dǫ g 0 (ǫ ′ ) g 1 (ǫ) + (1 − P 0 )P 1 . The normalization factor A is actually the acceptance rate of the Monte Carlo dynamics: as it was defined in 7 .
Using the same notation we can write the energy evolution as The r.h.s. of this equation can be equivalently obtained following the procedure presented in appendix A. Indeed, by definition of energy density, is Inserting eq. (24) in eq. (86) we get eq. (86) back.