Slow dynamics in glass-forming materials at α - and β -relaxation stages based on time-convolutionless mode-coupling theory

Physica A 526 (2019) 121074

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Slow dynamics in glass-forming materials at α - and β -relaxation stages based on time-convolutionless mode-coupling theory Michio Tokuyama Institute of Multidisciplinary Research for Advanced Materials, Tohoku University, Sendai 980-8577, Japan

highlights • alpha- and beta-relaxation processes are analyzed based on time-convolutionless mode-coupling theory. • A simple relation between time exponents in alpha and beta stages is found analytically. • Calculated values of time exponents are checked by simulations performed on various glass-forming materials.

article

info

Article history: Received 19 January 2019 Available online 11 April 2019 Keywords: α -relaxation process β -relaxation process Glass transition Supercooled liquids Time-convolutionless mode-coupling theory

a b s t r a c t The slow dynamics of glass-forming materials in α and β stages are analyzed based on the time-convolutionless mode-coupling theory recently proposed by the present author. Similarly to the analyses used in the mode-coupling theory, the two-step relaxation processes, so called critical decay with a time exponent a and von Schweidler decay with a time exponent b, are first discussed in β stage. In α stage, the Kohlrausch–Williams– Watts type function, often called stretched exponential with a time exponent β , is also investigated in the manner similar to that in β stage. The exponents a, b, and β are then shown to satisfy the relations

Γ [1 + b]2 Γ [2β − 1] 1 1 Γ [1 − a]2 = = = λ, + = γ, Γ [1 − 2a] Γ [1 + 2b] Γ [2β + 1]Γ [β + 1]Γ [−β − 1] 2a 2b where λ is an exponent parameter, γ a power-law exponent of the α -relaxation time, and Γ [x] an usual Γ -function. Thus, the numerical values of those exponents calculated under a given value of γ are compared with the simulation results for three types of glass-forming liquids with three different values of γ and are shown to describe them well within error. © 2019 Elsevier B.V. All rights reserved.

1. Introduction For past decades, a considerable study has been made to understand the slow dynamics of supercooled liquids and the mechanism of the liquid-glass transition. Many experiments and simulations performed on various glass-forming materials have thus confirmed that there exist two types of long-known relaxation phenomena near the glass transition, such as a power-law decay of von Schweidler type in a β -relaxation stage and a stretched exponential decay in an α -relaxation stage. The two-step relaxation processes in β stage have been revealed by the mode-coupling theory E-mail address: [email protected]. https://doi.org/10.1016/j.physa.2019.121074 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

M. Tokuyama / Physica A 526 (2019) 121074

(MCT) [1–3] and the power-laws decays in each time stage have been successfully found. On the other hand, however, our understanding is still incomplete for slow relaxation processes in α stage, although several empirical laws [4–7] and some theoretical approaches [8–10] have been proposed. In the present paper, therefore, we employ the time-convolutionless mode-coupling theory (TMCT) [11–14] recently proposed to study such long-known relaxation phenomena. Similarly to MCT, TMCT has been also shown to cause an ergodic to a non-ergodic transition at a critical point χ = χc , below which the intermediate scattering function reduces to a non-zero value, the so-called nonergodicity parameter, where χ is a control parameter, such as an inverse temperature 1/T and a volume fraction φ . Very recently, the TMCT equation has been solved numerically for monodisperse hard spheres [15–17] with the Percus–Yevick static structure factor S(q) [18]. Thus, the critical volume fraction φc has been shown to coincide with that obtained by the simulations. It has been also pointed out in Ref. [17] that the TMCT equation must be numerically solved consistently under given values of S(q, χ ) and the characteristic length ℓ(χ ) at each value of χ to obtain reasonable solutions consistent with the simulations. In general, however, it is difficult to do it since both physical quantities S(q) and ℓ are not known, except for monodisperse hard spheres. In the present paper, therefore, we focus only on the relaxation dynamics in each time stage near χc and then compare the theoretical results with the simulations. Depending on a time scale, there exist two types of relaxation stages in glass-forming materials near the glass transition, a β -relaxation stage and an α -relaxation stage. In β stage, two kinds of relaxation dynamics have been revealed by MCT, the so-called critical decay with a time exponent a and the von Schweidler decay with a time exponent b. In α stage, the Kohlrausch–Williams–Watts type function [4,5], often called stretched exponential with a time exponent β , is known to describe the dynamics near the glass transition well. By employing the formulation similar to that discussed in MCT, we analyze those different decays based on TMCT. The exponents a and b in β stage and the exponent β in α stage are then shown to be related to each other by an equation

λ=

Γ [1 + b]2 Γ [2β − 1] Γ [1 − a]2 = = , Γ [1 − 2a] Γ [1 + 2b] Γ [2β + 1]Γ [β + 1]Γ [−β − 1]

(1)

where λ is the so-called exponent parameter and is the only material-dependent quantity, and Γ [x] an usual Γ -function. This equation has the same form as that obtained by MCT, except the last term. If λ is given numerically, one can then find the values of a, b, and β . However, this is not the case because S(q) and ℓ are not known to calculate λ. Similarly to MCT, the power-law exponent γ of the α -relaxation time also satisfies the relation γ = (a + b)/(2ab). In the following, therefore, we use the value of γ to find a and b instead of λ. This can be done from a mean-field point of view by using the master curve for the diffusion coefficient [19–23], onto which not only the simulation results but also the experimental data are shown to be collapsed for 0 ≤ χ ≤ χg , where χg is a glass transition point and satisfies χg < χc . Once γ is known, one can first obtain a, b, and λ from Eq. (1) and then find β by using λ in Eq. (1). Those exponents are thus checked by the simulations performed on various glass-forming materials and are shown to describe them well within error. 2. Basic equations In this section, we first review the time-convolutionless mode-coupling theory (TMCT) equations briefly and then derive a starting equation to discuss the relaxation dynamics in α and β stages. 2.1. TMCT equation We consider the three-dimensional equilibrium glass-forming system, which consists of N particles with mass m and diameter σ in the total volume V at temperature T . The scaled intermediate scattering function fα (q, t) is given by fα (q, t) = ⟨ρα (q, t)ρα (q, 0)∗ ⟩/Sα (q)

(2)

with the collective density fluctuation

ρc (q, t) =

1 N 1/2

⎡ ⎤ N ∑ ⎣ ρs (q, t) − N δq,0 ⎦ ,

(3)

j=1

and the self-density fluctuation

ρs (q, t) = eiq·X j (t) ,

(4)

where ⟨ρs (q, t)⟩ = δq,0 , ⟨ρc (q, t)⟩ = 0, q = |q|, Sc (q) = S(q), and Ss (q) = 1, – S(q)(= ⟨|ρc (q, 0)| ⟩) – being the static structure factor. Here X j (t) indicates the position vector of a jth particle at time t and the brackets the average over an equilibrium ensemble. Since the density fluctuations ρα (q, t) are macroscopic physical quantities, one must set the magnitude of q as 0 < q ≤ qc . As discussed in the previous paper [17], the cutoff qc should be then chosen to be an inverse characteristic length ℓ(χ )−1 in a supercooled state, where χ is a control parameter, such as in inverse temperature 1/T and a volume fraction φ (= πσ 3 ρ/6), ρ (= N /V ) being the number density. It is convenient to introduce the cumulant function Kα (q, t) by 2

fα (q, t) = exp[−Kα (q, t)],

(5)

M. Tokuyama / Physica A 526 (2019) 121074

3

where Kα (q, 0) = 0 since fα (q, 0) = 1. As shown in the previous papers [12–14], one can then find the second-order differential equation for Kα (q, t) 2 q2 vth ∂ Kα (q, t) ∂ 2 Kα (q, t) = − ζα − ∂t2 Sα (q) ∂t

t

∫

∆ϕα (q, t − τ ) 0

∂ Kα (q, τ ) dτ ∂τ

(6)

with the nonlinear memory function

∆ϕα (q, t) =

2 ρvth

∫

2nα

dk <

(2π )3

vα (q, k)2 S(k)Sα (|q − k |)fc (k, t)fα (|q − k |, t), ∫

(7)

where ζα is a positive constant and < the sum over wave vectors k whose magnitudes are smaller than a cutoff qc (χ ) [17]. Here nc = 1, ns = 0, and vth = (kB T /m)1/2 . The vertex amplitude vα (q, k) is given by

vα (q, k) = qˆ · kc(k) + nα qˆ · (q − k)c(|q − k |),

(8)

where ρ c(k) = 1 − 1/S(k) and qˆ = q/q. Eq. (6) is an ideal TMCT equation, where the initial conditions are given by Kα (q, t = 0) = dKα (q, t)/dt |t =0 = 0. Here we note that since only a long-time relaxation process, such as α - and β relaxation processes, is discussed here, one can rather start from a simple equation given by Eq. (6), instead of a more precise one discussed in Ref. [14]. It is convenient to introduce the Laplace transform fα [q, z ] of fα (q, t) by fα [q, z ] = L[fα (q, t)][z ] =

∞

∫

e−zt fα (q, t)dt .

(9)

0

By taking a Laplace transform of Eq. (6), one then obtains zKα [q, z ] =

2 q2 vth

1

Sα (q) z 2 + ζα z + z ∆ϕα [q, z ]

.

(10)

For a long time (or z → 0) limit, from Eq. (10), one can easily find Kα (q, t) ≃ q2 Dα (q)t

(11)

with the q-dependent long-time diffusion coefficient 2 vth /Sα (q) ∫∞ . ζα + 0 ∆ϕα (q, s)ds

Dα (q) =

(12)

In the present paper, we only discuss two types of diffusion coefficients, Dα (qm ) and Ds (q = 0). 2.2. Ergodic to non-ergodic transition We first discuss an ergodic to non-ergodic transition at a critical point χ = χc , below which the long-time solution fα (q, t) (or Kα (q, t)) reduces to a non-zero value fα (q) (or Kα (q)), the so-called nonergodicity parameter. In the long-time t → ∞ (or z → 0) limit, Eq. (10) thus leads to [12–14] Kα (q) =

1

(13)

Fα (q)

with the long-time limit of the nonlinear memory function Fα (q) = lim

z →0

z ∆ϕα [q, z ]Sα (q) 2 q2 vth

=

1 2nα

∫

dk <

(2π )3

Vα (q, k)fc (k)fα (|q − k |),

(14)

where the vertex Vα is given by Vα (q, k) = (ρ/q2 )vα (q, k)2 Sα (q)S(k)Sα (|q − k |).

(15)

The existence of the critical point in Eq. (13) is mathematically confirmed since Kα (q) is a kind of the Lambert Wfunction. In fact, the critical volume fraction φc has been obtained by solving numerically Eq. (13) with the Percus–Yevick static structure factor as φc ≃ 0.582 [15–17] which agrees with that of simulations performed on monodisperse hard spheres [24–27] within error. 2.3. Starting equation near a critical point Once S(q, χ ) and ℓ(χ ) are given, one can then solve Eq. (6) (or Eq. (10)) numerically. As shown in Refs. [15,17], this is only possible for hard-sphere systems. In the present paper, therefore, we only discuss the asymptotic behavior of Kα (q, t) near the critical point χc in each time stage. From Eq. (10), the long-time dynamics is then determined by 1 zKα [q, z ]

= z L[Fα (q, fc (t), fα (t))][z ].

This is a starting equation to discuss the relaxation dynamics in α and β stages.

(16)

4

M. Tokuyama / Physica A 526 (2019) 121074

Fig. 1. A schematic plot of the scaled cumulant function Mi (t)/σ 2 versus time t /t0 , where t0 = σ /vth . The solid line indicates Mc (t), the dashed line Ms (t), and the dotted line M2 (t). Here tL = σ 2 /Ds (q = 0) and tLα = σ 2 /Dα (qm ).

3. Analyses of relaxation dynamics There exist multi-time stages near the critical point, where the different relaxation dynamics exist. In Fig. 1, we show in a schematic way typical time scales for two types of the scaled cumulant functions, M2 (t) = limq→0 (6/q2 )Ks (q, t) and Mα (t) = (6/q2m )Kα (qm , t), where qm is a first-peak position of S(q). The usual mean-square displacement M2 (t) shows two stages, the β stage for tγ ≤ t ≤ tβ and the late stage for tβ ≪ t, where tβ is a β -relaxation time discussed in the previous papers [19–23] and tγ = 1/γ . In β stage, each particle behaves as if it is trapped in a cage which is mostly formed by neighboring particles. The β stage is further split into two-time stages, a fast β -relaxation stage for tγ ≤ t ≤ tδ and a slow β -relaxation stage for tδ ≤ t ≤ tβ , where tδ is a crossover time to be determined. On a time scale of order tβ , the particles can escape their cages and on a time scale of order tL , they finally reach a long-time self-diffusion process described by M2 (t) ≃ 6Ds (q = 0)t. On the other hand, Mα (t) shows three stages, the β stage for tγ ≤ t ≤ tx , the α stage for tx ≤ t ≪ tLα , and the late stage for tα ≪ t, where tx is a crossover time from the β stage to the α stage and tα an α -relaxation time. Similarly to M2 (t), the β stage consists of a fast β -relaxation stage for tγ ≤ t ≤ tδ and a slow β -relaxation stage for tδ ≤ t ≤ tx . On a time scale of order tLα , the particles are then described by Mα (t) ≃ 6Dα (qm )t. Here we note that tx is physically different from tβ since the former is determined as a crossover time from the von Schweidler decay to the stretched exponential decay, while the latter is a crossover time from the von Schweidler decay to an exponential decay. In the following, we discuss the asymptotic behavior of Kα (q, t) in each stage. 3.1. A two-step relaxation in β stage We first discuss the power-law decays in β stage for tγ ≤ t ≤ tx . In β stage, one can use the formulation similar to that employed in MCT [3]. One can formally split Kα (q, t) into a trivial asymptotic part and a non-trivial part G; Kα (q, t) = Kαc (q)[1 − G(t)],

(17)

where Kαc (q) = Kα (q, χc ). There exist two-time stages; a fast β -relaxation stage for tγ ≤ t ≤ tδ and a slow β -relaxation stage for tδ ≤ t ≤ tx . At a fast β stage G(t) is described by the so-called critical decay G(t) = (ta /t)a ,

(18)

where a is a time exponent to be determined and ta a characteristic time to be determined. At a slow β stage G(t) is described by the so-called von Schweidler decay G(t) = −(t /tb′ )b ,

(19)

M. Tokuyama / Physica A 526 (2019) 121074

5

where b is a time exponent to be determined and tb′ a characteristic time to be determined. Since Kαc (q) ≪ 1 near χc , one can find, on a time scale of order tb′ , fc (k, t)fα (p, t) = fcc (k)fαc (p) exp {Kcc (k) + Kαc (p)}G(t)

[

]

≃ fcc (k)fαc (p) 1 + {Kcc (k) + Kαc (p)}G(t) ] 1 + {Kcc (k) + Kαc (p)}2 G(t)2 + . . . ,

[

(20)

2

c

where fαc (k) = e−Kα (k) and p = q − k. Taking the Laplace transform of Eq. (17) then leads to zKα [q, z ] = Kαc (q)[1 − zG[z ]].

(21)

Expanding the left hand side of Eq. (16) in powers of zG[z ] and using Eq. (20) in the right hand side of Eq. (16), one can obtain in a long time (or z → 0) limit

{zG[z ]}2 = λα z L[G(t)2 ][z ],

(22)

where the exponent parameter λα is given by

λα =

ℓ3

∑

2nα +1 V 2

Vα (q, k)fcc (k)fαc (p)Kαc (q) Kcc (k) + Kαc (p)

[

]2

.

(23)

|q|≤ℓ−1 ,|k |≤ℓ−1

Since

{ zG[z ] =

Γ [1 − a](zta )a , for a fast β stage, −Γ [1 + b](ztb′ )−b , for a slow β stage,

(24)

and

{

2

z L[G(t) ][z ] =

Γ [1 − 2a](zta )2a , for a fast β stage, Γ [1 + 2b](ztb′ )−2b , for a slow β stage,

(25)

use of Eq. (22) then leads to

Γ [1 − a]2 Γ [1 + b]2 = = λα . Γ [1 − 2a] Γ [1 + 2b]

(26)

This relation is the same as that obtained by MCT [3], except that the detailed formulation of λα is different from that of MCT. Here λα is the only material-dependent quantity which can be calculated from the vertices Vα . As is discussed later, however, it is just treated as an adjustable fitting parameter here. 3.2. A stretched exponential in α stage We now discuss the so-called Kohlrausch–Williams–Watts function (KWW) [4,5] in α stage, often also called stretched exponential fα (q, t) = e−Kα (q,t) with Kα (q, t) = Kαc (q) + (t /tα )β ,

(27)

where β is the KWW exponent and tα the so-called α -relaxation time. We should first mention that in order to find this stretched behavior, one has to solve the asymptotic Eq. (16) consistently for a whole-time region of α stage. This can be done only numerically. Hence we here discuss only the asymptotic solution of Eq. (16) to find the exponent β . In order to use the exponent parameter λα even in α stage, it is convenient to write Eq. (27) as Kα (q, t) = Kαc (q)[1 + (t /tα′ )β ],

(28)

1/β

c

tα′

tα′

where = Kα (q) tα (≪ tα ). Here we first show that the time is identical to the crossover time tx from the von Schweidler decay to the stretched exponential. In order to check this, we first rewrite Eq. (17) with Eq. (19) as Kα (q, t) = Kαc (q)[1 + (t /tb′ )b ] = Kαc (q) + (t /tb )b , where tb = tb /Kα (q) c

′

1/b

′

(29) ′

(≫ tb ). On a time scale of order tb , therefore, one can write the von Schweidler decay as

fα (q, t) ≃ fαc (q)[1 − (t /tb )b ].

(30)

By using Eqs. (28) and (30), one can then find the crossover time tx through the relation exp[−Kαc (q)(tx /tα′ )β ] = 1 − (tx /tb )b . When tx = tx =

tα′ , tα′

(31)

from Eq. (31) we obtain

= tb [1 − fαc (q)]1/b ,

(32)

6

M. Tokuyama / Physica A 526 (2019) 121074

Fig. 2. A plot of the collective-intermediate scattering function fc (q, t) versus time t for (a) a fragile liquid and (b) a strong liquid at q = qm . (a) The solid line indicates the simulation results for the binary mixtures A80 B20 with the Stillinger–Weber potential from Ref. [13] at T = 0.625, where qm σ = 7.25. The long-dashed line the von Schweidler decay with b = 0.3064 and tb /t0 = 104.278 and the dashed line the exponential decay with β = 0.6832 and tα /t0 = 102.681 , where fcc (qm ) = 0.86 and fc (qm , tα′ ) ≃ 0.7396. The symbol (□) indicates tx /t0 = tα′ /t0 = 101.478 , (•) tα /t0 , (△) tδ /t0 = 10−0.08 , and (▽) tγ /t0 = 10−1.08 . (b) The solid line indicates the simulation results for SiO2 with the Nakano–Vashishta potential from Ref. [13] at T = 2700 (K), where qm = 1.55 (Å−1 ). The long-dashed line the von Schweidler decay with b = 0.2779 and tb = 105.739 (ps) and the dashed line the exponential decay with β = 0.6809 and tα = 103.281 (ps), where fcc (qm ) = 0.922 and fc (qm , tα′ ) ≃ 0.850. The symbol (□) indicates tx = tα′ = 101.679 (ps), (•) tα , (△) tδ = 10−1.0 (ps), and (▽) tγ = 10−1.982 (ps).

and c

fα (q, tα′ ) = fαc (q)2 = e−2Kα (q) .

(33)

Since 1 − fαc (q) ≃ Kαc (q), one can thus find from Eq. (32) that tx = tα′ ≃ tb′ . Thus, it turns out that tα′ and tb′ are identical to the crossover time tx . These results are numerically checked later together with two different types of relaxation dynamics, von Schweidler decay and KWW decay, by employing the simulation results performed on different glass-forming materials (see Fig. 2 as simple examples).

M. Tokuyama / Physica A 526 (2019) 121074

7

Fig. 3. A plot of the scaled collective-intermediate scattering function fc (q, t)/fcc (q) versus scaled time tˆ for different temperatures at q = qm , where qm σ = 7.25. (a) α -relaxation stage with tˆ = t /tα , where the solid line indicates the stretched exponential given by y = exp[−tˆ β ]. (b) β -relaxation stage with tˆ = t /tb , where the solid line indicates the von Schweidler type decay given by y = 1 − tˆ b . The time exponents β and b are listed in Table 1. The dotted line indicates the simulation results for the binary mixtures A80 B20 with the SW potential from Ref. [23] at T = 0.5, the dot-dashed line at T = 0.556, the dashed line at T = 0.625, and the long-dashed line at T = 0.714.

We now discuss how to find the exponent β . The Laplace transform of Eq. (28) leads to zKα [q, z ] = Kαc (q)[1 + (ztα′ )−β Γ [1 + β]].

(34)

In the long-time (or z → 0) limit, the l.h.s. of Eq. (16) can be then written as 1 zKα [q, z ]

≃

(ztα′ )β

Kαc (q)Γ [β + 1]

.

(35)

Since Kαc (q) ≪ 1 near χc , one can also obtain, on a time scale of order tα′ , fc (k, t)fα (p, t) ≃ fcc (k)fαc (p) 1 − {Kcc (k) + Kαc (p)}(t /tα′ )β

[

+

1 2

] {Kcc (k) + Kαc (p)}2 (t /tα′ )2β + . . . .

(36)

8

M. Tokuyama / Physica A 526 (2019) 121074

Fig. 4. A plot of the scaled self-intermediate scattering function fs (q, t)/fsc (q) versus scaled time tˆ for the binary mixtures A80 B20 . The details are the same as in Fig. 3.

Using Eq. (36), one can thus write the r.h.s. of Eq. (16) as z L[Fα (q, fc (t), fα (t))][z ] ≃ λ′α (q)Γ [2β + 1](ztα′ )−2β ,

(37)

′ c q λα (q)Kc (q)

where = λα . In Eqs. (35) and (37), only the essential parts for long times are left. Hence the inverse Laplace transform of Eq. (35) does not necessarily coincide with that of Eq. (37) at t = tα′ . However, the tangent lines of them (or time derivatives of them) at t = tα′ should coincide with each other. Multiplying Eqs. (35) and (37) by z and taking the inverse Laplace transform of them, we then find

∑

(t /tα′ )−β t −2

Kαc (q)Γ [β + 1]Γ [−β − 1]

= λ′α (q)

Γ [2β + 1] (t /tα′ )2β t −2 . Γ [2β − 1]

(38)

At t = tα′ we thus obtain

Γ [2β − 1] = λα . Γ [β + 1]Γ [2β + 1]Γ [−β − 1]

(39)

M. Tokuyama / Physica A 526 (2019) 121074

9

Fig. 5. A plot of the scaled self-intermediate scattering function fs (q, t)/fsc (q) versus scaled time tˆ for different temperatures at q = qm , where qm = 4.25 Å−1 . The dotted line indicates the simulation results for Al2 O3 with the BM potential from Ref. [23] at T = 2200 K, the dot-dashed line at T = 2300 K, the dashed line at T = 2600 K, and the long-dashed line at T = 2800 K. The details are the same as in Fig. 3.

This is finally combined with Eq. (26) to find

Γ [1 + b]2 Γ [2β − 1] Γ [1 − a]2 = = = λα . Γ [1 − 2a] Γ [1 + 2b] Γ [2β + 1]Γ [β + 1]Γ [−β − 1]

(40)

This is an equation to obtain the exponent β together with the exponents a and b under a given value of λα . Next, we discuss how to find a numerical value of λα . 3.3. Time exponents We discuss how to find time exponents consistently. If the value of λα is known, one can easily calculate the time exponents, a, b, and β by using Eqs. (26) and (39). In general, however, it is difficult to do it because λα contains the static structure factors, whose analytic functions are not known, except the Percus–Yevick static structure factor. In the following, therefore, we take a different approach to find those exponents.

10

M. Tokuyama / Physica A 526 (2019) 121074

Fig. 6. A plot of the scaled collective-intermediate scattering function fc (q, t)/fcc (q) versus scaled time tˆ for different temperatures at q = qm , where qm = 1.55 Å−1 . The dotted line indicates the simulation results for SiO2 with the NV potential from Ref. [23] at T = 2600 K, the dot-dashed line at T = 2700 K, the dashed line at T = 2800 K, and the long-dashed line at T = 2900 K. The details are the same as in Fig. 3.

The formulation employed in MCT [3] can be also applied to TMCT. Then, one can find 1 −γ tα ∝ tb ∝ D− , α (qm ) ∝ |1 − χ/χc |

(41)

with the exponent

γ = (a + b)/(2ab),

(42)

where Dα (qm ) is a long-time diffusion coefficient given by Eq. (12). Once the value of γ is known, therefore, one can obtain the exponents a and b numerically from Eq. (40). Since γ has the same value both for a self case and for a collective case, the exponents a and b also have the same values for both cases. Thus, it turns out from Eq. (40) that the values of β and λα do not depend on cases, either. In the following, we show how γ is obtained. As discussed in Refs. [21–23], the mean-field analyses of the simulation results and the experimental data for various glass-forming materials suggest that the α - and β -relaxation times tα and tβ of the self-diffusion process of a tagged particle in a supercooled state obey the power laws

˜ s (q = 0)−(1+µ) , tα ∝ D

˜ s (q = 0)−(1−µ) tβ ∝ D

(43)

M. Tokuyama / Physica A 526 (2019) 121074

11

Fig. 7. A plot of the scaled self-intermediate scattering function fs (q, t)/fsc (q) versus scaled time tˆ for different volume fractions at q = qm , where qm σ = 7.0. The dotted line indicates the simulation results for the hard-sphere fluid with 10% size polydispersity from Ref. [28] at φ = 0.59, the dot-dashed line at φ = 0.58, and the long-dashed line at φ = 0.57. The details are the same as in Fig. 3.

˜ s (= qm Ds (q = 0)/vth ), where the value of the exponent µ depends with the scaled long-time self-diffusion coefficient D on types of glass-forming liquids; µ ≃ 2/10 for fragile liquids and 2/11 for strong liquids. Thus, it turns out that tα is physically different from tβ . Here Ds (q = 0) has the same critical point as that of Dα (qm ) and obeys a power law [12,13] Ds (q = 0) ∝ |1 − χ/χc |−γs ,

(44)

where γs is an exponent to be determined. Here the value of γs depends on the control parameter χ and also on types of liquids [21,22]. If χ is an extensive parameter such as volume fraction φ , one has γs = 2.0. On the other hand, if χ is an intensive parameter such as inverse temperature 1/T , one has γs ≃ 10/3 for fragile liquids and 11/3 for strong liquids. Use of Eqs. (41), (43) and (44) then leads to γ = (1 + µ)γs whose value is listed in Table 1. By using the value of γ , one can now obtain the numerical values of a, b, β , and λ from Eq. (40), which are also listed in Table 1. 4. Comparison with simulations In order to verify the time exponents listed in Table 1, we use the simulation results performed on different glassforming liquids [23,28,29]. For fragile liquids with the control parameter 1/T , we take the binary mixtures A80 B20 with

12

M. Tokuyama / Physica A 526 (2019) 121074

Fig. 8. A plot of the scaled self-intermediate scattering function fs (q, t)/fsc (q) versus scaled time tˆ for different volume fractions at q = qm , where qm σ = 7.0. The dotted line indicates the simulation results for the soft-sphere fluid with 15% size polydispersity from Ref. [29] at φ = 0.60, the dot-dashed line at φ = 0.59, the dashed line at φ = 0.58, and the long-dashed line at φ = 0.57. The details are the same as in Fig. 3. Table 1 Time exponents a, b, and β , and exponent parameter λ for different systems. System

γ

γs

a

b

β

λα

Control parameter 1/T Fragile liquids Strong liquids

12/3 13/3

10/3 11/3

0.2111 0.1973

0.3064 0.2779

0.6832 0.6809

0.8980 0.9130

Control parameter φ Hard spheres Soft spheres

2.4 2.4

2.0 2.0

0.3177 0.3177

0.6051 0.6051

0.710 0.710

0.7215 0.7215

the Stillinger–Weber (SW) potential [30] and Al2 O3 with the Born–Meyer (BM) potential [31]. For strong liquids with the control parameter 1/T , we take SiO2 with the Nakano–Vashishta (NV) potential [32] as a typical example of network glass formers. For fragile liquids with the control parameter φ , we take the hard-sphere fluids at 10% size polydispersity from Ref. [28] and the soft-sphere fluids with the potential U(r) ∼ r −36 at 15% size polydispersity from Ref. [29].

M. Tokuyama / Physica A 526 (2019) 121074

13

We first check two types of decays, the stretched exponential decay given by Eq. (27) and the von Schweidler decay given by Eq. (30), together with the crossover time tx = tα′ given by Eq. (32) and the height fα (qm , tα′ ) given by Eq. (33). As a simple example of fragile liquids, we take the simulation results [12,13,23] for the collective-intermediate scattering function of the Stillinger–Weber binary mixtures A80 B20 , while as an example of strong liquids, we take the simulation results [13,23] for the collective-intermediate scattering function of SiO2 with the NV potential. As is shown in Fig. 2, it turns out that the both simulation results are well described by two types of decays in each stage within error, whose crossover time satisfies Eqs. (32) and (33). These results also hold for other different glass-forming materials discussed above. In the following, therefore, we only discuss the universal behavior for the scaled intermediate scattering function fα (qm , t)/fαc (qm ) versus scale time tˆ , where (a) tˆ = t /tα and (b) tˆ = t /tb . We simply skip the critical decay given by Eq. (18) for simplicity. We next discuss a case where the control parameter χ is given by 1/T . We first investigate the fragile liquids. In Fig. 3, the simulation results on the SW binary mixtures A80 B20 for the scaled collective-intermediate scattering function fc (q, t)/fcc (q) are plotted versus scaled time tˆ for different temperatures at qm σ = 7.25, where (a) tˆ = t /tα and (b) tˆ = t /tb . They are shown to be well described by two types of decays near the critical temperature Tc ≃ 0.5 within error. In order to check whether the same values of a, b, and β as those obtained for a collective case holds even for a self case or not, the scaled self-intermediate scattering function fs (q, t)/fsc (q) are also plotted versus scaled time tˆ in Figs. 4 and 5. In Fig. 4, the simulation results on the SW binary mixtures A80 B20 are shown for different temperatures. They are shown to be well described by two types of decays near the critical temperature Tc ≃ 0.5 within error. In Fig. 5, the simulation results on Al2 O3 are shown for different temperatures at qm = 4.25 Å−1 . They are well described by two types of decays near Tc ≃ 1947.8 K. Here we note that the fluctuations in the simulation results on Al2 O3 are rather larger than those in the SW mixtures. This is mainly due to the fact that since the potential of the former is a long-range interaction, one needs more times to perform the simulations, while that of the latter is a short-range interaction. We next investigate the strong liquids. In Fig. 6, the simulation results on SiO2 with NV potential for fc (q, t)/fcc (q) are shown for different temperatures at qm = 1.55 Å−1 . They are well described by two types of decays near Tc ≃ 2544.2 K within error. We should mention here that the simulation results on the BKS potential [33] also satisfy the same decays as those for the NV potentials. We finally discuss a self case where χ is given by φ . In Fig. 7, the simulation results on the hard-sphere fluids with 10% size polydispersity [28] are shown for different volume fraction at qm σ = 7.0. They are well described by two types of decays near φc ≃ 0.5842. In Fig. 8, the simulation results on the soft-sphere fluids with 15% size polydispersity [29] are also shown for different volume fraction at qm σ = 7.0. They are also well described by two types of decays near φc ≃ 0.5952. 5. Summary In this paper, we have investigated two types of relaxation dynamics in α and β stages near the critical point based on TMCT. Employing the approach similar to that discussed in MCT and using the asymptotic TMCT equation given by Eq. (16), we have first shown how to find the time exponents a and b in the β -relaxation stage. Starting from Eq. (16), we have also explored how to obtain the time exponent β consistently in the α -relaxation stage. Thus, we have derived Eq. (26) to find a and b and also derived Eq. (39) to find β . Hence the exponent parameter λα turns out to just play a role in combining Eq. (26) with Eq. (39), leading to Eq. (40). In general, however, it is difficult to calculate it because the static structure factor is unknown. Instead of calculating λα , therefore, we have first used Eqs. (40) and (42) to find a and b, where the exponent γ are assumed to be given by the mean-field analyses. Once the value of γ is given, λα is also calculated through Eq. (40) together with a and b. Then, β is finally calculated through Eq. (40) under a given value of λα . All the values are thus listed in Table 1, where the value of γ depends not only on types of glass-forming materials, fragile liquids and strong liquids, but also on kinds of control parameters, an intensive parameter and an extensive parameter. In order to investigate the validity of those numerical values, we have then used the simulation results performed on various glass-forming materials. Thus, we have confirmed that all the simulation results are well described by the von Schweidler decay and the KWW function with the calculated values of a, b, and β within error. Finally, we should point out here that if TMCT equation is numerically solved under given values of S(q) and ℓ at each value of χ consistently, the value of γ must coincide with the mean-field value within error. In fact, this has been shown to be true for monodisperse hard-sphere fluids [15–17]. Even for other systems, therefore, this situation is expected to hold. This will be discussed elsewhere. References [1] [2] [3] [4] [5] [6] [7] [8]

E. Leutheusser, Phys. Rev. A 29 (1984) 2765. U. Bengtzelius, W. Götze, A. Sjölander, J. Phys. C 17 (1984) 5915. W. Götze, in: J.P. Hansen, D. Levesque, J. Zinn-Justen (Eds.), Liquids, Freezing and Glass Transition, North-Holland, Amsterdam, 1991. R. Kohlrausch, Ann. Phys. 12 (1847) 393. G. Williams, D.C. Watts, Trans. Faraday Soc. 66 (1980) 80. J. Wong, C.A. Angell, Glass Structure by Spectroscopy, Marcel Dekker, New York, 1976. K.L. Ngai, Comment. Solid State Phys. 9 (1979) 127. W. Götze, L. Sjögren, Z. Phys. B 65 (1987) 415.

14 [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

M. Tokuyama / Physica A 526 (2019) 121074 L. Sjögren, W. Götze, J. Non-Cryst. Solids 131 (1991) 153. W. Götze, L. Sjögren, Rep. Progr. Phys. 55 (1992) 241. M. Tokuyama, Physica A 395 (2014) 31. M. Tokuyama, Physica A 430 (2015) 156. M. Tokuyama, Physica A 465 (2017) 229. M. Tokuyama, Physica A 484 (2017) 453. Y. Kimura, M. Tokuyama, IL Nuovo Cimento 39 C (2016) 300. T. Narumi, M. Tokuyama, Phys. Rev. E 95 (2017) 032601. M. Tokuyama, T. Narumi, Physica A 514 (2019) 533. J.K. Percus, G.J. Yevick, Phys. Rev. 110 (1958) 1. M. Tokuyama, Physica A 364 (2006) 23. M. Tokuyama, Physica A 378 (2007) 157. M. Tokuyama, J. Phys. Chem. B 115 (2011) 14030. M. Tokuyama, AIP CP 1518 (2013) 47. M. Tokuyama, S. Enda, J. Kawamura, Physica A 442 (2016) 1. M. Tokuyama, H. Yamazaki, Y. Terada, Phys. Rev. E 67 (2003) 062403. M. Tokuyama, H. Yamazaki, Y. Terada, Physica A 328 (2003) 367. M. Tokuyama, Y. Terada, J. Chem. Phys. B 109 (2005) 21357. M. Tokuyama, T. Narumi, E. Kohira, Physica A 385 (2007) 439. S.K. Kumar, G. Szamel, J.F. Douglas, J. Chem. Phys. 124 (2006) 214501. S. Nakanishi, T. Narumi, Y. Terada, M. Tokuyama, Rep. Inst. Fluid Sci. 19 (2008) 1. T.A. Weber, F.H. Stillinger, Phys. Rev. B 31 (1985) 1954. V.V. Hoang, Phys. Rev. B 70 (2004) 134204. A. Nakano, L.S. Bi, R.K. Kalia, P. Vashishta, Phys. Rev. B 49 (1994) 9441. B.W.H. van Beest, G.J. Kramer, R.A. van Santen, Phys. Rev. Lett. 64 (1990) 1955.