We report an extensive and systematic investigation of the multi-point and multi-time correlation functions to reveal the spatio-temporal structures of dynamic heterogeneities in glass-forming liquids. Molecular dynamics simulations are carried out for the supercooled states of various prototype models of glass-forming liquids such as binary Kob–Andersen, Wahnström, soft-sphere, and network-forming liquids. While the first three models act as fragile liquids exhibiting super-Arrhenius temperature dependence in their relaxation times, the last is a strong glass-former exhibiting Arrhenius behavior. First, we quantify the length scale of the dynamic heterogeneities utilizing the four-point correlation function. The growth of the dynamic length scale with decreasing temperature is characterized by various scaling relations that are analogous to the critical phenomena. We also examine how the growth of the length scale depends upon the model employed. Second, the four-point correlation function is extended to a three-time correlation function to characterize the temporal structures of the dynamic heterogeneities based on our previous studies [K. Kim and S. Saito, Phys. Rev. E79, 060501

R
(2009); 
K.
Kim
and
S.
Saito
, J. Chem. Phys.133, 044511 (2010)]
. We provide comprehensive numerical results obtained from the three-time correlation function for the above models. From these calculations, we examine the time scale of the dynamic heterogeneities and determine the associated lifetime in a consistent and systematic way. Our results indicate that the lifetime of the dynamical heterogeneities becomes much longer than the α-relaxation time determined from a two-point correlation function in fragile liquids. The decoupling between the two time scales is remarkable, particularly in supercooled states, and the time scales differ by more than an order of magnitude in a more fragile liquid. In contrast, the lifetime is shorter than the α-relaxation time in tetrahedral network-forming strong liquid, even at lower temperatures.

Various liquids form disordered and amorphous solids if temperatures are reduced below their melting points while avoiding crystallizations. This transition to a disordered solid is known as the glass transition.1–4 A remarkable feature of supercooled states and glasses is the drastic increase in the viscosity and the structural relaxation time that accompanies non-exponentially observed in various time correlation functions.5–11 Understanding the universal mechanism of the slow dynamics in glass transitions is a challenging problem for condensed matter physics.

To tackle this problem, the notion of “spatially heterogeneous dynamics” or “dynamic heterogeneity” in glass-forming liquids has attracted much attention in recent decades and has been considered central to understanding the slow dynamics of glasses.12–18 Many theoretical, computational, and experimental efforts have been devoted to understanding of dynamic heterogeneities in glassy systems.

Experimentally, various nuclear magnetic resonance (NMR) and other spectroscopic techniques have revealed “heterogeneous” relaxation of the non-exponential decays in glassy systems.19–24 In such systems, the non-exponential relaxation can be explained as the superposition of individual particle contributions with different relaxation rates.13,14

A large number of molecular simulations have also provided information by allowing for the visualization of microscopic details regarding the molecular dynamics.25–47 These simulations have demonstrated direct evidences of the dynamic heterogeneities, i.e., molecular motions accompany correlated domains and, to some extent, exceed the molecular scale in a heterogeneous manner in both time and space. Experiments to directly visualize these dynamic heterogeneities have also been performed in colloidal glasses.48–55 

Those results have required characterizing and quantifying the length and time scales to determine the physical role of dynamic heterogeneities in the underlying mechanism of the glassy slow dynamics. Furthermore, such information regarding spatio-temporal structures would be indispensable to an assessment of the theoretical scenarios and hypotheses proposed thus far, such as the Adam–Gibbs,56 random first-order transition,57–59 dynamic facilitation,15 potential energy landscape,60 mode-coupling,61,62 replicated liquid,63,64 and frustration-limited domain65 approaches. Recent attention has focused on determining the physical origin of dynamic heterogeneities by proposing various length scales including those governed by cooperative rearranging regions,66,67 mosaic length,57,68 bond-orientational order (BOO),69–73 bond-breakage correlations,28,30,47 icosahedral order,74–76 locally preferred structures (LPSs),77,78 geometrical frustration,79–81 patch correlations,82,83 non-local viscoelasticity,84–86 Fickian diffusion,39,87,88 and point-to-set (PTS) correlations.89–95 

Recently, progress has been made in characterizing the dynamic length scale via multi-point correlations. Extractions of the spatial correlations between particle displacements during a typical time interval can lead to the four-point correlation function.36–39,96–105 Theoretical treatments have also been provided to analyze related multi-point susceptibilities based on the mode-coupling approach.104–108 Alternative multi-point correlations104,105,109–111 and nonlinear susceptibilities112–115 have also been proposed for experimental measurements.

In addition, a deeper understanding of the time scale of heterogeneous dynamics has been provided by nonlinear responses such as those studies by multidimensional NMR, hole-burning and photo-bleaching techniques,19–21,116–118 and single molecule measurements.119–124 In these experiments, the lifetime of the dynamic heterogeneity is evaluated from the exchange time between the slow- and fast-moving regions, which is much longer than the structural relaxation time near the glass transition temperature.

In contrast, much less attention has been paid to theoretical and computational explorations of the characteristic time scale of dynamic heterogeneities and its temperature dependence. Note that a characterization of the time scale of dynamic heterogeneities requires an analysis of the duration of the heterogeneous motions, which essentially requires a multi-time extension of the four-point correlation function. Some studies have introduced relevant multi-time correlation functions to investigate the time scale of heterogeneities.31,32,125–137 However, several time intervals are fixed at a characteristic time scale, and thus only limited information is available regarding the time scale of the heterogeneous dynamics.

In our previous studies,138,139 we have emphasized the importance of analyzing a four-point, three-time correlation function for various time intervals to systematically quantify the temporal structures of dynamic heterogeneities and their lifetimes. From the three-time correlation function, we have found that the lifetime increases and becomes much longer than the α-relaxation time τα of the two-point correlation function for a binary soft-sphere supercooled liquid. The observed decoupling between the two time scales at low temperatures is in agreement with the previously mentioned experimental results. The exploited strategy is analogous to the multi-time correlations used in recent multi-dimensional spectroscopy and other related optical techniques.140–148 These spectroscopic techniques have become powerful tools for examining the change from heterogeneous dynamics to homogeneous dynamics in various condensed phase systems.149–156 Furthermore, such multi-dimensional techniques have recently been applied to supercooled and glassy states.157–159 

The aim of the present paper is to investigate the spatio-temporal structures of dynamic heterogeneities by numerically calculating the multi-point and multi-time correlation functions. In particular, we extensively investigate the manner in which the specifics of the model influence the extracted length and time scales. Recently, the effects of the model details on the static correlation and the dynamics have been critically examined for the binary Kob–Andersen model and its Weeks–Chandler–Andersen modification.160–162 While the pair correlations of these models are almost identical, the structural relaxation times differ significantly. These results necessitate an investigation of the model dependence on the length and time scales of the dynamic heterogeneities extracted from multi-point and multi-time correlation functions. Thus, in the present study, we employ a frequently used binary mixture of Kob–Andersen, Wahnström, soft-sphere models, which exhibit a super-Arrhenius temperature dependence in their structural relaxation times, and a network-forming strong liquid model that exhibits Arrhenius behavior. We provide comprehensive numerical results for the four-point correlation function and the three-time correlation function extended from the four-point correlation. Using our extensive numerical results, we systematically determine the length and time scales of dynamic heterogeneities for various glass models.

The paper is organized as follows. In Sec. II, we introduce the simulation details of the glass-forming models used in this paper. In Sec. III, we present numerical calculations of the multi-point and multi-time correlation functions to characterize the spatio-temporal structures of the dynamic heterogeneities. In Sec. III A, we briefly summarize numerical results using conventional two-point correlation functions and determine the α-relaxation time τα. Then, in Sec. III B, we evaluate the growing length scales of the dynamics heterogeneities in the models by analyzing the four-point correlation functions. Finally, in Sec. III C, the lifetimes of the dynamic heterogeneities are quantified from the three-time correlation functions, and their dependence on the model and fragility is examined. In Sec. IV, we summarize our results and provide concluding remarks.

In this study, we carry out extensive molecular dynamics simulations for three-dimensional binary mixtures in the microcanonical ensemble. The system contains N1 particles of component 1 and N2 particles of component 2 under periodic boundary conditions. The total number density is fixed at ρ = N/V with the total number of particles N = N1 + N2 and a system volume V.

The models examined are the well-known prototype models of glass-forming liquids: the binary Kob–Andersen Lennard–Jones (KALJ) liquids,163,164 the binary Wahnström (WAHN) liquids,165 and binary soft-sphere (SS) liquids.166,167 In addition, we also study a model of a network-forming (NTW) liquid that mimics SiO2 with short-range spherical potentials.168 

The binary Lennard–Jones mixture is the most frequently utilized model for the study of glass transitions.163,164 The pair potentials are given by

\begin{equation}v_{\alpha \beta }(r)=4\epsilon _{\alpha \beta }\left[\left(\frac{\sigma _{\alpha \beta }}{r}\right)^{12} -\left(\frac{\sigma _{\alpha \beta }}{r}\right)^{6}\right],\end{equation}
vαβ(r)=4εαβσαβr12σαβr6,
(1)

in which α, β ∈ {1, 2} are the indexes of the particle species. The energy and size ratios are ε1211 = 1.5, ε2211 = 0.5 and σ1211 = 0.8, σ2211 = 0.88, respectively. The masses of the two particle species are equal, m1 = m2 = 1. The interaction is truncated at r = 2.5σαβ. The reduced units σ11, ε11/kB, and

$\sqrt{m_1\sigma _{11}^2/\epsilon _{11}}$
m1σ112/ε11 are used in the model for length, temperature, and time, respectively. The time step is Δt = 0.001 in the reduced time units. The total number density is fixed at ρ = 1.2 with N1 = 800 and N2 = 200. The temperatures investigated are T = 0.7, 0.65, 0.6, 0.55, 0.5, and 0.47.

The KALJ model is a non-additive mixture and thus disobeys the so-called Lorentz–Berthelot combining rules due to the strong attraction between components 1 and 2. Alternatively, the prototype model of the additive and equimolar binary Lennard–Jones mixture is introduced by Wahnström.165,169 The interaction potentials are the same as in Eq. (1), in which the size, mass, and energy ratios are given as σ12 = 1/1.2, m1/m2 = 1/2, and ε12 = 1, respectively. The Lorentz–Berthelot rules,

\begin{eqnarray}\sigma _{\alpha \beta } =\frac{\sigma _\alpha + \sigma _\beta }{2}, \quad \epsilon _{\alpha \beta } =\sqrt{\epsilon _\alpha \epsilon _\beta },\end{eqnarray}
σαβ=σα+σβ2,εαβ=εαεβ,
(2)

are obeyed in this model. Simulation results will be described in terms of the reduced units σ1, ε1/kB, and

$\sqrt{m_1\sigma _{1}^2/\epsilon _{1}}$
m1σ12/ε1 for length, temperature, and time, respectively. The system consists of N1 = 500 and N2 = 500 particles with a fixed number density ρ = 0.75. A time step of Δt = 0.001 is used. The temperatures investigated are T = 0.8, 0.75, 0.7, 0.65, 0.6, and 0.58.

We also study an equimolar binary mixture of soft-sphere particles.166,167 Particles interact via the soft-core potentials

\begin{equation}v_{\alpha \beta }(r)=\epsilon _{\alpha \beta }\left(\frac{\sigma _{\alpha \beta }}{r}\right)^{12},\end{equation}
vαβ(r)=εαβσαβr12,
(3)

with the cubic smoothing function vαβ(r) = B(ar)3 + C for distances

$r > r_c=\sqrt{3}$
r>rc=3⁠. The values of a, B, and C are determined by the continuity conditions up to the second derivative of vαβ(r). The size, mass, and energy ratios are the same as those of the WAHN model: σ12 = 1/1.2, m1/m2 = 1/2, and ε12 = 1, respectively. Thus, this model can be regarded as a purely repulsive interacting system of the WAHN model. Simulation results will be described in terms of the reduced units σ1, ε1/kB, and
$\sqrt{m_1\sigma _{1}^2/\epsilon _{1}}$
m1σ12/ε1
for length, temperature, and time, respectively.

The thermodynamic state of this model is usually characterized by the following non-dimensional parameter,

\begin{equation}\Gamma = \rho \left(\frac{\epsilon _1}{k_BT}\right)^{1/4} {l_0}^3,\end{equation}
Γ=ρε1kBT1/4l03,
(4)

in which l0 represents the effective particle size defined by 4l03 = (2σ1)3 + 2(σ1 + σ2)3 + (2σ2)3. In the simulation, the total number density is given as ρ = l0−3 with N1 = N2 = 500. The investigated states are Γ = 1.30, 1.35, 1.38, 1.42, 1.45, and 1.47. The corresponding temperatures are T = 0.350, 0.301, 0.276, 0.246, 0.226, and 0.214 with a time step of Δt = 0.005.

In addition, we study a model of network-forming liquids interacting via spherical short-ranged potentials.168 This model is simple model and imitates SiO2 glasses, in which tetrahedral networks strongly dominate the dynamics with an Arrhenius behavior for the structural relaxation time, even near the glass transition temperature. The interaction potentials are given as

\begin{equation}v_{\alpha \beta }(r)=\epsilon _{\alpha \beta }\left[\left(\frac{\sigma _{\alpha \beta }}{r}\right)^{12} -(1-\delta _{\alpha \beta })\left(\frac{\sigma _{\alpha \beta }}{r}\right)^{6}\right],\end{equation}
vαβ(r)=εαβσαβr12(1δαβ)σαβr6,
(5)

in which δαβ is the Kronecker delta. The interaction is truncated at r = 2.5σαβ. The size and energy units are determined as follows:

\begin{eqnarray}\sigma _{12} / \sigma _{11} = 0.49,\quad \sigma _{22} / \sigma _{11} = 0.85,\end{eqnarray}
σ12/σ11=0.49,σ22/σ11=0.85,
(6)
\begin{eqnarray}\epsilon _{12} / \epsilon _{11} = 24, \quad \epsilon _{22} / \epsilon _{11} = 1.\end{eqnarray}
ε12/ε11=24,ε22/ε11=1.
(7)

The mass ratio is determined as m2/m1 = 0.57 from the same ratio of O and Si. The units of length, time, and temperature are given as σ11, ε11/kB, and

$\sqrt{m_1\sigma _{11}^2/\epsilon _{11}}$
m1σ112/ε11⁠, respectively. These parameters are adjusted to reproduce the radial distribution functions of the SiO2 amorphous states. Tetrahedral networks are found to be formed due to the highly asymmetric size ratio and the strong attraction between the different components.168 The density of the investigated system is ρ = 1.655, with particle numbers N1 = 1, 000 and N2 = 2, 000. This value corresponds to the density ρ = 2.37 gÅ−3 of the so-called van Beest–Kramer–van Santen (BKS) model for the silica glass.170–172 The simulations are carried out at T = 0.42, 0.4, 0.38, 0.36, 0.34 and 0.32 with a time step of Δt = 0.0005.

First, we study the conventional two-point density correlation function and determine the structural α-relaxation time τα. The self-part of the intermediate scattering function of the component 1 particles

\begin{equation}F_s(k, t)=\left\langle \frac{1}{N_1}\sum _{j=1}^{N_1} \exp [i\bm {k}\cdot \Delta \bm {r}_j(0, t)]\right\rangle ,\end{equation}
Fs(k,t)=1N1j=1N1exp[ik·Δrj(0,t)],
(8)

is calculated for various glass-forming models. Here,

$\Delta \bm {r}_j(0, t)\equiv \bm {r}_j(t)-\bm {r}_j(0)$
Δrj(0,t)rj(t)rj(0) is the jth particle displacement vector at times 0 and t, and k is the wave vector.

The behavior of Fs(k, t) at various temperatures is demonstrated in Fig. 1. The wave vector

$k=|\bm {k}|$
k=|k| is chosen so that the static structure factor of component 1, S11(k), marks the main peak as (a) k = 7.25, (b) 6.65, (c) 6.55, and (d) 8.0 for the KALJ, WAHN, SS, and NTW models, respectively. From these calculations, we determine the α-relaxation time τα as Fs(k, τα) = e−1. In Fig. 2(a), the temperature dependence of τα is plotted as a function of the inverse temperature T0/T. Here T0 denotes the onset temperature introduced in the earlier work.160 We set T0 as T0 = 0.8, 0.7, 0.3, and 0.5 for the KALJ, WAHN, SS, and NTW models, respectively. Below this onset temperature, Fs(k, t) begins to develop a two-step relaxation in each model. Furthermore, the temperature dependence of τα exhibits super-Arrhenius behavior in the KALJ, WAHN, and SS models. This behavior is typical for the fragile glass-forming liquids. In contrast, in the NTW model the tetrahedral networks begin to develop strongly below T0. Correspondingly, the temperature dependence of the α-relaxation time obeys the Arrhenius law.168 In addition, we summarize the power-law behavior as τα ∝ |TTc|−Δ, as predicted by the mode-coupling theory.61,62 The results of the power-law fittings for the four model liquids are demonstrated in Fig. 2(b). We obtain the values of the exponent Δ and the mode-coupling temperature Tc as (Δ, Tc) ≈ (2.4, 0.435), (1.8, 0.56), (2.2, 0.198), and (2.8, 0.31) for KALJ, WAHN, SS, and NTW liquids, respectively. Similar results have been reported for the KALJ, SS, and NTW models.173 We also note that the power-law behavior of the NTW model is reliable over the limited temperature range, as investigated in the simulations for the BKS model.172 

FIG. 1.

Self-part of the intermediate scattering function Fs(k, t) of the component 1 particles for various glass-forming liquid models: (a) KALJ, (b) WAHN, (c) SS, and (d) NTW.

FIG. 1.

Self-part of the intermediate scattering function Fs(k, t) of the component 1 particles for various glass-forming liquid models: (a) KALJ, (b) WAHN, (c) SS, and (d) NTW.

Close modal
FIG. 2.

(a) α-relaxation time τα as a function of the inverse temperature T0/T with the onset temperature T0. (b) α-relaxation time τα as a function of the temperature difference (TTc)/Tc from the mode-coupling divergence temperature Tc. Each dashed line represents the power-law fits with exponent Δ = 2.4, 1.8, 2.2, and 2.8 for the KALJ, WAHN, SS, and NTW models, respectively.

FIG. 2.

(a) α-relaxation time τα as a function of the inverse temperature T0/T with the onset temperature T0. (b) α-relaxation time τα as a function of the temperature difference (TTc)/Tc from the mode-coupling divergence temperature Tc. Each dashed line represents the power-law fits with exponent Δ = 2.4, 1.8, 2.2, and 2.8 for the KALJ, WAHN, SS, and NTW models, respectively.

Close modal

To characterize the growth of the dynamic heterogeneities in supercooled states, the four-point correlation function is introduced to measure the correlation of the particle mobility field at a given time interval. There are several choices for the mobility field such as the particle displacement amplitude,33 the overlap function,38 and the intermediate scattering function.39 Physically, the choice of the mobility field does not alter the fundamental meaning of the four-point correlation function. We here calculate the four-point dynamical susceptibility χ4(k, t), which is defined as the variance of the self-part of the intermediate scattering function Fs(k, t):

\begin{equation}\chi _4(k, t) = N_1 [ \langle \hat{F}_s(\bm {k}, t)^2 \rangle - \langle \hat{F_s}(\bm {k}, t)\rangle ^2 ].\end{equation}
χ4(k,t)=N1[F̂s(k,t)2Fŝ(k,t)2].
(9)

We utilize

$\hat{F}_s(\bm {k}, t)$
F̂s(k,t) expressed as

\begin{equation}\hat{F}_s(\bm {k}, t)= \frac{1}{N_1}\sum _{j=1}^{N_1} \cos [\bm {k}\cdot \Delta \bm {r}_j(0, t)],\end{equation}
F̂s(k,t)=1N1j=1N1cos[k·Δrj(0,t)],
(10)

with

$F_s(k, t)=\langle \hat{F}_s(\bm {k}, t) \rangle$
Fs(k,t)=F̂s(k,t)⁠. The total value of χ4(k, t) is approximately proportional to the extension of the spatial correlations in dynamics with k at a given time interval t because χ4(k, t) investigates the increasing deviation of the two-point correlation function Fs(k, t) from the mean behavior, As is well documented in various studies101–104 and as demonstrated in Fig. 3, χ4(k, t) typically has its maximum value at the time scale of τα, which increases as the temperature decreases. The growth of χ4(k, t) for the strong NTW liquid is suppressed and is smaller than those of other fragile liquids. This behavior implies that the dynamic heterogeneity is less pronounced and plays a minor role in the strong liquid, as revealed in previous studies.104,174,175 Therefore, the dynamics of the strong liquid can be interpreted as mostly occurring through the rearrangement of strongly connecting tetrahedral networks. Interestingly, a similar suppression of χ4(k, t) in strong liquids, for longer times, has been observed in polydispersed systems,73,176 colloidal gels,177,178 and confined systems in random media.179,180

FIG. 3.

Four-point dynamical susceptibility χ4(k, t) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models.

FIG. 3.

Four-point dynamical susceptibility χ4(k, t) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models.

Close modal

The study of the spatial correlation of the four-point correlation function is also essential to extract the growing length scales ξ of the dynamic heterogeneities. In this study, instead of the self-intermediate scattering function, we utilize the frequently used four-point correlation function using the overlap function, which is defined as follows:

\begin{equation}S_4(q, t) = \frac{1}{N_1}\langle Q(\bm {q}, t) Q(-\bm {q}, t)\rangle ,\end{equation}
S4(q,t)=1N1Q(q,t)Q(q,t),
(11)
\begin{equation}Q(\bm {q}, t) =\sum _{j=1}^{N_1} W_j(a, t) \exp [-i \bm {q}\cdot \bm {r}_j(0)],\end{equation}
Q(q,t)=j=1N1Wj(a,t)exp[iq·rj(0)],
(12)

with the wave vector

$q = |\bm {q}|$
q=|q|⁠. Here,
$W_j(a, t) = \Theta (a\break - |\bm {r}_j(t) - \bm {r}_j(0)|)$
Wj(a,t)=Θ(a|rj(t)rj(0)|)
is the overlapping function with the Heaviside step's function Θ(x).36–38 The function Wj(a, t) selects a particle that moves further than the distance a during the time interval t. The value a = 0.3 is typically chosen. As studied in the previous study,38 the relaxation profile of the overlap function
$Q(t)= (1/N_1)\langle \sum _{j=1}^{N_1}W_j(a, t) \rangle$
Q(t)=(1/N1)j=1N1Wj(a,t)
with this probe length scale a = 0.3 approximately corresponds to that of the self-intermediate scattering function Fs(k, t) with the wave number marking the main peak of the static structure factor. We also note that any significant differences in S4(q, t) are not observed if we choose Fs(k, t) as the mobility field. Although recent reports have described the results of S4(q, t) using large-scale simulations,137,181–184 we revisit the identification of the correlation length ξ from the four-point correlation function, Eq. (11). Accordingly, we use much larger systems with N = 100, 000 for the KALJ, WAHN, and SS liquids and N = 90,000 for the NTW liquid to calculate S4(q, t). The other parameters described in Sec. II are unchanged throughout the simulations. Correspondingly, the linear dimension of the system, in termed of the unit length, L = V1/3, is given as L = 43.68, 51.09, 51.27, and 37.89 for the KALJ, WAHN, SS, and NTW models, respectively.

Figure 4 shows the wave-number dependence of S4(q, t) on the time scale τα. As indicated in Fig. 4, S4(q, t) at the α-relaxation time grows in small wave-numbers of q, particularly at low temperatures. This observation indicates that the mobile (or immobile) particles become highly correlated in space when the system undergoes supercooling. The behavior of S4(q, t) at small wave-numbers can be described by the Ornstein–Zernike (OZ) form,

\begin{equation}S_4(q, t) = \frac{\chi (t)}{1 + (q\xi (t))^{\alpha }},\end{equation}
S4(q,t)=χ(t)1+(qξ(t))α,
(13)

in which ξ(t) is regarded as the length scale of the dynamic heterogeneity and χ(t) is the dynamic susceptibility at q → 0. The OZ form with α = 2 has been used in previous studies.30,38,39 However, Fig. 4 shows that the exponent α depends on the details of the simulation model. Figure 5 displays the scaled function S4(q, τα)/χ(τα) as a function of qξ(τα). The results are in good agreement with α ≈ 2.4, 2.0, 2.4, and 1.5 for the KALJ, WAHN, SS, and NTW models, respectively. The apparent difference in the exponent α between the fragile (KALJ, WAHN, and SS) and strong (NTW) glass-forming liquids may be related to the change in geometrical characteristics of the heterogeneous motions if we employ an analogy to the critical phenomena.185 A similar difference in the exponent has been found in kinetically constrained models (KCMs), in which a snapshot of the dynamic heterogeneity in strong KCM model exhibits a smoother cluster structure than in fragile KCM model.41 

FIG. 4.

Four-point static structure factor S4(q, τα) at various temperatures as a function of the wave-number q for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models.

FIG. 4.

Four-point static structure factor S4(q, τα) at various temperatures as a function of the wave-number q for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models.

Close modal
FIG. 5.

Scaled four-point static structure factor S4(q, τα)/χ(τα) as a function of qξ(τα) at various temperatures for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The dashed line represents the Ornstein–Zernike form 1/(1 + (qξ(τα))α) with (a) α = 2.4, (b) 2.0, (c) 2.4, and (d) 1.5.

FIG. 5.

Scaled four-point static structure factor S4(q, τα)/χ(τα) as a function of qξ(τα) at various temperatures for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The dashed line represents the Ornstein–Zernike form 1/(1 + (qξ(τα))α) with (a) α = 2.4, (b) 2.0, (c) 2.4, and (d) 1.5.

Close modal

The values of the length scales ξ and the susceptibility χ at the α-relaxation time are determined from the fitting to Eq. (13). Figure 6(a) shows the temperature dependence of the correlation length ξ(τα). We find that the increasing length scale ξ(τα) with decreasing temperature can be described by the mode-coupling-like power-law ξ(τα) ∼ |TTc|−ν at the investigated temperatures. The exponent is ν ≈ 0.5 for the fragile KALJ, WAHN, and SS models, whereas ν decreases to ν ≈ 0.25 for the strong NTW model. However, the length ξ at lowest temperature T = 0.32 is apparently deviated from the power-law behavior of the NTW model, which is a strong liquid exhibiting the Arrhenius behavior. This limited power-law behavior is also observed in the temperature dependence of the α-relaxation time (see Fig. 2(b)). In addition, we obtained the scaling relationship τα ∼ ξ(τα)z found in Fig. 6(b). In analogy to the dynamical critical phenomena, the exponent z is approximately determined by the relationship z = Δ/ν, i.e., the exponent z becomes z ≈ 4.4, 4.8, 3.6, and 11.0 for the KALJ, WAHN, SS, and NTW models, respectively.

FIG. 6.

(a) Correlation length ξ(τα) as a function of (TTc)/Tc. The dashed line represents the power-law behavior, ξ(τα) ∼ |TTc|−ν with ν = 0.5 (red) and ν = 0.25 (orange), respectively. (b) The relationship between the correlation length ξ(τα) and the α-relaxation time τα. The dashed line represents the power-law behavior, τα ∼ ξ(τα)z with z = 4.4 (red) and z = 11.0 (orange), respectively. (c) The dynamic susceptibility χ(τα) as a function of (TTc)/Tc. The dashed line represents the power-law behavior, χ(τα) ∼ |TTc|−γ with ν = 1.2 (red) and ν = 0.37 (orange), respectively. (d) The relationship between the correlation length ξ(τα) and the dynamic susceptibility χ(τα). The dashed line represents the power-law behavior, χ(τα) ∼ ξ(τα)2−η with 2 − η = 2.4 (red) and 2 − η = 1.5 (orange), respectively.

FIG. 6.

(a) Correlation length ξ(τα) as a function of (TTc)/Tc. The dashed line represents the power-law behavior, ξ(τα) ∼ |TTc|−ν with ν = 0.5 (red) and ν = 0.25 (orange), respectively. (b) The relationship between the correlation length ξ(τα) and the α-relaxation time τα. The dashed line represents the power-law behavior, τα ∼ ξ(τα)z with z = 4.4 (red) and z = 11.0 (orange), respectively. (c) The dynamic susceptibility χ(τα) as a function of (TTc)/Tc. The dashed line represents the power-law behavior, χ(τα) ∼ |TTc|−γ with ν = 1.2 (red) and ν = 0.37 (orange), respectively. (d) The relationship between the correlation length ξ(τα) and the dynamic susceptibility χ(τα). The dashed line represents the power-law behavior, χ(τα) ∼ ξ(τα)2−η with 2 − η = 2.4 (red) and 2 − η = 1.5 (orange), respectively.

Close modal

We also examine the relationship between the susceptibility χ(τα) and the correlation length ξ(τα) via the scaling χ(τα) ∼ ξ(τα)2−η, as demonstrated in Figure 6(d). These results indicate that the scaling exponent of each model is correlated to the value of the exponent α utilized in the OZ function, as shown in Eq. (13). In addition, the dynamic susceptibility χ(τα) is approximated by the scaling relation χ(τα) ∼ |TTc|−γ with γ = (2 − η)ν, as shown in Fig. 6(c). Similar dynamic criticality has been investigated in Refs. 41 and 100.

Here, we remark that a recent theoretical treatment has been presented based on mode-coupling approach, referred to as the inhomogeneous mode-coupling theory (IMCT).107,108 The IMCT predicts the scaling exponents to be ν = 1/4 and 2 − η = 4 at the α-relaxation time scale. To examine the validity of these predictions, intensive large-scale molecular simulations have been carried out.181–184 Our exponents, z ≈ 4.4 and 2 − η ≈ 2.4, of the KALJ model are similar to the values obtained in Refs. 104,182. In these findings, the molecular simulations disagree with the theoretical predictions. However, the IMCT does not investigate a four-point correlation function such as Eq. (11) but treats the three-point correlation defined as the response function of the two-point correlation function due to the inhomogeneous external field. Further analysis is needed to understand the relationship between the four-point correlations utilized here and the inhomogeneous three-point susceptibilities.

In this subsection, we explore and accentuate the characteristic time scale of the dynamic heterogeneity. To this end, we search for the multi-time extension of the four-point correlation function χ4(k, t) by introducing

\begin{equation}\Delta F_4(k, t_1, t_2, t_3) = \biggl \langle \frac{1}{N_1}\sum _{j=1}^{N_1} \delta F_j(\bm {k}, \tau _2, \tau _3) \delta F_j(\bm {k}, 0, \tau _1))\biggr \rangle ,\end{equation}
ΔF4(k,t1,t2,t3)=1N1j=1N1δFj(k,τ2,τ3)δFj(k,0,τ1)),
(14)

in which

\begin{equation}\delta F_j(\bm {k}, 0, t) = \cos [\bm {k}\cdot \Delta \bm {r}_j(0, t)] - F_s(k, t)\end{equation}
δFj(k,0,t)=cos[k·Δrj(0,t)]Fs(k,t)
(15)

provides the individual fluctuations in the two-point correlation function at times 0 and t. This three-time correlation function examines the correlations at four different times, 0, τ1, τ2, and τ3. In practice, ΔF4(k, t1, t2, t3) reveals the correlation of the two-point correlation function Fs(k, t) between the two time intervals, t1 = τ1 and t3 = τ3 − τ2. Furthermore, the progressive changes in the waiting time t2 = τ2 − τ1 of the three-time correlation function ΔF4(k, t1, t2, t3) allow for an investigation of how the correlated motions change with time t2. In fact, the first time-interval portion

$\delta F_j(\bm {k}, 0, \tau _1)$
δFj(k,0,τ1) enables a distinction between the subensemble of slow- and fast-moving particles in the dynamic heterogeneities. In addition, the total function
$\delta F_j(\bm {k}, 0, \tau _1)\delta F_j(\bm {k}, \tau _2, \tau _3)$
δFj(k,0,τ1)δFj(k,τ2,τ3)
can reveal how long the correlations in the dynamics of the subensembles remain over the waiting time t2.138,139 Therefore, this multi-time correlation function can provide the temporal structures and the associated characteristic time scale of the dynamic heterogeneity. In other words, the time scale extracted from the three-time correlation function can be regarded as the lifetime of the dynamic heterogeneity, which should be associated with the length scale ξ determined from the four-point correlation function.

Figure 7 shows the time evolutions of the three-time correlation functions ΔF4 at the lowest temperature for each glass-forming liquid. Note that the wave-number k is chosen as the same value used in the calculation of Fs(k, t) in Fig. 1. We also provide the diagonal portions of the time evolutions along t = t1 = t3 for various t2 in Fig. 8. In previous studies,138,139 we have examined the three-time correlation functions for the SS model. This work confirms that the basic features of ΔF4 are similar and that the time evolutions at t2 occur similarly among all of the glass models studied: (i) As the temperature decreases, the intensity of the three-time correlation ΔF4(k, t1, t2, t3) gradually increases (see the results for other temperatures in the supplementary material186). This result is due to the correlations between particles that move slower (faster) during the first time interval t1 and remain slower (faster) during the second time interval t3. (ii) The correlations of ΔF4(k, t1, t2, t3) at t2 = 0 between the first time interval t1 and the subsequent time interval t3 are noticeable over wide time scales. This result implies that the particle motions are coupled not only at the (diagonal) α- and α-relaxation time scales, but also at the (off-diagonal) α- and β-relaxation time scales. (iii) With increasing the waiting time t2, ΔF4 gradually decays to zero, indicating that the dynamics change from heterogeneous to homogeneous because of the memory loss in the correlated motions between the two time intervals t1 and t3 for the given waiting time t2. However, it should be noted that the line shape of the NTW model, which exhibits the strong (Arrhenius) glass behavior, differs from those of the KALJ, WAHN, and SS models, which exhibit fragile glass (super-Arrhenius) behavior, particularly at t2 = 0. We observe the strong correlations of ΔF4 at small values of t1 and t3 in the NTW model, which approximately corresponds to the time scale of the early β-relaxation at which Fs(k, t) exhibits damped oscillations on a plateau as shown in Fig. 1(d). This well-known finite-size effect in silica glasses168,171 would be smaller if simulations were carried out for larger systems. Furthermore, the relaxation time scale of ΔF4(k, t1, t2, t3) depends on the model. As demonstrated in Fig. 7(d), ΔF4 of the NTW model decays rapidly. This time scale is clearly smaller than the α-relaxation time determined by the two-point correlation function. In contrast, the relaxation of ΔF4 occurs on a time scale comparable to (or exceeding) τα in the fragile KALJ, WAHN, and SS models, as depicted in Figs. 7(a)–7(c).

FIG. 7.

Two-dimensional correlation maps of the three-time correlation functions ΔF4(k, t1, t2, t3) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The state is chosen at the lowest temperature for each model as (a) T = 0.47, (b) T = 0.58, (c) Γ = 1.47, and (d) T = 0.32. The waiting times t2 normalized by τα are increased from left to right in each model. The vertical (horizontal) dashed line denotes the α-relaxation time as t1 = τα (t3 = τα). Note that the waiting time is different in each figure.

FIG. 7.

Two-dimensional correlation maps of the three-time correlation functions ΔF4(k, t1, t2, t3) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The state is chosen at the lowest temperature for each model as (a) T = 0.47, (b) T = 0.58, (c) Γ = 1.47, and (d) T = 0.32. The waiting times t2 normalized by τα are increased from left to right in each model. The vertical (horizontal) dashed line denotes the α-relaxation time as t1 = τα (t3 = τα). Note that the waiting time is different in each figure.

Close modal
FIG. 8.

The diagonal portions of the three-time correlation functions ΔF4(k, t, t2, t) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The state is chosen at the lowest temperature for each model as (a) T = 0.47, (b) T = 0.58, (c) Γ = 1.47, and (d) T = 0.32. The waiting times t2 are normalized by τα at each temperature.

FIG. 8.

The diagonal portions of the three-time correlation functions ΔF4(k, t, t2, t) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The state is chosen at the lowest temperature for each model as (a) T = 0.47, (b) T = 0.58, (c) Γ = 1.47, and (d) T = 0.32. The waiting times t2 are normalized by τα at each temperature.

Close modal

To quantitatively distinguish between the models we quantify the average lifetime of the dynamic heterogeneities using the waiting-time t2 dependence of ΔF4. To this end, we define the volume of the heterogeneities by integrating over the two time intervals t1 and t3,138,139

\begin{equation}\Delta _{\rm hetero}(t_2) = \int _0^{\infty }dt_3\int _0^{\infty }dt_1 \Delta F_4(k, t_1, t_2, t_3).\end{equation}
Δ hetero (t2)=0dt30dt1ΔF4(k,t1,t2,t3).
(16)

This integration resembles the underlying strategy of nonlinear responses such as NMR, hole-burning, and photo-bleaching techniques.13,19,21,117,118 In simulations, similar procedures have been utilized to analyze relevant multi-time correlations.133,134 Figure 9 illustrates the t2 dependence of Δhetero(t2) normalized by Δhetero(0) for the KALJ, WAHN, SS, and NTW models. The waiting time t2 is normalized by τα at each temperature. Figure 9(d), clearly shows that Δhetero of the NTW model decays rapidly and that the relaxation time is smaller than τα at any temperature. In contrast, Figs. 9(a)–9(c) demonstrate that in the fragile liquid models, the relaxation of Δhetero is slower than τα with decreasing temperature. Remarkably, the characteristic time scale of Δhetero for the WAHN model at the lowest temperature exceeds τα by more than an order of magnitude.

FIG. 9.

Waiting time t2 dependence of the integrated three-time correlation function Δhetero(k, t2)/Δhetero(k, 0) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The waiting times are normalized by τα for each temperature. The solid curve is determined by a fitting with the stretched-exponential form exp [−(t2hetero)c] with c ≈ 0.6, 0.5, 0.5, and 0.3 for the KALJ, WAHN, SS, and NTW models, respectively.

FIG. 9.

Waiting time t2 dependence of the integrated three-time correlation function Δhetero(k, t2)/Δhetero(k, 0) for the (a) KALJ, (b) WAHN, (c) SS, and (d) NTW models. The waiting times are normalized by τα for each temperature. The solid curve is determined by a fitting with the stretched-exponential form exp [−(t2hetero)c] with c ≈ 0.6, 0.5, 0.5, and 0.3 for the KALJ, WAHN, SS, and NTW models, respectively.

Close modal

The dependence of the normalized volume Δhetero(t2)/Δhetero(0) on the waiting time t2 can be fitted by the stretched-exponential function exp [−(t2hetero)c], as demonstrated in Fig. 9. The exponent c is approximately c ≈ 0.6, 0.5, 0.5, and 0.3 for the KALJ, WAHN, SS, and NTW models, respectively. From this analysis, we determine the average lifetime of the dynamic heterogeneity as τhetero at various temperatures for each model. We plot τhetero as a function of the inverse temperature T0/T in Fig. 10(a). As shown in Fig. 10(b), we obtain a relationship between the two time scales τhetero and τα that follows the power-law-like behavior, τhetero ∼ ταζ with ζ ≈ 1.25, 1.9, 1.5, and 0.9 for the KALJ, WAHN, SS, and NTW models, respectively. From the analysis, we find that τhetero of the network-forming strong liquid (NTW) is not greater than τα and tends to decrease with decreasing temperature T. That relationship is responsible for the minor role of the dynamic heterogeneities in strong liquids, as previously discussed when we examined the four-point correlations and the associated length scale ξ in Sec. III B. In contrast, the ratio τheteroα in the fragile liquid models exhibits the opposite temperature dependence, i.e., the lifetime τhetero exceeds the α-relaxation time τα with decreasing temperature. However, the increase in τhetero upon supercooling is considerably different among the fragile KALJ, WAHN, and SS models. The ratio τheteroα markedly exceeds the unity at lower temperatures in the WAHN and SS models, whereas τhetero in the KALJ model remains on the time scale of τα, even at the lowest temperature. We note that this feature of the KALJ model was also found in a previous study using a similar multi-time correlation function, in which its lifetime is comparable to τα.131 

A discussion of the significant observed dependence of the lifetime τhetero on the model details is meaningful, even for the fragile models. Recently, the fragility indexes of the simulated models such as the KALJ and WAHN models have been critically investigated from the perspective of the many-body static correlations hidden in the two-point correlations such as the usual radial distribution function. It has been found that the slow- and long-lived correlated domains correspond to the LPSs are characterized by the Voronoi polyhedra.77,78 These studies have also revealed that the non-additive KALJ mixture is less fragile than the additive WAHN mixture. This difference in fragility can be explained in terms of the spatial extent of the LPS domains. In fact, the growth of the LPS domains is significant in the WAHN model that develops icosahedral order upon supercooling. In contrast, the LPS domains in the KALJ model formed by a bicapped prismatic order are found to be smaller than that in the WAHN model. Given these findings, one can conclude that the overall temperature tendency of the lifetime τhetero shown in Fig. 10 is correlated and intensely sensitive to the fragility of the model, i.e., more fragile liquids tend to exhibit longer dynamic heterogeneity lifetimes τhetero.

FIG. 10.

(a) Average lifetime of DH τhetero normalized by the α-relaxation τα as a function of the inverse temperature T0/T with the onset temperature T0. The dashed line represents τhetero = τα as a viewing guide. (b) The relationship between two time scales, τhetero and τα. The dashed line is the power-law behavior τhetero ∼ ταζ with a slope of ζ ≈ 1.9, 1.5, 1.25, and 0.9 from top to bottom.

FIG. 10.

(a) Average lifetime of DH τhetero normalized by the α-relaxation τα as a function of the inverse temperature T0/T with the onset temperature T0. The dashed line represents τhetero = τα as a viewing guide. (b) The relationship between two time scales, τhetero and τα. The dashed line is the power-law behavior τhetero ∼ ταζ with a slope of ζ ≈ 1.9, 1.5, 1.25, and 0.9 from top to bottom.

Close modal

We have examined the four-point correlation function and its three-time correlation extension to systematically characterize the length and time scales of dynamic heterogeneities in prototype fragile (KALJ, WAHN, and SS) and strong (NTW) glass models. Analyses such as the extensive investigation performed herein to determine not only the length scale but also the time scale of various glass-forming models have not, to the best of our knowledge, been previously reported.

First, we quantified the growing length scale of the dynamic heterogeneities upon supercooling as determined by the wave-number dependence of the four-point correlation function using the overlap function. The scaling relationships of the extracted length scale ξ, which are analogous to the dynamical scaling obtained for critical phenomena, were consistently explored for the employed glass models. We observed that the length scale increases with decreasing temperature depending on the fragility of glass. In particular, the increase in the dynamic length scale of the strong liquid is suppressed compared with those of the fragile glass-forming liquids, indicating that the dynamic heterogeneity is less pronounced and plays only a minor role in the strong liquid. We also commented on the comparisons of our numerical results with the theoretical predictions of the IMCT.

Second, we investigated the time scale of the dynamic heterogeneities from determining how long the heterogeneous dynamics survive. Comprehensive numerical results of the three-time correlation function were demonstrated via two-dimensional correlation maps with an analogy to the multi-dimensional spectroscopic methods, as outlined in the Introduction. From the progressive changes in the second time interval of the three-time correlation function, we quantified the characteristic time scale of the dynamic heterogeneities and the associated lifetime τhetero. The lifetime τhetero exceeds the α-relaxation time τα, particularly for highly supercooled states in fragile glass-forming liquids. In contrast, τhetero is smaller than τα even at lower temperatures in the strong liquid, indicating that the dynamic heterogeneities play a minor role. Furthermore, we observed that the temperature dependence of τhetero depends significantly on the fragility, i.e., more fragile liquids exhibit long-lived dynamic heterogeneities with a time scale that exceeds the α-relaxation time. The two time scales differ by more than an order of magnitude in the WAHN model.

Finally, we remark that it is of important to investigate the relationship between the length and time scales of dynamic heterogeneities and their model dependence. In fact, we find that the length scale ξ of fragile liquids increases with decreasing temperature in a similar manner, as seen in Fig. 6(a). In contrast, the time scale τhetero is more sensitive to the fragility and becomes noticeably longer than τα as the fragility index increases. To clarify it, further work that utilizes not only multi-point and multi-time correlations but also other measurements including configurational entropy, LPS, PTS, and BOO is required.

The authors thank Ryoichi Yamamoto, Hideyuki Mizuno, Kunimasa Miyazaki, Takeshi Kawasaki, Hayato Shiba, and Atsushi Ikeda for helpful comments. We also gratefully acknowledge the information on the NTW model from Daniele Coslovich. K.K. is grateful to David Reichman and Glen Hocky for valuable discussions. This work was partially supported by Grant-in-Aid for Young Scientists (A) (Grant No. 23684037) (K.K.) and Scientific Research (B) (Grant No. 22350013) (S.S.) from Japan Society for the Promotion of Science (JSPS). This work was also supported by the Strategic Programs for Innovative Research (SPIRE) and the Computational Materials Science Initiative (CMSI) of the Ministry of Education, Culture, Sports, Science and Technology, Japan (MEXT). The computations were performed at Research Center for Computational Science, Okazaki Research Facilities, National Institutes of Natural Sciences, Japan.

1.
P. G.
Debenedetti
,
Metastable Liquids
(
Princeton University
,
Princeton, NJ
,
1996
).
2.
E.
Donth
,
The Glass Transition
(
Springer
,
Berlin
,
2001
).
3.
K.
Binder
and
W.
Kob
,
Glassy Materials and Disordered Solids
(
World Scientific
,
Singapore
,
2005
).
4.
P. G.
Wolynes
and
V.
Lubchenko
,
Structural Glasses and Supercooled Liquids
(
Wiley
,
USA
,
2012
).
5.
M. D.
Ediger
,
C. A.
Angell
, and
S. R.
Nagel
,
J. Phys. Chem.
100
,
13200
(
1996
).
6.
J.
Jackle
,
Rep. Prog. Phys.
49
,
171
(
1999
).
7.
K. L.
Ngai
,
J. Non-Cryst. Solids
275
,
7
(
2000
).
8.
C. A.
Angell
,
K. L.
Ngai
,
G. B.
McKenna
,
P. F.
McMillan
, and
S. W.
Martin
,
J. Appl. Phys.
88
,
3113
(
2000
).
9.
P. G.
Debenedetti
and
F. H.
Stillinger
,
Nature (London)
410
,
259
(
2001
).
11.
M. D.
Ediger
and
P.
Harrowell
,
J. Chem. Phys.
137
,
080901
(
2012
).
12.
H.
Sillescu
,
J. Non-Cryst. Solids
243
,
81
(
1999
).
13.
M. D.
Ediger
,
Annu. Rev. Phys. Chem.
51
,
99
(
2000
).
14.
R.
Richert
,
J. Phys.: Condens. Matter
14
,
R703
(
2002
).
15.
D.
Chandler
and
J. P.
Garrahan
,
Annu. Rev. Phys. Chem.
61
,
191
(
2010
).
16.
17.
Dynamical Heterogeneities in Glasses, Colloids, and Granular Media
, edited by
L.
Berthier
,
G.
Biroli
,
J.-P.
Bouchaud
,
L.
Cipelletti
, and
W.
van Saarloos
(
Oxford University Press
,
USA
,
2011
).
18.
L.
Berthier
and
G.
Biroli
,
Rev. Mod. Phys.
83
,
587
(
2011
).
19.
K.
Schmidt-Rohr
and
H. W.
Spiess
,
Phys. Rev. Lett.
66
,
3020
(
1991
).
20.
A.
Heuer
,
M.
Wilhelm
,
H.
Zimmermann
, and
H. W.
Spiess
,
Phys. Rev. Lett.
75
,
2851
(
1995
).
21.
R.
Böhmer
,
G.
Hinze
,
G.
Diezemann
,
B.
Geil
, and
H.
Sillescu
,
Europhys. Lett.
36
,
55
(
1996
).
22.
R.
Böhmer
,
R. V.
Chamberlin
,
G.
Diezemann
,
B.
Geil
,
A.
Heuer
,
G.
Hinze
,
S. C.
Kuebler
,
R.
Richert
,
B.
Schiener
,
H.
Sillescu
,
H. W.
Spiess
,
U.
Tracht
, and
M.
Wilhelm
,
J. Non-Cryst. Solids
235–237
,
1
(
1998
).
23.
U.
Tracht
,
M.
Wilhelm
,
A.
Heuer
,
H.
Feng
,
K. S.
Rohr
, and
H. W.
Spiess
,
Phys. Rev. Lett.
81
,
2727
(
1998
).
24.
E.
Vidal Russell
and
N. E.
Israeloff
,
Nature (London)
408
,
695
(
2000
).
25.
M. M.
Hurley
and
P.
Harrowell
,
Phys. Rev. E
52
,
1694
(
1995
)
26.
T.
Muranaka
and
Y.
Hiwatari
,
Phys. Rev. E
51
,
R2735
(
1995
).
27.
W.
Kob
,
C.
Donati
,
S. J.
Plimpton
,
P. H.
Poole
, and
S. C.
Glotzer
,
Phys. Rev. Lett.
79
,
2827
(
1997
).
28.
R.
Yamamoto
and
A.
Onuki
,
J. Phys. Soc. Jpn.
66
,
2545
(
1997
).
29.
C.
Donati
,
J. F.
Douglas
,
W.
Kob
,
S. J.
Plimpton
,
P. H.
Poole
, and
S. C.
Glotzer
,
Phys. Rev. Lett.
80
,
2338
(
1998
).
30.
R.
Yamamoto
and
A.
Onuki
,
Phys. Rev. E
58
,
3515
(
1998
).
31.
R.
Yamamoto
and
A.
Onuki
,
Phys. Rev. Lett.
81
,
4915
(
1998
).
32.
D. N.
Perera
and
P.
Harrowell
,
J. Chem. Phys.
111
,
5441
(
1999
).
33.
C.
Donati
,
S. C.
Glotzer
, and
P. H.
Poole
,
Phys. Rev. Lett.
82
,
5064
(
1999
).
34.
K.
Kim
and
R.
Yamamoto
,
Phys. Rev. E
61
,
R41
(
2000
).
35.
B.
Doliwa
and
A.
Heuer
,
Phys. Rev. E
61
,
6898
(
2000
).
36.
S. C.
Glotzer
,
J. Non-Cryst. Solids
274
,
342
(
2000
).
37.
N.
Lačević
,
F. W.
Starr
,
T. B.
Schrøder
,
V. N.
Novikov
, and
S. C.
Glotzer
,
Phys. Rev. E
66
,
030101
(
2002
).
38.
N.
Lačević
,
F. W.
Starr
,
T. B.
Schrøder
, and
S. C.
Glotzer
,
J. Chem. Phys.
119
,
7372
(
2003
).
39.
L.
Berthier
,
Phys. Rev. E
69
,
020201
R
(
2004
).
40.
L.
Berthier
and
J. P.
Garrahan
,
J. Phys. Chem. B
109
,
3578
(
2005
).
41.
A. C.
Pan
,
J. P.
Garrahan
, and
D.
Chandler
,
Phys. Rev. E
72
,
041106
(
2005
).
42.
A.
Widmer-Cooper
and
P.
Harrowell
,
Phys. Rev. Lett.
96
,
185701
(
2006
).
43.
C.
Brito
and
M.
Wyart
,
J. Stat. Mech.
2007
,
L08003
.
44.
A.
Widmer-Cooper
,
H.
Perry
,
P.
Harrowell
, and
D. R.
Reichman
,
Nat. Phys.
4
,
711
(
2008
).
45.
R.
Candelier
,
A. W.
Cooper
,
J. K.
Kummerfeld
,
O.
Dauchot
,
G.
Biroli
,
P.
Harrowell
, and
D. R.
Reichman
,
Phys. Rev. Lett.
105
,
135702
(
2010
).
46.
A. S.
Keys
,
L. O.
Hedges
,
J. P.
Garrahan
,
S. C.
Glotzer
, and
D.
Chandler
,
Phys. Rev. X
1
,
021013
(
2011
).
47.
H.
Shiba
,
T.
Kawasaki
, and
A.
Onuki
,
Phys. Rev. E
86
,
041504
(
2012
).
48.
A. H.
Marcus
,
J.
Schofield
, and
S. A.
Rice
,
Physs Rev. E
60
,
5725
(
1999
).
49.
W. K.
Kegel
and
A.
van Blaaderen
,
Science
287
,
290
(
2000
).
50.
E. R.
Weeks
,
J.
Crocker
,
A. C.
Levitt
,
A.
Schofield
, and
D.
Weitz
,
Science
287
,
627
(
2000
).
51.
E. R.
Weeks
and
D. A.
Weitz
,
Phys. Rev. Lett.
89
,
095704
(
2002
).
52.
E. R.
Weeks
,
J. C.
Crocker
, and
D. A.
Weitz
,
J. Phys.: Condens. Matter
19
,
205131
(
2007
).
53.
V.
Prasad
,
D.
Semwogerere
, and
E. R.
Weeks
,
J. Phys.: Condens. Matter
19
,
113102
(
2007
).
54.
T.
Narumi
,
S. V.
Franklin
,
K. W.
Desmond
,
M.
Tokuyama
, and
E. R.
Weeks
,
Soft Matter
7
,
1472
(
2011
).
55.
G. L.
Hunter
and
E. R.
Weeks
,
Rep. Prog. Phys.
75
,
066501
(
2012
).
56.
G.
Adam
and
J. H.
Gibbs
,
J. Chem. Phys.
43
,
139
(
1965
).
57.
T. R.
Kirkpatrick
,
D.
Thirumalai
, and
P. G.
Wolynes
,
Phys. Rev. A
40
,
1045
(
1989
).
58.
J. P.
Bouchaud
and
G.
Biroli
,
J. Chem. Phys.
121
,
7347
(
2004
).
59.
V.
Lubchenko
and
P. G.
Wolynes
,
Annu. Rev. Phys. Chem.
58
,
235
(
2007
).
60.
A.
Heuer
,
J. Phys.: Condens. Matter
20
,
373101
(
2008
).
61.
D. R.
Reichman
and
P.
Charbonneau
,
J. Stat. Mech.
2005
,
P05013
.
62.
W.
Götze
,
Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory
(
Oxford University Press
,
USA
,
2009
).
63.
M.
Mézard
and
G.
Parisi
,
J. Chem. Phys.
111
,
1076
(
1999
).
64.
G.
Parisi
and
F.
Zamponi
,
Rev. Mod. Phys.
82
,
789
(
2010
).
65.
G.
Tarjus
,
S. A.
Kivelson
,
Z.
Nussinov
, and
P.
Viot
,
J. Phys.: Condens. Matter
17
,
R1143
(
2005
).
66.
S.
Karmakar
,
C.
Dasgupta
, and
S.
Sastry
,
Proc. Natl. Acad. Sci. U.S.A.
106
,
3675
(
2009
).
67.
S.
Sengupta
,
S.
Karmakar
,
C.
Dasgupta
, and
S.
Sastry
,
Phys. Rev. Lett.
109
,
095705
(
2012
).
68.
X.
Xia
and
P. G.
Wolynes
,
Proc. Natl. Acad. Sci. U.S.A.
97
,
2990
(
2000
).
69.
H.
Shintani
and
H.
Tanaka
,
Nat. Phys.
2
,
200
(
2006
).
70.
T.
Kawasaki
,
T.
Araki
, and
H.
Tanaka
,
Phys. Rev. Lett.
99
,
215701
(
2007
).
71.
T.
Kawasaki
and
H.
Tanaka
,
J. Phys.: Condens. Matter
22
,
232102
(
2010
).
72.
H.
Tanaka
,
T.
Kawasaki
,
H.
Shintani
, and
K.
Watanabe
,
Nature Mater.
9
,
324
(
2010
).
73.
T.
Kawasaki
and
H.
Tanaka
,
J. Phys.: Condens. Matter
23
,
194121
(
2011
).
74.
M.
Dzugutov
,
S. I.
Simdyankin
, and
F. H. M.
Zetterling
,
Phys. Rev. Lett.
89
,
195701
(
2002
).
75.
U. R.
Pedersen
,
T. B.
Schrøder
,
J. C.
Dyre
, and
P.
Harrowell
,
Phys. Rev. Lett.
104
,
105701
(
2010
).
76.
M.
Leocmach
and
H.
Tanaka
,
Nat. Commun.
3
,
974
(
2012
).
77.
D.
Coslovich
and
G.
Pastore
,
J. Chem. Phys.
127
,
124504
(
2007
).
78.
D.
Coslovich
,
Phys. Rev. E
83
,
051505
(
2011
).
79.
F.
Sausset
,
G.
Tarjus
, and
P.
Viot
,
Phys. Rev. Lett.
101
,
155701
(
2008
).
80.
F.
Sausset
and
G.
Tarjus
,
Phys. Rev. Lett.
104
,
065701
(
2010
).
81.
B.
Charbonneau
,
P.
Charbonneau
, and
G.
Tarjus
,
Phys. Rev. Lett.
108
,
035701
(
2012
).
82.
J.
Kurchan
and
D.
Levine
,
J. Phys. A
44
,
035001
(
2010
).
83.
F.
Sausset
and
D.
Levine
,
Phys. Rev. Lett.
107
,
045501
(
2011
).
84.
A.
Furukawa
and
H.
Tanaka
,
Phys. Rev. Lett.
103
,
135703
(
2009
).
85.
A.
Furukawa
and
H.
Tanaka
,
Phys. Rev. E
84
,
061503
(
2011
).
86.
A.
Furukawa
and
H.
Tanaka
,
Phys. Rev. E
86
,
030501
R
(
2012
).
87.
L.
Berthier
,
D.
Chandler
, and
J. P.
Garrahan
,
EPL
69
,
320
(
2005
).
88.
G.
Szamel
and
E.
Flenner
,
Phys. Rev. E
73
,
011504
(
2006
).
89.
A.
Cavagna
,
T. S.
Grigera
, and
P.
Verrocchio
,
Phys. Rev. Lett.
98
,
187801
(
2007
).
90.
G.
Biroli
,
J. P.
Bouchaud
,
A.
Cavagna
,
T. S.
Grigera
, and
P.
Verrocchio
,
Nat. Phys.
4
,
771
(
2008
).
91.
W.
Kob
,
S.
Roldan-Vargas
, and
L.
Berthier
,
Nat. Phys.
8
,
164
(
2012
).
92.
G. M.
Hocky
,
T. E.
Markland
, and
D. R.
Reichman
,
Phys. Rev. Lett.
108
,
225506
(
2012
).
93.
L.
Berthier
and
W.
Kob
,
Phys. Rev. E
85
,
011102
(
2012
).
94.
S.
Karmakar
,
E.
Lerner
, and
I.
Procaccia
,
Physica A
391
,
1001
(
2012
).
95.
C.
Cammarota
and
G.
Biroli
,
Proc. Natl. Acad. Sci. U.S.A.
109
,
8850
(
2012
).
96.
C.
Dasgupta
,
A. V.
Indrani
,
S.
Ramaswamy
, and
M. K.
Phani
,
Europhys. Lett.
15
,
307
(
1991
).
97.
S.
Franz
and
G.
Parisi
,
J. Phys.: Condens. Matter
12
,
6335
(
2000
).
98.
S. C.
Glotzer
,
V. N.
Novikov
, and
T. B.
Schroder
,
J. Chem. Phys.
112
,
509
(
2000
).
99.
C.
Donati
,
S.
Franz
,
S. C.
Glotzer
, and
G.
Parisi
,
J. Non-Cryst. Solids
307–310
,
215
(
2002
).
100.
S.
Whitelam
,
L.
Berthier
, and
J. P.
Garrahan
,
Phys. Rev. Lett.
92
,
185705
(
2004
).
101.
C.
Toninelli
,
M.
Wyart
,
L.
Berthier
,
G.
Biroli
, and
J. P.
Bouchaud
,
Phys. Rev. E
71
,
041505
(
2005
).
102.
D.
Chandler
,
J. P.
Garrahan
,
R. L.
Jack
,
L.
Maibaum
, and
A. C.
Pan
,
Phys. Rev. E
74
,
051501
(
2006
).
103.
G.
Szamel
and
E.
Flenner
,
Phys. Rev. E
74
,
021507
(
2006
).
104.
L.
Berthier
,
G.
Biroli
,
J. P.
Bouchaud
,
W.
Kob
,
K.
Miyazaki
, and
D. R.
Reichman
,
J. Chem. Phys.
126
,
184503
(
2007
).
105.
L.
Berthier
,
G.
Biroli
,
J. P.
Bouchaud
,
W.
Kob
,
K.
Miyazaki
, and
D. R.
Reichman
,
J. Chem. Phys.
126
,
184504
(
2007
).
106.
G.
Biroli
and
J. P.
Bouchaud
,
Europhy. Lett.
67
,
21
(
2004
).
107.
G.
Biroli
,
J. P.
Bouchaud
,
K.
Miyazaki
, and
D. R.
Reichman
,
Phys. Rev. Lett.
97
,
195701
(
2006
).
108.
G.
Szamel
and
E.
Flenner
,
Phys. Rev. E
81
,
031507
(
2010
).
109.
L.
Berthier
,
G.
Biroli
,
J. P.
Bouchaud
,
L.
Cipelletti
,
E. D.
Masri
,
D.
L’Hôte
,
F.
Ladieu
, and
M.
Pierno
,
Science
310
,
1797
(
2005
).
110.
C.
Dalle-Ferrier
,
C.
Thibierge
,
C.
Alba-Simionesco
,
L.
Berthier
,
G.
Biroli
,
J. P.
Bouchaud
,
F.
Ladieu
,
D.
L’Hôte
, and
G.
Tarjus
,
Phys. Rev. E
76
,
041510
(
2007
).
111.
G.
Brambilla
,
D. E.
Masri
,
M.
Pierno
,
L.
Berthier
,
L.
Cipelletti
,
G.
Petekidis
, and
A. B.
Schofield
,
Phys. Rev. Lett.
102
,
085703
(
2009
).
112.
J. P.
Bouchaud
and
G.
Biroli
,
Phys. Rev. B
72
,
064204
(
2005
).
113.
M.
Tarzia
,
G.
Biroli
,
A.
Lefèvre
, and
J. P.
Bouchaud
,
J. Chem. Phys.
132
,
054501
(
2010
).
114.
C.
Crauste-Thibierge
,
C.
Brun
,
F.
Ladieu
,
D.
L’Hôte
,
G.
Biroli
, and
J. P.
Bouchaud
,
Phys. Rev. Lett.
104
,
165703
(
2010
).
115.
G.
Diezemann
,
Phys. Rev. E
85
,
051502
(
2012
).
116.
M. T.
Cicerone
and
M. D.
Ediger
,
J. Chem. Phys.
103
,
5684
(
1995
).
117.
C. Y.
Wang
and
M. D.
Ediger
,
J. Phys. Chem. B
103
,
4177
(
1999
).
118.
C. Y.
Wang
and
M. D.
Ediger
,
J. Chem. Phys.
112
,
6933
(
2000
).
119.
L. A.
Deschenes
and
D. A.
Vanden Bout
,
J. Phys. Chem. B
106
,
11438
(
2002
).
120.
A. N.
Adhikari
,
N. A.
Capurso
, and
D.
Bingemann
,
J. Chem. Phys.
127
,
114508
(
2007
).
121.
R.
Zondervan
,
F.
Kulzer
,
G. C. G.
Berkhout
, and
M.
Orrit
,
Proc. Natl. Acad. Sci. U.S.A.
104
,
12628
(
2007
).
122.
S. A.
Mackowiak
,
T. K.
Herman
, and
L. J.
Kaufman
,
J. Chem. Phys.
131
,
244513
(
2009
).
123.
S. A.
Mackowiak
,
L. M.
Leone
, and
L. J.
Kaufman
,
Phys. Chem. Chem. Phys.
13
,
1786
(
2011
).
124.
D.
Bingemann
,
R. M.
Allen
, and
S. W.
Olesen
,
J. Chem. Phys.
134
,
024513
(
2011
).
125.
A.
Heuer
and
K.
Okun
,
J. Chem. Phys.
106
,
6176
(
1997
).
126.
127.
B.
Doliwa
and
A.
Heuer
,
Phys. Rev. Lett.
80
,
4915
(
1998
).
128.
J.
Qian
and
A.
Heuer
,
Eur. Phys. J. B
18
,
501
(
2000
).
129.
R.
van Zon
and
J.
Schofield
,
Phys. Rev. E
65
,
011106
(
2001
).
130.
B.
Doliwa
and
A.
Heuer
,
J. Non-Cryst. Solids
307–310
,
32
(
2002
).
131.
E.
Flenner
and
G.
Szamel
,
Phys. Rev. E
70
,
052501
(
2004
).
132.
Y. J.
Jung
,
J. P.
Garrahan
, and
D.
Chandler
,
Phys. Rev. E
69
,
061205
(
2004
).
133.
Y. J.
Jung
,
J. P.
Garrahan
, and
D.
Chandler
,
J. Chem. Phys.
123
,
084509
(
2005
).
134.
S.
Léonard
and
L.
Berthier
,
J. Phys.: Condens. Matter
17
,
S3571
(
2005
).
135.
L. O.
Hedges
,
L.
Maibaum
,
D.
Chandler
, and
J. P.
Garrahan
,
J. Chem. Phys.
127
,
211101
(
2007
).
136.
H.
Mizuno
and
R.
Yamamoto
,
Phys. Rev. E
82
,
030501
R
(
2010
).
137.
H.
Mizuno
and
R.
Yamamoto
,
Phys. Rev. E
84
,
011506
(
2011
).
138.
K.
Kim
and
S.
Saito
,
Phys. Rev. E
79
,
060501
R
(
2009
).
139.
K.
Kim
and
S.
Saito
,
J. Chem. Phys.
133
,
044511
(
2010
).
140.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
USA
,
1999
).
141.
Ultrafast Infrared and Raman Spectroscopy
, edited by
M. D.
Fayer
(
Marcel Dekker
,
2001
).
142.
M.
Khalil
,
N.
Demirdoven
, and
A.
Tokmakoff
,
J. Phys. Chem. A
107
,
5258
(
2003
).
143.
Y.
Tanimura
,
J. Phys. Soc. Jpn.
75
,
082001
(
2006
).
144.
R. M.
Hochstrasser
,
Proc. Natl. Acad. Sci. U.S.A.
104
,
14190
(
2007
).
145.
M.
Cho
,
Two-Dimensional Optical Spectroscopy
(
CRC
,
USA
,
2009
).
146.
E. N.
Senning
and
A. H.
Marcus
,
Annu. Rev. Phys. Chem.
61
,
111
(
2010
).
147.
P.
Hamm
and
M.
Zanni
,
Concepts and Methods of 2D Infrared Spectroscopy
(
Cambridge University Press
,
Cambridge
,
2011
).
148.
M. A.
Berg
,
Adv. Chem. Phys.
150
,
1
(
2012
).
149.
J. B.
Asbury
,
T.
Steinel
,
C.
Stromberg
,
S. A.
Corcelli
,
C. P.
Lawrence
,
J. L.
Skinner
, and
M. D.
Fayer
,
J. Phys. Chem. A
108
,
1107
(
2004
).
150.
J. R.
Schmidt
,
S. A.
Corcelli
, and
J. L.
Skinner
,
J. Chem. Phys.
123
,
044513
(
2005
).
151.
J. J.
Loparo
,
S. T.
Roberts
, and
A.
Tokmakoff
,
J. Chem. Phys.
125
,
194522
(
2006
).
152.
D.
Kraemer
,
M. L.
Cowan
,
A.
Paarmann
,
N.
Huse
,
E. T. J.
Nibbering
,
T.
Elsaesser
, and
R. J. D.
Miller
,
Proc. Natl. Acad. Sci. U.S.A.
105
,
437
(
2008
).
153.
A.
Paarmann
,
T.
Hayashi
,
S.
Mukamel
, and
R. J. D.
Miller
,
J. Chem. Phys.
128
,
191103
(
2008
).
154.
S.
Garrett-Roe
and
P.
Hamm
,
J. Chem. Phys.
128
,
104507
(
2008
).
155.
T.
Yagasaki
and
S.
Saito
,
J. Chem. Phys.
128
,
154521
(
2008
).
156.
T.
Yagasaki
and
S.
Saito
,
Acc. Chem. Res.
42
,
1250
(
2009
).
157.
T.
Yagasaki
and
S.
Saito
,
J. Chem. Phys.
135
,
244511
(
2011
).
158.
F.
Perakis
and
P.
Hamm
,
J. Phys. Chem. B
115
,
5289
(
2011
).
159.
J. T.
King
,
M. R.
Ross
, and
K. J.
Kubarych
,
Phys. Rev. Lett.
108
,
157401
(
2012
).
160.
L.
Berthier
and
G.
Tarjus
,
Phys. Rev. Lett.
103
,
170601
(
2009
).
161.
L.
Berthier
and
G.
Tarjus
,
Phys. Rev. E
82
,
031502
(
2010
).
162.
L.
Berthier
and
G.
Tarjus
,
J. Chem. Phys.
134
,
214503
(
2011
).
163.
W.
Kob
and
H. C.
Andersen
,
Phys. Rev. E
51
,
4626
(
1995
).
164.
W.
Kob
and
H. C.
Andersen
,
Phys. Rev. E
52
,
4134
(
1995
).
165.
G.
Wahnström
,
Phys. Rev. A
44
,
3752
(
1991
).
166.
B.
Bernu
,
Y.
Hiwatari
, and
J. P.
Hansen
,
J. Phys. C
18
,
L371
(
1985
).
167.
B.
Bernu
,
J. P.
Hansen
,
Y.
Hiwatari
, and
G.
Pastore
,
Phys. Rev. A
36
,
4891
(
1987
).
168.
D.
Coslovich
and
G.
Pastore
,
J. Phys.: Condens. Matter
21
,
285107
(
2009
).
169.
T. B.
Schrøder
,
S.
Sastry
,
J. C.
Dyre
, and
S. C.
Glotzer
,
J. Chem. Phys.
112
,
9834
(
2000
).
170.
B. W. H.
van Beest
,
G. J.
Kramer
, and
R. A.
van Santen
,
Phys. Rev. Lett.
64
,
1955
(
1990
).
171.
J.
Horbach
,
W.
Kob
,
K.
Binder
, and
C. A.
Angell
,
Phys. Rev. E
54
,
R5897
(
1996
).
172.
J.
Horbach
and
W.
Kob
,
Phys. Rev. B
60
,
3169
(
1999
).
173.
L.
Berthier
,
G.
Biroli
,
D.
Coslovich
,
W.
Kob
, and
C.
Toninelli
,
Phys. Rev. E
86
,
031502
(
2012
).
174.
M.
Vogel
and
S. C.
Glotzer
,
Phys. Rev. E
70
,
061504
(
2004
).
175.
L.
Berthier
,
Phys. Rev. E
76
,
011507
(
2007
).
176.
S. E.
Abraham
and
B.
Bagchi
,
Phys. Rev. E
78
,
051501
(
2008
).
177.
T.
Abete
,
A.
de Candia
,
E.
Del Gado
,
A.
Fierro
, and
A.
Coniglio
,
Phys. Rev. Lett.
98
,
088301
(
2007
).
178.
A.
Fierro
,
E.
Del Gado
,
A.
de Candia
, and
A.
Coniglio
,
J. Stat. Mech.
2008
,
L04002
.
179.
K.
Kim
,
K.
Miyazaki
, and
S.
Saito
,
EPL
88
,
36002
(
2009
).
180.
K.
Kim
,
K.
Miyazaki
, and
S.
Saito
,
J. Phys.: Condens. Matter
23
,
234123
(
2011
).
181.
R. S. L.
Stein
and
H. C.
Andersen
,
Phys. Rev. Lett.
101
,
267802
(
2008
).
182.
S.
Karmakar
,
C.
Dasgupta
, and
S.
Sastry
,
Phys. Rev. Lett.
105
,
015701
(
2010
).
183.
E.
Flenner
and
G.
Szamel
,
Phys. Rev. Lett.
105
,
217801
(
2010
).
184.
E.
Flenner
,
M.
Zhang
, and
G.
Szamel
,
Phys. Rev. E
83
,
051501
(
2011
).
185.
A.
Onuki
,
Phase Transition Dynamics
(
Cambridge University Press
,
Cambridge
,
2007
).
186.
See supplementary material at http://dx.doi.org/10.1063/1.4769256 for more numerical results of three-time correlation functions ΔF4(k, t1, t2, t3) at various temperatures.

Supplementary Material