Glasslike dynamical behavior in hierarchical models submitted to continuous cooling and heating processes

Glasslike dynamical behavior in hierarchical models submitted to continuous cooling and heating processes 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 共Received 15 March 2001; published 24 September 2001兲 The dynamical behavior of a kind of models with hierarchically constrained dynamics is investigated. The models exhibit many properties resembling real structural glasses. In particular, we focus on the study of time-dependent temperature processes. In cooling processes, a phenomenon analogous to the laboratory glass transition appears. The residual properties are analytically evaluated, and the concept of fictive temperature is discussed on a physical basis. The evolution of the system in heating processes is governed by the existence of a normal solution of the evolution equations, which is approached by all the other solutions. This trend of the system is directly related to the glassy hysteresis effects shown by these systems. The existence of the normal solution is not restricted to the linear regime around equilibrium, but it is defined for any arbitrary, far-fromequilibrium, situation. DOI: 10.1103/PhysRevE.64.041505 PACS number共s兲: 64.70.Pf, 05.70.Ln, 45.70.⫺n I. INTRODUCTION In recent years, there has been quite a considerable amount of work in models in which glassy behavior is generated not by quenched disorder, but by kinetic constraints. The kinetic restrictions are responsible for the slow relaxation, since the state of a particle or a group of particles can only change if some condition of its environment is fulfilled. In particular, ‘‘facilitated’’ models have been considered, both for structural glasses 关1–5兴and for granular systems 关6兴. The characteristic feature of facilitated models is that a particle 共spin兲can only change its state if a certain number of its neighbors is in an excited state. Also, hierarchically constrained models have been used to study stretched exponential relaxation in glasses 关7兴. In these models, the system is structured in levels and a particle in a given level can only make a transition if a given cluster of particles in the lower level is in a certain subset of configurations. Then, the dynamics of the several levels are coupled, and the characteristic relaxation times increase with the level index. Hierarchically constrained dynamics may be relevant for those complex systems in which the time evolution of the slowest modes is controlled by the relaxation of the fastest ones. This qualitative picture is adequate to describe, among other problems, proteins relaxation 关8兴and the densification of powders and structural glasses at high pressure 关9,10兴. Very recently, a kind of hierarchically constrained dynamics has been shown to exhibit, in quite a natural way, logarithmic relaxation 关11兴. This kind of ‘‘anomalous,’’ highly nonexponential, decay is observed in a wide variety of complex systems, including spin glasses 关12,13兴, granular materials 关9,10,14,15兴, structural glasses 关16–18兴, and protein models 关8,19兴. The aim of this work is to study the dynamical behavior of a general class of hierarchically constrained models when submitted to more complicated processes. In particular, we are interested in the behavior of a system with hierarchical constraints when the temperature changes in time, which makes the coupling between the levels time dependent. The consideration of time-dependent temperature processes requires an extension of the original model as formulated by Palmer et al. 关7兴. This will be done in a very simple, but natural, way: the coupling between the levels varies in time because the probability of a cluster configuration allowing relaxation of a particle in the next level depends on the temperature. On the other hand, neither the number of particles in a given level nor the length of the clusters ‘‘facilitating’’ the relaxation depends on the temperature. They are considered as quantities defined in the coarse-grained description of the system introduced to model the physical problem at hand. Let us note that hierarchical models can also be applied to the analysis of nonthermal systems, such as granular materials in the dense regime. For those materials, thermal energy is not enough to make the system explore the phase space of configurations. Then, the system must be externally excited — for instance, vibrating it — in order to be able to evolve. In these situations, the role of temperature is played by the intensity of the external driving. If the stationary state reached by the system in the long-time limit can be described by Edward’s theory 关20,21兴, the compactivity X, which is analogy of the temperature in thermal systems, will be a function of the intensity of the external force. Then, by exploiting the analogies of Edward’s theory, i.e., substituting volume by energy and compactivity by temperature, it is possible to incorporate nonthermal systems in our formulation. Processes in which the temperature is time dependent are physically relevant because they can be used to study some characteristic dynamical aspects of glasses. For instance, when a supercooled liquid is cooled down to very low temperatures, a laboratory glass transition is observed. A dramatic change in the behavior of the system takes place, and it departs from the equilibrium curve, getting ‘‘frozen’’ in a far-from-equilibrium state. This transition appears as a consequence of the fast increase of the relaxation time with decreasing temperature. In order to characterize the cooling process, experimental physicists often use the residual value of the relevant physical properties, i.e., the difference between their actual values over the cooling curve and the value obtained by extrapolation of the equilibrium curve to PHYSICAL REVIEW E, VOLUME 64, 041505 1063-651X/2001/64共4兲/041505共15兲/$20.00 ©2001 The American Physical Society64 041505-1 the very-low-temperature region 关22,23兴. If the system is reheated from the nonequilibrium state, it returns to equilibrium for high temperatures, but it follows a different curve from the cooling one, giving rise to hysteresis effects. This phenomenon is related to the ‘‘nonlinearity’’ of glassy relaxation: the approach towards the equilibrium curve depends on the configuration of the system, measured by the so-called fictive temperature 关22–24兴. Narayanaswami’s theory provides a phenomenological explanation of this behavior 关22– 25兴. Interestingly, a similar behavior has been found in vibrated granular materials when the tapping intensity is varied in a cyclic way 关26兴, although the hysteresis effects are more evident when the heating process begins in a loosely packed state, referred as to the ‘‘irreversible’’ branch in the experiments. We will start from a very general hierarchical spin model, in which pseudospins are organized into levels, labeled by an index n. The pseudospins are assumed to correspond to some coarse-grained description of the system. They can take only two values, representing, for instance, two possible densities of a certain small subvolume of the system. One pseudospin in level n⫹1 can only flip between the two possible values if a cluster of ␮ nspins in level nis in a given subset of configurations. This is the basic characteristic of hierarchically constrained models as introduced in Ref. 关7兴, and it slows down the relaxation in level n⫹1, as compared with that of level n. Here we will consider the simple choice that all the spins in the cluster must be in the up 共excited兲state. The exact dynamical equation for the evolution of the pseudospins involves very complicated moments of the probability distribution. To get an exactly solvable model, a ‘‘meanfield’’ approximation will be introduced. Then, the characteristic relaxation time ␶ nof level nis seen to increase both with the index label n, due to the hierarchical constraint, and also with decreasing temperature, since the configurations allowing the system to relax become less probable when the temperature is lowered. Some exact dynamical results for systems described with master equations with time-dependent transition rates are known. In particular, the existence of a ‘‘normal’’ solution, i.e., a solution of the master equation that is approached by all the others, has been proved on a very general basis 关27兴. The main required conditions are the irreducibility of the Markov process for long enough times and that the transition rates be externally controlled, so that they do not depend on the probability distribution of the system. Moreover, it has been established that the normal solution tends to the equilibrium curve for very high temperatures in continuous heating processes 关27兴. The linear correction of the normal solution with respect to the equilibrium curve has been computed using Hilbert’s method 关28兴. It is not evident whether these results still hold when approximations are introduced in the dynamics of the system. In particular, in ‘‘mean-field’’-type approximations, the demonstrations in Refs. 关27,28兴are not valid, since the transition rates become functionals of the probability distribution function. Nevertheless, we will show that all the above mentioned properties apply to our simplified model. This is a good test of the plausibility of the approximations carried out and perhaps an indication that the results obtained here are more general than the derivations in Refs. 关27,28兴. The organization of the paper is as follows. In Sec. II the hierarchical model is introduced, and the exact evolution equation for the average spin is obtained. By introducing a mean-field approximation, this equation can be closed. Afterwards, a specific, but quite general, choice for the functions defining the model is made. This allows us to introduce a continuous limit in which the relaxation of the system at constant temperature is solved in Sec. III. Time-dependent temperature processes are considered in Sec. IV, where the general solution for the evolution of the probability distribution and the average energy are obtained. The general solution is similar to the expression proposed by Narayanaswami on a phenomenological basis 关22,24,25兴. Section IV is devoted to the analysis of Hilbert’s expansion, which is valid in the very-high-temperatures regime. Cooling processes are addressed in Sec. V where, for the sake of simplicity, a concrete cooling law is studied, for which the residual properties are analytically calculated. A qualitative analysis of the glasslike transition is presented in Sec. V. It allows us to give very good estimates of the residual properties and leads to the introduction of the concept of fictive temperature in a very natural way. The behavior of the system when it is reheated from low temperatures is considered in Sec. VI. The main role played by the normal solution for the understanding of the hysteresis effects shows up. Moreover, the analysis clearly indicates that the relevance of the normal solution is not restricted to near equilibrium situations. Finally, a discussion of the main points in this work is given in Sec. VII. II. DYNAMICS OF HIERARCHICALLY CONSTRAINED MODELS In this section a general kind of spin model with hierarchically constrained dynamics will be introduced. We will focus on the evolution of the average value of the spin, which is supposed to be the relevant variable. For instance, in a thermal system it will be directly related to the mean energy. Then, let us consider a system whose degrees of freedom can be classified into levels, labeled by an index n ⫽0,1,2,...,nmax . The degrees of freedom in level nwill be represented by Nnpseudospins, ␴ i (n)⫽⫾1, i⫽1,2,...,Nn. The Hamiltonian of the system is assumed to have the form H⫽h兺 n⫽0 nmax 兺 i⫽1 Nn mi (n),mi (n)⫽1⫹ ␴ i (n) 2.共2.1兲 Note that mi (n)is the occupation number of the ‘‘up’’ 共⫹1兲 state of the corresponding site. In order to write Eq. 共2.1兲we have supposed that there is no interaction between the pseudospins, but there is an ‘‘external field’’ h. For a thermal system, Hgives the energy of a given microstate of the system, while for a nonthermal system, like a powder, it could be interpreted as the volume of a given, mechanically stable, configuration of ‘‘grains’’ 关20,21兴. Using the terminology for thermal systems, the average value of the dimensionless energy per spin over the ensemble of systems considered is A. PRADOS AND J. J. BREY PHYSICAL REVIEW E 64 041505 041505-2 ␧⫽ 具 H 典 Nh ⫽1 N兺 n⫽0 nmax 兺 i⫽1 Nn pi (n),共2.2兲 where pi (n)⫽ 具 mi (n) 典 ⫽1⫹ 具 ␴ i (n) 典 2共2.3兲 is the probability that the ith spin of level be in the up state, N⫽兺n⫽0 nmaxNnis the total number of pseudospins, and the angular brackets denote statistical ensemble average. The system is considered to be in contact with a heat bath at temperature T, so that the equilibrium average value of the pseudospins does not depend either on inor on n, and it is given by 具 ␴ 典 e⬅ 具 ␴ i (n) 典 e⫽⫺tanh 冉 1 T* 冊 .共2.4兲 Here T*is a dimensionless temperature, and T*⫽2kBT/h, kBbeing Boltzmann’s constant. For the sake of concision, we will drop the asterisk in the following. The above average value of spin follows from the equilibrium probability for the ‘‘up’’ state of any spin, pe⬅pi,e (n)⫽e⫺1/T e⫺1/T⫹e1/T.共2.5兲 In the limit of infinite temperature or zero external field, both states of the pseudospins are equiprobable, pe⫽1/2, and 具 ␴ 典 e⫽0. From Eq. 共2.2兲, it follows that the equilibrium value of ␧is ␧e⫽pe.共2.6兲 For granular materials, the role of the temperature Tis played by the compactivity 关20,21兴, which is linked to the intensity of the perturbation, allowing the system to explore the configuration space. The dynamics of the model is formulated by means of a master equation with single-spin-flip Glauber transition rates 关29兴. Let us consider the flip of a given spin ␴ i (n). This transition connects a given configuration ␴ of the whole system with the configuration Ri (n) ␴ , where Ri (n)is the operator which rotates the spin ␴ i (n), keeping all the other spins the same. The transition rate for the flip of the spin ␴ i (n)in configuration ␴ is Wi (n)共 ␴ 兲⫽1 2 ␣ i (n)共 ␴ 兲 冋 1⫹ ␴ i (n)tanh 冉 1 T 冊 册 .共2.7兲 The characteristic relaxation rate ␣ i (n)of the spin ␴ i (n)depends on the configuration ␴ of the system through the hierarchical constraint ␣ i (n)共 ␴ 兲⫽关 ␣ ki (n⫺1)共 ␴ 兲 ␣ ki⫹1 (n⫺1)共 ␴ 兲••• ␣ ki⫹ ␮ n⫺1⫺1 (n⫺1) ⫻共 ␴ 兲兴1/ ␮ n⫺1 兿 j⫽ki ki⫹ ␮ n⫺1⫺1 ␦ ␴ j (n⫺1) ,⫹1,共2.8兲 where ␦ ij is the Kronecker delta. This expression implies that the spin ␴ i (n)needs, in order to flip, that all the spins belonging to a cluster of length ␮ n⫺1starting at a given spin kiof level n⫺1 must be in the up 共⫹1兲state. Besides, the characteristic flip rate of the spin is the average of the characteristic flip rates of the spins belonging to the cluster determining its possibility of change. This restricting condition is schematically depicted in Fig. 1. Note that the possibility of a given spin in level nto flip is restricted by the state of a set of clusters in all levels n⬘⬍n, the number of clusters involved in each level increasing as n⬘decreases. The hierarchical constraint implies that the configuration with all the spins in the down (⫺1) state is completely absorbent; i.e., the system does not evolve in time from that configuration. We are interested in the time evolution of the average spin ␴ i (n), which is given in Glauber dynamics by 关29兴 d dt 具 ␴ i (n) 典 ⫽⫺2 具 ␴ i (n)Wi (n)共 ␴ 兲 典 ,共2.9兲 and substitution of Eqs. 共2.7兲and 共2.8兲into this expression yields d dt 具 ␴ i (n) 典 ⫽⫺ 冓 关 ␣ ki (n⫺1)共 ␴ 兲 ␣ ki⫹1 (n⫺1)共 ␴ 兲••• ⫻ ␣ ki⫹ ␮ n⫺1⫺1 (n⫺1) 共 ␴ 兲兴1/ ␮ n⫺1共 ␴ i (n)⫺ 具 ␴ 典 e兲 ⫻ 兿 j⫽ki ki⫹ ␮ n⫺1⫺1 ␦ ␴ j (n⫺1) ,⫹1 冔 ,共2.10兲 where we have taken into account that ( ␴ i (n))2⫽⫹1 for all i, n. This equation is rather involved, since it couples the evolution of 具 ␴ i (n) 典 to moments of the probability distribution containing an increasing number of spins of all the levels n⬘ such that 0⭐n⬘⬍n. The levels 0⭐n⬘⭐n⫺2 enter into the equation through the rates ␣ j (n⫺1)( ␴ ). Then, we introduce at this stage a sort of ‘‘mean-field’’ approximation for the tranFIG. 1. The framed spin in level n, ␴ i (n), has a nonvanishing probability of changing its state only if the framed cluster of spins in level n⫺1 are in the state shown in the figure, i.e. all of them up 共⫹1兲. In this example, we have taken ki⫽i⫺2 and ␮ n⫺1⫽5. GLASSLIKE DYNAMICAL BEHAVIOR IN . . . PHYSICAL REVIEW E 64 041505 041505-3 sition rates, which is closely related, in spirit, to the seminal work of Palmer et al. 关7兴. Within this coarsening dynamics, it is reasonable to expect that spins in level nevolve over a time scale quite larger than that characteristic of level n⫺1. This means that spins in level n⫺1 change many times their state before a transition in level ntakes place. Thus, we replace the product of the Kronecker deltas in Eq. 共2.8兲by its average value, i.e., by the probability P(n⫺1)( ␴ ki (n⫺1) ⫽⫹1,..., ␴ ki⫹ ␮ n⫺1⫺1 (n⫺1) ⫽⫹1) that all ␮ n⫺1spins of the given cluster are in the up state, ␣ i (n)共 ␴ 兲⫽关 ␣ ki (n⫺1)共 ␴ 兲 ␣ ki⫹1 (n⫺1)共 ␴ 兲••• ␣ ki⫹ ␮ n⫺1⫺1 (n⫺1) 共 ␴ 兲兴1/ ␮ n⫺1 ⫻P(n⫺1)共 ␴ ki (n⫺1) ⫽⫹1,..., ␴ ki⫹ ␮ n⫺1⫺1 (n⫺1) ⫽⫹1兲.共2.11兲 Moreover, we will restrict ourselves to situations where there is spatial homogeneity within each of the levels, so that the dependence on the specific site considered in a given level can be dropped, obtaining ␣ (n)共 ␴ 兲⫽ ␣ (n⫺1)共 ␴ 兲P(n⫺1)共 ␴ 1 (n⫺1)⫽⫹1,..., ␴ ␮ n⫺1 (n⫺1)⫽ ⫹1兲.共2.12兲 Iteration of the above relation gives ␣ (n)共 ␴ 兲⫽ ␣ (0) 兿 j⫽0 n⫺1 P(j)共 ␴ 1 (j)⫽⫹1,..., ␴ ␮ j (j)⫽⫹1兲, 共2.13兲 with ␣ (0) being a constant that characterizes the relaxation rate of the spins belonging to level n⫽0, whose dynamics is not constrained. As ␣ (0) determines the basic time scale, which is arbitrary, we will take ␣ (0)⫽1 in the following. In the mean-field approximation just introduced, the time evolution of the average value of the spin, Eq. 共2.10兲, takes the form d dt 具 ␴ (n) 典 ⫽⫺ 兿 j⫽0 n⫺1 P共 ␴ 1 (j)⫽⫹1,..., ␴ ␮ j (j)⫽⫹1兲共 具 ␴ (n) 典 ⫺ 具 ␴ 典 e兲.共2.14兲 Now, a new approximation will be made. The probability P(j)in Eq. 共2.14兲will be substituted by its equilibrium value. This would be exact in linear response around equilibrium, but it will be taken here as an approximation leading to the basic equation of our hierarchically constrained model, namely, d dt 具 ␴ (n) 典 ⫽⫺ 兿 j⫽0 n⫺1 pe ␮ j共 具 ␴ (n) 典 ⫺ 具 ␴ 典 e兲,共2.15兲 where peis the equilibrium probability of any spin being in the up state, given by Eq. 共2.5兲. Again, on physical grounds, this is a sensible approximation due to the separation of the characteristic time scales of the different levels. Because of the hierarchically constrained dynamics, spins in level n⫺1 reach equilibrium over a time scale in which spins in level n have not begun to evolve. However, as a consequence of the last approximation, the configuration with all the spins in the down state is no longer absorbent. From Eq. 共2.15兲, an equivalent equation can be written for the evolution of the probability p(n)of the up state in level n, defined in Eq. 共2.3兲, i.e., d dtp(n)⫽⫺ ␣ n共p(n)⫺pe兲,共2.16兲 ␣ nbeing the characteristic relaxation rate of level n, ␣ n⫽pe gn,gn⫽兺 j⫽0 n⫺1 ␮ j.共2.17兲 Equation 共2.16兲implies that, due to the hierarchical constraints, the spin relaxation slows down with increasing level n, since ␣ nis a decreasing function of n, because pe⬍1. Equation 共2.16兲is the main result in this section. In the following, we will explore its implications, considering first processes at constant temperature in the next section, and cooling and heating processes in the remainder of the paper. III. RELAXATION AT CONSTANT TEMPERATURE For the case of constant temperature T, and therefore constant ␣ n, Eq. 共2.16兲is easily solved, p(n)共t兲⫺pe⫽关p(n)共0兲⫺pe兴e⫺ ␣ nt.共3.1兲 Then, each spin relaxes exponentially to equilibrium with the rate characteristic of its level. For the homogenous situations within each level we are considering, the dimensionless mean energy per spin 关30兴 defined in Eq. 共2.2兲simplifies to ␧共t兲⫽兺 n⫽0 nmax wnp(n)共t兲,共3.2兲 where wn⫽Nn/Nis the fraction of spins in level n, verifying 兺n⫽0 nmaxwn⫽1. Putting Eq. 共3.1兲into Eq. 共3.2兲yields ␧共t兲⫽␧e共T兲⫹兺 n⫽0 nmax wn关p(n)共0兲⫺pe兴e⫺ ␣ nt,共3.3兲 for the relaxation of the energy at constant temperature. In order to proceed, we will consider the simple case in which the initial probability distribution p(n)(0) does not depend on the index level n. This will be the situation, for instance, when the initial state corresponds to equilibrium at a different temperature T⫹⌬T. Thus, the relaxation function of the physical property described by the Hamiltonian of the system is given by ␾ 共t兲⬅␧共t兲⫺␧e ␧共0兲⫺␧e ⫽兺 n⫽0 nmax wne⫺ ␣ nt.共3.4兲 A. PRADOS AND J. J. BREY PHYSICAL REVIEW E 64 041505 041505-4 This equation is a generalization of the result derived by Palmer et al. in their pioneering work in hierarchically constrained dynamics 关7兴, which corresponds to the choice pe ⫽1/2. This is formally equivalent to the particularization of Eq. 共3.4兲for T→⬁. If other positive values of the temperature are considered, the effect is an increase of the relaxation times ␶ n⫽ ␣ n ⫺1,共3.5兲 since peis a decreasing function of the temperature. A main advantage of the formulation of the hierarchical models as presented here, aside from its larger generality, is that it allows analysis of processes in which the temperature of a thermal system 共or the vibration intensity in a granular system兲changes in time. This kind of processes will be addressed in the next section. A mean relaxation time ␶ can be defined as ␶ ⫽ 冕 0 ⬁dt ␾ 共t兲⫽兺 n⫽0 nmax wn ␶ n,共3.6兲 providing a quantitative measure of the time it takes the system to relax to equilibrium at temperature T. Let us consider that the fraction of spins in level n,wn, and the number of ‘‘facilitating’’ spins in level n, ␮ n, depend very smoothly on n; i.e., they can be expressed as functions of the form wn⫽w共n ␩ 兲 ␮ n⫽ ␮ 共n ␩ 兲,共3.7兲 where ␩ Ⰶ1. These seem to be sensible conditions when modeling a real system, in which the introduction of the levels and the pseudospins is associated to some coarsegrained description. By defining xn⫽n ␩ ,共3.8兲 which is a continuous variable in the limit ␩ →0, the sums over ncan be replaced by integrals. The relaxation rate of level n, given by Eq. 共2.17兲, becomes a function of the continuous variable x, ␣ 共x兲⫽pe g(x),共3.9兲 with g共x兲⫽ 兺 k⫽0 n⫺1 ␮ k ␩ 兺 k⫽0 nmax wk ␩ ⫽ 冕 0 xdx⬘ ␮ 共x⬘兲 冕 0 xmaxdx⬘w共x⬘兲 ,共3.10兲 where xmax⫽nmax ␩ , and the normalization of the weights wn has been used. Therefore, the relaxation function ␾ (t), given by Eq. 共3.4兲, becomes ␾ 共t兲⫽ 兺 n⫽0 nmax ␩ wne⫺ ␣ nt 兺 n⫽0 nmax ␩ wn ⫽ 冕 0 xmaxdxw共x兲e⫺ ␣ (x)t 冕 0 xmaxdxw共x兲 , 共3.11兲 in the continuous limit. In general, Eq. 共3.11兲is mathematically rather involved, since it depends both on the functions ␮ (x) and w(x). The simplest possibility appears to be ␮ (x) proportional to w(x), i.e., ␮ nproportional to wnor, equivalently, to Nn. In other words, the simplest kind of hierarchically constrained models shows up when the number of ‘‘facilitating’’ spins at a level is an extensive function of the number of spins at the same level 关11兴. This condition is expressed as ␮ 共x兲⫽aw共x兲,共3.12兲 with abeing a constant, independent of x. In this case, it is useful to define the new variable u⫽ 冕 0 xdx⬘w共x⬘兲 冕 0 xmaxdx⬘w共x⬘兲 ,共3.13兲 measuring the fraction of spins belonging to levels up to n ⫽x/ ␩ . In terms of u, the relaxation rates of Eq. 共3.9兲are given by ␣ 共u兲⫽pe g(u),g共u兲⫽au,共3.14兲 the relaxation function is expressed as ␾ 共t兲⫽ 冕 0 1due⫺ ␣ (u)t,共3.15兲 and the mean relaxation time reads ␶ ⫽ 冕 0 ⬁dt ␾ 共t兲⫽ 冕 0 1du ␶ 共u兲⫽pe ⫺a⫺1 a 兩 ln pe 兩 ,共3.16兲 with ␶ (u)⫽ ␣ ⫺1(u). It is interesting to consider situations for which pe aⰆ1, so that the minimum relaxation rate ␣ (1) is much smaller than the maximum one ␣ (0)⫽1 in our dimensionless time scale. In this case, the relaxation function ␾ (t) is linear in lntover an intermediate time window, 1 ⰆtⰆpe ⫺a, namely 关11兴, ␾ 共t兲⬃1⫺1 a 兩 ln pe 兩 共 ␥ ⫹lnt兲,共3.17兲 where ␥ stands for Euler’s constant, ␥ ⯝0.577. This kind of linear logarithmic behavior is characteristic of a great variety of complex systems, including spin glasses 关12,13兴, granular materials 关9,10,14,15兴, structural glasses 关16–18兴, and protein models 关8,19兴. In the present context, the condition pe a Ⰶ1 corresponds to a ‘‘low’’-temperature limit, in which the GLASSLIKE DYNAMICAL BEHAVIOR IN . . . PHYSICAL REVIEW E 64 041505 041505-5 mean relaxation time, given by Eq. 共3.16兲, is very large; i.e., the relaxation of the system becomes very slow. It is worth noting that, in terms of the uvariable, the continuous limit is formally obtained by changing the functions of the index level nby the corresponding functions of the variable uand by making the replacement 兺 n⫽0 nmax wn→ 冕 0 1du.共3.18兲 It follows that, in the continuous limit, the dynamical behavior of the system does not depend explicitly on the level populations Nn, but only on the relaxation rates ␣ expressed as functions of u, as given by Eq. 共3.14兲. IV. TIME-DEPENDENT TEMPERATURE PROCESSES In this section processes in which the temperature changes in time will be studied. The evolution equation 共2.16兲is now d dtp(n)共t兲⫽⫺ ␣ n共T兲关p(n)共t兲⫺pe共T兲兴,共4.1兲 where T⫽T(t). The general solution of Eq. 共4.1兲is p(n)共t兲⫽关p(n)共t0兲⫺pe共T0兲兴 ␹ n共t,t0兲⫹pe共T兲 ⫺ 冕 t0 tdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ␹ n共t,t⬘兲.共4.2兲 Here T0⫽T(t0) is the initial value of the temperature, T ⫽T(t), T⬘⫽T(t⬘), T⬙⫽T(t⬙), and we have introduced the function ␹ n共t1,t2兲⫽exp 冉 ⫺ 冕 t2 t1dt ␣ n共t兲 冊 .共4.3兲 The above equation is valid for any law of variation for the temperature. Taking into account Eq. 共3.2兲, the average energy per spin, ␧, is given by 关30兴 ␧共t兲⫽兺 n⫽0 nmax wn关p(n)共t0兲⫺pe共T0兲兴 ␹ n共t,t0兲⫹␧e共T兲 ⫺ 冕 t0 tdt⬘d␧e共T⬘兲 dT⬘ dT⬘ dt⬘兺 n⫽0 nmax wn ␹ n共t,t⬘兲,共4.4兲 where ␧eis the average equilibrium energy, defined in Eq. 共2.6兲. Let us assume that the system is initially at equilibrium with T⫽T0. Then the first term on the right-hand side 共RHS兲 of Eq. 共4.4兲vanishes and ␧共t兲⫽␧e共T兲⫺ 冕 t0 tdt⬘d␧e共T⬘兲 dT⬘ dT⬘ dt⬘M共t,t⬘兲,共4.5兲 where M共t,t⬘兲⫽兺 n⫽0 nmax wn ␹ n共t,t⬘兲共4.6兲 is a memory function. In the case of constant temperature, M(t,t⬘) is equal to the relaxation function ␾ (t⫺t⬘), as seen by comparing Eq. 共4.6兲with Eq. 共3.4兲. The structure of Eq. 共4.5兲is the same as that of Narayanaswami’s phenomenological theory of glasses 关22–25兴. A similar result was obtained some years ago for the one-dimensional Ising model with Glauber dynamics 关31兴. High-temperature limit: Hilbert’s method We are going to look for a solution of Eq. 共4.1兲by means of Hilbert’s method. A special solution pH (n)共t兲⫽兺 k⫽0 ⬁ pH (n),k共t兲共4.7兲 is constructed in an iterative way as follows. We take pH (n),0共t兲⫽pe共T兲,共4.8兲 while for k⭓1 pH (n),k共t兲⫽⫺ ␣ n ⫺1共T兲dp(n),k⫺1共t兲 dt .共4.9兲 Equation 共4.8兲shows that Hilbert’s expansion agrees with the equilibrium distribution to the lowest order. Besides, for k⫽1 we get from Eq. 共4.9兲 pH (n),1共t兲⫽⫺ ␶ n共T兲dpe共T兲 dT dT dt .共4.10兲 This equation indicates the main limitation of Hilbert’s method. Due to the divergence of the relaxation times ␶ n ⫽ ␣ n ⫺1in the low-temperature limit 关see Eq. 共2.17兲兴, also p(n),1 diverges in that limit. As a consequence, Hilbert’s solution is only accurate in the high-temperature regime, in which an expansion around equilibrium provides a good approximation. Restricting ourselves to high temperatures, we approximate pH (n)共T兲⯝pe共T兲⫺ ␶ n共T兲dpe共T兲 dT dT dt 共4.11兲 and, from Eq. 共3.2兲, ␧H共T兲⯝␧e共T兲⫺d␧e共T兲 dT dT dt 兺 n⫽0 nmax wn ␶ n共T兲.共4.12兲 Taking into account the definition of the average relaxation time ␶ , Eq. 共3.6兲, the above expression is seen to be equivalent to ␧H共T兲⯝␧e共T兲⫺d␧e共T兲 dT dT dt ␶ 共T兲,共4.13兲 A. PRADOS AND J. J. BREY PHYSICAL REVIEW E 64 041505 041505-6 which agrees with the high-temperature behavior of Eq. 共4.5兲. In Eqs. 共4.11兲and 共4.13兲,pH (n)and ␧Hdepend on time only through the temperature T(t). Thus, Hilbert’s method provides a ‘‘normal’’ solution, in the sense often used in kinetic theory. The validity of a expression identical to Eq. 共4.13兲, also in the high-temperature limit, has been established for a quite general class of systems whose dynamics is described by a master equation 关28兴. Although we have made here several drastic approximations in order to get a closed equation for the average spin, the high-temperature limit of the solution, given by Hilbert’s method, remains formally the same as that of the exact solution of the original model. Certainly, this is a good property of those approximations. What is the physical meaning of the failure of Hilbert’s expansion for low temperatures? Due to the divergence of the characteristic relaxation times, the system does not have enough time to relax to the equilibrium curve at very low temperatures, and it gets ‘‘frozen’’ in a far-from-equilibrium state. Since Hilbert’s method is an expansion around equilibrium, it fails in the low-temperature region. In fact, the RHS of Eq. 共4.11兲becomes negative for low enough temperatures. On the other hand, Hilbert’s expansion is useful to estimate the values of the physical properties in the ‘‘frozen’’ state 关28兴. Also, Hilbert’s method provides a qualitative understanding of the hysteresis effects appearing in thermal cycles 共cooling and reheating兲. In cooling processes (dT/dt⬍0), it is ␧H⭓␧e, while in heating processes (dT/dt⬎0), it is ␧H ⭓␧e. Then, ␧Hlies to opposite sides of the equilibrium curve for cooling and heating processes, and hysteresis effects show up in thermal cycling experiments, as will be discussed in more detail in Sec. VI. V. COOLING PROCESSES Next, we are going to study the continuous cooling of the system down to very low temperatures. The origin of time is taken at the beginning of the cooling process. The initial condition will be the equilibrium configuration at a ‘‘high’’ temperature T0, i.e., p(n)共0兲⫽pe共T0兲.共5.1兲 Then, the first correction in Hilbert’s expansion, p(n),1(t), is very small as compared with pe(T0) for T→T0. Particularization of Eq. 共4.2兲for the above initial condition gives p(n)共t兲⫽pe共T兲⫺ 冕 0 tdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ␹ n共t,t⬘兲.共5.2兲 Since pe(T) is an increasing function of the temperature, for continuous cooling processes it is p(n)共t兲⭓pe共T兲for all tand n.共5.3兲 The possible deviations from the equilibrium distribution always lead to an increase of the probability of the spin being in the excited state. Moreover, Eq. 共5.2兲directly implies that ␧共t兲⫽␧e共T兲⫺ 冕 0 tdt⬘d␧e共T⬘兲 dT⬘ dT⬘ dt⬘兺 n⫽0 nmax wn ␹ n共t,t⬘兲⭓␧e共T兲, 共5.4兲 where we have used Eq. 共3.2兲. This inequality has been experimentally observed in glass-forming liquids 关22,23兴. Since the reported experiments were made at constant pressure, the quantity ␧considered here must be interpreted as the enthalpy in that context. In order to proceed further in our analysis, the continuous limit introduced in the study of the relaxation at constant temperature in Sec. III will be considered. As already mentioned, this continuous limit is expected to be closer to the description of real systems than the discrete level picture. Besides, for the sake of concreteness, we will restrict ourselves to those models verifying Eq. 共3.12兲. The index level nis substituted by the continuous variable u, defined in Eq. 共3.13兲, representing the fraction of the total number of spins up to level n. With an obvious change of notation, Eq. 共5.2兲 becomes p共t;u兲⫽pe共T兲⫺ 冕 0 tdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ␹ 共t,t⬘;u兲,共5.5兲 where ␹ 共t,t⬘;u兲⫽exp 冉 ⫺ 冕 t⬘ tdt⬙ ␣ 共T⬙;u兲 冊 ,共5.6兲 with ␣ (T;u) given by Eq. 共3.14兲, i.e., ␣ 共T;u兲⫽pe共T兲au.共5.7兲 Also, using Eq. 共3.18兲, it is found that ␧共t兲⫽␧e共T兲⫺ 冕 0 tdt⬘d␧e共T⬘兲 dT⬘ dT⬘ dt⬘ 冕 0 1du ␹ 共t,t⬘;u兲. 共5.8兲 The time evolution of the probability p(t;u) depends on the explicit form of the cooling law. In experiments, linear cooling is usually employed, dT dt ⫽⫺rc,共5.9兲 where rc⬎0 is the cooling rate determining the time scale rc ⫺1over which the temperature changes. Linear cooling implies that pe(T) depends on time in a rather involved way. From Eqs. 共2.5兲and 共5.9兲one gets dpe dt ⫽⫺ 1 2rcpe共1⫺pe兲 冉 ln pe 1⫺pe 冊 2 .共5.10兲 We are interested in cooling processes for which the temperature changes slowly in time, rcⰆ1, so that the system departs from the equilibrium curve for very low temperatures, where peⰆ1. Then, a law equivalent to linear cooling, aside from logarithmic corrections, is GLASSLIKE DYNAMICAL BEHAVIOR IN . . . PHYSICAL REVIEW E 64 041505 041505-7 dpe dt ⫽⫺rcpe,共5.11兲 having the advantage that analytical calculations are much more simple with Eq. 共5.12兲关31–33兴. In the following, we will use the notation pe共T兲⫽pe,pe共T⬘兲⫽pe ⬘,pe共T0兲⫽pe0.共5.12兲 The time integrals in Eq. 共5.5兲can be transformed into integrals over peby means of the cooling law 共5.11兲, with the result p共t;u兲⫽pe⫹ 冕 pe pe0dpe ⬘exp 冉 ⫺pe ⬘au⫺pe au rcau 冊 .共5.13兲 The second term in this expression is dominant in the lowtemperature region, where peis very small. Therefore, it follows that spins in any level ufall in a nonequilibrium state for low enough temperatures. The details of this glasslike transition will be analyzed below, in a separate subsection. One of the main quantities characterizing a given cooling process is the residual value fres of a relevant property f. The residual value measures the excess with respect to the equilibrium curve, extrapolated to very low temperatures, fres⫽lim T→0 共f⫺fe兲.共5.14兲 In particular, for the probability p(t;u) we have from Eq. 共5.13兲 pres共u兲⫽ 冕 0 pe0dpe ⬘exp 冉 ⫺pe ⬘au rcau 冊 ,共5.15兲 which is easily transformed into pres共u兲⫽1 au共rcau兲1 au 冕 0 x0dxx 1 au ⫺1e⫺x,共5.16兲 where x0⫽pe0 au/(rcau). The slow cooling limit is defined by the residual properties being independent of the initial conditions and determined univocally by the cooling rate 关31– 34兴. In our case, slow cooling means that the upper integration limit in Eq. 共5.16兲can be substituted by infinity for all u, i.e., pe0 au rcauⰇ1. 共5.17兲 As the uvariable varies in the interval 0⭐u⭐1, the slow cooling condition is rcaⰆpe0 a⬍1. 共5.18兲 Then, with an exponentially small error, pres共u兲⬃1 au共rcau兲1/au⌫ 冉 1 au 冊 ⫽共rcau兲1/au⌫ 冉 1⫹1 au 冊 . 共5.19兲 This expression gives the probability that the spins in level u be in the up state at very low temperatures. In Fig. 2, the residual probability pres(u⫽1) is plotted as a function of the cooling rate rc, for a⫽1. The asymptotic result, given by Eq. 共5.19兲, is compared with the numerical integration of Eq. 共5.16兲with pe0⫽1/2; i.e., the system is taken initially at infinite temperature. The agreement is quite good up to rc ⯝0.1, which is not very small. The evolution of the average energy, for the cooling law 共5.11兲, is given by ␧共t兲⫺␧e共T兲⫽ 冕 0 1du关p共t;u兲⫺pe兴 ⫽ 冕 0 1du 冕 pe pe0dpe ⬘exp 冉 ⫺pe ⬘au⫺pe au rcau 冊 . 共5.20兲 The residual energy can be easily computed by particularizing this expression for T→0, ␧res⫽ 冕 0 1du pres共u兲.共5.21兲 The integrand pres(u), given by Eq. 共5.19兲, vanishes exponentially in the limit u→0. A standard Laplace analysis can be made, with the result ␧res⬃⌫ 冉 1⫹1 a 冊 共rca兲1/a 兩 ln共rca兲1/a 兩 .共5.22兲 The leading behavior is potential with rc, since ln␧res⬃1 aln共rca兲,共5.23兲 FIG. 2. Residual probability for the slowest modes pres(u⫽1) as a function of the dimensionless cooling rate rcdefined in the main text. The circles correspond to the numerical integration of Eq. 共5.16兲, while the solid line is the prediction of the asymptotic calculation, Eq. 共5.19兲. A good agreement is observed up to rc⯝0.1. A. PRADOS AND J. J. BREY PHYSICAL REVIEW E 64 041505 041505-8 which comes from the upper limit of integration, u⫽1, corresponding to the largest relaxation time. The integral over the whole distribution function pres(u) gives a logarithmic correction 兩 ln(rca) 兩 ⫺1, which makes the residual energy smaller than the dominant term (rca)1/a. This is due to the increasing behavior of pres(u) with u,pres(u)⭐pres(u⫽1). In Fig. 3 the residual value of the energy is plotted. The asymptotic expression, Eq. 共5.22兲, is compared with the numerical results from Eqs. 共5.21兲and 共5.16兲. Good agreement is found up to rc⫽0.01. It is worth noting that, for the cooling law considered, the logarithmic correction to the potential in rcbehavior is not present for other simple models of structural glasses previously studied 关31–35兴. Demarcation mode, fictive temperature, and glass transition The existence of nonvanishing residual properties is an indication of the departure of the system from equilibrium at low temperatures. Due to the divergence of the characteristic relaxation times ␶ n, for low enough temperatures the system does not have enough time to relax towards equilibrium, and a kinetic phenomenon resembling the laboratory glass transition 关22–24兴shows up. Next, we will try to understand the physical origin of this kinetic transition. Let us consider again the time evolution of the probability distribution p(t;u) in a cooling process, as given by Eqs. 共5.5兲and 共5.6兲. The integral in ␹ (t,t⬘;u), I共t,t⬘;u兲⫽ 冕 t⬘ tdt⬙ ␣ 共T⬙;u兲,共5.24兲 is a measure of the average number of transitions occurring in level uin the time interval between t⬘and t. Consequently, a mathematical definition for the limit of slow cooling is that the condition I共t*,0;u兲Ⰷ1共5.25兲 hold for all u, where t*is the time for which the temperature vanishes if extrapolated accordingly to the prescribed cooling law, i.e., T(t*)⫽0. Equation 共5.25兲guarantees that the system experiments a large number of transitions before getting eventually frozen, so that it has enough time to forget the details of the initial condition. For the cooling law defined in Eq. 共5.11兲it is easily verified that Eq. 共5.25兲is equivalent to Eq. 共5.17兲. If we are dealing with a slow cooling process, Eq. 共5.25兲 implies that there is a time window over which I共t,0;u兲Ⰷ1. 共5.26兲 This is the time regime we are interested in. Let us analyze the behavior of ␹ (t,t⬘;u) as a function of t⬘,t0⭐t⬘⭐t, for a given time tsuch that Eq. 共5.26兲holds. The function I(t,t⬘;u) changes from a very large value to zero when t⬘ goes from 0 to t. Consequently, ␹ (t,t⬘;u) increases from practically zero to unity when t⬘moves in the above time interval. Let us define a time tf(t;u), prior to t,by I共t,tf;u兲⫽1, 共5.27兲 so that the average number of transitions taking place in level uin the time interval between tfand tequals unity. Then, ␹ (t,tf;u)⫽e⫺1, and the function ␹ (t,t⬘;u) changes from zero to unity in a certain time interval around tf.In order to proceed, we will assume that this change takes place in the vicinity of tfvery rapidly, as compared with the variation of the rest of the integrand of Eq. 共5.5兲. Decomposing Eq. 共5.5兲in the form p共t;u兲⫽pe共T兲⫺ 冕 0 tfdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ␹ 共t,t⬘;u兲 ⫺ 冕 tf tdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ␹ 共t,t⬘;u兲,共5.28兲 the first integral is subdominant with respect to the second one and, moreover, 冕 tf tdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ␹ 共t,t⬘;u兲⯝ 冕 tf tdt⬘dpe共T⬘兲 dT⬘ dT⬘ dt⬘ ⫽pe共T兲⫺pe关Tf共t;u兲兴, 共5.29兲 where Tf(t;u) is the temperature of the system at time tf(t;u), i.e., Tf共t;u兲⫽T关tf共t;u兲兴.共5.30兲 Note that Tf(t;u)⬎T, since the time instant tf(t;u)⬍t. Substitution of Eq. 共5.29兲into Eq. 共5.28兲and use of the above approximations yields p共t;u兲⫽pe关Tf共t;u兲兴.共5.31兲 The arguments leading from Eq. 共5.5兲to Eq. 共5.31兲are formally equivalent to assume that FIG. 3. Dimensionless residual energy ␧res as a function of the cooling rate rc. As in Fig. 2, the circles are from the numerical integration of Eq. 共5.16兲and the solid line is the prediction of the asymptotic analysis, given by Eq. 共5.22兲. The agreement is good for rcⱗ0.01. GLASSLIKE DYNAMICAL BEHAVIOR IN . . . PHYSICAL REVIEW E 64 041505 041505-9