## Abstract

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.

## INTRODUCTION

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 (*2*–*9*). 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 (*12*–*16*). The SE violation is one of the important anomalies that has been widely reported for supercooled water (*17*–*25*). 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 (*32*–*35*). 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*, *38*–*44*), although this scenario is currently controversial (*45*–*47*). 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.

## RESULTS

### 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 (*49*–*52*). Furthermore, the structural relaxation of supercooled water is identified by the incoherent intermediate scattering function *F*_{s}(*k*,*t*) (see Materials and Methods). The time evolution of *F*_{s}(*k*,*t*) at various temperatures is illustrated in fig. S1C. As outlined in previous simulation studies (*53*–*57*), the behavior of *F*_{s}(*k*,*t*) of supercooled water is characterized by a two-step and nonexponential relaxation below the onset temperature *T*_{A} ≈ 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 *T*_{L} ≈ 220 K. Around this crossover temperature *T*_{L}, 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 *T*_{L} ≈ 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*).

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 *T*_{X} ≈ 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 *T*_{L} ≈ 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 *F*_{s}(*k*,*t*) = exp(−*Dk*^{2}*t*). 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 *T*_{X} ≈ 240 K.

### 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 *N*_{HB}(*t*) was calculated in the time interval *t*, and then the average number fraction *C*_{HB}(*t*) was calculated. The results are shown in fig. S1D [see the detailed definitions of the HB and *C*_{HB}(*t*) in Materials and Methods]. The HB lifetime τ_{HB} was then determined from *C*_{HB}(*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 *G*_{p} exp [−(*t*/τ_{η})^{β}] (see fig. S1B). *G*_{p} and τ_{η} denote the plateau modulus and the stress relaxation time, respectively. The exponent β (< 1) is the degree of nonexponentiality. The temperature dependence of *G*_{p} 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 *G*_{p} 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 *G*_{p}. This relationship will be clarified later. We additionally obtained another preservation of the SE relation, 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 *F*_{s}(*k*,*t*) is described by the Gaussian approximation using the MSD , the relation τ_{α} = (*Dk*^{2})^{−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 is plotted in fig. S1C. We introduce the peak time of Δ*F*_{s}(*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 [see the definition of α_{2}(*t*) in Materials and Methods and fig. S3A]. As demonstrated in Fig. 2C, the time scale is coupled with *D* even at supercooled states. From the definition, the first correction of cumulant expansion of Δ*F*_{s}(*k*,*t*) is given by α_{2}(*t*). Thus, τ_{NG} and exhibit similar temperature dependence. A similar observation has been reported in Lennard-Jones supercooled liquids (*62*); however, the SE ratio 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 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 , during the time interval *t*. We calculate the MSD due to the jp O atoms, 〈Δ*r*_{jp}(*t*)^{2}〉 = (1/*N*)∑_{i∈jp}〈Δ*r*_{i}(*t*)^{2}〉. The summation is over the jp molecule number *N*_{jp}(*t*) at time *t*. In Fig. 3A, the jp component of the MSD, 〈Δ*r*_{jp}(*t*)^{2}〉, is plotted at several temperatures. Because of the jp molecules, this restricted MSD exhibits the diffusive behavior 6*Dt* 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 . In practice, the value of is adjusted to the long-time regimes of the full MSD at each temperature. For the temperature *T* = 190 K, Å is chosen corresponding to the position at the first shoulder of the van Hove function 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, *G*_{s}(*r*,*t*) is largely deviated from the Gaussian form . This deviation implies that the spatial distribution of single-molecular displacement becomes heterogeneous. In particular, a double-peaked structure for *G*_{s}(*r*,*t*) indicates two distinct contributions due to jumping and nonjumping molecules. This non-Gaussianity can be clarified by the decomposition of *G*_{s}(*r*,*t*) due to the number of HBs broken, during the time *t* for the molecule *i* [see the definition of in Materials and Methods]. The molecules having more than three broken HBs [], which destroy the molecules’ local tetrahedral structures, are entirely subjected to the jumping motions. The displacements of these molecules exceed the cutoff length Å at 1 ns. As demonstrated in the study by Kawasaki and Onuki (*63*), this cutoff length 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*) ≡ 〈*N*_{jp}(*t*)〉/*N*, exhibiting activation jumps increases linearly over time. The jump rate is approximately given by τ_{HB}^{−1}, that is, ϕ_{jp}(*t*) ≃ *t*/τ_{HB}, which is demonstrated in Fig. 3B. If the mean jump length is assumed, then the jp component of the MSD 〈Δ*r*_{jp}(*t*)^{2}〉 increases linearly with time as from short time intervals. As demonstrated in Fig. 3A, 〈Δ*r*_{jp}(*t*)^{2}〉 exhibits 6*Dt*. 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 Å.

### 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 Δ*F*_{s}(*k*,*t*) at τ_{NG}, which is represented by . By using the SE preservation *D* ~ τ_{NG}^{−1}, we can express the degree of SE violation in *D*τ_{α} by the temperature dependence of (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 *G*_{s}(*r*,*t*) at lower temperatures is the main feature of dynamic heterogeneities (see fig. S4, A to C). Note that strongly characterizes the contribution of the mobile molecules that move faster than the Gaussian distribution (*65*). The peak time of α_{2}(*t*) becomes smaller than the structural relaxation time τ_{α}, particularly at low temperatures. Up to the time scale , 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 , 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 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 *F*_{s}(*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 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}, , 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 upon supercooling can be interpreted by the viscoelasticity and nonexponentiality in the stress relaxation function *G*_{η}(*t*). The viscosity η is mainly determined by *G*_{p} and τ_{η} according to the long-time behavior of *G*_{η}(*t*) ≃ *G*_{p} exp[−(*t*/τ_{η})^{β}] (see fig. S1B). This dependence of *G*_{η} on *G*_{p} and τ_{η} leads to the approximation of η as , where Γ(⋯) is the gamma function. Figure S2A shows that the plateau modulus *G*_{p} increases, whereas the stretched exponent β decreases with decreasing temperature. The clear correlation between η and *G*_{p}τ_{η}Γ(1/β)/β is demonstrated in fig. S2B except for high temperatures. The plateau moduli are well developed below *T*_{A} ≈ 260 K, which is correlated with the onset of the two-step relaxation in *F*_{s}(*k*,*t*). By combining it with *D* ~ τ_{HB}^{−1}, we obtain the relationship *D*η/*T* ~ [*G*_{p}Γ(1/β)/*T*β] × (τ_{η}/τ_{HB}). The linear relationship between τ_{HB} and τ_{η} provides an alternative representation for the SE violation as *D*η/*T* ~ *G*_{p}Γ(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 . This decoupling is in accord with the development of *G*_{p}, that is, the emergence of solid-like regions. Therefore, the increase in the non-Gaussianity is directly relevant to the increase in *G*_{p} (attained solidity) and to the decrease in β (increase in the nonexponentiality for the stress relaxation), resulting in the SE violation with lowering *T*.

## DISCUSSION

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 *T*_{X} ≈ 240 K is slightly below *T*_{A} ≈ 260 K, which is the onset temperature of the two-step relaxations exhibited in *F*_{s}(*k*,*t*) and *G*_{η}(*t*). These temperatures are above the FSC temperature *T*_{L} ≈ 220 K observed in the temperature dependence of η and τ_{α}. A similar observation, *T*_{L}< *T*_{X} ≲*T*_{A}, 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 Å, 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, 〈Δ*r*_{jp}(*t*)^{2}〉, was characterized by the diffusive behavior 6*Dt*, 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 , τ_{γ}, and τ_{NG}, and four-point dynamic susceptibility ) into the SE violation and preservation. Here, the time scales of τ_{η}, , and τ_{NG} characterize the mobile molecules within dynamic heterogeneities and are coupled with the diffusion constant *D* even for supercooled states. In contrast, τ_{γ} and 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 . Simultaneously, the SE relation is represented by *D*η/*T* ~ *G*_{p}Γ(1/β)/Tβ, where *G*_{p} 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 *G*_{p}) 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 *T*_{A} ≈ 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 *G*_{p}. 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 *G*_{p}. 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 *G*_{p} 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.

## MATERIALS AND METHODS

### Simulations

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 by(1)where **r**_{i}(*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 *F*_{s}(*k*,*t*) with , where *f*_{c}, τ_{s}, τ_{α}, and β_{α} are the fitting parameters. The exponent β_{α} is the degree of nonexponentiality of *F*_{s}(*k*,*t*).

The MSD of the O atom(2)was also calculated. The translational diffusion constant *D* was determined from the long-time behavior of the MSD using the Einstein relation *D* = lim_{t → ∞}〈Δ*r*(*t*)^{2}〉/6*t*.

### HB breakage and its lifetime

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

First, to characterize the local configuration change, we defined the number of HBs broken during time *t* for molecule *i* as . 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 *N*_{HB}(0). At time *t*, the number of remaining HBs, , was less than the initial value *N*_{HB}(0) due to HB breakages (*34*, *35*). The average fraction of HB bonds as a function of time *t* was then defined as(3)

The average HB lifetime τ_{HB} was determined by fitting *C*_{HB}(*t*) with , where the exponent β_{HB} is the degree of nonexponentiality of *C*_{HB}(*t*).

Furthermore, the present scheme is identical to the bond-breakage method applied to various supercooled liquids (*31*, *63*, *73*–*76*). 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 by(4)where *V* is the volume of the system, σ_{αβ} represents the αβ(α,β=*x*,*y*,*z*) components of the off-diagonal stress tensor, and *k*_{B} is the Boltzmann constant. The average stress correlation function is defined as *G*_{η}(*t*) = [*G*_{xy}(*t*) + *G*_{xz}(*t*) + *G*_{yz}(*t*)]/3. The shear viscosity η was determined from the integral of *G*_{η}(*t*) as 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 by(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 as(6)where . The difference is then given by .

This α_{2}(*t*) is mainly dominated by mobile components in the distribution of the single-molecular displacement *G*_{s}(*r*, *t*). To emphasize immobile and slower components, another type of a non-Gaussian parameter is given by(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 *F*_{s}(*k*,*t*)(8)where δ*F*_{i}(**k**,*t*) = cos{**k** ⋅ [**r**_{i}(*t*) − **r**_{i}(0)]} − *F*_{s}(*k*,*t*) is the *i*th molecular fluctuation in the real part of the density correlator (*66*).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/8/e1700399/DC1

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.

## REFERENCES AND NOTES

**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.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).