Identifying time scales for violation/preservation of Stokes-Einstein relation in supercooled water

See allHide authors and affiliations

Science Advances  18 Aug 2017:
Vol. 3, no. 8, e1700399
DOI: 10.1126/sciadv.1700399


The violation of the Stokes-Einstein (SE) relation D ~ (η/T)−1 between the shear viscosity η and the translational diffusion constant D at temperature T is of great importance for characterizing anomalous dynamics of supercooled water. Determining which time scales play key roles in the SE violation remains elusive without the measurement of η. We provide comprehensive simulation results of the dynamic properties involving η and D in the TIP4P/2005 supercooled water. This enabled the thorough identification of the appropriate time scales for the SE relation Dη/T. In particular, it is demonstrated that the temperature dependence of various time scales associated with structural relaxation, hydrogen bond breakage, stress relaxation, and dynamic heterogeneities can be definitely classified into only two classes. That is, we propose the generalized SE relations that exhibit “violation” or “preservation.” The classification depends on the examined time scales that are coupled or decoupled with the diffusion. On the basis of the classification, we explain the physical origins of the violation in terms of the increase in the plateau modulus and the nonexponentiality of stress relaxation. This implies that the mechanism of SE violation is attributed to the attained solidity upon supercooling, which is in accord with the growth of non-Gaussianity and spatially heterogeneous dynamics.


For simple liquids, the Stokes-Einstein (SE) relation between the shear viscosity η and the translational diffusion constant D is an important characteristic of their transport properties (1). Specifically, this relation implies that D ~ (η/T)−1, where T is the temperature. However, when liquids are supercooled below their melting temperatures, the SE relation is remarkably violated (SE violation), particularly near the glass transition temperature (29). Despite extensive efforts, the origins of the SE violation in supercooled liquids remain elusive.

Generally, transport coefficients such as D and η are mostly coupled at high temperatures. The characteristic time scale is associated with the structural α-relaxation time τα. By contrast, at supercooled states, the SE violation implies that D and η are determined by different time scales. Structural relaxations in supercooled liquids become spatially heterogeneous, which is a different behavior than the homogeneous dynamics observed in normal liquids (6, 10, 11). Thus, the physical implication of SE violation is relevant to the question regarding which time scales determine the transport coefficients in glass-forming liquids. Alternative types of the SE relation D ~ τα or D ~ τα/T have been controversially tested by assuming that τα is proportional to η/T (analogous to the Gaussian approximation) or η (analogous to the Maxwell model), respectively (7, 8).

For liquid water, various anomalies in both its thermodynamics and dynamics have been observed upon supercooling (1216). The SE violation is one of the important anomalies that has been widely reported for supercooled water (1725). In the previous studies on supercooled water, either Dτα or Dτα/T was tested for SE violation. However, the original SE relation Dη/T has not been widely studied because of the high computational costs for calculating η, particularly at low temperatures. Therefore, to determine the origin of the SE violation, obtaining η is important. Hence, the central aims of the present study are to obtain η and to identify the time scales associated with η and D to reveal the origin of the SE violation in supercooled water.

The outline of the present study is as follows. First, the SE violation in supercooled liquid water is examined using molecular dynamics simulations of the TIP4P/2005 model (26, 27). In particular, comprehensive numerical calculations with respect to shear viscosity are performed on the basis of the shear stress correlation function, which are comparable with recent studies for supercooled water using SPC/E (simple point charge/extended) (28) and TIP4P/2005f (29). Our results provide a more systematic examination of the SE violation in supercooled water. The justification of the scenario η/T ~ τα is demonstrated, which is consistent with the previous studies in simple liquids (7, 8, 30, 31).

Second, the role of the time scale associated with hydrogen bond (HB) dynamics in the SE relation is investigated. The rearrangement of the HB network in water is expected to play a critical role in determining its dynamical properties (3235). In addition, the tetrahedrality due to the HB network increases considerably with decreasing temperature (21, 36, 37). This highly structured tetrahedral network is associated with the hypothesized liquid-liquid transition between a high-density liquid and a low-density liquid (14, 3844), although this scenario is currently controversial (4547). Thus, these facts necessitate an investigation of the role of HB dynamics in the SE relation. We show that the SE relation is preserved (SE preservation) when we use the HB breakage time scales instead of τα, that is, the strong coupling between the diffusion constant D and the HB lifetime τHB at any temperature. This preservation is attributed to the activated jumps of mobile molecules that characterize the translational diffusion.

Third, the origin of the observed SE violation [that is, the decoupling between the diffusion constant D(~τHB−1) and the α-relaxation time τα] is elucidated. For this, non-Gaussian parameters and four-point dynamic correlations are examined to probe the degree of dynamic heterogeneities in supercooled water. Here, the SE violation/preservation is additionally demonstrated in terms of other significant time scales, such as the stress relaxation and the mobile/immobile contributions of the dynamic heterogeneities. From these classifications of various time scales, the degree of the SE violation is explained by the increase in the plateau modulus and the nonexponentiality of the stress correlation function upon supercooling. This elucidation for the SE violation is also correlated with the growing of non-Gaussianity and dynamic heterogeneities.


SE violation

The translational mean square displacement (MSD) was calculated at different temperatures (see Materials and Methods). The results are shown in fig. S1A. The diffusion constant was quantified from the long-time behavior of the MSD (see Materials and Methods). In Fig. 1A, we plot the temperature dependence of D. The overall behavior is in good agreement with the previously reported result of the TIP4P/2005 model (48). The shear viscosity η in the TIP4P/2005 supercooled water was investigated from the stress correlation function Gη(t) (see Materials and Methods and fig. S1B). The shear viscosity η was determined from the Green-Kubo formula (see Materials and Methods). The temperature dependence of the viscosity η is plotted in Fig. 1A along with that of D. At T = 300 K, the estimated value is η ≈ 0.78 centipoise (cP), which is approximately the same as the reported value for TIP4P/2005 (4952). Furthermore, the structural relaxation of supercooled water is identified by the incoherent intermediate scattering function Fs(k,t) (see Materials and Methods). The time evolution of Fs(k,t) at various temperatures is illustrated in fig. S1C. As outlined in previous simulation studies (5357), the behavior of Fs(k,t) of supercooled water is characterized by a two-step and nonexponential relaxation below the onset temperature TA ≈ 260 K. Figure 1B shows the temperature dependence of the α-relaxation time τα (see the definition of τα in Materials and Methods). In our calculations, the fragile-to-strong crossover (FSC) weakly occurs at approximately TL ≈ 220 K. Around this crossover temperature TL, the temperature dependence of η and τα changes from non-Arrhenius to Arrhenius behavior, as shown in Fig. 1 (A and B). The FSC is expected as a sign of the compressibility maximum locus (“Widom line”) originating from the liquid-liquid transition (58, 59). The observed TL ≈ 220 K is in accord with the crossing temperature at 1 g cm−3 of the Widom line determined in recent TIP4P/2005 simulations (44, 60, 61).

Fig. 1 Dynamical properties in the TIP4P/2005 supercooled water.

(A) Temperature dependence of viscosity η and translational diffusion constant D. The blue dashed curve is the fitting of the Vogel-Fulcher-Tammann law η ∝ exp [(BT0/(TT0)] with T0 = 170 K and B = 1. 79. The blue dotted line is the Arrhenius law for η ∝ exp (EA/T) at lower temperatures with an activation energy of EA = 52. 1 kJ/mol. Arrhenius behaviors D−1 ∝ exp (EA/T) in both the high and the low temperature ranges are also shown as two red dotted lines, with activation energies of EA = 19. 0 and 38. 6 kJ/mol, respectively. (B) Temperature dependence of the α-relaxation time τα and the HB lifetime τHB. The blue dashed curve is the fitting of the Vogel-Fulcher-Tammann law τα ∝ exp [BT0/(TT0)] with T0 = 175 K and B = 1. 87. The blue dotted line is the Arrhenius law for τα ∝ exp (EA/T) at lower temperatures with an activation energy of EA = 47. 9 kJ/mol. Arrhenius behaviors τHB ∝ exp (EA/T) in both the high and the low temperature ranges are also shown as two red dotted lines, with activation energies EA = 26. 1 and 41. 2 kJ/mol, respectively. (C) Relationship between η/T and τα. The direct proportional relation η/T ∝ τα is obtained. The dashed line is a guide to the eye.

The relationship between η/T and D is presented in Fig. 2A. The SE relation D ~ (η/T)−1 holds at high T but obeys the fractional formula of the SE relation D ~ (η/T)−ζ with ζ ≈ 0.8 below TX ≈ 240 K. The crossover from ζ = 1 to ζ = 0.8 in the fractional SE relation is similar to the recent experimental result (25). This onset temperature appears to be above the FSC TL ≈ 220 K. As noted in Introduction, the alternative expressions for the SE relation are conventionally examined via D ~ τα−1 or D ~ (τα/T)−1. The former formula uses the Gaussian approximation Fs(k,t) = exp(−Dk2t). If τα is characterized by η/T, then Dτα can play the role of the SE relation. Figure 1C shows the proportional relationship η/T ~ τα, which is consistent with the previous results in simple liquids (7, 8, 30, 31). The temperature dependence of Dτα is illustrated together with Dη/T in Fig. 2B. This shows that Dτα is a good indicator of the SE violation Dη/T below its onset temperature TX ≈ 240 K.

Fig. 2 SE violation and preservation.

(A) Translational diffusion constant D versus viscosity scaled by the temperature η/T. Dashed and dotted lines present the fractional SE relation D ∝ (η/T)−ζ. As T decreases, the crossover of the power law exponent from ζ ≈ 1. 0 (satisfying the SE relation) to 0.8 (SE violation) is obtained by the fitting of the data. (B) Inverse temperature dependence of the SE ratios Dη/T, Dτα, Embedded Image, Dτγ, and GpΓ(1/β)/Tβ scaled by their values at TA = 260 K. All SE ratios exhibit the SE violation in the lower-T regime. Note that the data of Embedded Image and GpΓ(1/β)/Tβ above TA = 260 K are omitted. (C) Inverse temperature dependence of the SE ratios DτHB, Dτη, DτNG, and Embedded Image scaled by their values at TA = 260 K. All SE ratios satisfy the SE preservation even at lower T. Note that the data of Dτη above TA = 260 K are omitted.

SE preservation

We introduce the generalized SE ratio Dτ with other significant time scales in supercooled water. First, we focus on the dynamics of HB breakage. The number of the nonbroken HBs for all molecules NHB(t) was calculated in the time interval t, and then the average number fraction CHB(t) was calculated. The results are shown in fig. S1D [see the detailed definitions of the HB and CHB(t) in Materials and Methods]. The HB lifetime τHB was then determined from CHB(t) (see the definition of τHB in Materials and Methods). Its temperature dependence is displayed in Fig. 1B along with that of τα. Remarkably, both D−1 and τHB exhibit a similar Arrhenius temperature dependence, which is different from that of η or τα exhibiting the FSC. Thus, we obtain a marked preservation of the SE relation (SE preservation) D ~ τHB−1 at any temperature, as evident in Fig. 2C. This SE preservation D ~ τHB−1 implies that the appropriate time scale associated with the translational diffusion D is not τα but should instead be the HB lifetime τHB. The HB breakage is commonly speculated to occur intermittently, inducing a markedly large number of jumping water molecules, particularly in supercooled states. Exploring how the HB dynamics are related to the translational diffusion via the jumping molecules is worthwhile. This issue will be discussed later in the paper.

Next, we examine the stress relaxation time τη. The long-time behavior of Gη(t) is well fitted by the stretched exponential function Gp exp [−(tη)β] (see fig. S1B). Gp and τη denote the plateau modulus and the stress relaxation time, respectively. The exponent β (< 1) is the degree of nonexponentiality. The temperature dependence of Gp and β is illustrated in fig. S2A. Note that the stress relaxation time τη differs from the relaxation time of the Maxwell model τM = η/G with the instantaneous shear modulus G = Gη(t = 0). If the temperature dependence of G is negligible, then the viscosity η is identified by τM. Furthermore, provided that τM equals τα, the linear relationship η ~ τα is obtained. However, as seen in fig. S1B, G and Gp increase slightly with decreasing temperature. Instead of the Maxwell model, the viscosity η is determined not only by the stress relaxation time τη but also by the plateau modulus Gp. This relationship will be clarified later. We additionally obtained another preservation of the SE relation, Embedded Image at any temperature, as evident in Fig. 2C. This observation implies that the HB breakage is correlated with the relaxation process of the local stress.

As mentioned in Introduction, the SE violation is possibly attributed to the heterogeneous dynamics, that is, the coexistence of correlated mobile and immobile motions. In this case, the distribution of the single-molecular displacement becomes non-Gaussian at supercooled states. When Fs(k,t) is described by the Gaussian approximation using the MSD Embedded Image, the relation τα = (Dk2)−1 at the diffusive regime is obtained (1). Therefore, the non-Gaussian behavior is directly linked with the SE violation. Analogous to the previous study (54), the degree of the non-Gaussianity Embedded Image is plotted in fig. S1C. We introduce the peak time of ΔFs(k,t) as τNG, which characterizes the time scale of the maximum deviation from the Gaussian behavior. As shown in Fig. 2C, the ratio DτNG represents the SE preservation at any temperature. We also calculated the conventional non-Gaussian parameter α2(t) and determined the peak time of α2(t) as Embedded Image [see the definition of α2(t) in Materials and Methods and fig. S3A]. As demonstrated in Fig. 2C, the time scale Embedded Image is coupled with D even at supercooled states. From the definition, the first correction of cumulant expansion of ΔFs(k,t) is given by α2(t). Thus, τNG and Embedded Image exhibit similar temperature dependence. A similar observation has been reported in Lennard-Jones supercooled liquids (62); however, the SE ratio Embedded Image was used, contrary to our results.

The relationship between the non-Gaussianity and the HB breakage is discussed next. The physical implication of the SE preservation Embedded Image is also given. Furthermore, the effects of characteristic time scales of dynamic heterogeneities on the SE violation/preservation are examined.

Relationship between translational diffusion and HB breakage

Let us examine how HB breakages are coupled with diffusion. To this end, we introduce the jumping (jp) molecules with large displacements. Here, the jp molecules undergoing jumping motions are defined as those O atoms that moved farther than an arbitrary cutoff length Embedded Image, Embedded Image during the time interval t. We calculate the MSD due to the jp O atoms, 〈Δrjp(t)2〉 = (1/N)∑i∈jp〈Δri(t)2〉. The summation is over the jp molecule number Njp(t) at time t. In Fig. 3A, the jp component of the MSD, 〈Δrjp(t)2〉, is plotted at several temperatures. Because of the jp molecules, this restricted MSD exhibits the diffusive behavior 6Dt even at short time regimes (t ≳ 1 ps). After a longer time, the jp contributions to the MSD asymptotically reach the full MSD curves at each temperature because all O atoms eventually move a distance greater than Embedded Image. In practice, the value of Embedded Image is adjusted to the long-time regimes of the full MSD at each temperature. For the temperature T = 190 K, Embedded Image Å is chosen corresponding to the position at the first shoulder of the van Hove function Embedded Image with r = |r|, which represents the distribution of single-molecular displacement (see fig. S4, A to C). At the time scale of τNG ≈ 1 ns, Gs(r,t) is largely deviated from the Gaussian form Embedded Image. This deviation implies that the spatial distribution of single-molecular displacement becomes heterogeneous. In particular, a double-peaked structure for Gs(r,t) indicates two distinct contributions due to jumping and nonjumping molecules. This non-Gaussianity can be clarified by the decomposition of Gs(r,t) due to the number of HBs broken, Embedded Image during the time t for the molecule i [see the definition of Embedded Image in Materials and Methods]. The molecules having more than three broken HBs [Embedded Image], which destroy the molecules’ local tetrahedral structures, are entirely subjected to the jumping motions. The displacements of these molecules exceed the cutoff length Embedded Image Å at 1 ns. As demonstrated in the study by Kawasaki and Onuki (63), this cutoff length Embedded Image enables the selection of irreversible jumps as a result of an activation process analogous to nucleation (64). The average number fraction of the jp molecules, ϕjp(t) ≡ 〈Njp(t)〉/N, exhibiting activation jumps increases linearly over time. The jump rate is approximately given by τHB−1, that is, ϕjp(t) ≃ tHB, which is demonstrated in Fig. 3B. If the mean jump length Embedded Image is assumed, then the jp component of the MSD 〈Δrjp(t)2〉 increases linearly with time as Embedded Image from short time intervals. As demonstrated in Fig. 3A, 〈Δrjp(t)2〉 exhibits 6Dt. Thus, these results clarify the correlation between translational diffusion and HB breakage and agree with the demonstrated SE preservation D ~ τHB−1 (see again Fig. 2C). Furthermore, the mean jump length can be estimated by Embedded Image Å.

Fig. 3 Diffusive properties of jump molecules.

(A) Translational MSD 〈Δr(t)2〉 for the O atom (solid curves), the MSD due to the jp O atoms 〈Δrjp(t)2〉 (points), and the Einstein relation 6Dt (black dashed lines). Here, D is determined by the long-time asymptotic value of 〈Δr(t)2〉/6t at each temperature. For temperatures T = 190, 200, 220, 240, and 260 K, Embedded Image is adjusted to 1.9, 1.8, 1.7, 1.5, and 1.4 Å, respectively. (B) Average number fraction of the jp molecules ϕjp(t). Dashed lines represent the linear growth relations tHB at each temperature.

Mechanism of SE violation

The demonstrated SE violation indicates that the translational diffusion constant D is not characterized by the α-relaxation time τα. The SE violation is explained in terms of the peak height of ΔFs(k,t) at τNG, which is represented by Embedded Image. By using the SE preservation D ~ τNG−1, we can express the degree of SE violation in Dτα by the temperature dependence of Embedded Image (see text S1 for details). That is, the increase in the degree of non-Gaussianity is in accord with the degree of the SE violation in Dτα.

As mentioned above, the non-Gaussianity is directly relevant with dynamic heterogeneities. The observed double-peak structure of Gs(r,t) at lower temperatures is the main feature of dynamic heterogeneities (see fig. S4, A to C). Note that Embedded Image strongly characterizes the contribution of the mobile molecules that move faster than the Gaussian distribution (65). The peak time Embedded Image of α2(t) becomes smaller than the structural relaxation time τα, particularly at low temperatures. Up to the time scale Embedded Image, a tagged molecule is trapped by the surrounding cage, which is observed as the plateau of MSD (see fig. S1A). The cage eventually breaks at Embedded Image, and then the tagged molecule begins to escape from the original position due to the jump motion. This physical implication is consistent with the demonstrated SE preservation DτHB.

The non-Gaussianity is additionally quantified by a new non-Gaussian parameter γ(t), which emphasizes the immobile and slower contribution of dynamic heterogeneities [see the definition of γ(t) in Materials and Methods and fig. S3B] (65). The peak time τγ of γ(t) becomes slower than Embedded Image with decreasing temperature. This indicates the decoupling between mobile and immobile molecules in supercooled states. As demonstrated in Fig. 2B, the SE ratio Dτγ exhibits the SE violation, following the similar temperature dependence of Dτα. Another quantity to examine the dynamic heterogeneities is the four-point correlation function χ4(k,t) that is defined by the variance of Fs(k,t) [see the definition of χ4(k,t) in Materials and Methods and fig. S3C] (66). The value of χ4(k,t) is related to the correlation length of dynamic heterogeneities at the time scale t. As demonstrated in fig. S3C, χ4(k,t) exhibits the peak value at τα, which increases as the temperature decreases. Figure 2B shows that the peak time Embedded Image of χ4(k,t) also acts as the SE violation. These results indicate that the immobile and slower component of non-Gaussianity is characterized by the time scales τα and τγ presenting the SE violation. In contrast, the time scales τNG, Embedded Image, and τHB are coupled with the diffusion constant D, which is markedly governed by the mobile and jumping molecules.

Furthermore, the increase in the degree of the non-Gaussianity Embedded Image upon supercooling can be interpreted by the viscoelasticity and nonexponentiality in the stress relaxation function Gη(t). The viscosity η is mainly determined by Gp and τη according to the long-time behavior of Gη(t) ≃ Gp exp[−(tη)β] (see fig. S1B). This dependence of Gη on Gp and τη leads to the approximation of η as Embedded Image, where Γ(⋯) is the gamma function. Figure S2A shows that the plateau modulus Gp increases, whereas the stretched exponent β decreases with decreasing temperature. The clear correlation between η and GpτηΓ(1/β)/β is demonstrated in fig. S2B except for high temperatures. The plateau moduli are well developed below TA ≈ 260 K, which is correlated with the onset of the two-step relaxation in Fs(k,t). By combining it with D ~ τHB−1, we obtain the relationship Dη/T ~ [GpΓ(1/β)/Tβ] × (τηHB). The linear relationship between τHB and τη provides an alternative representation for the SE violation as Dη/T ~ GpΓ(1/β)/Tβ, as demonstrated in Fig. 2B. Additionally, the SE violation is attributed to the immobile molecules within dynamic heterogeneities, whose time scales are τα, τγ, and Embedded Image. This decoupling is in accord with the development of Gp, that is, the emergence of solid-like regions. Therefore, the increase in the non-Gaussianity Embedded Image is directly relevant to the increase in Gp (attained solidity) and to the decrease in β (increase in the nonexponentiality for the stress relaxation), resulting in the SE violation with lowering T.


In summary, we reported comprehensive numerical results concerning the SE relation in the TIP4P/2005 supercooled water. In particular, the temperature dependence of the shear viscosity was quantified from the stress correlation function in a wide temperature range (190 to 300 K). Thus, the SE relation in supercooled liquid water was systematically examined as follows.

We reported that the violation of the SE relation is characterized by the fractional form D ~ (η/T)−ζ with ζ ≈ 0.8. The onset temperature of the SE violation TX ≈ 240 K is slightly below TA ≈ 260 K, which is the onset temperature of the two-step relaxations exhibited in Fs(k,t) and Gη(t). These temperatures are above the FSC temperature TL ≈ 220 K observed in the temperature dependence of η and τα. A similar observation, TL< TXTA, has been reported in numerical results using ST2 water model (40). Furthermore, the degree of the SE violation was identified by Dτα from the proportional relation η/T ~ τα. We also explored the role of HB breakage on the SE relation. The results revealed that the time scale associated with the translational diffusion constant D should be the HB lifetime τHB, in accordance with the preservation of the SE relation D ~ τHB−1 even for supercooled states. We observed that both D and τHB exhibit an Arrhenius temperature dependence with a similar activation energy. This SE preservation proposes the temperature-independent length scale Embedded Image Å, which has no relation with the Widom line and the possible liquid-liquid transition.

We quantitatively confirmed that the observed preservation of the SE relation D ~ τHB−1 was attributed to the effect of the activated jumping of mobile molecules on the translational diffusion. The distinction between jumping and nonjumping molecules in supercooled states is a manifestation of spatially heterogeneous dynamics, that is, the dynamic heterogeneities in supercooled water (67, 68). In particular, the MSD from the jp molecules, 〈Δrjp(t)2〉, was characterized by the diffusive behavior 6Dt, even on short time scales. The jumping rate was characterized by the inverse of the HB lifetime τHB. An analogous result showing the SE preservation between D and τHB has already been obtained in both binary soft-sphere mixtures (fragile liquids) (63) and silica-like network-forming liquids (strong liquids) (31). In these studies, the bond-breakage method characterizing the changes in local particle connectivity was used, which is essentially the same as the current analysis regarding the HB network in liquid water.

Furthermore, we categorized other time scales (such as stress relaxation time τη, time scales of the non-Gaussianity Embedded Image, τγ, and τNG, and four-point dynamic susceptibility Embedded Image) into the SE violation and preservation. Here, the time scales of τη, Embedded Image, and τNG characterize the mobile molecules within dynamic heterogeneities and are coupled with the diffusion constant D even for supercooled states. In contrast, τγ and Embedded Image exhibit the temperature dependence similar to that of the α-relaxation time τα. These time scales are governed by the immobile and slower molecules and are decoupled with D when the temperature decreases, leading to the SE violation.

Finally, we revealed that the SE violation was attributed to the increase in the degree of the non-Gaussianity Embedded Image. Simultaneously, the SE relation is represented by Dη/T ~ GpΓ(1/β)/Tβ, where Gp and β denote the plateau modulus and the stretched exponent in the stress relaxation function, respectively. Here, the proportional relationship between the stress relaxation time and the HB lifetime τη ~ τHB was used. Therefore, the time scales supporting the violation or preservation of the SE relation were thoroughly identified; attained solidity (increasing Gp) and increasing nonexponentiality (decreasing β) give rise to the SE violation with decreasing temperature. Note that the nonexponentiality in the stress relaxation is also a significant hallmark of the dynamic heterogeneities (69). In our simulations, the plateau modulus and the nonexponentiality develop largely below TA ≈ 260 K. Correspondingly, the growths of the non-Gaussianity and the dynamic susceptibility are noticeable, as demonstrated in fig. S3 (A to C).

There are other implications in developing the plateau modulus Gp. The SE violation with decreasing temperature will be relevant with the decoupling between translational and rotational motions in supercooled water. It is expected that translational relaxations become slower, whereas molecules undergo rotational motions even inside immobile solid-like regions (20). The mechanism of this decoupling will be clarified in terms of the attained solidity Gp. In addition, a recent theoretical study has shown that the spatially heterogeneous dynamics is attributed to the thermal excitation between the different metabasins of the free energy landscape (70). In the framework, the value of the plateau modulus Gp is determined by the curvature of the local metabasin. Considering these investigations, the demonstrated SE preservation DτHB will provide deeper insight into the activated jump events occurring between different metabasins, not only in supercooled water but also in various glassy systems, although further investigations are required to confirm it.



The molecular dynamics simulations of liquid water were performed using the LAMMPS package (71). The TIP4P/2005 model was used for the water molecules (26, 27). The NVT ensemble for N = 1000 water molecules was first simulated at various temperatures (T = 300, 280, 260, 250, 240, 230, 220, 210, 200, and 190 K) with a fixed density ρ = 1 g cm−3. The corresponding linear dimension of the system is L = 31.04 Å. After equilibration for a sufficient time at each temperature, the NVE ensemble simulations were completed, yielding five independent 100-ns trajectories from which the various physical quantities were calculated. The simulations were performed with a time step of 1 fs. The total CPU (central processing unit) time approximated about 20 years of single core time.

Incoherent intermediate scattering function and MSD

The incoherent intermediate scattering function is given byEmbedded Image(1)where ri(t) is the position vector of the O atom of the water molecule i at time t. The bracket indicates an average over the initial time t = 0. The wave number k = |k| was chosen as k = 3.0 Å̊−1, which corresponds to the first peak position of the static structure factors of the O atom. The α-relaxation time τα was determined by the fitting Fs(k,t) with Embedded Image, where fc, τs, τα, and βα are the fitting parameters. The exponent βα is the degree of nonexponentiality of Fs(k,t).

The MSD of the O atomEmbedded Image(2)was also calculated. The translational diffusion constant D was determined from the long-time behavior of the MSD using the Einstein relation D = limt → ∞〈Δr(t)2〉/6t.

HB breakage and its lifetime

The dynamics of HB was investigated by using r definition (72), where only the intermolecular O–H distance rOH is involved. An HB bond is present at the initial time if the rOH is less than 2.4 Å, corresponding to the first minimum of the radial distribution function gOH(r). At a later time t, the HB is broken when the distance rOH becomes larger than 2.4 Å̊, which was determined from the second minimum position of gOH(r).

First, to characterize the local configuration change, we defined the number of HBs broken during time t for molecule i as Embedded Image. Next, the characteristic time scale (that is, the HB lifetime τHB) was determined. The number of HBs was calculated at the initial time t = 0 and denoted as NHB(0). At time t, the number of remaining HBs, Embedded Image, was less than the initial value NHB(0) due to HB breakages (34, 35). The average fraction of HB bonds as a function of time t was then defined asEmbedded Image(3)

The average HB lifetime τHB was determined by fitting CHB(t) with Embedded Image, where the exponent βHB is the degree of nonexponentiality of CHB(t).

Furthermore, the present scheme is identical to the bond-breakage method applied to various supercooled liquids (31, 63, 7376). These previous studies have demonstrated that the bond-breakage method is more remarkable when the collective motions and dynamic heterogeneities peculiar to supercooled states are characterized.

Stress correlation function and shear viscosity

The autocorrelation function of the off-diagonal stress tensor is given byEmbedded Image(4)where V is the volume of the system, σαβ represents the αβ(α,β=x,y,z) components of the off-diagonal stress tensor, and kB is the Boltzmann constant. The average stress correlation function is defined as Gη(t) = [Gxy(t) + Gxz(t) + Gyz(t)]/3. The shear viscosity η was determined from the integral of Gη(t) as Embedded Image using the Green-Kubo formula.

Characterizations of dynamic heterogeneities

The non-Gaussian parameter for the displacements of the molecules is the conventional quantity to characterize dynamic heterogeneities in various glass-forming liquids. The equation is given byEmbedded Image(5)that represents the degree of the deviation from the Gaussian approximation in the density correlation function, which is revealed by the cumulant expansion such asEmbedded Image(6)where Embedded Image. The difference is then given by Embedded Image.

This α2(t) is mainly dominated by mobile components in the distribution of the single-molecular displacement Gs(r, t). To emphasize immobile and slower components, another type of a non-Gaussian parameter is given byEmbedded Image(7)which is referred to as the new non-Gaussian parameter (65).

Furthermore, the four-point dynamic susceptibility χ4(k,t) was used to identify the magnitude of dynamic heterogeneities. The equation is defined from the variance of Fs(k,t)Embedded Image(8)where δFi(k,t) = cos{k ⋅ [ri(t) − ri(0)]} − Fs(k,t) is the ith molecular fluctuation in the real part of the density correlator (66).


Supplementary material for this article is available at

fig. S1. Time correlation functions.

fig. S2. Properties of stress-relaxation function.

fig. S3. Characterizations of dynamic heterogeneities.

fig. S4. van Hove correlation functions.

text S1. SE violation evaluated by non-Gaussianity.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank A. Onuki for fruitful discussions. K.K. is grateful to T. Nakamura for helpful information on computation. The numerical calculations were performed at the Research Center of Computational Science, Okazaki, Japan. Funding: This work was partially supported by the Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research (KAKENHI) (grants nos. JP15H06263 and JP16H06018 to T.K. and JP26400428 and JP16H00829 to K.K.). Author contributions: T.K. and K.K. designed the research, analyzed the data, and wrote the paper. K.K. performed the simulation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article