PHYSICAL REVIEW E 88, 042317 (2013)
Formation of double glass in binary mixtures of anisotropic particles: Dynamic heterogeneities in rotations and displacements Kyohei Takae and Akira Onuki Department of Physics, Kyoto University, Kyoto 606-8502, Japan (Received 11 February 2013; revised manuscript received 3 September 2013; published 31 October 2013) We study glass behavior in a mixture of elliptic and circular particles in two dimensions at low temperatures using an orientation-dependent Lennard-Jones potential. The ellipses have a mild aspect ratio (∼1.2) and tend to align at low temperatures, while the circular particles play the role of impurities disturbing the ellipse orientations at a concentration of 20%. These impurities have a size smaller than that of the ellipses and attract them in the homeotropic alignment. As a result, the coordination number around each impurity is mostly 5 or 4 in glassy states. We realize double glass, where both the orientations and the positions are disordered but still hold mesoscopic order. We find a strong heterogeneity in the flip motions of the ellipses, which sensitively depends on the impurity clustering. In our model, a small fraction of the ellipses still undergo flip motions relatively rapidly even at low temperatures. In contrast, the nonflip rotations (with angle changes not close to ±π ) are mainly caused by the cooperative configuration changes involving many particles. Then, there arises a long-time heterogeneity in the nonflip rotations closely correlated with the dynamic heterogeneity in displacements. DOI: 10.1103/PhysRevE.88.042317
PACS number(s): 64.70.Q−, 61.20.Lc, 61.43.Fs, 64.70.P−
I. INTRODUCTION
Much attention has been paid to various types of glass transitions, where the structural relaxations become extremely slow with lowering the temperature T [1,2]. In experiments, colloidal particles can be spherical, but real molecules are mostly nonspherical. The translational and rotational diffusion constants have thus been measured in molecular systems near the glass transition [3]. Using generalized mode-coupling theories, some authors [4–6] have studied the coupled translation-rotation dynamics to predict translational glass and orientational glass. Theoretically, for double glass [6], the translational and orientational degrees of freedom can be simultaneously arrested at the same temperature. In real systems, the molecular rotations sensitively depend on many parameters including the molecular shapes, the density, and the concentration (for mixtures). Mixtures of anisotropic particles with mild differences in sizes and shapes such as (KCN)x (KBr)1−x form a cubic crystal without orientational order (plastic solid) at relatively high T . With further lowering T , they undergo a structural phase transition for small concentrations of impurities (Br− in KCNKBr mixtures) and become orientational glass in nondilute cases [7], where the particles are still in a crystal state. On the other hand, if the two species have significantly different sizes or shapes, translational glass without crystal order can emerge from liquid at low T . For rodlike molecules with relatively large aspect ratios, liquid crystal phase transitions occur with lowering T , but their glass transitions have not yet been well understood. In recent experiments on colloidal ellipsoids in monolayers, the aspect ratio was 6 [8] or 2.1 [9] with considerable size dispersities. In glassy states, these ellipsoids exhibited mesoscopic nematic or smectic order. Molecular dynamics simulations have also been performed on glass-forming fluids composed of anisotropic particles. They can be one-component fluids with a complex internal structure. Examples are methanol [10], orthoterphenyl (OTP) [11,12], and fluids of asymmetric dumbbells [13,14]. There are various kinds of two-component 1539-3755/2013/88(4)/042317(9)
glass formers. The simplest example is a mixture of two species of symmetric dumbbells [15–17]. Recently, we studied a mixture of spheroidal and spherical particles to examine the orientational glass using an orientation-dependent potential [18]. The physical picture of double glass is thus very complex. To give a clear example, we consider a mixture of elliptic particles with a mild aspect ratio and smaller circular particles (impurities). We assume orientation-dependent repulsive and attractive interactions, where the attractive part is between the ellipses and the impurities. Then, the impurities can strongly disturb the orientations and the positions of the surrounding ellipses. This is analogous to hydration of ions by surrounding water molecules [19,20]. If the impurity concentration c is increased from zero, orientational domains and crystalline grains of the ellipses are gradually fragmented and disordered [18]. In this paper, we realize double glass at low T at an impurity concentration of 20%. To produce glassy states, we slowly quench the mixture from liquid. In this situation, we encounter impurity clustering or aggregation at low T , which often results in small crystalline domains of impurities [21,22]. In our model, this tendency is considerably suppressed by the above-mentioned impurityellipse attractive interaction. Nevertheless, the impurity distribution is still mesoscopically heterogeneous, leading to a mesoscopic heterogeneity in the rotational motions. We shall see that some fraction of the ellipses still rotate under weak constraints even at low T . Furthermore, if anisotropic particles have the elliptic symmetry (the spheroidal one in three dimensions), they can undergo flip (turnover) motions with ±π angle changes [10,13,14]. These flip motions can occur thermally for mild aspect ratios, while they are sterically hindered by the surrounding particles for large aspect ratios. Thus, we expect a wide range of the rotational activity for mild aspect ratios. We shall find marked orientational and positional heterogeneities on mesoscopic scales in glass. Such heterogeneous patterns have been visualized in various model systems [23–25]. First, there arises a mesoscopic heterogeneity of
042317-1
©2013 American Physical Society
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013)
the flip motions correlated with the impurity clustering. Second, the positional configuration changes cause nonflip rotations of the ellipses, which are the origin of the longtime decay of the rotational correlation functions G (t) of even [10,12–17]. It follows a dynamic heterogeneity of the long-time nonflip rotations correlated with the dynamic heterogeneity in displacements or bond breakage [26–29]. The organization of this paper is as follows. In Sec. II, we will explain our simulation model and method. In Sec. III A, we will present simulation results on the heterogeneities in the orientations and the positions. In Sec. III B, the timecorrelation functions will be examined. In Sec. III C, the angular and translational mean-square displacements will be calculated. In Sec. III D, we will introduce the flip number for each ellipse in a time interval and study its heterogeneity. In Sec. III E, time development of a configurational change with large displacements and/or large angle changes will be illustrated. II. MODEL AND NUMERICAL METHOD
In two dimensions, we consider mixtures of anisotropic and circular particles with numbers N1 and N2 , where N = N1 + N2 = 4096. The concentration of the circular species is c = N2 /N. The particle positions are written as r i (i = 1, . . . ,N). The orientation vectors of the anisotropic particles are expressed as ni = (cos θi , sin θi ) in terms of angles θi (i = 1, . . . ,N1 ). The pair potential Uij between particles i ∈ α and j ∈ β (α,β = 1,2) is a modified Lennard-Jones potential given by [18] 6 12 σαβ σαβ Uij = 4 (1 + Aij ) 12 − (1 + Bij ) 6 , (1) rij rij where rij is the particle distance and is the interaction energy. In terms of characteristic lengths σ1 and σ2 , we set σαβ = (σα + σβ )/2. The potential is truncated at rij = 3σ1 . The particle anisotropy is accounted for by the anisotropic factors Aij and Bij , which depend on the angles between ni , nj , and the relative direction rˆ ij = rij−1 (r i − r j ). In this paper, we set Aij = χ [δα1 (ni · rˆ ij )2 + δβ1 (nj · rˆ ij )2 ],
(2)
Bij = ζ [δα1 δβ2 (ni · rˆ ij )2 + δα2 δβ1 (nj · rˆ ij )2 ],
(3)
where δαβ is the Kronecker δ, χ is the anisotropy strength of repulsion, and ζ is that of attraction between the two species. The Newton equations for r i (t) and θi (t) are given by ∂U d2 ri = − , (4) dt 2 ∂ ri d2 ∂U I 2 θi = − , (5) dt ∂θi where U = i 6, k = 6, and k < 6 and those of the impurities with k > 5, k = 5, and k < 5. These six fractions are denoted by φ1> , φ16 , φ1< , φ2> , φ25 , and φ2< , respectively, as functions of T . At low T , k is mostly 6 or 7 for the ellipses and is mostly 4 or 5 for the impurities. In fact, for the data in Figs. 1 and 2 at T = 0.05, we have (φ16 ,φ1> ) ∼ = (0.66,0.33) and (φ2< ,φ25 ) ∼ = (0.28,0.71) for ζ = 1, while these sets are (0.73,0.26) and (0.05,0.94), respectively, for ζ = 0.5. In the right panel of Fig. 3, φ1< and φ2> are very small at low T and increase with increasing T . Thus, the ellipses with k < 6 and the impurities with k > 5 represent liquidlike defects [35]. Hentschel et al. [35] studied the positional disorder using the Voronoi graphs for a mixture of circular particles with the soft-core potential in two dimensions. In their simulation, small (large) particles enclosed by heptagons (pentagons) form liquidlike defects decreasing at low T . B. Time-correlation functions
For strong short-range anchoring, the rotational dynamics sensitively depends on whether the anisotropic particles are close or far from the impurities. In Fig. 4, we show time
042317-3
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013) 0.8
0.15
0.8
0.6
0.4
0.12
φ5
0.6
0.06 0.2
0
0.03
0
0.2
0.4
0.6
0.8
1
0
ζ=1.0 T=0.2
0.4
0.2 0
0.2
0.4
T
0.6
0.8
1
T φ1> ,
0 -1
φ16 ,
-2/3
-1/3
evolution of the angle changes, θi (t0 ,t + t0 ) = θi (t + t0 ) − θi (t0 ),
(11)
where we pick up a rapidly rotating ellipse, a rarely flipping one, and an inactive one. We can see instantaneous flip motions by ±π . In previous simulations in glassy states, they observed flips for rodlike molecules [10,13,14] and large angle jumps for ortho-terphenyl (OTP) [11,12]. We introduce the distribution of the angle changes, 1 δ {[θi (t)]2π − θ }, N1 i∈1
(12)
where θi (t) = θi (t0 ,t + t0 ) and −π θ < π. For any angle ϕ, we define [ϕ]2π = ϕ − 2pπ , which is in the range [−π,π ] with an integer p. Furthermore, we consider the th
1/3
0
θ/
φ1
6, k = 6, and k < 6, respectively, and those φ2> , φ25 , and φ2< of the impurities with k > 5, k = 5, and k < 5, respectively, as functions of T with ζ = 1, where k is the number of the surrounding triangles in the Delaunay diagrams. Here, φ16 , φ1> , φ2< , and φ25 are large at any T , but φ1< and φ2> decrease at low T . Right: Fractions φ1< and φ2> vs T , which are the fractions of liquidlike defects [35] decreasing at low T .
G(t,θ ) =
t = 40 400 4000
σ2 / σ1 = 0.6
0.09
G(t, θ)
φ6
φ5
2/3
1
π
FIG. 5. (Color online) Time-dependent angle distribution function G(t,θ ) in Eq. (12) at t = 40, 400, and 4000 for ζ = 1 and T = 0.2. Peaks emerge at θ = ±π due to flip motions on the time scale of τ1 = 400. Afterwards, G(t,θ ) → 1/2π on the time scale of τ2 = 24 000.
moments of G(t,θ ) given by π G (t) = dθ G(t,θ ) cos(θ ) −π
1 cos[θi (t)]. = N1 i∈1
(13)
We calculated G(t,θ ), G1 (t), and G2 (t) by taking the average · · · over the initial time t0 and over five runs. In Fig. 5, we show time evolution of G(t,θ ), where flip motions give rise to peaks at θ = ±π growing on the time scale of τ1 . Thus, these flip motions cause the decay of G1 (t) in Fig. 6(a). However, G2 (t) is unchanged by the turnovers and
c=0.2, χ=1.2, σ2 / σ1 = 0.6, ζ=1.0 1
1
(a)
(b) 0.8
0.8
G1(t) 6
0.6
0.6
5
0.4
0.4
G2(t) σ2 / σ1 = 0.6
4
ζ=1.0 T=0.2
Δθ / π
3 2
0.2
T=0.05 0.10 0.15 0.20 0.25 0.30
0.2
T=0.05 0.10 0.15 0.20 0.25 0.30
0 -2 0 -2 -1 0 1 2 3 4 5 -1 0 1 2 3 4 5 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
1
1
t
t
3
6
1
0
(c)
5
(d)
0.8
-1
4
2
-2
0.6
-3 0.4
-4 0
200
400
600
t
0.2
FIG. 4. (Color online) Time evolution of angle changes θi (t0 ,t0 + t) in Eq. (11) for (1) a frequently flipping ellipse, (2) an infrequently flipping one, and (3) an inactive one for ζ = 1 and T = 0.1. Flip events occur at points (◦) on the curves (see the Appendix). These jumps are very different from thermal vibrations but occur thermally.
3
Fs (q,t)
log10(τ)
2
T=0.05 0.10 0.15 0.20 0.25 0.30
τ1 τ2 τα
1 0
-1 0 -2 10 10-1 100 101 102 103 104 105
t
2
4
6
8
10
1/T
FIG. 6. (Color online) (a) G1 (t), (b) G2 (t), (c) Fs (q,t) at q = 2π for ellipses at six T . (d) Relaxation times τ1 , τ2 , and τα in Eqs. (14)– (16) vs 1/T . Here, ζ = 1 and time t is in units of τ0 in Eq. (8).
042317-4
FORMATION OF DOUBLE GLASS IN BINARY MIXTURES . . .
G1 (τ1 ) = 1/e,
(14)
G2 (t) ∝ exp[−(t/τ2 )β ] (t > 1),
(15)
Fs (q,t) ∝ exp[−(t/τα )γ ]
(t > 1),
(16)
where the exponents β and γ are about 0.4 for T 0.2. We here determine τ2 and τα from the long-time relaxations of G2 (t) and Fs (q,t), respectively. In Fig. 6(d), we plot them, where τ1 τα ∼ τ2 . In all the T range in Fig. 6(d), τ1 may be nicely fitted to the Arrhenius form, ln(τ1 ) = 1.4/T − 1.2. On the other hand, τ2 and τα exhibit a changeover at T ∼ 0.2 and can be fitted to the Arrhenius forms as ln(τ2 ) = 2.6/T − 3.0 and ln(τα ) = 2.6/T − 4.5, for T 0.2, so τ2 /τα ∼ = 4 for T 0.2. Thus, G1 (t) decays mainly due to thermally activated flip motions, which are nearly decoupled from the translational motions. On the other hand, G2 (t) and Fs (q,t) decay at longer times due to irreversible configuration changes involving at least several particles. In three dimensions, the distribution of the angles cos−1 [ni (t0 + t) · ni (t0 )] was calculated for OTP [11] and for dumbbells [13,14]. In these papers, this distribution exhibited peaks due to orientational jumps. Also in the rotational time-correlation functions in three dimensions, the Legendre polynomials P [ui (t)] with ui (t) = ni (t0 + t) · ni (t0 ) were used [10,12–17], where G (t) with even decayed slower than those with odd at low T . These previous findings are in accord with our results. C. Mean-square displacements
In the literature, the angular mean-square displacement has been calculated to study the rotational diffusion [12–15,17]. In two dimensions, it is defined by 1 |θi (t0 ,t0 + t)|2 . (17) Mθ (t) = |θ |2 = N1 i∈1 We also introduce the usual positional mean-square displacement for the ellipses by 1 M(t) = |r|2 = |r i (t0 ,t0 + t)|2 , (18) N1 i∈1 where r i (t0 ,t0 + t) = r i (t0 + t) − r i (t0 ). At very short times, these quantities exhibit the ballistic behavior (∝ t 2 ). At long times, they grow linearly in time as (19) Mθ (t) ∼ = 2DR t, M(t) ∼ = 4Dt,
(20)
where DR and D are the rotational and translational diffusion constants, respectively. If the rotational activity is strongly heterogeneous, the so-called Stokes-Einstein-Debye relation DR ∼ kB T /ηa does not hold [3,12], where η is the viscosity and a is the radius of the diffusing particle. In our case, Mθ (t)
ζ=1.0, T=0.2, t f =8000 100
(a) 102 10
10-1
(b)
=1
n > 100 total
pe
slo
10-2
< |Δr | > 2
total 1
104
0
sl o pe =
decays more slowly after the initial relaxation in Fig. 6(b). Notice that G(t,θ ) tends to 1/2π on the time scale of τ2 . In Fig. 6(c), we also show the self-part of the density timecorrelation function Fs (q,t) at q = 2π for the ellipses, which closely resembles G2 (t). We define the relaxation times τ1 , τ2 , and τα as
PHYSICAL REVIEW E 88, 042317 (2013)
10-3 10-2 10
n= 0 n < 10
10-4
-4
10-2
100
t
102
104
10
|Δr | > 0.6
-5
10-2
100
t
102
104
FIG. 7. (Color online) Angular and positional mean-square displacements for ζ = 1 and T = 0.2. (a) Angular one Mθ (t) in Eq. (17) and contributions from those with ni 100, ni 10, and ni = 0, where ni is the flip number for tf = 20τ1 = 8000 (see the Appendix). The contribution from ni 100 approaches Mθ (t) for t 1, leading to DR = 0.14. (b) Positional one M(t) in Eq. (18) and contribution from those with ri > 0.6 in Eq. (21), where the latter grows linearly for t 20 with D = 1.4 × 10−5 .
is greatly increased by rapidly flipping ellipses and DR from it does not correspond to any experimentally observed relaxation rates at low T . In Fig. 7(a), we plot Mθ (t) for ζ = 1 and T = 0.2. Here, it attains the diffusion behavior with DR = 0.14 for t 1, while G1 (t) decays slower with τ1 = 400. We also display the contributions to the sum in Mθ (t) in Eq. (17) from the ellipses with ni = 0, ni 10, and ni 100, where ni is the flip number of ellipse i in a time interval with width tf = 8000 (see the Appendix). The fractions of these three groups are 0.16, 0.34, and 0.39, respectively. Remarkably, the contribution from ni 100 approaches Mθ (t) for t 1, while that from ni 10 behaves diffusively as 0.7 × 10−3 × 2t for t τ1 . Thus, the effective rotational diffusion constant of the ellipses with n 10 is 0.7 × 10−3 /0.34 = 2 × 10−3 . The Mθ (t) itself exhibits the plateau behavior at much lower temperatures (say, T = 0.05), while the contributions from n = 0 and n 10 exhibit it at T = 0.2. In Fig. 7(b), for ζ = 1 and T = 0.2, M(t) still in the course of plateau-to-diffusion crossover even at t = 104 . To obtain small D, we also plot the contribution from the ellipses with large displacements [36], 1 M > (t) = [ri (t) − c ]|ri (t)|2 , (21) N1 i∈1 where ri (t) is an abbreviation of |r i (t0 ,t0 + t)| and (u) is the step function being equal to 1 for u 0 and to 0 for u < 0. The threshold length c is set equal to 0.6. In this restricted sum, the thermal vibrational motions within transient cages are excluded, so it picks up the thermally activated jumps only. As a result, we have the linear growth M > (t) ∼ = 4Dt with D = 1.4 × 10−5 ∼ 10−4 DR from the early stage t 20. This behavior of M > (t) is insensitive to a small change of c [36]. For example, almost the same results followed for c = 0.8. As a similar observation, Chong and Kob found that DR τ2 grows strongly with lowering T for a mixture of rigid dumbbell molecules [17]. See more discussions for other molecular systems in item (i) in the Summary.
042317-5
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013)
We further introduce the n-dependent mean-square displacement among the ellipses with n flips as 1 |θ |2 n (t) = δnni |θi (t0 ,t0 + t)|2 . (23) f N1 φn i∈1
ζ=1.0, T=0.2, t f =8000 1
(b)
(a)
0.8
2
0.6
10nφnf 1
10
-2
2
10 0
1
50
(c)
0.4
< Δθ >n /π 2
2
0.2
φnf
G150 G120
It follows the sum relation Mθ (t) = φnf |θ |2 n (t).
G110
100
(d)
t
102
In our case, most ellipses undergo +π flips and −π flips equally on long times, so i∈1 δnni θi (t0 ,t0 + t) = 0. f f In Fig. 8(a), we plot φn , nφn , and |θ |2 n (tf ) (where t = tf ) for ζ = 1 and T = 0.2. Here, the fraction of the ellipses f with no flip is given by φ0 = 0.16. We find that the flip number distribution is very broad as
104
φnf ∼ 0.1n−1 ,
n < 10
n > 200
D. Distribution of flip numbers
The large size of DR is due to the presence of ellipses frequently undergoing flip motions. In the Appendix, we will give a method of determining the flip number ni for each ellipse i in a time interval [t0 ,t0 + tf ]. We should choose a sufficiently large width tf to detect a wide range of ni . In the following, tf = 8000 at T = 0.2 in Figs. 7 and 8 and tf = 105 at T = 0.05 in Fig. 9. The curves in Figs. 8(a), 8(b), and 9(a) are the averages over six runs. For a given time interval [t0 ,t0 + tf ], the fraction of the ellipses with n flips is written as δnni /N1 . (22) φnf = i∈1
1
ζ=1.0, T=0.05, t f =105
(a)
(b)
φ0f = 0.86
G1(t f )=0.82
0.8 0.6 0.4
10 φnf 2
10
-2
< Δθ >n /πf 2
2
10nφn
0.2 0
50
n
100
(25)
in the range 1 n < nmax , where nmax is an upper bound f about 103 . In the present case, the sums of φn in the ranges 1 n 10, 11 n 99, and n 100 are 0.18, 0.27, and 0.39, respectively. Furthermore, we find
FIG. 8. (Color online) (a) 102 φnf , 10nφnf for n 1, and −2 10 |θ |2 n /π 2 in Eq. (23), where tf = 20τ1 = 8000 (averages over six runs). (b) G1 (t) in Eq. (13) and Gn1 (t) in Eq. (28) with n = 0, 10, 20, and 50 at t = tf , where the latter approach the former at long times. Snapshots of ellipses with ni 10 in (c) and those with ni 200 in (d), whose heterogeneities are correlated with the impurity clustering.
1.2
(24)
n0
G10
0 100 10-2
n
G1
n < 10
FIG. 9. (Color online) 102 φnf , 10nφnf , and 10−2 |θ |2 n /π 2 (averages over six runs) in (a) and snapshot of ellipses with n 10 in f (b), where ζ = 1, T = 0.05, and tf = 105 . At this low T , φ0 = 0.86 and the fraction of the depicted ellipses in (b) is 0.92.
|θ |2 n (t) ∼ π 2 nt = π 2 nt/tf
(26)
for t 1 and n 1. In a general time width t, the ellipses with n flips in the reference time width tf should flip nt = nt/tf times on the average, where t 1 and n 1. Then, together with the sentence below Eq. (24), Eq. (26) is a natural relation. From Eqs. (24)–(26) we find DR ∼ nmax /tf ,
(27)
which means that DR is determined by rapidly rotating ellipses. To be self-consistent, nmax should be proportional to tf ; then, DR is independent of tf . In Fig. 8(b), we compare G1 (t) and the restricted sums, 1 Gn1 (t) = (n − ni ) cos[θi (t0 ,t0 + t)], (28) N1 i∈1 where we set n = 0, 10, 20, and 50. We here pick up the ellipses with flip numbers not exceeding n owing to the step function . We can see that these Gn1 (t) are nearly constant for some time and become nearly equal to G1 (t) after long times. For t 103 , G1 (t) is composed of the contributions from the ellipses with n 10. Note that the ellipses with n flips for tf = 20τ1 should undergo n/20 flips for tf = τ1 on the average. Thus, the ellipses which have undergone more than a few flips in the time interval [t0 ,t0 + τ1 ] do not appreciably contribute to G1 (τ1 ) = 1/e. In Figs. 8(c) and 8(d), we show snapshots of the ellipses with n 10 and n 200, respectively. The distributions of these rotationally inactive and active ellipses are highly heterogeneous. This marked feature is rather natural in view of the mild aspect ratio 1.23 and the significant impurity clustering. In fact, the impurities are nearly absent in the red regions in Fig. 8(b). In Fig. 9, we also show that the flip motions still remain even at T = 0.05. In this case, we find G1 (t) ∼ 0.8 at t = 105 in Fig. 6(a), but we estimate τ1 ∼ 1012 from the extrapolation of the Arrhenius form [see the sentences below Eq. (16)]. f f In Fig. 9(a), we find φ0 = 0.86 and φn ∼ 0.02n−1 and again
042317-6
FORMATION OF DOUBLE GLASS IN BINARY MIXTURES . . .
obtain Eq. (26) for ζ = 1, T = 0.05, and tf = 105 . In Fig. 9(b), displayed is a snapshot of the ellipses with ni 10, whose f fraction is n10 φn = 0.92. Even at this low T , 2% ellipses have ni > 200. We have DR = 1.2 × 10−3 due to these rapidly flipping ellipses. This snapshot was produced by the initial particle configuration common to that in Fig. 8(b). Most of the ellipses in the red regions in Fig. 8(b) are now inactive, since their orientation alignment increases with lowering T . E. Dynamic heterogeneities in nonflip rotations and displacements
Now, we examine the long-time structural relaxation caused by collective configuration changes, where large displacements induce large nonflip rotations. These events should occur repeated in the same fragile regions on time scales longer than τ2 [29,36]. In Fig. 10, we illustrate time development of a configuration change at successive times t0 + t with 0 t 10. Depicted are the ellipses with ri (t) = |r i (t0 ,t0 + t)| > 0.6 or π
θ 0
t=0
ζ =1.0, T=0.2
t=2
PHYSICAL REVIEW E 88, 042317 (2013)
ci (t) < 0.2, where ci (t) = cos[2θi (t0 ,t0 + t)].
The condition ci (t) < 0.2 means 0.22π < |θi (t0 ,t0 + t)| < 0.78π in the range [−π,π ]. From Eq. (13) we have G2 (t) = c i∈1 i (t)/N1 . In the narrow region in Fig. 10, the particle configuration was nearly stationary for t 0, but large particle motions started for t > 0 and continued on a time scale of 10. We can see circulating particle motions at t = 6 and 8 and stringlike ones at t = 10 [27–29,36]. The orientations of these ellipses are largely changing with their movements. For t > 10, the subsequent displacements became small, but considerable orientational motions persisted until t ∼ 20. We examine the nonflip motions in terms of ci (t) in Eq. (29), since it is invariant with respect to turnovers. In Fig. 11(a), we visualize the correlation between ci = ci (t) and ri = ri (t) for ζ = 1 and T = 0.2. We set t = tf = 8000, which is onethird of τ2 ∼ 24 000 (∼4τα ). We present a snapshot of the ellipses with (a) ci < 0.2 and ri > 0.6, (b) ci < 0.2 and ri < 0.6, and (c) ci > 0.2 and ri > 0.6. Here, we exclude the ellipses with ni > 200(∼30%), because they do not exhibit the glassy behavior. The fractions of these depicted groups are (a) 0.08, (b) 0.10, and (c) 0.09, while the fraction of the ellipses with ci > 0.2, ri < 0.6, and ni 200 is 0.41. Thus, if we consider the ellipses with ci < 0.2 and ni < 200, half of them have undergone displacements with ri > 0.6. Also, if we consider the ellipses with ri > 0.6 and ni < 200, half of them have undergone large angle changes with ci < 0.2. In Fig. 11(b), displayed is a snapshot of the ellipses with (a) ri > 0.6 and ni < 200 (18%) and (b) ri > 0.6 and ni > 200 (11%). The group (a) here consists of the groups (a) and (b) in Fig. 11(a). Here, the ellipses in these two groups form clusters, indicating collective displacements. In addition, clusters of one group are adjacent to those of another group. Thus, rotationally active ellipses with large ni tend to be translationally active also. (a)
t=4
ζ=1.0, T=0.2, t=8000
(b)
t=6
n200, |Δr | > 0.6
FIG. 11. (Color online) Comparison of ellipses with ci ≡ cos[2θi (t0 ,t0 + t)] < 0.2 and particles with large displacement ri ≡ |r i (t0 ,t0 + t)| > 0.6 for ζ = 1 and T = 0.2 at t = tf = 8000. In the left, depicted are three groups of ellipses with ci < 0.2 and ri > 0.6 (in red), ci < 0.2 and ri < 0.6 (in green), ci > 0.2 and ri > 0.6 (in yellow), whose fractions are 0.08, 0.10, and 0.09, respectively. In the right, depicted are 0.18N1 ellipses with ri > 0.6 and ni < 200 (in red) and 0.11N1 ones with ri > 0.6 and ni > 200 (in green).
042317-7
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013)
IV. SUMMARY AND REMARKS
We have performed simulation of a mixture of elliptic particle with a mild aspect ratio (=1.23) and smaller circular impurities with σ2 /σ1 = 0.6 at 20%. We have assumed an angle-dependent attractive interaction between the ellipses and the impurities (∝ζ ), which leads to the homeotropic anchoring of the ellipses around each impurity. We summarize our main simulation results. (1) We have shown snapshots of the orientations and the positions in Figs. 1 and 2, which are mesoscopically heterogeneous. From the Delaunay triangulation in the right panels of Fig. 2, we have found that the number of surrounding triangles (the coordination number) is 6 or 7 for the ellipses and 5 or 4 for the impurities in glassy states, as plotted in Fig. 3. A majority of the impurities (∼70%) are surrounded by five ellipses, analogously to the case of the Shintani-Tanaka model [24]. (2) We have calculated the distribution function of the angle changes G(t,θ ) in Eq. (12), which exhibits peaks at θ = ±π for large t due to flip motions as in Fig. 5. We have found that the rotational time-correlation functions G1 (t) and G2 (t) of the ellipses relax very differently at long times in Fig. 6, because G1 (t) decays due to flip motions and G2 (t) due to configuration changes. (3) We have found that the angular mean-square displacement Mθ (t) in Eq. (17) behaves as 2DR t rapidly for t 1 with very large DR in Fig. 7. This is in marked contrast to the slow time evolution of the translational mean-square displacement M(t). However, the contribution to M(t) from the largely displaced ellipses (|ri | > 0.6) has exhibited the diffusion behavior with D = 10−4 DR , because the diffusion is governed by the activation dynamics [36]. f (4) We have displayed the fractions φn of the ellipses with n flips in a time interval [t0 ,t0 + tf ], where tf = 8000 at T = 0.2 in Figs. 7 and 8 and tf = 105 at T = 0.05 in Fig. 9. We have f found a very broad distribution φn (∝ n−1 ) for 1 n < nmax . The angular mean-square displacement |θ |2 n (t) among the ellipses with n flips behaves as π 2 nt/tf . Then DR ∼ nmax /tf due to rapidly flipping ellipses. We have also shown that the long-time decay of G1 (t) is determined by rotationally inactive ellipses. (5) The distributions of the rotationally inactive (n 10) and active (n 200) ellipses at T = 0.2 have been presented in Figs. 8(c) and 8(d), while that of the inactive ones at T = 0.05 has been given in Fig. 9(b). (6) We have illustrated time development of a configuration change in Fig. 10, where large displacements and nonflip rotations are coupled. We have demonstrated close correlation between nonflip rotations and large displacements in Fig. 11. We make remarks as follows. (i) Our potential energy is invariant with respect to turnovers. This means that our model has the quadrupolar symmetry [7]. This is also the case of diatomic molecules or dumbbells [16,17], for which the flip motions should be highly heterogeneous for mild aspect ratios. For methanol, Sindzingre and Klein [10] found flip motions near the glass transition. For OTP, Lewis and Wahnstr¨om [11] found translation-free orientational jumps, while Lombardo et al. [12] found an enhancement in the rotational motions relative to the translation motions at low T . For these systems, the
heterogeneity of orientational jumps should be examined in more detail. However, the physical picture for one-component glass-formers should be different from that for binary mixtures in this paper. (ii) We have chosen a mild aspect ratio (=1.23) to find significant flip motions. However, an increase in the aspect ratio leads to a decrease in the flip frequency, on which we will report shortly. (iii) We have suppressed the clustering of small impurities by the angle-dependent attractive interaction. If we consider large impurities (say, σ2 /σ1 = 1.4), we may realize double glass by adding a repulsive interaction among the impurities suppressing crystal formation. (iv) The phase behavior of mixtures of two species of anisotropic particles should be studied in the future, where we expect nematic or smectic glass. (v) The spatial scales of the structural and orientational heterogeneities depend on various parameters. If the oriented domains are not too small, there arises a large orientationstrain coupling, leading to soft elasticity and a shape-memory effect [18]. These effects were observed for Ti-Ni alloys [37] (where atomic displacements within unit cells cause structural changes). When anisotropic particles have electric dipoles [7], mesoscopic polar domains appear as in ferroelectric glass (relaxors) [38]. Including such metallic alloys also, we point out relevance of the compositional heterogeneity in the glass transition coupled with some orientational phase transition. ACKNOWLEDGMENTS
This work was supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan. K.T. was supported by the Japan Society for Promotion of Science. The present numerical calculations were carried out on SR16000 at YITP in Kyoto University. APPENDIX: FLIP EVENTS IN NUMERICAL ANALYSIS
In our numerical analysis, we determine a series of flip times, ti1 ,ti2 ,ti3 , . . . for each ellipse i. We write the angle change as θi (t) = θi (t0 ,t + t0 ), suppressing t0 . (i) At the first flip time ti1 we set |θi (ti1 )| = 2π/3.
(A1)
For t > ti1 we introduce θi1 (t) = θi (t) ± π,
(A2)
where +π or −π is chosen such that |θi1 (ti1 + 0)| < π/2. (ii) At the second flip time ti2 we set |θi1 (ti2 )| = 2π/3.
(A3)
For t > ti2 we again introduce θi2 (t) = θi1 (t) ± π,
(A4)
where +π or −π is chosen such that |θi2 (ti2 + 0)| < π/2. (iii) We repeat these procedures to obtain the successive flip times. See Fig. 4 for examples of the flip time series.
042317-8
FORMATION OF DOUBLE GLASS IN BINARY MIXTURES . . .
PHYSICAL REVIEW E 88, 042317 (2013)
Note that the threshold 2π/3 in Eqs. (A1) and (A3) may be changed to another angle, say 5π/6. However, the resultant
flip time series is rather insensitive to this choice as long as it is in the range [π/4,π/2].
[1] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, and S. W. Martin, J. Appl. Phys. 88, 3113 (2000). [2] K. Binder and W. Kob, Glassy Materials and Disordered Solids (World Scientific, Singapore, 2005). [3] F. Fujara, B. Geil, H. Sillescu, and G. Fleischer, Z. Phys. B 88, 195 (1992); M. T. Cicerone and M. D. Ediger, J. Chem. Phys. 104, 7210 (1996). [4] A. Winkler, A. Latz, R. Schilling, and C. Theis, Phys. Rev. E 62, 8004 (2000). [5] S.-H. Chong and W. G¨otze, Phys. Rev. E 65, 041503 (2002). [6] R. Zhang and K. S. Schweizer, J. Chem. Phys. 133, 104902 (2010); 136, 154902 (2012). [7] U. T. H¨ochli, K. Knorr, and A. Loidl, Adv. Phys. 39, 405 (1990). [8] Z. Zheng, F. Wang, and Y. Han, Phys. Rev. Lett. 107, 065702 (2011). [9] C. K. Mishra, A. Rangarajan, and R. Ganapathy, Phys. Rev. Lett. 110, 188301 (2013). [10] P. Sindzingre and M. L. Klein, J. Chem. Phys. 96, 4681 (1992). [11] L. J. Lewis and G. Wahnstr¨om, Phys. Rev. E 50, 3865 (1994); J. Non-Cryst. Solids 172, 69 (1994). [12] T. G. Lombardo, P. G. Debenedetti, and F. H. Stillinger, J. Chem. Phys. 125, 174507 (2006). [13] S. K¨ammerer, W. Kob, and R. Schilling, Phys. Rev. E 56, 5450 (1997). [14] C. De Michele and D. Leporini, Phys. Rev. E 63, 036702 (2001). [15] S.-H. Chong, A. J. Moreno, F. Sciortino, and W. Kob, Phys. Rev. Lett. 94, 215701 (2005). [16] A. J. Moreno, S.-H. Chong, W. Kob, and F. Sciortino, J. Chem. Phys. 123, 204505 (2005). [17] S.-H. Chong and W. Kob, Phys. Rev. Lett. 102, 025702 (2009). [18] K. Takae and A. Onuki, Europhys. Lett. 100, 16006 (2012). [19] J. N. Israelachvili, Intermolecular and Surface Forces (Academic, London, 1991).
[20] C. A. Angell and E. J. Sare, J. Chem. Phys. 49, 4713 (1968); M. Kobayashi and H. Tanaka, J. Phys. Chem. B 115, 14077 (2011). [21] G. A. Vliegenthart, A. van Blaaderen, and H. N. W. Lekkerkerker, Faraday Discuss. 112, 173 (1999). [22] D. Antypov and D. J. Cleaver, J. Chem. Phys. 120, 10307 (2004). [23] T. Hamanaka and A. Onuki, Phys. Rev. E 74, 011506 (2006); 75, 041503 (2007). [24] H. Shintani and H. Tanaka, Nat. Phys. 2, 200 (2006). [25] T. Kawasaki, T. Araki, and H. Tanaka, Phys. Rev. Lett. 99, 215701 (2007). [26] R. Yamamoto and A. Onuki, J. Phys. Soc. Jpn. 66, 2545 (1997); Phys. Rev. E 58, 3515 (1998). [27] C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 80, 2338 (1998). [28] S. C. Glotzer, J. Non-Cryst. Solids 274, 342 (2000). [29] H. Shiba, T. Kawasaki, and A. Onuki, Phys. Rev. E 86, 041504 (2012). [30] S. Nos´e, Mol. Phys. 52, 255 (1984). [31] J. G. Gay and B. J. Berne, J. Chem. Phys. 74, 3316 (1981). [32] A. Ben-Naim, J. Chem. Phys. 54, 3682 (1971); K. A. T. Silverstein, A. D. J. Haymet, and K. A. Dill, J. Am. Chem. Soc. 120, 3166 (1998). [33] J. M. Drouffe, A. C. Maggs, and S. Leibler, Science 254, 1353 (1991); H. Noguchi, J. Chem. Phys. 134, 055101 (2011). [34] D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979). [35] H. G. E. Hentschel, V. Ilyin, N. Makedonska, I. Procaccia, and N. Schupper, Phys. Rev. E 75, 050404(R) (2007). [36] T. Kawasaki and A. Onuki, Phys. Rev. E 87, 012312 (2013); this paper shows that Dτb is nearly independent of T in a supercooled fragile binary mixture in three dimensions, where τb is the bond breakage time [26]. [37] S. Sarkar, X. Ren, and K. Otsuka, Phys. Rev. Lett. 95, 205702 (2005); Y. Wang, X. Ren, and K. Otsuka, ibid. 97, 225703 (2006). [38] R. A. Cowley, S. N. Gvasaliya, S. G. Lushnikov, B. Roessli, and G. M. Rotaru, Adv. Phys. 60, 229 (2011).
042317-9
Formation of double glass in binary mixtures of anisotropic particles: Dynamic heterogeneities in rotations and displacements Kyohei Takae and Akira Onuki Department of Physics, Kyoto University, Kyoto 606-8502, Japan (Received 11 February 2013; revised manuscript received 3 September 2013; published 31 October 2013) We study glass behavior in a mixture of elliptic and circular particles in two dimensions at low temperatures using an orientation-dependent Lennard-Jones potential. The ellipses have a mild aspect ratio (∼1.2) and tend to align at low temperatures, while the circular particles play the role of impurities disturbing the ellipse orientations at a concentration of 20%. These impurities have a size smaller than that of the ellipses and attract them in the homeotropic alignment. As a result, the coordination number around each impurity is mostly 5 or 4 in glassy states. We realize double glass, where both the orientations and the positions are disordered but still hold mesoscopic order. We find a strong heterogeneity in the flip motions of the ellipses, which sensitively depends on the impurity clustering. In our model, a small fraction of the ellipses still undergo flip motions relatively rapidly even at low temperatures. In contrast, the nonflip rotations (with angle changes not close to ±π ) are mainly caused by the cooperative configuration changes involving many particles. Then, there arises a long-time heterogeneity in the nonflip rotations closely correlated with the dynamic heterogeneity in displacements. DOI: 10.1103/PhysRevE.88.042317
PACS number(s): 64.70.Q−, 61.20.Lc, 61.43.Fs, 64.70.P−
I. INTRODUCTION
Much attention has been paid to various types of glass transitions, where the structural relaxations become extremely slow with lowering the temperature T [1,2]. In experiments, colloidal particles can be spherical, but real molecules are mostly nonspherical. The translational and rotational diffusion constants have thus been measured in molecular systems near the glass transition [3]. Using generalized mode-coupling theories, some authors [4–6] have studied the coupled translation-rotation dynamics to predict translational glass and orientational glass. Theoretically, for double glass [6], the translational and orientational degrees of freedom can be simultaneously arrested at the same temperature. In real systems, the molecular rotations sensitively depend on many parameters including the molecular shapes, the density, and the concentration (for mixtures). Mixtures of anisotropic particles with mild differences in sizes and shapes such as (KCN)x (KBr)1−x form a cubic crystal without orientational order (plastic solid) at relatively high T . With further lowering T , they undergo a structural phase transition for small concentrations of impurities (Br− in KCNKBr mixtures) and become orientational glass in nondilute cases [7], where the particles are still in a crystal state. On the other hand, if the two species have significantly different sizes or shapes, translational glass without crystal order can emerge from liquid at low T . For rodlike molecules with relatively large aspect ratios, liquid crystal phase transitions occur with lowering T , but their glass transitions have not yet been well understood. In recent experiments on colloidal ellipsoids in monolayers, the aspect ratio was 6 [8] or 2.1 [9] with considerable size dispersities. In glassy states, these ellipsoids exhibited mesoscopic nematic or smectic order. Molecular dynamics simulations have also been performed on glass-forming fluids composed of anisotropic particles. They can be one-component fluids with a complex internal structure. Examples are methanol [10], orthoterphenyl (OTP) [11,12], and fluids of asymmetric dumbbells [13,14]. There are various kinds of two-component 1539-3755/2013/88(4)/042317(9)
glass formers. The simplest example is a mixture of two species of symmetric dumbbells [15–17]. Recently, we studied a mixture of spheroidal and spherical particles to examine the orientational glass using an orientation-dependent potential [18]. The physical picture of double glass is thus very complex. To give a clear example, we consider a mixture of elliptic particles with a mild aspect ratio and smaller circular particles (impurities). We assume orientation-dependent repulsive and attractive interactions, where the attractive part is between the ellipses and the impurities. Then, the impurities can strongly disturb the orientations and the positions of the surrounding ellipses. This is analogous to hydration of ions by surrounding water molecules [19,20]. If the impurity concentration c is increased from zero, orientational domains and crystalline grains of the ellipses are gradually fragmented and disordered [18]. In this paper, we realize double glass at low T at an impurity concentration of 20%. To produce glassy states, we slowly quench the mixture from liquid. In this situation, we encounter impurity clustering or aggregation at low T , which often results in small crystalline domains of impurities [21,22]. In our model, this tendency is considerably suppressed by the above-mentioned impurityellipse attractive interaction. Nevertheless, the impurity distribution is still mesoscopically heterogeneous, leading to a mesoscopic heterogeneity in the rotational motions. We shall see that some fraction of the ellipses still rotate under weak constraints even at low T . Furthermore, if anisotropic particles have the elliptic symmetry (the spheroidal one in three dimensions), they can undergo flip (turnover) motions with ±π angle changes [10,13,14]. These flip motions can occur thermally for mild aspect ratios, while they are sterically hindered by the surrounding particles for large aspect ratios. Thus, we expect a wide range of the rotational activity for mild aspect ratios. We shall find marked orientational and positional heterogeneities on mesoscopic scales in glass. Such heterogeneous patterns have been visualized in various model systems [23–25]. First, there arises a mesoscopic heterogeneity of
042317-1
©2013 American Physical Society
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013)
the flip motions correlated with the impurity clustering. Second, the positional configuration changes cause nonflip rotations of the ellipses, which are the origin of the longtime decay of the rotational correlation functions G (t) of even [10,12–17]. It follows a dynamic heterogeneity of the long-time nonflip rotations correlated with the dynamic heterogeneity in displacements or bond breakage [26–29]. The organization of this paper is as follows. In Sec. II, we will explain our simulation model and method. In Sec. III A, we will present simulation results on the heterogeneities in the orientations and the positions. In Sec. III B, the timecorrelation functions will be examined. In Sec. III C, the angular and translational mean-square displacements will be calculated. In Sec. III D, we will introduce the flip number for each ellipse in a time interval and study its heterogeneity. In Sec. III E, time development of a configurational change with large displacements and/or large angle changes will be illustrated. II. MODEL AND NUMERICAL METHOD
In two dimensions, we consider mixtures of anisotropic and circular particles with numbers N1 and N2 , where N = N1 + N2 = 4096. The concentration of the circular species is c = N2 /N. The particle positions are written as r i (i = 1, . . . ,N). The orientation vectors of the anisotropic particles are expressed as ni = (cos θi , sin θi ) in terms of angles θi (i = 1, . . . ,N1 ). The pair potential Uij between particles i ∈ α and j ∈ β (α,β = 1,2) is a modified Lennard-Jones potential given by [18] 6 12 σαβ σαβ Uij = 4 (1 + Aij ) 12 − (1 + Bij ) 6 , (1) rij rij where rij is the particle distance and is the interaction energy. In terms of characteristic lengths σ1 and σ2 , we set σαβ = (σα + σβ )/2. The potential is truncated at rij = 3σ1 . The particle anisotropy is accounted for by the anisotropic factors Aij and Bij , which depend on the angles between ni , nj , and the relative direction rˆ ij = rij−1 (r i − r j ). In this paper, we set Aij = χ [δα1 (ni · rˆ ij )2 + δβ1 (nj · rˆ ij )2 ],
(2)
Bij = ζ [δα1 δβ2 (ni · rˆ ij )2 + δα2 δβ1 (nj · rˆ ij )2 ],
(3)
where δαβ is the Kronecker δ, χ is the anisotropy strength of repulsion, and ζ is that of attraction between the two species. The Newton equations for r i (t) and θi (t) are given by ∂U d2 ri = − , (4) dt 2 ∂ ri d2 ∂U I 2 θi = − , (5) dt ∂θi where U = i 6, k = 6, and k < 6 and those of the impurities with k > 5, k = 5, and k < 5. These six fractions are denoted by φ1> , φ16 , φ1< , φ2> , φ25 , and φ2< , respectively, as functions of T . At low T , k is mostly 6 or 7 for the ellipses and is mostly 4 or 5 for the impurities. In fact, for the data in Figs. 1 and 2 at T = 0.05, we have (φ16 ,φ1> ) ∼ = (0.66,0.33) and (φ2< ,φ25 ) ∼ = (0.28,0.71) for ζ = 1, while these sets are (0.73,0.26) and (0.05,0.94), respectively, for ζ = 0.5. In the right panel of Fig. 3, φ1< and φ2> are very small at low T and increase with increasing T . Thus, the ellipses with k < 6 and the impurities with k > 5 represent liquidlike defects [35]. Hentschel et al. [35] studied the positional disorder using the Voronoi graphs for a mixture of circular particles with the soft-core potential in two dimensions. In their simulation, small (large) particles enclosed by heptagons (pentagons) form liquidlike defects decreasing at low T . B. Time-correlation functions
For strong short-range anchoring, the rotational dynamics sensitively depends on whether the anisotropic particles are close or far from the impurities. In Fig. 4, we show time
042317-3
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013) 0.8
0.15
0.8
0.6
0.4
0.12
φ5
0.6
0.06 0.2
0
0.03
0
0.2
0.4
0.6
0.8
1
0
ζ=1.0 T=0.2
0.4
0.2 0
0.2
0.4
T
0.6
0.8
1
T φ1> ,
0 -1
φ16 ,
-2/3
-1/3
evolution of the angle changes, θi (t0 ,t + t0 ) = θi (t + t0 ) − θi (t0 ),
(11)
where we pick up a rapidly rotating ellipse, a rarely flipping one, and an inactive one. We can see instantaneous flip motions by ±π . In previous simulations in glassy states, they observed flips for rodlike molecules [10,13,14] and large angle jumps for ortho-terphenyl (OTP) [11,12]. We introduce the distribution of the angle changes, 1 δ {[θi (t)]2π − θ }, N1 i∈1
(12)
where θi (t) = θi (t0 ,t + t0 ) and −π θ < π. For any angle ϕ, we define [ϕ]2π = ϕ − 2pπ , which is in the range [−π,π ] with an integer p. Furthermore, we consider the th
1/3
0
θ/
φ1
6, k = 6, and k < 6, respectively, and those φ2> , φ25 , and φ2< of the impurities with k > 5, k = 5, and k < 5, respectively, as functions of T with ζ = 1, where k is the number of the surrounding triangles in the Delaunay diagrams. Here, φ16 , φ1> , φ2< , and φ25 are large at any T , but φ1< and φ2> decrease at low T . Right: Fractions φ1< and φ2> vs T , which are the fractions of liquidlike defects [35] decreasing at low T .
G(t,θ ) =
t = 40 400 4000
σ2 / σ1 = 0.6
0.09
G(t, θ)
φ6
φ5
2/3
1
π
FIG. 5. (Color online) Time-dependent angle distribution function G(t,θ ) in Eq. (12) at t = 40, 400, and 4000 for ζ = 1 and T = 0.2. Peaks emerge at θ = ±π due to flip motions on the time scale of τ1 = 400. Afterwards, G(t,θ ) → 1/2π on the time scale of τ2 = 24 000.
moments of G(t,θ ) given by π G (t) = dθ G(t,θ ) cos(θ ) −π
1 cos[θi (t)]. = N1 i∈1
(13)
We calculated G(t,θ ), G1 (t), and G2 (t) by taking the average · · · over the initial time t0 and over five runs. In Fig. 5, we show time evolution of G(t,θ ), where flip motions give rise to peaks at θ = ±π growing on the time scale of τ1 . Thus, these flip motions cause the decay of G1 (t) in Fig. 6(a). However, G2 (t) is unchanged by the turnovers and
c=0.2, χ=1.2, σ2 / σ1 = 0.6, ζ=1.0 1
1
(a)
(b) 0.8
0.8
G1(t) 6
0.6
0.6
5
0.4
0.4
G2(t) σ2 / σ1 = 0.6
4
ζ=1.0 T=0.2
Δθ / π
3 2
0.2
T=0.05 0.10 0.15 0.20 0.25 0.30
0.2
T=0.05 0.10 0.15 0.20 0.25 0.30
0 -2 0 -2 -1 0 1 2 3 4 5 -1 0 1 2 3 4 5 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
1
1
t
t
3
6
1
0
(c)
5
(d)
0.8
-1
4
2
-2
0.6
-3 0.4
-4 0
200
400
600
t
0.2
FIG. 4. (Color online) Time evolution of angle changes θi (t0 ,t0 + t) in Eq. (11) for (1) a frequently flipping ellipse, (2) an infrequently flipping one, and (3) an inactive one for ζ = 1 and T = 0.1. Flip events occur at points (◦) on the curves (see the Appendix). These jumps are very different from thermal vibrations but occur thermally.
3
Fs (q,t)
log10(τ)
2
T=0.05 0.10 0.15 0.20 0.25 0.30
τ1 τ2 τα
1 0
-1 0 -2 10 10-1 100 101 102 103 104 105
t
2
4
6
8
10
1/T
FIG. 6. (Color online) (a) G1 (t), (b) G2 (t), (c) Fs (q,t) at q = 2π for ellipses at six T . (d) Relaxation times τ1 , τ2 , and τα in Eqs. (14)– (16) vs 1/T . Here, ζ = 1 and time t is in units of τ0 in Eq. (8).
042317-4
FORMATION OF DOUBLE GLASS IN BINARY MIXTURES . . .
G1 (τ1 ) = 1/e,
(14)
G2 (t) ∝ exp[−(t/τ2 )β ] (t > 1),
(15)
Fs (q,t) ∝ exp[−(t/τα )γ ]
(t > 1),
(16)
where the exponents β and γ are about 0.4 for T 0.2. We here determine τ2 and τα from the long-time relaxations of G2 (t) and Fs (q,t), respectively. In Fig. 6(d), we plot them, where τ1 τα ∼ τ2 . In all the T range in Fig. 6(d), τ1 may be nicely fitted to the Arrhenius form, ln(τ1 ) = 1.4/T − 1.2. On the other hand, τ2 and τα exhibit a changeover at T ∼ 0.2 and can be fitted to the Arrhenius forms as ln(τ2 ) = 2.6/T − 3.0 and ln(τα ) = 2.6/T − 4.5, for T 0.2, so τ2 /τα ∼ = 4 for T 0.2. Thus, G1 (t) decays mainly due to thermally activated flip motions, which are nearly decoupled from the translational motions. On the other hand, G2 (t) and Fs (q,t) decay at longer times due to irreversible configuration changes involving at least several particles. In three dimensions, the distribution of the angles cos−1 [ni (t0 + t) · ni (t0 )] was calculated for OTP [11] and for dumbbells [13,14]. In these papers, this distribution exhibited peaks due to orientational jumps. Also in the rotational time-correlation functions in three dimensions, the Legendre polynomials P [ui (t)] with ui (t) = ni (t0 + t) · ni (t0 ) were used [10,12–17], where G (t) with even decayed slower than those with odd at low T . These previous findings are in accord with our results. C. Mean-square displacements
In the literature, the angular mean-square displacement has been calculated to study the rotational diffusion [12–15,17]. In two dimensions, it is defined by 1 |θi (t0 ,t0 + t)|2 . (17) Mθ (t) = |θ |2 = N1 i∈1 We also introduce the usual positional mean-square displacement for the ellipses by 1 M(t) = |r|2 = |r i (t0 ,t0 + t)|2 , (18) N1 i∈1 where r i (t0 ,t0 + t) = r i (t0 + t) − r i (t0 ). At very short times, these quantities exhibit the ballistic behavior (∝ t 2 ). At long times, they grow linearly in time as (19) Mθ (t) ∼ = 2DR t, M(t) ∼ = 4Dt,
(20)
where DR and D are the rotational and translational diffusion constants, respectively. If the rotational activity is strongly heterogeneous, the so-called Stokes-Einstein-Debye relation DR ∼ kB T /ηa does not hold [3,12], where η is the viscosity and a is the radius of the diffusing particle. In our case, Mθ (t)
ζ=1.0, T=0.2, t f =8000 100
(a) 102 10
10-1
(b)
=1
n > 100 total
pe
slo
10-2
< |Δr | > 2
total 1
104
0
sl o pe =
decays more slowly after the initial relaxation in Fig. 6(b). Notice that G(t,θ ) tends to 1/2π on the time scale of τ2 . In Fig. 6(c), we also show the self-part of the density timecorrelation function Fs (q,t) at q = 2π for the ellipses, which closely resembles G2 (t). We define the relaxation times τ1 , τ2 , and τα as
PHYSICAL REVIEW E 88, 042317 (2013)
10-3 10-2 10
n= 0 n < 10
10-4
-4
10-2
100
t
102
104
10
|Δr | > 0.6
-5
10-2
100
t
102
104
FIG. 7. (Color online) Angular and positional mean-square displacements for ζ = 1 and T = 0.2. (a) Angular one Mθ (t) in Eq. (17) and contributions from those with ni 100, ni 10, and ni = 0, where ni is the flip number for tf = 20τ1 = 8000 (see the Appendix). The contribution from ni 100 approaches Mθ (t) for t 1, leading to DR = 0.14. (b) Positional one M(t) in Eq. (18) and contribution from those with ri > 0.6 in Eq. (21), where the latter grows linearly for t 20 with D = 1.4 × 10−5 .
is greatly increased by rapidly flipping ellipses and DR from it does not correspond to any experimentally observed relaxation rates at low T . In Fig. 7(a), we plot Mθ (t) for ζ = 1 and T = 0.2. Here, it attains the diffusion behavior with DR = 0.14 for t 1, while G1 (t) decays slower with τ1 = 400. We also display the contributions to the sum in Mθ (t) in Eq. (17) from the ellipses with ni = 0, ni 10, and ni 100, where ni is the flip number of ellipse i in a time interval with width tf = 8000 (see the Appendix). The fractions of these three groups are 0.16, 0.34, and 0.39, respectively. Remarkably, the contribution from ni 100 approaches Mθ (t) for t 1, while that from ni 10 behaves diffusively as 0.7 × 10−3 × 2t for t τ1 . Thus, the effective rotational diffusion constant of the ellipses with n 10 is 0.7 × 10−3 /0.34 = 2 × 10−3 . The Mθ (t) itself exhibits the plateau behavior at much lower temperatures (say, T = 0.05), while the contributions from n = 0 and n 10 exhibit it at T = 0.2. In Fig. 7(b), for ζ = 1 and T = 0.2, M(t) still in the course of plateau-to-diffusion crossover even at t = 104 . To obtain small D, we also plot the contribution from the ellipses with large displacements [36], 1 M > (t) = [ri (t) − c ]|ri (t)|2 , (21) N1 i∈1 where ri (t) is an abbreviation of |r i (t0 ,t0 + t)| and (u) is the step function being equal to 1 for u 0 and to 0 for u < 0. The threshold length c is set equal to 0.6. In this restricted sum, the thermal vibrational motions within transient cages are excluded, so it picks up the thermally activated jumps only. As a result, we have the linear growth M > (t) ∼ = 4Dt with D = 1.4 × 10−5 ∼ 10−4 DR from the early stage t 20. This behavior of M > (t) is insensitive to a small change of c [36]. For example, almost the same results followed for c = 0.8. As a similar observation, Chong and Kob found that DR τ2 grows strongly with lowering T for a mixture of rigid dumbbell molecules [17]. See more discussions for other molecular systems in item (i) in the Summary.
042317-5
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013)
We further introduce the n-dependent mean-square displacement among the ellipses with n flips as 1 |θ |2 n (t) = δnni |θi (t0 ,t0 + t)|2 . (23) f N1 φn i∈1
ζ=1.0, T=0.2, t f =8000 1
(b)
(a)
0.8
2
0.6
10nφnf 1
10
-2
2
10 0
1
50
(c)
0.4
< Δθ >n /π 2
2
0.2
φnf
G150 G120
It follows the sum relation Mθ (t) = φnf |θ |2 n (t).
G110
100
(d)
t
102
In our case, most ellipses undergo +π flips and −π flips equally on long times, so i∈1 δnni θi (t0 ,t0 + t) = 0. f f In Fig. 8(a), we plot φn , nφn , and |θ |2 n (tf ) (where t = tf ) for ζ = 1 and T = 0.2. Here, the fraction of the ellipses f with no flip is given by φ0 = 0.16. We find that the flip number distribution is very broad as
104
φnf ∼ 0.1n−1 ,
n < 10
n > 200
D. Distribution of flip numbers
The large size of DR is due to the presence of ellipses frequently undergoing flip motions. In the Appendix, we will give a method of determining the flip number ni for each ellipse i in a time interval [t0 ,t0 + tf ]. We should choose a sufficiently large width tf to detect a wide range of ni . In the following, tf = 8000 at T = 0.2 in Figs. 7 and 8 and tf = 105 at T = 0.05 in Fig. 9. The curves in Figs. 8(a), 8(b), and 9(a) are the averages over six runs. For a given time interval [t0 ,t0 + tf ], the fraction of the ellipses with n flips is written as δnni /N1 . (22) φnf = i∈1
1
ζ=1.0, T=0.05, t f =105
(a)
(b)
φ0f = 0.86
G1(t f )=0.82
0.8 0.6 0.4
10 φnf 2
10
-2
< Δθ >n /πf 2
2
10nφn
0.2 0
50
n
100
(25)
in the range 1 n < nmax , where nmax is an upper bound f about 103 . In the present case, the sums of φn in the ranges 1 n 10, 11 n 99, and n 100 are 0.18, 0.27, and 0.39, respectively. Furthermore, we find
FIG. 8. (Color online) (a) 102 φnf , 10nφnf for n 1, and −2 10 |θ |2 n /π 2 in Eq. (23), where tf = 20τ1 = 8000 (averages over six runs). (b) G1 (t) in Eq. (13) and Gn1 (t) in Eq. (28) with n = 0, 10, 20, and 50 at t = tf , where the latter approach the former at long times. Snapshots of ellipses with ni 10 in (c) and those with ni 200 in (d), whose heterogeneities are correlated with the impurity clustering.
1.2
(24)
n0
G10
0 100 10-2
n
G1
n < 10
FIG. 9. (Color online) 102 φnf , 10nφnf , and 10−2 |θ |2 n /π 2 (averages over six runs) in (a) and snapshot of ellipses with n 10 in f (b), where ζ = 1, T = 0.05, and tf = 105 . At this low T , φ0 = 0.86 and the fraction of the depicted ellipses in (b) is 0.92.
|θ |2 n (t) ∼ π 2 nt = π 2 nt/tf
(26)
for t 1 and n 1. In a general time width t, the ellipses with n flips in the reference time width tf should flip nt = nt/tf times on the average, where t 1 and n 1. Then, together with the sentence below Eq. (24), Eq. (26) is a natural relation. From Eqs. (24)–(26) we find DR ∼ nmax /tf ,
(27)
which means that DR is determined by rapidly rotating ellipses. To be self-consistent, nmax should be proportional to tf ; then, DR is independent of tf . In Fig. 8(b), we compare G1 (t) and the restricted sums, 1 Gn1 (t) = (n − ni ) cos[θi (t0 ,t0 + t)], (28) N1 i∈1 where we set n = 0, 10, 20, and 50. We here pick up the ellipses with flip numbers not exceeding n owing to the step function . We can see that these Gn1 (t) are nearly constant for some time and become nearly equal to G1 (t) after long times. For t 103 , G1 (t) is composed of the contributions from the ellipses with n 10. Note that the ellipses with n flips for tf = 20τ1 should undergo n/20 flips for tf = τ1 on the average. Thus, the ellipses which have undergone more than a few flips in the time interval [t0 ,t0 + τ1 ] do not appreciably contribute to G1 (τ1 ) = 1/e. In Figs. 8(c) and 8(d), we show snapshots of the ellipses with n 10 and n 200, respectively. The distributions of these rotationally inactive and active ellipses are highly heterogeneous. This marked feature is rather natural in view of the mild aspect ratio 1.23 and the significant impurity clustering. In fact, the impurities are nearly absent in the red regions in Fig. 8(b). In Fig. 9, we also show that the flip motions still remain even at T = 0.05. In this case, we find G1 (t) ∼ 0.8 at t = 105 in Fig. 6(a), but we estimate τ1 ∼ 1012 from the extrapolation of the Arrhenius form [see the sentences below Eq. (16)]. f f In Fig. 9(a), we find φ0 = 0.86 and φn ∼ 0.02n−1 and again
042317-6
FORMATION OF DOUBLE GLASS IN BINARY MIXTURES . . .
obtain Eq. (26) for ζ = 1, T = 0.05, and tf = 105 . In Fig. 9(b), displayed is a snapshot of the ellipses with ni 10, whose f fraction is n10 φn = 0.92. Even at this low T , 2% ellipses have ni > 200. We have DR = 1.2 × 10−3 due to these rapidly flipping ellipses. This snapshot was produced by the initial particle configuration common to that in Fig. 8(b). Most of the ellipses in the red regions in Fig. 8(b) are now inactive, since their orientation alignment increases with lowering T . E. Dynamic heterogeneities in nonflip rotations and displacements
Now, we examine the long-time structural relaxation caused by collective configuration changes, where large displacements induce large nonflip rotations. These events should occur repeated in the same fragile regions on time scales longer than τ2 [29,36]. In Fig. 10, we illustrate time development of a configuration change at successive times t0 + t with 0 t 10. Depicted are the ellipses with ri (t) = |r i (t0 ,t0 + t)| > 0.6 or π
θ 0
t=0
ζ =1.0, T=0.2
t=2
PHYSICAL REVIEW E 88, 042317 (2013)
ci (t) < 0.2, where ci (t) = cos[2θi (t0 ,t0 + t)].
The condition ci (t) < 0.2 means 0.22π < |θi (t0 ,t0 + t)| < 0.78π in the range [−π,π ]. From Eq. (13) we have G2 (t) = c i∈1 i (t)/N1 . In the narrow region in Fig. 10, the particle configuration was nearly stationary for t 0, but large particle motions started for t > 0 and continued on a time scale of 10. We can see circulating particle motions at t = 6 and 8 and stringlike ones at t = 10 [27–29,36]. The orientations of these ellipses are largely changing with their movements. For t > 10, the subsequent displacements became small, but considerable orientational motions persisted until t ∼ 20. We examine the nonflip motions in terms of ci (t) in Eq. (29), since it is invariant with respect to turnovers. In Fig. 11(a), we visualize the correlation between ci = ci (t) and ri = ri (t) for ζ = 1 and T = 0.2. We set t = tf = 8000, which is onethird of τ2 ∼ 24 000 (∼4τα ). We present a snapshot of the ellipses with (a) ci < 0.2 and ri > 0.6, (b) ci < 0.2 and ri < 0.6, and (c) ci > 0.2 and ri > 0.6. Here, we exclude the ellipses with ni > 200(∼30%), because they do not exhibit the glassy behavior. The fractions of these depicted groups are (a) 0.08, (b) 0.10, and (c) 0.09, while the fraction of the ellipses with ci > 0.2, ri < 0.6, and ni 200 is 0.41. Thus, if we consider the ellipses with ci < 0.2 and ni < 200, half of them have undergone displacements with ri > 0.6. Also, if we consider the ellipses with ri > 0.6 and ni < 200, half of them have undergone large angle changes with ci < 0.2. In Fig. 11(b), displayed is a snapshot of the ellipses with (a) ri > 0.6 and ni < 200 (18%) and (b) ri > 0.6 and ni > 200 (11%). The group (a) here consists of the groups (a) and (b) in Fig. 11(a). Here, the ellipses in these two groups form clusters, indicating collective displacements. In addition, clusters of one group are adjacent to those of another group. Thus, rotationally active ellipses with large ni tend to be translationally active also. (a)
t=4
ζ=1.0, T=0.2, t=8000
(b)
t=6
n200, |Δr | > 0.6
FIG. 11. (Color online) Comparison of ellipses with ci ≡ cos[2θi (t0 ,t0 + t)] < 0.2 and particles with large displacement ri ≡ |r i (t0 ,t0 + t)| > 0.6 for ζ = 1 and T = 0.2 at t = tf = 8000. In the left, depicted are three groups of ellipses with ci < 0.2 and ri > 0.6 (in red), ci < 0.2 and ri < 0.6 (in green), ci > 0.2 and ri > 0.6 (in yellow), whose fractions are 0.08, 0.10, and 0.09, respectively. In the right, depicted are 0.18N1 ellipses with ri > 0.6 and ni < 200 (in red) and 0.11N1 ones with ri > 0.6 and ni > 200 (in green).
042317-7
KYOHEI TAKAE AND AKIRA ONUKI
PHYSICAL REVIEW E 88, 042317 (2013)
IV. SUMMARY AND REMARKS
We have performed simulation of a mixture of elliptic particle with a mild aspect ratio (=1.23) and smaller circular impurities with σ2 /σ1 = 0.6 at 20%. We have assumed an angle-dependent attractive interaction between the ellipses and the impurities (∝ζ ), which leads to the homeotropic anchoring of the ellipses around each impurity. We summarize our main simulation results. (1) We have shown snapshots of the orientations and the positions in Figs. 1 and 2, which are mesoscopically heterogeneous. From the Delaunay triangulation in the right panels of Fig. 2, we have found that the number of surrounding triangles (the coordination number) is 6 or 7 for the ellipses and 5 or 4 for the impurities in glassy states, as plotted in Fig. 3. A majority of the impurities (∼70%) are surrounded by five ellipses, analogously to the case of the Shintani-Tanaka model [24]. (2) We have calculated the distribution function of the angle changes G(t,θ ) in Eq. (12), which exhibits peaks at θ = ±π for large t due to flip motions as in Fig. 5. We have found that the rotational time-correlation functions G1 (t) and G2 (t) of the ellipses relax very differently at long times in Fig. 6, because G1 (t) decays due to flip motions and G2 (t) due to configuration changes. (3) We have found that the angular mean-square displacement Mθ (t) in Eq. (17) behaves as 2DR t rapidly for t 1 with very large DR in Fig. 7. This is in marked contrast to the slow time evolution of the translational mean-square displacement M(t). However, the contribution to M(t) from the largely displaced ellipses (|ri | > 0.6) has exhibited the diffusion behavior with D = 10−4 DR , because the diffusion is governed by the activation dynamics [36]. f (4) We have displayed the fractions φn of the ellipses with n flips in a time interval [t0 ,t0 + tf ], where tf = 8000 at T = 0.2 in Figs. 7 and 8 and tf = 105 at T = 0.05 in Fig. 9. We have f found a very broad distribution φn (∝ n−1 ) for 1 n < nmax . The angular mean-square displacement |θ |2 n (t) among the ellipses with n flips behaves as π 2 nt/tf . Then DR ∼ nmax /tf due to rapidly flipping ellipses. We have also shown that the long-time decay of G1 (t) is determined by rotationally inactive ellipses. (5) The distributions of the rotationally inactive (n 10) and active (n 200) ellipses at T = 0.2 have been presented in Figs. 8(c) and 8(d), while that of the inactive ones at T = 0.05 has been given in Fig. 9(b). (6) We have illustrated time development of a configuration change in Fig. 10, where large displacements and nonflip rotations are coupled. We have demonstrated close correlation between nonflip rotations and large displacements in Fig. 11. We make remarks as follows. (i) Our potential energy is invariant with respect to turnovers. This means that our model has the quadrupolar symmetry [7]. This is also the case of diatomic molecules or dumbbells [16,17], for which the flip motions should be highly heterogeneous for mild aspect ratios. For methanol, Sindzingre and Klein [10] found flip motions near the glass transition. For OTP, Lewis and Wahnstr¨om [11] found translation-free orientational jumps, while Lombardo et al. [12] found an enhancement in the rotational motions relative to the translation motions at low T . For these systems, the
heterogeneity of orientational jumps should be examined in more detail. However, the physical picture for one-component glass-formers should be different from that for binary mixtures in this paper. (ii) We have chosen a mild aspect ratio (=1.23) to find significant flip motions. However, an increase in the aspect ratio leads to a decrease in the flip frequency, on which we will report shortly. (iii) We have suppressed the clustering of small impurities by the angle-dependent attractive interaction. If we consider large impurities (say, σ2 /σ1 = 1.4), we may realize double glass by adding a repulsive interaction among the impurities suppressing crystal formation. (iv) The phase behavior of mixtures of two species of anisotropic particles should be studied in the future, where we expect nematic or smectic glass. (v) The spatial scales of the structural and orientational heterogeneities depend on various parameters. If the oriented domains are not too small, there arises a large orientationstrain coupling, leading to soft elasticity and a shape-memory effect [18]. These effects were observed for Ti-Ni alloys [37] (where atomic displacements within unit cells cause structural changes). When anisotropic particles have electric dipoles [7], mesoscopic polar domains appear as in ferroelectric glass (relaxors) [38]. Including such metallic alloys also, we point out relevance of the compositional heterogeneity in the glass transition coupled with some orientational phase transition. ACKNOWLEDGMENTS
This work was supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology of Japan. K.T. was supported by the Japan Society for Promotion of Science. The present numerical calculations were carried out on SR16000 at YITP in Kyoto University. APPENDIX: FLIP EVENTS IN NUMERICAL ANALYSIS
In our numerical analysis, we determine a series of flip times, ti1 ,ti2 ,ti3 , . . . for each ellipse i. We write the angle change as θi (t) = θi (t0 ,t + t0 ), suppressing t0 . (i) At the first flip time ti1 we set |θi (ti1 )| = 2π/3.
(A1)
For t > ti1 we introduce θi1 (t) = θi (t) ± π,
(A2)
where +π or −π is chosen such that |θi1 (ti1 + 0)| < π/2. (ii) At the second flip time ti2 we set |θi1 (ti2 )| = 2π/3.
(A3)
For t > ti2 we again introduce θi2 (t) = θi1 (t) ± π,
(A4)
where +π or −π is chosen such that |θi2 (ti2 + 0)| < π/2. (iii) We repeat these procedures to obtain the successive flip times. See Fig. 4 for examples of the flip time series.
042317-8
FORMATION OF DOUBLE GLASS IN BINARY MIXTURES . . .
PHYSICAL REVIEW E 88, 042317 (2013)
Note that the threshold 2π/3 in Eqs. (A1) and (A3) may be changed to another angle, say 5π/6. However, the resultant
flip time series is rather insensitive to this choice as long as it is in the range [π/4,π/2].
[1] C. A. Angell, K. L. Ngai, G. B. McKenna, P. F. McMillan, and S. W. Martin, J. Appl. Phys. 88, 3113 (2000). [2] K. Binder and W. Kob, Glassy Materials and Disordered Solids (World Scientific, Singapore, 2005). [3] F. Fujara, B. Geil, H. Sillescu, and G. Fleischer, Z. Phys. B 88, 195 (1992); M. T. Cicerone and M. D. Ediger, J. Chem. Phys. 104, 7210 (1996). [4] A. Winkler, A. Latz, R. Schilling, and C. Theis, Phys. Rev. E 62, 8004 (2000). [5] S.-H. Chong and W. G¨otze, Phys. Rev. E 65, 041503 (2002). [6] R. Zhang and K. S. Schweizer, J. Chem. Phys. 133, 104902 (2010); 136, 154902 (2012). [7] U. T. H¨ochli, K. Knorr, and A. Loidl, Adv. Phys. 39, 405 (1990). [8] Z. Zheng, F. Wang, and Y. Han, Phys. Rev. Lett. 107, 065702 (2011). [9] C. K. Mishra, A. Rangarajan, and R. Ganapathy, Phys. Rev. Lett. 110, 188301 (2013). [10] P. Sindzingre and M. L. Klein, J. Chem. Phys. 96, 4681 (1992). [11] L. J. Lewis and G. Wahnstr¨om, Phys. Rev. E 50, 3865 (1994); J. Non-Cryst. Solids 172, 69 (1994). [12] T. G. Lombardo, P. G. Debenedetti, and F. H. Stillinger, J. Chem. Phys. 125, 174507 (2006). [13] S. K¨ammerer, W. Kob, and R. Schilling, Phys. Rev. E 56, 5450 (1997). [14] C. De Michele and D. Leporini, Phys. Rev. E 63, 036702 (2001). [15] S.-H. Chong, A. J. Moreno, F. Sciortino, and W. Kob, Phys. Rev. Lett. 94, 215701 (2005). [16] A. J. Moreno, S.-H. Chong, W. Kob, and F. Sciortino, J. Chem. Phys. 123, 204505 (2005). [17] S.-H. Chong and W. Kob, Phys. Rev. Lett. 102, 025702 (2009). [18] K. Takae and A. Onuki, Europhys. Lett. 100, 16006 (2012). [19] J. N. Israelachvili, Intermolecular and Surface Forces (Academic, London, 1991).
[20] C. A. Angell and E. J. Sare, J. Chem. Phys. 49, 4713 (1968); M. Kobayashi and H. Tanaka, J. Phys. Chem. B 115, 14077 (2011). [21] G. A. Vliegenthart, A. van Blaaderen, and H. N. W. Lekkerkerker, Faraday Discuss. 112, 173 (1999). [22] D. Antypov and D. J. Cleaver, J. Chem. Phys. 120, 10307 (2004). [23] T. Hamanaka and A. Onuki, Phys. Rev. E 74, 011506 (2006); 75, 041503 (2007). [24] H. Shintani and H. Tanaka, Nat. Phys. 2, 200 (2006). [25] T. Kawasaki, T. Araki, and H. Tanaka, Phys. Rev. Lett. 99, 215701 (2007). [26] R. Yamamoto and A. Onuki, J. Phys. Soc. Jpn. 66, 2545 (1997); Phys. Rev. E 58, 3515 (1998). [27] C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 80, 2338 (1998). [28] S. C. Glotzer, J. Non-Cryst. Solids 274, 342 (2000). [29] H. Shiba, T. Kawasaki, and A. Onuki, Phys. Rev. E 86, 041504 (2012). [30] S. Nos´e, Mol. Phys. 52, 255 (1984). [31] J. G. Gay and B. J. Berne, J. Chem. Phys. 74, 3316 (1981). [32] A. Ben-Naim, J. Chem. Phys. 54, 3682 (1971); K. A. T. Silverstein, A. D. J. Haymet, and K. A. Dill, J. Am. Chem. Soc. 120, 3166 (1998). [33] J. M. Drouffe, A. C. Maggs, and S. Leibler, Science 254, 1353 (1991); H. Noguchi, J. Chem. Phys. 134, 055101 (2011). [34] D. R. Nelson and B. I. Halperin, Phys. Rev. B 19, 2457 (1979). [35] H. G. E. Hentschel, V. Ilyin, N. Makedonska, I. Procaccia, and N. Schupper, Phys. Rev. E 75, 050404(R) (2007). [36] T. Kawasaki and A. Onuki, Phys. Rev. E 87, 012312 (2013); this paper shows that Dτb is nearly independent of T in a supercooled fragile binary mixture in three dimensions, where τb is the bond breakage time [26]. [37] S. Sarkar, X. Ren, and K. Otsuka, Phys. Rev. Lett. 95, 205702 (2005); Y. Wang, X. Ren, and K. Otsuka, ibid. 97, 225703 (2006). [38] R. A. Cowley, S. N. Gvasaliya, S. G. Lushnikov, B. Roessli, and G. M. Rotaru, Adv. Phys. 60, 229 (2011).
042317-9
