Multiple length and time scales of dynamic heterogeneities in model glass-forming liquids: A systematic analysis of multi-point and multi-time correlations

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,⁵⁶ random first-order transition,^57–59 dynamic facilitation,¹⁵ potential energy landscape,⁶⁰ mode-coupling,^61,62 replicated liquid,^63,64 and frustration-limited domain⁶⁵ 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 correlations^{104,105,109–111} and nonlinear susceptibilities^112–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 N₁ particles of component 1 and N₂ particles of component 2 under periodic boundary conditions. The total number density is fixed at ρ = N/V with the total number of particles N = N₁ + N₂ 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,¹⁶⁵ and binary soft-sphere (SS) liquids.^166,167 In addition, we also study a model of a network-forming (NTW) liquid that mimics SiO₂ with short-range spherical potentials.¹⁶⁸ 

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

in which α, β ∈ {1, 2} are the indexes of the particle species. The energy and size ratios are ε₁₂/ε₁₁ = 1.5, ε₂₂/ε₁₁ = 0.5 and σ₁₂/σ₁₁ = 0.8, σ₂₂/σ₁₁ = 0.88, respectively. The masses of the two particle species are equal, m₁ = m₂ = 1. The interaction is truncated at r = 2.5σ_αβ. The reduced units σ₁₁, ε₁₁/k_B, and

$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 N₁ = 800 and N₂ = 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 σ₁/σ₂ = 1/1.2, m₁/m₂ = 1/2, and ε₁/ε₂ = 1, respectively. The Lorentz–Berthelot rules,

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

$m1σ12/ε1$ for length, temperature, and time, respectively. The system consists of N₁ = 500 and N₂ = 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

with the cubic smoothing function v_αβ(r) = B(a − r)³ + C for distances

$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: σ₁/σ₂ = 1/1.2, m₁/m₂ = 1/2, and ε₁/ε₂ = 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 σ₁, ε₁/k_B, 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,

in which l₀ represents the effective particle size defined by 4l₀³ = (2σ₁)³ + 2(σ₁ + σ₂)³ + (2σ₂)³. In the simulation, the total number density is given as ρ = l₀⁻³ with N₁ = N₂ = 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.¹⁶⁸ This model is simple model and imitates SiO₂ 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

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

The mass ratio is determined as m₂/m₁ = 0.57 from the same ratio of O and Si. The units of length, time, and temperature are given as σ₁₁, ε₁₁/k_B, and

$m1σ112/ε11$⁠, respectively. These parameters are adjusted to reproduce the radial distribution functions of the SiO₂ amorphous states. Tetrahedral networks are found to be formed due to the highly asymmetric size ratio and the strong attraction between the different components.¹⁶⁸ The density of the investigated system is ρ = 1.655, with particle numbers N₁ = 1, 000 and N₂ = 2, 000. This value corresponds to the density ρ = 2.37 gÅ⁻³ 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

is calculated for various glass-forming models. Here,

$Δ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 F_s(k, t) at various temperatures is demonstrated in Fig. 1. The wave vector

$k=|k|$ is chosen so that the static structure factor of component 1, S₁₁(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 F_s(k, τ_α) = e⁻¹. In Fig. 2(a), the temperature dependence of τ_α is plotted as a function of the inverse temperature T₀/T. Here T₀ denotes the onset temperature introduced in the earlier work.¹⁶⁰ We set T₀ as T₀ = 0.8, 0.7, 0.3, and 0.5 for the KALJ, WAHN, SS, and NTW models, respectively. Below this onset temperature, F_s(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 T₀. Correspondingly, the temperature dependence of the α-relaxation time obeys the Arrhenius law.¹⁶⁸ In addition, we summarize the power-law behavior as τ_α ∝ |T − T_c|^−Δ, 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 T_c as (Δ, T_c) ≈ (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.¹⁷³ 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.¹⁷² 

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,³³ the overlap function,³⁸ and the intermediate scattering function.³⁹ 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 χ₄(k, t), which is defined as the variance of the self-part of the intermediate scattering function F_s(k, t):

We utilize

$F̂s(k,t)$ expressed as

with

$Fs(k,t)=⟨F̂s(k,t)⟩$⁠. The total value of χ₄(k, t) is approximately proportional to the extension of the spatial correlations in dynamics with k at a given time interval t because χ₄(k, t) investigates the increasing deviation of the two-point correlation function F_s(k, t) from the mean behavior, As is well documented in various studies^101–104 and as demonstrated in Fig. 3, χ₄(k, t) typically has its maximum value at the time scale of τ_α, which increases as the temperature decreases. The growth of χ₄(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 χ₄(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

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:

with the wave vector

$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 W_j(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,³⁸ 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 F_s(k, t) with the wave number marking the main peak of the static structure factor. We also note that any significant differences in S₄(q, t) are not observed if we choose F_s(k, t) as the mobility field. Although recent reports have described the results of S₄(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 S₄(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 = V^1/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 S₄(q, t) on the time scale τ_α. As indicated in Fig. 4, S₄(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 S₄(q, t) at small wave-numbers can be described by the Ornstein–Zernike (OZ) form,

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 S₄(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.¹⁸⁵ 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.⁴¹ 

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 ξ(τ_α) ∼ |T − T_c|^−ν 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.

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 χ(τ_α) ∼ |T − T_c|^−γ 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 χ₄(k, t) by introducing

in which

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, τ₁, τ₂, and τ₃. In practice, ΔF₄(k, t₁, t₂, t₃) reveals the correlation of the two-point correlation function F_s(k, t) between the two time intervals, t₁ = τ₁ and t₃ = τ₃ − τ₂. Furthermore, the progressive changes in the waiting time t₂ = τ₂ − τ₁ of the three-time correlation function ΔF₄(k, t₁, t₂, t₃) allow for an investigation of how the correlated motions change with time t₂. In fact, the first time-interval portion

$δ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 t₂.^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 ΔF₄ 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 F_s(k, t) in Fig. 1. We also provide the diagonal portions of the time evolutions along t = t₁ = t₃ for various t₂ 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 ΔF₄ are similar and that the time evolutions at t₂ occur similarly among all of the glass models studied: (i) As the temperature decreases, the intensity of the three-time correlation ΔF₄(k, t₁, t₂, t₃) gradually increases (see the results for other temperatures in the supplementary material¹⁸⁶). This result is due to the correlations between particles that move slower (faster) during the first time interval t₁ and remain slower (faster) during the second time interval t₃. (ii) The correlations of ΔF₄(k, t₁, t₂, t₃) at t₂ = 0 between the first time interval t₁ and the subsequent time interval t₃ 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 t₂, ΔF₄ 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 t₁ and t₃ for the given waiting time t₂. 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 t₂ = 0. We observe the strong correlations of ΔF₄ at small values of t₁ and t₃ in the NTW model, which approximately corresponds to the time scale of the early β-relaxation at which F_s(k, t) exhibits damped oscillations on a plateau as shown in Fig. 1(d). This well-known finite-size effect in silica glasses^168,171 would be smaller if simulations were carried out for larger systems. Furthermore, the relaxation time scale of ΔF₄(k, t₁, t₂, t₃) depends on the model. As demonstrated in Fig. 7(d), ΔF₄ 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 ΔF₄ 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).

To quantitatively distinguish between the models we quantify the average lifetime of the dynamic heterogeneities using the waiting-time t₂ dependence of ΔF₄. To this end, we define the volume of the heterogeneities by integrating over the two time intervals t₁ and t₃,^138,139

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 t₂ dependence of Δ_hetero(t₂) normalized by Δ_hetero(0) for the KALJ, WAHN, SS, and NTW models. The waiting time t₂ 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.

The dependence of the normalized volume Δ_hetero(t₂)/Δ_hetero(0) on the waiting time t₂ can be fitted by the stretched-exponential function exp [−(t₂/τ_hetero)^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 T₀/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 τ_α.¹³¹ 

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.

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.