Molecular Dynamics Study of Orientational Cooperativity in Water

Recent experiments on liquid water show collective dipole orientation fluctuations dramatically slower then expected (with relaxation time $>$ 50 ns) [D. P. Shelton, Phys. Rev. B {\bf 72}, 020201(R) (2005)]. Molecular dynamics simulations of SPC/E water show large vortex-like structure of dipole field at ambient conditions surviving over 300 ps [J. Higo at al. PNAS, {\bf 98} 5961 (2001)]. Both results disagree with previous results on water dipoles in similar conditions, for which autocorrelation times are a few ps. Motivated by these recent results, we study the water dipole reorientation using molecular dynamics simulations in bulk SPC/E water for temperatures ranging from ambient 300 K down to the deep supercooled region of the phase diagram at 210 K. First, we calculate the dipole autocorrelation function and find that our simulations are well-described by a stretched exponential decay, from which we calculate the {\it orientational autocorrelation time} $\tau_{a}$. Second, we define a second characteristic time, namely the time required for the randomization of molecular dipole orientation, the {\it self-dipole randomization time} $\tau_{r}$, which is an upper limit on $\tau_{a}$; we find that $\tau_{r}\approx 5 \tau_{a}$. Third, to check if there are correlated domains of dipoles in water which have large relaxation times compared to the individual dipoles, we calculate the randomization time $\tau_{\rm box}$ of the site-dipole field, the net dipole moment formed by a set of molecules belonging to a box of edge $L_{\rm box}$. We find that the {\it site-dipole randomization time} $\tau_{\rm box}\approx 2.5 \tau_{a}$ for $L_{\rm box}\approx 3$\AA, i.e. it is shorter than the same quantity calculated for the self-dipole. Finally, we find that the orientational correlation length is short even at low $T$.

When water is cooled, the cooperativity of water molecules increases. Recent experiments on water show large correlated domains of dipoles at ambient conditions which have a relaxation time much larger than the autocorrelation time of individual dipoles [21]. MD studies of water models also show the possibility of formation of large correlated domains of dipoles in bulk as well as interfacial water [35] (where these correlated patterns of dipoles are pinned to solvated amino acids). These two studies are the principal motivation for the present investigation of the rotational cooperativity of water molecules.
A challenging problem is to develop methods of describing molecular motion in water that are better able to interpret experimental results, such as incoherent quasielastic neutron scattering, light scattering, dielectric, and nuclear magnetic resonance experiments [2,18]. Several approximation proposals have been made for various autocorrelation functions describing both rotational and translational motion [20,27]. These methods usually assume the Kohlrausch-Williams-Watts stretched exponential for the long time relaxation behavior of autocorrelation functions φ(t), as predicted by mode coupling theory (MCT), The relaxation time τ a , the exponent β, and the non-ergodicity factor A are fitting parameters that depend on temperature T and density ρ [22,23,24,25,27,28,29,30,31].
Our interest here is to study the orientational dynamics of water by simulating SPC/E water. First we calculate the orientational autocorrelation time as the fitting parameter τ a appearing in Eq. (1) [22,23]. Other definitions are possible, e.g., based on other fitting functions for the orientational autocorrelation function decay, such as the biexponential [26,36] or the von Schweidler law [33]. In all cases, the orientational autocorrelation times are the result of multi-parameter fitting procedures and roughly correspond to the characteristic time over which the orientational autocorrelation function decays by a factor of e ≈ 2.7.
To find an upper limit of the orientational autocorrelation time τ a , we will introduce a new quantity, the dipole randomization time τ r , as the time after which the fluctuations of the dipoles resemble an uncorrelated random variable [37] (Sec. IV A). We find τ r > τ a , and that τ r and τ a are linearly related (Sec. IV B), which is consistent with the MCT predictions that: (i) The autocorrelation times of all the autocorrelation functions of any fluctuation coupled to density fluctuations diverge at the same temperature T MCT with the same power law exponent; (ii) All the characteristic times of a supercooled liquid are proportional to one another.
where r i,j is the distance between molecules i and j, ǫ = 0.65 kJ/mol and σ = 0.3166 nm.
We perform MD simulations for a system of N = 1728 molecules at density ρ = 1.0 g/cm 3 , 210 K ≤ T ≤ 300 K, with periodic boundary conditions and a simulation time step of 1 fs.
The temperature is controlled by the Berendsen method of rescaling the velocities [40]. The long-range Coulombic interactions [41] are treated with the reaction field technique with a cutoff of 0.79 nm. For each state point, we run two independent simulations to improve statistics.
To estimate the orientational autocorrelation time of water molecules in the supercooled regime, we average the scalar product of the normalized dipole vectors µ i of each water molecule i in the system, where θ i (t) is the angle between µ i (t) and µ i (0). This function corresponds to the average of the Legendre polynomial P 1 (cos θ i (t)) evaluated for each molecule and can be directly measured by dielectric experiments.  Table I. Both parameters in Eq. (1), A and β, show weak dependences on T . The resulting values of these parameters are consistent with previous simulations of a smaller system of SPC/E molecules [23].
The estimated autocorrelation times τ a agree ( Fig. 2) with the power law behavior predicted by the MCT, We estimate T MCT = (194 ± 4) K and γ a = 2.9 ± 0.3, in agreement with previous results for similar densities and temperatures [24].
The estimated values of τ a , verify well the the von Schweidler law (Appendix) and the time-temperature superposition principle predicted by MCT, i.e. that the autocorrelation functions in the α-relaxation regime at different temperatures follow the same master curve if the time is rescaled by the autocorrelation time ( Fig. 1b) [29].

A. Definition and Methods
Here we define the randomization time τ r , a new quantity that we propose to characterize the orientational autocorrelation time. We consider the normalized dipole µ i of molecule i over a time interval ∆t = N δt,μ whereμ i is a function of δt and ∆t, t k ≡ kδt, and δt is the time interval between two consecutive samples of µ i .
If δt is greater than the autocorrelation time of µ i , then two consecutive samples µ i (t) average over all the molecules N in the system. Hence because ( µ i ) 2 = 1 for any t k , and This is the result of a freely jointed chain of N bonds of the same length, for which the mean square end-to-end distance is N 2 μ 2 i = N [42]. Therefore, if δt is larger than the orientational autocorrelation time for µ i , the µ rms decreases as 1/ √ ∆t.
If, instead, δt is shorter than the orientational autocorrelation time, consecutive elements in the sum in Eq. (6) are correlated µ i (t) · µ i (t + δt) = z, resulting in a smaller fluctuation.
This can be formally understood by considering the freely rotating chain model [42], where consecutive bonds in the chain are free to rotate, each around the axis of the previous bond, at an angle θ, such that cos(θ) = z. With this assumption, the resulting mean square end-to-end distance for n bonds of unit length is In the case of small θ, we have z = 1 − ǫ + O(ǫ 2 ), with ǫ = θ 2 /2 ≪ 1 and z n ≃ exp(−nǫ).
Then, from Eq. (8), we obtain In our problem, the bonds are dipole vectors sampled at time intervals δt, and n = ∆t/δt = N . Therefore Eq. (9) becomes The right-hand side of this equation behaves as 1/ √ ∆t for N → ∞, i.e., the random case behavior is recovered for large ∆t/δt. Therefore, if we define τ r as the time at which the correlation goes to zero as 1/ √ ∆t, it is possible to see that If we consider the fluctuation of any observable, the relation (11) defines the randomization time τ r for that observable [37] and τ r is equal to the smallest δt such that µ rms ∼ 1/ √ ∆t for any ∆t. The evaluation of τ r from a plot such as in Fig. 3 could be problematic, since it depends critically on the data errors. Therefore, to define in a clear way τ r , we fit the first eight points (∆t = δt, 2δt, ...., 8δt) using where λ = λ(δt). In this way we study how the deviation from the asymptotic regime decreases by increasing δt. We find that the exponent λ increases toward the asymptotic value 1/2 for increasing δt, and λ = 1/2 for any δt ≥ τ r (Fig. 4). We therefore define τ r as the extrapolated values of δt at which λ = 1/2. We find that λ approaches 1/2 as 1/δt, to the leading order, for low temperatures (Fig. 4).
The resulting values of τ r are presented in Fig.5a as functions of T − T MCT , showing that the power law behavior Eq. (4) is well satisfied by τ r . In this case our estimates are T MCT = (191.5 ± 2.5) K and γ r = 3.3 ± 0.2, both consistent within the errors with the estimates based on τ a (Fig.2). Therefore, the prediction (i) of MCT is verified.
By plotting τ r against τ a , we verify the MCT prediction (ii). We find (Fig. 6) that τ r and τ a are linearly related and that τ r is approximately five times larger than τ a .
The large value of τ r with respect to τ a is consistent with the fact that the latter measures the decay of the self-dipole correlation to a finite value, while the first measures the time needed for the self-dipole autocorrelation to decay to zero. This result is also reminiscent of the recent MD analysis in bulk water for the site-dipole field, a measure of the average orientation of the molecules passing through each spatial position, recently introduced in Ref. [35]. Higo et al. [35] find coherent patterns for the site-dipole field, at ambient pressure and T = 298 K and T = 300 K, that persist for more than 100 ps, a time much larger than the single molecule orientational relaxation time τ a of approximately 5 ps (Table I). It is, therefore, interesting to calculate the randomization time τ box and to find its relation with the autocorrelation time τ a for T → T MCT .

V. THE SITE-DIPOLE FIELD
To check if there are large correlated domains of dipoles in water which have large relaxation times compared to the individual dipole correlation time, we next study site-dipole field introduced by Higo et. al. [35]. We define the instantaneous coarse-grained site-dipole as the average of dipoles µ i of all the molecules n i (t) at time t belonging to box i of edge L box , volume v = L 3 box and centered at r i . If n i (t) = 0, then d v i = 0 by definition [43,44]. We chose vectors r i in such a way that the corresponding boxes do not overlap [45]. The time averaged v i over an interval ∆t, is defined analogously to Eq. (5). The rms average d v rms , is defined in analogy to Eq. (6) and (7), but instead of summation over all molecules we perform a summation over all boxes.
Since the argument presented for µ rms is also valid for d v rms , the relation (11) also holds for d v rms and allows us to estimate the randomization time τ box for d v rms . We find that τ box , calculated for L box = 3.33Å , diverges at T box = (194 ± 2) K with a power law with exponent γ box = 3.2 ± 0.2, consistent with our estimates of γ a and T MCT , respectively (Fig. 7).
If we compare τ box with τ r (Fig. 8), we again find a linear relation, as in Fig. 6 for τ a , consistent with the MCT statement (ii). The proportionality factor is approximately 2.5 [46], smaller than the factor ≈ 5 found for τ r in Fig. 6. Therefore, we conclude that in bulk water the coarse-grained site-dipole randomization time τ box is larger than the self-dipole autocorrelation time τ a , but smaller than τ r . Thus we do not find a significant increase in the box dipole autocorrelation time compared to the autocorrelation time τ a . These results do not support the results of Refs. [21,35].
To test the existence of cooperative domains in the SPC/E model, we perform coarse-

VI. DISCUSSION
Considerable numerical evidence shows that MCT predictions apply to orientational dynamics of water, despite the fact that MCT has been developed for particles interacting through spherically symmetric potentials [47]. However, recent extensions of MCT to liquids of linear molecules [48,49], and single solute molecules in a simple solvent liquid [50], confirm the main MCT predictions about the orientational autocorrelation functions [33].
Our study of supercooled water confirms the validity of MCT predictions for the orientational autocorrelation time τ a , estimated through a stretched exponential of the dipole autocorrelation function, for the temperature range 210 K≤ T ≤ 300 K at density ρ = 1 g/cm 3 .
Our results agree with the time-temperature superposition principle and the power law Eq. (4), with T MCT = (194 ± 4) K and γ a = 2.9 ± 0.3.
By evaluating the randomization time τ r , defined as the time needed to randomize the molecular dipoles, we verify the MCT prediction that all the characteristic times of quantities coupled to density fluctuations of a supercooled liquid are proportional to each other and follow the same power law Eq. (4). We find τ r ≈ 5 τ a , with T MCT = (191.5 ± 2.5) K and γ r = 3.3 ± 0.2, consistent with the estimates based on the calculation of τ a .
We also calculate the randomization time τ box for the box dipole field, a quantity introduced in Ref. [35] to measure the local orientational memory of molecules passing through a given spatial position. Our results for L box = 3.33Åshow that τ box diverges at T box = T MCT , following a power law with exponent γ box = γ a , and that τ box ≈ τ r /2. As a consequence, the local memory is lost faster than the self-dipole orientational memory.
Our results also show the existence of domains of correlated dipoles of short spatial range, with a correlation length comparable to the dipole-dipole correlation length at ambient T [44]. Whether this conclusion is specific to the SPC/E model with reaction field is an open question, and requires further investigation using other models of water, e.g. polarizable models.

Acknowledgments
We thank N. Giovambattista for his helpful collaboration on initial phase of this work.
We thank the NSF Chemistry Program for support. GF thanks the Spanish Ministerio de Educación y Ciencia (Programa Ramón y Cajal and Grant No. FIS2004-03454), and M.
Sasai for his hospitality during a visit to Nagoya University.

APPENDIX A: THE VON SCHWEIDLER LAW
The MCT predicts that the autocorrelation function departs from the plateau A as a power law with exponent b, known as the von Schweidler law, where the von Schweidler exponent b does not depend on T . We verify that at lower temperatures Eq. (A1) holds for roughly two decades in time (Fig. 10) and we find a clear deviation only for T ≥ 260 K at short times, possibly due to the fact that for T ≥ 260 K it is more difficult to estimate the plateau A. The estimated value of b is 0.6 ± 0.1, consistent with previous results [29] and with the MCT prediction that γ a , a, and b are related by the equation Here a is the exponent of the power law that describes the short-time approach to the plateau , and a is related to b by the transcendental equation