Glass Polymorphism in TIP4P/2005 Water: A Description Based on the Potential Energy Landscape Formalism

TL;DR

This study uses the potential energy landscape formalism to analyze pressure-induced transformations between amorphous ice forms in water, revealing signatures of first-order phase transition behavior and identifying distinct PEL regions for LDA and HDA.

Contribution

It applies the PEL formalism to TIP4P/2005 water simulations, uncovering anomalous landscape features associated with LDA-HDA transformations and linking simulation results to experimental conditions.

Findings

PEL properties show anomalous behavior during LDA-HDA transition.

IS pressure exhibits van der Waals-like loop.

Transformation pressure extrapolates to experimental values.

Abstract

The potential energy landscape (PEL) formalism is a statistical mechanical approach to describe supercooled liquids and glasses. Here we use the PEL formalism to study the pressure-induced transformations between low-density amorphous ice (LDA) and high-density amorphous ice (HDA) using computer simulations of the TIP4P/2005 molecular model of water. We find that the properties of the PEL sampled by the system during the LDA-HDA transformation exhibit anomalous behavior. In particular, at conditions where the change in density during the LDA-HDA transformation is approximately discontinuous, reminiscent of a first-order phase transition, we find that (i) the inherent structure (IS) energy, $e_\text{IS}(V)$, is a concave function of the volume, and (ii) the IS pressure, $P_\text{IS}(V)$, exhibits a van der Waals-like loop. In addition, the curvature of the PEL at the IS is anomalous, a non-monotonic function of $V$. In agreement with previous studies, our work suggests that conditions (i) and (ii) are necessary (but not sufficient) signatures of the PEL for the LDA-HDA transformation to be reminiscent of a first-order phase transition. We also find that one can identify two different regions of the PEL, one associated to LDA and another to HDA. Our computer simulations are performed using a wide range of compression/decompression and cooling rates. In particular, our slowest cooling rate (0.01 K/ns) is within the experimental rates employed in hyperquenching experiments to produce LDA. Interestingly, the LDA-HDA transformation pressure that we obtain at $T=80$ K and at different rates extrapolates remarkably well to the corresponding experimental pressure.

Tables

Table 1. Table 1: Combinations of compression/decompression temperature T 𝑇 T , cooling-rate q c subscript 𝑞 c q_{\text{c}} , and compression/decompression rate q p subscript 𝑞 p q_{\text{p}} studied here. For the case T = 280 𝑇 280 T=280  K the liquid was equilibrated at 280 K and 1 bar before compression/decompression. For the case q c = 0.01 subscript 𝑞 c 0.01 q_{\text{c}}=0.01  K/ns some cooling runs were started from the liquid equilibrated at 200 K and others from 240 K (all at 1 bar). For all other cases the liquid was initially equilibrated at 240 K and 1 bar. The slowest compressions, at q p = 0.1 subscript 𝑞 p 0.1 q_{\text{p}}=0.1  MPa/ns, were started from the configurations obtained at 400 MPa during the compression runs at q p = 1 subscript 𝑞 p 1 q_{\text{p}}=1  MPa/ns (decompressions at the slowest rate considered are not included).

Changing Parameter	Constant Parameters
$T=20$, 80, 160, 200, 280 K	$q_{\text{c}}=30$ K/ns	$q_{\text{p}}=300$ MPa/ns
$q_{\text{c}}=0.01$, 0.1, 1, 30, 100 K/ns	$T=80$ K	$q_{\text{p}}=10$ MPa/ns
$q_{\text{p}}=0.1$, 1, 10, 300, 1000 MPa/ns	$T=80$ K	$q_c=0.1$ K/ns

Equations

$$Z(T,V)=\int \Omega (e_{\text{IS}})\text{d}{e}_{\text{IS}}{\text{e}}^{-\beta F_{\text{basin}}(e_{\text{IS}},T,V)},$$

$$F_{\text{basin}}(e_{\text{IS}},T,V)=e_{\text{IS}}+F_{\text{vib}}(e_{\text{IS}},T,V),$$

$$\beta F_{\text{vib}}(e_{\text{IS}},T,V)\approx {⟨\sum _{i=1}^{6N-3}\ln \left(\beta \hslash {\omega }_i(e_{\text{IS}},V)\right)⟩}_{e_{\text{IS}}}.$$

$$\beta F_{\text{vib}}(e_{\text{IS}},T,V)=(6N-3)\ln \left(\beta A_0\right)+\mathcal{S}(e_{\text{IS}},V),$$

$$\mathcal{S}(e_{\text{IS}},V)={⟨\sum _{i=1}^{6N-3}\ln \left(\frac{\hslash {\omega }_i(e_{\text{IS}},V)}{A_0}\right)⟩}_{e_{\text{IS}}}.$$

$$F=F_{\text{IS}}(E_{\text{IS}},T,V)+F_{\text{vib}}(E_{\text{IS}},T,V),$$

$$F_{\text{IS}}(E_{\text{IS}},T,V)=E_{\text{IS}}(N,T,V)-T{S}_{\text{conf}}(E_{\text{IS}}).$$

$$S_{\text{conf}}(E_{\text{IS}})\equiv k_{\text{B}}\ln [\Omega (E_{\text{IS}})\text{d}{E}_{\text{IS}}].$$

$$\frac{{\partial }^{2}F}{\partial V^2}>0.$$

$$\frac{\partial P}{\partial V}<0.$$

$$\frac{{\partial }^2{F}_{\text{IS}}}{\partial V^2}+\frac{{\partial }^2{F}_{\text{vib}}}{\partial V^2}>0.$$

$$\frac{\partial P_{\text{IS}}}{\partial V}+\frac{P_{\text{vib}}}{\partial V}<0,$$

$$P_{\text{IS}}$$

$$P_{\text{vib}}$$

$${\Delta }_P={-\frac{\partial P}{\partial V_{\text{m}}}|}_{V_{\text{m,mid}}}$$

$${\Delta }_{P_{\text{IS}}}={-\frac{\partial P_{\text{IS}}}{\partial V_{\text{m}}}|}_{V_{\text{m,mid}}}.$$

$${\Delta }_{e_{\text{IS}}}=\min \left[\frac{{\partial }^2{e}_{\text{IS}}}{\partial V_{\text{m}}^2}\right],$$

Full text

††preprint: APS/123-QED

I Introduction

Water is a prototypical complex substance. This complexity is manifested in numerous anomalous properties present in the liquid state Errington and Debenedetti (2001); Ludwig (2001); Debenedetti (2003), such as the maximum in density upon isobaric cooling ($4\text{\,}\mathrm{\SIUnitSymbolCelsius}$ at 1 bar), and the maximum in diffusivity upon isothermal compression ($\approx 200$ MPa at $0\text{\,}\mathrm{\SIUnitSymbolCelsius}$). In the solid state, water can exist in a surprisingly large number of crystalline polymorphs (17 distinct ices have been identified so far Salzmann (2019)) and in at least two different glassy states (amorphous ices) Mishima and Stanley (1998); Angell (2004); Loerting and Giovambattista (2006); Loerting et al. (2011); Handle et al. (2017). The most common forms of glassy water are low- (LDA) and high-density (HDA) amorphous ice. Amorphous ices can be obtained utilizing several thermodynamic paths Burton and Oliver (1935); Brüggeller and Mayer (1980); Mishima et al. (1984); Handle and Loerting (2015) and the behavior of these systems is well documented Mishima et al. (1985); Mishima (1994); Loerting et al. (2001); Gromnitskaya et al. (2001); Klotz et al. (2005); Loerting et al. (2006); Winkel et al. (2008); Handle and Loerting (2018a, b). Remarkably, most experimental studies indicate that, when properly annealed Perakis et al. (2017); Handle and Loerting (2018a), LDA and HDA can be interconverted by sharp and reversible transformations reminiscent of first-order phase transitions between equilibrium states Mishima et al. (1985); Mishima (1994); Gromnitskaya et al. (2001); Klotz et al. (2005); Koza et al. (2005); Winkel et al. (2008); Handle and Loerting (2018a).

The puzzling behavior of water has spawned several potential theoretical scenarios Angell (2008); Gallo et al. (2016); Handle et al. (2017); Anisimov et al. (2018). According to the liquid-liquid critical point (LLCP) scenario, water is hypothesized to exist in two different liquid states at low temperature, a low- (LDL) and a high-density (HDL) liquid. In addition, LDL and HDL are separated by a first-order phase transition line ending at an LLCP at higher temperatures  Poole et al. (1992). One relevant advantage of the LLCP scenario, relative to other available theoretical explanations Angell (2008); Gallo et al. (2016); Handle et al. (2017); Anisimov et al. (2018), is that the LLCP scenario naturally rationalizes the experimental phenomenology present in the amorphous ices. Specifically, in the LLCP scenario, LDA and HDA are the glass counterparts of LDL and HDL, respectively. Accordingly, the sharp LDA$\rightleftarrows$HDA transformation is a result of extending the liquid-liquid phase transition into the glass domain Poole et al. (1992); Mishima and Stanley (1998), explaining the sharpness of the LDA$\rightleftarrows$HDA transformation found in experiments. Experimental evidence for the connection between LDA and LDL, and between HDA and HDL, can be found in Refs. McMillan and Los, 1965; Johari et al., 1987; Smith and Kay, 1999; Handle et al., 2012; Amann-Winkel et al., 2013; Perakis et al., 2017. However, the true nature of LDA and HDA Tse et al. (1999); Johari (2000) and their relationship with the liquid state are still a matter of debate Johari (2014, 2015); Stern et al. (2015); Handle et al. (2016); Shephard and Salzmann (2016); Fuentes-Landete et al. (2019); Stern et al. (2019).

In this work, we study the LDA$\rightleftarrows$HDA transformation in water using molecular dynamics (MD) simulations in conjunction with the potential energy landscape (PEL) approach Goldstein (1969); Stillinger and Weber (1982); Debenedetti and Stillinger (2001); Sciortino (2005); Stillinger (2015). The PEL approach is a powerful theoretical framework within statistical mechanics that has been used extensively to study the dynamic and thermodynamic behavior of liquids at low temperatures Sastry (2001); Mossa et al. (2002); La Nave et al. (2006); Heuer (2008), including water Sciortino et al. (2003); La Nave and Sciortino (2004); Handle and Sciortino (2018a, b). In particular, it allows one to express the Helmholtz free energy $F(N,T,V)$ of a liquid [and hence, the corresponding equation of state (EOS)] in terms of statistical properties of the PEL. In the case of water, the PEL formalism has been successfully applied to obtain the EOS for the SPC/E Berendsen et al. (1987) and TIP4P/2005 Abascal and Vega (2005) water models Sciortino et al. (2003); Handle and Sciortino (2018a). Such an EOS can be used to extrapolate the behavior of the supercooled liquid to low temperatures. Interestingly, the PEL-EOS for TIP4P/2005 predicts the existence of an LLCP at $T_{\text{c}}=175$ K, $P_{\text{c}}=175$ MPa, and $\rho_{\text{c}}=0.997$ g/cm3 Handle and Sciortino (2018a), consistent with other predictions for this model Abascal and Vega (2010); Sumi and Sekino (2013); Singh et al. (2016); Biddle et al. (2017). In the case of SPC/E water, the LLCP is estimated to be located below the Kauzman temperature 111 Within the PEL formalism, the Kauzman temperature $T_{\text{K}}$ is defined as the temperature where the liquid has access to only one basin, i.e., $S_{\text{conf}}=0$ Sciortino (2005). This means that below $T_{\text{K}}$ the system cannot undergo structural changes and hence, it cannot show a phase separation.

and hence, it is not accessible to the liquid state Scala et al. (2000a, b); Sciortino et al. (2003).

Besides liquids, the PEL approach has also been applied to study several atomic and molecular glasses Shell et al. (2003); Shell and Debenedetti (2004); Stillinger (2015). In particular, it was used to study the LDA$\rightleftarrows$HDA transformation in SPC/E and ST2 Stillinger and Rahman (1974) water Giovambattista et al. (2003, 2016, 2017). In the case of SPC/E water, where an LLCP is not accessible Scala et al. (2000a, b); Sciortino et al. (2003), the changes in the PEL properties sampled by the system during the LDA$\rightarrow$HDA transformation are rather smooth and change monotonically with density Giovambattista et al. (2003). These are the expected results for normal glasses, such as for a system of soft-spheres Sun et al. (2018); Sciortino (2005). Instead, in the case of ST2 water, where the LLCP is accessible Poole et al. (1992, 2005); Cuthbertson and Poole (2011); Liu et al. (2012); Palmer et al. (2014); Smallenburg and Sciortino (2015); Palmer et al. (2018a, b), the PEL properties sampled by the system during the LDA$\rightarrow$HDA transformation exhibit anomalous behavior consistent with a first-order like phase transition between the two glass states Giovambattista et al. (2016, 2017). This conclusion was also supported by a PEL study of a water-like monatomic model that exhibits liquid and glass polymorphism Sun et al. (2018).

One limitation of the ST2 water model to study glassy water is its inability to reproduce the structure of HDA Chiu et al. (2013, 2014). Instead, the TIP4P/2005 water model reproduces relatively well the structure of LDA and HDA Wong et al. (2015). Thus, it is a natural question whether TIP4P/2005 water also exhibits PEL anomalies during the pressure-induced LDA$\rightarrow$HDA transformation as found in ST2 water. In this work we address this question and study the PEL of TIP4P/2005 water during the pressure-induced LDA$\rightleftarrows$HDA transformation. We note that the TIP4P/2005 model is presently regarded as one of the most realistic (rigid) models to study liquid and crystalline water Vega and Abascal (2011). Moreover, several studies Abascal and Vega (2010); Wikfeldt et al. (2011); Sumi and Sekino (2013); Yagasaki et al. (2014); Russo and Tanaka (2014); Singh et al. (2016); Biddle et al. (2017); Handle and Sciortino (2018a) are consistent with the presence of an LLCP in TIP4P/2005 water. TIP4P/2005 water also displays an apparent first-order transition between LDA and HDA Wong et al. (2015), as observed experimentally.

A peculiar property of glasses is that their properties depend on the preparation process considered, i.e., glasses are history-dependent materials Binder and Kob (2011). This implies that the LDA$\rightleftarrows$HDA transformation can be sensitive to the cooling and compression rates employed as was shown in MD simulation studies Giovambattista et al. (2003, 2016, 2017). Accordingly, in this work we pay particular attention to the effects of cooling ($q_{\text{c}}$) and/or compression ($q_{\text{P}}$) rates. Specifically, we are able to reach, for the first time, cooling rates as slow as $q_{\text{c}}=0.01$ K/ns which is comparable to experimental cooling rates necessary to avoid crystallization in hyperquenching techniques Brüggeller and Mayer (1980); Dubochet and McDowall (1981); Mayer and Brüggeller (1982); Mayer (1985); Kohl et al. (2000). The slowest compression rate employed here ($q_{\text{P}}=0.1$ MPa/ns) expands beyond the slowest compression/decompression rates studied so far in MD simulations but it is still about three orders of magnitude faster than the fastest experimental rate we are aware of ($$65\text{,}\mathrm{GPa}\text{,}{\mathrm{s}}^{-1}$=6.5\cdot 10^{-5}$ MPa/ns) Chen and Yoo (2011) and more than seven orders of magnitude faster than the rates commonly used experimentally ($\approx 10^{-1}$–101 MPa/s) Mishima (1994); Winkel et al. (2008); Handle and Loerting (2018a).

The structure of this work is as follows. In Sec. II, we discuss briefly the PEL formalism and its application to the study of liquids. In Sec. III we describe the computer simulation details and the numerical methods employed. The results are presented in Sec. IV where we discuss the PEL properties of TIP4P/2005 water during the preparation of LDA (Sec. IV.1) and during the pressure-induced LDA$\rightleftarrows$HDA transformation (Sec. IV.2). A summary and discussion is included in Sec. V.

II The PEL formalism for Liquids

The PEL approach, as introduced by Stillinger and Weber Stillinger and Weber (1982), is a powerful tool to describe the properties of liquids. For a system of $N$ rigid water molecules, the PEL is a hypersurface embedded in a $(6N+1)$-dimensional space. It is defined by the potential energy of the system, $U(\vec{r}^{N},\phi^{N},\theta^{N},\psi^{N})$, as function of the $3N$ coordinates of the molecules’ center of mass $\vec{r}$, and the corresponding $3N$ Euler angles $\phi$, $\theta$, $\psi$. At any given time $t$, the system is represented by a single point on the PEL with coordinates given by the values of $\vec{r}^{N}(t),\phi^{N}(t),\theta^{N}(t),\psi^{N}(t)$. It follows that, as time evolves, the representative point of the system describes a trajectory on the PEL, as it moves from one basin of the PEL to another. A basin is defined as the set of points of the PEL that lead to the same local minimum by potential energy minimization. The local minimum associated to a given basin is called inherent structure (IS) and its associated energy is denoted $e_{\text{IS}}$. Depending on the temperature considered, the representative point of the system may or may not be able to overcome potential energy barriers separating different basins. Accordingly, different regions of the PEL may be accessible to the system depending on the temperature considered.

In the PEL framework, the canonical partition function can be formulated as a one-dimensional integral Debenedetti and Stillinger (2001); Sciortino (2005)

[TABLE]

where $\Omega(e_{\text{IS}})\text{d}e_{\text{IS}}$ is the number of basins with IS energy between $e_{\text{IS}}$ and $e_{\text{IS}}+\text{d}e_{\text{IS}}$, $F_{\text{basin}}(e_{\text{IS}},T,V)$ is the average basin free energy of basins with IS energy $e_{\text{IS}}$, $\beta=1/k_{\text{B}}T$ and $k_{\text{B}}$ is Boltzmann’s constant. All basins containing a significant amount of crystalline order are by definition excluded in the integration over phase space in Eq. 1 Sciortino (2005). The basin free energy can further be written as

[TABLE]

where $F_{\text{vib}}$ accounts for the vibrational motion of the system around the IS with energy $e_{\text{IS}}$. If we assume that the PEL around an IS can be approximated by a quadratic function (harmonic approximation), we can write

[TABLE]

Here, the values $\omega_{i}(e_{\text{IS}},V)$ are the $6N-3$ normal mode frequencies and $\hbar$ is Planck’s constant in its reduced form. To separate the $T$ and $e_{\text{IS}}$ dependence of $F_{\text{vib}}$, we write

[TABLE]

where

[TABLE]

The latter is called the basin shape function and it quantifies the average local curvature of the PEL around the IS. In Eq. 5, the average is taken over all IS with energy $e_{\text{IS}}$ and $A_{0}=1$  kJ/mol is a constant that ensures the arguments of the logarithm to have no units.

It can further be shown that the system free energy $F(N,T,V)$ can be expressed as

[TABLE]

where

[TABLE]

Here $E_{\text{IS}}=\left<e_{\text{IS}}(N,T,V)\right>$ is the average energy of the IS sampled by the system and $S_{\text{conf}}$ is the configurational entropy. The latter is defined as

[TABLE]

In summary, the PEL formalism allows the free energy of the system $F(N,T,V)$ to be expressed in terms of three basic properties of the PEL at constant $N,T,V$ Sciortino (2005); Sastry (2001); Debenedetti and Stillinger (2001):

• (i)

the average energy of the IS sampled by the system, i.e., $E_{\text{IS}}$;

• (ii)

the number of IS with energy between $E_{\text{IS}}$ and $E_{\text{IS}}+\text{d}E_{\text{IS}}$, i.e., $\Omega(E_{\text{IS}})\text{d}E_{\text{IS}}$;

• (iii)

the average curvature of the PEL at the IS, as quantified by $\mathcal{S}$.

For systems at constant $N,T,V$, thermodynamic arguments show that the system is in stable or metastable equilibrium if and only if Stanley (1971)

[TABLE]

Here and in the following, the partial derivatives are evaluated at constant $T$ and $N$. Eq. 9 states that $F(N,T,V)$ must be a convex function of $V$ along an isotherm at constant $N$. Alternatively, since $P=-\partial F/\partial V$, Eq. 9 can be rewritten as

[TABLE]

Eq. 10 implies that the isothermal compressibility of the system must be positive. If Eqs. 9 or 10 are violated then the system is unstable and exhibits a phase transition.

Within the PEL formalism for liquids, Eq. 9 can be re-written in terms of Eq. 6 as

[TABLE]

Similarly, Eq. 10 can be re-written as

[TABLE]

where $P=P_{\text{IS}}+P_{\text{vib}}$ and

[TABLE]

Eq. 11 shows that, for a system at constant $T$ and $N$, a phase transition may occur due to a concavity in $F_{\text{IS}}(V)$, $F_{\text{vib}}(V)$, or both. Alternatively, Eq. 12 implies that a phase transition can occur due to a positive slope in $P_{\text{IS}}(V)$, $P_{\text{vib}}(V)$, or both.

Eqs. 11 and 12 are thermodynamic stability conditions, in terms of PEL properties, that apply only to equilibrium systems. Accordingly, they are not suitable to define phase transitions between glasses. However, for the case of two different polyamorphic systems Giovambattista et al. (2016, 2017); Sun et al. (2018), it was shown that during the first-order-like phase transition between LDA and HDA, $\partial^{2}E_{\text{IS}}/\partial V^{2}>0$ and $\partial P_{\text{IS}}/\partial V<0$. For a system in equilibrium, this could lead to a violation of Eqs. 11 and 12. Indeed, if the basin shape function is approximately independent of the volume of the system then it can be shown that, within the harmonic approximation, $\partial P_{\text{vib}}/\partial V\approx 0$ and $\partial F_{\text{vib}}/\partial V\approx 0$. In this case, the liquid exhibits a phase transition if and only if $\partial P_{\text{IS}}/\partial V<0$, or alternatively, if and only if $F_{\text{IS}}$ is a concave function of $V$ (at constant $N$ and $T$).

III Simulation Details

The basis for our study of the LDA$\rightleftarrows$HDA transformation in TIP4P/2005 water are MD simulations starting from LDA configurations. We first prepare LDA at 1 bar by quenching the equilibrium liquid from temperature $T_{0}$ down to $T\leq 80$ K, using cooling rates $q_{\text{c}}$ in the range $0.01$–100 K/ns (step (i); vertical arrow in Fig. 1). $T_{0}=240$ K for $q_{\text{c}}\geq 0.1$ K/ns and $T_{0}=200$ or 240 K for $q_{\text{c}}=0.01$ K/ns. During this cooling from the liquid, we save configurations at different intermediate $T$. Configurations so prepared are then compressed up to $P=3$ GPa to produce HDA, using compression rates $q_{\text{P}}$ in the range 1–1000 MPa/ns (step (ii), horizontal right arrows in Fig. 1). Compressions at the slowest rate, $q_{\text{P}}=0.1$ MPa/ns, are started from $P=400$ MPa, using the respective configurations obtained during compression with rate $q_{\text{P}}=1$ MPa/ns as initial configurations. HDA is then decompressed at the same rate starting from 2 GPa (step (iii), horizontal left arrows in Fig. 1).

MD simulations to study glassy water have been criticized in the past due to the fast cooling and compression rates employed. Accordingly, in this work, we explore in detail the influence of the compression/decompression temperature $T$ as well as rates $q_{\text{c}}$ and $q_{\text{P}}$ on the LDA$\rightleftarrows$HDA transformation observed in our MD simulations. All combinations of $T$, $q_{\text{c}}$, and $q_{\text{P}}$ studied are listed in Table 1. We stress that the smallest cooling rate studied here ($q_{\text{c}}=0.01$ K/ns) is comparable to the rates used in hyperquenching experiments of liquid water Brüggeller and Mayer (1980); Dubochet and McDowall (1981); Mayer and Brüggeller (1982); Mayer (1985); Kohl et al. (2000). With such a slow cooling rate, we need to simulate 16 $\mu$s to cool the system from $240$ to $80$ K.

The systems studied consist of $N=1728$ TIP4P/2005 water molecules in a cubic box with periodic boundary conditions. All our MD simulations are performed at constant $N$, $P$, $T$ using the GROMACS 5.1.4 and 2016.5 Van Der Spoel et al. (2005) simulation packages. Simulations use the leap-frog integrator with a time step of 2 fs. Temperature is controlled using a Nosé-Hoover thermostat Nosé (1984); Hoover (1985) and pressure is controlled using a Parinello-Rahman barostat Parrinello and Rahman (1981). For the Coulomb interactions, we use a particle mesh Ewald treatment Essmann et al. (1995) with a Fourier spacing of 0.1 nm. For both the Lennard-Jones (LJ) and the real space Coulomb interactions, a cutoff of 0.85 nm is used. Lennard-Jones interactions beyond 0.85 nm have been included assuming a uniform fluid density. Water molecules are treated as rigid by using the LINCS (Linear Constraint Solver) algorithm Hess (2008) of 6th order with one iteration to correct for rotational lengthening.

At $T=240$ K the system was simulated for 10 ns. After an initial equilibration period of 1 ns we extracted ten independent configurations separated by 1 ns. These ten configurations served as the starting points for ten independent simulations [steps (i) to (iii)] for every set of compression/decompression temperature $T$ ($\leq 200$ K), $q_{\text{c}}$, and $q_{\text{p}}$. The length of these cooling and compression/decompression runs is determined by $q_{\text{c}}$ and $q_{\text{P}}$, respectively.

In order to calculate the properties of the PEL, i.e., $e_{\text{IS}}$, $P_{\text{IS}}$, ${\mathcal{S}}$, throughout the LDA$\rightleftarrows$HDA transformation, we obtain the IS during both the compression and decompression runs by minimizing the potential energy of the system every 10 MPa. The minimization directly yields $e_{\text{IS}}$ and the Virial Frenkel and Smit (2001) at the IS is used to calculate $P_{\text{IS}}$. The basin shape function $\mathcal{S}$ is obtained from Eq. 5 using the normal mode frequencies given by the eigenvalues of the Hessian at the IS.

IV Results

IV.1 Preparation of LDA

As discussed previously, to generate LDA, the system is equilibrated at $T_{0}=240$ K and $P=1$ bar, and then cooled isobarically to different temperatures $20\leq T\leq 200$ K, using different cooling rates $q_{\text{c}}$. The behavior of the density of the system upon cooling is shown in Fig. 2. In all cases the density first decreases, reaches a minimum around $200$ K and then increases linearly. The influence of the cooling rate is clearly visible in Fig. 2, which shows that the LDA samples obtained after cooling are less dense as $q_{\text{c}}$ is decreased. As is typical for glasses, fast cooling increases the glass transition temperature, leaving the system trapped in higher density and higher energy regions of the PEL. The very similar slope displayed by all samples in the low $T$ part (i.e., below the density minimum) results from the decrease in vibrations around the IS of the basins the systems are trapped in.

We note that the densities of the LDA forms obtained with $q_{\text{c}}=0.01$ K/ns (i.e., the experimental rate) are very similar to the density of TIP4P/2005 ice I${}_{\text{h}}$. This suggests that the LDA so produced consists of an almost fully developed, highly tetrahedral hydrogen-bond network. Consistent with this assessment is the structure factor of LDA at 80 K reported in Fig. 3. It is visible that the pre-peak ($Q<20$ nm*-1*) grows as $q_{\text{c}}$ is decreased. At the same time it moves to lower $Q$ and separates more clearly from the main peak ($Q\approx 30$ nm*-1*). This also indicates that lower $q_{\text{c}}$ yield an LDA with a more developed tetrahedral hydrogen-bond network. The main peak in $S(Q)$ grows only slightly as $q_{\text{c}}$ is lowered and the features at higher $Q$ coincide within the noise of our data for all cooling rates studied. The structural data in Fig. 3 indicate further, that our samples have not crystallized during cooling. We confirm this by calculating the local order parameter as defined in Ref. Russo et al., 2014 and find that $>99\%$ of the molecules are classified as liquid (cf. also Refs. Engstler and Giovambattista, 2017; Martelli et al., 2018).

IV.2 Pressure-Induced LDA-HDA Transformation

Next, we discuss the properties of the system during the pressure induced LDA$\rightleftarrows$HDA transformation. In Sec. IV.2.1, we study the effects of varying $T$ at constant ($q_{\text{c}}$, $q_{\text{P}}$). The effects of varying $q_{\text{c}}$ at constant ($q_{\text{P}}$, $T$), and $q_{\text{P}}$ at constant ($q_{\text{c}}$, $T$) are addressed in Secs. IV.2.2 and IV.2.3, respectively.

IV.2.1 Temperature Dependence

Samples of LDA obtained at $T=20$, 80, 160, 200, and 280 K and $P=0.1$ MPa, using a cooling rate $q_{\text{c}}=30$ K/ns, were compressed at constant temperature from $P=0.1$ MPa to 3 GPa with $q_{\text{P}}=300$ MPa/ns (see Fig. 4). At $T\leq 160$ K, i.e., below the estimated LLCP temperature $T_{\text{c}}\approx 175-193$ K Abascal and Vega (2010); Sumi and Sekino (2013); Singh et al. (2016); Biddle et al. (2017); Handle and Sciortino (2018a), the system exhibits a sharp increase in density which signals the LDA$\rightarrow$HDA transformation. This density increase is sharper and shifts towards lower pressures as $T\rightarrow T_{\text{c}}$. At 280 K ($T>T_{\text{c}}$), the system is in the liquid state and the density increases smoothly and monotonically during compression (inset, purple lines). The behavior of $\rho(P)$ also overlaps with the corresponding equilibrium liquid isotherm obtained from the PEL-EOS in Ref. Handle and Sciortino, 2018a (inset, solid orange line). At 200 K, however, $\rho(P)$ (inset, blue lines) does not follow the corresponding PEL-EOS isotherm (inset, dashed orange line) and also shows a relatively sharp density step although no LDA$\rightarrow$HDA transformation is expected ($T>T_{\text{c}}$). This indicates that the used compression rate is too large relative to the relaxation time of the liquid at 200 K and the system cannot reach equilibrium during compression.

The HDA configurations obtained at 2 GPa were used as the starting configurations for the decompression runs conducted at the same $T$ and rate $q_{\text{P}}$ (see Fig. 4). At very low and negative pressures, all amorphous ices fracture. This corresponds to the sudden (almost vertical) density drop at $P_{\text{fract}}<-500$ MPa in Fig. 4. Interestingly, $P_{\text{fract}}$ becomes more negative as $T$ decreases implying that the amorphous ices are stronger under tension at lower $T$.

The HDA$\rightarrow$LDA transformation, at $P_{\text{fract}}<P<0$, is very weak for the present model and at the studied rates. It is observable at 160 K ($T<T_{\text{c}}$) where $\rho(P)$ exhibits a change of slope at $\approx-85$ MPa and an LDA-like state can be identified at $-600<P<-300$ MPa. However, at $T\leq 80$ K, the HDA state seems to expand continuously until it finally fractures (cf. Refs Wong et al., 2015; Engstler and Giovambattista, 2017). We also note that, a density step during the decompression of HDA is visible at $T=200$ K in Fig. 4. However, $$200\text{,}\mathrm{K}$>T_{\text{c}}$ and hence, no HDA$\rightarrow$LDA transformation can exist at this temperature. As mentioned previously, for the compression rate employed, the system is not able to constantly accommodate to the change in pressure at $T=200$ K, since the equilibration time is too long relative to the $q_{\text{P}}$ employed. Accordingly, the densities during compression and decompression do not coincide with each other and also deviate from the equilibrium liquid isotherm obtained from the PEL-EOS (inset, dashed orange line) Handle and Sciortino (2018a). We note that the results shown in Fig. 4 are in full agreement with previous studies of glassy TIP4P/2005 water Wong et al. (2015); Engstler and Giovambattista (2017) and in qualitative agreement with glassy ST2 water Poole et al. (1992); Giovambattista et al. (2016, 2017).

At $T=280$ K the system is in the liquid state and $\rho(P)$ during the decompression coincides with the density during compression. In addition, $\rho(P)$ coincides with the density of the equilibrium liquid obtained from the PEL-EOS in Ref. Handle and Sciortino, 2018a (inset, solid orange line).

Fig. 5 shows that the effect of increasing $T$ is to shift (i) the LDA$\rightarrow$HDA transformation (blue symbols) to lower pressures, and (ii) the HDA$\rightarrow$LDA (red symbols) transformation to higher pressures. This results in a narrower hysteresis during the LDA$\rightleftarrows$HDA transformation as $T\rightarrow T_{\text{c}}$. We stress that these results are qualitatively consistent with experiments (see purple lines). However, the LDA$\rightleftarrows$HDA transformation pressures obtained in our MD simulations (at the employed rates) are off relative to experimental values. As we will show in Secs. IV.2.2 and  IV.2.3, lowering $q_{\text{c}}$ and/or $q_{\text{P}}$ reduces the hysteresis, improving the MD results relative to experiments. We also note that the pressure-interval associated to the LDA$\rightarrow$HDA transformation (blue error bars) shrinks with increasing temperatures (i.e., the transformation becomes sharper). We expect a similar $T$-effect on the HDA$\rightarrow$LDA transformation. However, the HDA$\rightarrow$LDA transformation is not clearly observable in Fig. 4 at low temperatures. Accordingly, in Fig. 5, the HDA$\rightarrow$LDA transformation is indicated only for $T=160$ K. We also report the density steps for 200 K, although they are clearly not related to HDA$\rightleftarrows$LDA transformations.

Included in Fig. 5 are also the liquid spinodal lines associated to the LLCP predicted by the PEL-EOS Handle and Sciortino (2018a). Within the LLCP scenario, the LDA$\rightleftarrows$HDA transformations are nothing else but the extensions of the LDL$\rightleftarrows$HDL spinodal lines into the glass domain (cf., e.g., Ref. Mishima and Stanley, 1998). It follows from Fig. 5 that the liquid-liquid spinodal lines predicted by the PEL-EOS cannot be used to estimate the LDA$\rightleftarrows$HDA transformation lines obtained either from MD simulations nor experiments (see purple lines). This discrepancy is not due to a deficiency of the PEL formalism, but due to the fact that the PEL-EOS is parameterized based on IS sampled by the equilibrium liquid at $200<T<270$ K Handle and Sciortino (2018a). As we show in the supplementary material, the IS sampled by the system during compression at low $T$ are not explored by the equilibrium liquid. It is probably this difference between the PEL regions sampled by the liquid and the glass that makes the extension of the PEL-EOS into the glass domain of very limited applicability.

Next, we discuss the PEL properties of the system corresponding to the runs shown in Fig. 4. Fig. 6 shows $P_{\text{IS}}$, $e_{\text{IS}}$, and ${\mathcal{S}}$ as a function of density during the compression/decompression processes. Consistent with studies based on the SPC/E and ST2 models Giovambattista et al. (2003, 2016, 2017), we find that, at low temperatures, all PEL properties are different along the compression and decompression paths. This implies that the system explores different regions of the PEL during the LDA$\rightarrow$HDA and HDA$\rightarrow$LDA transformations. This is not the case for $T=280$ K, because, as explained previously, the system is in the equilibrium liquid state along both the compression and decompression processes.

The behavior of $P_{\text{IS}}(\rho)$ shown in Fig. 6(a) is rather complex. During compression at $T=280$ K, in the equilibrium liquid state, $P_{\text{IS}}(\rho)$ is a monotonic function of density, as expected. However, at $T=160$ K ($T<T_{\text{c}}$) and even at $T=200$ K ($T>T_{\text{c}}$), $P_{\text{IS}}(\rho)$ exhibits a van der Waals-like loop [i.e., a section of negative slope in $P_{\text{IS}}(\rho)$]. At these temperatures, the behavior of $P_{\text{IS}}(\rho)$ shown in Fig. 6(a) is reminiscent of the behavior expected for equilibrium systems during a phase transition (cf. Sec. II). Interestingly, there is no van der Walls-like loop at very low temperatures (80 and 20 K).

It follows from Figs. 4 and 5 that, as $T$ increases from 20 to 160 K, the LDA$\rightarrow$HDA transformation becomes sharper, more reminiscent of a first-order phase transition, and accordingly, $P_{\text{IS}}(\rho)$ develops a van der Waals-like loop [see Fig. 6(a)]. To make this point clear, we compare the slopes of $P(\rho)$ and $P_{\text{IS}}(\rho)$ at the mid-point of the LDA$\rightarrow$HDA transformation in Fig. 7. We use the following notation (cf. Ref. Giovambattista et al., 2017):

[TABLE]

and

[TABLE]

Here $V_{\text{m}}$ denotes the molar volume and $V_{\text{m,mid}}$ denotes the molar volume at the midpoint of the transformation. We note that $\Delta_{P}=0$ corresponds to a discontinuous density jump. Hence, the closer $\Delta_{P}$ is to zero, the sharper the LDA$\rightarrow$HDA transformation. We also note that negative values for $\Delta_{P_{\text{IS}}}$ indicate a van der Waals-like loop in $P_{\text{IS}}(\rho)$. It follows from Fig. 7 that, as $\Delta_{P}$ decreases, $\Delta_{P_{\text{IS}}}$ also decreases and becomes negative. That is, the sharper the transformation, the more pronounced is the van der Walls-like loop in $P_{\text{IS}}(\rho)$. Consistent with Fig. 4, we find that during the decompression runs, $P_{\text{IS}}(\rho)$ is very smooth showing no van der Waals-like loop, at least for the rates considered [see Fig. 6(a)].

The behavior of $e_{\text{IS}}(\rho)$ is shown in Fig. 6(b). During compression at $T=280$ K, the system is able to reach the equilibrium liquid state and hence, it samples the same values of $e_{\text{IS}}$ accessed by equilibrium liquid (see supplementary material). However, at all lower temperatures, a concavity in $e_{\text{IS}}$ develops during the LDA$\rightarrow$HDA transformation, a feature again reminiscent of the behavior expected for equilibrium systems during phase transitions (cf. Sec. II). To clarify this point, we show in Fig. 7 the minimum curvature (i.e., maximum concavity) of $e_{\text{IS}}(\rho)$ during the LDA$\rightarrow$HDA transformation,

[TABLE]

as a function of $\Delta_{P}$. Interestingly, $e_{\text{IS}}(\rho)$ exhibits a mild concavity at all $T<280$ K, even at 20 and 80 K where $P_{\text{IS}}(\rho)$ shows no van der Walls-like loop. It follows that, a concavity in $e_{\text{IS}}(\rho)$ is not a sufficient (anomalous) property of the PEL for a glass to exhibit a first-order-like transition. Indeed, as argued in Ref. Sun et al., 2018, the van der Waals-like loop in $P_{\text{IS}}(\rho)$ and a concavity in $e_{\text{IS}}(\rho)$ seem to be necessary (but not sufficient) conditions for a glass to exhibit a first-order-like phase transition. In addition, we note that the concavity in $e_{\text{IS}}(\rho)$ is rather $T$-independent (at least for $q_{\text{c}}=30$ K/ns and $q_{\text{P}}=300$ MPa/ns).

During decompression, $e_{\text{IS}}(\rho)$ decreases monotonically with decreasing density until the system reaches its limit of stability at $\rho\approx 0.85-0.90$ g/cm3, where the amorphous ices are prompt to fracture [see Fig. 6(b)]. An exception is the case of $T=200$ K, were a very mild concavity seems to develop at $\approx 1.1$ g/cm3. Interestingly, the amorphous ices obtained at $\approx 0.9$ g/cm3 after decompression have a very large IS energy relative to the starting LDA configurations at $\approx 0.95$ g/cm3. This implies that the recovered LDA-like states are stressed glasses located in high regions of the PEL, within the corresponding LDA domain (cf. also Ref. Giovambattista et al., 2016).

$\mathcal{S}(\rho)$ behaves similarly to $P_{\text{IS}}(\rho)$ [see Figs. 6(a) and (c)]. Specifically, during compression at 160 and 200 K, $\mathcal{S}(\rho)$ shows a van der Waals-like loop. In other words, the sampled basins become narrower during compression up to $\approx 1.05$ g/cm3, then become wider during compression up to $\approx 1.2$ g/cm3 and then become narrower again. At 20 and 80 K (glass state), and at 280 K (liquid state), $\mathcal{S}(\rho)$ increases monotonically upon compression. During the decompression runs, $\mathcal{S}(\rho)$ decays monotonically with decreasing $\rho$ at all $T$ studied, i.e., the basins sampled by the system become wider as density decreases. Fig. 6(c) is fully consistent with previous studies on glassy ST2 water Giovambattista et al. (2016, 2017). We note, however, that in the case of a water-like monoatomic system that exhibits an LDA$\rightleftarrows$HDA transformation, $\mathcal{S}(\rho)$ shows no van der Waals-like loop during the LDA$\rightleftarrows$HDA transformation Sun et al. (2018).

We conclude this section with a brief discussion on the phenomenology found for the compressions and decompressions at 200 K. As noted above this temperature is higher than all $T_{\text{c}}$ estimates for TIP4P/2005 Abascal and Vega (2010); Sumi and Sekino (2013); Singh et al. (2016); Biddle et al. (2017); Handle and Sciortino (2018a). Hence, it may be surprising that we found relatively sharp density steps during compression and decompression (including a hysteresis), as well as a van der Waals-like loop in $P_{\text{IS}}$ and a concavity in $e_{\text{IS}}$. These are signatures that we also find for the LDA$\rightarrow$HDA transformation, and which are similar to what is expected for first-order phase transitions. Since $$200\text{,}\mathrm{K}$>T_{\text{c}}$, a phase transition is ruled out as the cause of this phenomenology. Instead, the anomalous properties of the PEL at $T=200$ K can be rationalized by noticing that the compression rate used is large enough so that the system is not able to reach equilibrium during compression and decompression. This reminds us that phenomena observed in non-equilibrium systems should be interpreted with caution. We note that, as we will show in Sec. IV.2.3, a decrease in the compression rate $q_{\text{P}}$ increases the sharpness of the LDA$\rightleftarrows$HDA transformation at 80 K (and leads to a more pronounced van der Waals-like loop in $P_{\text{IS}}$). Instead, at $T=$200\text{,}\mathrm{K}$>T_{\text{c}}$, reducing $q_{\text{P}}$ must bring the system to equilibrium, as we find for the case $T=280$ K. Accordingly, a slower rate $q_{\text{P}}$ must decrease the sharpness of the apparent of the LDA/LDL-HDA/HDL transformation at $T=200$ K.

IV.2.2 Cooling-Rate Dependence

In order to study the cooling rate effects on the LDA$\rightleftarrows$HDA transformation, we consider LDA samples prepared at 1 bar and 80 K using cooling rates $q_{\text{c}}=0.01$–100 K/ns. All samples are then compressed/decompressed at $T=80$ K using $q_{\text{P}}=10$ MPa/ns. Again, we point out that the smallest cooling rate studied here (0.01 K/ns $=10^{7}$ K/s) corresponds to the estimated rate reached in experimental hyperquenching techniques Brüggeller and Mayer (1980); Dubochet and McDowall (1981); Mayer and Brüggeller (1982); Mayer (1985); Kohl et al. (2000).

Fig. 8 reports $\rho(P)$ during the compressions/decompressions cycles starting from LDA forms prepared using different $q_{\text{c}}$. In close similarity to our discussion in Sec. IV.2.1, Fig. 8 shows that all LDA samples experience a sudden densification that signals the LDA$\rightarrow$HDA transformation during compression. Instead, the HDA$\rightarrow$LDA transformation is rather smooth with no evident LDA-like state recovered at negative pressure. HDA evolves continuously until it fractures at $P<-0.7$ GPa.

Fig. 8 clearly shows that, as $q_{\text{c}}$ decreases, the densification step associated to the LDA$\rightarrow$HDA transformation becomes sharper and shifts slightly to higher $P$. We stress that, the LDA$\rightarrow$HDA transformation at $q_{\text{c}}=0.1$–0.01 K/ns is remarkably similar to the experimental results Mishima (1994); Handle and Loerting (2018a). Indeed, based on the slope of the transformation, it is difficult to deny the first-order phase transition nature of the LDA$\rightarrow$HDA transformation. Instead, during decompression of HDA, $\rho(P)$ decreases monotonically, with no evident LDA-like state (at the present conditions). Interestingly, $\rho(P)$ behaves identically during decompression for all samples considered. This indicates that, once HDA forms, the system looses memory of the process followed to prepare LDA. Accordingly, the HDA$\rightarrow$LDA transformation is unique for $T=80$ K and $q_{\text{P}}=10$ MPa/ns.

We summarize these results in Fig. 9, where the pressure of the LDA$\rightarrow$HDA transformation is shown as function of $1/q_{\text{c}}$. For comparison, we also include available experimental data. It follows from Fig. 9 that as $1/q_{\text{c}}$ increases and $q_{\text{c}}$ approaches the experimental rate, the LDA$\rightarrow$HDA transformation pressure seems to reach an asymptotic value of $\approx 950$ MPa. Although this pressure is larger than the corresponding experimental pressure of $\approx 550$ MPa, one should note that the compression rates in experiments and MD simulation are very different (the role of $q_{\text{P}}$ on our systems is discussed in the next section). We note that, the pressure associated to the HDA$\rightarrow$LDA transformation is not shown in Fig. 9, because, at the present $T$ and $q_{\text{P}}$, this transformation is very smooth.

The PEL properties sampled by the system during the compression/decompression cycles in Fig. 8 are shown in Fig. 10. During compression, $P_{\text{IS}}(\rho)$ increases monotonically for LDA forms prepared with fast cooling rates, $q_{\text{c}}\geq 30$ K/ns [Fig. 10(a)]. However, as $q_{\text{c}}$ decreases, a clear van der Waals-like loop develops in $P_{\text{IS}}(\rho)$ during the LDA$\rightarrow$HDA transformation. In particular, this van der Waals-like loop becomes more pronounced as the slope of $\rho(P)$ in Fig. 8 becomes sharper. This is clearly shown in Fig. 11 where the $\Delta_{P_{\text{IS}}}$ is plotted as function of $\Delta_{P}$. Interestingly, during the decompression of HDA, $P_{\text{IS}}(\rho)$ decreases monotonically (i.e., it exhibits no van der Waals-like loop), which is also consistent with the smooth behavior of $\rho(P)$ shown in Fig. 8 along the decompression paths. We also note that the behavior of $P_{\text{IS}}(\rho)$ during decompression of HDA is independent of $q_{\text{c}}$. This, again, is consistent with the assessment that the system looses memory once it reaches the HDA state.

The behavior of $e_{\text{IS}}(\rho)$ [Fig. 10(b)], is fully consistent with the evolution of $P_{\text{IS}}(\rho)$ during the compression/decompression cycles. Specifically, during compression $e_{\text{IS}}(\rho)$ exhibits a concavity region during the LDA$\rightarrow$HDA transformation that becomes more pronounced as $q_{\text{c}}$ decreases. In other words, as the van der Waals-like loop in $P_{\text{IS}}(\rho)$ becomes more pronounced [see Fig. 10(a)], and the LDA$\rightarrow$HDA transformation becomes sharper (see Fig. 8), $e_{\text{IS}}(\rho)$ becomes increasingly a more concave function of $\rho$. This is also visible in Fig. 11 where the $\Delta_{e_{\text{IS}}}$ is plotted as function of $\Delta_{P}$. Instead, during decompression of HDA, $e_{\text{IS}}(\rho)$ shows no concavity which is consistent with the lack of a van der Waals-like loop in $P_{\text{IS}}(\rho)$ [see Fig. 10(a)] and the smooth behavior of $\rho(P)$ (see Fig. 8) along the corresponding path.

A subtle point follows from Fig. 10(b). Specifically, the $e_{\text{IS}}(\rho)$ corresponding to the starting LDA forms at $\rho\approx 0.95$–$0.97$ g/cm3 become more negative with decreasing $q_{\text{c}}$. In other words, as the cooling rate decreases during the preparation of LDA, the system is able to relax for longer times during the cooling process, and the final LDA state is able to reach deeper regions of the PEL. Accordingly, our results imply that in order to observe a first-order-like phase transition during the compression of LDA, one should start with LDA forms located deep within the LDA region of the PEL (cf. Ref. Giovambattista et al., 2017 for the case of glassy ST2 water). From a microscopic point of view, the LDA forms with lower $e_{\text{IS}}(\rho)$ are characterized by higher tetrahedral order as is obvious from the structure factor shown in Fig. 3 (cf. also Refs. Engstler and Giovambattista, 2017; Martelli et al., 2018). This high degree of tetrahedrally makes the hydrogen-bond network of water stronger and more resistant to collapse under pressure during the LDA$\rightarrow$HDA transformation. Thus, it is the high tetrahedral order in LDA what is required to observe a sharp first-order like transition.

Regarding the curvature of the IS sampled by the system [Fig. 10(c)], we find that the shape function $\mathcal{S}(\rho)$ follows closely the behavior of $P_{\text{IS}}(\rho)$ during the compression/decompression cycles. For example, when $P_{\text{IS}}(\rho)$ shows a van der Waals-like loop, $\mathcal{S}(\rho)$ does it as well. When $P_{\text{IS}}(\rho)$ is a monotonic function of $\rho$, $\mathcal{S}(\rho)$ is also a monotonic function. This parallel behavior of $P_{\text{IS}}(\rho)$ and $\mathcal{S}(\rho)$ was also noted in Sec. IV.2.1.

IV.2.3 Compression Rate Dependence

Next, we study the effects of varying $q_{\text{P}}$ on the LDA$\rightleftarrows$HDA transformation. All LDA samples considered in this section are prepared by isobaric cooling with $q_{\text{c}}=0.1$ K/ns and the compression/decompression runs are performed at $T=80$ K. The behavior of $\rho(P)$ during these runs ($q_{\text{p}}=0.1$–1000 MPa/ns) is shown in Fig. 12. It is visible that decreasing $q_{\text{P}}$ makes the LDA$\rightarrow$HDA transformation much sharper and it shifts it to lower $P$. This is a reasonable behavior since the system has more time to relax during the transformation as $q_{\text{P}}$ is reduced. Similarly, during the decompression process, the HDA$\rightarrow$LDA transformation also becomes slightly more evident and shifts to higher (less negative) $P$ with decreasing values of $q_{\text{P}}$.

The results from Fig. 12 are summarized in Fig. 13 where we also include available experimental data. It is again visible that the LDA$\rightarrow$HDA transformation pressure decreases and the density steps become steeper as $q_{\text{P}}$ decreases. In particular, the MD data extrapolates reasonably well to the experimental data from Ref. Mishima, 1994. Unfortunately there is no systematic experimental study of the $q_{\text{P}}$ dependence of the LDA$\rightarrow$HDA transformation at 80 K. At 125 K, however, experiments found no significant change in the transformation pressure as the rate was increased from $0.1$ MPa/s to $100$ MPa/s Loerting et al. (2006). During decompression we find no HDA$\rightarrow$LDA transition at positive pressures for all rates studied (please note that the smallest $q_{\text{P}}$ used for the decompressions is 1 MPa/ns), a finding consistent with experiments Mishima (1994). Interestingly, the MD simulation data in Fig. 12 shows that the hysteresis in $~{}\rho(P)$ during the LDA$\rightleftarrows$HDA transformation becomes smaller as $q_{\text{P}}$ decreases. It seems, however, that changes in $q_{\text{P}}$ affect the LDA$\rightarrow$HDA more than the HDA$\rightarrow$LDA transformation.

The behavior of $P_{\text{IS}}(\rho)$, $e_{\text{IS}}(\rho)$ and $\mathcal{S}(\rho)$ during the compression/decompression cycles in Fig. 12 are shown in Fig. 14. During compression, a van der Waals-like loop develops in $P_{\text{IS}}(\rho)$ [see Fig. 14(a)] as well as in $\mathcal{S}_{\text{IS}}(\rho)$ [see Fig.14(c)] for decreasing values of $q_{\text{p}}$. In particular, we note that the anomalous behavior in $P_{\text{IS}}(\rho)$ becomes more pronounced as the LDA$\rightarrow$HDA transformation becomes sharper. In addition, $e_{\text{IS}}(\rho)$ shows a concavity region during the LDA$\rightarrow$HDA transformation [see Fig.14(b)]. These findings are also evident in Fig. 15, where we plot $\Delta_{P_{\text{IS}}}$ and $\Delta_{e_{\text{IS}}}$ as a function of $\Delta_{P}$. Here we note a slight increase of $\Delta_{e_{\text{IS}}}$ at low $\Delta_{P}$ although it is unclear if this increase is indeed significant given the variance of the data.

During decompression, there is no anomalous behavior in any of the PEL properties studied. $P_{\text{IS}}(\rho)$ and $\mathcal{S}(\rho)$ decay monotonically during decompression and $e_{\text{IS}}(\rho)$ has positive curvature at all densities. This is consistent with the smooth decrease in $\rho(P)$ during the decompression of HDA (see Fig. 12). We note, however, that as the HDA$\rightarrow$LDA transformation becomes sharper with decreasing $q_{\text{P}}$, the slope of $P_{\text{IS}}(\rho)$ at the transformation seems to approach zero. One may expect that as $q_{\text{P}}$ decreases below 1 MPa/ns, $P_{\text{IS}}(\rho)$ may also exhibit a van der Waals-like loop during the HDA$\rightarrow$LDA transformation. A consistent trend follows from Fig. 14(b) for the case of $e_{\text{IS}}(\rho)$. During the decompression path, the curvature at $\rho\approx 1.1$ g/cm3 decreases with decreasing $q_{\text{P}}$ and, for $1$ MP/ns, the curvature of $e_{\text{IS}}(\rho)$ is practically zero. Hence we expect that, $e_{\text{IS}}(\rho)$ should develop a concavity during the HDA$\rightarrow$LDA transformation for $q_{\text{P}}<1$ MPa/ns. We note that the sharpness of the HDA$\rightarrow$LDA transformation during the decompression path is sensitive to $q_{\text{P}}$, but not $q_{\text{c}}$. Not surprisingly, only reducing $q_{\text{P}}$ allows the system to increasingly relax during the HDA$\rightarrow$LDA transformation. This is indicated by the IS energies explored by the system during decompression in Fig. 14(b). As $q_{\text{P}}$ decreases, the $e_{\text{IS}}(\rho)$ of the system at $\approx 0.95$ g/cm3 decreases meaning that the system accesses LDA configurations that are located deeper in the LDA state and closer to the starting LDA. Instead, as shown in Fig. 10(b), the system reaches the same $e_{\text{IS}}$ at $\rho\approx 0.9$ g/cm3 during the decompression of HDA for all cooling rates $q_{\text{c}}$ considered.

V Conclusions

In this report, we explored the PEL of TIP4P/2005 water during the pressure induced LDA$\rightarrow$HDA and HDA$\rightarrow$LDA transformations. The initial LDA form for the compression runs was produced by quenching the liquid at $P=0.1$ MPa with cooling rates as low as $0.01$ K/ns. This cooling rate coincides with rates reached in hyperquenching experiments Brüggeller and Mayer (1980); Dubochet and McDowall (1981); Mayer and Brüggeller (1982); Mayer (1985); Kohl et al. (2000). Reducing the cooing rate from 100 K/ns to 0.01 K/ns allows the system to access deeper and deeper regions of the PEL. From a microscopic point of view, the deeper the system is within the LDA region of the PEL, the more tetrahedral it is. At our slowest cooling rate, the density of LDA is slightly above the density of TIP4P/2005 ice I${}_{\text{h}}$.

During the compression-induced LDA$\rightarrow$HDA transformation, a pronounced density increase occurs. The sharpness of this density increase was found to be strongly dependent on the LDA preparation process (i.e., cooling rate $q_{\text{c}}$) as well as the compression rate $q_{\text{P}}$ and temperature $T$. At $T<80$ K, the LDA$\rightarrow$HDA transformation is rather smooth, due to the slow kinetics of the transformation and relatively fast compression rates employed. However, as $T\rightarrow T_{\text{c}}$, the LDA$\rightarrow$HDA transformation becomes sharp, reminiscent of a first-order phase transition, as observed in experiments. By studying the compression-induced LDA$\rightarrow$HDA transformation at fixed temperature and compression rate ($T=80$ K, $q_{\text{P}}$=10 MPa/ns), we find that reducing $q_{\text{c}}$ leads to a sharper transformation between LDA and HDA. In other words, as the tetrahedrality of LDA increases, the LDA$\rightarrow$HDA transformation becomes more reminiscent of a first-order phase transition. Similarly, by studying the compression-induced LDA$\rightarrow$HDA transformation at fixed temperature and cooling rate ($T=80$ K, $q_{\text{c}}$=0.1 K/ns), we find that reducing $q_{\text{P}}$ leads to a sharper transition between LDA and HDA. Remarkably, our LDA$\rightarrow$HDA transformation pressures, obtained at different compression rates, extrapolate fairly well to the experimental transformation pressure.

During the decompression of HDA, we can only observe an HDA$\rightarrow$LDA transformation at $T=160$ K, close to the estimated LLCP. At lower $T$, the transformation is rather smooth. However, we also find that as the decompression rate $q_{\text{P}}$ decreases, an HDA$\rightarrow$LDA transformation becomes more apparent in the behavior of $\rho(P)$. Interestingly, the HDA$\rightarrow$LDA transformation is insensitive to the cooling rate employed in the preparation of LDA. Consistent with the case of ST2 water Giovambattista et al. (2016, 2017), our results indicate that once HDA forms (at $P=2$ GPa), the system seems to completely loose memory of its history.

In agreement with previous studies Giovambattista et al. (2016, 2017); Sun et al. (2018), we find that at those conditions ($T$, $q_{\text{c}}$, $q_{\text{P}}$) where the LDA$\rightarrow$HDA transformation is reminiscent of a first-order-like phase transition, the PEL properties sampled by the system during the transformation become anomalous. Specifically, during the LDA$\rightarrow$HDA transformation, (i) $e_{\text{IS}}(\rho)$ becomes a concave function of $\rho$ and (ii) a van der Waals-like loop develops in $P_{\text{IS}}(\rho)$. In addition, and in agreement with results obtained for ST2 water Giovambattista et al. (2016, 2017), we also find that ${\cal S}(\rho)$ is anomalous, exhibiting also a van der Waals-like loop. These features are very weak or absent during smooth HDA$\rightarrow$LDA transformations.

Our studies at $T=80$ K and different rates ($q_{\text{c}}$, $q_{\text{P}}$), show that the anomalous van der Waals-like loop in $P_{\text{IS}}(\rho)$ becomes more pronounced as the LDA$\rightarrow$HDA transformation becomes sharper, i.e., more reminiscent of a first-order phase transition. The case of $e_{\text{IS}}(\rho)$ is less clear but our MD simulations show that $e_{\text{IS}}(\rho)$ remains a concave function of $\rho$ at all conditions studied. We argue that the PEL anomalies (i) and (ii) are necessary but not sufficient conditions for a system to exhibit a fist-order-like phase transition between LDA and HDA forms. Indeed, as discussed in Sec. II, for the case of a supercooled liquid (e.g., close to the glass transition temperature), these anomalies of the PEL may originate a first-order phase transition between two liquid states.

Previously similar PEL studies were conducted for the SPC/E and ST2 models of water Giovambattista et al. (2003, 2016, 2017). The main difference between SPC/E and ST2 water is that an LLCP is accessible in (metastable) equilibrium simulations of ST2 water Poole et al. (1992, 2005); Cuthbertson and Poole (2011); Liu et al. (2012); Smallenburg and Sciortino (2015); Palmer et al. (2014, 2018a, 2018b), while an LLCP is not accessible in SPC/E water Scala et al. (2000a, b); Sciortino et al. (2003). Consistent with this difference the PEL of SPC/E water shows smooth changes during the LDA$\rightleftarrows$HDA transformations including a very weak, concavity in $e_{\text{IS}}$  Giovambattista et al. (2003). In ST2 on the other hand, van der Waals-like loops in $P_{\text{IS}}(\rho)$ and $\mathcal{S}(\rho)$ as well as a concavity in $e_{\text{IS}}(\rho)$ are present Giovambattista et al. (2016, 2017). Even a maximum in $e_{\text{IS}}(\rho)$ was reported, consistent with the presence of two megabasins Giovambattista et al. (2016, 2017). That is, the signs expected for a first-order-like phase transition in the PEL are significantly weaker in SPC/E water than in ST2 water. This suggests that the glass phenomenology observed in SPC/E water can be thought of as “supercritical”, analogous to a liquid$\rightleftarrows$gas transformation at $T>T_{\text{c}}$. For ST2 water, the glass phenomenology resembles more a “subcritical” first-order phase transition. In our study of TIP4P/2005 water we find PEL features similar to ST2 water, including van der Waals-like loops in $P_{\text{IS}}(\rho)$ and $\mathcal{S}(\rho)$ as well as a concavity in $e_{\text{IS}}(\rho)$. However, we find no maximum in $e_{\text{IS}}(\rho)$ for TIP4P/2005 water. That is, in the case of TIP4P/2005 water, the two distinct regions of the PEL associated to LDA and HDA are not necessarily two different megabasins of the PEL. Similar conclusions were found in Ref. Sun et al., 2018 for a water-like model system. We stress that our study of TIP4P/2005 water is consistent with studies indicating an LLCP in this model Abascal and Vega (2010); Wikfeldt et al. (2011); Sumi and Sekino (2013); Yagasaki et al. (2014); Russo and Tanaka (2014); Singh et al. (2016); Biddle et al. (2017); Handle and Sciortino (2018a).

We conclude by noticing that the amorphous ices sample regions of the PEL that are different from the ones sampled by the equilibrium liquid (see supplementary material). This has profound implications in the relationship between the liquid and glass state. Specifically, this implies that the PEL regions sampled by the equilibrium liquid and the glass state differ. It follows that it may not be possible to predict quantitatively the behavior of the glass state based only on properties of the equilibrium liquid. In this work, we found that, indeed, there are no direct correlations between the LDL$\rightleftarrows$HDL spinodal lines of TIP4P/2005 liquid water obtained from the PEL-EOS of Ref. Handle and Sciortino, 2018a and the corresponding transformation pressures between LDA and HDA. However, in the LLCP scenario, the LDA$\rightleftarrows$HDA transformation lines are extensions of the LDL$\rightleftarrows$HDL spinodal lines into the glass domain.

Supplementary Material

In the supplementary material we compare the PEL regions sampled during the pressure induced LDA$\rightleftarrows$HDA transformations to the PEL regions sampled by the equilibrium liquid. We show that the regions of the PEL sampled by the equilibrium liquid and the amorphous ices are, in general, different.

Acknowledgements.

PHH thanks the Austrian Science Fund FWF (Erwin Schrödinger Fellowship J3811 N34) and the University of Innsbruck (NWF-Project 282396) for support. FS acknowledges support from MIUR-PRIN grant 2017Z55KCW. The computational results presented have been achieved using the HPC infrastructure LEO of the University of Innsbruck and the HPCC of CUNY. The CUNY HPCC is operated by the College of Staten Island and funded, in part, by grants from the City of New York, State of New York, CUNY Research Foundation, and National Science Foundation Grants CNS-0958379, CNS-0855217 and ACI 1126113.