Time Scales of Fickian Diffusion and the Lifetime of Dynamic Heterogeneity

Dynamic heterogeneity, believed to be one of the hallmarks of the dynamics of supercooled liquids, is expected to affect the diffusion of the particles comprising the liquid. We have carried out extensive molecular dynamics simulations of two model glass forming liquids in three dimensions to study the time scales at which Fick's law of diffusion starts to set in. We have identified two different Fickian time scales: one at which the mean squared displacement begins to show a linear dependence on time, and another one at which the distribution of particle displacements becomes Gaussian. These two times scales are found to be very different from each other and from the α relaxation time in both systems. An interesting connection is found between one of these Fickian time scales and the time scale obtained from the bond-breakage correlation function. We discuss the relation among these different time scales and their connection to dynamic heterogeneity in the system.

1. Introduction

One of the important unsolved problems in condensed matter physics is to understand the complex dynamics of supercooled liquids near the glass transition. In the last several decades, immense efforts have been made to understand the rapid rise of the viscosity and the slowing down of the dynamics of supercooled liquids as the glass transition is approached [1–12]. It is now well-accepted that the slow dynamics become increasingly heterogeneous near the glass transition. Numerous studies have been performed to understand this behavior [12–14]. While the origin of dynamic heterogeneity is still not fully understood, the lifetime of these heterogeneities are also debated. In a dynamically heterogeneous system, there are regions in which the local dynamics are considerably faster or slower than the dynamics averaged over the whole sample. One can consider any quantity that characterizes the local dynamics over a time scale τ and define the lifetime of dynamical heterogeneity as the minimum value of τ for which this quantity becomes homogeneous in space. In an earlier study [15], it has been shown that the lifetime of dynamic heterogeneity grows faster than the α-relaxation time, whereas in some other works [16–20], the ratio of this time scale with the α-relaxation time is found to depend sensitively on the temperature.

It is well-known that the mean squared displacement (MSD), 〈Δr²〉, computed for a supercooled liquid shows an initial ballistic growth with time and then tends to plateau at an intermediate time scale. The plateau is followed by a diffusive region. The plateau in the MSD appears because of the hindrance of the motion of a particle caused by the “cage” formed by its neighbors. The diffusive regime sets in when the particles manage to break out of the cage, but the dependence of the MSD on time does not become linear until a time scale that is larger than the α-relaxation time. Although the MSD shows diffusive behavior at this time scale, the dynamics of the system remain spatially heterogeneous. In [15], it has been shown that the distribution of the displacements of the particles becomes Gaussian at a time scale that is at least 30 − 40 times larger than the α-relaxation time. True Fickian diffusion starts after this time scale. Two different time scales related to diffusion therefore exist. One of these is the time at which the slope of the MSD vs. time in a log-log plot becomes unity. We call this time scale τ_D. The other one is the time at which the distribution of the displacements of the particles becomes Gaussian. We call this time τ_F. A similar time scale, τ_H, is obtained via the Binder cumulant of the van Hove correlation function.

It is interesting to inquire whether the Fickian time scales defined above are related to the lifetime of dynamic heterogeneity. In recent studies [21, 22], it has been shown that the short-time β-relaxation in supercooled liquid is a cooperative process. It is also shown that the temperature dependence of the length scale associated with the cooperative β-relaxation process is the same as that of the dynamic heterogeneity length scale obtained at the α-relaxation time. Note that the α-relaxation time can be many orders of a magnitude larger than the β-relaxation time, especially at low temperatures. Thus, the dynamics of a supercooled liquid is heterogeneous at time scales that can be as short as the β-relaxation time, and these heterogeneities persist up to time scales that can be much longer. In the experimental studies of Cicerone and Ediger [16, 17], the lifetime of dynamic heterogeneity was defined as the time scale over which the local dynamics in a region where it is slower than the dynamics averaged over the whole system continues to remain so. The difference between this time scale and the α-relaxation time was found [16, 17] to increase with increasing supercooling. Numerical studies [4, 23–25] have suggested that dynamic heterogeneity in glass-forming liquids of binary hard-sphere systems can survive up to a time scale that is larger than the α-relaxation time by a factor of a few tens. In [15], it has been suggested that the time scale obtained from the distribution of single particle displacements may provide a lower bound for the lifetime of dynamic heterogeneities present in the system, but the actual lifetime of the heterogeneities may be much larger. The violation of the Stokes-Einstein relation between the diffusion coefficient and the viscosity [26–28], another poorly understood phenomenon in glassy dynamics, is also believed to be closely related to dynamic heterogeneity.

In this paper, we have addressed some of the issues mentioned above via extensive molecular dynamics simulations of two well-known model glass forming liquids in three spatial dimensions. The primary objective of this study is to estimate the time scales τ_D, τ_F, and τ_H in different generic glass-forming liquids and to compare them with other well-known important time scales in the system, such as the α-relaxation time τ_α and the bond-breakage time scale τ_BB (to be defined later) We have looked for the existence of correlations among these different time scales and tried to figure out their consequence in understanding the dynamics of glass forming liquids approaching the glass transition. We find that the time scales of diffusion, τ_D, τ_F, and τ_H, are longer than τ_α in the temperature range considered here. The temperature dependence of τ_F is similar to that of τ_H and their growth with decreasing temperature is faster than that of τ_α. The growth of τ_D as the temperature is reduced but is, on the other hand, slower than that of τ_α. We find that the relations between pairs of these time scales are described by power laws at low temperatures, so that the glass transition temperatures obtained from Vogel-Fulcher-Tammann [29, 30] (VFT) or mode coupling [31, 32] (power law) fits to the temperature dependence of these time scales and are the same within error bars. We also find an intriguing connection between τ_D and the bond-breakage time scale τ_BB. The implications of these findings for the violation of the Stokes-Einstein relation and the kinetic fragility are discussed.

The rest of the paper is organized as follows. We present the details of the models studied and the simulation methodology in section 2. We then introduce different correlation functions that are computed to estimate different time scales and present our results for two different generic glass-forming liquids. Finally, we conclude with a discussion of the importance of these different time scales for understanding the intricate role played by dynamic heterogeneity in the dynamics of glass-forming liquids.

2. Models and Methods

We have studied two different model glass forming liquids in three dimensions.

• 3dKA Model [33]: This is a well-known 80:20 binary Lennard-Jones mixture in three dimensions. Particles interact via the following pairwise interaction

$\begin{array}{ll}{V}_{\alpha \beta }(r)=4{ϵ}_{\alpha \beta }\left[{\left(\frac{{\sigma }_{\alpha \beta }}{r}\right)}^{12}-{\left(\frac{{\sigma }_{\alpha \beta }}{r}\right)}^6\right],& (1)\end{array}$

where α, β ∈ {A, B} and ϵ_AA = 1.0, ϵ_AB = 1.5, ϵ_BB = 0.5, σ_AA = 1.0, σ_AB = 0.80, σ_BB = 0.88. We cut off the interaction potential at 2.50σ_αβ. A quadratic polynomial is used to make the potential and its first two derivatives smooth at cutoff distance. The temperature range for the simulations is T ∈ [0.450 − 1.000] and the system size is N = 8, 000.

• 3dR10 Model [34]: It is a 50:50 binary mixture where the particles interact via

$\begin{array}{ll}{V}_{\alpha \beta }(r)={ϵ}_{\alpha \beta }{\left(\frac{{\sigma }_{\alpha \beta }}{r}\right)}^n,& (2)\end{array}$

with n = 10. The potential is cut off at a distance 1.38σ_αβ and again, a quadratic polynomial is used to make potential and its first two derivative smooth at the cutoff distance. Here ϵ_αβ = 1.0, σ_AA = 1.0, σ_AB = 1.22, and σ_BB = 1.40. The temperature range for the simulations is T ∈ [0.520 − 1.000] and N = 1, 0000.

For all the studied models we have performed NVT molecular dynamics (MD) simulations with a modified leap-frog algorithm. We used the Berendsen thermostat to keep the temperature constant. Our integration time step is dt = 0.005. The number of MD steps are 10⁷ − 10⁸, depending on the temperature.

3. Correlation Functions

3.1. Overlap Correlation Function

We characterize the dynamics by calculating the well-known two-point density-density correlation function Q(t). It measures the overlap between two configurations separated by time t. The self-part of this correlation function is defined as

$\begin{array}{ll}{Q}_s(t)=\frac{1}{N}\left[\langle {\displaystyle \sum _{i=1}^N}w\left(|{\overrightarrow{r}}_i(t)-{\overrightarrow{r}}_i(0)|\right)\rangle \right],& (3)\end{array}$

where w(x) is a window function with

$\begin{array}{ll}w(x)=\{\begin{array}{cc}1.0\hfill & \text{if }x\le a\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}\hfill & (4)\end{array}$

The value of a is chosen to be 0.30 which is close to the value at which the MSD forms a plateau. The window function is chosen to remove the possible decorrelation due to the vibrational motion of the particles inside their cage. A different choice of the parameter a does not change the results qualitatively. The α-relaxation time τ_α is obtained from Q_s(t) using the condition Q_s(t = τ_α) = exp(−1).

3.2. van Hove Correlation Function

The van Hove correlation function measures the probability that a particle has displacement x after time t. The self-part of the van Hove correlation function [35] is defined as

$\begin{array}{ll}{G}_s(x,t)={\displaystyle \sum _{i=1}^N}\delta (x-x_i(t)+x_i(0)),& (5)\end{array}$

where x_i(t) is the x-component of the position vector of particle i at time t.

3.3. Bond-Breakage Correlation Function

The bond-breakage correlation function is defined in the following way [36–38]. At t = t₀ a pair of particle i (of type α) and j (of type β) is considered to be bonded if

$\begin{array}{ll}{r}_{ij}(t_0)=|{\overrightarrow{r}}_i(t_0)-{\overrightarrow{r}}_j(t_0)|\le {\sigma }_{\alpha \beta }.& (6)\end{array}$

If r_ij(t) ≤ 1.35σ_αβ, the bond is said to have survived at time t. We calculate the number of bonds that survive at time t and the bond-breakage correlation function is defined as the ratio of this number to the initial number of bonds. The bond-breakage time scale τ_BB is the time at which the correlation function goes to 1/e of its initial value.

4. Results

We show below the results for only one type of particle (type A). The results for the other type of particle (type B) are qualitatively similar [15]. In the 3dKA model (80:20 mixture), the number of B-type particles is substantially smaller, resulting in larger statistical fluctuations. So, we prefer to show the results for A-type particles.

4.1. Time Scale From the Mean Square Displacement

As mentioned before, the mean square displacement (MSD) shows three distinct regimes for a supercooled liquid: initial ballistic growth i.e., 〈Δr²〉 ~ t², then a plateau followed by a diffusive regime i.e., 〈Δr²〉 ~ t. So, a plot of the logarithmic derivative of the MSD with respect to time starts from 2, goes to a minimum and then eventually goes to 1, corresponding to the above three distinct regimes. The time scale where it goes to 1 is of importance as this is the time at which the system begins to show diffusive behavior. We have calculated this time scale, τ_D, for different temperatures for both the model systems. In the top panel of Figure 1, we show the logarithmic derivative of the MSD as a function of time for all the studied temperatures for the 3dR10 model. We have taken $\frac{dln\langle \Delta r^2\rangle }{dlnt}=0.97$ to define the time scale τ_D. The temperature dependence of this time scale is shown in the bottom panel of the same figure. The temperature dependence can be fitted very well using the VFT form with the divergence temperature T_VFT≃0.406, which is close to ~0.404, the value obtained using τ_α. We consider two other characteristic temperatures—the mode coupling critical temperature T_c obtained from fitting the temperature dependence of a time scale to a power-law divergence, and the “glass transition temperature” T_g defined as the temperature at which the time scale being considered reaches the value 10⁶ (this is similar to the definition of the experimental glass transition temperature). The values of T_c and T_g obtained from the temperature dependence of τ_D are ~0.495 and ~0.484, respectively. These values are close to those (~0.491 and ~0.481, respectively) obtained using τ_α.

FIGURE 1

Figure 1. (Top) Log-derivative of MSD is plotted vs. time for the 3dR10 model. The horizontal line corresponds to 0.97. (Bottom) The dependence of the time scale τ_D on the temperature. The line passing through the data points is a fit of the data to the VFT form, ${\tau }_{\alpha }={\tau }_{0}exp\left[\frac{A}{T-T_{VFT}}\right]$.

To understand the mutual relationship between τ_D and τ_α, we have plotted the ratio of these two time scales (τ_D/τ_α) as a function of temperature in the top panel of Figure 2. As the temperature is decreased, this ratio first increases, reaches a maximum near T = 0.8, and then decreases as the temperature is lowered further. The temperature at which τ_D/τ_α peaks is close to the onset temperature T_on [39] for this model. Thus, at temperatures lower than T_on which corresponds to the temperature below which the effects of the potential energy landscape becomes important, the diffusion time scale τ_D grows more slowly than the structural relaxation time τ_α with decreasing temperature.

FIGURE 2

Figure 2. The ratio τ_D/τ_α is plotted as a function of T for the 3dR10 model. The inset shows a comparison of this ratio with the ratio of τ_BB/τ_α for the same model.

We have also calculated the bond-breakage time scale τ_BB following the method described in section 3. It is the time scale where most of the particles have gone though some amount of reshuffling of their neighbors, as measured using the breakage of nearest neighbor bonds. The temperature dependence of τ_BB/τ_α is similar to that of τ_D/τ_α, as shown in the inset of Figure 2. In particular, these two ratios exhibit nearly identical temperature dependence for values of T lower than that at which they exhibit a peak. This observation suggests an intriguing connection between diffusion and bond breakage processes at temperatures lower than the onset temperature. Possible significance of this observation is discussed in section 8.

To calculate τ_D, we have chosen the cutoff value of the logarithmic derivative of the MSD to be 0.97. Ideally, τ_D should be calculated where this value actually goes to unity. This is quite challenging as for low temperatures, the derivative does not reach the value of unity in our simulation time scale. Even if the derivative goes to unity within the simulation time scale, the noise in its measurement makes it difficult to accurately determine the time at which it first reaches this value. To check whether the non-monotonic behavior of the temperature dependence of τ_D/τ_α depends on how we calculate τ_D, we have done a systematic study of this ratio by changing the cutoff. Figure 3 shows the cutoff dependent ratio for the 3dKA model systems. From the plot, it is clear that the non-monotonic behavior is a feature that is independent of the choice of the cutoff: the non-monotonicity increases systematically as the cutoff approaches unity.

FIGURE 3

Figure 3. Cutoff dependence of τ_D/τ_α for the 3dKA model. It is clear that although the value of this ratio depends on the cutoff, the non-monotonic temperature dependence is present for all values of the cutoff.

4.2. Time Scale From the Distribution of Single-Particle Displacements

Next, we estimate the time scale of the onset of Fickian diffusion by following the procedure discussed in [15]. To study the onset of Fickian diffusion, it is important to calculate the distribution of single particle displacements, P(log₁₀(δr), t) at time t. As discussed in [15], this distribution is connected to the self part of the van Hove correlation function as

$\begin{array}{ll}P({\text{log}}_{10}(\delta r),t)=\text{ln }104\pi \delta r^3{G}_s(\delta r,t).& (7)\end{array}$

If the motion of the particles is governed by Fick's law of diffusion, the distribution of single particle displacements i.e., P[log₁₀(δr, t)], should become Gaussian and remain so at longer times. Moreover, the peak value of the distribution should reach a value ≈2.13. Any deviation from this behavior would indicate non-Fickian diffusion as well as dynamic heterogeneity [15, 40]. Figure 4 shows plots for P(log₁₀(δr, t)) at different times for two different temperatures T = 0.450 and 1.000 for the 3dKA model. For the lower temperature (T = 0.450), deviations from Gaussian behavior are clearly visible at smaller times, whereas the distribution approaches a Gaussian shape for longer times and the peak value of the distribution also approaches ≈2.13.

FIGURE 4

Figure 4. (Top) Distribution of the logarithm of single particle displacements at different times at T = 0.450 for the 3dKA model. For relatively short times, the distribution shows signs of bimodality, whereas at longer times it tends to be more Gaussian. (Bottom) Similar plot at T = 1.000.

At the higher temperature, T = 1.000, deviations from the Gaussian shape are not very clear, but the peak value reaches the value ~2.13 only after a few τ_α. This suggests that even at high temperatures, the dynamics are affected by the presence of spatial heterogeneity at time scales larger than τ_α. We define the time scale of onset of Fickian diffusion as the time at which the peak value of the distribution goes to ≈1.92, similar to the criterion adopted in [15] and refer to this time scale as τ_F. In the left panel of Figure 5, we show plots of the peak value of P(log₁₀(δr), t) with increasing time for all the studied temperatures. In the right panel we show the temperature dependence of the time scale τ_F. The temperature dependence of this time scale can be fitted very well to the VFT formula. The line passing through the data points in the right panel shows the VFT fit with T_VFT≃0.289.

FIGURE 5

Figure 5. (Left) Value of the maximum of P(log₁₀(Δr), t) is plotted against time t for different temperatures for the 3dKA model. The time where it reaches 1.92 is defined as the Fickian time scale τ_F. (Right) Variation of the logarithm of the time scale τ_F with temperature T. The line passing through the data points is a fit to the VFT form.

To see the mutual relationship between τ_α and τ_F, in Figure 6 we have plotted the ratio (τ_F/τ_α) for all the studied temperatures. From the plot, one can see that τ_F increases much more rapidly than τ_α when the temperature is decreased. The ratio changes from 5 to 25 in the temperature window T ∈ [1.000 − 0.450] as shown in Figure 6.

FIGURE 6

Figure 6. The ratio of τ_F/τ_α is plotted as a function of T for the 3dKA model. Note that τ_F increases much faster than τ_α with decreasing temperature.

Ideally, one should use a cutoff at 2.13 in order to calculate τ_F. We could not do that because reaching this value within our simulation time scale is not possible. We have, however, checked that the difference between τ_F and the other time scale τ_D related to diffusion is not an artifact of using a cutoff lower than 2.13 in the calculation of τ_F and a cutoff lower than 1.0 in the calculation of τ_D. It is evident from the data that τ_F is larger than τ_D and it increases faster than τ_D as the temperature is reduced. There are indications that these two time scales approach each other at temperatures higher than 3.0.

The van Hove function is expected to have a Gaussian form for Fickian diffusion. In the top panel of Figure 7, we have plotted the van Hove function for the 3dKA model at T = 0.450 for different times. From the plots, one can see evidence for non-Gaussian behavior for smaller times and an approach to a Gaussian form for longer times. The bottom panel shows similar plots for T = 1.000. To get a second measure of the degree of non-Gaussianity, we have calculated the Binder cumulant of the van Hove function for different times. The Binder cumulant, which provides a measure of the excess kurtosis of the distribution, is defined as

$\begin{array}{ll}B(T,t)=1-\frac{\langle x^4\rangle }{3{\langle x^2\rangle }^2}.& (8)\end{array}$

If the van Hove function is perfectly Gaussian, the Binder cumulant should be zero. For that to happen, one needs to run the simulations for very long times. Here we have defined a time scale τ_H as the time where the Binder cumulant of the van Hove function becomes ≈0.07. In the top panel of Figure 8, we show the Binder cumulant for all studied temperatures at different times for the 3dKA model. The horizontal line in the figure corresponds to the value of 0.07. The bottom panel shows the ratio of τ_H and τ_α as a function of temperature. The temperature dependence of τ_H/τ_α is similar to that of τ_F/τ_α. This is not surprising because both τ_F and τ_H correspond to the time at which the distribution of particle displacements becomes Gaussian. Note that τ_H is larger than τ_F for our choice of the cutoffs used for measuring these time scales. The values of the ratio τ_H/τ_α lie in the range ~6 − 35 for temperatures in the range T ∈ [1.000 − 0.450]. The characteristic temperatures, T_VFT, T_c and T_g, obtained from the different time scales in the 3dKA model are listed in Table 1. The values of T_VFT and T_c obtained from the α-relaxation time are found to be close to those obtained from the other time scales. The value of T_g, which is sensitive to the absolute value of the time scale, depends of the time scale used to define it. At this point it is worth mentioning that another time scale τ₄, obtained from the four-point susceptibility χ₄(t) [3] (τ₄ is the time scale where the χ₄(t) attains its maximum value), is proportional to τ_α with proportionality constants close to one [3]. Therefore, this time scale does not provide any new information.

FIGURE 7

Figure 7. (Top) The van Hove correlation function at different times at T = 0.450 for the 3dKA model. (Bottom) Similar plot for T = 1.000.

FIGURE 8

Figure 8. (Top) The Binder Cumulant of the van Hove correlation function for different times and all studied temperatures for the 3dKA model. The horizontal line corresponds to B(T, t) = 0.07. The time scale τ_H is defined as the time at which B(T, t) attains this value. (Bottom) The temperature dependence of the ratio τ_H/τ_α.

TABLE 1

Table 1. The values of T_VFT, T_c, and T_g for all the three time scales for 3dKA model.

As suggested in [15], the time scales τ_F and τ_H may be thought of as lower bounds of the lifetime of dynamic heterogeneity. As reported in [41, 42], the distribution of diffusion constants is a good measure of dynamic heterogeneity in the system. The distribution of diffusion constants is obtained from the self-part of the van Hove function using Lucy iteration [43]. We have followed the method described in [41, 42]. In the top and bottom panels of Figure 9, we show the distribution of diffusion constants for T = 0.450 and T = 1.000, respectively, for different times for the 3dKA model. For a high temperature e.g., T = 1.000, the distribution does not have multiple peaks but the broadness of the distribution reflects the underlying heterogeneity. For T = 0.450, one can see that there are multiple peaks in the distributions even at time scales longer than τ_F. It approaches unimodal behavior at much longer time scales. So, heterogeneity survives even after the system starts to follow Fick's law of diffusion, which is in agreement with the conclusion of previous studies [15].

FIGURE 9

Figure 9. (Top) Distribution of diffusion constants for T = 0.450 for 3dKA model. (Bottom) Similar plot for T = 1.000.

5. Relation Between Different Time Scales

In the previous section, we found that the temperature dependence of the different time scales (τ_α, τ_D, τ_F, and τ_H) can be well-fitted by the VFT form with very similar values of the divergence temperature T_VFT. This suggests that the different time scales are related to one another via power laws. In Figure 10, we have shown such a relationship between different time scales. The linear dependence in a log-log plot confirms that a power-law relation describes the data quite well, especially at low temperatures. The exponent m that describes the power-law relation between one of the time scales of diffusion and the α-relaxation time is ≃0.83 for τ_D, implying that τ_D increases more slowly than τ_α with decreasing temperature. The values of m for τ_F and τ_H obtained from the fits, m≃1.21 and m≃1.17, respectively, are close to one another. This is consistent with the observation, mentioned above, that these two time scales have a very similar temperature dependence. The larger than unity values of m for these time scales imply that their growth, with decreasing temperature, is faster than the growth of τ_α. Since τ_F and τ_H are believed to provide lower bounds to the lifetime of dynamic heterogeneity, this result implies that dynamic heterogeneities exist over time scales that are much longer than the α-relaxation time at low temperatures. The values of m for τ_D, τ_F and τ_H also imply that τ_F and τ_H increase faster than τ_D as the temperature is decreased, indicating that at low temperatures, there is a large time interval over which the MSD is linear in time, but the distribution of particle displacements is not Gaussian. Yet another consequence of the observed power-law dependence of the diffusion-related time scales on τ_α is that the transition temperature of mode-coupling theory obtained by fitting the temperature dependence of τ_D, τ_F or τ_H to a power law would be the same as that obtained from a power-law fit to the temperature dependence of τ_α. However, the exponent of the power-law would depend on which time scale is considered.

FIGURE 10

Figure 10. Different relaxation times, τ_D, τ_F, and τ_H, are plotted against the α-relaxation time τ_α for 3dKA model. Reasonable linear dependence in a log-log plot suggests that they are related to one another via power laws. The lines passing through the data points are best fits to power laws.

6. Stokes-Einstein Relation and Kinetic Fragility

The violation of the Stokes-Einstien relation [26–28] between the diffusion coefficient D and the viscosity (or a suitable time scale τ) in deeply supercooled liquids [41, 42] is believed to be closely related to dynamic heterogeneity. It is interesting to examine the extent of violation of the SE relation when one of the time scales studied above is considered to be the relevant time scale τ. It is known from earlier work [41] that D∝1/τ_α (i.e., the SE relation is satisfied) at high temperatures, but a modified relation, $D\propto {\tau }_{\alpha }^{-1+\omega }$, where ω is called the fractional SE exponent, is satisfied at lower temperatures. The value of the exponent (1 − ω) is 0.78 and 0.75 for 3dKA and 3dR10 models, respectively.

The power-law dependence of the time scales of diffusion on τ_α implies that these time scales also satisfy a modified SE relation with values of the fractional SE exponent given by (1 − ω)/m. Since m for τ_D for the 3dKA model is close to the value of (1 − ω), the SE relation is expected to be valid for this time scale (i.e., D∝1/τ_D) at low temperatures. On the other hand, the fractional SE exponent is expected to be small (close to 0.65) for the two other time scales. The plots of the diffusion coefficient vs. different time scales shown in Figure 11 illustrate this behavior. The SE relation is found to be satisfied for τ_D, whereas the fractional SE exponent for τ_F is 0.65. Similar behavior is found for the 3dR10 model, as illustrated in the bottom panel of Figure 11. The SE relation between D and τ_D≃τ_BB was also reported in an earlier study [44].

FIGURE 11

Figure 11. (Top) The diffusion constant D is plotted against different time scales for the 3dKA model system. The lines passing through the data points are power-law fits. (Bottom) Similar plot for the 3dR10 model.

As noted above, τ_α satisfies the SE relation at temperatures higher than the onset temperature. Since the temperature dependence of the diffusion-related time scales is different from that of τ_α, even at these higher temperatures, these time scales do not satisfy the SE relation at temperatures near the onset temperature. The SE relation may be restored for these time scales at temperatures substantially higher than the onset temperature.

The power-law dependence of the diffusion-related time scales on τ_α also implies that the kinetic fragility parameter κ, obtained from a fit of the temperature dependence of a time scale τ to the VFT form, τ(T) = τ₀exp[1/(κ(T/T_VFT − 1)], for these time scales is different from that for τ_α. The value of κ for one of these time scales would be 1/m times the value for τ_α, implying that the fragility for τ_D is higher than that for τ_α, whereas the temperature dependence of τ_F and τ_H corresponds to a smaller value of the fragility.

7. System-Size Dependence of Time Scales

In Figure 12, we show the system-size dependence of τ_D (left panel) at different temperatures for the 3dKA model. Almost no system-size dependence is observed at the higher temperatures, and only a small dependence is found at the lowest temperature for system sizes in the range N ∈ [500 − 10, 000]. We have done simulations for smaller system sizes in the range N ∈ [100 − 1, 000] for the 3dR10 model. Note that for this model, one can simulate systems with smaller sizes as the inter-particle potential has shorter range. The results are shown in the right panel of Figure 12 where one can clearly see the evidence for an initial increase of the time scale with system size at low temperatures. Note that τ_α has a very different system size dependence—its value decreases with increasing system size at low temperatures [7]. This observation is puzzling and warrants further investigation.

FIGURE 12

Figure 12. (Left) System size dependence of τ_D at different temperatures for the 3dKA model system. The observed system size dependence for N ∈ [500 − 10, 000] is very weak and only small dependence is observed at the low temperature. (Right) Zoomed in version of similar plot but for small system sizes N ∈ [100 − 1, 000] for the lowest four temperatures for the 3dR10 model. Note for the 3dR10 model, one can simulate smaller sizes as the range of the inter particle potential is small. See text for details.

8. Summary and Conclusions

In summary, we have estimated different time scales related to diffusion for two well-known glass-forming liquids. One of these time scales, τ_D, is the time at which the dependence of the MSD on time becomes linear. The other, τ_F or τ_H, represents the time at which the distribution of particle displacements becomes Gaussian. The actual values of these time scales depend on the cutoff or tolerance factor used in their measurement from MD simulations. However, our investigation of the dependence of the time scales on the cutoff suggests that the qualitative behavior we find is independent of the choice of the cutoff. The results obtained for the two model systems are similar. We show that both τ_D and τ_F are larger than the α-relaxation time τ_α in the temperature range considered in our simulations. At low temperatures (T lower than the onset temperature T_on), the growth of τ_D with decreasing temperature is slower, whereas the growth of τ_F is faster than that of τ_α. Between the two diffusion-related time scales, τ_F is substantially larger than τ_D, indicating that there is a long time interval over which the MSD is linear in time, but the distribution of the displacements is not Gaussian. If both linear dependence of the MSD on time and Guassian distribution of displacements are considered to be essential features of Fickian diffusion, then the longer time scale τ_F should be identified as the time at which Fickian diffusion sets in. Our results then imply that Fickian diffusion sets in at times much longer than the α-relaxation time in deeply supercooled liquids. This is in agreement with earlier results [15].

We find that the ratio τ_D/τ_α is a non-monotonic function of temperature and it peaks at a temperature close to the onset temperature T_on at which the dynamics starts being influenced by the energy landscape. This observation suggests an intriguing connection between the behavior of the MSD as a function of time and the influence of the energy landscape on the dynamics. Another interesting observation is that τ_D is close to the bond-breakage time scale τ_BB for temperatures lower than T_on.

The diffusion-related time scales τ_D and τ_F and the α-relaxation time τ_α are mutually correlated and the relation between pairs of them is well-described by a power law at low temperatures. This implies that if any one of them diverge at some temperature, then the others also diverge at the same temperature. Thus, the VFT temperature and the critical temperature of mode-coupling theory are the same for all these time scales. The power-law relation also implies that the fractional SE exponent and the kinetic fragility associated with the diffusion-related time scales are simply related to those obtained for τ_α. An interesting result in this context is that the time scale τ_D and the diffusion coefficient D satisfy the SE relation at temperatures lower than T_on. Similar results, τ_D≃τ_BB∝1/D, have been found in simulations of supercooled water [45, 46] and supercooled silica [47]. These results suggest that the observed relations among τ_D, τ_BB, and D are probably universal in the sense that they are satisfied in all supercooled liquids.

On the other hand, the Fickian time scale τ_F exhibits strong violation on the SE relation, with a fractional SE exponent close to 0.65. This is puzzling because one would naively expect the SE relation to be valid when the diffusion is Fickian. However, the violation of the SE relation is believed to be caused by dynamic heterogeneity and the lifetime of dynamic heterogeneity may be longer than τ_F. We have found some evidence for this from a calculation of the distribution of local diffusion constants (Figure 9) which exhibits multiple peaks for time scales longer than τ_F, indicating that the lifetime of the heterogeneity is longer than the Fickian time scale. This is consistent with the suggestion [15] that τ_F is a lower bound to the lifetime of dynamic heterogeneity.

At present, there is no clear understanding of the origin of most of the observations described above. A few conclusions can be drawn from these observations. The equality of the time scales τ_D and τ_BB at temperatures lower than the onset temperature indicates that the onset of linear diffusion (MSD proportional to time) depends crucially on bond-breaking events that lead to the escape of a particle from the “cage” formed by its neighbors. The fact that the structural relaxation time τ_α increases faster than τ_BB at low temperatures (the ratio τ_BB/τ_α decreases as the temperature is reduced below the onset temperature) suggests that structural relaxation involves additional processes that are slower than bond-breaking. The observations that the time scale τ_F of Fickian diffusion is much longer than τ_D and the lifetime of dynamic heterogeneity obtained from the distribution of local diffusion constants is even longer indicate that dynamic heterogeneities persist in the system and continue to influence diffusive behavior over time scales that are much longer than the time scale of linear diffusion.

These observations suggest the following picture for the dynamics over time scales longer than τ_D but shorter than τ_F: the dynamically heterogeneous system consists of groups of particles with the property that the particles belonging to a particular group exhibit linear diffusion with diffusion constant D, but the value of D is different for different groups. Under these circumstances, the mean square displacement averaged over all the particles will be linear in time with a diffusion constant that is a properly weighted average of the diffusion constants of individual groups. However, the distribution of the displacements of all the particles will not be Gaussian, even if the distribution is Gaussian for the particles in a particular group—a Gaussian distribution will be obtained only after a longer time scale τ_F at which the heterogeneity in the values of the diffusion constants disappears. This time scale τ_F is clearly related to the lifetime of dynamic heterogeneity, but as far as we know, a theoretical understanding of the magnitude and temperature dependence of τ_F is not yet available.

All these observations raise several interesting questions about the characteristics of the dynamics of supercooled liquids. Answers to these questions will be important in understanding the role of dynamic heterogeneity in glassy relaxation.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

SK and CD designed the project. RD carried out the research. RD, CD, and SK analyzed the data and wrote the paper. All authors contributed to the article and approved the submitted version.

Funding

This project was funded by intramural funds at TIFR Hyderabad from the Department of Atomic Energy (DAE). Support from Swarna Jayanti Fellowship grants DST/SJF/PSA-01/2018-19 and SB/SFJ/2019-20/05 are also acknowledged.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Cavagna A. Supercooled liquids for pedestrians. Phys Rep. (2009) 476:51–124. doi: 10.1016/j.physrep.2009.03.003

CrossRef Full Text | Google Scholar

2. Berthier L, Biroli G. Theoretical perspective on the glass transition and amorphous materials. Rev Mod Phys. (2011) 83:587–645. doi: 10.1103/RevModPhys.83.587

CrossRef Full Text | Google Scholar

3. Karmakar S, Dasgupta C, Sastry S. Growing length scales and their relation to timescales in glass-forming liquids. Annu Rev Condens Matter Phys. (2014) 5:255–84. doi: 10.1146/annurev-conmatphys-031113-133848

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Berthier L, Biroli G, Bouchaud JP, Cipelletti L, Masri DE, L'Hôte D, et al. Direct experimental evidence of a growing length scale accompanying the glass transition. Science. (2005) 310:1797–800. doi: 10.1126/science.1120714

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Biroli G, Bouchaud JP, Miyazaki K, Reichman DR. Inhomogeneous mode-coupling theory and growing dynamic length in supercooled liquids. Phys Rev Lett. (2006) 97:195701. doi: 10.1103/PhysRevLett.97.195701

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Biroli G, Bouchaud JP, Cavagna A, Grigera TS, Verrocchio P. Thermodynamic signature of growing amorphous order in glass-forming liquids. Nat Phys. (2008) 4:771. doi: 10.1038/nphys1050

CrossRef Full Text | Google Scholar

7. Karmakar S, Dasgupta C, Sastry S. Growing length and time scales in glass-forming liquids. Proc Natl Acad Sci USA. (2009) 106:3675–9. doi: 10.1073/pnas.0811082106

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Karmakar S, Dasgupta C, Sastry S. Length scales in glass-forming liquids and related systems: a review. Rep Prog Phys. (2015) 79:016601. doi: 10.1088/0034-4885/79/1/016601

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kirkpatrick TR, Thirumalai D, Wolynes PG. Scaling concepts for the dynamics of viscous liquids near an ideal glassy state. Phys Rev A. (1989) 40:1045–54. doi: 10.1103/PhysRevA.40.1045

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Lubchenko V, Wolynes PG. Theory of structural glasses and supercooled liquids. Annu Rev Phys Chem. (2007) 58:235–66. doi: 10.1146/annurev.physchem.58.032806.104653

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Biroli G, Bouchaud J. Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications. Hoboken, NJ: Wiley publication (2012).

Google Scholar

12. Ediger MD. Spatially heterogeneous dynamics in supercooled liquids. Annu Rev Phys Chem. (2000) 51:99–128. doi: 10.1146/annurev.physchem.51.1.99

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Richert R. Heterogeneous dynamics in liquids: fluctuations in space and time. J Phys. (2002) 14:R703. doi: 10.1088/0953-8984/14/23/201

CrossRef Full Text | Google Scholar

14. Andersen HC. Molecular dynamics studies of heterogeneous dynamics and dynamic crossover in supercooled atomic liquids. Proc Natl Acad Sci USA. (2005) 102:6686–91. doi: 10.1073/pnas.0500946102

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Szamel G, Flenner E. Time scale for the onset of Fickian diffusion in supercooled liquids. Phys Rev E. (2006) 73:011504. doi: 10.1103/PhysRevE.73.011504

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Cicerone MT, Ediger MD. Relaxation of spatially heterogeneous dynamic domains in supercooled orthoterphenyl. J Chem Phys. (1995) 103:5684–92. doi: 10.1063/1.470551

CrossRef Full Text | Google Scholar

17. Wang CY, Ediger MD. How long do regions of different dynamics persist in supercooled o-Terphenyl? J Phys Chem B. (1999) 103:4177–84. doi: 10.1021/jp984149x

CrossRef Full Text | Google Scholar

18. Deschenes LA, Vanden Bout DA. Single-molecule studies of heterogeneous dynamics in polymer melts near the glass transition. Science. (2001) 292:255–8. doi: 10.1126/science.1056430

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Böhmer R, Hinze G, Diezemann G, Geil B, Sillescu H. Dynamic heterogeneity in supercooled ortho-terphenyl studied by multidimensional deuteron NMR. Europhys Lett. (1996) 36:55.

Google Scholar

20. Böhmer R, Diezemann G, Hinze G, Sillescu H. A nuclear magnetic resonance study of higher-order correlation functions in supercooled ortho-terphenyl. J Chem Phys. (1998) 108:890–9. doi: 10.1063/1.475452

CrossRef Full Text | Google Scholar

21. Das R, Tah I, Karmakar S. Possible universal relation between short time β-relaxation and long time α-relaxation in glass-forming liquids. J Chem Phys. (2018) 149:024501. doi: 10.1063/1.5033555

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Karmakar S, Dasgupta C, Sastry S. Short-time beta relaxation in glass-forming liquids is cooperative in nature. Phys Rev Lett. (2016) 116:085701. doi: 10.1103/PhysRevLett.116.085701

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Flenner E, Szamel G. Dynamic heterogeneity in a glass forming fluid: susceptibility, structure factor, and correlation length. Phys Rev Lett. (2010) 105:217801. doi: 10.1103/PhysRevLett.105.217801

CrossRef Full Text | Google Scholar

24. Lačević N, Starr FW, Schrøder TB, Glotzer SC. Spatially heterogeneous dynamics investigated via a time-dependent four-point density correlation function. J Chem Phys. (2003) 119:7372–87. doi: 10.1063/1.1605094

CrossRef Full Text | Google Scholar

25. Crauste-Thibierge C, Brun C, Ladieu F, L'Hôte D, Biroli G, Bouchaud JP. Evidence of growing spatial correlations at the glass transition from nonlinear response experiments. Phys Rev Lett. (2010) 104:165703. doi: 10.1103/PhysRevLett.104.165703

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Hansen JP, McDonald JR. Theory of Simple Liquids. 2nd ed. London: Academic Press (1986).

Google Scholar

27. Einstein A. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann Phys. (1905) 322:549–60. doi: 10.1002/andp.19053220806

CrossRef Full Text | Google Scholar

28. Landau LD, Lifshitz EM. Chapter VI: Diffusion. In: Landau LD, Lifshitz EM, editors. Fluid Mechanics. 2nd edition ed. Pergamon (1987). p. 227–37. Available online at: http://www.sciencedirect.com/science/article/pii/B9780080339337500143

Google Scholar

29. Fulcher GS. Analysis of recent measurements of the viscosity of glasses. J Am Ceram Soc. (1925) 8:339–55. doi: 10.1111/j.1151-2916.1925.tb16731.x

CrossRef Full Text | Google Scholar

30. Vogel DH. Das temperaturabhaengigkeitsgesetz der viskositaet von fluessigkeiten. Phys Zeitsch. (1921) 22:645.

Google Scholar

31. Götze W. Liquids, Freezing and the Glass Transition. Amsterdam: North-Holland (1991).

Google Scholar

32. Gotze W, Sjogren L. Relaxation processes in supercooled liquids. Rep Prog Phys. (1992) 55:241–376. doi: 10.1088/0034-4885/55/3/001

CrossRef Full Text | Google Scholar

33. Kob W, Andersen HC. Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: the van Hove correlation function. Phys Rev E. (1995) 51:4626–41. doi: 10.1103/PhysRevE.51.4626

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Karmakar S, Lerner E, Procaccia I. Direct estimate of the static length-scale accompanying the glass transition. Phys A. (2012) 391:1001–8. doi: 10.1016/j.physa.2011.11.020

CrossRef Full Text | Google Scholar

35. Van Hove L. Correlations in space and time and born approximation scattering in systems of interacting particles. Phys Rev. (1954) 95:249–62. doi: 10.1103/PhysRev.95.249

CrossRef Full Text | Google Scholar

36. Shiba H, Yamada Y, Kawasaki T, Kim K. Unveiling dimensionality dependence of glassy dynamics: 2D infinite fluctuation eclipses inherent structural relaxation. Phys Rev Lett. (2016) 117:245701. doi: 10.1103/PhysRevLett.117.245701

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Yamamoto R, Onuki A. Kinetic heterogeneities in a highly supercooled liquid. J Phys Soc Jpn. (1997) 66:2545–8. doi: 10.1143/JPSJ.66.2545

CrossRef Full Text | Google Scholar

38. Tah I, Sengupta S, Sastry S, Dasgupta C, Karmakar S. Glass transition in supercooled liquids with medium-range crystalline order. Phys Rev Lett. (2018) 121:085703. doi: 10.1103/PhysRevLett.121.085703

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Sastry S, Debenedetti PG, Stillinger FH. Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid. Nature. (1998) 393:554–7. doi: 10.1038/31189

CrossRef Full Text | Google Scholar

40. Ikeda A, Miyazaki K. Slow dynamics of the high density Gaussian core model. J Chem Phys. (2011) 135:054901. doi: 10.1063/1.3615949

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Sengupta S, Karmakar S. Distribution of diffusion constants and Stokes-Einstein violation in supercooled liquids. J Chem Phys. (2014) 140:224505. doi: 10.1063/1.4882066

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Bhowmik BP, Das R, Karmakar S. Understanding the stokes–Einstein relation in supercooled liquids using random pinning. J Stat Mech. (2016) 2016:074003. doi: 10.1088/1742-5468/2016/07/074003

CrossRef Full Text | Google Scholar

43. Lucy LB. An iterative technique for the rectification of observed distributions. Astron J. (1974) 79:745.

PubMed Abstract | Google Scholar

44. Kawasaki T, Onuki A. Slow relaxations and stringlike jump motions in fragile glass-forming liquids: breakdown of the Stokes-Einstein relation. Phys Rev E. (2013) 87:012312. doi: 10.1103/PhysRevE.87.012312

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Kawasaki T, Kim K. Identifying time scales for violation/preservation of Stokes-Einstein relation in supercooled water. Sci Adv. (2017) 3:e1700399. doi: 10.1126/sciadv.1700399

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Kawasaki T, Kim K. Classification of mobile- and immobile-molecule timescales for the Stokes–Einstein and Stokes–Einstein–Debye relations in supercooled water. J Stat Mech. (2019) 2019:084004. doi: 10.1088/1742-5468/ab3114

CrossRef Full Text | Google Scholar

47. Kawasaki T, Kim K, Onuki A. Dynamics in a tetrahedral network glassformer: vibrations, network rearrangements, and diffusion. J Chem Phys. (2014) 140:184502. doi: 10.1063/1.4873346

PubMed Abstract | CrossRef Full Text | Google Scholar