Glasslike relaxation in a one-dimensional bonded-fluid model

PHYSICAL REVIEW BVOLUME 44, NUMBER 17 Glasslike relaxation in aone-dimensional bonded-fluid model 1NOVEMBER 1991-I J.J.Brey and M. J.Ruiz-Montero F&sica Teorica, Universidad de Seuilla, Apartado Correos 1065, Sector Sur, E-41080Seuilla, Spain (Received 5March 1991) AMonte Carlo investigation of the linear and nonlinear relaxation of aone-dimensional bonded fiuid is presented. The simulation shows that the model possesses many of the qualitative properties of the glass-forming liquids. The linear response function exhibits nonexponential relaxation, and it can be fitted quite accurately with the Kohlrausch-Williams-Watts expression. The average relaxation time has anon-Arrhenius temperature dependence. Simulations of temperature jumps have also been performed, and we find that in some cases the relaxation ofthe model is not monotonic. I. INTRODUCTION Glass-forming liquids display aquite characteristic relaxation behavior. Although the quantitative aspects seem to depend strongly on the details of the specific material under consideration, the general picture that appears from the experimental results corresponds to some kind of universal behavior. In fact, several general phenomenological theories have been proposed that have proved to be very successful in describing agreat variety of viscous Auids and glasses. ' From aphysical point of view, little is known about the relevant processes that control relaxation in glassforming liquids. Alarge number of models showing at least some of the observed experimental properties have been analyzed, but the physical mechanisms invoked to account for the glasslike behavior are very different, and even contradictory. 'Afundamental question in this respect is whether it is meaningful to talk about glasses, in general, or one has to consider each of them at an (almost) individual level. Lacking ageneral well-established theory, it seems worth while to keep looking for specific models, as simple as possible, that mimic the dynamical properties of real supercooled liquids and glasses. In particular, we believe that it is interesting to study models having the following two features. First, they must have well defined and consistent thermodynamic and kinetic properties. Second, they Inust be simple enough to allow arigorous investigation. Aprimary goal should be to identify the conditions required for amodel in order to be able to describe glassy dynamical properties. In this paper we use Monte Carlo simulation to discuss the linear and nonlinear relaxation properties of aonedimensional bonded Quid. It is alattice model whose thermodynamic properties were studied by Bell. Very recently, the model was also given kinetic properties by defining a master equation with transition probabilities that are consistent with the Hamiltonian. Previous Monte Carlo simulation has shown that the model presents, under continuous cooling, aphenomenon that is similar to the laboratory glass transition in real glasses. In Sec. II, the one-dimensional bonded Quid is defined. Since it has already been discussed elsewhere, 'only a brief description will be given. The same also applies to the technical details of the Monte Carlo simulation to be discussed at the beginning of Sec. III. The simulation is used to analyze the equilibrium time correlation function characterizing the response in energy to an infinitesimal perturbation in the limit of avery long wavelength. In Sec. IV, the nonequilibrium relaxation following afinite temperature jump is studied. Section Vcontains some comments and discussions. II. BOND-FLUID MODEL We consider acollection of Mparticles in aonedimensional lattice of Xsites. The interaction energy between two nearest-neighbor particles is — c. Two particles separated by asingle empty site may establish abond of energy — (a+co) between them. In this case it is convenient to think of the empty site as being occupied by a bond. No other interactions are possible in the system. To specify aconfiguration, the position of the particles on the lattice and also the existing bonds must be given. The energy of agiven configuration xof the lattice is, assuming periodic boundary conditions, E(x)=— (M — Z)E— Bco, where 8is the number of bonds and 2Z is the number of particle-hole contacts in the configuration. Ahole is an empty site that is not being occupied by abond. Since we are dealing with aone-dimensional model with shortrange interactions, the system will not present thermodynamic phase transition. The equilibrium properties of the system are easily derived by constructing the isothermal-isobaric partition function, 'or by using the transfer-matrix method. Their expressions will not be given here. The dynamics of the model is defined by means of amaster equation for the probability distribution, P(x,t), of configuration xat time t, aP(x,t)= y[W( 'x~ )xP(x', t) at — W( ~ xx) (Pxt)] . The transition rates W(x— +x ') from configuration xto 9259 1991 The American Physical Society 9260 J.J.BREY AND M. J.RUIZ-MONTERO configuration x' are defined in the following way. Aparticle can jump to anext empty site if and only if it is not bonded to another particle. Bonds between particles can be destroyed, and abond can also be created between two particles separated by asingle empty site. For the allowed transitions, the transition rate is given by 1AE fV(x~x ')=1— tanh 7QB(2.3) where 7"p is anatural attempt frequency, kB is the Boltzm ann constant, Tthe temperature, and b,E=E(x')— E(x). The above definition of the transition probabilities guarantees that the Markov process is irreducible and, therefore, any initial distribution P(x,0) will approach in time the canonical form P,„(x)~exp E(x) B(2.4) As the temperature is lowered, the number of holes and particle-particle contacts decreases. To relax aperturbation, the system needs to destroy an increasing nurnber of bonds before ahole finds acouple of particles, it jumps between them, and afterwards creates anew bond. It follows that we can expect our system to show a temperature-dependent activation energy. For comparison purposes it must be realized that the dynamics of the bonded-Quid model does not involve any kind of external disorder or any geometrical constraints restricting the number of allowed transitions. III. LINEAR RESPONSE The response in energy of the system to an infinitesimal homogeneous perturbation is characterized by the correlation function &E(t)E(0))— &E )' &E') — &E)' (3.1) where the angular brackets denote the equilibrium average at the given temperature. The function P(t) is defined such that P(0)=1. We have investigated the behavior of the response function P(t) of our model through Monte Carlo simulations with periodic boundary conditions. The details of the algorithm have been presented in Ref. 7. We started the simulations with the system in equilibrium at the temperature T=O .For asystem with number of sites twice the number of particles, there is aunique equilibrium configuration corresponding to such temperature: a structure in which every two particles are separated by a hole, without bonds between particles. From this state of zero energy, an instantaneous quench to agiven temperature was carried out. Afterwards, the system was allowed to relax to equilibrium before we start recording the energy E(t). The lattice used in the simulation had 200 sites and 100 particles, and the interaction energies were given values c, =200 and co =500, in arbitrary units. For all the ranges of temperatures investigated, 200~ kBT~700, the equilibrium correlation length, estimated by using the 500 1000 1500 FIG. 1. Energy correlation function at k&T=300 for three different sizes of the system. The curves are the average over 500 independent runs for N=200, over 50 for N=400, and over 40 for N= 1000. I LA C) tnt FICr. 2. KWW fit of the linear relaxation function for kz T=280. The solid line corresponds to values of the parameters ~„=42.87 and y=0.368. transfer-matrix method, is much smaller than the length of the system. Anyway, we have checked that the results did not change significantly as the size of the system is increased. As an example, Fig. 1shows the results obtained for P(t) on the systems of sizes N=200, 400, and 1000, at kBT=300. Due to computer time limitations, the number of independent runs performed was smaller the larger the system was. For N=200, the results were averaged over 500 runs, for N=400 the number of runs was 50, and 40 for N=1000. From the comparison of the three curves in Fig. 1it follows that the best compromise choice is %=200, because it gives enough accuracy and allows one to decrease the Auctuations by considering alarger number of independent runs. All the results to be presented in the following correspond to a system with 100particles, averaged over 500 trajectories. For each of the studied temperatures, the correlation function P(t) was fitted to the Kohlrausch-WilliamsWatts (KWW) function: GLASSLIKE RELAXATION IN AONE-DIMENSIONAL. . . 9261 (3.2) IIIII with rand ybeing adjustable parameters. Only the values in the interval 0.03 ~/(t) ~0.5were used in the fitting of the curves. This region accounts for most of the area under P(t). In Fig. 2we have plotted ln( — in') versus lnt for k~T=280. The solid line is the KWW fit with parameters ~=42.87 and y=0.368. In general, the agreement is better the lower the temperature, as can be seen from Figs. 3and 4, where the correlation functions at kz T=280 and 500 are shown. Also notice that systematic deviations from the KWW behavior are observed in the short-time region. The values of ~and yobtained from the fittings are given in Table I. From the values of yit follows that the correlation function has ahighly nonexponential decay, showing abroad relaxation spectrum. Besides, the temperature dependence of yimplies the lack of thermorheological simplicity of the model, i.e.,the temperature dependence of Pcannot be described by means of asingle time scaling. Nevertheless, it must be noticed that, for not too large temperature intervals, the model could give the appearance of thermorheological simplicity. Although there is aclear tendency of yto decrease as the temperature decreases, the dependence of yon the temperature does not appear to be monotonous in the low-temperature region. In this respect our results differ from those obtained by Fredrickson for facilitated Ising models, 'and agree with the behavior found by Stillinger and Weber for tiling models. ' We have also investigated the temperature dependence of the average relaxation time defined by ~= fdt P(t) .(3.3) 0 To evaluate this expression we have used the KWW expression with the values of the parameters obtained in the previous fit. The results are also listed in Table I. An expression that is often used to fit the experimental data for the temperature dependence of the average relax100 ~~nnnnqnnnnnnnnnnnn~ I~I 'I 50 200 FIG. 4. Linear relaxation function for k~ T=500. The circles are the results of the simulation and the solid line the KWW fit with r=4.17 and y=0.413. ation time in glasses is the so-called Vogel-TammanFulcher (VTF) form Eo p(3.4) d d(1/k~ T) (3.5) is seen to increase as the temperature decreases. This is also acharacteristic of many glass-forming liquids. " It must be noticed that the value of To obtained from the fit of Eq. (3.4) is well below the range of temperatures of the Monte Carlo simulation data. This prevents us from extrapolating the validity of Eq. (3.4) to the neighborhood of To, and to conclude that the model presents a kinetic singularity at kz T=74. In fact, the results to be discussed in the next section show that there is no such singularity. Therefore, the VTF expression must only be where Eo and To are adjustable parameters. Figure 5 shows afit of Eq. (3.4) to our data. The parameters are kg To =74 and k~Eo =1098. The curvature of the plot clearly indicates that the kinetics of the model does not obey the Arrhenius equation. More precisely, the apparent activation energy TABLE I. Values of the K%'W fitting parameters of the linear relaxation function for several values of the temperature. Also, the values of the average relaxation times are given. k, T C) 0500 1000 t1500 FIG. 3. Linear relaxation function for k&T=280. The circles are the results of the simulation and the solid line the best KWW fit with the values of the parameters given in Fig. 2. 650 550 500 400 350 300 280 250 0.479 0.429 0.413 0.370 0.356 0.345 0.368 0.377 3.1 3.48 4.17 7.15 11.81 25.11 42.87 117.20 6.73 9.64 12.73 29.928 56.05 132.82 182.47 462.63 9262 J.J.BREY AND M. J.RUIZ-MONTERO IV. RESPONSE TO INSTANTANEOUS QUENCHES We have carried out a series of experiments to study the relaxation of the system after an instantaneous quench from T=O to alow (positive) temperature. It must be kept in mind that, in this kind of model, negative temperatures correspond to states with alarger energy. The relaxation of the energy of the model can be characterized by the nonlinear response function E(t)— (E) 1E(t) E(0)— «) «) (4.1) 1.522.5(c+tu)y FIG. 5. Temperature dependence of the average relaxation time ofthe linear relaxation function. The solid line is aVTF fit with k~ To =74 and km Eo =1098. considered as auseful expression to describe the behavior of the average relaxation time in the range of temperatures we have studied. The temperature To can, at the most, be associated with achange in the qualitative behavior of ~. It is also interesting to attempt to fit our data to the expression proposed by Adam and Gibbs' 7=exp C (3.6) where E, is atemperature-independent parameter and S, is the configurational entropy. The latter can be identified in our model with the total equilibrium entropy that is easily computed. 'In Fig. 6the logarithm of the average relaxation time is plotted versus (TS,)'. The apparent curvature of the plot indicates that our model significantly deviates from the Adam-Gibbs expression. More will be said about this point in the last section. Here E(t) is the nonequilibrium average energy at time t, and (E)the equilibrium average energy corresponding to the temperature after the quench. The results we present in this section correspond again to alattice of 200 sites with 100 particles, and they have been averaged over 1000 independent trajectories. For some of the lowest temperatures considered, we have checked that the results are not significantly altered when the number of particles is increased to 1000 (2000 sites). From Eq. (4.1) it follows that g(0)=1, and that the relaxation of the system towards equilibrium implies 1t (~)=0. At all the temperatures investigated (SO~ k~T ~1000), it was observed that g(t) decays to zero for large enough times. This is important with regard to the discussion in the previous section. If there were adynamical transition at kz To — — 74, one should expect that the system would not be able to equilibrate below that temperature. Since that was not the case, we conclude that the singularity present in Eq. (3.4) is artificial and, therefore, that expression does not apply when one approaches To from above. Nevertheless, although the system always goes to equilibrium, important qualitative differences appear in the relaxation. For Ez T)150 the relaxation is not rnonotonous. There is afast initial decay of the average energy that is overshot below its equilibrium value. An example is shown in Fig. 7that corresponds to k&T=1000. The final relaxation takes place through negative values of g. This behavior can be easily explained in terms of the C3 2x10 3x10 1/yg 4x10 C FIG. 6. Plot of the average relaxation time vs (TS,)', where S, is the equilibrium configurational entropy at temperature T. If the Adam and Gibbs expression would apply, the points should be accurately fitted by astraight line. 50 100 150 FIG. 7. Nonmonotonous relaxation following an instantaneous quench from k& T=1000. GLASSLIKE RELAXATION IN AONE-DIMENSIONAL. . . 9263 elementary processes taking place in the system. Because of the initial configuration we have chosen, there is a great tendency to create bonds, as those are the most favorable transitions from the energetic point of view. In this way, the system builds up anumber of bonds larger than the equilibrium one. This has been verified by following the evolution of the configurations in the simulation. At lower temperatures (k~ T(150), the relaxation is monotonous, but highly nonexponential. There is an almost instantaneous decay of the response function to values of the order of 0.2. Afterwards, the relaxation is much slower. At the lowest temperatures studied, k~T=50, runs of approximately 3X10 MCS were required for P(t) to decay to values of the order of 0.05. Following asimilar study by Fredrickson and Brawer' for afacilitated Ising model, we have fitted the relaxation data at low temperatures to aKWW function r' g(t) =exp (4.2) Of course, the meaning of these fits is quite different from the ones carried out for the linear relaxation function. Here, Eq. (4.2) must be understood only as aconvenient way of analyzing the results. In obtaining these fits, all the data points for P(t) were used. The fits turn out to be reasonably good beyond the initial fast decay. In Fig. 8we have plotted the simulation data and the fit for the lowest studied temperature, kz T=50. The values of the parameters obtained are y' =0.16 and z=9.0X 10, showing that the relaxation cannot be deqQ scribed by asingle exponential. In Table II we present the values of y' and ~obtained for all the low temperatures (ks T~ 120) investigated. Also, the average relaxation times, defined by the area under g(t), are given. Comparison of Tables Iand II shows that the exponents y' are smaller than the yones, while the relaxation times 7are much larger than the 7 ones. Although it can be associated to the fact that the TABLE II. Values of the KWW fitting parameters and of the average relaxation times following an instantaneous quench from k& T=O k, r 120 100 90 80 70 65 50 0.163 0.146 0.139 0.137 0.138 0.139 0.134 4.306 9.52 18.03 48.58 210.15 486.41 8.99x10' Tq 4X 10 3.55 x10' 1.35 x10' 4.50x10' 1.75 X10 3.64 X10 1.17x10' temperatures are smaller in Table II than in Table I, the different physical meaning ofboth sets of relaxation times must be noticed. The exponent y' seems to depend very weakly on the temperature, although it appears to decrease slightly as the temperature is lowered. On the other hand, rand wq are strongly temperature dependent, and show afast increase as the temperature decreases. Figure 9is an rrhenius plot of the average relaxation times. Althoug the data points cannot be accurately fitted by astraig t line atendency towards Arrhenius behavior is observed line, aten ency owa for temperatures below k&T=90. This reminds the behavior of some fragile Auids, for which areturn to an Arrhenius temperature dependence of the relaxation time has been observed, when approaching the laboratory glass transition temperatures from above. V. DISCUSSION The results presented in this paper and in Ref. 7show that the one-dimensional bonded-Quid model presents many of the dynamical properties of real glass-forming liquids. It has been found that the linear relaxation function for the model has ahighly nonexponential decay that can be fitted reasonably well with the KWW expression. The values of the KWW exponent are of the same order as in real glasses, and they show atendency to decrease as C3 05x10 10 I 1.5x10 7 2x10 2.5x10 t3x10 0 FIG. 8. KWW fit of the nonlinear response function following an instantaneous quench from k& T=O kT=O to 50. The values of the fitting parameters are given in Table II. 0.01 0.0151ykT 0.02 FIG. 9. Arrhenius plot of the average relaxation times after an instantaneous quench from k~ T=O J.J.BREY AND M. J.RUIZ-MONTERO the temperature is lowered. The average relaxation times increase very fast when one approaches the lowtemperature region, and their behavior can be accurately described by the empirical VTF law in all the ranges of temperatures we have studied. Apoint that deserves some comments is the validity of the Adam-Gibbs relation, Eq. (3.6), for this kind of lattice model. Once one identifies the structural entropy with the equilibrium entropy of the system, Eq. (3.6) establishes arelation between an equilibrium property and a pure dynamical property, the average relaxation time. Nevertheless, the statics of the model does not determine its dynamics when the latter is formulated in terms of a master equation. There are many difFerent ways of defining aset of transition probabilities that are compatible with agiven equilibrium distribution, even if detailed balance is required. Of course, the situation is not the same when acomplete and consistent Hamiltonian description of the system is used. From the above discussion it follows that, unless some additional condition is introduced in the definition of the transition rates, it is amatter of chance whether the model will obey the Adam-Gibbs relation or not. Therefore, we believe that nothing physically deep is behind the fact that the facilitated Ising models accurately obey the Adam-Gibbs relation, whereas the tiling models' and the bonded-Quid model show aclear departure therefrom. An important feature of our model is the nonmonotonous relaxation after some instantaneous jumps of temperature. At this moment we do not know whether this is only an artificial property of our model or it is related to some actual behavior of real glasses. We have now several models that successfully describe many of the characteristics of real glass-forming liquids. These models allow quite eKcient computer simulations and, therefore, can be used to test approximate methods. Besides, they prove that very difFerent mechanisms are able to produce glassy behavior and, therefore, one of the key questions remains. Can we really talk about ageneral kind of systems named glasses? ACKNOWLEDGMENTS We acknowledge partial support from the Direccion General de Investigacion Cienttfica yTecnica (Spain) through Grant No. PB89-0618. O. S.Narayanaswamy, J.Am. Ceram. Soc.54, 491 (1971). G. W. Scherer, Relaxation in Glasses and Composites (Wiley, New York, 1986). K. L. Ngai, R. W. Rendell, A. K. Rajagopol, and S. Teitler, Ann. N.Y.Acad. Sci.484, 150 (1986). 4G. H. Fredri. ckson, Annu. Rev. Phys. Chem. 39, 181 (1988). 5A. K. Rajagopol, K. L. Ngai, R. W. Rendell, and S. Teitler, Physica A149 358 (1988). G. M. Bell, J.Math. Phys. 10 1753 (1968). 7J.J.Brey and M. J.Ruiz-Montero, Phys. Rev. B43, 585 (1991). 8N. G. Van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981). G. H. Fredrickson, Ann. N.Y.Acad. Sci. 484, 185 (1986). T. A. Weber and F. H. Stillinger, Phys. Rev. B36, 7043 (1987). S.Brawer, J.Chem. Phys. 81,954 (1984). G. Adam and J.H. Gibbs, J.Chem. Phys. 43, 139 (1965). G. H. Fredrickson and S.A. Brawer, J.Chem. Phys. 84, 3351 (1986). i4T. A. Weber, G. H. Fredrickson, and F. H. Stillinger, Phys. Rev. B34, 7641 (1986).