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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\boldsymbol{\ell }}$$\end{document}ℓ-th degree Legendre polynomials \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\tau }}}_{{\boldsymbol{\ell }}}$$\end{document}τℓ, and viscosity divided by temperature, η/T. Experiments on supercooled liquids are frequently performed to measure the SED relations, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\tau }}}_{{\boldsymbol{\ell }}}$$\end{document}τℓkBT/η and Dt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\tau }}}_{{\boldsymbol{\ell }}}$$\end{document}τℓ, 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\tau }}}_{{\boldsymbol{\ell }}}$$\end{document}τℓ values exhibit a temperature dependence similar to that of η/T, whereas the lowest-order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\tau }}}_{{\boldsymbol{\ell }}}$$\end{document}τℓ 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.

Characterization of the translational and rotational motions of molecules in liquid states is of great significance [1][2][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][5][6][7][8][9][10][11] .
The Stokes-Einstein (SE) relation is one of the important characteristics of the translational diffusion constant, D t , in many liquid state systems, = D k T t B / πηR (6 ), where k B , 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, D r , and η as = 3 . Based on the Debye model, D r 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 . τ 1 and τ 2 are the most-commonly investigated; they are characterized by dielectric relaxation and NMR spectroscopies, respectively. Note that a deviation from τ 1 /τ = 3 2 has been reported in supercooled molecular liquids, which is regarded as a sign of the breakdown of the Debye model [15][16][17] . Those two equations result in the Stokes-Einstein-Debye (SED) relation, t 2 by combining further with the SE relation between D t and η/T. This SED relation is proportional to the quotient D t /D r , which accounts for the coupling between the translational and rotational diffusion dynamics at any temperature. The violation of the SE relation between D t and η has been intensively observed in various glass-forming liquids, such as o-terphenyl [18][19][20][21][22][23] . In particular, the quantity η D t /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][25][26][27][28][29][30][31][32] . MD simulations have also been variously performed to address their molecular mechanisms [33][34][35][36][37][38][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, τ  T /η, τ  D t , and D t /D r , that need to be quantified. Recently, the D r of supercooled molecular liquids has been calculated using MD simulations following the Einstein relation for rotational Brownian motions [41][42][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 τ D t 2 and D t /D r are completely different in o-terphenyl liquids 42 and diatomic molecular liquids 43 ; D t /D r significantly decreases with decreasing temperature, indicating the translational-rotational decoupling. In contrast, τ D t 2 exhibits the opposite temperature dependence, i.e., increases in τ 2 exceed the time scale of the translational diffusion constant, 1/D t , 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 τ 2 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][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67] . In particular, the violation of the SED relation and the translational-rotational decoupling in supercooled water have been reported by calculating D t and D r , 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][61][62][63] . The SED relation has been investigated by calculating the η of SPC/E supercooled water, during which D r was not calculated 62 . Under these conditions, the SED relation, particularly for the  dependence of τ  T /η, has not yet been thoroughly investigated, while only one experimental data analysis has been conducted for τ T 2 /η 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, τ  T /η, τ  D t , and D t /D r . In particular, we aim to demonstrate that the  dependence of τ  T /η is an important factor in exploring the inherent translational-rotational dynamics in supercooled water.

Results
Here we examine the translational and rotational SE relations, η ∝ D t /T and η ∝ D r /T, respectively. We determined D t and D r 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 D t and D r 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, D t /D t , 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.
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, R t = k B T/ πηD (6 ) t , and the rotational counterpart, . Figure 1(c) shows the temperature dependences of R t and R r . At T = 300 K, R t and R r 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 R t and R r increases with decreasing the temperature, implying that the translational and rotational diffusions are decoupling. The relevance of the decoupling D t /D r will be discussed below.
www.nature.com/scientificreports www.nature.com/scientificreports/ Next, we investigate τ  as determined from the -th order Legendre polynomials, and explore its relationship with D r . Figure 2(a) shows τ  (for =  1, 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 τ 1 /τ 2 , τ 1 /τ 3 , and τ 1 /τ 6 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), τ  D r 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 R t and R r in supercooled states suggest that the SED relations, τ  T /η, τ  D t , and D t /D r , 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 .
Here, we address the SED relation τ  T /η, which is the counterpart to recent experimental data 59 . Figure 3(a) shows the relationship between η/T and τ  for =  1 and =  6. Note that the results of =  2 and =  3 are omitted from the plot to improve its clarity. As observed in Fig. 3(a), τ 1 deviates from the value predicted by the SED relation, particularly at lower temperatures (T < 250 K), instead exhibiting the fractional form τ η ∝ − . T ( / ) 1 0 8 . In contrast, τ  with at higher-order =  6 follows the SED relation, τ η ∝ 6 /T. Figure 3(b) shows the temperature dependence of η/ τ  T ( ) (for =  1 and =  6), in comparison with that of the translational SE relation, η D t /T. We observe that the temperature dependence of η/ τT ( ) 1 is analogous to that of η D t /T, suggesting the violation of the /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, D r /D t . As T decreases, this ratio increases, indicating the translational-rotational diffusion decoupling. (c) Temperature dependence of translational and rotational hydrodynamic radii, R t and R. Both R t and R r decrease significantly upon cooling, accompanied with violation of SE relations. In particular, upon cooling, R r decreases at a higher rate than that of R t in response to decreasing temperature. (a) Temperature dependence of -th order rotational relaxation times τ  for =  1, 2, 3, and 6. The temperature dependences of hydrogen-bond lifetime, τ HB , and α-relaxation time, τ α , are included for comparison. (b) Temperature dependence of ratio τ 1 /τ 2 , τ 1 /τ 3 , and τ 1 /τ 6 . Each quantity is scaled by the value at T = 260 K. (c) Assessments of the Debye model, made by plotting the temperature dependence of τ  D r for =  1, 2, 3, and 6.
www.nature.com/scientificreports www.nature.com/scientificreports/ SE relation, whereas η/ τ T ( ) 6 exhibits a weaker temperature dependence. We previously demonstrated the relation- 63 . Here, τ α denotes the α-relaxation time that was determined from the incoherent intermediate scattering function F s (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, S OO (k). This implies that τ α D t is a good proxy for the translational SE relation η D t /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 η/ τT ( ) 1 from this value at high temperatures superimposes the violation of the translational SE relation, η D t /T or τ α D t . We next examine the second SED relation, τ  D t . Figure 4(a and b) display the relationship between D t and τ  and the temperature dependence of τ  D t , respectively. As τ 6 can serve as a proxy of η/T as observed in Fig. 3, τ D t 6 exhibits a comparable temperature dependence with the SE ratio η D t /T (see Fig. 4(b)). In contrast, τ 1 exhibits a temperature dependence similar to that of D t . Furthermore, we show an alternative quantity, τ D t HB , with the hydrogen-bond lifetime, τ HB . Here, τ HB represents the time scale characterizing the irreversible hydrogen-bond breakage process, which is determined from the hydrogen-bond correlation function [69][70][71][72][73] . As demonstrated in a previous study and displayed in Fig. 4(b), the SE relation is preserved as . 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 ∼ f 1/τ HB at investigated temeperatues. Correspondingly, the translational diffusion constant is modeled as /τ HB , where  jump denotes a typical jump length (~Å). Therefore, τ D t HB 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 τ 1 and τ HB , which is demonstrated in Fig. 4.
To elucidate the molecular mechanism of the demonstrated relationship between D t 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, Here, Δ → r t ( ) j and Δθ t ( ) j 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 π θ r G r t 4 ( , ; ) s 2 with = | → | r r for = . t 0 1 ps, 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  C t ( ) (see Supplementary Fig. S2). At longer time scales, however, θ G r t ( , ; ) s 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 θ G r t ( , ; ) s . A tagged molecule is trapped by a cage composed of neighbor molecules for longer times in supercooled regime. The rotational relaxation time τ 1 , 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 www.nature.com/scientificreports www.nature.com/scientificreports/ cage, utilizing large translational and rotational mobilities. Thus, τ 1 is regarded as the time scale coupled with D t . In addition, the time scale of τ 1 is similar to the hydrogen-bond lifetime τ HB , as demonstrated in Fig. 2(a).
In contrast, the relaxation time τ 6 corresponds to a molecular reorientation with an angle of 0.37 rad, which lies near to the dominant peak of θ G r t ( , ; ) s 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, χ t ( ) t,r (see Methods and Supplementary Fig. S3(a)). We found that the time scale τ 6 is akin to the peak time of χ t ( ) t,r , which shows that its temperature dependence is similar to that of τ α (see Supplementary Fig. S3(b)). Consequently, the similar temperature dependences between τ 6 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 D r needs particular care due to the limitation of the angular Brownian motion scnenario. Furthermore, it has been revealed that D r is superfacial for describing the reorientational motion in supercooled molecular liquids 43 . The angular mean-square displacement φ 〈∆ 〉 t ( ) 2 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  C t ( ), particularly at lower temperatures (See Supplementary Figs S1 and S2). In  www.nature.com/scientificreports www.nature.com/scientificreports/ contrast, the plateau of  C t ( ) 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 D t and D r has no direct relevance to the real translational-rotational coupling, τ ∝ − D t 1 1 . 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: η ∝ − , 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.
(ii) We test the validity of the Debye model, , for the orders =  1, 2, 3, and 6 of the Legendre polynomials, demonstrating that the rotational relaxation time τ  is entirely inconsistent with the rotational diffusion constant D r . (iii) Furthermore, we systematically examine the SED relations τ  T /η and τ  D t . 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 τ 1 , but is instead satisfied at the higher-order time scale of τ 6 . In contrast, we find that τ D t 6 deviates from values at high temperatures, similarly to the violation of the translational SE relation η D t /T, while τ D t 1 superficially satisfies the SED relation.
(iv) We observe the coupling between the translational diffusion constant, D t , and the lowest rotational relaxation time, τ 1 . We characterize the correlation between large translational and rotational mobilities using from the van Hover correlation function θ G r t ( , ; ) s . Furthermore, we find that τ 1 exhibits the temperature dependence similar to that of the hydrogen-bond lifetime τ HB , which is consistent with the previously demonstrated result, (v) On the contrary, we show that the higher-order rotational relaxation time τ 6 is analogous with the α-relaxation time τ α , rather than with τ HB . This time scale is characterized by immobile molecular mobilities showing dynamic heterogeneities, which we investigate using the four-point correlation functions χ t ( ) t and χ t ( ) r . It is also of essential to examine the role of the length scale of dynamic heterogeneities ξ 4 on the violations of the SE/SED relations in supercooled water. The ξ 4 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. (vi) 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 D t show anomalous density dependence, particularly at low temperaures 65 . Further investigations along this line are necessary to clarify this issue.

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][66][67][79][80][81][82][83][84][85][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 D t , and the hydrogen-bond lifetime τ HB , 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 D r , the rotational relaxation time of the -th degree Legendre polynomials τ  . Furthermore, the fourpoint 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.
We obtained the angular displacement vector φ ∆ → t ( ) j of the molecule j is obtained through the time integration of the angular velocity vector, . Note that Δt is chosen by a sufficiently small time interval; this was 10 fs in our calculations. We determined the rotational diffusion constant D r was determined from the long-time limit of φ 〈∆ → 〉 t ( ) 2  , where Ω → t ( ) j 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 D r as . We confirmed that the D r values obtained from these two methods are consistent at any temperature.
The rotational correlation function  C t ( ) is defined by the autocorrelation function of the normalized polarization vector → e t ( ) [ ] is the -th order Legendre polynomial as a function of x (see Supplementary Fig. S2).  C t ( ) decays from 1 to 0 as t increases. We obtained the -th ( =  1, 2, 3, and 6) order rotational relaxation time τ  by fitting  C t ( ) 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  . We previously calculated χ t ( ) t using the wave-number k = 3.0 Å −1 , and quantified the peak time τ t (note that the same quantity was denoted by τ χ 4 in ref. 63

Data Availability
The data supporting the findings of this study are available from the corresponding authors upon reasonable request.