Glassy behavior in a simple model with entropy barriers

Glassy behavior in a simple model with entropy barriers A. Prados and J. J. Brey Fı ´sica Teo ´rica, Facultad de Fı ´sica, Universidad de Sevilla, Apartado de Correos 1065, E-41080 Sevilla, Spain B. Sa ´nchez-Rey Escuela Polite ´cnica Superior, Universidad de Huelva, E-21819 La Ra ´bida, Huelva, Spain ~Received 13 September 1996! We study the dynamical behavior of a system with a variable number of particles n. The empty state n50 is the ground state, while all the other states n.0 are degenerate in energy. In equilibrium, the mean number of particles is equal to unity, independently of the temperature. The static properties are the same as for the Backgammon model recently proposed by Ritort @Phys. Rev. Lett. 75, 1190 ~1995!#, while a variation of the kinetics is considered. The elementary dynamical processes are the arrival and departure of a particle. The rate of the departure process is constant, while the arrival rate is obtained from the detailed balance condition. Thus, there is no energy barrier separating the ground state n50. Nevertheless, glassy behavior appears due to the presence of effective entropy barriers. At low temperatures, the response functions are shown to obey f (t).exp@2(t/ t ) g #. In thermal cycles of cooling and reheating from low temperatures, the system shows hysteresis, which follows from the trend of the system to approach the normal curve characterizing the heating program. @S0163-1829~97!03710-7# I. INTRODUCTION The study of glassy behavior has been quite an active field in recent years. A review of the main features observed in real glasses, and several microscopic models showing similarity with them, can be found in Refs. 1 and 2. In relaxation experiments, the linear response functions show nonexponential behavior. In particular, a Kohlrausch-Williams-Watts ~KWW!decay is usually found. In cooling experiments a laboratory glass transition, in which the properties defining the state of the system become frozen, is observed. The transition is associated to a fast increase of the relaxation time as the temperature is lowered. During reheating, hysteresis effects show up, with the system returning to equilibrium following a path which is different from the cooling one. A more detailed discussion of the rich phenomenology of glasses is available in Refs. 3 and 4. There is a great variety of models trying to explain glassy behavior. The simplest one is a two-level system ~TLS!, where an energy barrier must be surpassed in order to go from the excited to the ground state.5,6 In some models, the increase of the relaxation time is associated to the introduction of cooperativity in the dynamics of the system,7but there is also an energy barrier separating the ground state from the excited ones. This barrier plays an essential role in the divergence of the relaxation time at low temperatures. On the other hand, entropy is known to play an important role in the description of glassy behavior since the pioneering work by Adam and Gibbs.8 Recently, Ritort9,10 has proposed a model without energy barriers, in the sense that the system can always reach the ground state without any energy-activated process. The dynamical study of the model has focused on thermal cycles of cooling and reheating, and zero-temperature properties as aging. The system displays glassylike behavior, despite the absence of energy barriers, and the divergence of the relaxation time is due to the entropic contribution to free energy barriers. These appear because of the small number of directions in phase space along which the energy decreases. Slow relaxation shows up because the system has to explore a wide phase space region before reaching the ground state. Because of the rules governing its dynamics, the model has been referred to as the Backgammon ~BG!model. It can be visualized in several different, although equivalent, ways. Here we present one of them, while another one is discussed in the final section. Suppose we have a two-dimensional lattice with a particle at each site. Then, an external mechanism is introduced such that particles tend to aggregate in the direction perpendicular to the lattice. Particles remaining on the lattice have a larger energy than those which are aggregate to them, so that the minimum energy is reached when all particles form a unique aggregate at a given site. All sites and particles being equivalent, this state has a degeneration given by the number of sites ~or particles!. The dynamics of the system is defined by means of a Markov process in which each particle can move to any other site, with transitions rates given by Metropolis dynamics. Since the spatial arrangement of the sites in the plane plays no role at all, the model is of a mean-field type. Mean-field approximations are not accurate to describe relaxation through energy barriers in real structural glasses, because of the nucleation processes taking place in them. Nevertheless, as pointed out by Ritort,9 the effect of entropy barriers should not depend very strongly on the range of the interactions and the information obtained from this kind of model is expected to be relevant also in the case of short-ranged interactions. In this work we introduce a model that keeps the main characteristic of Ritort’s model, namely the absence of energy barriers for transitions to the ground state, and allows an analytical treatment of the dynamics. We consider a system with a variable number of particles n, in which the ground state has no particles, n50, while all the states with n.0 PHYSICAL REVIEW B 1 MARCH 1997-IIVOLUME 55, NUMBER 10 55 0163-1829/97/55~10!/6343~13!/$10.00 6343 © 1997 The American Physical Society are degenerate. The statics of this system is equivalent to the BG model with indistinguishable particles.11 The dynamics is formulated by means of a master equation with transition rates verifying the detailed balance condition. The equation can be exactly solved for constant temperature processes, allowing the identification of the mechanisms leading to nonexponential relaxation and to the divergence of the relaxation time. For cooling processes we show the relevance of the relaxation modes of the master equation and of the energy relaxation time to characterize the laboratory glass transition and the freezing temperature, respectively. Along heating, the dynamical behavior of the model is understood from the trend of the system towards a ‘‘normal’’ curve.12 In particular, the hysteresis effect, which is so characteristic of glasses, is directly related to the approach to this normal curve. The existence of such a curve is a quite strong prediction of models bases on a master equation formulation of the dynamics. Whether there is a normal curve also for real structural glasses remains an open question. The results obtained will be compared to other previously considered models and, in particular, the one-dimensional Ising model with Glauber dynamics.13 Let us mention that the Ising model may be relevant in the context of structural glasses, since it has been proved to accurately describe the evolution of the configuration of a one-dimensional system of particles with anharmonic and competing interactions.14 Although energy barriers exist in the model studied in Ref. 13, glassy behavior appears in both cases for similar reasons. Probably, this is also the case for any model showing glassy behavior, as long as its dynamics is described by a master equation. The plan of the paper is the following. In Sec. II the model is formulated, and the master equation describing its dynamics is solved for the constant temperature case. Relaxation properties are considered in Sec. III, focusing on the stretched exponential decay found at low temperatures. Section IIIA is devoted to the study of the equilibrium time autocorrelation function of the energy, while linear relaxation after a temperature perturbation is the subject of Sec. IIIB. Thermal cycles are studied in Sec. IV, and cooling processes are considered in Sec. IVA, where the laboratory glass transition is analyzed in detail. Section IVB deals with heating processes. The normal curve associated with a given heating program is defined, and its relation to the observed hysteresis effect is discussed. Finally, the main conclusions of the paper are summarized in Sec. V. II. THE MODEL The model we consider has a variable number of particles n. This number completely specifies the state of the system. The empty state, n50, has zero energy, e 050, and all the states with n.0 are degenerate, with energy e n5 e . The system is in contact with a heat and particle bath characterized by a temperature Tand fugacity z [exp(2 a ). Therefore, the equilibrium probability of finding the system in state nis pn ~0!5Ce2 be ne2 a n,~2.1! where b 5(kBT)21,kBbeing Boltzmann’s constant. The constant Cis determined from the normalization condition, and it is given by C512e2 a 12e2 a ~12e2 be !.~2.2! Now, we assume that the bath is such that the equilibrium average number of particles is unity, independently of the temperature, i.e., ^ n & 05( n50 ` npn ~0!51. ~2.3! This provides a relationship between the fugacity and the temperature, namely a 5ln~11e2 be /2!.~2.4! We notice that expressions of this kind are typical when passing from a canonical description to a grand-canonical one, and the latter is required to correctly reproduce the number of particles in the system. Using this relation, Eq. ~2.2! reduces to C5exp(2 a ), and the equilibrium distribution can be written p0 ~0!5e2 a ,~2.5a! pn ~0!5e2 be 2 a ~n11!,n>1. ~2.5b! The introduction of a bath verifying Eq. ~2.4!has been stimulated by the work carried out in Refs. 9–11, where two variations of the BG model are studied. In these models, N particles can occupy Ndifferent ‘‘abacuses’’ r51,...,N. While in one of the models9,10 the particles are considered as distinguishable, in the other one11 they are treated as indistinguishable. This is the only difference between both models. Except for an additive constant, the energy of a given configuration is proportional to the number of occupied abacuses. There is no limitation in the number of particles nr being in a particular abacus r, except the one following from the total number of particles, (rnr5N. Since all the abacuses are equivalent, the average number of particles in each of them must be unity at equilibrium. The model described above mimics the equilibrium properties of the BG model with indistinguishable particles. A brief discussion of this is given in Appendix A. The idea is to focus on one of the abacuses, considering the remainder of them as a bath in the limit N→`. The condition given by Eq. ~2.4!guarantees that this limit is taken keeping the same both the number of abacuses and the number of particles. From Eq. ~2.5!it is straightforward to obtain the equilibrium properties of the system, as functions of the temperature. The average energy is ^ E & 05( n50 ` e npn ~0!5 e ~12p0 ~0!!5 e e2 be /2 11e2 be /2 ,~2.6! and its fluctuations are given by s E 25 ^ E2 & 02 ^ E & 0 25 e 2e2 be /2 ~11e2 be /2!2.~2.7! Fluctuations in the number of particles are s N 25 ^ n2 & 02 ^ n & 0 252e be /2.~2.8! 6344 55A. PRADOS, J. J. BREY, AND B. SA ´NCHEZ-REY This quantity diverges in the low-temperature limit, where most of the probability corresponds to the ground state. Due to the condition of the mean number of particles being equal to unity, the probability distribution has a long tail as a function of n. Therefore, there is an effective correlation length associated to the divergence of the fluctuations of the number of particles. Finally, the equilibrium entropy reads S kB 52 ( n50 ` pn ~0!lnpn ~0!52ln~11e2 be /2!1 be e2 be /2 11e2 be /2 . ~2.9! This expression coincides with the entropy per abacus in the BG model with indistinguishable particles. It is free of the pathological behavior shown by the entropy in the case of considering the particles as distinguishable, where it becomes negative at low temperatures.9,11 Next, we proceed to formulate the kinetics of the model. The elementary dynamical processes we will consider are the arrival or the departure of one particle, and therefore the dynamical evolution of the system will be given by a onestep process15 master equation, dpn dt 5rn11pn111gn21pn212~rn1gn!pn,~2.10! where pn(t) is the probability that the system has nparticles at time t,rnis the transition rate from state nto state n21 ~loss of one particle!, and gnis the transition rate from state nto state n11~gain of a particle!. Of course, the state n50 is a reflecting boundary, r050. ~2.11! As we do not want to introduce any energy barrier obstructing the relaxation of the system towards the ground state, we will take rn5 n ,n.0, ~2.12! where n is a constant parameter with dimensions of frequency. The transition rates gnare chosen in order to verify the detailed balance condition, i.e., g05 n e2 be 2 a ,~2.13a! gn5 n e2 a ,n.0. ~2.13b! Since the ground state can be reached at any temperature from any other state without surmounting any energy barrier, possible divergence of the characteristic relaxation time and glassy behavior can only appear in the model due to the presence of entropy barriers. In fact, glassy behavior is to be expected, because at low temperatures a →0 according to Eq. ~2.4!, and the leading behavior of the transition rates is given by gn.rn5 n ,n.0, ~2.14a! g0. n e2 be .~2.14b! One can argue where the model, as formulated here, incorporates the entropy barriers which are so evident in the original BG model. At low temperatures, the probability of finding the system in an excited state far from n51, which is the only one from which the energy can decrease, is of the same order as p1up to n5O( a 21). This reflects the equivalence of all the abacuses in the original BG model. The relaxation slows down because the random walk performed by namong the excited states contributing to the energy is symmetric, and it takes a very large time to the system to relax from states n5O( a 21). Very recently, some random walk models have been proposed to mimic the zero-temperature dynamics of the BG model.16,17 To put our work in a proper context, it is important to note that, first, we will study here the finite temperature kinetics of our model and, secondly, that we are using a grand-canonical ensemble description. For this reason our random walk is not symmetric, except in the limit T→0. Our aim is not to propose a model that exactly reproduces the dynamics of the BG model. Instead, we want to retain its main features in a solvable model, in order to identify the relevant mechanisms leading from entropy barriers to glassy behavior. The solution of the master equation in the case of time independent temperature can be obtained by using standard procedures.15 The constant n in the transition rates will be used to set up the time scale, and thus it will be taken equal to unity in the following. We look for the eigenvalues and eigenvectors of the problem. The former are given by l~q!511e2 a 22e2 a /2cosq,~2.15! where qruns in the interval @0, p #. Besides, there is the eigenvalue l50, whose eigenvector is the equilibrium distribution given by Eq. ~2.5!. All the eigenvalues lother than l50 are strictly positive, as it must be the case for a master equation with transition rates verifying detailed balance. The eigenvector associated to l(q)is j 0 ~ q ! 5 S 2 p D 1/2 e~ be 2 a !/2cos h ~q!,~2.16a! j n~q!5 S 2 p D 1/2 e2[ be 1 a ~11n!]/2cos@nq1 h ~q!#,n>1. ~2.16b! Here h (q) is a real function defined by e2i h ~q!52e2iq2 be 2 a /22e2 be 2 a 111e2 a 22e2 a /2cosq eiq2 be 2 a /22e2 be 2 a 111e2 a 22e2 a /2cosq, ~2.17! and h ~0!5 p /2. ~2.18! It has the property h (2q)52 h (q)1 p . The eigenvectors j n(q) verify the closure relation pn ~0!1 E 0 p dq j n~q! j m~q! pm ~0!5 d nm .~2.19! By using the above equation, any initial condition can be expressed as a sum over the eigenvectors, 55 6345GLASSY BEHAVIOR IN A SIMPLE MODEL WITH . . . Dn~t50![pn~t50!2pn ~0!5 E 0 p dqg~q! j n~q!, ~2.20! with g~q!5( m50 `Dm~t50! j m~q! pm ~0!.~2.21! Now, it is trivial to write the time evolution of the deviation from equilibrium Dn(t), Dn~t![pn~t!2pn ~0!5 E 0 p dqg~q! j n~q!e2tl~q!. ~2.22! This provides the general solution of the master equation, for the time independent temperature case. Since l(q) is strictly positive for all q,0<q< p ,D n(t) goes to zero in the infinite time limit, as expected. As already discussed, at zero temperature the model reduces to a symmetric random walk with an absorbing boundary at n50. Therefore, the probability distribution tends to a stationary state with pn5 d n,0 . The decay to this state is very slow and aging effects occur, even if the system was initially in equilibrium at low temperatures. Since it is easily seen that at T50 our model becomes equivalent to ‘‘model B’’ studied in detail in Ref. 17, we will not discuss the aging effects here. III. RELAXATION PROPERTIES In this section we are going to study the relaxation properties of the model at a given constant temperature. Attention will be focused on ~a!the time autocorrelation function of energy in equilibrium and ~b!the linear relaxation of energy after a temperature perturbation. It must be stressed that both quantities does not coincide, because the ensemble description of the model does not correspond to the canonical one. A. Energy time autocorrelation function The time autocorrelation function of the energy in equilibrium is given by ^ E~0!E~t! & 05( n50 ` ( m50 ` e n e mp1 u 1~n,t u m,0!pm ~0!,~3.1! where p1 u 1(n,t u m,0) is the conditional probability of finding the system in state nat time tgiven it was initially in state m. Let us introduce the response function f ~t!5 ^ E~0!E~t! & 02 ^ E & 0 2 ^ E2 & 02 ^ E & 0 2,~3.2! that verifies f ~0!51, lim t→` f ~t!50. ~3.3! The conditional probability p1 u 1(n,t u m,0) is the solution of the master equation ~2.10!with the initial condition p1 u 1~n,0 u m,0!5 d nm .~3.4! By making use of Eq. ~2.19!it is easy to see that p1 u 1~n,t u m,0!5pn ~0!1 E 0 p dq j m~q! pm ~0! j n~q!e2tl~q!, ~3.5! since pn (0) corresponds to the null eigenvalue, and j n(q)to l(q). Therefore, it is ^ E~0!E~t! & 05( n,m50 ` e n e mpn ~0!pm ~0! 1( n,m50 ` e n e m E 0 p dq j n~q! j m~q!e2tl~q! 5 ^ E & 0 21 E 0 p dqa2~q!e2tl~q!,~3.6! where we have introduced the function a~q!5( n50 ` e n j n~q!.~3.7! Substitution of Eq. ~3.6!into Eq. ~3.2!yields f ~t!5 * 0 p dqa2~q!e2tl~q! * 0 p dqa2~q!.~3.8! It follows that f (t) decays monotonically from its initial value, f (0)51, to zero. This could have been foreseen, since it is a general property for equilibrium autocorrelation functions in models whose dynamics is described by means of master equations with the transition rates verifying the detailed balance condition. The problem has been reduced to calculate the function a(q), defined by Eq. ~3.7!, that can be written as a~q!5( n51 ` e j n~q!52 e j 0~q!,~3.9! because ( n50 ` j n~q!50, ~3.10! due to the orthogonality of the eigenvectors j (q) with respect to the equilibrium distribution. From Eqs. ~3.9!and ~2.16a!we obtain a~q!}cos h ~q!.~3.11! The proportionality constant in the above relation is irrelevant for the calculation of the response function, given by Eq. ~3.8!. The function h (q) defined in Eq. ~2.17!is rather involved, but simple expressions are derived both in the limits of short and long times. For short times, t!1, it is f ~t!;e2tlM,~3.12! where lM[ t S 215 * 0 p dql~q!a2~q! * 0 p dqa2~q!5e a 21. ~3.13! 6346 55A. PRADOS, J. J. BREY, AND B. SA ´NCHEZ-REY Thus, the relaxation in the short time regime is exponential, as it is the usual case in systems described by master equations.18 In the limit of long times, a Laplace’s analysis of Eq. ~3.8!gives f ~t!;e a 21 2 p 1/2e9 a /4~11e2 a /22e2 a !2 e2t~12e2 a /2!2 ~12e2 a /2!4t3/2 . ~3.14! Aside from slow algebraic corrections, the relaxation is again exponential, but with a characteristic time t L5~12e2 a /2!2,~3.15! which is different from the one of the short time regime. This is also the most common case in models described by master equations. This fact, together with the monotonic decay of the equilibrium autocorrelation function, leads to a nonexponential relaxation regime at intermediate times.18 This regime is expected to be more relevant as the time scales separation becomes larger. This is the case when e a →1. Then, both characteristic times diverge, but t L@ t S.~3.16! Taking into account Eq. ~2.4!for a , it follows that e a →1 is equivalent to b →`or T→0. In this limit, both Eqs. ~3.12!and ~3.14!become much simpler. For short times it is f ~t!;e2 a t,~3.17! whereas in the long time region f ~t!;1 p 1/2 e2 a 2t/4 ~ a 2t/4!3/2 .~3.18! This latter equation shows that relaxation takes place over a time scale s5 a 2t 4,~3.19! which is much longer than the defined by the initial exponential. Thus, separation of time scales comes up, and nonexponential relaxation is to be expected in an intermediate time window. The picture we have obtained is similar to the one found in the low-temperature relaxation of Glauber’s Ising model.19–21 Therefore, we make use of the same techniques to derive the behavior of the correlation function in the intermediate time regime in the low-temperature limit. To begin with, we obtain an expression which is valid in the time scale defined by Eq. ~3.19!. We introduce a new variable u through q5 a u/2. ~3.20! Then, a simple analysis gives f ~t![ f ¯ ~s!54 p E 0 `du u2 ~11u2!2e2s~11u2!,~3.21! where terms of order a have been neglected. For very long times, s@1, Eq. ~3.18!is of course recovered. However, in the region s!1 we do not get the low-temperature version of the short time behavior, as given by Eq. ~3.17!, since it corresponds to the much shorter time scale defined by a 21. Over the time scale s, that behavior collapses onto the point s50. Actually, s!1 corresponds to an intermediate time window where tis large but s5 a 2t/4 is small ( a !1). It is easy to see that ln f ¯ ~s!;24 p 1/2 s1/2,s!1. ~3.22! Therefore, from Eqs. ~3.19!and ~3.21!we get ln f ~t!;2 S 4 a 2t p D 1/2 ,~3.23! which is a stretched exponential or Kohlrausch-WilliamWatts ~KWW!function, ln f ~t!52 S t t D g ,~3.24! with g 51/2, ~3.25a! t 5 p 4 a 2; p 4e be .~3.25b! Thus, at low temperatures the relaxation time t obeys the Arrhenius law, with an ‘‘activation’’ energy e . One may ask himself which is the physical origin of this behavior, since the system does not have to surmount any energy barrier to reach the ground state. In our model, as in the one proposed by Ritort,9there is an entropy barrier. At low temperatures a →0 and a symmetric random walk is performed by the system among all the excited states. The characteristic relaxation time will be dominated by the diffusion process from the mean position in the excited region to the state n50. The mean position in the excited states n ¯ exc is given by n ¯ exc5(n51 `npn ~0! (n51 `pn ~0!5 ^ n & 0 12p0 ~0!5~12e2 a !21,~3.26! where we have made use of Eq. ~2.5a!. This quantity must not be confused with the average number of particles in excited states. In the limit of low temperatures n ¯ exc; a 21@1. ~3.27! Then, an estimation to the time needed to diffuse until n50 will be t dif5O~n ¯ exc 2!5O~ a 22!.~3.28! The above equation can be considered as a qualitative explanation of the relaxation time t dependence on the temperature shown by Eq. ~3.25b!, since it is reasonable to expect that t 5O( t dif), the mean time taken by the system to get to the ‘‘bottleneck’’ in the configuration space. A simplified picture of the evolution of the equilibrium time autocorrelation function f (t) of the energy can be given in terms of the three time regimes we have found, 55 6347GLASSY BEHAVIOR IN A SIMPLE MODEL WITH . . . f ~t!5 H e2 a tt!1, e2~4 a 2t/ p !1/2 1!t!4 a 22, p 21/2~ a 2t/4!23/2e2 a 2t/4 t@4 a 22. ~3.29! A similar behavior has been previously obtained for relaxation in different models.18,19,22 It must be noticed that the scheme described by Eq. ~3.29!is consistent with both empirical and numerical results for glassy systems, where nonexponential relaxation and KWW behavior is usually found over an intermediate time window.23 It is possible to estimate roughly the range of validity of the KWW function. One can determine the time intersections tiand tfof the KWW function with the short and long time exponentials, respectively. It is found that ti54/ p and a 2tf.6.28. In the time interval (ti,tf) the KWW function is expected to hold, and the relaxation function verifies exp(24 a / p )> f (t)>0.06. Although this is a very crude estimation, we conclude that most of the relevant part of the relaxation of f (t) at low temperatures is given by the stretched exponential in Eq. ~3.29!, because a !1. In Fig. 1 we have plotted f (t) for be 510, which corresponds to a 56.731023. The solid line is the KWW function given by Eq. ~3.23!. As discussed in the paragraph above, it is valid over an intermediate time window corresponding to the relevant part of the relaxation. For very long times, relaxation is exponential, and the KWW function is not a good approximation. For very short times, relaxation is also exponential, but the difference with the KWW function is negligible over the scale of the figure. Finally, it must be remarked once more that the KWW decay found at intermediate times follows from the existence of two exponential regimes valid at very short and very long times with a clear separation of their respective time scales. A detailed discussion can be found in Ref. 18 for any system whose dynamics is described by a master equation. The main point is whether most of the relevant part of the relaxation can be described by a KWW function as a consequence of a clear time scale separation. This happens in our model because the relaxation spectrum becomes very broad at low temperatures. This is not a general property for all master equations, and KWW relaxation may not show up for a given choice of the transition rates, if the relaxation spectrum associated to them remains narrow at low temperatures. B. Linear relaxation of the energy The energy relaxation after a temperature perturbation is characterized by the response function c ~t!5 ^ E~t! & 2 ^ E & 0 ^ E~0! & 2 ^ E & 0,~3.30! where ^ E~t! & 5( n50 ` e npn~t!5 e @12p0~t!#.~3.31! Using the definition of Dnin Eq. ~2.22!, we have c ~t!5D0~t! D0~0!.~3.32! Substitution of the exact solution of the master equation for constant temperature obtained in Sec. II, Eqs. ~2.21!and ~2.22!, leads to c ~t!5 * 0 p dqg~q!cos h ~q!e2tl~q! * 0 p dqg~q!cos h ~q!.~3.33! We have to calculate g(q), from the initial conditions Dn(0). We will consider that the system was in equilibrium at a temperature b 1D b . Then, the temperature was instantaneously changed to b at t50. In the linear response approximation, Dn~0!5pn ~0!~ b 1D b !2pn ~0!~ b !5dpn ~0! d b D b ,~3.34! and the function g(q) in Eq. ~2.21!reads FIG. 2. Energy relaxation in the low-temperature region, for a temperature value corresponding to e /kBT510. The diamonds are the numerical evaluation of Eq. ~3.33!, while the solid line corresponds to the KWW function of Eq. ~3.43!. In this logarithmic scale, we have restricted ourselves to positive values of the response function. FIG. 1. Plot of the equilibrium autocorrelation function of energy, for a temperature value corresponding to e /kBT510. The diamonds are the numerical evaluation of Eq. ~3.8!, and the solid line is the stretched exponential of Eq. ~3.29!. 6348 55 A. PRADOS, J. J. BREY, AND B. SA ´NCHEZ-REY g~q!5D b ( n50 ` j n~q!d d b lnpn ~0!.~3.35! The expression of the equilibrium distribution, Eq. ~2.5!, is equivalent to lnpn ~0!52 be ~12 d n0!2 a ~n11!,;n>0, ~3.36! and substitution of Eqs. ~2.16!and ~3.36!into Eq. ~3.35!, together with the relation d a d b 52 e 2 e2 be /2 11e2 be /2 52 e 2e2 a 2 be /2,~3.37! leads, after some algebra, to g~q!5 S 2 p D 1/2 e D b F e~ be 2 a !/2cos h ~q!11 2e2 be 22 a cos@q1 h ~q!#22e2 a /2cos h ~q!1e2 a cos@ h ~q!2q# ~11e2 a 22e2 a /2cosq!2 G .~3.38! The above expression for g(q) is rather involved for arbitrary temperature. In the low-temperature limit, a →0, introducing again the time scale sdefined by Eq. ~3.19!and the variable uof Eq. ~3.20!, one gets c ~t![ c ¯ ~s!58 p E 0 `duu2~u221! ~11u2!3e2s~11u2!.~3.39! The relaxation of the energy takes place over a time scale of order a 22, as it was the case of the equilibrium energy autocorrelation. In the stime scale, the initial exponential relaxation does not show up, because the short time behavior of c (t) is given by c ~t!.e22tsinh a ,~3.40! and at low temperatures its characteristic time scale (2 a )21collapses onto the point s50. For very long times, s@1, a Laplace analysis of Eq. ~3.39!yields c ¯ ~s!;22 p 1/2 e2s s3/2 .~3.41! The previous equation tells us that energy relaxation is not monotonic. In fact, it is proved in Appendix B that E 0 `dt c ~t!50, ~3.42! implying that c (t) is negative in a time region. However, in the intermediate time window s!1 a stretched exponential decay is again obtained, though the general argument developed in Ref. 18 cannot be directly applied. For s!1, it is easy to show from Eq. ~3.39!that ln c ~t!;2 S 16 a 2t p D 1/2 .~3.43! Therefore, a simplified picture of the energy relaxation at low temperatures is obtained, which is similar to the one found before for the energy autocorrelation. In terms of the three relevant time regimes that have arisen in our discussion, c ~t!5 H e22 a tt!1, e2~16 a 2t/ p !1/2 1!t!4 a 22 22 p 21/2~ a 2t/4!23/2e2 a 2t/4 t@4 a 22. ~3.44! At very long times, the relaxation function c (t) crosses the taxis and decays to zero from negative values. This is quite a small effect, since a numerical estimation of the minimum of c (t) gives c min.20.05. Therefore, the KWW function in Eq. ~3.44!also gives a relevant information about the energy relaxation at low temperatures. In particular, its relaxation time t E5 p 16 a 22; p 16e be ,~3.45! can be used to characterize the relaxation of energy after a homogenous perturbation in temperature. It must be remarked that t Ealso follows an Arrhenius law at low temperatures. A qualitative explanation of this behavior, in terms of the diffusive motion of the system, can be given along the same way as in the previous section. Obviously, the stretched exponential approximation is not able to explain the crossing of the taxis that takes place at very long times, but it accurately fits most of the relevant part of energy relaxation, namely up to c .0.1. In Fig. 2 the energy relaxation function obtained numerically is compared with the KWW function in Eq. ~3.44!. The value of the parameter is the same as in Fig. 1, i.e., be 510 ( a 56.731023). In the variables used in Fig. 2, exponential relaxation corresponds to a straight line of unity slope, while KWW relaxation is represented by a straight line of slope equal to the parameter g in Eq. ~3.24!. The logarithm scale used amplifies the discrepancies, especially for short times, where the difference between the KWW function and the initial exponential is in fact negligible. IV. THERMAL CYCLES Here we are interested in studying the behavior of the model when it is continuously cooled down from high to low temperatures, and afterwards reheated. This is usually called a thermal cycle. Upon describing it, the system may deviate from equilibrium while being cooled, leading to the kinetic phenomenon known as the laboratory glass transition. In the heating process, equilibrium is approached again at high 55 6349GLASSY BEHAVIOR IN A SIMPLE MODEL WITH . . . temperatures, but the system follows a different curve from the cooling one, and hysteresis shows up. The kinetic behavior just discussed is shown by a wide class of materials,2,4 and also by some simple models.9,11–13,21,24 Nevertheless, analytical results are scarce,13,21 although quite a general explanation of the hysteresis phenomenon has been given.12 It can be understood as the monotonic approach to a ‘‘normal’’ curve, different from the equilibrium one, characterizing heating processes. As the proof in Ref. 12 was made for the canonical ensemble, a generalization for the case considered here is presented in Sec. IVB. The remainder of this section is organized as follows. First, we study cooling processes, and the existence of the laboratory glass transition. Secondly, heating processes are considered, paying special attention to the appearance of hysteresis, and relating it to the trend of the system to approach the normal curve. Let us point out that we have not been able to solve exactly the master equation for the case of time-dependent temperature, that implies that the transition rates are also time dependent. The procedure developed in Ref. 13 is valid when the eigenvectors of the master equation do not depend on temperature. This is not the case here, since the eigenvectors of the master equation, given by Eq. ~2.16!, are temperature dependent through the function h (q) in Eq. ~2.17!. Therefore, we have performed a Monte Carlo simulation of the master equation, using a generalization of the Bortz-KalosLebowitz algorithm25 for master equations with timedependent transition rates.26 Nevertheless, some analytical estimations can be done, and they will be compared with numerical results. A. Cooling processes and laboratory glass transition Now we are going to study the continuous cooling of the system to low temperatures. In order to analyze the deviation from equilibrium values of the properties of the system, we will follow a reasoning similar to that used in Ref. 13. We start from the relaxation spectrum of the master equation, Eq. ~2.15!, and notice that the modes ldepend on temperature through a , and therefore they are time dependent in a given cooling program T(t). In this spectrum, the relaxation rates of the system vary with their label q, from the minimum value, corresponding to q50, l1511e2 a 22e2 a /25~12e2 a /2!2,~4.1! to the maximum one, for q5 p , l2511e2 a 12e2 a /25~11e2 a /2!2.~4.2! Given a cooling law, to each of the relaxation modes we can associate a characteristic time scale s~q!5 E t t0dt8l~q;T8!,~4.3! where t0is the extrapolated time for which the temperature would vanish according to the cooling program, and T8[T(t8). The time s(q) is roughly proportional to the effective number of transitions left to the mode l(q;T) before reaching T50. For times longer than the one t(q) making s(q)51, one can consider that the mode will not experiment any more transitions. Thus, for temperatures lower than the one corresponding to t(q) the contribution of the mode will not evolve in time and can be considered as ‘‘frozen.’’ In this way, we can determine a freezing temperature T(q) for each value of q. Equivalently, one can introduce the notion of a ‘‘demarcation’’ mode qD(T), such that modes with q<qD(T) are frozen, while modes with q.qD(T) are still relaxing at the given temperature.13,27 The laboratory glass transition begins at the temperature T1[T(t1) given by the relation E t1 t0dt8l1~T8!51~4.4! or qD~T1!50, ~4.5! i.e., only the slowest relaxation rate is frozen, and the deviation from equilibrium starts off. On the other hand, the system will be completely frozen at a temperature T2[T(t2) for which the fastest relaxation mode does not evolve any more, namely, E t2 t0dt8l2~T8!51~4.6! or qD~T2!5 p .~4.7! A global image of the freezing phenomenon can be obtained by means of the time scale s5 E t t0dt81 t ~T8!,~4.8! where t (T) is the time characterizing the relaxation of the property Pwe are interested in after a temperature perturbation. For instance, in our model t would be the KWW relaxation time t Ein Eq. ~3.45!, if we want to describe the energy evolution during the cooling process. An estimation of the ‘‘global’’ freezing temperature Tffor the property Pis obtained by making s51, i.e., 15 E tf t0dt 1 t ~T!,~4.9! and then Tf5T(tf). Since the laboratory glass transition is very narrow in temperature, at least when the system is slowly cooled, an approximation to the frozen value of the property Punder consideration would be P0(Tf), i.e., the equilibrium value at its freezing temperature Tf. It is important to note that the temperatures T1,T2, and Tfdepend both on the cooling rate rcand the cooling law f(T) defining the cooling program, dT dt 52rcf~T!.~4.10! This is also the case in other simple models whose dynamics is described in terms of master equations. For some choices of the cooling law f(T), the system remains in equilibrium at all temperatures.6,13 We are not going to discuss this problem 6350 55A. PRADOS, J. J. BREY, AND B. SA ´NCHEZ-REY here, but focus our attention on the behavior of the system when it is being linearly cooled, dT dt 52rc,~4.11! i.e., f(T)51, which is the most usual cooling program in real experiments4and also in theoretical studies of model systems.9,11,24,28 The relaxation modes l1and l2, Eqs. ~4.1!and ~4.2!, and the time characterizing the energy relaxation t Eare written as functions of a , defined by Eq. ~2.4!. Then, it is useful to transform the time integral in the definition of the sscales into an integral over a with the aid of d a dt 52 r c~12e2 a !@ln~e a 21!#2,~4.12! where Eq. ~4.11!has been taken into account, and r c52kBrc e ~4.13! is an adimensional cooling rate, giving the time scale over which a evolves. As discussed above, the beginning of the laboratory glass transition is estimated to take place at a time t1such that T(t1)5T1, being T1the temperature in Eq. ~4.5!, i.e., the one at which the slowest relaxation mode freezes. By using Eqs. ~4.1!and ~4.12!, we can write 151 r c E 0 a 1d a ~12e2 a /2!2 ~12e2 a !@ln~e a 21!#2,~4.14! where a 1[ a (t1). In the limit of slow cooling, r c!1, and it follows that a 1!1. For this case, Eq. ~4.14!simplifies to 151 4 r c E 0 a 1d aa ~ln a !2.~4.15! To solve this relation for a 1, we make the change of variable a 5 a 1x, E 0 a 1d aa ~ln a !25 a 1 2 ~ln a 1!2 E 0 1dx x @11~lnx/ln a 1!#2 ; a 1 2 2~ln a 1!2.~4.16! The last integral can be done by dividing the interval (0,1) into the two subintervals (0, u ln a 1 u 21) and ( u ln a 1 u 21,1). In the first interval, the integrand is bounded by unity, and the integral is negligible. In the second interval, it is u lnx u ! u ln a 1 u , giving rise to the result in Eq. ~4.16!. Substitution into Eq. ~4.15!yields 151 8 r c a 1 2 ~ln a 1!2.~4.17! By making use of the slow cooling condition, a 1!1, we have 2ln a 1;ln~8 r c!.~4.18! Now, we take into account that a 1;exp(2 b 1 e /2), to get T1; e kB 1 u ln~8 r c! u .~4.19! In order to calculate the fictive temperature Tf, we start from Eq. ~4.9!, with the relaxation time of energy t Egiven by Eq. ~3.45!, 15 E 0 a fd a dt d at E 21~ a !516 pr c E 0 a fd aa 2 ~e a 21!@ln~e a 21!#2. ~4.20! As before, a fis the value of a corresponding to Tf. For slow cooling, it is a f!1, since a f, a 1. Then, the above equation reduces to 1516 pr c E 0 a fd aa ~ln a !2.~4.21! In this way, we have arrived at an expression similar to Eq. ~4.15!for a 1. Therefore, E 0 a fd aa ~ln a !2; a f 2 2~ln a f!2~4.22! and 158 pr c a f 2 ~ln a f!2.~4.23! Again, a reasoning along the line of the one above Eq. ~4.19! gives us Tf; e kB 1 u ln~ pr c/8! u .~4.24! A similar dependence on the cooling rate is obtained in real experiments,4and has also been found in Glauber’s Ising model.13 Taking into account the comment below Eq. ~4.9! one can estimate the residual value of the energy, i.e., er[lim T→0 @ ^ E & ~T!2 ^ E & 0~T!#5 ^ E & 0~Tf!.~4.25! In the limit of slow cooling, r c!1, Eq. ~4.25!leads to a potential dependence on the cooling rate of the residual energy, er} r c 1/2 . A similar behavior of the residual properties has been obtained in some models of glasses.5,13,6 In Fig. 3 the evolution of the mean energy for the cooling program in Eq. ~4.11!, with an adimensional cooling rate r c50.02, is plotted. The departure from equilibrium roughly begins at the temperature obtained from Eq. ~4.19!, kBT1/ e .0.55. The estimation of the freezing temperature, obtained by using Eq. ~4.24!is kBTf/ e .0.21, in good agreement with the numerical result. The frozen value of the energy given by the Monte Carlo simulation is ^ E & / e 50.076, while the value obtained from Tfis ^ E & 0(Tf)/ e 50.081. Again, the approximated theory provides a reasonable estimation of the actual value. 55 6351GLASSY BEHAVIOR IN A SIMPLE MODEL WITH . . .