Predictions of Dynamic Behavior Under Pressure for Two Scenarios to Explain Water Anomalies

Using Monte Carlo simulations and mean field calculations for a cell model of water we find a dynamic crossover in the orientational correlation time $\tau$ from non-Arrhenius behavior at high temperatures to Arrhenius behavior at low temperatures. This dynamic crossover is independent of whether water at very low temperature is charaterized by a ``liquid-liquid critical point'' or by the ``singularity free'' scenario. We relate $\tau$ to fluctuations of hydrogen bond network and show that the crossover found for $\tau$ for both scenarios is a consequence of the sharp change in the average number of hydrogen bonds at the temperature of the specific heat maximum. We find that the effect of pressure on the dynamics is strikingly different in the two scenarios, offering a means to distinguish between them.

• The singularity-free (SF) scenario hypothesizes the presence of a line of temperatures of maximum density T MD (P ) with negative slope in the (T, P ) plane. As a consequence, K T and |α P | have maxima that increase upon increasing P , as shown using a cell model of water. The maxima in C P do not increase with P , suggesting that there is no singularity [4] [ Fig. 2(b)].
Above the homogeneous nucleation line T H (P ) where data are available, the two scenarios predict roughly the same equilibrium phase diagram. Here we show that dynamic measurements should reveal a striking difference between the two scenarios. Specifically, the low-T dynamics depends on local structural changes, quantified by the variation of the number of hydrogen bonds, that are affected by pressure differently for each scenario. We find this result by studying-using Monte Carlo (MC) simulations and mean field calculations-a cell model which has the property that by tuning a parameter its predictions conform to those of either the LLCP or the SF scenario. This cell model is based on the experimental observations that on decreasing P at constant T , or on decreasing T at constant P , (i) water displays an increasing local tetrahedrality [5], (ii) the volume per molecule increases at sufficiently low P or T , and (iii) the O-O-O angular correlation increases [6].
The entire system is divided into cells i ∈ [1, . . . , N], each containing a molecule with volume v ≡ V /N, where V ≥ Nv hc is the total volume of the system, and v hc is the hardcore volume of one molecule. The cell volume v is a continuous variable that gives, in d dimensions, the mean distance r ≡ v 1/d between molecules. The van der Waals interaction is represented by a potential with attractive energy ǫ > 0 between nearest-neighbor (n.n.) molecules and a hard-core repulsion at For a regular square lattice, each molecule i has four bond indices σ ij ∈ [1, . . . , q], corresponding to the four n.n. cells j, giving rise to q 4 different molecular orientations. Bonding and intramolecular (IM) interactions are accounted for by the two Hamiltonian terms where the sum is over n.n. cells, 0 < J < ǫ is the bond energy, δ a,b = 1 if a = b and δ a,b = 0 otherwise, and where (k,ℓ) i denotes the sum over the IM bond indices (k, l) of the molecule i and J σ > 0 is the IM interaction energy with J σ < J, which models the angular correlation between the bonds on the same molecule. The total energy of the system is the sum of the van der Waals interaction and Eqs. (1) and (2).
At constant P , the density of water decreases for T < T MD (P ) which the model takes into account by increasing the total volume by an amount v B > 0 for each bond formed.
Hence the total molar volume v of the system is where v free is a variable for the molar volume without taking into account the bonds, p B = N B /(2N) is the fraction of bonds formed and N B is the number of bonds [4,7].
We perform simulations in the NP T ensemble [7] for q = 6, v B /v hc = 0.5, J/ǫ = 0.5, and for two different values of J σ /ǫ: (i) J σ /ǫ = 0.05, which gives rise to a phase diagram with a LLCP [ Fig. 1(a)], and (ii) J σ = 0, which leads to the SF scenario [4]. We study two square lattices with 900 and 3600 cells, and find no appreciable size effects. We collect statistics over 10 6 MC steps after equilibrating the system for all P and T .
For J σ /ǫ = 0.05, |α P | for P < P C ′ displays a maximum, α max P [ Fig. 1(b)]. As P increases, α max P increases and shifts to lower T , converging toward T W (P ) [ Fig. 1(a)]. We find that the number of bonds, N B , increases on decreasing T , and at constant T decreases for increasing P , and is almost constant at T W (P ) [8]. This is consistent with trends seen both in experiments [5] and in simulations [9], suggesting that for T > T W (P ) the liquid is less structured and more HDL-like, while for T < T W (P ) it is more structured and more LDL-like.
We find that |dp B /dT | shows a clear maximum for all P < P C ′ which shifts to lower T upon increasing P [ Fig. 1(c)]. Remarkably, we also find that the locus of |dp B /dT | max coincides with the Widom line T W (P ) [ Fig. 1(a)] and that the value of |dp B /dT | max increases on approaching P C ′ . This is the same qualitative behavior as |α P (T )| max and C P (T ) max , which are used to locate T W (P ) [Figs. 1(b) and 2(a)]. The relation of |dp B /dT | with the fluctuations is revealed by its proportionality to |α P (T )| and to the fluctuation of the number of bonds where k B is the Boltzmann constant.
For J σ = 0 (SF scenario) we observe no difference for the behavior of N B and |dp B /dT |.
We further verify the prediction of the SF scenario [4] that C max P remains a constant upon Next, we study how this different behavior affects the dynamics. Previous simulations [10] found a crossover from non-Arrhenius to Arrhenius dynamics for the diffusion constant of models that display a LLCP, and showed the temperature of this crossover coincided with T W (P ). We calculate, for both scenarios, the relaxation time τ of S i ≡ j σ ij /4, which quantifies the degree of total bond ordering for site i. Specifically, we identify τ as the time For both scenarios we find a dynamic crossover (Fig. 3). At high T , we fit τ with the Vogel-Fulcher-Tamman (VFT) function where τ VFT 0 , T 1 , and T 0 are three fitting parameters. We find that τ has an Arrhenius T where τ 0 is the relaxation time in the high-T limit, and E A is a T -independent activation energy. We find that for J σ /ǫ = 0.05 the crossover occurs at T W (P ) for P < P C ′ [ Fig. 3(a)], and that for J σ = 0 the crossover is at Fig. 3(b)]. We note that for both scenarios the crossover is isochronic, i.e. the value of the crossover time τ C is approximately independent of pressure.
We next calculate the Arrhenius activation energy E A (P ) from the low-T slope of log τ vs. 1/T [ Fig. 4(a)]. We extrapolate the temperature T A (P ) at which τ reaches a fixed macroscopic time τ A ≥ τ C . We choose τ A = 10 14 MC steps > 100 sec [11] [ Fig. 4(b)]. We find that E A (P ) and T A (P ) decrease upon increasing P in both scenarios, providing no distinction between the two interpretations. Instead, we find a dramatic difference in the P dependence of the quantity E A /(k B T A ) in the two scenarios, increasing for the LLCP scenario and approximately constant for the SF scenario [ Fig. 4(c)].
We can better understand our findings by developing an expression for τ in terms of thermodynamic quantities, which will then allow us to explicitly calculate E A /(k B T A ) for both scenarios. For any activated process, in which the relaxation from an initial state to a final state passes through an excited transition state, ln(τ where ∆(U + P V − T S) is the difference in free energy between the transition state and the initial state. Consistent with results from simulations and experiments [12,13], we propose that at low T the mechanism to relax from a less structured state (lower tetrahedral order) to a more structured state (higher tetrahedral order) corresponds to the breaking of a bond and the simultaneous molecular reorientation for the formation of a new bond. The transition state is represented by the molecule with a broken bond and more tetrahedral IM order. Hence, where p B and p IM , the probability of a satisfied IM interaction, can be directly calculated.
To estimate ∆S, the increase of entropy due to the breaking of a bond, we use the mean We next test that the expression of ln(τ /τ 0 ), in terms of ∆S and Eq. (6), describes the simulations well, with minor corrections at high T . Here τ 0 ≡ τ 0 (P ) is a free fitting parameter equal to the relaxation time for T → ∞. From Eq. (7) we find that the ratio E A /(k B T A ) calculated at low T increases with P for J σ /ǫ = 0.05, while it is constant for J σ = 0, as from our simulations [ Fig. 4(d)].
In summary, we have seen that both the LLCP and SF scenarios exhibit a dynamic crossover at a temperature close to T (C max P ), which decreases for increasing P . We interpret the dynamic crossover as a consequence of a local breaking and reorientation of the bonds for the formation of new and more tetrahedrally oriented bonds. Above T (C max P ), when T decreases, the number of hydrogen bonds increases, giving rise to an increasing activation energy E A and to a non-Arrhenius dynamics. As T decreases, entropy must decrease. A major contributor to entropy is the orientational disorder, that is a function of p B , as described by the mean field expression for ∆S. We find that, as T decreases, p B -hence the orientational order -increases. We find that the rate of increase has a maximum at T (C max P ), and as T continues to decrease this rate drops rapidly to zero -meaning that for T < T (C max P ), the local orientational order rapidly becomes temperature-independent and the activation energy E A also becomes approximately temperature-independent, for the Eq.(6). Corresponding to this fact the dynamics becomes approximately Arrhenius.
We find that the crossover is approximately isochronic (independent of the pressure) consistent with our calculations of an almost constant number of bonds at T (C max P ). In both scenarios, E A and T A decrease upon increasing P , but the P dependence of the quantity   that the same behavior is found using the mean field approximation. In all the panels, where not shown, the error bars are smaller than the symbol sizes.