Continuous cooling of a one-dimensional bonded fluid: A Monte Carlo simulation study

PHYSICAL REVIE%' 8VOLUME 43, NUMBER 11JANUARY 1991 Continuous cooling of aone-dimensional bonded fluid: AMonte Carlo simulation study J.J.Brey and M. J.Ruiz-Montero Departamento de I'Isica Teorica, Uniuersidad de Seuilla, Apartado de Correos 1065, Sector Sur, E-41080Seuilla, Spain (Received 23 July 1990) Akinetic one-dimensional bonded Quid model is formulated by extending amodel for water due to Bell. The behavior of the system under continuous cooling is analyzed by use of aMonte Carlo simulation method. It shows aphenomenon similar to the laboratory glass transition. The dependences of the glass transition temperatures and of the residual energy on the cooling rate are discussed. I. INTRODUCTION In spite of the intense scientific effort carried out in recent years, the nature of the liquid-glass transition is far from being we11 understood. Although it is clear that the glass transition as observed in real experiments is akinetic phenomenon associated with the departure of the liquid from the equilibrium, little is known about the physical processes controlling the behavior of the system through the transition. 'In some theories and models the laboratory glass transition is the manifestation of an underlying thermodynamic or kinetic phase transition, where one or more relaxation times diverge. Nevertheless, it has also been shown that an underlying phase transition is not a necessary condition for amodel to present glass transition phenomena similar to those observed in real glasses. Atypical way of preparing aglass is by the fast cooling of aliquid. In fact, it seems that any liquid can be led to aglassy state if the cooling rate is large enough. Besides, the properties of the resulting glass depend on the cooling procedure. Therefore, the study of the time evolution of systems whose temperature is externally changed shows up as acritical step for the understanding of the glass transition. This poses us with aformidable task from a theoretical point of view. It is not even clear how to formulate, in an appropriate way, akinetic equation incorporating the mechanism controlling the temperature. Of course, this can be done, in principle, by considering time-dependent boundary conditions in the kinetic equation for an isolated system. Nevertheless, this introduces inhomogeneities into the system and renders the problem extremely di%cult. Asimpler, although unjustified, possibility is to assume that the equation describing the evolution of asystem at constant temperature can be directly extended by formally replacing the constant value of the temperature by the given time function. Most of the models invented to analyze the dynamics of glasses are formulated in this way. Nevertheless, even with this choice, the explicit analytical solution for the case of continuous cooling is known only for afew very simple models in some asymptotic limits. ' Recently, several models have been proposed that mimic the dynamics of real fiuids near the glass transition. Although no fully satisfactory theory has been developed to explain their kinetic properties yet, their simplicity allows very efFicient computer simulations. As Fredickson has pointed out, ageneral characteristic of these models is that they contain some dynamical constraints that are spatially short ranged, but lead to ahigh cooperativity in the dynamics at larger distances. The purpose of this paper is to investigate, using Monte Carlo simulation methods, the time evolution under continuous cooling of amodel of one-dimensional bonded fiuid. This model was introduced by Bell as a model for liquid water. It is alattice model for which only the statics was specified. The free Gibbs energy can be easily computed and, from it, analytical expressions for the thermodynamic properties are derived. Here we introduce adynamics of the bonded fiuid by amaster equation with transition probabilities obeying detailed balance. It must be noticed that Bell's model has been previously used to study the properties of glasses by Singh and Kovac. Nevertheless, their approach is quite different from the one we present here. They assume that, in the quenched fiuid, there is some kind of disorder associated with agiven distribution ofbond energies. No attempt to describe the dynamics of the glass transition is made and, therefore, the kinetics of the model was not defined. Section II defines our model. Explicit expressions for some of the thermodynamical properties are given and the kinetics of the system is established. The Monte Carlo simulation method employed is outlined in Sec. III. The results obtained for the evolution of the system are presented in Sec. IV. There, it is seen that the system exhibits abehavior similar to alaboratory glass transition. The dependence of the residual energy and the glass transition temperature on the cooling rate is analyzed. Section Vcontains some comments and conclusions. II. DESCRIPTION OF THK MODEL We consider aone-dimensional lattice of Nequally spaced sites. Each site can be either occupied by aparticle or it can be empty. The number of particles is M. Two nearest-neighbor particles have an interaction ener43 1991 The American Physical Society 586 J.J.BREY AND M. J.RUIZ-MONTERO 43 E(N, M, Z,B)=— (M — Z)E— Bc@ .(2.1) The number of possible configurations with given values of Zand 8is g(N, M, Z,B) (M — 1)!N(N — M— B— 1)! Z! (Z — 1)!(M — Z— B)!(N — M— Z— B)! gy — c. Also, two particles separated by asingle empty lattice site may establish between them abond of energy — (8+co). Both Eand co are taken as positive. These are the only interactions in the system. Following Bell, we will refer to asite between two bonded particles as occupied by abond. The empty sites are considered as holes. Let us think of agiven configuration of the system. We will assume periodic boundary conditions. The energy of the system can be expressed in terms of the number of bonds 8and the number of particle-hole contacts, denoted by 2Z, as Refs. 3and 4, Eq. (2.6) does. In fact, the latter can be directly derived by using the transfer matrix method. By the differentiation of the Gibbs free energy, all the equilibrium thermodynamical properties of the system can be derived. In particular, the equation of state reads aG N=— ADJ '+J(J— 1) ADJ '+ A+(J— 1) (2.7) and the average energy is BG —AE+ADJ '(E+co) ADJ-'+ A+(J— 1)-' (2.8) Other properties can be derived by using the isothermisobaric distribution. The average number of bonds is given by (2.2) 8ink BD ADJ ADJ '+ A+(J— 1) (2.9) Since this expression differs from those reported by Bell and Singh and Kovac, aderivation of it is presented in the Appendix. Equation (2.2) is useful in estimating finite-size effects in acomputer simulation of the model. In our case it has been checked that these effects are negligible for the size of the system we have considered. Now the partition function is constructed from Eqs. (2.1) and (2.2). It is convenient to work in the isothermisobaric ensemble where the partition function for our model is given by M— Zco 6(M,p, T)=gggg(N, M, Z,B) Z=O B=0 N=M+B+Z Xexp[ 13(E+pNl) j, — (2.3) where Tis the temperature, P=(k~T) ', pis the pressure, and Iis the distance between neighboring sites. Substitution of Eqs. (2.1) and (2.2) into Eq. (2.3) yields, after carrying out the summations, 6(M,p, T)=(J— 1) 'J +'[(J— 1) +DAJ '+ Aj Since we are dealing with aone-dimensional model with short-range interactions, the system does not present any equilibrium phase transition. In the following we will consider asystem with aconstant length, i.e.,number of sites X. The equilibrium pressure at each temperature is determined by Eq. (2.7). Once the thermodynamics of our system is known, we specify its dynamics. It is assumed to be governed by a master equation. Taking into account that the number of sites and particles are kept fixed, the allowed transitions consist either in the creation (or destruction) of abond, or in the motion of aparticle to anext empty site. Of course, abond can only be created starting from a configuration in which two particles are separated by a single hole. Also, aparticle cannot jump to asite that is being "occupied" by abond. Before that, the bond must be destroyed. This implies that at low temperatures„ when the concentration of bonds is large, the relaxation of the model requires the system to jump over high energy barriers. If xand x' denote two configurations of the system connected by an allowed single transition, the transition rate used is X[J '+2AD(J — 1)J AE I— tanh 2kB T(2.10) +AJ '(J — 1)+(J— 1) 'j, with A=e~', D=e~, and J=e~~. The Gibbs free energy 6is 6=— kBT ink and, in the limit ofMgoing to infinity, one gets 6=MkBT ln J ADJ '+ A+(J— 1) (2.4) (2.5) (2.6) Although Eq. (2.4) does not agree with those derived in where QE =Q(x') E(x) and zo is a — constant that fixes the natural time unit of our system. The dynamics defined in this way satisfies the detailed balance condition. Besides, it is clear that all the configurations can be reached from agiven one by astring of transitions with nonzero probability, i.e.,the Markov process defined by our equation is irreducible. It follows that any initial distribution of configurations will tend towards the equilibrium canonical distribution. It must be noticed that we have not introduced in the dynamics any constraint other than the purely energetic 43 CONTINUOUS COOLING OF AONE-DIMENSIONAL BONDED. .. ones associated with the destruction of bonds. In this sense our model differs from the facilitated Ising models 'and also from the tiling models. ' For convenience we have introduced an adimensional cooling rate as III. MONTE CARLO SIMULATION r~0 5= (E+co) (4.2) The simulation of the system has been carried out by means of aMonte Carlo algorithm similar to that formulated by Bortz, Kalos, and Lebowitz. 'For computational reasons, this method was preferred to the standard Metropolis Monte Carlo method, where the probability of rejection of amovement increases very fast as the temperature is lowered. Given the configuration of the system at atime t, the transition rates for all the possible transitions are computed according to Eq. (2.10). The total transition rate from the given state O(x) =gW(x ~x') (3.1) — lna Q(x) (3.2) where ais another random number between 0and 1. As indicated before, the time was measured in Monte Carlo steps, taking so= 1in Eq (2.10. ). For the interaction energies we chose c.=200 and co =500. The data presented in this paper have been obtained with alattice of 200 sites and 100 particles, and periodic boundary conditions have been used. For the highest cooling rates considered, we have checked that the results are not significantly altered when the number of particles is increased to 1000 (2000 sites). Besides, by using the transfer matrix method we have calculated the normalized equilibrium density correlation function. Of course, it is an oscillating function whose amplitude monotonously decays with distance. For kz T=400 it is smaller than 10 for distances greater than 10 lattice spacings. For kz T=200, that is the smallest temperature at which the system remains near equilibrium in our simulation, the correlation has already decayed to 10 for distances of 30 lattice spacings. For each cooling rate, 100 independent cooling runs were performed and the results were averaged. In spite of the rather small size of our system, the runs at the lowest cooling rates required more than 35 hin aCONVEX-220 computer. IV. CONTINUOUS COOLING We have cooled the system following alinear law for the variation of the temperature, i.e.,with the constant cooling rate d(k~ T) dt (4.1) is also calculated. Arandom number uniformly distributed between 0and 1is then generated to determine which of the possible transitions takes place. This is done by dividing the interval (0,1) into anumber of intervals equal to the number of possible transitions. To each of them it corresponds aportion of length W'(x~x )/Q(x). The time interval associated with this change in the configuration of the system is C3 C3 C3 r—0.1 (D ~= a.o 00 e,.mo0 "i 00 0 0 0 0 200 300 400 500 T FIG. 1. Average energy as afunction of the temperature in continuous cooling experiments for r=1, 0.1, and 0.01. Also, the equilibrium values have been plotted (circles). We started with the equilibrium configuration at P= — ca that, for the system considered here, is an open structure with ahole between every two particles. As is the case for the facilitated Ising model, the fact that the temperature is negative does not present any physical problem since the model does not show pathological behavior in passing from negative to positive temperatures. Besides, the state at 13= — ~is only used in our simulation as a convenient initial state, from which we carry out an instantaneous quench to alow temperature. This temperature was chosen in such away that significant departures from the equilibrium behavior when cooling the system at the considered rate show up well below it. Its value ranged between k~T=500 for r=1 and k~T=300 for r=0.005. The aim of the quench was to save computer time. Then, the system was equilibrated at aconstant temperature until the relative difference between the averaged energy over 100 runs and the equilibrium value given by Eq. (2.8) was of the order of 1%. This required runs of typically 10 Monte Carlo time steps. After that, the system was continuously cooled. Ten difterent values of the constant rate in the range 2.5X10 ~r~1have been investigated. That interval corresponds to 3.6X10 ~6 &1.4X10 '. Figure 1shows the time evolution of the average energy as afunction of the temperature in some of the continuous cooling experiments. For the sake of clarity only the results corresponding to three values of the cooling rate are represented, namely r=1,0.1, and 0.01. Nevertheless, the following comments apply to all the cases studied. It is seen that the system has abehavior similar to the laboratory glass transition observed in real glasses. For agiven value of the cooling rate, the evolution of the average energy follows the equilibrium curve at high temperatures, but the system falls out of equilibrium when J.J.BREY AND M. J.RUIZ-MONTERO 43 e„,(5)=5 (4.3) with a=0.199. For 0.05 &r&1the fit gives o,"=0.139. To see whether alogarithm dependence on the cooling rate is more accurate, in Fig. 3we have plotted the logarithm of the residual energy against ln( — ln5). Again, two linear regions show up quite clearly. Each of them can be fitted to abehavior of the form e„,(5)=(— ln5) (4.4) with y=2.184 for r&0.05 and y' =1.087 for 0.05 &r&1. Although the fitting in Fig. 3is slightly better than in Fig. 2, our results are not at all conclusive about which of the laws given by Eqs. (4.1) and (4.2) describes best the 0 0 FIG. 2. Double logarithm representation of the residual energy as afunction ofthe cooling rate. sufBciently low temperatures are reached, and the energy departs from its equilibrium value. If we define the glass transition temperature Tas that corresponding to the point where the system falls out of equilibrium, then Tdecreases as the cooling rate rdecreases. Amore precise definition of Twill be given later on. Subsequent cooling leads the energy to avalue practically constant; that is, not affected by afurther decrease of the temperature. The difference between the value of the energy and the energy of the ground state is the so-called residual energy. In our model the energy per particle in the ground state is — (E+co). Figure 1also shows that the residual energy becomes smaller as the cooling rate is decreased. In an attempt to identify the dependence of the residual energy on the cooling rate, two different representations have been used. They correspond to the two laws usually considered in the bibliography. In Fig. 2adouble logarithm representation is shown. Although one could try to fit all the data to asingle straight line, it seems that two different regions can be identified, with acrossover around r=0.05 (5=7.1X10 ). Asimilar behavior of the residual energy was found by Kob and Schilling" for adifferent one-dimensional configurational glass model. For r&0.05, the results fit quite well to the power law 2.22.4in(-i&b) FIG. 3. Logarithm of the residual energy vs ln (— ln6). p(T)=p, (T~)+m„(T— TI), (4.5) where mis the variation ofpwith the temperature associated with the glass relaxation. The fictive temperature depends, in general, on the property under consideration. In acontinuous cooling experiment T& reaches astationary value that can be identified as the laboratory glass transition T.Its value is easily obtained by extending the glass behavior until reaching the equilibrium curve. The temperature of the intersection point is the glass transition temperature. In our case, the system becomes frozen at low temperatures. Therefore, Tis given by the equilibrium temperature of the system when it has a value of the macroscopic property equal to its residual value in the glass. From the above discussion it is clear that different properties can lead to different values of Tfor agiven value of the cooling rate. We have evaluated two glass transition temperatures for each of our cooling experiments by using the residual values of the energy and the number of bonds, respectively. The results are plotted in Fig. 4. Is is seen that there is aclear correlation between the values of both transition temperatures. Besides, their difference seems to go to zero as the cooling rate decreases. It must be noticed that, according to Eq. (2.1), the energy is not determined by the number of bonds and vice versa. Figure 4shows that, from apractical point of view, Tg is quite accurately given by alinear function of the logarithm of the cooling rate. Asimilar behavior has been found from phenomenological theories and also observed in some real glasses. dependence of the residual energy on the cooling rate. Aconcept that has proved to be very useful in phenomenological theories of glassy relaxation is that of fictive temperature. The fictive temperature T& of a glass is the temperature at which the equilibrium liquid has approximately the same structure as the glass. Usually, it is defined with reference to agiven property pof the system by 43 CONTINUOUS COOLING OF AONE-DIMENSIONAL BONDED. ..589 C3 O0 0 have tried to identify in our model the energy of apossible configurational defect and the corresponding barrier height leading to the values of aand a' reported below Eq. (4.3). Nevertheless, we have not found agreement, even approximate, for any of the possible simplest defects. The relation between the low-temperature relaxation properties of our model and an ensemble of two-level systems is an important point that deserves further study. ACKNOWLEDGMENTS — 12 We acknowledge partial support from the Direccion General de Investigacion Cientifica yTecnica (Spain) through Grant No. PB86-0205. FIG. 4. Glass transition temperature obtained from the energy (circles) and from the number ofbonds (asterisks). V. DISCUSSION The one-dimensional bonded fluid model discussed in this paper exhibits aphenomenon analogous to the laboratory glass transition observed in real glasses. The model is ergodic and its thermodynamic properties are known. It does not present any kind of thermodynamic or kinetic ideal glass transition. In this sense, it is similar to the facilitated Ising models and differs from the tiling models. ' Here we have restricted ourselves to continuous cooling computer experiments. The linear and nonlinear relaxation properties at constant temperature will be discussed elsewhere, both analytically and by means of a Monte Carlo simulation. Nevertheless, it can be anticipated that, at low temperatures, the relaxation of the model is dominated by adefect diffusion picture. At such temperatures, the concentration of bonds is very high, and the concentration of holes and particle-particle contacts is very low. In order to decrease the energy, the holes must move until they find aparticle-particle contact, and then create anew bond. The process requires the destruction of all the bonds separating the initial positions of the hole and the particle-particle contact. Of course, the bonds can be created again after the hole (or the particle) has jumped. We arrive then to apicture of defects diffusing in aperiodic potential. The same kind of picture was found for the one-spin facilitated Ising model. The residual energy of the model can be fitted quite well to both power laws and logarithm laws in certain ranges of the cooling rates. No clear conclusion has been established about which of the two laws leads to abetter fit of the data. Both of them have been justified theoretically by considering an ensemble of two-level systems. "' The logarithm law corresponds to asmooth distribution of low-excitation energies, while the power law is obtained when the excitation energy is the same for all the two-level systems. In this latter case the exponent is given by the ratio between the excitation energy and the potential barrier. In an attempt to elucidate between the two laws, we APPENDIX In this appendix we present aderivation of Eq. (2.2). We start from Mparticles without holes or bonds between them and let us insert Zholes and 8bonds with the condition that every two of them must be separated by at least one particle. It is convenient to analyze the six possible following cases separately. (a) Configurations having one particle to the left and another particle to the right. Its number is easily seen to be (M — 1)! Z!8!(M — Z8— 1)!— (Al) (b) Configurations having one hole to the left and, of course, ending by one particle. There are (M — 1)! (Z — 1)!8!(M — Z— 8)! (A2) of them. (c) Configurations having one hole to the right and, of course, beginning by one particle. This number is also given by Eq. (A2). (d) Configurations having one bond to the left and, of course, one particle to the right: (M — 1)! Z!(8 — 1)!(M — Z— 8)! (A3) (N — M8— 1)!— (Z — 1)!(N — M8— Z)!— (A4) (b) Now holes can also be introduced to the right of the lattice without changing the value of Z. The number of different ways of introducing the new holes is (e) Configurations having one bond to the right and one particle to the left are given by Eq. (A3). The next step is to introduce in each of the above configurations N— (M+8+Z) holes in all the possible ways, but without creating any new particle-hole contacts. This implies that the new holes can only be inserted next to one of the Zholes already present. In the process, the periodic boundary conditions must be taken into account. Let us consider each of the above group of configurations. (a) It is seen that the number of possibilities for every configuration in this group is 590 J.J.BREYAND M. J.RUIZ-MONTERO 43 (N — M— 8)! Z!(N — M— 8— Z)! (N — M— 8— l)! (Z — 1)!(N — M— 8— Z)! (A6) In this way we generate all the configurations having at least ahole to the left. (c) Only those final configurations with one particle to the left and one or more holes to the right have to be considered here. Then, for each of the starting configurations given by Eq. (A2) we have possibilities. (d) and (e) Equation (A4) also holds. The final step is to collect all the above results, namely multiply Eq. (Al) by Eq. (A4), Eq. (A2) by Eq. (AS), and so on, adding all the resulting products. This 1eads very easily to Eq. (2.2). 'C. A. Angell and M. Goldstein, Dynamical Aspects ofStructur al Change in Liquids and Glasses (Annals of the New York Academy of Sciences, New York, 1986). 2For arecent review see, for instance, G. H. Fredickson, Annu. Rev. Phys. Chem. 39, 181 (1988). See also, J. Jackie, Rep. Prog. Phys. 49, 171 (1986). An interesting and speculative discussion can be found in J. P. Sethna, Europhys. Lett. 6, 529 (1988). G. W. Scherer, Relaxation in Glasses and Composl'tes (Wiley, New York, 1986). 4S. A. Langer, J. P. Sethna, and E. R. Grannan, Phys. Rev. B 41, 2261 (1990). 5R. Schilling, J.Stat. Phys. 53, 1277 (1988). G. M. Bell, J.Math. Phys. 10, 1753 (1968). 7R. R. Singh and J.Kovac, J.Chem. Phys. 89, 1083 (1988). 8G. H. Fredickson and H. C. Andersen, Phys. Rev. Lett. 53, 1244 (1984); G. H. Fredickson and S. A. Brawer, J. Chem. Phys. 84, 3351 (1986). T. A. Weber, G. H. Fredickson, and F. H. Stillinger, Phys. Rev. B34, 7641 (1986); T. A. Weber and F. H. Fredickson, ibid. 36, 7043 (1987). Amore elaborated tiling model is discussed in J. G. Harris and F. H. Stillinger, ibid. 41, 519 (1990). ~oA. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comput. Phys. 17, 10 (1975). 'W. Kob and R.Schilling, Z. Phys. B68, 245 (1987). D. A. Huse and D. S.Fisher, Phys. Rev. Lett. 57, 2203 (1986).