Spurious violation of the Stokes-Einstein-Debye relation in supercooled water.

The theories of Brownian motion, the Debye rotational diffusion model, and hydrodynamics together provide us with the Stokes-Einstein-Debye (SED) relation between the rotational relaxation time of the [Formula: see text]-th degree Legendre polynomials [Formula: see text], and viscosity divided by temperature, η/T. Experiments on supercooled liquids are frequently performed to measure the SED relations, [Formula: see text]kBT/η and Dt[Formula: see text], where Dt is the translational diffusion constant. However, the SED relations break down, and its molecular origin remains elusive. Here, we assess the validity of the SED relations in TIP4P/2005 supercooled water using molecular dynamics simulations. Specifically, we demonstrate that the higher-order [Formula: see text] values exhibit a temperature dependence similar to that of η/T, whereas the lowest-order [Formula: see text] values are decoupled with η/T, but are coupled with the translational diffusion constant Dt. We reveal that the SED relations are so spurious that they significantly depend on the degree of Legendre polynomials.

Introduction

Characterization of the translational and rotational motions of molecules in liquid states is of great significance[1-3]. For this purpose, various transport properties, such as shear viscosity, translational diffusion constant, and rotational relaxation time have been measured both experimentally, and through molecular dynamics (MD) simulations. These properties also play crucial roles in the understanding of the detailed mechanism of hydrogen-bond network dynamics in liquid water[4-11].

The Stokes–Einstein (SE) relation is one of the important characteristics of the translational diffusion constant, Dt, in many liquid state systems, /, where kB, T, η represent the Boltzmann constant, the temperature, and the shear viscosity, respectively. This SE relation is derived originally from the theories of hydrodynamics and Brownian motion, where a rigid spherical particle with a radius R is assumed to be perfectly suspended in a Stokes flow of a constant shear viscosity η under the stick boundary condition[12]. Thus, R is conventionally regarded as the effective hydrodynamic radius of the molecule when applying the SE relation to molecular liquids[13].

Analogous to translational motion, the rotational Brownian motion leads to another SE relation between the rotational diffusion constant, Dr, and η as /. Based on the Debye model, Dr can also be determined by solving the rotational diffusion equation for the reorientation of the molecular dipole as /, where is the rotational relaxation time of the -th order Legendre polynomials[14]. and are the most-commonly investigated; they are characterized by dielectric relaxation and NMR spectroscopies, respectively. Note that a deviation from / has been reported in supercooled molecular liquids, which is regarded as a sign of the breakdown of the Debye model[15-17]. Those two equations result in the Stokes–Einstein–Debye (SED) relation,

The SED relation can also be expressed asby combining further with the SE relation between Dt and η/T. This SED relation is proportional to the quotient Dt/Dr, which accounts for the coupling between the translational and rotational diffusion dynamics at any temperature.

The violation of the SE relation between Dt and η has been intensively observed in various glass-forming liquids, such as o-terphenyl[18-23]. In particular, the quantity /T increases towards the glass transition temperature, but exhibits a constant value at high temperatures. These experiments indicate that the translational diffusion occurs in a more enhanced manner than estimations using shear viscosity. Many theoretical efforts have therefore been devoted to explaining the violation of the SE relation in glass-forming liquids[24-32]. MD simulations have also been variously performed to address their molecular mechanisms[33-39]. It is commonly argued that the violation of SE relation is a sign of spatially heterogeneous dynamics and of the non-Gaussian property of the particle displacement distribution[22,40].

The validity of the SED relation is still highly controversial, because there are three possible candidates, /η, , and Dt/Dr, that need to be quantified. Recently, the Dr of supercooled molecular liquids has been calculated using MD simulations following the Einstein relation for rotational Brownian motions[41-43]. Experimental analogs have also been reported using optical spectroscopy in colloidal glasses[44,45]. In particular, it has been shown that the temperature dependences of and Dt/Dr are completely different in o-terphenyl liquids[42] and diatomic molecular liquids[43]; Dt/Dr significantly decreases with decreasing temperature, indicating the translational-rotational decoupling. In contrast, exhibits the opposite temperature dependence, i.e., increases in exceed the time scale of the translational diffusion constant, 1/Dt, as the temperature decreases. This discrepancy is thus attributed to the inconsistency between the two expressions, Eqs (1) and (2). However, the direct measurement of R is impractical for molecular liquids both in experiments and MD simulations. More practically, the breakdown of the Debye model, i.e., the dependence of , prevents us from making a precise assessment of the SED relation, whichever one of three quantities is utilized.

For liquid water, it has been widely accepted that the validity of the Debye model for molecular reorientation is limited even in normal states, although that is frequently used when analyzing experimental data. Instead, various large-amplitude rotational jump models have been developed to give an accurate prediction of the rotational relaxation time [7,8,11]. Particularly for supercooled water, the appropriate description for the violation of the SED relation becomes more complicated. Recent MD simulations have demonstrated that the translational and rotational dynamics become spatially heterogeneous upon cooling[46,47]. Furthermore, the violations of the SE and SED relations have been intensely characterized through both experiments and simulations[11,48-67]. In particular, the violation of the SED relation and the translational-rotational decoupling in supercooled water have been reported by calculating Dt and Dr, while has not been calculated[49,52,57], despite the fact that plays an essential role in the precise assessment of the SE relation[60-63]. The SED relation has been investigated by calculating the of SPC/E supercooled water, during which Dr was not calculated[62]. Under these conditions, the SED relation, particularly for the dependence of /η, has not yet been thoroughly investigated, while only one experimental data analysis has been conducted for /η[59].

The purpose of this study is to shed light on the controversy regarding the violation of the SED relation, specifically through the numerical calculations of three quantities, /η, , and Dt/Dr. In particular, we aim to demonstrate that the dependence of /η is an important factor in exploring the inherent translational-rotational dynamics in supercooled water.

Results

Here we examine the translational and rotational SE relations, /T and /T, respectively. We determined Dt and Dr from the long-time behaviors of the translational and rotational mean-square displacements, respectively (see Methods). We calculated η using the Green–Kubo formula for the shear stress correlation function, as detailed in a previous study[63]. Figure 1(a) shows both Dt and Dr as a function of η/T. Comparing these with the dashed line representing the linear relationship, we find that both the translational and the rotational SE relations are invalid in supercooled regimes, particularly at T < 250 K. Note that the rotational SE relation is violated to a greater extent than the translational SE relation. Figure 1(b) shows the ratio of the translational and rotational diffusion constants, Dt/Dt, as a function of the scaled inverse of temperature. The substantial decoupling displayed between the two diffusion constants indicates that the translational and rotational dynamics are decoupling, which is comparable with the previously reported results on ST2[49], SPC/E[52], and TIP4P/2005[57] models. Furthermore, similar results are also demonstrated in o-terphenyl liquids[42] and diatomic molecular liquids[41,43] using MD simulations.

Figure 1

(a) Assessments of the translational and rotational SE relations, /T and /T, made by plotting the relationships between the translational diffusion constant Dt or the rotational diffusion constant Dr and the shear viscosity divided by the temperature η/T. The dashed line represents the linear relation /T, which represents the SE relation. Neither the translational nor the rotational SE relations are satisfied in supercooled region (T < 250 K). (b) Temperature dependence of the ratio of rotational and translational diffusion constants, Dr/Dt. As T decreases, this ratio increases, indicating the translational-rotational diffusion decoupling. (c) Temperature dependence of translational and rotational hydrodynamic radii, Rt and R. Both Rt and Rr decrease significantly upon cooling, accompanied with violation of SE relations. In particular, upon cooling, Rr decreases at a higher rate than that of Rt in response to decreasing temperature.

The observed decoupling of translational-rotational diffusion is directly related to the inconsistency regarding the effective hydrodynamic radius observed when using the SE relations. We quantified the hydrodynamic radius for the translational degree of freedom, Rt = kBT/, and the rotational counterpart, . Figure 1(c) shows the temperature dependences of Rt and Rr. At T = 300 K, Rt and Rr are approximated by 1.2 Å and 1.0 Å, respectively. These values are slightly smaller than the van der Waals radius of the TIP4P/2005 model. As seen in Fig. 1(c), these two radii sharply decrease upon supercooling, accompanied by violation of the translational and rotational SE relations. Moreover, the difference between Rt and Rr increases with decreasing the temperature, implying that the translational and rotational diffusions are decoupling. The relevance of the decoupling Dt/Dr will be discussed below.

Next, we investigate as determined from the -th order Legendre polynomials, and explore its relationship with Dr. Figure 2(a) shows (for , 2, 3, and 6) as a function of the scaled inverse of the temperature. increases for all values as the temperature decreases. Interestingly, we observe that values with higher-order degrees exhibit stronger temperature dependence than those of the lowest order. In other words, the ratios /, /, and / notably decrease as the temperature decreases (see Fig. 2(b)). A similar result was found using MD simulations of the SPC/E supercooled water[62]. As evident in Fig. 2(b), exhibits strong temperature dependence, indicating the breakdown of the Debye model. The observed deviation increases for higher-order values with decreasing temperatures. The breakdown of the Debye model and the inconsistency between Rt and Rr in supercooled states suggest that the SED relations, /η, , and Dt/Dr, are likely spurious quantities. More precisely, these quantities cannot represent real translational and rotational dynamics in supercooled water, regardless of whether they exhibit anomalous deviations from values at high temperatures. We below demonstrate ambiguities of the SED relations, of which results markedly depend on the order .

Figure 2

(a) Temperature dependence of -th order rotational relaxation times for , 2, 3, and 6. The temperature dependences of hydrogen-bond lifetime, , and α-relaxation time, , are included for comparison. (b) Temperature dependence of ratio /, /, and /. Each quantity is scaled by the value at T = 260 K. (c) Assessments of the Debye model, made by plotting the temperature dependence of for , 2, 3, and 6.

Here, we address the SED relation /η, which is the counterpart to recent experimental data[59]. Figure 3(a) shows the relationship between η/T and for and . Note that the results of and are omitted from the plot to improve its clarity. As observed in Fig. 3(a), deviates from the value predicted by the SED relation, particularly at lower temperatures (T < 250 K), instead exhibiting the fractional form . In contrast, with at higher-order follows the SED relation, /T. Figure 3(b) shows the temperature dependence of / (for and ), in comparison with that of the translational SE relation, /T. We observe that the temperature dependence of / is analogous to that of /T, suggesting the violation of the SE relation, whereas / exhibits a weaker temperature dependence. We previously demonstrated the relationship /T in TIP4P/2005 supercooled water[63]. Here, denotes the α-relaxation time that was determined from the incoherent intermediate scattering function F(k, t). The wave-number, k, was chosen as k = 3.0 Å−1, which corresponds to the main peak of the static structure factor of oxygen, SOO(k). This implies that is a good proxy for the translational SE relation /T. Similar results have also been reported for other supercooled liquid systems[34,38,39,68]. Accordingly, the temperature dependence of with higher-order resembles the coupling with that of (see Fig. 2(a)). On the other hand, the deviation of η/ from this value at high temperatures superimposes the violation of the translational SE relation, /T or .

Figure 3

(a) Assessments of the SED relation, /η, made by plotting the relationships between rotational relaxation time for and , and the shear viscosity divided by the temperature, η/T. Both the dotted line and the long-dashed line represent the SED relation, /T. The short-dashed line represents . (b) Assessments of the SED relation, /η, made by plotting the temperature dependence of η/ for and . The violation of the SE relation, /T, is also plotted for comparison. Each quantity is scaled by the value at T = 260 K. The SED ratio η/ exhibits a temperature dependence similar to the violation of the SE relation, /T, whereas the plot of η/ resembles the SED relation, Eq. (1). The background color (white region) indicates the onset temperature of the SE violation.

We next examine the second SED relation, . Figure 4(a and b) display the relationship between Dt and and the temperature dependence of , respectively. As can serve as a proxy of η/T as observed in Fig. 3, exhibits a comparable temperature dependence with the SE ratio /T (see Fig. 4(b)). In contrast, exhibits a temperature dependence similar to that of Dt. Furthermore, we show an alternative quantity, , with the hydrogen-bond lifetime, . Here, represents the time scale characterizing the irreversible hydrogen-bond breakage process, which is determined from the hydrogen-bond correlation function[69-73]. As demonstrated in a previous study and displayed in Fig. 4(b), the SE relation is preserved as if the time scale is replaced with [63]. Similar observations have also been reported in binary soft-sphere supercooled liquids[74] and silica-like network-forming supercooled liquids[68]. This SE preservation is understood by a possible “jump model”. As outlined in the ref.[63], the frequency of the jump motion can be represented as / at investigated temeperatues. Correspondingly, the translational diffusion constant is modeled as /, where denotes a typical jump length (~Å). Therefore, becomes constant. This SE preservation indicates that irreversible hydrogen-bond breakages destroy the local tetrahedral structures, and lead to the translational and rotational molecular jumps with high mobility. It also implies the the coupling between and , which is demonstrated in Fig. 4.

Figure 4

(a) Assessments of the SED relation, , made by plotting the relationship between the translational diffusion constant, Dt, and the rotational relaxation time, , for and . Both the dotted line and the long-dashed line represent the SED relation, . The short-dashed line represents . (b) Assessments of the SED relation, , made by plotting the temperature dependence of the SED ratios . The violation of SE relation, Dt/T, is also plotted for comparison. , with hydrogen-bond lifetime , is also shown. Each quantity is scaled by the value at T = 260 K. Note that shows the preservation of the SE relation[63]. The SED ratio η/ exhibits a temperature dependence similar to the preservation of the SE relation, , whereas bears a certain resemblance to the violation of the SE relation, /T. The background color (white region) indicates of the onset temperature of the SE violation.

To elucidate the molecular mechanism of the demonstrated relationship between Dt and , we analyze the generalized van Hove self correlation function, i.e., the joint probability distribution function for the translational displacement and the rotational angle of the molecule, [75]. Here, Δ and Δ are the translational displacement vector of oxygen and the rotational angle of the dipole moment of a molecule j during a time t, respectively. Figure 5 shows the contour maps of with for , 1 ns, and 10 ns at T = 190 K. For the shorter time interval, t = 0.1 ps, the distribution is stretched towards the rotational angle direction, θ, which is caused by the libration motion of the molecule. This observation corresponds to the oscillations of (see Supplementary Fig. S2). At longer time scales, however, shows the coupling between the translational displacement and the rotational angle, which is consistent with the previously reported results of ref.[75]. Furthermore, the broad ridge separated from the main peak denotes the non-Gaussianity of . A tagged molecule is trapped by a cage composed of neighbor molecules for longer times in supercooled regime. The rotational relaxation time , of which the characteristic angle is π/2 rad, is governed by this large rotational mobility. The single molecule eventually begins diffusion by escaping from the cage, utilizing large translational and rotational mobilities. Thus, is regarded as the time scale coupled with Dt. In addition, the time scale of is similar to the hydrogen-bond lifetime , as demonstrated in Fig. 2(a).

Figure 5

Joint probability distribution functions (generalized van Hove correlation function) for t = 0.1 ns, 1 ns, and 10 ns at T = 190 K. The value of the color bar is normalized by Å2. Remarkable positive correlations between translational displacement and rotational angle Δ are observed, particularly for large r and θ values. This indicates that a large translational motion correlates with a large rotational motion of molecule.

In contrast, the relaxation time corresponds to a molecular reorientation with an angle of 0.37 rad, which lies near to the dominant peak of at t = 10 ns (see Fig. 5(c)). The higher-order mostly highlights immobile molecules both for translational and rotational motions, which will contribute to the dynamical heterogeneities. To investigate this, we characterize the dynamic heterogeneity of translational and rotational motions using the four-point correlation functions, (see Methods and Supplementary Fig. S3(a)). We found that the time scale is akin to the peak time of , which shows that its temperature dependence is similar to that of (see Supplementary Fig. S3(b)). Consequently, the similar temperature dependences between and are demonstrated in Fig. 2(a).

Finally, we discuss the strong decoupling behavior of the translational and rotational diffusion constants, as demonstrated in Fig. 1. As already pointed out in ref.[8], the use of the rotational diffusion constant Dr needs particular care due to the limitation of the angular Brownian motion scnenario. Furthermore, it has been revealed that Dr is superfacial for describing the reorientational motion in supercooled molecular liquids[43]. The angular mean-square displacement is largely influenced by the accumulation of the libration motion, which has a time scale of 0.1 ps. Each molecule can rotate, despite being trapped by the cage, at this short time scale, as indicated in Fig. 5(a). Accordingly, the angular mean-square displacement exhibits a plateau, but its persistent time is much smaller than that of , particularly at lower temperatures (See Supplementary Figs S1 and S2). In contrast, the plateau of after the time scale of libration motion indicates the occurrence of the cage effect (see Supplementary Fig. S2). These findings imply that the decoupling between Dt and Dr has no direct relevance to the real translational-rotational coupling, . This translational-rotational coupling scenario is in accord with the observasion in supercooled molecular liquids[43].

Discussion

In this paper, we report the numerical results of MD simulations of the relationship between the translational and rotational dynamics in TIP4P/2005 supercooled liquid water. Our contributions to the assessment of translational and rotational SE relations and the Debye model can be summarized as follows:

Both translational and rotational SE relations, and , are significantly violated in supercooled states. In particular, the rotational SE relation is violated stronger to a greater extent than that of translational SE relation. Correspondingly, the rotational hydrodynamic radius becomes significantly smaller than translational one with decreasing temperature.

We test the validity of the Debye model, , for the orders , 2, 3, and 6 of the Legendre polynomials, demonstrating that the rotational relaxation time is entirely inconsistent with the rotational diffusion constant Dr.

Furthermore, we systematically examine the SED relations /η and . We reveal that these SED relations strongly depend on the order of , leading to the following spurious argument: The SED relation /T is violated with the lowest-order rotational relaxation time , but is instead satisfied at the higher-order time scale of . In contrast, we find that deviates from values at high temperatures, similarly to the violation of the translational SE relation /T, while superficially satisfies the SED relation.

We observe the coupling between the translational diffusion constant, Dt, and the lowest rotational relaxation time, . We characterize the correlation between large translational and rotational mobilities using from the van Hover correlation function . Furthermore, we find that exhibits the temperature dependence similar to that of the hydrogen-bond lifetime , which is consistent with the previously demonstrated result,  [63].

On the contrary, we show that the higher-order rotational relaxation time is analogous with the α-relaxation time , rather than with . This time scale is characterized by immobile molecular mobilities showing dynamic heterogeneities, which we investigate using the four-point correlation functions and . It is also of essential to examine the role of the length scale of dynamic heterogeneities on the violations of the SE/SED relations in supercooled water. The value is conventionally quantified by the wave-number dependence of the four-point correlation functions[76]. This calculation necessitates MD simulations using more substantial large systems, which are currently undertaken.

In conclusion, in this paper we provide significant and unprecedented insights into the appropriate assessment of SE, Debye, and SED relationships, in doing so clarifying previously awkward and confusing contradictions. Finally, it is worth mentioning the importance of the density dependence on the SE/SED relations in supercooled water. In fact, both and Dt show anomalous density dependence, particularly at low temperaures[65]. Further investigations along this line are necessary to clarify this issue.

Methods

Molecular dynamics simulations

We performed MD simulations of liquid water using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)[77], and used the TIP4P/2005 model to calculate the water molecule potentials[78]. Other MD simulations were also carried out to investigate various properties in supercooled states of this model[65-67,79-86]. The comparison with other rigid non-poralizable models was also made in the recent review[87]. Remark that recent ab initio MD simulations provide a more realistic behavior of the dynamics in supercooled regime[88]. We used a Coulombic cutoff 1 nm. The particle-particle particle-mesh solver was utilized to calculate long-range Coulomb interactions, and the SHAKE algorithm was also used for bond and angle constraints. Periodic boundary conditions were used, and the time step of simulation was 1 fs. First, we employed the NVT ensembles for N = 1,000 water molecules was employed at various temperatures (T = 300, 280, 260, 250, 240, 230, 220, 210, 200, and 190 K) with a fixed mass density of ρ = 1 g cm−3. The corresponding system size was L = 31.04 Å. We conducted the NVE ensemble simulations after the equilibration with a sufficient long time at each temperature, The dynamical quantities including, the α-relaxation time , the translational diffusion constant Dt, and the hydrogen-bond lifetime , and shear viscosity used in here are reported in a previous study[63]. In this study, we newly calculated time correlation functions for characterzing the rotational diffusion constant Dr, the rotational relaxation time of the -th degree Legendre polynomials . Furthermore, the four-point correlation function was also calculated for characterizing dynamic heterogeneities of rotational molecular motions. The trajectories for the calculations of various quantities were for 10 ns (T ≥ 220 K) and 100 ns (T ≤ 210 K). We average over five independent simulation runs for the calculations.

Rotational diffusion constant and rotational relaxation time

We calculated the angular mean-square displacement , following ref.[52] (see Supplementary Fig. S1(a)). We obtained the angular displacement vector of the molecule j is obtained through the time integration of the angular velocity vector, , where the angular velocity vector of the molecule j is given by the cross-product of the normalized polarization vector as, /, with the magnitude . Note that Δt is chosen by a sufficiently small time interval; this was 10 fs in our calculations. We determined the rotational diffusion constant Dr was determined from the long-time limit of as /4t. Furthermore, we independently calculated the angular velocity time correlation function, , where denotes the angular velocity vector of the molecule j in the world reference frame following ref.[57] (see Supplementary Fig. S1(b)). We also used the Green–Kubo formula to obtain Dr as . We confirmed that the Dr values obtained from these two methods are consistent at any temperature.

The rotational correlation function is defined by the autocorrelation function of the normalized polarization vector as , where is the -th order Legendre polynomial as a function of x (see Supplementary Fig. S2). decays from 1 to 0 as t increases. We obtained the -th (, 2, 3, and 6) order rotational relaxation time by fitting to the Kohlrausch–Williams–Watts function .

Rotational four-point correlation functions

We used the four-point correlation function to elucidate the degree of dynamic heterogeneity in supercooled liquids[76]. The four-point correlation for translational motions is defined by the variance of the intermediate scattering function as, , with . We previously calculated using the wave-number k = 3.0 Å−1, and quantified the peak time (note that the same quantity was denoted by in ref.[63]). The rotational four-point correlation function can be analogously defined as , with . The peak time of is represented by .

supplementary information