Insight into liquid polymorphism from the complex phase behaviour of a simple model

Polyamorphism, the occurrence of a material in more than one non-crystalline form, such as amorphous solids or distinctly different liquid phases, keeps fascinating materials scientists. While such phases do not exhibit long range order, their local atomic level structure can be markedly different, resulting in different physical characteristics, such as the density or electronic properties. The existence of multiple disordered phases, especially that of liquid phases with different densities, and the first order phase transitions between them have been long suspected to have a far reaching influence on materials properties, for example potentially explaining the high pressure melting line maximum observed in certain metals.[1]

Due to the substantial technical challenges, experimental observation of the peculiar liquid-liquid polyamorphism is rare in one-component systems. Over the past two decades its existence has been suggested only for a few such materials: phosphorus[2] silicon,[3, 4, 5], nitrogen[6], cerium[7] or triphenyl phosphite,[8] as well as the most widely studied material long suspected to have a liquid-liquid phase transition, water.[9, 10, 11]

Besides a surprisingly complex phase diagram (with so far 19 identified solid phases [12, 13]), water is also known to display a wealth of anomalous properties compared to other simple liquids. These include the non-monotonous temperature dependence of the density and the viscosity, as well as the increasing isobaric specific heat and isothermal compressibility when the temperature is decreased in the deep supercooled regime. This latter phenomenon can be interpreted as indicative of a potential phase transition in the metastable liquid state upon further cooling.[14, 15] The unconventional behavior of water is understood to stem from the fact that hydrogen bonding, which has a heavily directional character, is a major contributor to the cohesive interaction between molecules. This has two important consequences for the structure of water: bonding distances are very short and the first neighbour shell comprises, on average, only 4 molecules, in contrast to 12 in case of most other simple atomic or molecular liquids. At ordinary conditions, this dominance of hydrogen bonds gives rise to a strongly bound open tetrahedral structure of relatively low density, but it was suggested that larger pressures lead to the partial collapse of the hydrogen bonding network, resulting in higher coordination numbers.[16, 17] As a consequence, the equilibrium between these two competing local arrangements can bring about the existence of the two distinct liquid phases: the low- and high-density liquid (LDL and HDL),[18, 19] explaining numerous aspects of water’s anomalous behaviour,[20] and consistent with further experimental observations of two amorphous arrested glassy states showing structural similarities with them.[21, 22].

However, the suspected liquid-liquid critical point (LLCP), specifically in water, lies so deep in the metastable, supercooled region of the phase diagram, that it is highly challenging to access not just experimentally, but also in computer simulations[16]. While several studies indicate the existence of a HDL phase of a range of water potential models,[23, 24, 25] so far unequivocal proofs have only been presented in the case of the ST2[26][27] and TIP5P models[28].

In order to provide more insight into the properties of such short-lived liquid phases, the attention turned towards designing and studying model potentials, where the LLCP lies above the crystallization boundary, allowing the detailed examination of structural properties and thermodynamic response functions in the stable region of the phase diagram.

It was found that the existence of multiple liquid phases is a generic feature of model tetrahedral network fluids (patchy particles) with an LLCP[29] which can be brought to the stable regime upon increasing the angular flexibility of the model, which lowers the energy cost of the disruption of the network structure.[30]

Jagla showed in their seminal works that directionality is not a necessary feature for an interatomic potential model to display liquid-liquid polyamorphism. They studied the soft-core double step potential[31, 32] and the hard-core two-scale ramp model[33, 34], which are both monoatomic, isotropic pair potentials. The latter model, which was originally proposed, and its second critical point examined in the case of the one-dimensional fluid, by Hemmer and Stell [33], became known as the Jagla model. The Jagla or Hemmer-Stell model has since received increased attention, as it was found that its LLCP appears in the stable liquid regime, and thus has been used as a toy model to understand the polyamorphism of more complex liquids, such as water.[35, 36] A range of other non-directional soft-core potential models has been since studied,[37] several of which display anomalous behaviour and suspected multiple liquid or amorphous phases. Several of these are known to form not only close packed, but low-coordinated crystals as well,[38, 39, 40, 41, 42] moreover, it is possible to design such models using inverse statistical mechanical methods and targeting specific ground state structures.[43] However, the thermodynamic stability of these phases, especially at finite temperature, is not always well understood nor established unequivocally. In the current work we demonstrate the importance of exhaustive phase space search in clarifying phase behaviour on the example of the Jagla model.

The Jagla model is the combination of a hard-sphere core and two linear functions accounting for the repulsive and the attractive ramps, as depicted in the inset of Fig. 1. We note that the locations of phase transitions show a strong dependence on the parameterization and the implementation particulars, although the qualitative features of the phase diagram remain unchanged[44, 45, 46, 47, 48, 49, 50, 51, 52]. For details of the Jagla potential used in this work and for the derived Jagla units, we refer the reader to the inset of Fig. 1 and to the Supplementary Material[53].

In the original work of Jagla[34], both liquid phases were presumed to be thermodynamically stable, because no sign of crystallization was seen in the simulations, and no further effort had been made to identify possible crystalline structures. The heuristic criteria of Hansen-Verlet[54] and Lowen[55] were also applied, using the structure factor to conclude that the LLCP is in the stable region.

In studies where the structure of the solid phase is identified, the hexagonal-close-packed (hcp) crystal,[48, 50] or at certain potential parameters and at very low pressures, face-centred-cubic (fcc)[50, 36] have been reported. More exhaustive enumerations of further potential structures of the Jagla model have only been discussed in the two-dimensional system.[56, 57]

Lomba et al performed detailed calculations on the location of the melting line and the two triple points between fcc solid-LDL-vapour and fcc solid-LDL-HDL.[36] However, most other studies provide only an estimate to the temperature at which the solid-liquid first-order phase transition occurs, either by giving an upper limit to the melting line by gradually heating a system initially consisting of a crystalline configuration, or estimating conditions where both the liquid and solid phases remain stable within the same simulation cell.[48]

In order to fill this gap in our knowledge of the phase behaviour of the Jagla model, we performed an exhaustive study of its potential energy landscape, employing the nested sampling (NS) method,[58, 59, 60, 61] which allows unbiased exploration without the need of prior knowledge of the stable structures. NS simulations were carried out in the pressure range between $p=0.003-0.25~{}U_{0}/r_{0}^{3}$, allowing us to determine thermodynamical phase stability at arbitrary temperatures. We also performed thermodynamic integration and grand canonical ensemble simulations to accurately determine phase transition conditions[53]. Our results are summarised in Fig. 1 showing the revised pressure-temperature phase diagram of the Jagla model.

At pressures below $p=0.05$, our results show excellent agreement with previous findings,[50] confirming that the liquid freezes into a hcp solid, which has a lower density than that of the liquid, thus the slope of the melting line is found to be negative. It has been found that the density of the Jagla liquid, similarly to that of liquid water, reaches a maximum when heated isobarically[50, 36], a behaviour accurately reproduced by our NS simulations.

Unexpectedly, NS identified a novel, thermodynamically stable, high-pressure solid phase, which has a density at least 30% higher than the liquid. This structure has a primitive cubic unit cell of 8 atoms, and has an $Ia\overline{3}d$ symmetry. We provide an overview of the energetics of the stable crystalline phases of the Jagla solid in Fig. 2.

We have found that the Jagla $Ia\overline{3}d$ crystal is dynamically unstable, as illustrated in the inset Fig. 2 for the displacement of an atom in the $[100]$ direction. As a result, at temperatures above $T=0.22$, the average structure is $Ia\overline{3}d$, but at lower temperatures the structure undergoes two types of distinct symmetric distortions, transforming into lower energy crystalline phases with either $Fddd$ or $C2/c$ symmetries. The density of the $Fddd$ and $C2/c$ phases is slightly lower than that of the $Ia\overline{3}d$, and upon compression the energy of these phases decreases monotonically until the first neighbour shell reaches the hard sphere limit at the interatomic separation of $r_{0}$, where no further compression is possible.

Regarding the low-temperature region of the phase diagram, we found that by isotropically compressing the close-packed hcp and fcc structures, two distinct stability regions can be identified. At low densities, the first neighbour-shell is positioned at the distance of the potential minimum, at $b=1.72r_{0}$, whereas at very high densities the interaction energy is further lowered, reaching the global minimum structure when the interatomic distances at the first neighbour shell are at the hard-sphere limit. Such iso-structural transitions are known for potential models with multiple characteristic distances.[37, 56] We expect that the high-density closed pack phases (HDfcc and HDhcp) are only stable at very low temperatures, as the phase space associated with these highly compressed structures is restricted by the hard-particle repulsive interaction. Though we saw no sign of the fcc stacking or the extreme high density phases in our calculations, the scaling of the structures suggest that both the HDhcp and HDfcc are much lower in enthalpy, as seen in Figure 2. The relative stability of hcp and fcc, or potentially other stacking variants, depend on the choice of potential parameters[62, 63], leading to HDfcc being the ground state structure for the parameterization employed in this work.

Having identified the novel $Ia\overline{3}d$ crystalline structure, we carried out thermodynamic integration and coexistence simulations to accurately determine its region of stability. Both of these methods suggest that from the pressure $p=0.07$, the $Ia\overline{3}d$ structure becomes more stable than the LDhcp phase, with the triple point between the three phases (LDL, LDhcp, $Ia\overline{3}d$) estimated in the temperature region $T=0.36-0.37$. According to our calculations, above the triple point pressure the melting temperature of the $Ia\overline{3}d$ phase increases steeply, and freezing occurs at significantly higher temperatures than the LDL-HDL binodal in the entire pressure range. We have explicitly compared the Gibbs free energy of the HDL to that of the $Ia\overline{3}d$ phase and found that the HDL does not become more stable, even at higher temperatures. Hence the LLCP, in contrary with previous assumptions, is buried deep in the metastable region of the phase diagram.

The microscopic structure of the $Ia\overline{3}d$ phase of the Jagla solid, depicted in the left panels of Fig. 3, is characterised by each atom having three nearest neighbors, which form two interpenetrating networks, each of which is identical to the $K_{4}$ structure proposed for a range of materials, such as boron[64], phosphorus[65] and hydrogenated carbon[66]. A strikingly similar structure has been observed in the ice VI phase of water[67], as shown in the right panels Fig. 3, where the water molecules also form an intertwined hydrogen bonding network. Interestingly, ice VI is one of the two stable crystalline structures which are suspected to hide the HDL and VHDL phases of water on the phase diagram.[68, 69]

We argue that the existence of a stable high density crystalline phase highlights yet another curious qualitative similarity between the Jagla matter and water. Their radial distribution functions (RDF), after rescaling, are compared in Fig. 4. In the case of the Jagla model, there are two competing local structural features, one characterized by a closed-packed first neighbor shell located around the potential minimum, while the other features a more open structure with three neighbors almost at the hard-core contact limit. These are observed individually in the corresponding solid phases, LDhcp and $Ia\overline{3}d$, but appear simultaneously in the liquid phases. This further highlights the analogy of the local structures observed in water and the Jagla matter, and we also note the signature of the intercalated networks appearing as a second peak in the RDF, around $r=1.3$ of Ice VI and the Jagla $Ia\overline{3}d$ phases.

The fact that a new, structurally complex phase was found for such a simple model potential that is thermodynamically stable in a large pressure-temperature region demonstrates that relying on just physical intuition may be insufficient in predicting phase diagrams, and without an exhaustive search of the phase space by appropriate computational tools, crucial properties of the system can remain hidden. Our simulations have identified a previously unexplored crystalline phase, refuting the commonly held presumption that both the high and low density liquid phases of the Jagla matter, and the transition between them are thermodynamically stable. This raises a series of further questions on whether reparametrization of the model could lower the free energy of the high density liquid and move the LLCP in the stable region of the phase diagram. This new, hitherto missing piece of the liquid-liquid polyamorphism puzzle provides additional support, and a potential route of generalization of the two-state model[72] that has been suggested as a possible explanation for the phenomena.