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 timescales: 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 $\alpha$ 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 relationships between these different timescales and their connection to dynamic heterogeneity in the system.

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 times cales: 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 times cales and their connection to dynamic heterogeneity in the system.

I. 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 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][2][3][4][5][6][7][8][9][10][11][12]. It is now well-accepted that the slow dynamics becomes increasingly heterogenious near the glass transition. Numerous studies have been performed to understand this behavior [12][13][14]. While the origin of dynamic heterogeneity is still not fully understood, the lifetime of these heterogeneities is also debated. 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, this time scale has been shown to be comparable to the αrelaxation time [16][17][18][19][20].
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 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 survive up to a time scale that can be much larger than the α-relaxation time. The difference between the survival time of dynamic heterogeneity and * Electronic address: rajsekhard@tifrh.res.in † Electronic address: cdgupta@iisc.ac.in ‡ Electronic address: smarajit@tifrh.res.in the α-relaxation time increases with increasing supercooling. In earlier studies [23][24][25][26], it has been 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.
It is well-known that the mean squared displacement (MSD), ∆r 2 , computed for a supercooled liquid shows an initial ballistic growth with time and then tends to a 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 neighbours. 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 behaviour at this time scale, the dynamics of the system remain spatially heterogeneous. In Ref. [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. Thus there exists two different time scales, related to diffusion. One of these is the time at which the slope of the MSD versus 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 Ref. [15], it has been suggested that the time scale obtained via 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 re-lation between the diffusion coefficient and the viscosity [27][28][29], 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 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 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 [30,31](VFT) or mode coupling [32,33] (power law) fits to the temperature dependence of these time scales 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 the Model & Method section. 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.

II. MODELS & METHODS
We have studied two different model glass forming liquids in three dimensions.
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 = 8000.
• 3dR10 Model [35] : It is a 50 : 50 binary mixture where the particles interacts via 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 = 10000.
For all the studied models we have performed NVT molecular dynamics (MD) simulations with a modified leap-frog algorithm. We use the Berendsen thermostat to keep the temperature constant. Our integration time step is dt = 0.005. The number of MD steps are 10 7 − 10 8 , depending on the temperature.

A. 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 where w(x) is a window function with 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.

B. 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 [36] is defined as where x i (t) is the x-component of the position vector of particle i at time t.

C. Bond-Breakage Correlation Function
The bond-breakage correlation function is defined in the following way [37][38][39]. At t = t 0 a pair of particle i (of type α) and j (of type β) is considered to be bonded if If r ij (t) ≤ 1.35σ αβ , the bond is said to have survived at time t. We calculate 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.

A. 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 2 ∼ t 2 , then a plateau followed by a diffusive regime i.e. ∆r 2 ∼ 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 Fig. 1 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 Fig. 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 [40] 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.
We have also calculated the bond-breakage time scale τ BB following the method described in Sec. III. It is the time scale where the 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 bottom panel of Fig. 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 at remperatures lower than the onset temperature. However, we do not have a clear understanding of this behavior.
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 determine accurately 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. Fig. 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.

B. Time Scale from the Distribution of Single-particle Displacements
Next, we estimate time at 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 10 (δr), t) at time t. As discussed in [15], this distribution is connected to the self part of the van Hove correlation function as P (log 10 (δr), t) = ln 10 4πδr 3 G s (δr, t) .
If the motion of the particles is governed by Fick's law of diffusion, the distribution of single particle displacements i.e., P (log 10 (δ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,41]. Fig. 4 shows plots for P (log 10 (δ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.
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 is 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 Ref. [15] and refer to this time scale as τ F . In the left panel of Fig. 5, we show plots of the peak value of P (log 10 (δr), t) with increasing time for several 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 V F T 0.289. To see the mutual relationship between τ α and τ F , in Fig. 6 we have plotted the ratio (τ F /τ α ) for all the stud-ied temperatures. Form 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 Fig. 6. 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 Fig 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 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 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 Fig. 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].
As suggested in Ref. [15], the time scales τ F and τ H may be thought of as lower bounds of the lifetime of dynamic heterogeneity. As reported in Refs. [42,43], 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 [44]. We have followed the method described in [42,43]. In the top and bottom panels of Fig. 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].

V. RELATION BETWEEN DIFFERENT TIME SCALES
In the previous section, we found that the temperature dependence of the different relaxation times (τ α , τ D , τ F and τ H ) can be well-fitted by the VFT form with very similar values of the divergence temperature T V F T . This suggests that the different time scales are related to one another via power laws. In Fig.10 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 powerlaw relation between one of the time scales of diffusion and the α-relaxation time is 0.83 for τ D , implying that τ D increases more slowly that τ α 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 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 heterogenieties 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.

VI. STOKES-EINSTEIN RELATION AND KINETIC FRAGILITY
The violation of the Stokes-Einstien [27][28][29]relation between the diffusion coefficient D and the viscosity (or a suitable time scale τ ) in deeply supercooled liquids [42,43] 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 [? ] that D ∝ 1/τ α (i.e. the SE relation is satisfied) at high temperatures, but a modified relation, D ∝ τ −1+ω α , 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 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 versus different time scales shown in Fig 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 Fig 11. 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 ) = τ 0 exp[1/(κ(T /T V F T − 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.

VII. SYSTEM-SIZE DEPENDENCE OF TIME SCALES
In Fig. 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 − 10000]. We have done simulations for smaller system sizes in the range N ∈ [100−1000] 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 Fig. 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 [45]. This observation is puzzling and warrants further investigation.

VIII. 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 . This result raises questions that merit further investigation.
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 diverges 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 powerlaw relation also implies that the fractional SE exponent and the kinetic fragility associated with the diffusionrelated 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 , implying that the time at which the MSD begins to be linear in time and the value of the diffusion coefficient obtained from this linear dependence are strongly related to each other (D ∝ 1/τ D ).
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 (Fig.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.
All these observations, raises several interesting questions about characteristics of the dynamics of supercooled liquids. Answers to these questions will be important in understanding the role of dynamic heterogeneity in glassy relaxation.