Signature of Dynamical Heterogeneity in Spatial Correlations of Particle Displacement and its Temporal Evolution in Supercooled Liquids

TL;DR

This study uses molecular dynamics simulations to analyze the spatial and temporal evolution of dynamic heterogeneity in supercooled liquids, revealing its persistence across multiple relaxation timescales and its temperature-dependent length scale.

Contribution

The paper provides detailed insights into the temporal evolution and temperature dependence of dynamic heterogeneity and its length scale in supercooled liquids through extensive simulations.

Findings

Dynamic heterogeneity peaks near the alpha-relaxation time.

Heterogeneity exists at the beta-relaxation timescale and persists up to several alpha-relaxation times.

Heterogeneity length scale's temperature dependence remains consistent over time.

Abstract

The existence of heterogeneity in the dynamics of supercooled liquids is believed to be one of the hallmarks of the glass transition. Intense research has been carried out in the past to understand the origin of this heterogeneity in dynamics and a possible length scale associated with it. We have done extensive molecular dynamics simulations of few model glass-forming liquids in three dimensions to understand the temporal evolution of the dynamic heterogeneity and the heterogeneity length scale. We find that although the strength of the dynamic heterogeneity is maximum at a timescale close to characteristic $\alpha$-relaxation time of the system, dynamic heterogeneity itself is well-developed at timescale as short as $\beta$-relaxation time and survives up to a timescale as long as few tens of $\alpha$-relaxation time. Moreover, we discovered that temperature dependence of heterogeneity length remains the same in the whole time window although its absolute value changes over time in a non-monotonic manner.

Equations

$$g_{uu}(r,\Delta t)=\frac{⟨\sum _{i,j=1,j\ne i}^N{\text{u}}_i(t,\Delta t){\text{u}}_j(t,\Delta t)\delta (r-|{\mathbf{r}}_{ij}(t)|)⟩}{4\pi r^2\Delta rN\rho ⟨u(\Delta t){⟩}^2},$$

$${\chi }_4^P(L_B,T)={\chi }_4^P(\infty ,T)f(L_B/{\xi }_d(T)),$$

Full text

glass transition — dynamic length scale

Dynamic heterogeneity (DH) is ubiquitous in a vast variety of natural processes spanning from molecular systems to biological cells and tissues. Existence, characterization and its role in different dynamical processes particularly in the dynamics of glass forming liquid approaching glass transition, is an active field of research. The dynamics of the supercooled liquid become increasingly heterogeneous Berthier and Biroli (2011) as the system approaches the putative glass transition point. Different regions of the system relax at time scales that differ from each other by several orders of magnitude. It is argued that the slowing down of dynamics arises mainly due to the cooperative motion of the particles. The typical size of these cooperatively rearranging regions (CRR) is believed to be of the order of few particles diameter within which motions of particles are spatially and temporarily correlated. It is already shown both in experiments Weeks et al. (2000); Ediger (2000); Gokhale et al. (2016); Cicerone et al. (1995) and in computer simulations Berthier et al. (2011); Yamamoto and Onuki (1997); Keys et al. (2007); Widmer-Cooper et al. (2004) that the mobile particles are non uniformly distributed in the system and tend to form a cluster. This inclination of clustering, as well as the size of the cluster increases as glass transition, is approached.

Extensive studies Berthier and Biroli (2011); Karmakar et al. (2014); Ediger (2000); Berthier et al. (2005) have been performed in past to understand the behaviour of DH at the characteristic long relaxation time scale or the $\alpha$-relaxation time scale, $\tau_{\alpha}$ (defined later) Fragiadakis and Roland (2017) and only a handful studies are done at shorter $\beta$-relaxation time scale Karmakar et al. (2016). In Karmakar et al. (2016), it has been shown that $\beta$-relaxation time, $\tau_{\beta}$ (defined later) has a strong finite size effect, which can be rationalized if one assumes the existence of a growing correlation length. It was surprisingly found that the temperature dependence of this growing correlation length at $\beta$-relaxation time is the same as that of the heterogeneity length scale obtained at $\alpha$-relaxation time. This observation is very surprising as these time scales can differ by many orders of magnitude, especially at low temperatures.

The main goals of this work are two-fold. The first goal is to find signatures of DH in the dynamics at the short time scale of the order of $\tau_{\beta}$. Then we would like to understand the subsequent growth and temporal evolution of DH at timescales ranging from $\tau_{\beta}$ to an order of magnitude larger than $\tau_{\alpha}$. As most of the research works have focused on the characterization of DH in the $\alpha$-relaxation timescale, it is very important to comprehend the time evolution of DH in the intermediate as well as long-time scale compare to $\tau_{\alpha}$, in order to understand the role of DH in glass transition. The main results of this work are the observation of signature of DH in the displacement fields of particles over $\beta$ time scale and the survival of DH at timescales that are larger than $\tau_{\alpha}$ by a factor of $30-40$ in the studied temperature range. We have also discovered that temperature dependence of the heterogeneity length scale remains the same throughout the studied time window, but the region of heterogeneity or the spatial extent of heterogeneity changes with time in a non-monotonic way with its maximum appearing at or near $\tau_{\alpha}$.

Although in Karmakar et al. (2016), it was shown that DH seems to be already developed in the system at timescale close to $\tau_{\beta}$, it was not immediately clear how particle motions at this short time scale get affected due to the presence of the heterogeneity, in other words, whether particle motions at $\tau_{\beta}$ are correlated over the length scale of dynamic heterogeneity, is not immediately clear. Following Karmakar et al. (2016), we define $\tau_{\beta}$ to be the time at which logarithmic derivative of MSD develops a minimum which corresponds to an inflection point in the $\log-\log$ plot of MSD (see SI for further details).

To measure the spatial correlation and to extract the associated length scale in the displacements of particles at $\tau_{\beta}$, we have implemented the procedure given in Ref. Poole et al. (1998); Kob et al. (1997); Donati et al. (1999). Note that this measure of the spatially correlated motion in super-cooled liquids do not depend on arbitrary cutoff parameters as already conclusively shown in Refs.Poole et al. (1998); Kob et al. (1997); Donati et al. (1999) for DH at $\tau_{\alpha}$. The spatial correlation of the particle displacements $g_{uu}(r,\Delta t)$ is defined as

[TABLE]

where $\textbf{u}_{i}(t,\Delta t)=\textbf{r}_{i}(t+\Delta t)-\textbf{r}_{i}(t)$ is the vector displacement of the particle between time $t$ and $t+\Delta t$. $\langle u(\Delta t)\rangle=\langle\frac{1}{N}\sum_{i=1}^{N}|\textbf{u}_{i}(t,\Delta t)|\rangle$, is the average displacement of particles within time interval $\Delta t$. ${\bf r}_{ij}(t)={\bf r}_{j}(t)-{\bf r}_{i}(t)$, is the distance between $i^{th}$ and $j^{th}$ particles. Note that our definition of displacement-displacement correlation is slightly different from the definition given in Ref. (Poole et al., 1998). In Ref. Poole et al. (1998), scalar displacement of the particles are considered, whereas we have considered the vector displacement of the particles (Puosi and Leporini, 2012). It captures the orientational as well as translational correlation in the particle displacements within the time of observation.

We have performed extensive computer simulation of four well-studied model glass formers in three dimensions with different inter particle potential over a wide range of temperatures. The model systems studied are the following (i) 3dKA Kob and Andersen (1995), (ii) 3dR10 Karmakar et al. (2010), (iii) 3dIPL Pedersen et al. (2010) and (iv) 3dHP O’Hern et al. (2002). The details of the models and simulations are given in the Supplemental Information (SI). We find that the growth of DH identified using displacement-displacement correlation function show strong system size dependence at least for at $\tau_{\beta}$. This was not the case in Kim and Saito (2009); Mizuno and Yamamoto (2011) when the correlation function was computed at $\tau_{\alpha}$. Thus we demonstrate the results of $g_{uu}(r,\Delta t)$ for the above models for two different system sizes $N=10000$ and $N=108000$.

In Fig. 1 (left panel) we show the $g_{uu}(r,\Delta t)$ for $N=108000$. It is observed that $g_{uu}(r,\Delta t)$ exhibits damped oscillation which is in agreement with previous numerical Donati et al. (1999); Bennemann et al. (1999) as well as experimental studies Weeks et al. (2007). The correlation function decays to zero exponentially as a function of distance $r$ and with decreasing temperature, the dynamics of the liquid becomes more heterogeneous as the correlation between the particles’ displacements in space extend up to a larger distance as shown in the left panel of Fig. 1. It physically means that particles in the liquids are moving in a co-operative fashion with a monotonically increasing size of the co-operative region as the temperature is lowered. We find a strong system size dependence in $g_{uu}(r,\Delta t)$ as shown in the inset of the left panel of Fig. 1. We have computed the correlation for $N=108000$ and $N=8000$ for the 3dKA model as shown in the inset of Fig. 1 which shows that at low temperature, relative correlation increases with increase in system size. For the robustness of our results, we have calculated the correlation for the other two model systems. We found the results are qualitatively similar (see SI for further details).

In the right panel of Fig. 1, we show a plot of heterogeneity length scale as a function of temperature for the 3dKA model. As expected we observe strong finite size effects on the temperature dependence of DH length obtained from displacement-displacement correlation function. Heterogeneity length obtained from very large system size $(N=108000)$ grow very similarly with dynamical length scale obtained from finite size scale (FSS) of $\tau_{\beta}$ Karmakar et al. (2016). For $N=8000$ system size one observes that the heterogeneity length grows mildly, thus studies on smaller system size would have lead to a conclusion that DH is not present at $\beta$-relaxation time. After conclusively establishing the existence of DH at short time scale, we now focus on the temporal evolution of DH and the associated length scale across the whole range of timescales that can be accessed in simulation.

In Mizuno and Yamamoto (2011), although DH was studied over different time scales but a systematic study on the temperature dependence of the DH length scale was not done. Equipped with the method of block analysis, introduced in Ref.Chakrabarty et al. (2017), a systematic study of the temperature dependence of the dynamical length scale across different time scale for different model glass-forming liquids became computationally feasible. Usually, four-point correlation function, $g_{4}(r,t)$ and the corresponding susceptibilities, $\chi_{4}(t)$ Dasgupta et al. (1991) are used to study DH.$\chi_{4}(t)$ is related to fluctuations in two-point function, $Q(t)$. The Fourier transform of $g_{4}(r,t)$ is known as four-point structure factor $S_{4}(q,t)$ and related to $\chi_{4}(t)$ as $\lim_{q\to 0}S_{4}(q,t)\equiv\chi_{0}(t)$. $\tau_{\alpha}$ is defined as $\langle Q(t=\tau_{\alpha})\rangle=1/e$, where $\langle\ldots\rangle$ denotes ensemble averages. (See SI for further details and definitions). To perform the block analysis, we equilibrate a large system of $N=108000$ particles and measure various quantities of interest by coarse-graining over different block sizes, $L_{B}$. We then obtain the dynamic length scale $\xi_{d}$ by FSS analysis of $\chi_{4}(L_{B},t)$.

In previous studies Chakrabarty et al. (2017); Karmakar et al. (2009), dependence of $\chi_{4}^{P}(L_{B})$, the maximum intensity (peak value) of dynamical susceptibility on block size was studied and the dynamic heterogeneity correlation length $\xi_{d}$ has been estimated by FSS analysis using the following scaling form

[TABLE]

where $\chi_{4}^{P}(\infty,T)$ is the $L_{B}\to\infty$ value of dynamical susceptibility at a temperature $T$. In this work we have done similar scaling analysis for the block size dependence of the intensity of dynamical heterogeneity at few particular time scales $t=\tau_{\alpha}/3,\tau_{\alpha}/2,3\tau_{\alpha}$ to understand how the characteristic length scale and the exponent associated with dynamical susceptibility and the heterogeneity length scale change with time at different temperatures.

In Fig. 2 we have plotted the results for the 3dKA model system. In top panels, we reported the block size dependence of $\chi_{4}(L_{B},T)$ (inset) and the scaling collapse of $\chi_{4}(L_{B},T)$ at $\tau_{\alpha}/3$ (left panel) and $\tau_{\alpha}/2$ (right panels) time scales respectively. A similar analysis also is shown in the bottom panels at time scales $\tau_{\alpha}$ (left panel) and $3\tau_{\alpha}$ (right panel). The scaling observed in all these four cases are indeed very good and the calculated length scales from FSS analysis of the block method is found to be in good agreement with the $\xi_{d}$ obtained from the wave vector dependence of $S_{4}(q,t)$ Lacevic et al. (2003) (discussed in subsequent paragraph). In FSS, $\xi_{d}(T)$ is known up to a multiplicative factor which is the same for all temperatures. In order to fix this uncertainty, $\xi_{d}$ obtain from FSS is scaled to match with $\xi_{d}$ obtained from $S_{4}(q,t)$ at one temperature.

By fitting the q dependence of $S_{4}(q,t)$ for small q values to the Ornstein-Zernike (OZ) form $S_{4}(q,t)\simeq\chi_{0}(t)/[1+(q\xi)^{2}]$, one can also obtain $\xi_{d}$. It has already been shown that heterogeneity length scale obtained from FSS of block method is in good agreement with the same obtained from $S_{4}(q,t)$. In top panels of Fig. 3 we plot the wave vector dependence of the inverse of four point structure factor $S_{4}(q,t)$ for the 3dKA model for two different times $\tau_{\alpha}/3$ (left), $3\tau_{\alpha})$ (right). One can clearly see that OZ form fits the data very well, thus the extracted length scale will be quite accurate. In the bottom panels of the same figure, the temperature dependence of the length scales computed by different methods are compared for $t=\tau_{\alpha}/3$ (left) and $t=3\tau_{\alpha}$ (right). The legend “$\xi_{d}$ from FSS” refers to the length scale obtained from finite size scaling at $t=\tau_{\beta}$ and taken from Karmakar et al. (2009). Note that $\xi_{d}$ from FSS are scaled at $T=0.80$. It is worth highlighting that these results suggests that temperature dependence of $\xi_{d}$ is same across timescales starting from $\tau_{\beta}$ to at least $3\tau_{\alpha}$. To check the robustness of our results, length scales for other two models (3dR10 and 3dIPL) are also computed and the temperature dependence of the length scale are found to be same over the time interval $(\tau_{\beta},3\tau_{\alpha})$ (see SI for further details). Although variation of heterogeneity with decreasing temperature remains same over studied timescales across model systems, the spatial extent of the heterogeneity is observed to increases up to some timescale and then starts to decrease. We are not able to estimate the dynamical length from FSS of block methods at very early time ($t<\tau_{\alpha}/3$) and very long time ($t>3\tau_{\alpha}$) as the intensity of $\chi_{4}(t)$, itself becomes extremely small to do any computation reliably.

Next we examine the power law relation between $\chi_{4}(T)$ and $\xi_{d}(T)$ across these studied timescales. According to Inhomogeneous mode coupling theory (IMCT) Biroli et al. (2006); Berthier et al. (2007a, b) there exists a power relation between $\chi_{3}(\tau_{\alpha})$ (a three point correlator which is believed to be similar to $\chi_{4}$ at least in the scaling behaviour) and $\xi_{d}(\tau_{\alpha})$ as $\chi_{3}(\tau_{\alpha})=\xi_{d}^{2-\eta}(\tau_{\alpha})$ where the theory predicts exponent $2-\eta=4$ in the $\alpha$ regime and $2$ for the $\beta$ regime Stein and Andersen (2008). Following Ref.Chakrabarty et al. (2017) we have done the scaling analysis of $\chi_{4}^{t^{*}}(T)$ to obtain the exponent $2-\eta$ at different times, $t^{*}$. In the large system size limit, $L_{B}>>\xi_{d}$, the $L_{B}$ dependence should disappear in the scaling relation in Eq.2 and it should approach a constant value for $x>>1$. On the other hand for $\xi_{d}\to\infty$ and $L_{B}$ remains finite the dependence of $\chi_{4}$ on $\xi_{d}$ should go away. This implies that scaling function $f(x)$ should be proportional to $x^{2-\eta}$ at $x\to 0$ and $\chi_{4}^{t}(L_{B},T)$ should grow as $L_{B}^{2-\eta}$. We show the results of such an analysis for the 3dKA model in Fig. 4.

The exponent value is found to be $2-\eta\simeq 2$ for both $t=\tau_{\alpha}/3$ as well as at $3\tau_{\alpha}$. Note this is completely different from the exponent ($2-\eta\simeq 4$) predicted by IMCT in the $\alpha$ regime but in very good agreement with the prediction at $\beta$ regime. This observation can be rationalized if one assumes that there will be less activated relaxation at a short time scale as compared to $\alpha$ relaxation time scale and MCT approximation should then be reasonable. Thus we can expect to have a reasonable agreement with MCT predictions at short timescales but not so good at longer timescales. Our results very nicely highlights this agreement with good quality data. We have done similar analysis for other models systems (3dIPL and 3dR10) and found that the exponent $2-\eta$ is very close to 2 (see SI for details).

Next, we look at the time evolution of the DH length scale. In Lacevic et al. (2003), it was shown that time dependence of $\xi_{d}(t)$ is same as $\chi_{4}(t)$, which is in contradiction with the results reported in Toninelli et al. (2005). In Toninelli et al. (2005), $\xi_{d}$ is found to increase monotonically with time, in partial agreement with the results reported in Flenner et al. (2011) for hard sphere systems. Moreover in Flenner et al. (2011) it is found that $\xi_{d}$ saturates to a plateau at a later time. On the other hand, IMCT predicts $\xi_{d}$ to remain constant in between $\tau_{\beta}$ and $\tau_{\alpha}$. We then look at the mutual time evolution of $\xi_{d}$ and $\chi_{4}$ for all these different model systems to understand the apparent contradiction in results reported in the literature.

In top panel of Fig. 5, the time dependence of $\xi_{d}(t)$ for the 3dR10 model is shown. It is clear that $\xi_{d}(t)$ grows up to $\tau_{\alpha}$ and decreases at later time, in agreement with Ref. Lacevic et al. (2003). We then looked at the results obtained for 3dHP model (see SI) O’Hern et al. (2002), which is one of the paradigmatic models in the context of jamming physics. For the 3dHP model, $\xi_{d}(t)$ shows a peak at timescale close to $4\tau_{\alpha}$ (bottom panel of Fig. 5), which implies that $\xi_{d}(t)$ increases even though the overall strength of the heterogeneity decreases after $\tau_{\alpha}$. This suggests that a hard sphere like model are probably different from those models where the particles are treated as point particles. These observations although corroborate the previous observations, the reason for it is not immediately clear. A different choice of the steepness of the potential does change the qualitative results (see SI for detailed discussion). Although for 3dHP model $\xi_{d}(t)\sim\log(t)$ in agreement with Ref. Flenner et al. (2011) for hard sphere systems, one does not see similar dependence for 3dR10 model. The mutual dependence between $\xi_{d}(t)$ and $\chi_{4}(t)$, as $\chi_{4}(t)\sim\xi_{d}(t)^{2-\eta}$, seems to have two different regime - up to $t\sim\tau_{\alpha}$ it is power law like with exponent $2-\eta\sim 4$ for all the models [$4.18$ (3dKA), $3.50$ (3dR10) and $3.39$ for 3dHP model] but a extremely different one with very large exponent for $t>\tau_{\alpha}$ as shown in the insets of Fig. 5. Thus it suggests that one will not be able to extract even approximately the dynamic heterogeneity length $\xi_{d}$ from the measurement of $\chi_{4}(t)$ alone as the exponent $2-\eta$ varies across model systems as well as over the time window of calculation.

Finally, to conclude, we have shown the presence of dynamic heterogeneity in the displacement field of particles at $\tau_{\beta}$ timescale and also highlighted the strong system size effect. Then we showed that although the absolute value of the dynamic heterogeneity length changes with time, the temperature dependence of this length scale across timescale spanning from $\tau_{\beta}$ to couple of times $\tau_{\alpha}$ remain same. This is really surprising and may become an important test for validation of different theories of glass transition. We also found that absolute value of $\xi_{d}$ reaches its maximum value at $t\sim\tau_{\alpha}$ for 3dKA, 3dIPL and 3dR10 models but it does so at $\sim 4\tau_{\alpha}$ for 3dHP models. These observations are indeed extremely intriguing and warrants further research as to why heterogeneity length scale continues to increase at timescales where $\chi_{4}$ already decays to small values for soft sphere and hard sphere model systems.