We use molecular dynamics simulations to study relations between thermodymamic, structural, and dynamical properties of TIP4P/2005 water models with systematically reduced partial charges and, thus, weaker hydrogen bonds. Observing a crossing of isochores in the P–T diagram, we show that these water-like models have a readily accessible liquid–liquid critical point (LLCP) associated with a transition between high-density liquid (HDL) and low-density liquid (LDL) forms and determine the dependence of the critical temperature *T*_{c}, pressure *P*_{c}, and density *ρ*_{c} on the charge-scaling factor from fits to a two-structure equation of states. The results indicate that the water-like models exhibit liquid polyamorphism in a wide range of interaction parameters. Considering elongated systems, we observe a decomposition into extended and stable HDL-like and LDL-like regions at appropriate pressures and low temperatures and analyze the respective structural and dynamical properties. We show that the diverse local order results in very different correlation times of local dynamics, while the fragility is hardly changed. The results yield insights into the origin of a dynamical crossover, which is observed when lowering the temperature along isobars and was previously interpreted in terms of a fragile-to-strong transition. Our findings imply that the effect does not involve two liquid phases with an exceptionally large difference of the fragility but rather a high temperature dependence near the LLCP results from a rapid conversion from HDL-like environments with faster dynamics to LDL-like ones with slower dynamics.

## I. INTRODUCTION

The anomalies of water are remarkable examples of competing many-particle effects. Most researchers agree that they result from the competition between hydrogen bonding and other intermolecular interactions and originate in the supercooled regime of the liquid,^{1} where thermodynamic response functions appear to diverge upon cooling below the homogeneous nucleation temperature,^{2} which was recently specified as *T*_{H} ≈ 231 K in studies on micrometer-sized water droplets.^{3} In particular, it was conjectured that the anomalies of water are caused by a second critical point,^{4} which terminates a liquid–liquid phase transition (LLPT) line at reduced temperatures *T* and elevated pressures *P*. The proposed liquid phases of water differ in density motivating the denotations as high-density liquid (HDL) and low-density liquid (LDL).^{5–9} In addition, it was argued that they show diverse structures, e.g., the tetrahedral order is lower in the HDL than that in the LDL phase, and distinguishable dynamics.^{10}

The rapid crystallization of water, however, interferes with a straightforward characterization of the structural and dynamical properties of the proposed HDL and LDL phases. The term no man’s land was coined to denote these difficulties pertaining to the studies on supercooled water in the temperature range of ∼150 K–231 K. Therefore, despite recent progress in the reduction of the no man’s land,^{11–14} the existence of the HDL and LDL states and, thus, of the LLPT line is still a subject of controversial discussions.^{15,16} Notwithstanding, experimental evidence for the liquid polyamorphism of water is provided, e.g., by the existence of amorphous ice phases, which are separated by a first-order-like transition and show glass transitions at different temperatures *T*_{g}.^{5–9,17}

A liquid polyamorphism is not unique to water. A recent review article provides a comprehensive overview.^{18} Experimental evidence for this phenomenon was also reported for, e.g., carbon,^{19} phosphorous,^{20,21} and sulfur.^{22} In simulation work, liquid polyamorphism was observed for various systems ranging from patchy particles with tetrahedral functionality^{23} to the spherically symmetric two length scale Jagla model.^{24} Moreover, even simple theoretical approaches can produce this feature.^{25,26} These studies imply that a LLPT line and the associated liquid–liquid critical point (LLCP) can occur when spacious low-energy and low-entropy local structures such as tetrahedral order exist. In such a situation, cooling promotes the LDL phase consisting of these structural motifs, while pressurizing favors the HDL phase with more compact particle arrangements.

Molecular dynamics (MD) simulations are well suited for investigations of liquid polyamorphism.^{27–31} In studies on water, it is advantageous that high cooling rates allow one to avoid crystallization. A possible drawback is that the LLCP, if existent, often occurs on the edge of the temperature–pressure region accessible to current computational approaches. Nevertheless, MD simulations were successfully used to scrutinize the existence of a LLCP for various water models.^{32} While evidence for a LLPT was reported for the TIP4P/2005,^{33–37} TIP4P/Ice,^{37} TIP5P,^{36} WAIL,^{38} E3B3,^{39} and ST2^{4,40,41} models, this phenomenon does not occur for the mW model, where the entropy of mixing is too high for a phase separation.^{42} For the SPC/E model, a LLCP was proposed to exist in a parameter range inaccessible to current MD simulations.^{43} Even if observed, appropriate equilibration of the systems may be an issue owing to the slow structural relaxation near the LLCP, and opposing views can exist relating to the location of the LLCP^{33–36} and the effects of ice-like clusters,^{44,45} as was reported for TIP4P/2005 water.

For the theoretical description of polyamorphism, the two-structure equation of states (TSEOS) formalism was successfully applied.^{36,41,42,46–52} The TSEOS approach assumes that the Gibbs energy of HDL–LDL mixtures is given by the contributions of the individual phases supplemented by a nonideal mixing term, which can be enthalpic^{35,41,53,54} or entropic^{42,48,55} in nature. Although it is based on a mean-field approximation, it often enables a good description of the data and a reliable determination of the location of the LLCP,^{52} i.e., of the critical temperature *T*_{c}, pressure *P*_{c}, and density *ρ*_{c}, including cases where a modest extrapolation to inaccessible state points is required.

This literature survey suggests that liquid polyamorphism is more abundant than anticipated. However, it is still elusive for which types of particle interactions the phenomenon occurs and at which state points the LLPT is located. A major advantage of MD simulations is that the interactions can be varied continuously and that temperature, pressure, and density changes are straightforward. These capabilities were exploited to ascertain the LLPT of various liquids. For the Stillinger–Weber potential of silicon,^{56} it was found that the LLPT is erased when the three-body repulsion parameter is reduced.^{57,58} For the ST2 model of water, it was observed that the LLPT becomes more prominent when the bond angle is varied.^{59} Such variations are of particular relevance for water, for which a unique model does not yet exist so that the dependence of the properties on the parameters is an important issue.

Another efficient way to alter the interactions in chemically realistic models is to vary the partial charges of the atoms. While the model of Woodcock, Angell, and Cheeseman (WAC) of silica did not show a LLCP in the accessible region of the phase diagram,^{28} a clear crossing of isochores in the P–T phase diagram and extensive free energy calculations indicated a LLCP when the partial charges of this model are decreased in magnitude.^{29,60} Moreover, charge scaling proved to be very useful to investigate the role of the particle interactions in the glass transition of water,^{61,62} silica,^{63} and an ionic liquid.^{64} In particular, this approach allows one to tune the strength of the hydrogen bonds and the degree of the tetrahedral order in water models.^{61,62}

The conjectured liquid polyamorphism is assumed to result also in peculiar dynamical behaviors of water. For example, the dynamics of water speeds up when the pressure is increased.^{1} To rationalize this anomalous effect, one may argue that the HDL replaces the LDL upon pressurizing and that the former shows faster dynamics than the latter. Consistent with this conjecture, it was reported that the HDL has a lower glass transition temperature *T*_{g} than the LDL,^{17} but the interpretation of the experimental observations is still controversial.^{65} In an extensive simulation study on the TIP4P/2005 water model, different glass transition temperatures were clearly identified, whereby the pressure dependence *T*_{g}(*P*) was found to be normal for the HDL, but anomalous for the LDL.^{66} Furthermore, a LLPT may be accompanied by a fragile-to-strong (FS) dynamical transition.^{67} Explicitly, assuming that the HDL and LDL are fragile and strong liquids, i.e., deviate from and conform to Arrhenius behaviors, respectively, a FS transition is expected in temperature-dependent studies, when, depending on the pressure, the LLPT line or the Widom line are crossed, where the latter marks the position of the maxima of thermodynamic response functions in the one-phase region above the LLCP.^{68–71} However, alternative explanations for crossovers in the temperature-dependent dynamics of water were also put forward.^{72} For example, the effect was attributed to an onset of hopping motions or secondary relaxations.^{70,73} MD simulations proved to be very useful to investigate the dynamical behavior of the ST2 water model near the LLPT.^{74}

Here, we perform MD simulations to ascertain the thermodynamic, structural, and dynamical properties of TIP4P/2005 water models with reduced partial charges and, thus, weaker hydrogen bonds. Observing an isochore crossing in the P–T phase diagram, it is shown that the studied water-like models have LLCPs. Performing a TSEOS analysis, the dependence of the critical parameters on the magnitude of the partial charges is determined. The analysis reveals that charge scaling shifts the LLCP to computationally better accessible state points. This effect is a major advantage of our approach and allows us to follow the evolution of the structure and dynamics of the charge scaled models across the LLCP. To rationalize the findings, we consider systems, which decompose into extended HDL-like and LDL-like spatial regions and determine the structure and dynamics in these domains. The findings show that a FS-like dynamical crossover is mainly related to a rapid conversion from HDL-like environments with faster dynamics to LDL-like ones with slower dynamics.

## II. METHODS

### A. Simulation details

We performed MD simulations for several charge scaled TIP4P/2005 water models. Specifically, we started from the original TIP4P/2005 model^{75} and scaled the partial charges of all atoms by a factor *q* ranging between 0.86 and 0.91, while all other interaction parameters remained unchanged. Thus, the water-like molecules of the charge scaled TIP4P/2005 models are less polar and form weaker hydrogen bonds.

The simulations were conducted using the GROMACS simulation package.^{76,77} Periodic boundary conditions were employed and the leapfrog integrator with a time step of 2 fs was used. The particle mesh Ewald method^{78} was utilized to calculate the Coulomb and Lennard-Jones interactions where the cutoff was set to 1.0 nm and the Fourier spacing to 0.16 nm^{−1}. All systems were carefully equilibrated prior to data acquisition. Specifically, the equilibration runs and the production runs were at least several ten times longer than the rotational correlation time at the given temperature.

The majority of the studied systems were cubic and comprised *N* = 2000 molecules, leading to box lengths of ∼4 nm. Most simulations of these systems were performed in the isothermal–isochoric (NVT) ensemble, applying the Nosé–Hoover thermostat.^{79,80} For simulations in the isothermal–isobaric (NPT) ensemble, we additionally employed the Parrinello–Rahman barostat.^{81} Moreover, we utilized elongated rectangular simulation boxes, which facilitate a decomposition into two phases because they allow for smaller interface-to-volume ratios than their cubic counterparts.^{29,82} These coexistence simulations were performed in the NVT ensemble for *q* = 0.86 and *q* = 0.88. We chose a ratio of the box dimensions of *L*_{x} : *L*_{y} : *L*_{z} ≡ 1 : 1 : 3 for systems comprising *N* = 2000 molecules and carried out additional simulations for *N* = 12 000 molecules and a ratio of 1:1:6. With these choices, the smaller box dimensions were *L*_{x} = *L*_{y} ≥ 2.7 nm and *L*_{x} = *L*_{y} ≥ 4.0 nm for the smaller and larger elongated systems, respectively.

### B. Two-structure equation of states

To evaluate the isochores and determine the LLCP, we follow previous studies^{35,36,41,51} and perform a TSEOS analysis. This approach assumes that the liquid consists of an interconverting mixture of a high-density structure A with a fraction 1 − *x* and a low-density structure B with a fraction *x*. The molar Gibbs energy is written as^{35,36,41,51,52}

Here, *G*^{BA} = *G*^{B} − *G*^{A} is the difference of the Gibbs energies *G*^{A} and *G*^{B} of the pure states. The third summand accounts for the entropy of mixing the structures A and B, where the first two terms of this contribution describe ideal mixing and the last one nonideal contributions. Because the determination of the LLCP is of interest, it is convenient to express temperature and pressure as dimensionless variables, $\Delta T^=T\u2212TcTc$ and $\Delta P^=P\u2212Pc\rho cRTc$. *G*^{A} is modeled as a third-order expansion around the LLCP,

with $\u0109mn$ as adjustable coefficients. *G*^{BA} determines the *A ⇌ B* equilibrium; explicitly, the equilibrium constant amounts to *K*(*T*, *P*) = exp(−*G*^{BA}/*RT*). It is approximated as

The meaning of the parameters *λ*, *a*, and *b* was discussed in detail in previous studies.^{35,41,51} The LLPT, LLCP, and the Widom line, i.e., the continuation of the LLPT line beyond the LLCP, are determined by the condition ln *K*(*T*, *P*) = 0.

The nonideal contribution to the entropy of mixing, *ωx*(1 − *x*), is assumed to be symmetric in *x*. The parameter *ω* is a measure for the relevance of this effect. While a temperature-independent value of *ω* indicates an entropic origin, $\omega \u221d1T$ is found when the nonideality is enthalpic in nature. In previous works on the TIP4P/2005 water model,^{35,51} it was shown that the data are well described by the form

with $T^=T/Tc$ and *ω*_{0} as the adjustable parameter. Therefore, we also apply this form to all modified TIP4P/2005 models.

The equilibrium fraction *x*_{e} is defined by the condition

This equation has to be treated numerically. When describing isochoric data, the molar volume can be obtained from the partial derivative of the dimensionless Gibbs energy $\u011c=GRTc$ with respect to $P^=P\rho cRTc$. In a recent work,^{51} the TSEOS analysis was augmented with contributions from the liquid–vapor spinodal to improve the performance in a wide range of state points. We found that this extension was not necessary in our study, which was restricted to the vicinity of the LLCP, in particular, to *T* < 1.25*T*_{c}. More details about our TSEOS analysis can be found in the supplementary material.

### C. Observables

#### 1. Structure

We use two local structure identifiers to distinguish between HDL-typical and LDL-typical environments. In the following, we will refer to these particle arrangements as high-density (HD) and low-density (LD) local structures, respectively. For their calculation, we consider a molecule *i* and obtain the distances *d*_{i,j} to its neighboring molecules *j* from the positions of the respective oxygen atoms. Moreover, we determine the number *n*(*i*) of neighbors within the cutoff distance *d*_{i,j} < 0.37 nm. Ordering the values of *d*_{i,j} according to increasing distance, the local structure index (LSI) of molecule *i* is defined as^{83}

Here, Δ(*i*, *j*) = *d*_{i,j+1} − *d*_{i,j} and $\Delta \xaf(i)$ is the average of these differences over all *n*(*i*) neighbors. For tetrahedral LD local structures, the difference between the *d*_{i,j} values is small for the first four neighbors but large between the fourth and fifth neighbors, resulting in large LSI values. For HD local structures, by contrast, the distribution of the Δ(*i*, *j*) is narrower, and thus, the LSI values are smaller. Moreover, we delineate a simple structure identifier *N*_{4}(*i*) based on the number *n*(*i*) of neighboring molecules. In detail, we define *N*_{4}(*i*) = 0 if *n*(*i*) > 4 and *N*_{4}(*i*) = 1 if *n*(*i*) ≤ 4. Thus, the average value of *N*_{4} is small in HDL-like regions with more than four neighbors in the next neighbor shell, while it is large in tetrahedral LDL-like neighborhoods. Likewise, the mere number of neighbors was used to define a structure identifier in a previous computational approach to water properties near the LLCP.^{74} In the supplementary material, we also analyze the distance to the fifth neighbor, *d*_{5} ≡ *d*_{i,5}, and exploit that low and high values are indicative of HD and LD local structures, respectively.

#### 2. Dynamics

To characterize the translational motion of the water-like molecules, we analyze the time-dependent positions $r\u2192i(t)$ of their oxygen atoms. Specifically, we calculate the mean-square displacement (MSD)

where the pointed brackets denote the averages over all molecules *i* and various time origins *t*_{0}, and determine the self-diffusion coefficient *D* by fitting the MSD in the diffusive regime at long times to *r*^{2}(*t*) = 6*Dt*.

For the analysis of the reorientation dynamics, we consider the O–H bond vectors of the water-like molecules and compute the rotational correlation function (RCF),

Here, the unit vectors $e\u2192i(t)$ describe the time dependence of the O–H bond orientations. To quantify the time scale of the molecular reorientation, we fit the RCF *F*_{1}(*t*) beyond the vibrational regime to a scaled Kohlrausch–Williams–Watts (KWW) function,

Here, the factor *A* is used to consider that vibrational motion causes a minor loss of correlation prior to the onset of structural relaxation. Moreover, *τ*_{kww} and *β*_{kww} are the time constant and the stretching parameter of the KWW function, respectively. From the fit results, we calculate the mean correlation time $\tau =1A\u222b0\u221eF1(t)dt$ according to *τ* = (*τ*_{kww}/*β*_{kww})Γ(1/*β*_{kww}), where Γ(*x*) denotes the gamma function.

#### 3. Spatially resolved analyses

In the simulations for elongated boxes near the LLCP, we observe a time-dependent decomposition into HDL-like and LDL-like fractions along the *z* dimension. To characterize the space–time characteristics of the structural and dynamical properties in these regions, we follow the respective behaviors in layers, which are parallel to the *xy* plane and have a width of Δ*z* = 1 nm, along the trajectories. To analyze the density fluctuations, we define the layer-averaged density *ρ*_{L}(*z*, *t*) based on the number of molecules in the layer centered about the position *z* at time *t*. Likewise, we average over the molecules in the respective layers to obtain *z*- and *t*-dependent values of layer-averaged structural identifiers, *I*_{L} and *N*_{4,L}. For *N* = 2000, these averages include about 250 molecules. To analyze the overall temporal aspect of the fluctuations, we calculate correlation functions of the layer-averaged quantities,^{84}

where *ξ*_{L} denotes the respective random variable and *σ*_{ξ} is its standard deviation. Moreover, we perform KWW fits to obtain mean time constants *τ*_{ξ} of these fluctuations. In these analyses, vibrational motions and numerical uncertainties lead to rapid changes in the numbers of particles within the layers, which can, however, be removed by averaging the coordinates over the time scale of the vibrational motion. To investigate the local dynamics–structure relations, we calculate *F*_{1}(*t*) separately for molecules in structurally different environments. Specifically, we distinguish between molecules characterized by various layer-averaged densities or structures at the time origin *t*_{0}, where the latter are obtained from the layers centered about the instant molecular positions, as described above.

## III. RESULTS

### A. TSEOS analysis

Commonly, liquids show a monotonic decrease in the pressure *P* when the temperature *T* is reduced along an isochore. However, when two liquid phases with different densities compete, the pressure can increase upon cooling. A crossing of isochores in the P–T diagram indicates a LLPT.^{28,29,60,85} Therefore, we calculate the pressure along isochores of the charge scaled water models from the virial stress. The P–T diagrams for the lowest (*q* = 0.86) and highest (*q* = 0.91) studied partial charges are shown in Fig. 1. For both *q* values, we observe that the pressure decreases when the temperature is reduced along isochores for sufficiently high densities, while it increases along those with sufficiently low densities in the studied temperature range. More importantly, the isochores for appropriate intermediate densities cross, indicating the existence of a LLPT. At low temperatures, the isochores end where the equilibration times exceed 2 *µ*s. We see that the metastable region and the expected spread of the crossing isochores below the LLPT can be explored to some degree. The highest temperature, for which an isochore crossing is observed, and the corresponding pressure and density mark the LLCP. Qualitatively, the same behavior is observed for all studied scaling factors between *q* = 0.86 and *q* = 0.91; see the supplementary material.

Taking the isochore crossing as evidence for the existence of two competing liquid states, we use the TSEOS approach to further analyze the simulation data. The results of this analysis are included as solid lines in Fig. 1. We see that the TSEOS approach well describes the simulation results for *q* = 0.86 and *q* = 0.91 in the vicinity of the LLCP. Similar agreement is observed for the other studied scaling factors, in harmony with previous findings for various model liquids with the tetrahedral order.^{35,36,41,51} The parameters resulting from the TSEOS fits for all studied scaling factors *q* are compiled in the supplementary material.

In Fig. 2, we depict the location of the LLCP for the studied water-like models, as obtained from the TSEOS analysis. We see that the critical temperature *T*_{c}, pressure *P*_{c}, and density *ρ*_{c} continuously decrease when the scaling factor *q* is reduced, consistent with previous findings for the WAC model of silica.^{60} Qualitatively, our findings can be rationalized by the fact that the electrostatic interactions and, thus, the hydrogen bonds are weaker for small values of *q*. Because of this *q* dependence of the interactions, the relevant thermal energy and, thus, *T*_{c} decrease. Moreover, weaker hydrogen bonds lead to a weaker tetrahedral order so that less pressure is required to push a fifth molecule into the next neighbor shell. Specifically, the LLCP is located at elevated pressures when the magnitude of the partial charges is lowered by less than 10% but at negative pressures when the molecular polarity is further reduced. Finally, the electrostatic interactions pull the molecules closer together than expected based on the Lennard-Jones parameters so that the density is reduced at small scaling factors *q*. To validate the outcome of the TSEOS analysis, we follow the evolution of the density *ρ* upon isobaric cooling at atmospheric pressure in the supplementary material. As expected, we find that *ρ* decreases very strongly near *T*_{c} for the higher scaling factors, for which *P*_{c} ≈ 1 bar, whereas such drops are not observed for the lower *q* values.

To estimate the location of the LLCP for the regular TIP4P/2005 model, we extrapolate the basically linear dependence of the critical parameters on the scaling factor to *q* = 1. We obtain *T*_{c} = 169 K, *P*_{c} = 136 MPa, and *ρ*_{c} = 1017 kg/m^{3}. In Table I, we observe that these values of *T*_{c} and *P*_{c} are smaller but still close to those reported for the TIP4P/2005 model in the literature, while *ρ*_{c} is in the range of the reported critical densities.^{33–35,87} We expect that the former differences result because the *q* range of the extrapolation is nearly twice as large as that of the analysis and the *q* dependence is not necessarily linear in a broader parameter range.

Study . | T_{c} (K)
. | P_{c} (MPa)
. | ρ_{c} (kg/m^{3})
. |
---|---|---|---|

Abascal and Vega^{33} | 193 | 135 | 1012 |

Sumi and Sekino^{34} | 182 | 158–162 | 1020 |

Yagasaki et al.^{86} | 185 | … | 1020 |

Singh et al.^{35} | 182 | 170 | 1017 |

Handle and Sciortino^{87} | 175 | 175 | 997 |

Debenedetti et al.^{37} | 177 | 175 | … |

Present extrapolation | 169 | 136 | 1017 |

### B. Structure and dynamics

Having determined the LLCP of the charge scaled TIP4P/2005 water models, we move onto their structure and dynamics. Figures 3 and 4 show results for *q* = 0.86 along the isochore with a density close to the critical one, *ρ* = 925 kg/m^{3} ≈ *ρ*_{c}.

In Fig. 3(a), we observe a weakly bimodal LSI distribution, in particular, at lower temperatures, where a peak at larger values *I* ≈ 0.0011 nm^{2} emerges and grows with respect to that at smaller values *I* ≈ 0.0002 nm^{2} upon cooling. Although the effect is less prominent than in an inherent structure analysis for the original TIP4P/2005 model,^{88} these findings indicate that two local structural motifs become preferred in the vicinity of the LLCP and that the HDL-like environments are transformed into LDL-like arrangements upon cooling. However, at none of the studied temperatures, we observe a decomposition into extended HDL and LDL phases for the studied cubic system.

To further analyze the interconversion between the HD and LD local structures, we determine their temperature-dependent fractions. For this purpose, we assume that a water-like molecule is in a LD local structure if *N*_{4} = 1 or *I* ≥ 0.0006 nm^{2}, where the latter value is roughly midway between the *p*(*I*) peaks. In Fig. 3(b), we see that the LD fractions *x*_{LD} obtained from the different structural identifiers and from the TSEOS analysis, respectively, gradually grow when the temperature is decreased along the studied isochore with *ρ* ≈ *ρ*_{c}. By contrast, a sharp rise in *x*_{LD} occurs near *T*_{c} when cooling under isobaric conditions with *P* ≈ *P*_{c}. In the supplementary material, we present the pressure dependence of *x*_{LD} at characteristic temperatures. We find that the LD fraction gradually decreases when the pressure is increased through *P*_{c} well above *T*_{c}, while a sharp drop occurs near the LLCP. Remarkably, well below *T*_{c}, the coexisting HD:LD compositions result in re-entrant *x*_{LD}(*P*) behaviors for all studied structural identifiers, also including *d*_{5}. This effect was not reported in an analogous study of the original TIP4P/2005 water model,^{35} where it is computationally more expensive to explore state points below the LLCP, once more demonstrating the added value of our charge-scaling approach.

In Fig. 4, the MSD and RCF indicate a prominent slowdown of the translational and rotational dynamics when the temperature is decreased. Typical of viscous liquids, the temporary caging of the molecules by their neighbors leads to a plateau of both quantities at times intermediate between the regimes of vibrational motion and viscous flow. To characterize the translational and rotational dynamics, we determine the self-diffusion coefficients *D* and the mean correlation times *τ* from *r*^{2}(*t*) and *F*_{1}(*t*), respectively. Analogous analyses are performed for various isochores and isobars. For none of the studied cubic systems, we find evidence for a clear bimodality of water dynamics due to coexisting HDL-like and LDL-like regions with different molecular mobilities.

Figures 5 and 6 present the diffusion coefficients and correlation times of the charge scaled TIP4P/2005 water model with *q* = 0.86 obtained along various isochores and isobars. We see that the translational and rotational dynamics strongly slow down when the density is decreased in the vicinity of *ρ*_{c}, in harmony with the dynamical anomalies reported for water and other tetrahedrally ordered liquids such as silica.^{1,60} When cooling along isochores, we do not find evidence for changes in the temperature dependence in the vicinity of the LLCP. In particular, there is no FS crossover for the isochore with the density *ρ* = 925 kg/m^{3} ≈ *ρ*_{c} in the studied temperature range. Rather, the diffusion coefficients and correlation times can be described by the Vogel–Fulcher–Tammann (VFT) equation, e.g.,

where the VFT parameters *τ*_{0}, *B*, and *T*_{0} can be described as second order polynomial functions of the density. On the other hand, we observe an intriguing temperature dependence, when we cool under isobaric conditions at pressures in the vicinity of *P*_{c} = 55 MPa. Here, isobaric conditions mean that we perform NPT simulations and apply a barostat to set a desired pressure *P* and to adjust the box volume *V* at each studied temperature *T* along the cooling path. Specifically, the temperature dependence of the diffusion coefficients *D*(*T*) at *P* = const. resembles that at *V* = const. at sufficiently high or low temperatures, while it is much stronger under isobaric conditions in the vicinity of the LLCP, resulting in a sigmoidal shape; see Fig. 5. Similar effects are found for the correlation times *τ*(*T*), as is demonstrated in the supplementary material. When comparing the results for various *P* values, we observe that the onset of the enhanced temperature dependence is consistent with the location of the Widom line. In harmony with our findings, the temperature dependence of the diffusion coefficients was found to be higher along isobars than along isochores near the LLCP,^{74} and the enhanced temperature dependence was reported to occur upon crossing the Widom line for the ST2 and TIP4P/2005 water models.^{4,89} Thus, the gradual and rapid interconversions between the HD and LD local structures observed above upon isochoric and isobaric cooling across *T*_{c}, respectively, are accompanied by corresponding changes in the translational and rotational dynamics. This is a first hint at a structural origin of the FS dynamical crossover for the studied charge scaled water models.

Next, we discuss the shape of the *F*_{1}(*t*) decays for the *q* = 0.86 system. For this purpose, we depict the KWW stretching parameters *β*_{kww} obtained along various isochores in Fig. 6. To compare the situation at isokinetic states, these data are not plotted as a function of the inverse temperature but of the correlation time. Along all isochores, the stretching parameter continuously decreases, and hence, the nonexponentiality of the decays increases when the motion slows down. In particular, we do not observe a minimum of *β*_{kww}(*τ*) near the LLCP, which would result if prominent structural inhomogeneity owing to coexisting HD and LD local structures at these state points caused enhanced dynamical heterogeneity, as was reported in recent work.^{90} At a given correlation time *τ*, *β*_{kww} is smaller for lower values of *ρ* than for higher ones. This observation may be unexpected because it implies that the dynamics is more heterogeneous when the tetrahedral order is, on average, higher at lower densities. However, the *β*_{kww}(*ρ*) dependence may be non-monotonous and show a weak minimum at *ρ* ≈ *ρ*_{c} owing to a moderate dynamical inhomogeneity associated with the LLCP.

The knowledge about the density dependence of the VFT parameters *τ*_{0}(*ρ*), *B*(*ρ*), and *T*_{0}(*ρ*) from the above analysis can be exploited to obtain the correlation time *τ*_{c} and, analogously, the diffusion coefficients *D*_{c} at the LLCP of the charge scaled water models from interpolations. In Fig. 2, we see that the translational and rotational dynamics at the critical point are significantly faster for smaller scaling factors *q*. This means that the LLCP is accessible more easily to simulations when the water-like molecules are less polar, which is an important advantage of our charge-scaling approach and will also be exploited in the following spatially resolved analysis. In addition, we find that the *q* dependence of the diffusion coefficients at the LLCP can be mapped onto that of the correlation times using the relation *a*^{2} = 6*D*_{c}*τ*_{c} and setting the jump length *a* to 0.325 nm, slightly larger than the next neighbor distance. Thus, short-range and long-range dynamics are related in the same way for all studied water-like models at the LLCP, while this is not necessarily true well inside the metastable region.

The findings from the isochoric and isobaric cooling schemes can be compared in more detail using the information from the above TSEOS and VFT analyses. Specifically, the TSEOS results yield the temperature-dependent densities *ρ* required to keep the pressure constant at a desired value upon cooling. Then, the knowledge about the density dependence of the VFT parameters can be employed to calculate the dynamics along an isobar from that along the isochores. In Fig. 5, we see that the, thus, calculated diffusion coefficients well agree with the above discussed NPT simulation results. Altogether, the results imply that the sigmoidal temperature dependence along isobars, which is particularly apparent when extrapolating the TSEOS–VFT based predictions to lower temperatures, cannot be interpreted in terms of a FS transition in the strict sense, i.e., as a result of a transition from an exceptionally fragile to a truly strong liquid. Rather, the diffusion coefficients and correlation times show a high temperature dependence in the vicinity of the LLCP because there is a rapid transition from a HDL-like liquid with fast dynamics to a LDL-like liquid with slow dynamics. Accordingly, when the pressure *P* of the isobar is reduced and the distance to the LLCP increases, the inflection points of 1/*D*(1/*T*) and *τ*(1/*T*) shift to higher temperatures and the sigmoidal shape becomes less prominent; see Fig. 5.

Similarly, previous studies used two-state models of water dynamics and argued that the FS transition results from a rapid change between HDL and LDL phases with different water mobilities. Specifically, Shi *et al.*^{91} proposed that there is no true FS transition between liquids with largely different fragilities but rather HDL and LDL waters obey Arrhenius laws albeit with different activation energies and a rapid crossover from one Arrhenius law to the other occurs. Likewise, Caupin and co-workers^{92,93} described experimental and computational results for the temperature- and pressure-dependent viscosity of supercooled water by a weighted superposition of two liquid behaviors and conjectured that the very high temperature dependence near the Widom line cannot be ascribed to either of the liquid phases but to their rapid interconversion. However, unlike Tanaka and co-workers, they reported that the high-temperature water phase is a fragile liquid and, hence, does not show Arrhenius behavior.

### C. Spatially resolved analyses

To scrutinize the liquid–liquid coexistence expected in the metastable region of the phase diagram and investigate the properties of both phases, we move from cubic to elongated systems and perform simulations at and below the LLCP. In doing so, we exploit that, for a given system size, the energetically unfavorable interface between the phases is minimized in elongated simulation boxes when it aligns parallel to the smallest box surface, in our case, the *xy* plane, facilitating phase separation. These simulations are restricted to the scaling factors *q* = 0.86 and, in particular, *q* = 0.88, where the LLCP is clearly detectable and the dynamics in the vicinity of this state point is relatively fast; see Fig. 2. It is ensured in the supplementary material that the pressure *P* at a given state point is the same within the error margins in the simulations for the cubic and elongated boxes.

In Fig. 7, we show the spatial variation of the density and structure for the charge-scaled model with *q* = 0.88 and *ρ* = 940 kg/m^{3} at three exemplary times during a trajectory at *T* = 125 K < *T*_{c}. Specifically, we depict the normalized density *ρ*_{L}/*ρ* and the structural identifiers *I*_{L} and *N*_{4,L} along the elongated *z* direction, as obtained from averaging over layers, which are Δ*z* = 1 nm thick and centered about *z*. Evidently, the system separates into two phases with significantly different properties. In particular, the layer-averaged density *ρ*_{L} is ∼10% lower in the center than at the edges of the simulation box, while the layer-averaged *I*_{L} and *N*_{4,L} parameters are higher in the former than the latter regions. The interface between both regions weakly changes during the course of the simulation. These findings indicate that the system decomposes into extended and stable LDL-like and HDL-like regions in the central and outer box regions. In the supplementary material, we present the results from analogous analyses at *T* = 132 K > *T*_{c}, which reveal that corresponding fluctuations of the structural properties also exist somewhat above the LLCP but their length and time scales are smaller so that a stable decomposition into extended LDL-like and HDL-like regions does not occur. We note in passing that, because of the periodic boundary conditions, the system is comprised of only one HDL-like region and one LDL-like region at 125 K so that the extension of these regions is artificially determined by the finite size of the simulation box.

Figure 8 compares the probability distribution of the layer-averaged local density *ρ*_{L} and structural identifiers *I*_{L} and *N*_{4,L} for the charge-scaled model with *q* = 0.88 at various temperatures. We see that all distributions have a Gaussian shape at *T* > *T*_{c} ≈ 131 K, while they are bimodal well below the LLCP. We emphasize that the considered distributions do not characterize the structural diversity of the immediate local environments, but rather the variation of the average properties in Δ*z* = 1 nm thick layers. Thus, for the observation of bimodal behaviors, it is not sufficient that HD and LD local structures exist but it is additionally required that the respective structural identity is largely preserved throughout the respective layer. Therefore, the observed crossovers from unimodal to bimodal behaviors indicate in our case that the length scale of structural fluctuations becomes larger than the thickness of the layers, and hence, the extension of HDL-like and LDL-like regions reaches the nanometer regime, in harmony with the above findings from the visual inspection of Fig. 7.

Next, we ascertain the relation between structure and dynamics for the charge scaled water models in more detail. To achieve this objective, one may tag the molecules with a structural identifier, e.g., the LSI, and characterize the dynamical behaviors of the, thus, defined subsets separately. If the used structural identifier changes rapidly in space and time, this approach will, however, not yield the true dynamics in HDL-like and LDL-like environments. Specifically, the occurring structural fluctuations affect the observed dynamical behaviors of the subsets, e.g., their correlation functions, to a degree, which depends on the relative length and time scales on which the structural identifier and the molecular dynamics are probed. Therefore, we exploit that extended and stable HDL-like and LDL-like regions exist in the elongated systems at lower temperatures. Moreover, we use layer-averaged quantities for tagging the molecules to ensure that the, thus, specified structural environments differ over sufficiently large regions, and we employ a local probe of the molecular dynamics, explicitly, the RCF *F*_{1}(*t*), to minimize the changes in the structural environments while assessing the respective dynamics. Therefore, we expect that our findings are characteristic for the rotational motions in regions with clearly distinguishable structures.

Figure 9 presents the results for the elongated *q* = 0.88 system at *T* = 125 K < *T*_{c}. Comparing *F*_{1}(*t*) for subsets of molecules characterized by different layer-averaged *N*_{4,L} values, it is evident that the loss of correlation shifts to longer times when *N*_{4,L} and, hence, the tetrahedral order increases. To quantify the dependence, the corresponding correlation times of the subsets *τ*_{L} are depicted as a function of *N*_{4,L}. We find that the rotational motion is nearly two orders of magnitude slower in the LDL-like regions with *N*_{4,L} > 0.7 than that in the HDL-like regions with *N*_{4,L} < 0.3. Consistently, a similar slowdown is observed when the layer-averaged density *ρ*_{L} decreases from HDL-like to LDL-like values. As a consequence of this substantial dependence of the dynamics on the structure, the ensemble averaged RCF shows a two-step signature, which is, however, smeared out to some extent because the interfacial regions lead to deviations from a perfect bimodality of the *N*_{4} distribution; see Fig. 8. This bimodal nature is quantified in the supplementary material. Comparing these results for *F*_{1}(*t*) with the decay of *C*_{N4}(*t*), we find that the molecular orientation changes faster than the layer-averaged structural identifier. Hence, the correlation times at low *N*_{4,L} (high *ρ*_{L}) and high *N*_{4,L} (low *ρ*_{L}), indeed, reflect the dynamical behaviors in extended and stable LDL-like and HDL-like regions, respectively, consistent with the observation that the dependence of *τ*_{L} on *N*_{4,L} or *ρ*_{L} levels off in these regions.

Finally, we investigate the temperature dependence of the rotational dynamics in HDL-like and LDL-like regions. Figure 10 shows the rotational correlation times *τ*_{L} in layers with various average structural identifiers *N*_{4,L} in an Arrhenius diagram. Independent of the value of *N*_{4,L}, we observe weak deviations from Arrhenius behavior across the critical temperature *T*_{c}, whereby the temperature dependence seems to be higher for larger *N*_{4,L} values. However, the time scale *τ*_{N4} of the *N*_{4,L} fluctuations becomes comparable to the correlation times *τ*_{L} for *N*_{4,L} > 0.7 when increasing the temperature toward *T*_{c}. This finding implies that the *τ*_{L} values for large *N*_{4,L} do no longer reflect the true LDL-related dynamics near the LLCP, because the correlation functions of these subsets decay prematurely when LDL-like environments are transformed into HDL-like environments with faster dynamics. In other words, the results for the LDL-related dynamics in the elongated system with *N* = 2000 are reliable only at sufficiently low temperatures, e.g., *T* = 125 K used in the above analysis.

Therefore, we increase the system size to *N* = 12 000 and the aspect ratio to 1:1:6 to further reduce the size and relevance of a potential HDL–LDL interface and, hence, to stabilize the phase separation. In Fig. 10, we see that the fluctuations of the layer-averaged *N*_{4,L} values are, indeed, slower for the larger than the smaller system. Moreover, the correlation times *τ*_{L} for high *N*_{4,L} values mildly increase, while those for small *N*_{4,L} remain essentially unchanged. On the other hand, long simulation times prevent us from exploring the temperature range well below the LLCP for *N* = 12 000. Thus, a reliable determination of the exact temperature dependence of the rotational dynamics in LDL-like regions is possible neither for the small nor for the large system. Still, the combined results for these systems indicate that distinct changes of LDL-like dynamics at *T*_{c} do not occur.

The correlation times *τ*_{L} extracted for *q* = 0.86 in regions with different layer-averaged densities *ρ*_{L} from the elongated *N* = 2000 system can be compared with those obtained along the isochores with the corresponding densities *ρ* from the cubic *N* = 2000 system. In Fig. 6, we observe that the temperature-dependent *τ*_{L} from the spatially resolved analysis of the “phase-separated” elongated system agrees with the correlation times *τ* extracted from the isochores for high densities *ρ* = *ρ*_{L}. This agreement confirms that these time constants characterize the dynamics in defined HDL-like regions. By contrast, the correlation times *τ*_{L} from the former analyses are smaller than the *τ* values from the latter for low densities *ρ* = *ρ*_{L} < *ρ*_{c} at temperatures above and near *T*_{c}. Again, this discrepancy can be attributed to the relatively fast HDL–LDL fluctuations in the elongated system at these temperatures. Accordingly, the *τ*_{L} values for low *ρ*_{L} tend to show a higher temperature dependence and approach the extrapolations of the isochore data with the corresponding densities *ρ* when the temperature is decreased well below the LLCP.

## IV. SUMMARY

We performed MD simulations to ascertain the thermodynamic, structural, and dynamical properties of various water-like models in the supercooled temperature range. Specifically, we started from the well established TIP4P/2005 model and systematically reduced the magnitude of the partial charges of the atoms and, thus, the strength of the hydrogen bonds and the degree of the tetrahedral order.^{61} This charge-scaling approach allowed us to investigate the role of local order for liquid polyamorphism and to determine the consequences of this interplay for the dynamical behavior. This knowledge is of use for all liquids with spacious low-energy and low-entropy structural motifs and, in particular, for water, for which a unique model is still lacking.

For such liquids with a spacious local order, the existence of a LLPT and a related LLCP was proposed to cause anomalous properties but is still the subject of a controversial scientific debate. Here, we exploited that charge scaling allows us to shift these phenomena to state points, which are conveniently accessible in state-of-the-art computer simulations. Similar results were reported for the WAC model of silica.^{60} For all studied charge-scaling factors, *q* = 0.86–0.91, we observed isochore crossing in the P–T diagram indicative of a LLPT. Despite its mean-field nature, the TSEOS formalism well describes the isochore data and enables the determination of the LLCP parameters for the water-like models. These results indicate that a LLCP exists for a significant range of interaction parameters. The finding that the phenomenon is not restricted to very specific particle interactions suggests that it can also occur for various real liquids, including water.

The TSEOS analysis showed that the critical temperature *T*_{c}, pressure *P*_{c}, and density *ρ*_{c} decrease when the magnitude of the partial charges is reduced, leading even to negative values of *P*_{c} at *q* < 0.9. While the critical parameters showed an approximately linear *q* dependence in the studied range *q* = 0.86–0.91, we expect a mildly nonlinear behavior in a broader range, because extrapolation to *q* = 1.00 underestimates the *T*_{c} and *P*_{c} values known from studies of the regular TIP4P/2005 model.^{33–35,37,86,87} When the charge-scaling factor and the hydrogen-bond strength are decreased, the isochores flatten out in the P–T diagram, suggesting that the liquid polyamorphism vanishes not too far below the covered *q* range. Studies on silica and silicon reported that the LLCP disappears below the liquid–vapor spinodal.^{31,60} Because the critical pressure takes strongly negative values for *q* = 0.86, one may speculate that the liquid–vapor spinodal is also close for water-like models with substantially reduced partial charges.

As for structure and dynamics, we found prominent changes when cooling in the vicinity of the LLCP. Specifically, distributions of local structure identifiers become bimodal, indicating a decomposition into HDL-like and LDL-like environments. Furthermore, the temperature dependence of the diffusion coefficients *D* and correlation times *τ* differ significantly under isochoric and isobaric conditions, consistent with previous results for the ST2 water model.^{74} Along the isochores, we observed moderate deviations from Arrhenius behavior, and the temperature dependence does not change even when crossing the LLCP for *ρ* ≈ *ρ*_{c}. Along the isobars, by contrast, the temperature dependence of *D* and *τ* is much higher near the LLCP than at *T* ≫ *T*_{c} and *T* ≪ *T*_{c}. These distinct differences between isochoric and isobaric conditions are confirmed when we use the temperature-dependent correlation times along the isochores to successfully predict those along the isobars. In addition, we observed that the dynamics of the water-like models slows down when the density or the pressure is decreased in the vicinity of the LLCP, consistent with the dynamical anomaly of water.

To ascertain the origin of the remarkable dynamical behaviors, we analyzed the respective dynamics in regions rich in HD or LD local structures, which span over a few nanometers and are stable for a few microseconds. For this purpose, we exploited that a decomposition into such regions occurs at sufficiently low temperatures in elongated systems, where the interface between these structurally different regions can be minimized along the smallest (*xy*) cross section of the simulation box. We characterized the variation of the structural and dynamical properties along the elongated *z* dimension by averaging over Δ*z* = 1 nm thick layers centered about the respective *z* value. In this way, we ensured that, except for the highest studied temperatures, the structural differences persist on the length and time scales on which the dynamical properties are probed. Contrarily, in previous approaches to the structure–dynamics relations of water,^{90,94} different subensembles were discriminated based on local structural identifiers, which fluctuate more rapidly in space and time so that partially averaged dynamical behaviors may be observed.

This spatially resolved analysis was restricted to the systems with the small scaling factors of *q* = 0.86 and *q* = 0.88, which exhibit relatively fast dynamics at the LLCP. It revealed that the rotational motion of the water-like molecules is about two orders of magnitude faster in HDL-like regions than that in LDL-like regions. Consistently, experimental work reported that the rotational correlation times differ by about this amount when heating high-density and low-density amorphous ice phases through their glass transition temperatures *T*_{g}.^{17} Considering that LDL-like configurations with slow dynamics are transformed into HDL-like ones with fast dynamics when the pressure is increased, these findings yield an explanation for the anomalous pressure dependence of the diffusion coefficients and correlation times. For neither the LDL-like nor the HDL-like regions, we found a change in the temperature dependence of the dynamics near the LLCP.

These findings yield insights into the nature of dynamical crossovers observed for water and other tetrahedral liquids near the LLCP upon cooling along isobars. This observation was often interpreted in terms of a FS transition related to a LLPT. Explicitly, one may argue that the HDL and LDL are very fragile and truly strong liquids, respectively, so that a FS transition occurs when the LLPT or the Widom line are crossed upon cooling. At variance with this conjecture, the present findings for water-like systems indicated that both HDL-like and LDL-like regions show only weak deviations from Arrhenius behavior. Thus, the dynamic crossover is not due to a transition between liquids with highly different fragilities. In other words, the HDL-like and LDL-like dynamic properties are observed only at a sufficient distance to the LLCP, and the effects in the vicinity of the LLCP cannot be attributed to either of the two liquid phases. Rather, the high temperature dependence near the LLCP is caused by the interconversion between fast HDL-like structures and slow LDL-like structures, i.e., it results from the crossover between the largely different *D*(*T*) or *τ*(*T*) curves of the respective regions. As a consequence, in the one-phase regime of the phase diagram, the temperature dependence along isobars becomes higher when the interconversion occurs more rapidly upon approaching *P*_{c}. While the dynamics in the NVT and NPT ensembles agree at the same state point, a clear dynamical crossover is not observed along the isochores where rapid HDL–LDL interconversions do not occur. This scenario has certain similarities with two-state models of water dynamics by Tanaka and co-workers^{91} and Caupin and co-workers,^{92,93} who proposed that the FS transition involves a crossover between water phases with dynamics on very different time scales.

## SUPPLEMENTARY MATERIAL

See the supplementary material for details on the TSEOS analysis, isochores for further charge-scaling factors, temperature-dependent densities along isobars at atmospheric pressure, pressure-dependent fractions of the low-density local structures, additional data for the temperature-dependent dynamics along isochores and isobars, a comparison of pressures for cubic and elongated simulation boxes, and spatially resolved analyses above *T*_{c}.

## ACKNOWLEDGMENTS

The authors thank the Deutsche Forschungsgemeinschaft (DFG) for funding through Forschungsgruppe FOR 1583 (Project no. Vo 905/9-2).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

Because ergodicity was not fully reached for the layer-averaged quantities at the lowest studied temperatures, it proved useful to calculate $\u27e8\xi L\u27e92$ as ⟨ξ_{L}(*t* + *t*_{0})⟩⟨ξ_{L}(*t*_{0})⟩ and to normalize the correlation functions by C(*t* = 0). In the ergodic limit, this is equivalent to Eq. (10).