{'fraction_non_alphanumeric': 0.02900769763198188, 'fraction_numerical': 0.011179394676678483, 'mean_word_length': 4.4389687428701805, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Validity of the zero-thermodynamic law in off-equilibrium coupled harmonic oscillators

In order to describe the thermodynamics of the glassy systems it has been recently introduced an extra parameter also called effective temperature which generalizes the fluctuation-dissipation theorem (FDT) to systems off-equilibrium and supposedly describes thermal fluctuations around the aging state. Here we investigate the applicability of a zero-th law for non-equilibrium glassy systems based on these effective temperatures by studying two coupled subsystems of harmonic oscillators with Monte Carlo dynamics. We analyze in detail two types of dynamics: 1) sequential dynamics where the coupling between the subsystems comes only from the Hamiltonian and 2) parallel dynamics where there is a further coupling between the subsystems arising from the dynamics. We show that the coupling described in the first case is not enough to make asymptotically the effective temperatures of two interacting subsystems coincide, the reason being the too small thermal conductivity between them in the aging state. This explains why different interacting degrees of freedom in structural glasses may stay at different effective temperatures without never mutually thermalizing.


I. INTRODUCTION
The dynamics of glassy systems has been a subject of intensive research 2 . Despite the fact that glassy systems are off-equilibrium systems, some regularities that allow the rationalization of the problem have been found. One of the most striking regularities is the presence of aging. This means that the correlation and response functions are not only functions of time-differences but also of the time elapsed since the system was prepared 6 . Thus, qualitatively, the longer one waits in the low temperature phase, the smaller the response to an external field will be. A salient feature of systems in equilibrium is the fact that the linear response functions and the equilibrium fluctuations are related by the well known fluctuation-dissipation theorem (FDT) 1 . This relation does not hold for off-equilibrium systems. Several studies of spin-glass mean-field models have shown that a generalization of the fluctuation-dissipation theorem is possible through the definition of the "fluctuation-dissipation ratio" (FDR) 3,4 : which is equal to 1 in equilibrium. It turns out that the behavior of the quantity X(t, s) is non trivial in the limit t, s → ∞. If the lowest time s is sent to infinity the quantity X(t, s) becomes a non-trivial function of the autocorrelation C(t, s). This a strong statement which has been proved to hold in the framework of mean-field spin glasses 3,4 . Moreover, it has been recently recognized that the quantity X is generally related to the Parisi order parameter P (q) which appears in equilibrium studies of spin-glasses providing a natural link between the static and dynamical properties 5 . What is the physical interpretation of X? According to relation (1) the fluctuation-dissipation relation would be satisfied if the temperature into the right hand side of (1) were T /X(t, s). This last ratio receives the name of effective temperature and it has been shown 11 that it has some of the good properties of a macroscopic temperature. In fact a proper thermometer coupled to the slow degrees of freedom can measure it. The value of T eff (t, s) = T /X(t, s) would then be different (and higher) than that of the thermal bath. The question about the convenience of this temperature to describe the non-equilibrium behavior has been a subject of controversy in the last years 7 . While there are some evidences (not only theoretical but also experimental 14,2 ) that the violation of FDT gives a good temperature in the thermodynamic sense, it is unclear what properties of standard (i.e. equilibrium) temperatures are common to the non-equilibrium ones.
The motivation of this paper is to answer to the following question: How effective temperatures equalize when two systems out of equilibrium are put in contact? In other words, does there exist a zeroth law for non-equilibrium systems? Let us imagine about a vitrified piece of silica quenched to the room temperature. Because the glass is off-equilibrium its effective temperature is higher than room temperature. But, if we touch the piece of glass it is not hotter than the room temperature. We must conclude that some degrees of freedom within the piece of silica are thermalized to the room temperature while other remain non-thermalized and still hotter. Touching the piece of silica we feel the fast modes, not the slow ones. This poses the question, how is that possible that different interacting degrees of freedom have not reached thermal equilibrium for sufficient long times? Despite of some considerations present in the literature 11,12 there are no clear answers to this question. We believe that some of them may require a more deep understanding through a detailed analysis of an illustrative example as a previous stage to offer more simple and generic considerations. It is our purpose here to follow this route trying to give a general answer to this question by deriving exact results in the framework of a solvable model.
The model is a set of harmonic oscillators evolving by Monte Carlo dynamics introduced in 8 (hereafter referred as BPR model). The importance of this model relies on the fact that it is exactly solvable and shows one of the main features of glasses, namely aging in correlation and response functions. Our interest will be in considering two coupled sets of harmonic oscillators. Thus, we can see how the main observables are affected by the coupling, in particular how the effective temperature evolves for the two sets of interacting degrees of freedom (represented by the two different sets of harmonic oscillators). The interaction may then appear through the Hamiltonian or through the Monte Carlo dynamics itself. We will discover that the effective temperature for the two sets of oscillators depends on how the coupling is done, and we will understand why in vitreous systems different degrees of freedom may stay at different temperatures without thermalising at very long times. The central idea is that interacting non-equilibrium systems each one with very different effective temperatures may not equalize because the conductivity in the aging state can be extremely small. In this sense the utility of the extension of the zeroth thermodynamic law to the non-equilibrium aging state is questioned due to the smallness of the non-equilibrium conductivities.
The paper is organized as follows. Section II describes the main aspects as well as the interest of the model. Section III describes the two classes of couplings we have considered. Section IV analyzes the case in which the main coupling is ruled by the Monte Carlo dynamics. Section V describes the case where coupling appears only in the Hamiltonian. Sections IV and V show how to solve the dynamics of the system. The reader who is not interested in technical issues can skip them. Section VI discusses the results and the physical consequences of our work. The last section presents the conclusions. Three appendices are devoted to some other technical issues.

II. A SIMPLE AND SOLVABLE MODEL OF GLASS
As a simple model of glass we will consider a system of uncoupled harmonic oscillators evolving with Monte Carlo dynamics. The Hamiltonian is: This model was introduced in 8 and was also reviewed in 9,10 . The low-temperature Monte Carlo dynamics of an ensemble of linear harmonic oscillators shows typical nonequilibrium features of glassy systems like aging in the correlation and response functions. The interest of this model is that the slow dynamics at low temperatures is a consequence of the entropy barriers generated by the low acceptance rate. The simplicity of this model makes it exactly solvable yielding a lot of results about the non-equilibrium behavior.
The Monte Carlo move consists on the following: the x i are moved to x i + r i / √ N where r i are random variables Gaussian distributed with zero average and variance ∆ 2 . The move is accepted according to the transition probability W (∆E) which satisfies detailed balance: W (∆E) = W (−∆E) exp(−β∆E), where ∆E is the change in the Hamiltonian. In Appendix A we show the computation of the correlation and response functions. Here we only quote the main results, 1. Slow decay of the energy. The evolution equation for the energy is Markovian. This simplicity allows for an asymptotic large-time expansion showing that the energy decays logarithmically E(t) ∼ 1/ log(t) and the acceptance ratio decays faster A(t) ∼ 1/(t log(t)).

2.
Aging in correlations and responses. The correlation function C(t, s) is defined by: The response function is calculated by applying an external field to the system. Then, the response function is the variation of the magnetization of the system when the field is applied: with the magnetization given by, Details on how to solve correlations and responses are given in Appendix A. The final results are equations (102,111). Both correlation and responses show dominant s/t scaling with logarithmic corrections. The asymptotic scaling behavior is given by, with C(s, s) = 2E(s) K , G(s, s) = f (s) K where the expression f (t) is given in equation (98) and L(t) ∼ t(log 2 (t)) 15 . The slow decay of the response function shows the presence of long-term memory which manifests as aging in the integrated response function 8 .
3. The effective temperature. As said in the introduction, the effective temperature is defined in terms of the FDR eq.(1): In equilibrium E(s) = T /2 and we recover the expected result T eff = T . Interestingly (7) yields a result for T eff which only depends on the smallest time s. The unique dependence of the effective temperature on the lowest time s is generally believed to be satisfied in the asymptotic large s limit for generic structural glasses and spin-glass models with a one step of replica symmetry breaking. This expectation holds here for all times. At zero temperature when slow motion sets in, the system never reaches the ground state and ages forever. In this regime the effective temperature verifies in the long-time limit (i.e s −→ ∞ ): This gives a thermodynamic relationship between the effective temperature and the dynamical energy in the offequilibrium regime showing how the equipartition theorem can be extended to the glassy regime. The effective temperature measures how a quasi-stationary or adiabatic hypothesis is exact for the present model suggesting that some features of equilibrium thermodynamics may be applied to the aging regime.

III. TWO COUPLED SYSTEMS
Now we consider the case in which we couple two systems of harmonic oscillators. In this case it is possible to compute analytically how one system affects the other without loosing the benefit of evaluating Gaussian integrals. The Hamiltonian we have to deal with is: where we take K 1 K 2 > ǫ 2 , otherwise the system has no bounded ground state. We define the following extensive quantities (per oscillator): where E 1 and E 2 are the energy of the bare systems while Q is the overlap between them. In this case we also consider Monte-Carlo dynamics, where the transition probability is performed by the Metropolis algorithm which satisfies detailed balance. The random changes in the degrees of freedom x i , y i are defined in the same way we have explained in the previous section for the case of a single system. But there are different ways to implement the dynamics in the model depending on the updating procedure of the variables x i , y i . Here we have analyzed two important and different procedures which yield quite different results: 1. Uncoupled or sequential dynamics. In this case the two sets of variables x and y are sequentially updated. First the x i variables are updated and the move is accepted according to the total change of energy ∆E = ∆E 1 −ǫ∆Q. Next, the variables y i are changed and the move accepted according to the energy change ∆E = ∆E 2 − ǫ∆Q. This procedure is then iterated. In this case, the dynamics does not affect simultaneously the two sets of variables but each set is updated independently from the other. The only coupling between the two sets of oscillators comes from the explicit coupling term ǫQ in the Hamiltonian. Note that for ǫ = 0 the dynamics becomes trivial because the dynamical evolutions are that of two independent sets of harmonic oscillators everything reducing to the original model described in section II.
2. Coupled or parallel dynamics. In this first case the x i , y i variables are updated in parallel according to the rule The transition probability for that move W (∆E) is determined by the change in the total energy ∆E = ∆E 1 + ∆E 2 − ǫ∆Q introducing, on top of the explicit coupling term ǫQ in the Hamiltonian, an additional coupling between the whole set of oscillators through the parallel updating dynamics. Contrarily to the uncoupled case, the ǫ = 0 case is interesting by itself because it shows how this kind of dynamical coupling strongly influences the glassy behavior. In fact, in the limiting case ǫ = 0, there will be some changes which make the energy of one of the two systems increase, this change being accepted because the total energy will decrease. Because of that, despite of the fact that there is no direct coupling in the Hamiltonian the dynamics turns out to be strongly coupled.
In what follows we describe the main set of quantities we are interested in. The solution of the dynamical equations for the coupled and uncoupled cases is very similar. The Appendix B shows in detail the derivation of the dynamical solution for the uncoupled case.

A. Correlation, overlaps and responses
On top of the time evolution of one-time quantities our interest will also focus on the behavior of two-times quantities such as correlations and responses. These quantities will refer to three classes of systems: the set of oscillators described by the x variables, the set of oscillators described by the y variables and the whole set of x and y variables. In the rest of the paper, as a rule, the subindex 1 will refer to quantities describing the set x of oscillators, the subindex 2 will refer to quantities describing the set y of oscillators and the subindex T will refer to quantities describing the whole set of oscillators x plus y. The main set of correlation and response functions we are interested in are: • Correlations. The correlation function for the sets x and y, as well as the global correlation C T (t, s) = 1 2 (C 1 (t, s) + C 2 (t, s)). • Overlaps. These are cross-correlations involving different sets of variables: with Q 1 (t, s) = Q 2 (s, t). As we will see later, it is useful to define these two functions Q 1 , Q 2 which essentially are the same overlap function but acting on different time sectors.
• Response functions. The response function for the sets x and y are defined in the following way. Define the magnetizations for the two sets of oscillators x and y, Consider also two external fields h 1 and h 2 conjugated respectively to M 1 and M 2 , We define four types of response functions G 1,2 , G ′ 1,2 . The G 1 (t, s), G 2 (t, s) functions measure the change in the magnetization M 1 (t), M 2 (t) induced by their respective conjugated field h 1 , h 2 applied at time s. These are defined by where the index i = 1, 2 represents each one of the systems. Apart from these two response functions we may define the global response function G T (t, s) as the change in the global magnetization M T = 1 2 (M 1 + M 2 ) induced by a field conjugate to the total magnetization, The primed response functions G ′ 1 (t, s), G ′ 2 (t, s) functions measure the change in the magnetization in each set of oscillators M 1 (t), M 2 (t) induced by a conjugated field (respectively h 2 , h 1 ) applied on the other set of oscillators at a time s: where the indices i, j = 1, 2 are different i = j. In the absence of a coupling term ǫQ in the Hamiltonian (9) the two response functions G ′ 1,2 vanish but for ǫ = 0 they enter into the solution of the dynamical equations. • Effective temperatures. From the correlation and response functions we may define three effective temperatures: T 1 eff for the system 1, T 2 eff for system 2 and T T eff for the global system. These are defined as follows, We will analyze in detail the three effective temperatures for the coupled and the uncoupled cases. From them we will learn whether the systems equalize their temperatures and how they do.

B. Equilibrium regime
Here we present the results for the statics for the general model (9). The equilibrium solution is the stationary state of the dynamics coinciding for both coupled and uncoupled dynamics. The results for the one-time quantities can be simply evaluated from the partition function, which involves simple Gaussian integration. By performing the appropriate partial derivatives we calculate the different thermodynamic quantities: where the parameters I, J are defined by, where the total energy E = E 1 + E 2 − ǫQ is given by the equipartition relation E = T . Note that K 1 K 2 − ǫ 2 > 0 in order for E eq 1 , E eq 2 to be positive. The equilibrium correlations C, overlaps Q and responses G, G ′ only depend on the time differences. While the precise form of these functions depends on the particular type of dynamics, the magnetic susceptibilities do not. These are given by: and are temperature independent as expected for oscillator systems. Nonetheless, in equilibrium the three effective temperatures (18) coincide with the bath temperature T .

IV. THE DYNAMICALLY UNCOUPLED (OR SEQUENTIAL) CASE
In this section we solve the dynamics of the thermodynamic relevant quantities for the case in which the two subsystems of oscillators are dynamically uncoupled. As explained in the previous section, in this case we make a sequential dynamics avoiding direct dynamical coupling effects coming from the Monte Carlo dynamics. The derivation of the dynamical equations is explained in the Appendix B. The equations for the energies and overlap (10) are written down in (126,127,133), with the following definitions: where the error function was defined in (97). Definitions (27) hold for i = 1, 2, each i representing one of the two systems. Note that the whole dynamics is contained in the function f Ri(t) . In what follows we will be especially interested in the zero-temperature case where relaxation time diverges and dynamics is slow and glassy. For T = 0 the function f Ri(t) in (27) becomes A. Asymptotic long-time expansion for the one-time quantities The asymptotic solution of equations (23,24,25) may be guessed from the behavior of the energy (113) for the uncoupled systems. Trying a solution of the type we can only solve the asymptotic behavior in the limit ǫ ≈ 0. This is a consequence of the fact that the quantities R 1 and R 2 are different in general and we have a system of four equations with three parameters. There is not any general solution for this system, but in the limit ǫ ≈ 0 the quantities R 1 and R 2 become identical to the corresponding energies to leading order yielding only three equations with three unknown parameters (a, b, c). In this limit the value of the coefficients a, b, c may be easily obtained yielding where J was defined in eq.(21).At first order in logarithmic corrections 1/ log(t) we find in the limit ǫ ≈ 0 : Note that, in the long time limit both energies E 1 , E 2 tend to zero logarithmically but their relative difference E1−E2 E1 stays finite. In the limit of small coupling constant we can do a more refined expansion yielding: where we have put explicitly the terms of order 1/ log 2 (t) as sub-dominant corrections. These terms come from the fact that the true expressions for the energies and the overlap should be, in order to match the coefficients in (23,24,25): which gives a correction of order 1/ log 2 (t) in expressions (31).
In Fig.1 we show the evolution for the energies and the overlap for two systems of harmonic oscillators with a small value of ǫ. We also show the asymptotic behavior (32,33,34). We can see that the two energies remain different even at long times. We will see that this feature is very important for describing the non-equilibrium state of the whole system. We can also see that the asymptotic expansions are in good agreement with the numerical solution of the dynamic equations. Nevertheless, there are systematic deviations at long times being consequence of the limited range of validity (ǫ ≪ 1) of the asymptotic solution (32,33,34). If K 1 ∆ 2 1 = K 2 ∆ 2 2 the energies of the two oscillators become identical (note that equations (96,97,98) only depend on the constant K∆ 2 ) and asymptotic dynamics also.

B. Correlations and responses.
The set of equations for the four correlation and overlap functions defined in (11,12) can be written as: with the subsidiary boundary conditions The equilibrium solution can be appropriately worked out because the matrix coefficients are time-independent. If we write the matrix equation in compact form The precise results for correlations and overlaps are reported in the Appendix B (formulae (151-159)), the initial conditions being given in (37). In the non-equilibrium case it is not possible to write down an exact solution for equation (38) for any value of ǫ. The formal solution of (38) is which can be worked out perturbatively up to any order around ǫ = 0 (T stands for the time ordered product). In the Appendix C, we give some details how to construct such expansion. Up to order ǫ 2 the solution for the components C 1 , C 2 of the four component vector C are Similar expansions are obtained from Q 1 (t, s), Q 2 (t, s). Here we do not report them because correlations are enough to analyze the effective temperatures. Similarly we can also obtained expressions for the responses as detailed in the Appendix B. The time evolution for the four possible response functions is given by As explained in Appendix B these equations must be solved with the subsidiary boundary conditions, The initial conditions for G 1 , G 2 come from the delta-terms in their equations. The other two initial conditions for G ′ 1 , G ′ 2 come from the fact that there is no discontinuous jump in the response function of one system when we apply the field to the other system. This result also holds in the framework of the Langevin dynamics and manifests in the equations for the magnetizations (see in Appendix B (143,145)) as the absence of a field h 2 in the equation for M 1 and the absence of a term h 1 in the equation for M 2 .
In equilibrium the expressions for the responses G 1 , G 2 are given in the Appendix B. Up to order ǫ 2 the expression for the off-equilibrium responses G 1 , G 2 can be solved analogously as done for the correlations and are given by (49)

V. THE DYNAMICALLY COUPLED (OR PARALLEL) CASE
For the dynamically coupled case the calculations proceed similarly as to the previous dynamically uncoupled case. The evolution equations for the overlap and the energies (10) are: I and J were defined in (21) and the new quantities E N , f EN and E T (E T is not the total energy) are given by and the expression (54) at zero-temperature is A. Asymptotic long-time expansion for the one-time quantities Proceeding similarly as done in the former section we can find the asymptotic expressions for the energies and overlaps. In this case we can find a solution for ǫ finite due to the fact that in this case we have only one dynamic function f EN . We find, to leading order in 1/ log(t) Note that, contrarily to results (31) for the dynamically uncoupled case, the energies E 1 , E 2 asymptotically coincide and the relative difference (E 1 − E 2 )/E 1 vanishes in the long-time limit. This difference of behaviors is not casual and has a physical interpretation that we will discuss later. The more precise expansion turns out to be, ) .
(59) Let us stress that, contrarily to the dynamically uncoupled or sequential case the previous expressions are valid to any order in ǫ. The origin of the 1/ log 2 (t) terms in previous expressions is the same as in the uncoupled case (32,33,34). In Fig. 2, we show the numerical solution for the evolution of the energies and the overlap as well as the asymptotic expansions (57,58,59).
We have said that the relative difference (E 1 − E 2 )/E 1 vanishes in the long-time limit. It is not difficult to see how this happens. The time-evolution for the quantity E 1 /E 2 is easy to derive from eqs.(50,51,52) in the asymptotic long-time limit E N → 0. One then finds the following expansion to leading order If K 1 ∆ 2 1 = K 2 ∆ 2 2 the correction is even smaller. Interestingly, this leading correction does not depend on ǫ showing that the two energies E 1 , E 2 approach each other at a rate determined by the fact that the whole dynamics of the model is coupled and not by the fact that the two oscillator systems are coupled by the presence of a term ǫQ in the Hamiltonian. The behavior of this quantity is shown in figure (4) for different values of ǫ together with the asymptotic expansion (60).

B. Correlations and responses.
Following similar methods as for the dynamically uncoupled case presented before we can write down the equations for correlations and overlaps with the subsidiary boundary conditions given by equations (37).
In matrix form these equations reduce to the equation (38). As explained for the uncoupled case, this set of equations can be exactly solved only in the equilibrium regime where the coefficients are time-independent. The solution is then given by the equation (39), the expressions for correlations and overlaps are the formulae (151,156) with a = − In the most general case where the coefficients of the matrix equation are time dependent the exact solution can be written in the closed form (40) which can be expanded to any order in ǫ as explained in the Appendix C. As in the uncoupled case we present here the results up to order ǫ 2 only for the correlations, Responses G 1,2 , G ′ 1,2 can be worked out in a similar way as shown in the Appendix B for the uncoupled case, with the subsidiary boundary conditions (note that for G 1 , G 2 they are different from those in (47)), It is a simple exercise to check in equilibrium whether these responses give the correct value of the susceptibility (22). In equilibrium responses only depend on the difference of times. As the susceptibility is just the integral of the response function we can integrate the equations (for simplicity we shall consider that s = 0). Then, the equilibrium susceptibility of every system is just: Integrating the equations for response functions we obtain: At very long times, ergodicity imposes G 1,2 (∞) = G ′ 1,2 (∞) = 0. These equations give the exact results for the equilibrium susceptibilities eqs.(22). In equilibrium the responses can be easily computed and one gets (to keep formulae at minimum we only report the results for G 1 and G 2 ), with the usual expression (155) for λ 1 , λ 2 . For G 1 the values of the constants are: while for G 2 the same results (79) are valid but interchanging the indices 1 and 2.
In the general off-equilibrium case the result for G 1 , G 2 , G ′ 1 , G ′ 2 can be worked out perturbatively. Here we only write the expression up to order ǫ 2

VI. RESULTS AND DISCUSSION
First of all we can check that the equilibrium results are the expected ones. It is easy to prove that, independent of the dynamics, the effective temperatures are just the temperature of the bath: Because in equilibrium the energies of the subsystems are the same (see (20)).

A. Sequential case
In the off-equilibrium case the results are more interesting. It is easy to verify the following expressions for (18) up to order ǫ 2 : The expression for the total effective temperature is just where the correlations are given by expressions (41,42) and the responses are given by (48,49). From the equations (83,84) we can see immediately that the effective temperatures are well defined in the regime in which the ratio t s is finite. Otherwise, the last term in the right hand side of (83) and (84) would diverge. At zero temperature and up to order ǫ 2 in the coupling constant, we have found that E 1 , E 2 , Q decrease logarithmically implying that both f R1 and f R2 decay like 1/t. Now let us consider t, s both large but t − s ≪ s. For a weak coupling (i.e ǫ ≈ 0) the value of the effective temperatures are, in the limit s → ∞ but with t s finite: This yields in the s → ∞ limit a non vanishing relative difference T 1 eff /T 2 eff − 1. This is a consequence of the fact that the two energies are different in the long-time regime. Note that each effective temperature verifies the equipartition theorem in the limit of long times as expected. The physical interpretation is clear: each system is relaxing towards its equilibrium state slowly and at any time we can consider that the systems are at "quasiequilibrium" at their corresponding effective temperatures. Obviously the concept of "quasi-equilibrium" is meaningful in a time window smaller than the characteristic time-scale in which the system relaxes (i.e. during this time-scale the effective temperatures do not change), hence we need to impose that t/s is finite.
Let's think now about the global system. As we have seen, the energies for the two systems remain different even at infinite times. This can be explicitly seen in figure (4) where we show how the relative difference between the energies (or the effective temperatures according to (86,87)) increases monotonically as a function of time (for late times) for any value of ǫ. We may then conclude that a coupling in the Hamiltonian is not enough to reach an equalization of effective temperatures. This difference of the two effective temperatures implies that there are some degrees of freedom hotter than others. One can then imagine that there is always some kind of heat transfer or current flow going from the "hot degrees" of freedom to the "cold" ones. Then, one may ask why the effective temperatures do not asymptotically equalize. The reason is that the off-equilibrium conductivity may vanish with time fast enough for the heat transfer not to be able to compensate such difference. In this situation, if we now compute the total effective temperature (85) for the whole system we see that in the off-equilibrium regime this temperature does not coincide with the sum of the energies of the systems. This fact fortifies the definition of the effective temperature using the FDR (1) in off-equilibrium systems. For two systems in "local" equilibrium at two different temperatures, despite the fact that each system verifies FDT, the sum of the two systems never verifies FDT unless the two temperatures are equal. In our case, we have two systems which are in "quasi-equilibrium" at two different effective temperatures, so the T T eff would never be the sum of the two energies unless the two effective temperatures T 1 eff , T 2 eff were the same. In other words, two systems thermodynamically stable at different temperatures are not globally stable when put in contact.

B. Parallel dynamics
The effective temperatures (18) can be exactly computed to order ǫ 2 . In the equilibrium regime, both the full expression derived from (151,156,77,78) and the general approximate solutions (65,66,80,81) up to order ǫ 2 computed in the equilibrium regime yield the bath temperature for the three effective temperatures. In the non-equilibrium case, up to order ǫ 2 , the results are: The expression for the total effective temperature is just: where the correlations are given by expressions (65,66), and the responses are given by (80,81). As in the case without coupling, the interesting dynamics is when the temperature of the bath is zero. In this case, the energies and the overlap decay to zero logarithmically which implies that f EN (t) vanishes like 1/t. A careful evaluation of the integrals contributing to the ǫ 2 term shown in equations (88,89) reveals that they are a function of t/s which stays finite provided that ratio is finite. As we discussed in the previous uncoupled or sequential case the effective temperatures (88,89) have full sense when we consider times t/s finite so no appreciable transfer of energy between the two systems has still occurred.
It is clear from the asymptotic expressions for the energies and the overlap that in the long-time limit (s → ∞) : While the energies themselves decay as 1/ log(s) the relative difference (Q(s)/E 1 (s) − ǫ 2 K1K2 ) decays like 1/ log(s). Up to order ǫ 2 we may write, in the limit s −→ ∞ (with t s finite): because the asymptotic values of the E 1 (s) and E 2 (s) are the same the effective temperatures for the subsystems become identical in the long-time limit. Note that the case with dynamic coupling or parallel dynamics is qualitatively different from the case without dynamic coupling or sequential, because now all the degrees of freedom are at the same effective temperature in the long-time limit. Moreover, if we consider the global system it is easy to prove that the total effective temperature defined in (18) is, in the limit s → ∞ with t s finite: where E = E 1 = E 2 and Q are given by (56). This is a consequence of the fact that the energies of the two systems equalize due to the dynamic coupling. Then, the whole system has the same effective temperature and we can define an effective temperature for the global system using FDT. The situation is the same as in equilibrium systems. If we have two systems in equilibrium at a certain temperature T , FDT not only holds for each subsystem but also holds for the whole system bringing the temperature of the bath T . At higher-orders in ǫ we expect that all terms with be subleading for t/s to be finite and asymptotically all three temperatures coincide.
If we restrict to the case in which the coupling constant vanishes, ǫ = 0, then the systems are still coupled only through the dynamics and we obtain the same qualitatively results: with E(s) ≈ J 8 log(s) . We conclude that the dynamic coupling does not allow the presence of more than one effective temperature in the whole system because even in the absence of explicit coupling in the Hamiltonian, the dynamics itself makes the energies to equalize in the long-time limit regime.

VII. CONCLUSIONS
In this paper we have solved exactly the dynamics of two systems of harmonic oscillators. We focused our attention on the concept of the effective temperature defined through the FDR eq.(1). The effective temperature, a parameter defined by a relation of the correlation and response functions, has been introduced in the context of glass theory in order to understand the physics behind the dynamic behavior of these out-off-equilibrium systems. In this paper we hope to have clarified some aspects behind the physical meaning of this effective temperature.
We have studied two types of couplings between the two subsystems of oscillators, both in an aging state, finding that the way we couple them is crucial for the validity of the zero-temperature law in the off-equilibrium regime to hold. The two cases we studied are the dynamically uncoupled or sequential case and the dynamically coupled or parallel case. In short, for the sequential case the coupling between the variables of the two subsystems in the resulting dynamics arises only through the Hamiltonian term ǫQ. For the parallel case, the variables of the two subsystems are simultaneously updated leading to further interaction between the two subsystems (on top of the ǫQ coupling term in the energy).
We have discovered that for the dynamically uncoupled or sequential case the two subsystems asymptotically reach different effective temperatures which never equalize. So the whole system is divided in two parts, each part characterized by its own effective temperature. The explanation for this odd behavior lies behind the time dependence of the off-equilibrium thermal conductivity which decays very quickly to allow for an asymptotic equalization of the two effective temperatures. This raises the question whether different interacting degrees of freedom do eventually reach the same effective temperature in the asymptotic regime, condition tightly related to the validity of the zeroth law for the off-equilibrium aging state. Our conclusion is that the zeroth law is probably valid but hardly effective due to the very small conductivity between the two subsystems in the aging state. A calculation of the thermal conductivity in this model will be shown elsewhere 16 and reveals that it decreases very quickly with time, the heat transfer being unable to compensate for the difference of the effective temperatures of the two subsystems.
For the dynamically coupled or parallel case, the two effective temperatures equalize and the two subsystems are in a sort of thermal equilibrium between them in the aging state. Consequently, the union of the two subsystems has an effective temperature which coincides with the temperature of each subsystem. In this case, the direct coupling of the two subsystems through the parallel dynamics makes the conductivity much larger than in the sequential case so in this case a zero-th law for the aging state is effective and holds. In fact, these results are also valid when we consider the particular case ǫ = 0 in which the dynamics in itself is enough to equalize the effective temperatures.
From these two type of couplings the first one is the only realistic. Dynamics in real structural glasses involves short scale motions of atoms and coupling between the different degrees of freedom occurs at the level of the energy or Hamiltonian and never at the level of the dynamics (at least in the classical regime). The results of this paper explain then why different degrees of freedom in structural glasses can stay at different effective temperatures forever. The off-equilibrium conductivity or heat transfer between the different degrees of freedom is small enough for the equalization of the effective temperatures associated to the different degrees to never occur. This explains why when we touch a piece of glass we feel it at the room temperature. The heat transfer coming from the hotter non-thermalized degrees of freedom is extremely small. Before finishing we must note one particular feature of our model. All the calculations were done at zero temperature where the energy vanishes asymptotically. The fact that the energy (and consequently the conductivity) of the system is exhausted in the asymptotic limit can lead to a pathological behavior not present in structural glasses at finite temperature. Nevertheless, the fact that the thermal conductivity vanishes much faster than the energy itself, suggests that the vanishing of the conductivity is not related to zero-temperature dynamics.
In the present calculation we have focused on the interaction between two subsystems, both in the aging state. When one of the subsystems is in an aging state and the other is in equilibrium the analysis proceeds similarly, the conclusion being that the non-thermalized subsystem determines the rate of heat transfer and hence the measurement of the effective temperature. The value of the effective temperature measured by a thermometer and other related questions can be analyzed in detail in the present model and will be presented elsewhere 16 .
To conclude, although a zero-th law for non-equilibrium glassy systems may hold, it is hardly effective because of the small energy transfer occurring between degrees of freedom at different effective temperatures. It would be very interesting to pursue this investigation further by studying other solvable examples and showing that what we have exemplified here is a generally valid for structural glasses as well as for other glassy systems.

APPENDIX A: A SHORT REVIEW OF THE BPR MODEL
In this appendix we show how derive the results for the correlation and the response functions in 8 in order to understand the techniques we will use throughout this paper. In that model, the system is constituted by N uncoupled harmonic oscillators which evolve with Monte Carlo dynamics. The energy of this system is The result for the dynamical evolution for the energy is: where we have defined the quantities The stationary solution is just E = 1 2 T as expected. Another important quantity is the acceptance rate which is the number of accepted Monte Carlo movements at a time t: In the same way we can compute the equation for the correlation function defined as: and the evolution of the correlation function is given by the equation where the quantity f (t) has been previously defined in (98). The solution for the correlation function (which depends explicitly on two times) is: where we have to add the initial condition C(s, s) = 2E(s) K . In order to compute the equation for the response function defined by, we have to consider the Hamiltonian perturbed by a small external field Then we compute the dynamical evolution for the magnetization, which in our model is defined by: where we have defined a 'new' energy as: The quantity f E ′ (t) is identically defined as in (98) but with the new energy E ′ . Note that in the case in which h = 0 the magnetization will always be zero because of the initial condition we consider, i.e M (0) = 0. Also, when we compute the evolution for the response function, the first term in the right hand side of the (105) is just the response.
Then we have to analyze carefully the second term (which is proportional to the external field) in the right hand side of (105). First of all, we write We consider the variation of the magnetization as follows and by keeping only the linear term in h in the second term of the r.h.s in (108) we get The first term in the r.h.s of (109) is obviously zero and only the the second term gives a non-vanishing contribution. So the evolution for the response function is whose solution is Now, we are in position to compute the effective temperature based on the violation of FDT Note that the effective temperature only depends on the smallest time s. This feature is due to the simplicity of the model. In this model, due to the finite amplitude of the Monte Carlo movements the system never reaches the ground state {x i = 0}. In fact, Monte Carlo dynamics induces entropic barriers which manifest as activated behavior for the relaxation time. The interesting dynamics is found when we study the relaxation of the system at zero temperature. To obtain the dynamical equations at zero temperature we have to consider only the negative changes in the energy. It can be seen that in the long time limit the relaxation of the energy is logarithmic moreover, we obtain the following asymptotic behavior of the function f (t) and the acceptance rate: For the long time behavior of the correlation and the response functions we obtain to leading order in log(s)/ log(t)

APPENDIX B: SOLUTION OF THE DYNAMICALLY UNCOUPLED OR SEQUENTIAL CASE
In this appendix we show explicitly the detailed calculations for the case in which we sequentially update the two subsystems. Note that each subsystem is updated in parallel but no simultaneous updating of the whole system is performed so there is no direct coupling of the two subsystems through the dynamics but only through an explicit coupling term ǫQ in the Hamiltonian. We have to take into account this fact when we compute the distribution probability for a change in the energy. The Hamiltonian we have to deal with is The main quantities we work with are where E 1 and E 2 are the energy of the bare systems while Q is the overlap between them. The Monte Carlo updating procedure is the following. First all the x i are moved to x i + r i / √ N where the r i are random variables Gaussian distributed with zero average and variance ∆ 2 1 . The move is accepted according to a rule defined by an acceptance probability W (∆E) which satisfies detailed balance: W (∆E) = W (−∆E) exp(−β∆E), where ∆E is the change in the Hamiltonian. Later all the y i are moved to y i + s i / √ N , where the s i are random variables Gaussian distributed with zero average and variance ∆ 2 2 . The same transition probability is now applied for the y i variables. This sequential updating of the x i and y i variables is then iterated. Note that the coupling in the dynamics only appears through the change ǫ∆Q of the total energy. Now we compute the distribution probability of a change in the energy of the first system. This probability distribution can be expressed and in the same way we can compute the probability for the other system: Using the integral representation of the delta function: we obtain P (δE 1 ) = 1 with the quantities Note that due to the explicit coupling ǫ the probability of a change in the energy of one system not only depends on their energy, but also on the energy of the other system and the overlap. Now, we can compute the evolution of the energies yielding We compute the equation for the evolution of the overlap in two steps. The first is the change in the overlap when the variables of the first system are moved; and the second one is when the variables of the second system are moved. So we must to compute two joint probability distributions Then we compute the evolution equation for the overlap in each step and sum the two equations ∂Q(1st)) ∂t = The solution of these equations is which yields the final equation with the quantities defined in (27,97). In the same way we can compute the equation for the correlation and overlap functions defined in (11,12). To compute their evolution equations we must evaluate the joint probability of a change in the energy and a change in the correlation function. Note that when we consider the change in the variables x i we have to consider the energy of the system one, and when we consider the change in the variables y i we have to take into account the energy of the other system. The joint probability can be decomposed into the probability distribution for a change in the energy multiplied by a conditional probability P (δE 1 , δC 1 ) = P (δE 1 )P (δC 1 |δE 1 ) , P (δE 2 , δC 2 ) = P (δE 2 )P (δC 2 |δE 2 ) , P (δE 2 , δQ 1 ) = P (δE 2 )P (δQ 1 |δE 2 ) , P (δE 1 , δQ 2 ) = P (δE 1 )P (δQ 2 |δE 1 ) .
Then the evolution for the correlation functions can be computed using yielding ∂ ∂t For the response functions we have to compute the dynamic evolution equations for the magnetizations. We consider an external field coupled to each system, so the new Hamiltonian is We define the magnetizations as follows Then, we perform the same steps as we did for the other quantities. First of all we have to compute the joint probability of a change in the magnetization and a change in the energy. For example, for computing the response function for the first system we make h 2 = 0 and h 1 = 0; then we compute the joint probability distribution for a change in M 1 and E 1 . After that we can obtain the evolution for the magnetization of this system Note that in this case we are considering h 2 = 0 but still the equation for M 1 depends on M 2 . For the sequential updating procedure we have to consider the evolution for M 2 (t) with h 1 = 0 and h 2 = 0 which is, by symmetry considerations We finally get the equations for the four different response functions using the same procedure we followed for the single system (see Appendix A). This yields