# Fast and slow dynamics of hydrogen bonds in liquid water

## Abstract

We study hydrogen-bond dynamics in liquid water at low temperatures using molecular dynamics simulations. We find that bond lifetime (“fast dynamics”) has Arrhenius temperature dependence. We also calculate the bond correlation function and find that the correlation time (“slow dynamics”) shows power-law behavior. This power-law behavior, as well as the decay of the bond correlations, is consistent with the predictions of the mode-coupling theory. The correlation time at the lowest temperature studied shows deviation from power-law behavior that suggests continuity of dynamic functions between the liquid and glassy states of water at low pressure.

Both experiments  and simulations  of water have focused on understanding various aspects of hydrogen bond dynamics, such as the average bond lifetime $`\tau _{HB}`$ and the structural relaxation time $`\tau _R`$. Experiments show that $`\tau _{HB}`$ follows a simple Arrhenius law, but that characteristic relaxation times commonly obey power laws on supercooling . Simulations are particularly useful for investigating water in the supercooled regime since nucleation does not occur on the time scale of the simulations and access to quantitative hydrogen bond information is available. Accordingly, we carry out simulations of normal and supercooled water and find Arrhenius behavior of $`\tau _{HB}`$, as expected from experimental results ; the activation energy associated with $`\tau _{HB}`$ provides a simple test of the criteria used in the simulations to define a bond. Further, we find power-law behavior of $`\tau _R`$ which breaks down at the lowest temperature. This non-trivial growth of $`\tau _R`$ is consistent with a continuous transition of the supercooled liquid to the glassy state.

To elucidate how simple temperature dependence of $`\tau _{HB}`$ \[Fig. 1\] can lead to complex dependence of $`\tau _R`$ \[Fig. 2\], we perform lengthy molecular-dynamics simulations (up to 70 ns, over 7 orders of magnitude) at seven temperatures between 200 K and 350 K . We use the extended simple point charge (SPC/E) potential  for water and calculate the bond dynamics by considering two definitions of an intact hydrogen bond: (i) an energetic definition , which considers two molecules to be bonded if their oxygen-oxygen separation is less than 3.5 Å and their interaction energy is less than a threshold energy $`E_{HB}`$, and (ii) a geometric definition , which uses the same distance criterion but no energetic condition, instead requiring that the O–H…O angle between two molecules must be less than a threshold angle $`\theta _{HB}`$.

We calculate $`\tau _{HB}`$ using both bond definitions over the entire temperature range simulated and find Arrhenius behavior of $`\tau _{HB}`$ \[Fig. 1\]. Measurements of $`\tau _{HB}`$ using depolarized light scattering techniques  find Arrhenius behavior . Previous simulations did not find Arrhenius behavior of $`\tau _{HB}`$ . The activation energy $`E_A`$ associated with $`\tau _{HB}`$ has been interpreted as the energy required to break a bond via librational motion, a “fast” motion . Comparison of experimental and simulated values of $`E_A`$ provides a primitive test of the bonding criteria in our simulations; we obtain reasonable agreement between experimental and simulated values of $`E_A`$ using thresholds of $`E_{HB}=10`$ kJ/mol for the energetic definition and $`\theta _{HB}=30^{}`$  for the geometric definition. We find better quantitative agreement with experiments for $`\tau _{HB}`$ values obtained from the geometric definition than for $`\tau _{HB}`$ values obtained from energetic definition – possibly because the geometric bond definition, like the depolarized light scattering experiments, is highly sensitive to the linearity of the bond. We also calculate $`\tau _{HB}`$ using the thresholds $`E_{HB}=0`$ kJ/mol  and $`\theta _{HB}=35^{}`$ and find Arrhenius behavior, but with $`E_A`$ roughly 30% smaller for the energetic definition, and roughly 10% smaller for the geometric definition.

The quantity $`\tau _{HB}`$ is one characteristic time – the mean – of the distribution of bond lifetimes $`P(t)`$, which measures the probability that an initially bonded pair remains bonded at all times up to time $`t`$, and breaks at $`t`$. $`P(t)`$ is obtained from simulations by building a histogram of the bond lifetimes for each configuration, and may be related to $`\tau _{HB}`$ by $`\tau _{HB}=_0^{\mathrm{}}tP(t)𝑑t`$ . For both bond definitions and all seven temperatures simulated, we calculate the behavior of $`P(t)`$ \[Fig. 3(a)\]. We observe neither power-law nor exponential behavior for either bond definition, unlike previous calculations using the ST2 potential which revealed power-law dependence of $`P(t)`$ . The difference in $`P(t)`$ is not surprising, since the previous study considered a different threshold value and also followed a path near a liquid-liquid critical point . The behavior of $`P(t)`$ for the two bond definitions is different, likely caused by differences in sensitivity of the definitions to librational motion.

We study $`\tau _R`$ (“slow” dynamics) for an initially bonded pair by calculating the bond correlation function $`c(t)`$, the probability that a randomly chosen pair of molecules is bonded at time $`t`$ provided that the pair was bonded at $`t=0`$ (independent of possible breaking in the interim time). To calculate $`c(t)`$, we define

$$c(t)h(0)h(t)/h^2,$$
(1)   
where $`h(t)`$ is a binary function for each pair of molecules $`\{i,j\}`$, and $`h(t)=1`$ if molecules $`\{i,j\}`$ are bonded at time $`t`$ and $`h(t)=0`$ if $`\{i,j\}`$ are not bonded at time $`t`$ . The angular brackets denote an average of all pairs $`\{i,j\}`$ and starting times. We define $`\tau _R`$ by $`c(\tau _R)e^1`$. Short time fluctuations and choice of bond definition do not strongly affect the long-time behavior of the “history-independent” quantity $`c(t)`$ (which does not depend on the continuous presence of a bond), but such fluctuations cause both qualitative and quantitative differences in the long-time behavior of the “history-dependent” $`P(t)`$ (which does depend on the continuous presence of a bond).

We can interpret the behavior of $`\tau _R`$ in terms of the mode-coupling theory (MCT) for a supercooled liquid approaching a glass transition . In accordance with MCT, we find – independent of the two bond definitions considered – power-law growth for $`T210`$ K

$$\tau _R(TT_c)^\gamma ,$$
(2a)   
with $`\gamma =2.7\pm 0.1`$ and $`T_c=197.5\pm 1.0`$ K, approximately 50 K less than the temperature of maximum density $`T_{MD}`$ of SPC/E \[Fig. 2(a)\] . Fits of experimental relaxation times to Eq. (2a) also find $`T_c`$ at a temperature 50 K less than the $`T_{MD}`$ of water . In MCT, $`T_c`$ is the temperature of the ideal kinetic glass transition and is larger than the glass transition temperature $`T_g`$ determined from, e.g., viscosity or relaxation time measurements. At $`T=200`$ K, $`\tau _R`$ is smaller than would be estimated by Eq. (2a), most likely because MCT does not account for activated processes which aid diffusion and reduce relaxation times at low temperatures . Typically, these activated processes become important near $`T_c`$, as observed here.

Our simulation results for $`\tau _R`$ can also be fit by the Vogel-Fulcher-Tammann (VFT) form  for $`T210`$ K

$$\tau _Re^{A/(TT_0)},$$
(2b)   
with $`T_0=160`$ K \[Fig. 2(b)\]. In the entropy theory of the glass transition , $`T_0`$ is associated with the Kauzmann temperature, the temperature where the extrapolated entropy of the supercooled liquid approaches the entropy of the solid. For typical liquids we expect $`T_0<T_g`$, so estimating $`T_0`$ and $`T_c`$ provides lower and upper bounds for $`T_g`$. For water, however, fits to Eq. (2b) of experimental data (which are far above $`T_g`$) yield $`T_0>T_g`$ ; hence we do not consider our $`T_0`$ value to be a lower bound for $`T_g`$ of the SPC/E model. $`T_g`$ is defined experimentally as the temperature where the viscosity reaches $`10^{12}`$ Pa$``$s or $`\tau =100`$ s. Experiments near $`T_g`$ often show a crossover from VFT behavior to “normal” Arrhenius behavior . While our simulations are still relatively far from $`T_g`$ (based on the value of $`\tau _R`$), a naive extrapolation assuming that temperatures $`T=210`$ K and $`200`$ K follow Arrhenius behavior yields $`T_g105`$ K . Applying the same shift to $`T_g`$ as is observed for the $`T_{MD}`$ of SPC/E relative to water (specifically 35 K), the speculated $`T_g`$ of SPC/E is consistent with experimental measurements of $`T_g140`$ K in water .

A continuous crossover to Arrhenius behavior in water might account for the fact that $`T_g<T_0`$ when experimental data are fit to Eq. (2b). VFT or power law behavior (“fragile”) for $`T220`$ K changing to Arrhenius behavior (“strong”) for $`T220`$ K could smoothly connect the structural relaxation times in the liquid with those of the glass. A fragile-to-strong transition in water has been previously suggested , and recent experimental results for the diffusion constant (which scales as $`\tau _R^1`$) at temperatures closer to $`T_g`$ may help to determine if such a transition occurs .

The reactive flux, defined by the derivative

$$k(t)dc(t)/dt,$$
(3a)   
measures the effective decay rate of an initial set of hydrogen bonds. At temperature $`T=300`$ K, non-exponential decay of $`k(t)`$ was found for the closely-related SPC model using the geometric bond definition . Our calculations of $`k(t)`$ \[Fig. 3(b)\] reveal a power-law region $`k(t)t^\zeta `$ for $`T250`$ K for both bond definitions, with an exponent $`\zeta =0.5\pm 0.1`$. The range of the power-law region increases from about one decade at 350 K to about two decades at 250 K.

The value of $`\zeta `$ can be interpreted using MCT, which predicts that $`c(t)`$ decays from a plateau value $`c_p`$ with power-law dependence

$$c_pc(t)t^b$$
(3b)   
in the range where $`k(t)`$ appears to be power-law . From Eq. (3b), $`k(t)t^{b1}`$, so $`b=1\zeta `$. We find $`b=0.5\pm 0.1`$, consistent with previous work , further suggesting that the bond behavior is consistent with MCT predictions for a glass transition. For T $`<250`$ K, the decay of $`k(t)`$ appears to be neither power-law nor exponential.

We note that even though the functional form of $`\tau _{HB}`$ does not appear to be strongly dependent on bond definition, $`P(t)`$ is different for the two definitions—suggesting that $`P(t)`$ may not be the best function for studying bond dynamics. In contrast $`k(t)`$—which includes bond reformation—appears to be largely independent of the bond definition at long times .

The Arrhenius behavior and short time scale of $`\tau _{HB}`$ indicates that librational motion is a thermally activated process that is largely independent of the structural slowing. The non-Arrhenius behavior of $`\tau _R`$ and the non-exponential relaxation of $`k(t)`$ are consistent with the presence of the proposed ideal kinetic glass transition for SPC/E approximately 50 K below the $`T_{MD}`$ of SPC/E water , which coincides with the temperature at which many experimental relaxation times appear to diverge . It was hypothesized that the apparent singular temperature of liquid water observed experimentally may be identified with the $`T_c`$ of MCT ; our results support this hypothesis. They also complement scenarios that account for the anomalous thermodynamic behavior, such as the liquid-liquid transition hypothesis  or the the singularity-free hypothesis .

We thank C.A. Angell, S.V. Buldyrev, S.-H. Chen, I. Große, S.T. Harrington, A. Luzar, O. Mishima, P. Ray, F. Sciortino, A. Skibinsky, D. Stauffer, and especially S. Havlin and S. Sastry for helpful discussions and/or comments on the manuscript. All simulations were performed using the Boston University 192-processor SGI/Cray Origin 2000 supercomputer. FWS is supported by a NSF graduate fellowship. The Center for Polymer Studies is supported by NSF grant CH9728854 and British Petroleum.