A simplified nonlinear memory function for the dynamics of glass-forming materials based on time-convolutionless mode-coupling theory
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n A simplified nonlinear memory functionfor the dynamics of glass-forming materialsbased on time-convolutionless mode-coupling theory
Michio Tokuyama
Institute of Multidisciplinary Research for Advanced Materials, Tohoku University, Sendai980-8577, Japan
Abstract
A simplified nonlinear memory function is proposed in the ideal time-convolutionlessmode-coupling theory equation to study the dynamics of glass-forming liquids. The nu-merical solutions are then compared with the simulation results performed on fragile liq-uids and strong liquids. They are shown to recover the simulation results in a supercooledstate well within error, except at a β -relaxation stage because of the ideal equation. A tem-perature dependence of the nonlinearity µ in the memory function then suggests that thesupercooled state must be clearly separated into two substates, a weakly supercooled statein which µ increases rapidly as T decreases and a deeply supercooled state in which µ becomes constant up to the glass transition as T decreases. On the other hand, it is shownthat in a glass state µ increases rapidly as T decreases, while it is constant in a liquid state.Thus, it is emphasized that the new model for the simplified memory function is muchmore reasonable than the conventional one proposed earlier by the present author not onlyqualitatively but also quantitatively. Key words:
Glass transition; Simplified nonlinear memory function; Supercooled liquids;Time-convolutionless mode-coupling theory; Universality
Understanding of the glass transition from first principles is one of the importantworks left in condensed matters physics. A considerable attention has been drawnto study the dynamics of glass-forming materials near the glass transition experi-mentally and theoretically [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15].In order to study the glass transition from a statistical-mechanical point of view, wehave recently proposed the time-convolutionless mode-coupling theory (TMCT)
Preprint submitted toElsevier 22January 2020
14] and then derived the second-order differential equation for the cumulant func-tion K ( q, t )(= − ln[ f ( q, t )]) [16], where f ( q, t ) is the scaled collective-intermediatescattering function. We call this an ideal TMCT equation. This equation containsthe nonlinear memory function ∆ ϕ ( q, t ) , which has the same form as that obtainedin the mode-coupling theory (MCT) [1,2,3]. The memory function is written interms of the static structure factor S ( k ) and is integrated over the wavevector k .Once the analytic form of S ( k ) is known, the ideal TMCT equation is easily solvednumerically. Here it is well-known that depending on temperature T , there existthree characteristic states in the glass-forming liquids, a liquid state [L], a super-cooled state [S], and a glass state [G]. It is also well-known that depending on atime scale, there are four characteristic stages in the dynamics of supercooled liq-uids, a first stage for a short time, a β -relaxation stage for an intermediate timewhere f ( q, t ) is described by the von Schweilder decay [17], an α -relaxation stagefor a later time where f ( q, t ) is well-known to obey the Kohlrausch-Williams-Watts(KWW) function (or a stretched exponential decay) [18,19], and a diffusion stagefor a long time where f ( q, t ) obeys a diffusion equation. In order to check whetherthe numerical solutions of the ideal TMCT equation obey such relaxation processesor not, we have recently solved it numerically for hard spheres [20,21,22] by us-ing the Percus-Yevick static structure factor [23]. Thus, it has been shown that thecritical volume fraction φ c ≃ . agrees with that predicted from the molecular-dynamics simulations [24,25,26,27] and also that the numerical solutions describethe simulation results well, except at β stage because of an ideal equation [22,28].Since S ( k ) is in general not known, however, one has to obtain its numerical valuesfrom either the simulations or the experiments. This is not an easy task to do. There-fore, the main purpose of the present paper is to propose a reasonably simplifiednonlinear memory function and thus to check whether the numerical solutions ofthe ideal TMCT equation with such a memory function can describe the simulationresults performed on (F) fragile liquids and (S) strong liquids well or not.In the previous paper [16], we have first simplified the memory function ∆ ϕ ( q, t ) around q = q m as ∆ ϕ ( q m , t ) = ∆ ϕ ( q m , f ( q m , t ) by just employing the same ap-proach as that used in MCT [2], where q m is a first peak position of S ( q ) . Then, wehave shown that the numerical solutions can roughly describe the simulation resultsonly in [L]. In order to describe the dynamics of supercooled liquids, therefore, wehave next proposed the renormalized simplified memory function [29] ∆ ϕ ( q m , t ) = ∆ ϕ ( q m , f ( q m , t ) ε = ∆ ϕ ( q m , e − εK ( q m ,t ) , (1)where ε is a nonlinear exponent to be determined. We call this a conventionalmodel. Thus, it has been shown that ε is constant to be 2.0 in [L], while it growsfrom 2.0 in [S] as temperature decreases. The numerical solutions have been shownto describe the simulation results performed on two types of glass-forming liquids(F) and (S) qualitatively but not quantitatively. Especially, this model fails in de-scribing the simulation results at α stage near the glass transition quantitatively. Infact, the stretched (or KWW) exponent β obtained from the numerical solutions2 ig. 1. (Color online) A plot of the nonlinear parameters µ and ε versus scaled inversetemperature T c /T for (a) a new simplified model given by Eq. (2) and (b) a conventionalmodel given by Eq. (1), where T c is a critical temperature. The symbols ( • ) indicate thefitting values for a fragile liquid and ( ✷ ) for a strong liquid. The symbols ( ✸ ) indicate theglass transition point T g , ( ⊙ ) the supercooled point T s , and ( + ) the deeply supercooledpoint T f , where the numerical values of T i are listed in Table 2. The solid lines are guidesfor eyes. The label [L] stands for a liquid state, [S w ] for a weakly supercooled state, [S f ]for a deeply supercooled state, and [G] for a glass state. does not coincide with that obtained from the simulation results near the glass tran-sition [30].In order to simplify the memory function, we have assumed in the previous papers[16,29] that the dominant contribution in ∆ ϕ ( q, t ) results from the first peak q m of S ( k ) and then simplified S ( k ) in terms of δ function as S ( k ) = 1 + Aδ ( k − q m ) [2], where A is a positive constant to be determined. However, the disagreementbetween the numerical solutions and the simulation results at α stage suggests thatthe integration over the wavevector k contained in ∆ ϕ ( q, t ) must plays an importantrole in finding a reasonable value of the KWW exponent β at α stage. Consideringsuch a wavevector dependence, therefore, we now propose the following simplifiednonlinear memory function: ∆ ϕ ( q m , t ) = ∆ ϕ ( q m , f ( q m , t )[1 + µK ( q m , t )] β/b , (2)where µ is a nonlinear parameter to be determined and b the von Schweidler expo-nent [17]. We call this a new model. In Ref. [28], we have pointed out theoreticallythat the solutions of the ideal TMCT equation with the original (non-simplified)memory function deviate from the simulation results at β stage because of the ap-proximation employed (see next section for details). In fact, this has been checkednumerically by using the PY static structure factor [22,28]. Hence we note that thecomparison of the numerical solutions with the simulation results must be done ex-3ept at β stage. Because of this reason, in this paper we solve the TMCT equationnumerically not only for the new model given by Eq. (2) but also for the con-ventional model given by Eq. (1) under the same initial conditions obtained fromthe simulation results and compare both solutions together with the simulation re-sults consistently. Thus, it turns out from the new model that the supercooled stateshould be further separated into two substates, a weakly supercooled state [S w ]and a deeply supercooled state [S f ]. As is shown in Fig. 1(a), this is because astemperature decreases, the nonlinear parameter µ becomes constant in a deeplysupercooled state [S f ], while it grows rapidly in [S w ] and [G], where in [L] µ isconstant to be 1.150 for (F) and 1.046 for (S). Similar behavior is seen for twotypes of glass-forming liquids, (F) and (S). It is also shown that the new modelenables us to recover the same value of β as that obtained for the simulation resultsnear the glass transition [30]. On the other hand, the conventional model shows inFig. 1(b) that the nonlinear exponent ε grows rapidly in [S] and [G] as temperaturedecreases, while ε is constant to be 2.564 in [L].In Section 2 we summarize the basic equations briefly, which are used in the presentpaper. In Section 3 we derive the universal TMCT equation, which has only one so-lution for different glass-forming liquids of type i , where i = (F) and (S). In Section4 we first briefly summarize the simulation results performed on two types of glass-forming liquids, (F) and (S). Then, we solve the universal TMCT equation for thenew model numerically under the initial conditions obtained from the simulationresults. Thus, we show that the numerical solutions agree with the simulation resultswell within error, except at β stage. We also compare them with those obtained forthe conventional model under the same initial conditions consistently. We concludein Section 5 with a summary. In this section, we briefly summarize the macroscopic equations for the scaledcollective-intermediate scattering function f ( q, t ) , which have been derived fromfirst principles by using a new formulation based on TMCT [14,16,28,29].We consider the three-dimensional equilibrium glass-forming system, which con-sists of N particles with mass m and diameter σ in the total volume V at tempera-ture T . Let ξ denote the control parameter, such as a volume fraction φ (= πρσ / and an inverse temperature /T , where ρ (= N/V ) is the number density. Then, thescaled collective-intermediate scattering function f ( q, t ) is given by f ( q, t ) = h ρ ( q , t ) ρ ( q , ∗ i /S ( q ) (3)with the collective-density fluctuation ρ ( q , t ) = N − / [ P Nj =1 e i q · X j ( t ) − N δ q , ] ,where X j ( t ) indicates the position vector of a j th particle at time t , the brackets4 ig. 2. (Color online) Classification of the basic equations discussed in the present paperinto three stages, [N], [K], and [H], depending on a space ( r )-time ( t ) scale, where ℓ i and τ i indicate the relevant length and time of interest, respectively. the average over an equilibrium ensemble, S ( q )(= h| ρ ( q , | i ) the static structurefactor, q = | q | , and f ( q,
0) = 1 . Depending on a space-time scale, there exist three characteristic stages, a micro-scopic stage [N], a kinetic stage [K], and a hydrodynamic stage [H] (see Fig. 2). Ina microscopic stage [N], the position X j ( t ) and the momentum P j ( t ) of j th particleat time t are described by the Newton (or Heisenberg) equations. Depending on aspace-time scale, the kinetic stage [K] further consists of two substages, a fast stage[K ] and a slow stage [K ]. In a fast stage [K ], the relevant variables are given bythe number density ρ ( q , t ) and the current density j ( q , t )( ∝ dρ ( q , t ) /dt ) . As dis-cussed in the previous papers [14,16,28,29], it is important to employ two typesof projection-operator methods to derive the basic equations for different relevantvariables. In fact, it is indispensable to use the time-convolutionless formalism forthe number densities to recover the cumulant expansion proposed by Kubo [31]. Alinear non-Markov stochastic diffusion equation for ρ ( q , t ) is then derived fromthe Heisenberg equation by employing the Tokuyama-Mori projection-operatormethod [33,34] (see a dotted arrow (T) in Fig. 2), where the memory terms areconvolutionless in time and are written in terms of correlation function of the fluc-tuating current. A linear non-Markov Langevin type equation for j ( q , t ) is also de-5ived from the Heisenberg equation by using the Mori projection-operator method[32] (see a dashed arrow (M) in Fig. 2), where the memory term is convolution intime and is written in terms of correlation function of the fluctuating forces. Byintroducing the cumulant function K ( q, t ) by K ( q, t ) = − ln[ f ( q, t )] , those cou-pled equations are then used to obtain the following closed nonlinear non-Markovsecond-order differential equation for K ( q, t ) (see a bold T arrow in Fig. 2): [28] ∂ ∂t K ( q, t ) = q v th S ( q ) − ζ ∂∂t K ( q, t ) − Z t ds Z s dτ ∆ ϕ ( q, s − τ ) f ( q, τ ) f ( q, s ) ∂ ∂τ K ( q, τ ) (4)with the nonlinear memory function ∆ ϕ ( q, t ) = ρv th Z < d k (2 π ) v ( q , k ) S ( k ) S ( | q − k | ) f ( k, t ) f ( | q − k | , t ) , (5)where ζ is a friction coefficient, R < the sum over wave vectors k whose magnitudesare smaller than a cutoff q c [22], and v th = ( k B T /m ) / . The vertex amplitude v ( q , k ) is given by v ( q , k ) = ˆ q · k c ( k ) + ˆ q · ( q − k ) c ( | q − k | ) , (6)where ρc ( k ) = 1 − /S ( k ) and ˆ q = q /q . This is a TMCT equation to discussthe dynamics of glass-forming equilibrium liquids. Here we note that the nonlin-ear memory function ∆ ϕ ( q, t ) has exactly the same form as that obtained in MCT[1,2,3], although TMCT equation is quite different from MCT equation. By usingthe Percus-Yevic static structure factor [23], one can solve this equation numeri-cally and show that it describes not only the α -relaxation process but also the β -relaxation process consistently [28]. However, we should mention here that becauseof the double time integrals in the memory function, such numerical calculationsare limited only for smaller volume fractions far from the glass transition. In gen-eral, therefore, it is difficult to obtain numerical solutions for higher values of ξ near the glass transition. In a slow stage [K ], however, one can avoid such a diffi-culty because one can apply the approximation (A2) given by f ( q, τ ) /f ( q, s ) ≃ [28] to Eq. (4) (see a upper bold arrow in Fig. 2). Thus, one can find the asymptoticsecond-order differential equation for K ( q, t ) ∂ K ( q, t ) ∂t = q v th S ( q ) − ζ ∂K ( q, t ) ∂t − Z t ∆ ϕ ( q, t − s ) ∂K ( q, s ) ∂s ds, (7)where the initial conditions are given by K α ( q, t = 0) = dK α ( q, t ) /dt | t =0 = 0 .This is an ideal TMCT equation to discuss the dynamics of glass-forming equilib-rium liquids. Here we note that because of the approximation (A2), the numericalsolutions of Eq. (7) coincide with those of Eq. (4), except at β stage.6n a hydrodynamic stage [H], the diffusion equation for K ( q, t ) is derived fromeither Eq. (4) or Eq. (7) in the long time limit as (see a lower bold arrow in Fig. 2) ∂K ( q, t ) ∂t = q D c ( q ) , or ∂f ( q, t ) ∂t = − q D c ( q ) f ( q, t ) (8)with the q -dependent collective diffusion coefficient D c ( q ) = v th /S ( q ) ζ + R ∞ ∆ ϕ ( q, s ) ds . (9) We now discuss an ergodic to non-ergodic transition at a critical point ξ = ξ c , abovewhich the long-time solution f ( q, t ) (or K ( q, t ) ) reduces to a non-zero value f c ( q ) (or K c ( q ) ), the so-called nonergodicity parameter. As shown in the previous papers[16,28,29], starting from either Eq. (4) or Eq. (7), in the long-time t → ∞ limit,one can find K c ( q ) = 1 F ( q ) (10)with the long-time limit of the nonlinear memory function F ( q, f c , f c ) = 12 Z < d k (2 π ) Θ( q , k ) f c ( k ) f c ( | q − k | ) , (11)where the vertex Θ is given by Θ( q , k ) = ( ρ/q ) v ( q , k ) S ( q ) S ( k ) S ( | q − k | ) . (12)The existence of the critical point in Eq. (10) is mathematically confirmed since K c ( q ) is a kind of the Lambert W-function. In fact, the critical volume fraction φ c has been obtained by solving numerically Eq. (10) with the Percus-Yevick staticstructure factor as φ c ≃ . [20,21,22], which agrees with that obtained fromthe simulation results performed on monodisperse hard spheres [24,25,26,27] withinerror. In the long-time limit we obtain K ( q, t ) = q D c ( q ) t for ξ < ξ c , while we find K ( q, t ) = K c ( q ) and D c ( q ) = 0 for ξ ≥ ξ c since the integral term of the memoryfunction ∆ ϕ ( q, s ) in Eq. (9) diverges at the critical point ξ c . Thus, the existenceof the nonergodicity parameter suggests a singularity as a function of the distance (1 − ξ/ξ c ) at ξ = ξ c . In fact, by employing the same mathematical approach as thatused in MCT [2,3], one can also write D c ( q ) in terms of a singular function as D c ( q, ξ ) ∝ (1 − ξ/ξ c ) γ , (13)where γ is a power exponent to be determined.7 .3 Characteristic decays We briefly review two types of characteristic decays near the glass transition, whichare used in this paper. One is the so-called von Schweidler (VS) decay [17] at theslow β stage given by f ( q, t ) ≃ f c ( q )[1 − ( t/t b ) b ] , (14)where b is a time exponent and t b a characteristic time. Another is the so-calledstretched exponential decay (or KWW function) [18,19] at α stage given by f ( q, t ) = f c ( q ) exp[ − ( t/t α ) β ] , (15)where β is the KWW exponent and t α the α -relaxation time. As shown in Ref. [30],near the critical point ξ c , the crossover time t x from the VS decay to the KWWdecay is given by t x ≃ K /βc t α ≃ K /bc t b . By using Eq. (7), one can also find asimple relation between the exponent parameter λ [3], b , and β as [30] λ = Γ[1 + b ] Γ[1 + 2 b ] = Γ[2 β − β + 1]Γ[ β + 1]Γ[ − β − , (16)where Γ[ x ] is an usual Γ -function. In the following, these results are used to discusshow the new model is different from the conventional model not only qualitativelybut also quantitatively. In this subsection, we discuss the universal equation to describe the dynamics ofglass-forming liquids of type i with the intensive control parameters such as an in-verse temperature, where i = (F) for fragile liquids and i = (S) for strong liquids.As shown in the previous papers [35,36,37,38,39,40,41,42], the mean- n th displace-ment in a liquid of type i coincides with the other mean- n th displacements in dif-ferent liquids of type i if the universal parameter u s in those liquids has the samevalue, where u s ( q = 0) = − log( D s ( q = 0) q m /v th ) . Here D s ( q = 0) is a long-time self-diffusion coefficient. However, we note that the displacement in (F) nevercoincides with that in (S), even if u s has the same value. In general, these situationsalso hold for the cumulant function K ( t )(= K ( q m , t )) . We discuss this next.Since the characteristic features between different liquids of type i can be comparedat a first peak of the static structure factor, we set q = q m in Eq. (7). Then, Eq. (7)depends on the physical quantities q m , v th , and S ( q m ) . In order to eliminate themfrom Eq. (7), it is convenient to introduce the time scale τ D as τ D = S ( q m ) / / ( q m v th ) . (17)8y using a scaled time τ = t/τ D , one can then find the following asymptoticsolution of Eq. (7): K ( τ ) ≃ τ / , τ ≪ ,τ /τ L , τ > τ α , (18)with the diffusion time τ L (= D − ) , where D ( q m ) is the scaled diffusion coefficientgiven by D ( q m ) = q m τ D D c ( q m ) = q m S ( q m ) / D c ( q m ) /v th . (19)For a short time the ballistic motion dominates the dynamics of the system, whilefor a long time the diffusion process dominates it. For such time regions there areno difference between the dynamics of (F) and that of (S). The clear differencebetween them appears at α and β stages through the nonlinear memory function ∆ ϕ ( q m , τ ) . One can then transform Eq. (7) into a dimensionless equation ∂ K ( τ ) ∂τ = 1 − ζ ∂K ( τ ) ∂τ − κ Z t M ( τ − s ) ∂K ( s ) ∂s ds, (20)where ζ = ζ τ D and M ( τ ) = ∆ ϕ ( q m , τ ) / ∆ ϕ ( q m , . Here the coupling parameter κ is given by κ = τ D ∆ ϕ ( q m , . This is an universal equation to describe thedynamic of glass-forming equilibrium liquids of type i .In order to solve Eq. (20) numerically, one needs to fix the values of two unknownparameters ζ and κ consistently. As shown in the previous papers [16,29], this isdone by using the simulation results at each value of D ( T ) . In fact, the frictioncoefficient ζ has been found to be constant from the short-time behavior of thesimulation results [29]. Since ζ depends on T through τ D , one can fix the value of ζ at each value of D . The coupling parameter κ is also found by using the simulationresults under the fixed value of D . In fact, from Eqs. (9) and (19), one can obtain κ = 1 − ζ DD R ∞ M ( τ ) dτ . (21)Thus, the coupled equations (20) and (21) can be solved self-consistently at a givenvalue of D . Hence the universal parameter u for K ( τ ) is now given by u = − log( D ( T )) . (22)Later, we also show that the scaled α -relaxation time τ α is related to D as τ α ∝ D − . In the previous paper [29], we have fixed a value of ε so that the plateau height ofthe numerical solution coincides with that of the simulation results. However, thisapproach is not appropriate because the numerical solutions of Eq. (7) with the PY9tatic structure factor are known to deviate from the simulation results at β stage[22,28]. In order to solve Eq. (20) with the simplified memory function, therefore,we start from the simulation results as an initial condition. Then, it is convenient touse the following formal solution of Eq. (20): K ( τ ) = K ( τ ) − R τ Σ( s )[ K ( τ − s ) − K ( τ )] ds R τ Σ( s ) ds (23)with Σ( s ) = κ Z s e − ζ ( s − s ′ ) M ( s ′ ) ds ′ , (24)where K ( τ ) = ( ζ τ − e − ζτ ) /ζ . This is a starting equation to find the numericalsolutions by iterations. The nonlinear exponents µ and ε are then found so that start-ing from the simulation results, the first iteration result coincides with the plateauheight of the simulation results around the crossover time τ x . Several iterations arethus done at the fixed values of µ and ε until κ reaches a constant value within errorof order − . In the following, we thus obtain the numerical solutions not only forthe new model but also for the conventional one and compare them together withthe simulation results from a unified point of view consistently. Finally, we notethat Eq. (23) is an exact equation derived from Eq. (20) and is different from theapproximate equation derived in the previous paper [29], where K ( τ − s ) has beenexpanded in powers of s/τ , up to order ( s/τ ) . The numerical solutions of Eq. (23) for the conventional model are shown to failin describing the α -relaxation process near the glass transition. This means thatthe integration over the wavevector k in the nonlinear memory function ∆ ϕ ( q m , t ) must play an important role in finding a reasonable value of β at α stage. In orderto perform the integration approximately, one may simply assume that S ( k ) obeysa Gaussian distribution with the peak position q m as S ( k ) ≃ a exp[ − a ( k − q m ) ] ,where a and a are positive constants. We also assume around k = q m that K ( k, τ ) ≃ ( k/q m ) K ( q m , τ ) and K ( | q m − k | , τ ) ≃ K ( q m , τ ) . Then, the inte-gration of Eq. (5) over k leads to an asymptotic function of (1 + µK ) − , where µ is a fitting parameter to be determined. Thus, we propose the following simplifiednonlinear memory function: M ( τ ) = f ( τ )[1 + µK ( τ )] ν , (25)where ν is a fitting parameter to be determined.This model contains two unknown parameters µ and ν . Since µK ( τ ) = − ln( f ( τ ) µ ) ,the parameter µ must play an role of a nonlinear exponent in f ( τ ) . On the other10 able 1Time exponents b and β , exponent parameter λ , and exponent ν for different systems [30].System b β λ ν (= β/b ) control parameter /T Fragile liquids 0.3064 0.6832 0.8980 2.230Strong liquids 0.2779 0.6809 0.9130 2.451control parameter φ Hard spheres 0.6051 0.710 0.7215 1.173Soft spheres 0.6051 0.710 0.7215 1.173 hand, the parameter ν must be found through the dynamics of α - and β -relaxationprocesses near the glass transition. In fact, at a slow β stage, f ( τ ) is described bythe von Schweidler decay given by Eq. (14) (or K ( τ ) ≃ K c + ( τ /τ b ) b ). As dis-cussed in Ref. [30], the VS decay dominates the system, up to the crossover time τ x . In α stage for τ ≥ τ x , therefore, the memory function M ( s ) in the time integralof Eq. (23) is dominated by the VS decay K ( s ) ∼ s b . From Eq. (25), we then obtain M ( s ) ∼ f c s − νb . On the time scale of order τ ( ≥ τ x ≥ s ) , one can expand K ( τ − s ) in powers of s/τ as K ( τ − s ) ≃ K ( τ ) + O ( s/τ ) . Thus, Eq. (23) reduces to K ( τ ) ∼ τ R τ s − νb ds ∼ τ νb . (26)This must be compared with the KWW function given by K ( τ ) ∼ ( τ /τ α ) β to find ν = β/b. (27)Here the numerical value of ν is found by using the numerical values of β and b obtained in Ref. [30], which are listed in Table 1.In order to check whether the new model is reasonable or not, one may use thescaled memory function ∆ ϕ ( q m , t ) / ∆ ϕ ( q m , directly obtained from the numer-ical calculation of the ideal TMCT equation given by Eq. (7) with the PY staticstructure factor [43]. In Fig. 3, the calculated scaled memory function is plottedversus (a) log( K ) and (b) log( t ) together with the new model given by Eq. (25) andthe conventional model given by Eq. (1). The calculated memory function is thenshown to be well described by the new model at ν = 1 . and µ = 162 . withinerror. Thus, the new model is expected to describe the simulation results well notonly qualitatively but also quantitatively. On the other hand, the conventional modelcan describe the calculated memory function only for a short time ( t < − ) anda long time ( t > t α ). In fact, one can write Eq. (25) as M ( t ) ≃ exp[ − (1 + µν ) K ] for a short time. From Eq. (1), one obtains ε = 1 + µν . The conventional model isthen shown to coincide with the calculated memory function only for a short and along times at ε = 191 . (see Fig. 3), where ν = 1 . and µ = 162 . . Thus, theconventional model is easily shown to disagree with the calculated memory func-11 ig. 3. (Color online) A plot of the scaled memory function ∆ ϕ ( q m , t ) / ∆ ϕ ( q m , versus(a) log( K ) and (b) log( t ) . The solid line indicates the scaled memory function obtainedfrom the numerical calculation of Eq. (7) at φ = 0 . , q m σ = 7 . , and q c σ = 60 [43]. Thelong-dashed line indicates the numerical values of Eq. (25) at ν = 1 . and µ = 162 . .The dashed line indicates the numerical value of Eq. (1) at ε = 191 . . The symbols ( • ) indicate the crossover time t x = 10 . ( K x = 10 − . ) and ( ✷ ) the α -relaxation time t α = 10 . ( K α = 1 ). tion at α and β stages even for any other values of ε . This is the main reason why itfails to describe the α -relaxation process. In the present section, we solve Eq. (23) based on the simulation results.
We first briefly summarize the simulation results. As shown in the previous papers[40,41,42], the molecular-dynamics simulations have been performed for two typesof glass-forming liquids, (F) fragile liquids and (S) strong liquids, where the con-trol parameter ξ is given by an inverse temperature /T . By using those simulationresults, one can then find the physical quantities, such as q m and ζ . As typical ex-amples, we here take the binary mixtures A B with the Stillinger-Weber poten-tial (SW) [44] for (F) and SiO with the Nakano-Vashishta potential (NV) [45] for(S). The detailed information about those simulations is found in Refs. [40,41,42].In SW, length, time, and temperature are scaled by σ , t (= σ/v ) , and ǫ/k B , re-spectively, where ǫ is energy unit, and v = ( ǫ/m ) / . In NV, length, time, andtemperature are measured by ˚A, ps, and K, respectively. The particle mass is givenby average mass over components. The friction constant ζ is then found from the12 ig. 4. (Color online) A log plot of the scaled diffusion coefficient D versus inverse temper-ature for (a) a fragile liquid and (b) a strong liquid. The symbols ( • ) indicate the simulationresults. The solid lines indicate the singular function given by Eq. (28). The symbols ( ✸ )indicate the glass transition point T g , ( ⊙ ) the supercooled point T s , and ( + ) the deeplysupercooled point T f , whose numerical values are listed in Table 2. simulations as ζ ( T ) ≃
12 (SW) and 96 (NV), while the peak position q m of thetotal static structure factor S ( q ) is chosen near the glass transition as q m σ = D and ζ are found from the simulation results. In Fig. 4, thesimulation results for D ( q m ) are then plotted versus inverse temperature togetherwith the singular function of D . As discussed in Refs. [29,39], the value of thepower exponent γ of D c depends on type i . In fact, we have γ = 4 . for (F) and4.333 for (S) [39]. Similarly to Eq. (13), one may also assume that D obeys thesingular function near T c as D ( T ) = A r (1 − T c /T ) γ r , (28)where γ r is a power exponent to be determined and A r a prefactor to be determined.From the simulation results, we thus find γ r ≃ . and A r ≃ . for (F) and γ r ≃ . and A r ≃ . for (S). Here we note that γ r is slightly differentfrom γ because D also depends on S ( q m ) and v th through Eq. (19). By using the Table 2Critical temperature T c , glass transition temperature T g , supercooled temperature T s ,deeply supercooled temperature T f , and corresponding universal parameters u g = u ( T = T g ) , u s = u ( T = T s ) , and u f = u ( T = T f ) for different systems.type system T c T g T f T s u g u f u s fragile SW 0.505 0.527 0.644 0.891 6.082 3.102 1.705strong NV 2490.4 2826.5 3317.9 3830.0 4.09 0 2.620 1.920 ig. 5. (Color online) Universal function of u . (a) A plot of the scaled friction coefficient ζ versus u . The symbols ( ✷ ) indicate the numerical values for SW, ( ⊙ ) for Al O , ( • ) forNV, and ( + ) for BKS. The solid lines are guides for eyes. (b) A plot of f ( q m , τ ) versus log( τ ) at u ≃ . and ζ ≃ . . The solid line indicates the simulation results for SW at T = 0 . and the dashed line for Al O at T = 2300 (K) [29]. simulation results, one can also find τ α in terms of D as τ α ( T ) = B r /D ( T ) = B r τ L ( T ) , (29)where B r ≃ . for (F) and 0.875 for (S). Hence we obtain τ L /τ α ≃ . for(F) and 1.182 for (S).In Fig. 5(a), the scaled friction coefficient ζ is plotted versus u . Although ζ isconstant, ζ increases as u increases because it depends on S ( q m ) and v th . Herewe note that ζ is an universal function of u in different systems of type i . In fact,the values of ζ obtained from the simulation results [29] for Al O with the Born-Meyer potential [46] and for SiO with the Beest-Kramer-Santen (BKS) potential[47] are also shown in Fig. 5(a), where ζ ≃
82 and q m σ =4.25 for Al O and ζ ≃
113 and q m σ =1.65 for BKS. In each type of liquids ζ is thus shown to coincide witheach other. Hence the dynamics in different systems of type i is expected to coincidewith each other if u has the same value in those systems [29]. As one of suchexamples, in Fig. 5(b) the simulation results for SW at T = 0 . are comparedwith those for Al O at T = 2300 K, where u ≃ . and ζ ≃ . in both systems.Both simulation results are thus shown to coincide with each other within error.Here we note that those are only available data to be compared at the same valueof u in [S]. This is because those simulations have been done independently. Inorder to confirm the universality in dynamics, therefore, the simulations in varioussystems should be done at the same value of u .As discussed in Ref. [30], the simulation results satisfy the relation given by Eq.(16) near the glass transition, where the relevant exponents are listed in Table 1.In Fig. 6, we show how the simulation results for SW and NV are well described14 ig. 6. (Color online) A plot of the scaled function f ( τ ) /f c versus log( τ /τ α ) near theglass transition. (a) The long-dashed and the dashed lines indicate the simulation resultsfor SW at T = f c = 0 . . Here log( τ α ) = 3 . for T = 0 . and 4.621 for 0.556. (b) The long-dashed and the dashed lines indicate thesimulation results for NV at T = f c = 0 . . Here log( τ α ) = 3 . at T = 2900 (K) and 4.101 at 2800 (K). In both figures, the solid linesindicate the scaled KWW decay and the dotted lines the scaled VS decay. by the VS decay and the KWW decay near the glass transition within error. Thoseresults are also used later to check from a unified point of view wether the numericalsolutions for both models can be described by those decays well or not. We now solve Eq. (23) for the new model by iterations numerically. In Figs. 7(a)and 7(b), the iteration results are plotted together with the simulation results nearthe glass transition for SW and NV, respectively. Starting from the simulation re-sults, the numerical solutions are shown to deviate from the simulation results onlyat β stage for . < τ ≤ τ x . As the iteration increases, the plateau height decreasesin SW, while it increases in NV. This difference is clearly seen in κ . As the iterationincreases, κ increases in SW and reaches a constant value around 11.595 within er-ror of order − after 16th iterations, while it decreases in NV, reaching a constantvalue around 18.430 within error of order − after 15th iterations. In general, thenumber of iterations increases as T decreases. At each temperature T (or u ) onecan thus fix the value of µ at a first iteration. In Fig. 8(a), the u dependence of thenonlinear parameter µ is then shown for SW and NV. As u increases, µ is expectedto grow in both systems since the magnitude of the nonlinear fluctuations becomeslarger. In fact, the increment of the nonlinearity δµ (= µ ( u ) − µ ( u s )) is found in[S]. However, the numerical results for µ suggest that the supercooled state [S]should be further separated into two substates, a weakly supercooled state [S w ] and15 ig. 7. (Color online) A log plot of the scattering function f ( τ ) versus τ near the glass tran-sition. (a) The dotted line indicates the simulation results for SW at T = 0 .
556 ( u ≃ . .The solid lines indicate the numerical solutions for the new model obtained by 16th it-erations (from top to bottom) at µ = 1 . . The symbol ( • ) indicates the crossover time log( τ x ) ≃ . . (b) The dotted line indicates the simulation results for NV at T = 3000 (K) ( u ≃ . . The solid lines indicate the numerical solutions for the new model ob-tained by 15th iterations (from bottom to top) at µ = 3 . . The symbol ( • ) indicates thecrossover time log( τ x ) ≃ . . The symbols ( ✷ ) in the insets indicate the numericalvalues of κ obtained by iterations. The solid lines in the insets are guides for eyes.Fig. 8. (Color online) A plot of the relevant parameters for the new model versus u , where(a) µ and (b) κ . In both figures, the symbols ( • ) indicate the fitting values for SW and ( ✷ )for NV. The symbols ( ✸ ) indicate the glass transition point u g , ( ⊙ ) the supercooled point u s , and ( + ) the deeply supercooled point u f , where the numerical values of u i are listed inTable 2. The solid lines are guides for eyes. The label [L] stands for a liquid state, [S w ] fora weakly supercooled state, [S f ] for a deeply supercooled state, and [G] for a glass state. ig. 9. (Color online) A plot of f ( τ ) versus log( τ ) for different temperatures for the newmodel. (a) The dotted lines indicate the simulation results for SW at T = [L] 5.00, 1.25,[S w ] 0.833, 0.714, [S f ] 0.625, and 0.556 from left to right. The symbols ( • ) indicate thecrossover time in [S f ] given by log( τ x ) ≃ . T=0.644 ) , . (0.625), 3.615 (0.556),and 4.931 (0.527), where f ( τ x ) ≃ . . The dashed lines indicate the numerical solutionsat T = T f ) and 0.527 ( T g ). (b) The dotted lines indicate the simulation results forNV at T = [L] 5500, 4300, [S w ] 3600, [S f ] 3300, 3000, 2900, [G] 2800, and 2600 (K) fromleft to right. The symbols ( • ) indicate the crossover time in [S f ] given by log( τ x ) ≃ . ( T = 3300 K), 2.192 (3000), 2.551 (2900), and 2.978 (2800), where f ( τ x ) ≃ . . Thesolid lines in both figures indicate the numerical solutions for the new model. a deeply supercooled state [S f ]. In fact, in each state the value of µ is given for (F)and (S) by µ ≃ . , [L] for u ≤ u s . u + 0 . , [S w ] for u s ≤ u < u f . , [S f ] for u f ≤ u ≤ u g (30) µ ≃ . [L] for u ≤ u s . u − . , [S w ] for u s ≤ u < u f . , [S f ] for u f ≤ u ≤ u g (31)respectively. In Fig. 8(b), the coupling parameter κ is shown to grow rapidly as u increases, except in [S f ] where it grows slowly.In Figs. 9(a) and 9(b), the numerical solutions at different temperatures are com-pared with the simulation results for SW and NV, respectively. For higher temper-atures in [L], the solutions are shown to coincide with the simulation results wellwithin error. For lower temperatures, the solutions are shown to agree with themwell within error, except at β stage. Since one can find a value of µ approximatelyin [S], one can obtain the numerical solutions at any temperatures in [S], eventhough there is no simulation result performed at such temperatures. As typical ex-amples, therefore, we also plot the solutions at T = T f and T g in Fig. 9(a). Thedeviation at β stage is clearly seen in a supercooled state [S] and a glass state [G],17 ig. 10. (Color online) A plot of the scaled function f ( τ ) /f c versus log( τ /τ α ) for thenew model near the glass transition. (a) The long-dashed and the dashed lines indicatethe numerical solutions for SW at T = f c = 0 . and log( τ b /τ α ) = 1 . . (b) The long-dashed and the dashed lines indicate the numericalsolutions for NV at T = f c = 0 . . The detailsare the same as in Fig. 6. while they are very small in a liquid state [L]. As pointed out before, such a de-viation just results from the ideal TMCT equation. Hence we should mention herethat such a technical deviation must disappear if one can solve the original TMCTequation given by Eq. (4) for the new model numerically.Similarly to Fig. 6, we also explore the dynamics near the glass transition carefullyby using the KWW function and the VS decay. The scaled function f ( τ ) /f c is plot-ted versus scaled time τ /τ α in Figs. 10(a) and 10(b) for SW and NV, respectively.Thus, the numerical solutions for the new model in both systems are shown to bedescribed by the KWW function well within error for τ x ≤ τ ≤ τ α . Here we notethat the numerical solutions start to deviate from the KWW function for τ > τ α .This is just because the long time dynamics is governed by the diffusion processgiven by f ( τ ) = e − τ/τ L .Finally, we compare the numerical results for the new model with those for theconventional model. Starting from the same initial conditions as those for the newmodel, one can also obtain the numerical solutions for the conventional model con-sistently. As shown in Fig. 1(b), the u (or T c /T ) dependence of ε is quite differentfrom that of µ . The most important difference appears in [S f ]. For a new model, µ is constant as T decreases, while ε grows monotonically. On the other hand, inother states the u dependence of both nonlinear parameters is qualitatively similarto each other but quantitatively not. As shown in Ref. [29], the nonlinear parameter ε also grows linearly in u for (F) and (S). In fact, in each state the value of ε is given18 ig. 11. (Color online) A plot of f ( τ ) versus log( τ ) for different temperatures. (a) Thedotted lines indicate the numerical solutions for the new model in SW at T = [L] 5.00, 1.25,[S w ] 0.833, 0.714, [S f ] 0.644 ( T f ), 0.625, 0.556, and [G] 0.527 ( T g ) from left to right. Thesymbols ( • ) indicate the crossover time in [S f ] given by log( τ x ) ≃ . ( T = 0 . ),2.306 (0.625), 3.615 (0.556), and 4.931 (0.527), where f ( τ x ) ≃ . . (b) The dotted linesindicate the numerical solutions for the new model in NV at T = [L] 5500, 4300, [S w ] 3600,[S f ] 3300, 3000, 2900, [G] 2800, and 2600 (K) from left to right. The symbols ( • ) indicatethe crossover time in [S f ] given by log( τ x ) ≃ . ( T = 3300 K), 2.192 (3000), 2.551(2900), and 2.978 (2800), where f ( τ x ) ≃ . . The solid lines in both figures indicate thenumerical solutions for the conventional model. for (F) and (S) by ε ≃ . , [L] for u ≤ u s . u + 2 . , [S w ] for u s ≤ u < u f . u + 2 . , [S f ] for u f ≤ u ≤ u g (32) ε ≃ . [L] for u ≤ u s . u − . , [S w ] for u s ≤ u < u f . u + 3 . , [S f ] for u f ≤ u ≤ u g (33)respectively. The coupling parameter κ is also shown to grow monotonically in SWand NV as u increases. The u dependence of κ seems to be rather similar to thatfor the new model. This is just because κ is determined by Eq. (21), where f ( τ ) changes, depending on u , whether µ is constant or not in [S f ].In Figs. 11(a) and 11(b), the numerical solutions for the conventional model arenow compared directly with those for the new model in SW and NV, respectively.For higher temperatures in [L], both solutions are shown to coincide with each otherwell within error. As T decreases, the difference between both solutions becomesclear. At β stage for . < τ ≤ τ x , the plateau height of f ( τ ) for the conventionalmodel is shown to be always lower than that for the new one. At α stage for τ x ≤ τ ≤ τ α , the distinct difference between both solutions becomes obvious. Here we19ote that one can also obtain the numerical solutions for the conventional model atany temperatures in [S] since one can find a value of ε approximately by using alinear equation for u . As typical examples, therefore, we also compare the solutionsfor SW at T = T f and T g in Fig. 11(a). In this paper, we have proposed the new simplified model given by Eq. (25) for thenonlinear memory function in the TMCT equation. This model has contained twounknown parameters µ and ν . In order to fix a value of ν , we have positively usedtwo types of decays, the VS decay and the KWW decay, near the glass transition.Then, we have checked whether the new model can describe the nonlinear memoryfunction obtained numerically by directly solving the ideal TMCT equation givenby Eq. (7) with the PY static structure factor or not. Thus, we have shown that thenew model can describe it very well, while the conventional model can not for anyvalues of ε . We have first transformed the ideal TMCT equation into an universalform given by Eq. (20) because there exists only one solution at each value of u forsame type of liquids. Hence we have chosen A B as a typical example of fragileliquids and SiO as a typical example of strong liquids. In order to fix a value of µ , we have used the simulation results performed on those systems. We have thensolved the universal equation numerically for both models under the same initialconditions. Then, we have compared both numerical solutions together with thesimulation results from a unified point of view consistently. Thus, we have shownthat the new model can describe the simulation results well not only qualitativelybut also quantitatively, except at β stage because of the ideal equation, while theconventional model can describe them approximately only at higher temperaturesin [L].The new model has suggested that the supercooled state should be further separatedinto two substates, a weekly supercooled state [S w ] and a deeply supercooled state[S f ]. In [S f ] the nonlinear parameter µ becomes constant, up to T g . Hence this situ-ation seems to be similar to that in [L] since µ is constant in [L], up to T s . However,there exists an essential difference between them. In [S f ] the magnitude of the equi-librium density fluctuation ρ ( x , t ) around ρ is large and u obeys a power law givenby Eq. (28), where x is a position vector. On the other hand, in [L] it is small and u is approximately described by a linear equation for T − . In [S w ] µ grows rapidlyas T decrease. This situation is similar to that in [G]. However, there also exists anessential difference between them. In [S w ] the system is equilibrium and the mag-nitude of the density fluctuations is larger as T decreases. On the other hand, in [G]the system is out of equilibrium and the magnitude of the nonequilibrium densityfluctuation around the average density ¯ n ( x , t ) is small, where ¯ n ( x , t ) obeys a de-terministic nonlinear equation and describes a nonequilibrium state. Because of thesmall fluctuations, therefore, u increases slowly, following the Arrhenius law given20y u ∼ T − . Those differences may be explored more clearly by using the so-called spatial heterogeneity [48,49,50,51,52,53,54]. As discussed in Refs. [49,50],the physical picture is as follows. In [S w ] the spatially heterogeneous regions with | ρ ( x , t ) | > ρ start to change in space and time at T . As T decreases, such re-gions grow larger, leading to an increment of µ . Among those regions, there arestill enough space for particles to diffuse, leading to medium values of a diffusioncoefficient D ( T ) . In [S f ] the heterogeneous regions seem to almost cover a wholespace and change slowly in space and time at T . As T decreases, such regionsgrow slowly, keeping µ constant. There exist small cages formed by surroundingheterogeneous regions, within which particles can diffuse, leading to smaller val-ues of D ( T ) . In [G] the heterogeneous regions with | ¯ n ( x , t ) | > ρ cover a wholespace and change very slowly in space and time at T . As T decreases, the hetero-geneous regions do not hardly grow, while the nonequilibrium state determined by ¯ n ( x , t = 0) depends on T sensitively, leading to a sharp increment of nonlinearity.There still exist smaller cages for particles to diffuse, following the Arrhenius law.This will be discussed elsewhere.Finally, we have mentioned that the u dependence of µ and κ for (F) is qualitativelysimilar to that for (S) but it is quantitatively different from that. This is reasonablebecause such a difference is originally based on the fact that the static structurefactor of SiO is structurally quite different from that of fragile liquids because theformer has a network structure [55,56,57]. Hence this confirms that although theuniversality in the dynamics of glass-forming liquids holds in the same type, thedynamics of different types does not coincides with each other even if u has thesame value in each type. In this paper, we have chosen two types of the simulationresults as simple examples, SW for (F) and NV for (S). In order to confirm theuniversal properties and also to obtain precise values of universal quantities, oneneeds to perform more extensive molecular-dynamics simulations on various glass-forming materials at the same value of u consistently.AcknowledgmentsThe author wishes to thank T. Narumi for supplying his unpublished numericalsolutions of ideal TMCT equation with the PY static structure factor. This workwas partially supported by Institute of Multidisciplinary Research for AdvancedMaterials, Tohoku University, Japan. References [1] E. Leutheusser, Phys. Rev. A (1984) 2765.[2] U. Bengtzelius, W. G ¨otze, A. Sj¨olander, J. Phys. C (1984) 5915.[3] W. G ¨otze, in: J.P. Hansen, D. Levesque, J. Zinn-Justen (Eds.), Liquids, Freezing andGlass Transition, North-Holland, Amsterdam, 1991.
4] E. W. Fischer, Physica A (1993) 183.[5] Relaxation Kinetics in Supercooled Liquids - Mode Coupling Theory and ItsExperimental Tests, edited by S. Yip, special issue of Trans. Theory Stat. Phys. (6-8) 1995.[6] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, and S. W. Martin, J. Appl.Phys. (2000) 3113.[7] P. G. Debenedetti and F. H. Stillinger, Nature (London) (2001) 259.[8] Proceedings of the 4th International Discussion Meeting on Relaxation in fragileSystems, edited by K. L. Ngai, J. Non-Cryst. Solids (2002) 1.[9] P. G. Debenedetti, J. Phys.: Condens. Matter (2003) R1669.[10] K. Binder and W. Kob, Glassy Materials and Disordered Solids (World Scientific,Singapore, 2005).[11] W. G ¨otze, Complex Dynamics of Glass-Forming Liquids (Oxford University Press,New York, 2009).[12] Proceedings of the International Symposium on Slow Dynamics in Complex Systems,edited by M. Tokuyama and I. Oppenheim (AIP, New York, 2004, 2013).[13] L. Berthier and G. Biroli, Rev. Mod. Phys. (2011) 587.[14] M. Tokuyama, Physica A (2014) 31.[15] M. Micoulaut, Rep. Prog. Phys. (2016) 066504.[16] M. Tokuyama, Physica A (2015) 156.[17] E. R. von Schweidler, Annalen der Physik (1907) 711.[18] R. Kohlrausch, Ann. Phys. (1847) 393.[19] G. Williams and D. C. Watts, Trans. Faraday Soc. (1980) 80.[20] Y. Kimura and M. Tokuyama, IL NUOVO CIMENTO C (2016) 300.[21] T. Narumi and M. Tokuyama, Phys. Rev. E (2017) 032601.[22] M. Tokuyama and T. Narumi, Physica A (2019) 533.[23] J. K. Percus and G. J. Yevick, Phys. Rev. (1958) 1.[24] M. Tokuyama, H. Yamazaki, Y. Terada, Phys. Rev. E (2003) 062403.[25] M. Tokuyama, H. Yamazaki, Y. Terada, Physica A (2003) 367.[26] M. Tokuyama, Y. Terada, J. Chem. Phys. B (2005) 21357.[27] M. Tokuyama, T. Narumi, E. Kohira, Physica A 385 (2007) 439.[28] M. Tokuyama, Physica A (2017) 453.
29] M. Tokuyama, Physica A (2017) 229.[30] M. Tokuyama, Physica A (2019) 121074.[31] R. Kubo, J. Phys. Soc. Japan (1962) 1100.[32] H. Mori, Prog. Theoret. Phys. (1965) 423.[33] M. Tokuyama, H. Mori, Progr. Theoret. Phys. (1975) 918.[34] M. Tokuyama, H. Mori, Progr. Theoret. Phys. (1976) 411.[35] M. Tokuyama, Physica A (2007) 157.[36] M. Tokuyama, Phys. Rev. E (2009) 031503.[37] M. Tokuyama, Phys. Rev. E (2010) 041501.[38] M. Tokuyama, J. Phys. Chem. B (2011) 14030.[39] M. Tokuyama, AIP CP (2013) 47.[40] M. Tokuyama, S. Enda, Physica A (2014) 31.[41] M. Tokuyama, S. Enda, and J. Kawamura, Physica A (2016) 1.[42] S. Enda, M. Tokuyama, AIP Conf. Proc. (2013) 324.[43] T. Narumi, unpublished.[44] T. A. Weber and F.H. Stillinger, Phys. Rev. B (1985) 1954.[45] A. Nakano, L. S. Bi, R. K. Kalia, P. Vashishta, Phys. Rev. B (1994) 9441.[46] V. V. Hoang, Phys. Rev. B bf70 (2004) 134204.[47] B. W. H. van Beest, G. J. Kramer, R. A. van Santen, Phys. Rev. Lett. (1990) 1955.[48] G. Adam and J. H. Gibbs, J. Chem. Phys. (1965) 139.[49] M. Tokuyama, Y. Enomoto, and I. Oppenheim, Phys. Rev. E (1997) 2302.[50] M. Tokuyama, Y. Enomoto, and I. Oppenheim, Physica A (1999) 380.[51] H. Sillescu, Journal of Non-Crystalline Solids (1999) 81.[52] M. D. Ediger, Annu. Rev. Phys. Chem. (2000) 99.[53] E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, Science (2000) 627.[54] R. Richert, J. Phys.: Condens. Matter (2002) R703.[55] R. Zallen, The Physics of Amorphous Materials (New York: Wiley, 1983)[56] U. Mohanty, Adv. Chem. Phys. (1995) 89.[57] W. Kob, J. Phys.: Condens. Matter (1999) R85.(1999) R85.