Potential energy landscape of the apparent first-order phase transition between low-density and high-density amorphous ice

Water is one of the most complex liquids, exhibiting many anomalous thermodynamic and dynamical properties (see, e.g., Refs. 1 and 2). In the glassy state, water is also a complex substance.^1,3–5 Notably, amorphous solid water can be found in at least two different forms, low-density amorphous ice (LDA) and high-density amorphous ice (HDA), having very different properties. For example, the densities of LDA and HDA differ by $≈ 20%$⁠.^6–10 A large number of experiments indicate that LDA and HDA can be interconverted via many thermodynamic paths, such as isobaric heating and isothermal compression/decompression processes (see, e.g., Refs. 1, and 3–5). The LDA-HDA transformation, between properly annealed LDA and HDA forms, is rather sharp and reversible, and exhibits many of the characteristics of a first-order phase transition.^7,11–14 The explanation of this unusual behavior requires answers to fundamental questions of statistical mechanics, such as how to properly define or interpret a “phase-transition” between out-of-equilibrium glassy states.

In this work, we address these questions by studying the LDA-HDA transformations in water, and the relationship of these glasses with the liquid state, using the potential energy landscape (PEL) approach (see, e.g., Ref. 15). The PEL is a statistical mechanical approach that formally separates the configurational contributions to the partition function into contributions from local energy minima (“inherent structures (IS)”), and vibrational excitations within the basins of attraction surrounding these minima. Consequently, the PEL formalism has been used extensively to study supercooled liquids and glasses,^16–18 and other equilibrium systems,¹⁹ where such a separation between vibrational and configurational degrees of freedom arises naturally.

Specifically, for a system of N particles, the PEL is the hyper-surface in $(3N+1)$-dimensional space defined by the potential energy of the system as a function of the $3N$ coordinates, $U(r→1,r→2,….r→N)$⁠. At any given time $t$⁠, the system is represented by a single point in the PEL defined by the particle coordinates ${r→1(t),r→2(t),….r→N(t)}$⁠. As time evolves, the representative point of the system moves, sampling different regions of the PEL. In the PEL approach, the thermodynamic^17,20 and dynamic^21–24 properties of the system can be defined in terms of the topography of the PEL regions being sampled. The topography of the PEL can be rather complex with comparatively shallow basins residing within deeper and broader “megabasins”.^15,25,26 In the case of supercooled liquids, the free energy can be expressed in terms of three basic properties of the PEL:^15–17,20 (i) the average energy $EIS$ of the IS associated with the basins sampled in equilibrium; (ii) the number of IS having a given value of $EIS$⁠; and (iii) the average curvature of the basins associated with each inherent structure.

The behavior of glassy and liquid water is necessarily related, and several scenarios have been proposed to explain their unusual properties within a common framework (see, e.g., Refs. 27–30). One of the more widely accepted explanations is based on the idea that LDA and HDA are the glass counterparts of two liquids, low-density liquid (LDL) and high-density liquid (HDL), that are separated by a first-order phase transition.^3,30,31 This liquid-liquid phase transition (LLPT) ends at a liquid-liquid critical point (LLCP), recently estimated to exist at a temperature $TLLCP≈223$ K and pressure $PLLCP≈50$ MPa, based on experiments and theory.^8,32 The LLPT hypothesis was originally proposed on the basis of computer simulations using the ST2 water model,³⁰ leaving the validation of the hypothesis for experiments. Unfortunately, crystallization makes experimental verification challenging, and thus the hypothesis remains a point of debate. However, there is experimental evidence supporting the existence of a LLPT in water.^33–35 Indeed, most of the evidence supporting the LLPT hypothesis is from studies involving glassy water. A LLPT has been directly observed in experiments on other substances such as phosphorus³⁶ and cerium,³⁷ demonstrating the possibility of such a scenario. Furthermore, computer simulations of atomistic (see, e.g., Refs. 38–41), nanoparticle⁴² and molecular systems including modified water models^43–45 have shown the possibility for a LLPT in the equilibrium (as opposed to metastable) region of the phase diagram, establishing that a LLPT can exist as a proper phase transition in the thermodynamic limit.

In order to probe the apparent first-order transitions between LDA and HDA states, we employ the ST2 model. This choice is important since, contrary to other models such as the SPC/E and mW models,^46–48 the ST2 model reproduces the sharp LDA-HDA transformation observed in experiment. The behavior of ST2 water in glassy states has been recently characterized in detail and it has been shown that it is in qualitative agreement with experiments.^49,50 Specifically, in the case of the ST2 model, the density ρ has been observed to change abruptly with little change in the pressure P [i.e., $(∂P/∂ρ)T≈0$] during the LDA-HDA transformation; see, e.g., Refs. 46, 49, and 51.

In the present work, we explore the LDA-HDA transformation in ST2 water using the PEL formalism in order to clarify the thermodynamic differences between the LDA and HDA forms, and to assess the extent to which it may be appropriate to refer to this transformation as a first-order-like phase transition. As we show below, our results are consistent with the view that LDA and HDA occupy distinct megabasins of the PEL, and that the transformations between LDA and HDA exhibit a number of behaviors observed in well-defined equilibrium first-order phase transitions.

The structure of this manuscript is as follows: In Sec. II we describe the computer simulation details and methods employed. In Sec. III we discuss the changes in the PEL properties of ST2 water during the LDA-HDA transformations. We compare our results for ST2 water to those reported previously for SPC/E water in Section IV. We compare the LDA-HDA transformation and the liquid-to-ice-VII first-order phase transition using the PEL formalism in Sec. V. The regions of the PEL sampled by the liquids and glasses are compared in Sec. VI. A summary and discussion, including a simple model of the PEL for polyamorphic water, is presented in Sec. VII. In the  Appendix, we study the effects of reducing the compression/decompression rates on our results. Additional material is included as supplementary material where results at different temperatures are compared.

We perform extensive out-of-equilibrium molecular dynamics (MD) simulations of water using the ST2 model.⁵² Long-range (electrostatic) interactions are treated using the reaction field technique.⁵³ The details of the MD procedure and our implementation of the ST2 model are identical to that described in Ref. 54. Reference 49 contains a complete analysis of the thermodynamic and structural changes accompanying the pressure-induced LDA-HDA transformations to be analyzed in this work using the PEL formalism.

To summarize the thermodynamic procedure briefly, an LDA configuration is prepared by cooling equilibrium liquid water at $P=0.1$ MPa from $T0=350$ K down to the chosen temperature T for compression/decompression, using a cooling rate of $qc=30$ K/ns. This preparation method corresponds to the experimental process used to produce the LDA form known as hyperquenched glassy water (HGW), although we use a faster cooling rate than experiments; see discussions in Refs. 49–51. Our HGW configuration is then compressed isothermally, producing a sample of HDA. The resulting HDA form is then decompressed (at the same temperature T) leading to a new LDA sample. When subjected to sufficiently negative pressure, this LDA sample fractures. In order to compare our HGW sample and the LDA form we obtained by decompression of HDA, we also subject the HGW configuration to increasingly negative pressure, starting from $P=0.1$ MPa, until it also fractures.

Our compression/decompression rate is $qP= 300$ MPa/ns, which leads to sharp LDA-HDA transformations, as observed experimentally using much slower rates.^46,49,51 For each compression/decompression temperature, we perform 10 runs starting from independently generated starting configurations to account for sample-to-sample variations in the non-equilibrium state. In all cases, we use a cubic box with $N=1728$ water molecules with periodic boundary conditions.

During the compression/decompression runs, configurations are saved every 10 MPa and the corresponding IS are obtained using the conjugate gradient algorithm.⁵⁵ The virial expression for the pressure at the IS configuration defines the IS pressure, $PIS$ (see, e.g., Ref. 56). The curvature of a basin near the IS minimum is quantified by the shape function $SIS$⁠. $SIS$ is obtained from the set of eigenvalues ${ωi2}$ (where $i= 1$ to $6N$⁠) of the Hessian matrix evaluated at the IS configuration (see, e.g., Ref. 17 and references therein),

Here $ωi$ is the frequency of vibrational mode $i$⁠, and $ℏ=h/2π$ where h is Planck’s constant. The constant $A0=1$ kJ/mol is included so that the argument of the $ln$ function is dimensionless.

The same definition for $SIS$ was employed in Ref. 57 for the case of SPC/E water, allowing for a direct comparison of the present results with those reported in Ref. 57. While $PIS$ and $EIS$ are calculated for 10 independent runs, the shape function $SIS$ is calculated for only 2 of these runs, owing to the computational expense of evaluating the eigenvalues ${ωi2}$⁠. Technical details on the evaluation of the Hessian can be found in Refs. 58 and 59 where the density of states of water in equilibrium is reported.

The PEL properties $EIS$⁠, $PIS$⁠, and $SIS$ are fundamental quantities in the PEL formalism.¹⁷ For example, for a low-temperature liquid in equilibrium (or metastable equilibrium), knowledge of $EIS$ and $SIS$ as a function of V and T is sufficient to quantify key contributions to the free energy. In this case, the Helmholtz free energy $F$ can be written as

where $Sconf$ and $Svib$ are the configurational and vibrational entropies and $Evib$ is the average energy of the system due to vibrational motion in the PEL around the IS. In the harmonic approximation,

where $f(N,T)$ is a function only of temperature and number of molecules,⁶⁰ and

In the case of water, there can also be significant anharmonic contributions.⁶¹ An extension of the PEL formalism to quasi-equilibrium liquids shows that $EIS$⁠, $PIS$⁠, and $SIS$ are sufficient to write an analogous expression for the free energy of out-of-equilibrium states, from which the thermodynamic properties of the liquid can be predicted.^62–64 In this case, however, knowledge of the PEL properties sampled in equilibrium is necessary.¹⁷ 

We focus on the three fundamental properties of the PEL (⁠$PIS$⁠, $EIS$⁠, and $SIS$⁠) sampled by the system during the compression-induced LDA-to-HDA transformation and the decompression-induced HDA-to-LDA transformation. Doing so allows us to characterize the PEL basins associated with these states, and to determine the overall topography of the PEL for this system.

First we examine the LDA-to-HDA transformation at very low temperatures. Figs. 1(a)–1(f) show $PIS$⁠, $EIS$⁠, and $SIS$ as a function of ρ at two representative temperatures in the glass state, $T=20$ and $T=80$ K, during the LDA-to-HDA transformation. The behavior of $PIS(ρ)$⁠, $EIS(ρ)$⁠, and $SIS(ρ)$ does not differ significantly between these two T. Hence we focus our discussion on the case $T=80$ K.

In Fig. 2(b) of Ref. 49, the behavior of the pressure $P(ρ)$ during the LDA-to-HDA transformation at $T=80$ K is shown for the same 10 compression runs for which $PIS(ρ)$ is shown here in Fig. 1(b). The LDA-to-HDA transformation during compression at 80 K occurs at a pressure P_LDA-to-HDA which ranges between 1050 and 1200 MPa. Within the LDA regime (P < P_LDA-to-HDA) and the HDA regime (P > P_LDA-to-HDA), we note that P varies rapidly with ρ and is approximately linear, indicating that in these regimes the system responds to volume changes like a stiff, elastic solid. These elastic regimes correspond to approximately $ρ<0.95$ for LDA and $ρ>1.45$ g/cm³ for HDA. As shown in Fig. 1(b), we find that the behavior of $PIS(ρ)$ in these density ranges is similar to that observed for $P(ρ)$ in Ref. 49.

A significant difference between $P(ρ)$ and $PIS(ρ)$ is that P is a monotonic function of ρ, while $PIS$ is not. That is, during the LDA-to-HDA transformation (occurring in the range $0.95<ρ< 1.45$ g/cm³), $PIS(ρ)$ exhibits a van der Waals-like loop reminiscent of a first-order phase transition; see Fig. 1(b). This van der Waals-like loop in $PIS(ρ)$ becomes more evident at $T=160$ K [see Fig. 3(a)], i.e., as the liquid phase is approached from the glass state, and it vanishes at $T>TLLCP$ [see Fig. 3(b)]. Since our system is glassy and therefore not in equilibrium, a van der Waals loop in $PIS(ρ)$ does not definitively indicate a thermodynamic instability. At the same time, it is interesting that during the LDA-to-HDA transformation, we find $(∂PIS/∂ρ)T<0$⁠, in analogy with the condition of instability for an equilibrium system, $(∂P/∂ρ)T<0$⁠. We note that for $T→0$ K, these two inequalities become identical, since in this limit $PIS→P$ because the vibrational contributions to the pressure vanish.¹⁷ 

We next examine the behavior of $EIS$ during the LDA-to-HDA transformation at $T=80$ K, shown in Fig. 1(d). In the following analysis we assume that $EIS$ contributes to $F$ as described by Eq. (2). In this case, however, $Sconf$ is considered to be an out-of-equilibrium configurational entropy that depends on the glass preparation. This implies that $F$ is assumed to be an additive function of $EIS$ and $−kBTSIS$⁠. We also note that the inequality $(∂2F/∂V2)T<0$ identifies conditions at which a system is thermodynamically unstable as a single homogeneous phase. As shown in Fig. 2, the curvature of $EIS$ as a function of V is positive in the density ranges associated with the elastic regimes of both LDA and HDA, consistent with the mechanical stability of the system in these two forms. However, during the transformation of LDA to HDA (in the range $0.95<ρ< 1.45$ g/cm³), $EIS$ exhibits a pronounced negative curvature, i.e., $(∂2EIS/∂V2)T<0$⁠. Our data for $EIS$ therefore indicate that the influence of the PEL on the thermodynamic behavior of the system is to introduce a region of instability between the LDA and HDA regimes.

We also note that $EIS$ exhibits a weak maximum during the LDA-HDA transformation. Hence, the present results are consistent with the PEL of ST2 water containing two megabasins, one for LDA and another for HDA, separated by potential energy barriers. We stress that such barriers separating the LDA and HDA megabasins are not necessary for the system to become unstable. In other words, as argued in the previous paragraph (based on Eq. (2)), the negative curvature of $EIS$ (which does not need to be accompanied by potential energy barriers) may be sufficient to introduce an instability $[(∂2F/∂V2)T<0]$ between the LDA and HDA states.

The behavior of $SIS$ during the LDA-to-HDA transformation at $T=80$ K is shown in Fig. 1(f). During the elastic compression of LDA and HDA (for $ρ<0.95$ and $ρ>1.45$ g/cm³), $SIS$ increases monotonically. In other words, in regimes of elastic deformation, the system explores basins in the PEL which become narrower as density increases. However, the overall behavior of $SIS$ as a function of V across the LDA-to-HDA transformation is non-monotonic. At the beginning of the LDA-to-HDA transformation, $SIS$ decreases abruptly, indicating that “wider” basins (relative to LDA) become available to the system in the transformation zone. We also note that $(∂2SIS/∂V2)T>0$ throughout much of the LDA-to-HDA transformation region. Within the harmonic approximation of the PEL formalism (see Eqs. (2) and (3), the positive curvature of $SIS$ will act to reduce the instability introduced by the negative curvature of $EIS$ noted above. However, in terms of their respective contributions to F according to Eq. (2), the relative variation of $EIS$ over the LDA-to-HDA transformation range is approximately 20 times larger than for $kBTSIS$ when $T=80$ K. Hence the curvature of $EIS$ dominates over that contributed by $SIS$ in determining the overall curvature of F.

To summarize, Figs. 1(a)–1(f) show three characteristic features of the LDA-to-HDA transformation in the PEL: (i) a van der Waals-like loop in $PIS$⁠; (ii) negative curvature in $EIS$⁠; and (iii) non-monotonic variation of $SIS$⁠. These features are all consistent with behavior analogous to a first-order phase transition in an equilibrium system, as we show in Sec. V. We note that these characteristic features, found here at $T≤80$ K, become more pronounced as T increases within the range of T in which the system remains glassy. For example, see the results for $T= 160$ K in the left column of Fig. 3 for $ρ<1.6$ g/cm³. See also the supplementary material where we include $EIS(ρ)$⁠, $PIS(ρ)$⁠, and $SIS(ρ)$ for all temperatures studied.

During decompression, the behavior of $PIS(ρ)$⁠, $EIS(ρ)$⁠, and $SIS(ρ)$ are consistent with the system moving from the HDA megabasin back to the LDA megabasin. For example, at $T= 160$ K [Figs. 3(a), 3(c), and 3(e)], all three of the characteristic features noted above [points (i), (ii), (iii)] are observed during the HDA-to-LDA transformation. However, these features are less evident, or absent, during the HDA-to-LDA transformations at $T≤ 80$ K. As shown in Fig. 1, the negative curvature of $EIS$ is less prominent during decompression, and non-monotonic behavior in $PIS$ and $SIS$ is barely observed. This demonstrates that the IS basins of the PEL visited by the system during the LDA-to-HDA and HDA-to-LDA transformations are quite different, that is, the system follows different trajectories on the PEL.

A natural question follows: Is the LDA-to-HDA transformation reversible? As shown in Fig. 1, the values of $PIS$⁠, $EIS$⁠, and $SIS$ at $ρ≈0.87$ g/cm³ for (a) the LDA (HGW) used to initiate the compression runs (black lines) and (b) the LDA recovered from HDA via decompression (red lines) are different. While these two LDA forms are not identical, as shown in Ref. 49, both are structurally similar based on the corresponding OO, OH, and HH radial distribution functions. Reference 49 argued that both of these LDA forms should be considered to be different members of the LDA “family.” The present results based on the PEL approach support this view. A region of negative curvature in $EIS$ (accompanied by a weak maximum) is encountered both during compression and decompression, consistent with the view that these trajectories take the system first from the LDA to the HDA megabasin, and then back again to the LDA megabasin. However, $EIS$ for the recovered LDA is much higher than for the starting LDA (HGW); see Figs. 1(c), 1(d), and 3(c) for $ρ=0.87$ g/cm³. This difference suggests that the recovered LDA form is a highly stressed glass residing higher in the LDA megabasin than HGW. This would not be surprising given the relatively fast compression rates that are accessible in simulations.^{46,47,49,51,65} To confirm this interpretation, we subjected both the HGW and the recovered LDA forms to large negative pressures, to test if their properties become more similar when they are both brought close to their limits of mechanical tensile stress. As shown by the green lines in Figs. 1 and 3, the PEL properties of HGW do approach more closely those of the recovered LDA at negative pressure. Also consistent with this view, we show in the  Appendix that $EIS(ρ)$ for HGW and recovered LDA become closer to one another when the compression/decompression rate is reduced.

In a previous work, a similar PEL study of the LDA-HDA transformation in water was carried out using the SPC/E model.⁵⁷ The computer simulation protocol followed in Ref. 57 is similar to the procedure used in the present work. In particular, both studies employ the same compression/decompression and cooling rates, which is crucial for a proper comparison. The main difference between the SPC/E and ST2 models is that a LLPT is accessible in (metastable) equilibrium simulations of ST2 water^{30,45,66–68} while it is not accessible in SPC/E water via ordinary molecular dynamics simulations. Earlier studies^69,70 suggested that a LLPT may be present in SPC/E water at very low T, but it is not accessible at the available computational time scales. In short, Ref. 57 describes the PEL properties explored during the LDA-HDA transformation if the LLPT is not accessible to simulations. The present study corresponds to the case when the LLPT is accessible.

The main difference between the LDA-HDA transformation results for SPC/E and ST2 water is that in the case of SPC/E water, the evidence (at, e.g., $T≈80$ K; see Fig. 2 in Ref. 57) for first-order phase transition-like behavior in the PEL is significantly weaker. Specifically, during the LDA-HDA transformations in SPC/E water, (i) $PIS$ does not exhibit van der Waals-like loops; (ii) $EIS$ exhibits only very weak, barely noticeable, negative curvature upon compression, and not at all during decompression; and (iii) non-monotonic behavior in $SIS$ is also weak during compression, and absent during decompression. This pattern of behavior is consistent with the absence of an accessible LLPT in SPC/E water.

We also note that at the compression rates studied, the LDA-HDA transformation in SPC/E water exhibits a more gradual change of density on compression/decompression, compared to ST2 water or experimental water.^46,49,51 This suggests that the glass phenomenology observed in SPC/E water can be thought of as a “supercritical” LDA-HDA transformation, analogous to a liquid-gas transformation at $T>Tc$⁠. In contrast, for ST2 water, and perhaps for real water, the glass phenomenology corresponds to a “subcritical” first-order phase transition between LDA and HDA.

A theory for amorphous ices has been presented in Ref. 47 based on theory/computer simulations using the mW model, a coarse-grain model for water.⁷¹ In the supplementary material of Ref. 47 it is shown that for a compression rate $qP= 500$ MPa/ns, similar to the compression rate employed here and in Ref. 57, the slope of $ρ(P)$ in the LDA state and during the LDA-HDA transformation is similar, i.e., there is not a sudden change in density associated with the LDA-to-HDA transformation. The authors conclude that there is a state at which LDA, HDA, and the liquid coexist at an out-of-equilibrium triple point. Since the mW model seems to lack the sharp change in density associated with the LDA-to-HDA transformation, its behavior is closer to the case of SPC/E than ST2 water. It is thus suggestive that the theory proposed in Ref. 47 may apply only for the case of  “supercritical” LDA-HDA transformations, such as the transformation observed in SPC/E water.

We finally note some of the characteristic features indicative of a LDA-HDA “first-order phase transition,” found in ST2 water at $T≤160$ K, vanish with increasing temperature as $T→TLLCP$⁠, where $TLLCP≈245$ K is the LLCP temperature of ST2 water.⁷² For example, at temperatures slightly above $TLLCP$⁠, we find that the curvature of $EIS(ρ)$ increases (becomes less negative) and $PIS(ρ)$ shows no van der Waals-like loop (see supplementary material). In particular, the behavior of $EIS(ρ)$ and $PIS(ρ)$ becomes closer to that reported in Fig. 2 of Ref. 57 for the case of SPC/E water. A similar phenomenology is found at high temperature, in the liquid state. For example, the right column of Fig. 3 shows $EIS(ρ)$⁠, $PIS(ρ)$⁠, and $SIS(ρ)$ for the case of $T=280$ K $>TLLCP$⁠. At these temperatures, $EIS(ρ)$ and $PIS(ρ)$ increase smoothly and monotonically with increasing density (at $ρ<1.45$ g/cm³, where crystallization does not interfere).

In this section, we study the PEL behavior during the pressure-induced HDA-to-ice (at low temperature) and liquid-to-ice transformation (at high temperature). References 46, 49, and 67 showed that ST2 water crystallizes spontaneously at high-pressure during isothermal compression, or isothermal heating of HDA at high pressure. The resulting crystal has the structure of ice VII, characterized by interpenetrating tetrahedral networks. The occurrence of spontaneous crystallization in our system provides the opportunity to identify the PEL behavior that is characteristic of an unambiguous first-order phase transition. We can thus compare the behavior of the PEL properties sampled during the liquid-to-ice first-order phase transition and with that occurring during the LDA-HDA transformation.

We first focus on the case of isothermal compression at $T=160$ K; see the left column of Fig. 3. At this temperature, we observe the successive LDA-to-HDA and HDA-to-ice transformations upon compression. The similarities in the behavior of $EIS(ρ)$⁠, $PIS(ρ)$⁠, and $SIS(ρ)$ during these two transformations are striking. In both the LDA-HDA and HDA-ice transformations, (i) $PIS(ρ)$ exhibits a van der Waals loop; (ii) $EIS(ρ)$ exhibits negative curvature; and (iii) a sharp change in $SIS(ρ)$ occurs as the transformation proceeds. In addition, the changes observed in the PEL properties during both transformations are consistent with the system evolving from one megabasin of the PEL to another. Specifically, at $T= 160$ K, a maximum in $EIS$ separates LDA, HDA, and ice indicating that at this temperature, there are well-separated LDA, HDA, and ice megabasins in the PEL. The only minor difference is that during the HDA-ice transformation, $SIS(ρ)$ increases monotonically (see Fig. 3(e) for $ρ> 1.5$ g/cm³), while for the LDA-HDA transformation $SIS(ρ)$ increases non-monotonically (see Fig. 3(e) for $0.85<ρ<1.7$ g/cm³). Yet, in both cases, the basins become narrower [i.e., they have a larger $SIS(ρ)$] in the higher-density “phase.”

The right column of Fig. 3 shows the PEL properties at $T=280$ K. At this temperature the system starts in the liquid state during compression. For the compression rate employed, the system at this temperature is essentially in equilibrium, as demonstrated by the fact that the PEL properties overlap during compression and decompression. The main effect of increasing the temperature from $T= 160$ K $<TLLCP$ to $T=280$ K $>TLLCP$ is to lose the signature of a distinct megabasin associated with the HDL-like liquid. However, all the PEL features of the HDA-ice transition observed in the left panels of Fig. 3 also occur in the liquid-ice transition shown in the right panels.

We note that the use of the term “van der Waals loop” to describe the shape of the $PIS$ isotherms for the HDA-to-ice transformations depicted in Figs. 3(a) and 3(b) deserves some qualification. During a glass-to-crystal or liquid-to-crystal transformation we should not necessarily expect a smooth maximum followed by a smooth minimum in $PIS$ as a function of $ρ$⁠, as is observed for the glass-to-glass transformation. Under the assumption that $PIS$ dominates the density dependence of $P$ at fixed $T$⁠, these smooth extrema correspond approximately to metastability limits (analogous to mean-field spinodals) at which the isothermal compressibility diverges. Such divergences are not surprising for a liquid-liquid transition that terminates in a critical point, but there is no thermodynamic requirement that such a spinodal should occur during a liquid-to-crystal transition. Indeed, the liquid-to-crystal transition depicted in Fig. 3(b) exhibits quite sharp cusp-like changes at the extrema of $PIS$⁠. However, since the liquid-liquid and liquid-to-crystal transitions considered here are both first-order phase transitions, both can display metastability, in which $PIS$ first overshoots and then undershoots the coexistence value of $PIS$ during the transformation. It is in this broader sense that we refer to the behavior of $PIS$ for both the LDA-HDA and HDA-ice transformations as being “van der Waals-like” and a “signature” of first-order phase transitions.

In summary, these results suggest that during a true first-order phase transition, $PIS(ρ)$ shows a van der Waals-like loop, and $EIS(ρ)$ exhibits a region of negative curvature. Interestingly, in the case of the liquid-to-ice transformation, the system also crosses a potential energy barrier in the PEL that separates two distinct (liquid and ice) megabasins. The similarities in the PEL properties sampled during the liquid-to-ice and LDA-HDA transformations provide support for the view that the LDA-HDA transformation is an out-of-equilibrium first-order phase transition, the origins of the LDA-HDA and LDL-HDL transformations being the system moving between the same two (LDL/LDA and HDL/HDA) megabasins of the PEL.

We now consider how the LDA and HDA configurations compare with the configurations sampled by LDL and HDL in equilibrium. This question has deep implications for our understanding of glass and liquid polymorphism. For example, if the PEL basins sampled by the LDA (HDA) glass are the same as the ones sampled by the equilibrium LDL (HDL) liquid at some fictive temperature $Tf$⁠, then it becomes possible to provide a thermodynamic modeling of the LDA-HDA transformation in terms of the LDL-HDL first-order phase transition. We addressed this question in Ref. 57 for the case of SPC/E water and found that, surprisingly, during the LDA-HDA transformation the system explores regions of the PEL never sampled by the equilibrium liquid. In this section, we show that this conclusion also holds for the case of ST2 water.

To compare the regions of the PEL associated with LDA, HDA, and the equilibrium liquid, in Figs. 4 and 5 we present these states in the $PIS$-$EIS$⁠, and $SIS$-$EIS$ planes at $T=80$ and $160$ K. See Fig. 3 of Ref. 57 for a comparison with the case of SPC/E water. The starting LDA form (i.e., HGW) is close to the low-energy end of the $ρ=0.9$ g/cm³ isochores for the liquid in the $PIS$-$EIS$ and $SIS$-$EIS$ planes. This is expected since the density of LDA in this temperature range is approximately $0.87$ g/cm³. However, as soon as the compression starts, the glass deviates abruptly from the liquid isochores, exploring regions of the PEL that are not accessed by the liquid, similar to the behavior in SPC/E water.⁵⁷ For example, at $T=80$ K the LDA form at $ρ=1.0$ g/cm³ is characterized by $PIS=1000$ MPa, $EIS=−54.5$ kJ/mol, and $SIS=8.7$ (see empty diamond on the black lines in Fig. 4), while the corresponding values for the equilibrium liquid at the same density and $EIS$ are $PIS=100$ MPa and $SIS<8.35$ (see blue solid diamonds). Interestingly, the LDA-HDA transformation occurs in a region of the PEL that is far from the accessible configurations explored by the liquid at all densities.

During decompression of HDA (red lines in Figs. 4 and 5), the system returns to regions of the PEL with properties similar to those of the liquid. However, at a given density, the LDA form exhibits values of $PIS$⁠, $EIS$⁠, and $SIS$ that correspond to the equilibrium liquid at a different density. For example, Figs. 4(a) and 4(b) show that the decompressed glass at $ρ=1.0$ g/cm³ has $PIS≈−600$ MPa, $EIS≈−53$ kJ/mol, and $SIS≈8.05$ (see empty diamond on the red lines), while the corresponding values for the equilibrium liquid at the same density and $EIS$ are approximately $PIS=0$ MPa and $SIS=8.2$ (blue solid diamonds).

The results shown in Fig. 5(a) for $T=160$ K are similar. The main difference between the cases $T= 80$ and $160$ K occur during the decompression paths. Specifically, Figs. 4(a) and 4(b) show decompression paths (red lines) where both $PIS$ and $SIS$ decrease monotonically with decreasing $EIS$⁠. Instead, Figs. 5(a) and 5(c) show a sharp kink at low $PIS$ and $SIS$ corresponding to the HDA-LDA transformation. Yet, at both temperatures the LDA form at a given density is in a different region of the PEL than the equilibrium liquid at the same density.

For comparison, we include in Figs. 5(a)–5(d) the location in the $PIS$-$EIS$ and $SIS$-$EIS$ planes of the high-pressure ice that forms during compression at $T= 160$ K. Clearly, the PEL region explored by the system in the ice state is very different from the regions associated with LDA, HDA, or the equilibrium liquid.

In summary, we find that the PEL properties in our simulations of amorphous solid water during compression and decompression support the view that LDA and HDA correspond to two distinct megabasins in the PEL. We also show that the PEL behavior we observe in the LDA-HDA transformation is qualitatively the same as the PEL signatures observed during unambiguous first-order phase transitions in the same system, specifically the liquid-ice transitions. This finding is consistent with the interpretation of the LDA-HDA transformation as a sub-glass-transition manifestation of the equilibrium LLPT that has been demonstrated to occur in ST2 water. Comparison of our ST2 results with those obtained using SPC/E water demonstrates that when the LLPT is inaccessible, the phase-transition-like characteristics of the LDA-HDA transformation are correspondingly weakened or lost entirely. Given the closer similarity between the behavior of amorphous solid ST2 water and real amorphous solid water, our results therefore support the possibility that a LLPT occurs in deeply supercooled liquid water.

An interesting result of this work is that, in ST2 water, the regions of the PEL sampled by the liquid (LDL and HDL) and the glass (LDA and HDA) differ. This is in agreement with previous simulations of SPC/E water. It follows that, at least for the compression/decompression rates explored here, the concept of fictive temperature cannot be used to associate LDA or HDA with “frozen” equilibrium configurations of LDL or HDL. As a consequence, it seems unlikely that a simple thermodynamic modeling of the LDA-HDA transformation in terms of the LDL-HDL first-order phase transition is possible.

However, we also note that the fact that LDL/HDL and LDA/HDA sample different regions of the PEL is not at odds with the view that there are two distinct megabasins in the PEL, one associated with both LDL and LDA, and another megabasin corresponding to both HDL and HDA. A schematic and highly idealized PEL for real water, based on our results for ST2 water, is shown in Fig. 6. Only the LDL/LDA and HDL/HDA megabasins are illustrated; for clarity we have omitted the individual IS basins that are distributed all over the PEL. In Fig. 6, the region of the PEL explored by a given realization of a glass (LDA or HDA) should be small compared to the region accessible to the liquid. The slower a liquid is cooled toward the glass state, the deeper the system gets in the PEL in the final glassy state. Accordingly, the well-annealed LDA and HDA samples prepared in experiments are expected to lie quite close to the minima of their respective megabasins. The region of the PEL explored by the equilibrium liquid will include many higher energy configurations. The path followed in the PEL when compressing at experimental rates from LDA to HDA might look something like the green path in Fig. 6. However, if a HDA sample is prepared by compression of LDA at a very fast rate compared to experiments, as is done in our simulations, the glass can be driven into regions of the PEL quite different, and at higher energy, than those explored by either the liquid or the well-relaxed HDA obtained in experiments; a possible realization of this scenario is shown by the blue path in Fig. 6. From this point of view, we should not be surprised that different regions of the PEL are sampled by the glass during compression, relative to the liquid. A similar argument has been used to explain the differences between HDA samples prepared under different procedures.⁷³ 

One of the limitations of simulation studies of glass polymorphism is the fast compression/decompression rates employed, relative to experiments. In particular, we have shown that the starting LDA form for our compression/decompression cycle (HGW) is not identical to the LDA form recovered after decompression of HDA. However, as shown in the  Appendix, the difference between these two LDA forms decreases upon reducing the compression/decompression rate. Hence the present results are consistent with HGW and the recovered LDA form belonging to the same LDL/LDA megabasin. Our work supports the interpretation that the LDL/LDA and HDL/HDA megabasins each host a family of glasses and liquids, all with similar properties. In a following work, we will expand on the reversibility of the pressure-induced LDA-HDA transformation in simulations using the PEL approach.

See supplementary material for a close comparison of the $EIS$⁠, $PIS$⁠, and $SIS$ isotherms at all temperatures studied (⁠$T=20, 80, 160, 220, 250$⁠, and 280 K).

This project was supported, in part, by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grants Nos. CNS-0855217, CNS-0958379, and ACI-1126113. P.H.P. thanks NSERC and ACEnet. We thank Wesleyan University for computational resources.

In this section we study the effects of reducing the compression/decompression rate on the behavior of $EIS(ρ)$⁠, $PIS(ρ)$⁠, and $SIS(ρ)$ during the LDA-HDA transformations. We find that the qualitative PEL behavior is not altered by reducing $qP$ from 300 to 30 MPa/ns. Interestingly, while the compression-induced LDA-to-HDA transformation is clear at both rates studied, the decompression-induced HDA-to-LDA transformation becomes much more evident when $qP$ is reduced. In addition, the recovered LDA becomes closer to the starting LDA (HGW) as $qP$ decreases.

Fig. 7 shows $EIS(ρ)$⁠, $PIS(ρ)$⁠, and $SIS(ρ)$ during the pressure-induced LDA-to-HDA transformation at $T= 80$ K and $T= 160$ K, for both $qP=30$ and $300$ MPa/ns. During compression, reducing $qP$ has the main effect of reducing the values of $EIS(ρ)$⁠, $PIS(ρ)$⁠, and $SIS(ρ)$ during the LDA-HDA transformation. Yet at both rates one observes the signatures expected during the transition from one megabasin (LDA) to another (HDA) upon compression: (i) $PIS(ρ)$ shows a van der Waals-like loop; (ii) $EIS(ρ)$ exhibits a region of negative curvature; and (iii) $SIS(ρ)$ decreases sharply.

We observe a larger effect during the decompression path when $qP$ is reduced. As discussed in Sec. III, the HDA-LDA transformation is rather smooth at $T= 80$ K. The van der Waals-like loop in $PIS(ρ)$ and the negative curvature in $EIS(ρ)$ are weaker effects during decompression, relative to the compression process. However, upon reducing the decompression rate to $qP=30$ MPa/ns, we observe that these features become significantly more prominent. This observation is consistent with our conclusion that during decompression, the system transitions from the HDA megabasin back to the starting LDA megabasin of the compression cycle. In particular, Figs. 7(c) and 7(d) show that the value of $EIS(ρ)$ for the recovered LDA at $ρ=0.87$ MPa decreases with decreasing $qP$ and becomes closer to the IS energy of the starting LDA (HGW). That is, the recovered LDA is deeper in the LDA megabasin with decreasing compression/decompression rate, as one would expect. This supports the view that the LDA-HDA transformation in ST2 water is indeed reversible.⁴⁹ 

Fig. 8 shows $PIS(EIS)$ and $SIS(EIS)$ at $T=80$ and $160$ K for $qP=30$ and $300$ MPa/ns. At both rates studied, these properties exhibit the same qualitative behavior during compression (black and green lines) and decompression (red and blue lines). The main point of Fig. 8 is that the system samples IS that are never visited by the equilibrium liquid, even if the compression/decompression rate is reduced to $qP=30$ MPa/ns.