As a liquid approaches the gas state, the properties of the potential energy landscape (PEL) sampled by the system become anomalous. Specifically, (i) the mechanically stable local minima of the PEL [inherent structures (IS)] can exhibit cavitation above the so-called Sastry volume, vS, before the liquid enters the gas phase. In addition, (ii) the pressure of the liquid at the sampled IS [i.e., the PEL equation of state, PIS(v)] develops a spinodal-like minimum at vS. We perform molecular dynamics simulations of a monatomic water-like liquid and verify that points (i) and (ii) hold at high temperatures. However, at low temperatures, cavitation in the liquid and the corresponding IS occurs simultaneously and a Sastry volume cannot be defined. Remarkably, at intermediate/high temperatures, the IS of the liquid can exhibit crystallization, i.e., the liquid regularly visits the regions of the PEL that belong to the crystal state. The model liquid considered also exhibits a liquid–liquid phase transition (LLPT) between a low-density and a high-density liquid (LDL and HDL). By studying the behavior of PIS(v) during the LLPT, we identify a Sastry volume for both LDL and HDL. The HDL Sastry volume marks the onset above which IS are heterogeneous (composed of LDL and HDL particles), analogous to points (i) and (ii) above. However, the relationship between the LDL Sastry volume and the onset of heterogeneous IS is less evident. We conclude by presenting a thermodynamic argument that can explain the behavior of the PEL equation of state PIS(v) across both the liquid–gas phase transition and LLPT.

The potential energy landscape (PEL) approach1–3 is a formalism within statistical mechanics that has been extensively used to study condensed matter.4 Systems studied in the past include low-temperature liquids,5–7 glasses,8–10 granular materials,11 clusters of atoms,12,13 and proteins.14 Briefly, for a system of N atoms, the PEL is the hypersurface in a (3N + 1)-dimensional space defined by the potential energy of system V as a function of 3N degrees of freedom, V(r1,r2,,rN). The PEL approach allows one to characterize quantitatively the dynamic and thermodynamic properties of a system in terms of the topography of V(r1,r2,,rN) {e.g., the distribution of local minima [inherent structures (IS)] in the PEL, curvature of the PEL at the IS, and IS energy and pressure}. In particular, it is possible within the PEL formalism to obtain the equation of state (EOS) P(v) of a given liquid5 based on the topography of V(r1,r2,,rN). Perhaps, one of the most insightful applications of the PEL formalism is in the study of out-of-equilibrium systems, including fast-cooled liquids and glasses.8,15 More recently, the PEL has been extended to the study of quantum liquids.16 

Interestingly, applications of the PEL approach to the study of phase transitions have been limited. Most studies have focused on the liquid–gas phase transition.17–19 In a series of publications, Sastry, Stillinger, Debenedetti, and collaborators studied the liquid–gas transition and show that, as the liquid approaches the gas state, the PEL equation-of-state PIS(v) exhibits a spinodal-like behavior.18–23 Specifically, as the liquid is decompressed at constant temperature, PIS(v) decreases and reaches a spinodal-like minimum at a volume vS called the Sastry volume. Importantly, vS occurs while the system remains in the liquid state, i.e., before the true liquid-to-gas spinodal limit is reached. In particular, it was shown that the IS sampled by the liquid at v < vS were homogeneous while, at v > vS, the IS could also exhibit cavitation (fractured, heterogeneous IS). Hence, at least for finite-size systems,21,22 the Sastry volume signals the volume above which the (stable) liquid is able to sample fractured IS anticipating the liquid-to-gas transition upon further decompression.

The PEL formalism has also been applied to the study of another class of first-order phase transitions, the liquid–liquid phase transition (LLPT). Experiments indicate that a LLPT exists in water24 and aqueous solutions25,26 and other polymorphic systems, such as silicon27 and sulfur;28 computer simulations also suggest that a LLPT exists in high-pressure hydrogen.29 As for the liquid–gas phase transition, it was found that the PEL-EOS PIS(v) across the LLPT is analogous to the corresponding EOS, P(v). For example, PIS(v) exhibits a van der Waals-like loop during the LLPT.15,30–32 In addition, the PEL approach has been shown to be extremely useful in providing a framework to understand the (out-of-equilibrium) thermodynamics of glass polymorphism and its close connection to the LLPT.15,30,31 It was shown that there are common features in the PEL-EOS of the equilibrium liquid across the (a) LLPT and (b) during the associated glass–glass first-order-like phase transition. Hence, the PEL formalism is a potential tool to extend the concept of phase transition, for equilibrium systems, to the case of phase transitions between glasses.

In this work, we explore the PEL of the water-like monatomic Fermi–Jagla (FJ) liquid, which, in addition to the standard liquid–gas phase transition, also exhibits a LLPT. The aim of this work is to (i) revisit the results of Refs. 21 and 22 for the case of the FJ liquid, which exhibits complex phase-behavior; and (ii) to explore whether the concept of Sastry volume associated with the liquid–gas phase transition can be generalized to other first-order phase transitions, specifically the LLPT. An ultimate goal of this work is (iii) to provide a common explanation of the behavior of the PEL-EOS during both the liquid–gas and LLPT. This work is organized as follows: The computer simulation details are provided in Sec. II. The results are presented in Sec. III. We show that (i′) our results at high-temperatures are fully consistent with Refs. 21 and 22. However, the PEL-EOS is different at low temperature and the Sastry volume cannot be identified due to cavitation in the liquid. In addition, we find that the IS of the FJ liquid can exhibit crystallization, even when the liquid is stable and hence, amorphous and homogeneous. (ii′) Our results also indicate that the PEL-EOS during the LLPT shares some similarities with the PEL-EOS during the liquid–gas transition. In particular, it is possible to identify “Sastry volumes” associated with the LLPT. Section III also (iii′) includes a brief discussion that relates equilibrium thermodynamics and the behavior of the PEL-EOS across first-order phase transitions. A summary is included in Sec. IV.

We perform Monte Carlo (MC) simulations of a monatomic system with isotropic pair interactions given by the Fermi–Jagla (FJ) pair potential, U(r).33–35 This is a core-softened potential characterized by two length-scales, a hard-core radius ra, and an attractive minimum at r = b = 1.97 a; see Fig. 1. The functional form of the FJ pair potential is given by

U(r)=ϵ0anrn+A01+expA1A0(r/aA2)B01+expB1B0(r/aB2).
(1)

The parameters Ai and Bi (i = 0, 1, 2) are provided in Table I (the parameters ϵ0 and a are irrelevant since they define, respectively, the units of energy and length). The FJ potential is obtained from the well-studied Jagla potential,36,37UJ(r), by (i) approximating the two linear ramps of UJ(r) with Fermi distributions and (ii) replacing the hard-core repulsion of UJ(r) by a power law repulsion. The main advantage of the pair potential given in Eq. (1), relative to the Jagla potential, is that U(r) is a smooth function of r with a smooth first-order derivative with respect to r. Hence, contrary to the Jagla potential, the FJ pair potential can be implemented in standard MD simulations and is suitable for PEL studies. The thermodynamic, structural, and dynamic properties of the FJ model liquid can be found in Refs. 33–35. Briefly, the FJ model system exhibits many water anomalous properties, such as the presence of an isothermal compressibility maximum upon isobaric cooling and a diffusivity maximum upon isothermal compression. In particular, this model liquid exhibits a first-order LLPT separating a low-density and a high-density liquid state (LDL and HDL) at low temperatures—a phenomenon due to the presence of two-length scales in U(r); see, e.g., Ref. 38. In the P-T plane, the LLPT of the FJ liquid extends toward high temperatures and ends at a liquid–liquid critical point (LLCP) located at Tc = 0.18 and Pc = 0.35 (vc = Vc/N = 2.9), in reduced units.33 Reduced units are defined by setting the atomic mass m = 1 and the Boltzmann constant kB = 1; energies and distances are given in units of ϵ0 and a.

FIG. 1.

Fermi–Jagla pair interaction potential, U(r). The FJ potential is characterized by a hard-core radius ra, a core-softened part at approximately a < r < b ≈ 2a and a weak attractive part at rb. For comparison, we include a Lennard-Jones pair potential with the same minimum depth and location as U(r).

FIG. 1.

Fermi–Jagla pair interaction potential, U(r). The FJ potential is characterized by a hard-core radius ra, a core-softened part at approximately a < r < b ≈ 2a and a weak attractive part at rb. For comparison, we include a Lennard-Jones pair potential with the same minimum depth and location as U(r).

Close modal
TABLE I.

Pair interaction potential parameters used in Eq. (1); see Fig. 1.

nA0A1A2B0B1B2
20 4.56 28.88 1.36 1.00 3.57 2.36 
nA0A1A2B0B1B2
20 4.56 28.88 1.36 1.00 3.57 2.36 

The FJ potential is truncated at the cutoff distance rc = 4.0, and a switching function is added, as implemented in Ref. 39. The switching function introduces minor modifications in the original FJ potential at 3.6 ≤ r ≤ 4.0 only, and makes the potential energy and the corresponding forces smooth functions of r. This is important for the calculation of the IS of the system. For a given system configuration, the corresponding IS is obtained by minimizing the total potential energy using the conjugate gradient method;40 see Ref. 39 for details.

The system is composed of N = 1000 particles located in a cubic box with periodic boundary conditions. MC simulations are performed at constant volume per particle, 2.0 ≤ v ≤ 6.5, and temperature, T = 0.70, 0.50, 0.35, 0.20 > Tc and T = 0.15 < Tc. In a single MC step, all the atoms of the system are subjected to a move try. The system is equilibrated for 106 MC steps followed by 106–5 × 106 steps for data production, depending on v and T. During data production, configurations are saved every 10 000 MC steps, and the configurations are then subjected to potential energy minimization. Hence, for a given (v, T) state, we collect 100–500 independent IS. We also perform molecular dynamics (MD) simulations at T = 0.20 and for large volumes (see below) to confirm that our results are not affected by any potential artifacts that may originate by using the MC simulation technique. Details on the MD simulations can be found in Refs. 30, 33, and 39.

The properties of the PEL of simple liquids [e.g., Lennard-Jones binary mixtures (LJBM)], in the proximity of the liquid–vapor phase transition, have been studied in detail by Sastry, Stillinger, Debenedetti, and collaborators.18,19,21,22 One of the main properties in the PEL formalism is the so-called PEL equation of state (PEL-EOS), PIS(v). PIS(v) is the average pressure of the system at the IS sampled at a given temperature and as a function of volume. It was shown that, as the volume of the liquid increases at constant temperature, approaching the gas state, PIS(v) becomes negative and reaches a minimum at the so-called Sastry volume, vS.18–20 In particular, it was shown that, at v < vS, the liquid can only sample homogeneous IS, while at v > vS, the liquid is able to sample IS that are fractured, i.e., they exhibit cavitation. Therefore, vS signals the onset of mechanically stable configurations (IS) accessible to the liquid that have cavities coexisting with amorphous solid domains.21,22 Importantly, the Sastry volume is smaller than the liquid-to-vapor spinodal vspl, at which the liquid becomes unstable. Hence, the cavitation in the IS occurs at volumes where the liquid is still stable and hence, homogeneous.

In this section, we follow Refs. 18, 19, 21, and 22 and study the PEL-EOS of the FJ liquid as the system expands from the liquid to the gas states. The aim of this section is to test whether the results presented in the binary mixtures of Refs. 21 and 22 hold for the rather complex (monodisperse) FJ model liquid. Briefly, we find that the picture presented in those works holds at high temperatures (T = 0.70) (Sec. III A 1). Surprisingly, it is also found that crystallization in the IS can occur while the system is in the liquid state (Sec. III A 2). The role of temperature on the nature of the IS sampled by the liquid is also discussed (Sec. III A 3).

1. High-temperatures

Most of the studies that focused on the PEL-EOS of simple liquids close to the liquid–vapor phase transition have been performed at high temperatures; e.g., T = 1.0 in the case of LJBM.21,22 Accordingly, we discuss first the PEL-EOS of the FJ liquid at T = 0.70; the role of temperature is discussed in Sec. III A 3. The results presented in this section are fully consistent with Refs. 21 and 22 and hence, we will only present a brief discussion; we refer the reader to Refs. 21 and 22 for a detailed discussion.

The equation of state of the FJ liquid at T = 0.70, P(v), is shown in Fig. 2(a). In the (large) range of volumes considered, the P(v) of the liquid decreases monotonically until the liquid–vapor spinodal volume is reached, vspl = 5.65, at which P(v) exhibits a minimum. As it is shown below, cavitation in the liquid occurs only at v > vspl. Figure 2(b) shows the PEL-EOS at T = 0.70. As for the case of P(v), PIS(v) also decreases monotonically with increasing volumes and reaches a minimum at the Sastry volume, vS = 5.23 ± 0.03. In particular, it follows that for the FJ liquid at T = 0.70, vS < vspl. The results in Fig. 2 are fully consistent with the findings in Refs. 21 and 22 for simple binary mixtures of similar size (N ≈ 1000 atoms).

FIG. 2.

(a) Pressure and (b) inherent structure pressure of the FJ liquid as a function of volume at T = 0.70. The liquid-to-vapor spinodal occurs at vspl = 5.65, above which the liquid becomes unstable and exhibits cavitation. The Sastry volume is vS = 5.23 ± 0.03 < vspl and hence, it occurs within the equilibrium liquid domain. Error bars in (b) represent the standard deviations in PIS(v). Fluctuations in PIS(v) become visible at volumes slightly smaller than vS.

FIG. 2.

(a) Pressure and (b) inherent structure pressure of the FJ liquid as a function of volume at T = 0.70. The liquid-to-vapor spinodal occurs at vspl = 5.65, above which the liquid becomes unstable and exhibits cavitation. The Sastry volume is vS = 5.23 ± 0.03 < vspl and hence, it occurs within the equilibrium liquid domain. Error bars in (b) represent the standard deviations in PIS(v). Fluctuations in PIS(v) become visible at volumes slightly smaller than vS.

Close modal

Next, we show that, consistent with Refs. 21 and 22, the Sastry volume signals approximately the onset above which the liquid IS may exhibit cavitation. To do this, we implement the methodology of Ref. 41 and identify liquid–vapor interfaces that may enclose cavities in a given instantaneous configuration (from the MC simulations) or IS. Then, we identify the atoms of the system that are within a distance d = 2.2 of such interfaces (this distance corresponds to the first minimum of the radial distribution function at large volumes). We define fLG as the fraction of atoms in the system that belongs to the interface of a cavity. The behavior of fLG(v) at T = 0.70, calculated from the instantaneous configurations and IS are shown in Fig. 3. Consistent with Refs. 21 and 22, we find that fLGIS(v)=0 for IS at approximately v < vS while fLGIS(v) increases rapidly as v increases above vS. Hence, practically, only the IS at v > vS may show cavitation. Instead, in the instantaneous configurations, fLG(v) ≈ 0 for all volumes v < vsp,meaning that the system is always homogeneous in the stable/metastable liquid states; fLG(v) > 0 at approximately v > vspl. A snapshot of the system at an IS obtained at v = 5.40, slightly above vS, is shown in Fig. 4.

FIG. 3.

Fraction of particles of the FJ liquid at the interface with empty cavities, fLG. The blue line corresponds to the fLG(v) calculated from the instantaneous configurations (MC simulation); the red line is the fLGIS(v) calculated from the corresponding IS. Cavitation in the IS occurs above the Sastry volume, vS ≈ 5.23 (dashed line), at which fLGIS(v)>0. Cavitation in the instantaneous liquid configurations occurs above the liquid–gas spinodal, vspl = 5.65 > vS.

FIG. 3.

Fraction of particles of the FJ liquid at the interface with empty cavities, fLG. The blue line corresponds to the fLG(v) calculated from the instantaneous configurations (MC simulation); the red line is the fLGIS(v) calculated from the corresponding IS. Cavitation in the IS occurs above the Sastry volume, vS ≈ 5.23 (dashed line), at which fLGIS(v)>0. Cavitation in the instantaneous liquid configurations occurs above the liquid–gas spinodal, vspl = 5.65 > vS.

Close modal
FIG. 4.

(a) Snapshot of an IS of the FJ liquids at T = 0.70 and v = 5.40 > vS (eIS = −4.97, Q6 = 0.06). At approximately v > vS, the liquid can sample IS that are fractured, i.e., with cavities (while the liquid remains homogeneous). (b) Cavity formed in the IS shown in (a); the cavity is split into four pieces due to periodic boundary conditions.

FIG. 4.

(a) Snapshot of an IS of the FJ liquids at T = 0.70 and v = 5.40 > vS (eIS = −4.97, Q6 = 0.06). At approximately v > vS, the liquid can sample IS that are fractured, i.e., with cavities (while the liquid remains homogeneous). (b) Cavity formed in the IS shown in (a); the cavity is split into four pieces due to periodic boundary conditions.

Close modal

It was also shown in Refs. 21 and 22 that at v < vS the distribution of IS energies, eIS, is unimodal, while at v > vS, the distribution of eIS splits into two separate subsets. The large-eIS subset of IS was found to correspond to IS that were homogeneous or liquid-like, while the low-eIS subset of IS corresponded to IS that were inhomogeneous or fractured, with amorphous domains coexisting with empty cavities (e.g., Fig. 4). We confirm that this finding also holds in the case of the FJ liquid. Surprisingly, we also find that a small subset of the IS shows crystallization (Sec. III A 2).

Figure 5 shows the distribution of IS energies Pb(eis) at selected volumes for v < vS. At first sight, Pb(eis) seems unimodal, but there are a few apparently residual IS that exhibit very small values of eIS. For example, at v = 3.90, Pb(eis) shows a large maximum at eIS ≈ −5.28 and a small subset of IS with energies about eIS ≈ −5.37 (see also Fig. S1 of the supplementary material). As we will show in Sec. III A 2, the small subset of IS with eIS ≈ −5.37 corresponds to IS that crystallized during the potential energy minimization [the corresponding instantaneous configurations are all liquid-like (amorphous)]. Here, we focus on the majority of the IS sampled by the liquid that correspond to the large maximum in Pb(eis) (e.g., at eIS ≈ −5.28 for v = 3.90).

FIG. 5.

Distribution probability of the IS energies sampled by the liquid, Pb(eIS), at T = 0.70. Upper panels: at v < vS, Pb(eIS) is practically unimodal, but there is a small subset of IS sampled by the liquid with an anomalous small IS energy. The high-eIS, large subset of these IS corresponds to liquid-like, homogeneous IS; the low-eIS, small subset of these IS corresponds to IS that exhibit crystallization (see also Fig. 7). Lower panels: at v > vS, Pb(eIS) is practically bimodal. The high-eIS subset of these IS corresponds to amorphous, homogeneous IS; the low-eIS subset of these IS are fractured (i.e., they exhibit cavitation; see also Fig. 4). As the liquid approaches the gas state (v > vS), more IS are fractured and less IS are homogeneous.

FIG. 5.

Distribution probability of the IS energies sampled by the liquid, Pb(eIS), at T = 0.70. Upper panels: at v < vS, Pb(eIS) is practically unimodal, but there is a small subset of IS sampled by the liquid with an anomalous small IS energy. The high-eIS, large subset of these IS corresponds to liquid-like, homogeneous IS; the low-eIS, small subset of these IS corresponds to IS that exhibit crystallization (see also Fig. 7). Lower panels: at v > vS, Pb(eIS) is practically bimodal. The high-eIS subset of these IS corresponds to amorphous, homogeneous IS; the low-eIS subset of these IS are fractured (i.e., they exhibit cavitation; see also Fig. 4). As the liquid approaches the gas state (v > vS), more IS are fractured and less IS are homogeneous.

Close modal

If the small set of IS that exhibits crystallization is removed, then Pb(eis) is unimodal. This maximum in Pb(eis) (associated with the homogeneous liquid-like IS) shifts smoothly toward negative values as v increases in the range v = 3.5–4.0, and then, it shifts toward larger values at v = 4.00–6.5. Consistent with Ref. 21, at v > vS, Pb(eis) develops a second maximum at lower energies. For example, at v = 5.30, the large maximum in Pb(eis) is located at eIS ≈ −4.90, while a small maximum develops at eIS ≈ −4.95. This extra subset of IS with eIS ≈ −4.95 corresponds to fractured IS. The subset of IS with eIS ≈ −4.90, corresponds to the liquid-like homogeneous IS. It follows from Fig. 5 that, as the volume increases (v > vS), the number of homogeneous IS decreases while the number of fractured IS increases. Accordingly, at very large volumes, only those IS with low-eIS remain, implying that all IS sampled by the liquid exhibit cavitation. The average IS energy of the FJ liquid as a function of volume, EIS(v), is shown in Fig. 6. The fact that Pb(eis) becomes bimodal at approximately v > vS is reflected in the sudden increase in the error bars of EIS(v) at vvS.

FIG. 6.

Average energy of the IS sampled by the FJ liquid at T = 0.70 and as a function of volume. The dashed-line indicates the Sastry volume above which the IS may show cavitation. The presence of fractured IS leads to the enhanced fluctuations in EIS(v) at vvS.

FIG. 6.

Average energy of the IS sampled by the FJ liquid at T = 0.70 and as a function of volume. The dashed-line indicates the Sastry volume above which the IS may show cavitation. The presence of fractured IS leads to the enhanced fluctuations in EIS(v) at vvS.

Close modal

2. Crystallization in the IS

Next, we focus on that small subset of IS that exhibits very small IS energy (eIS ≈ −5.37 at v = 3.90; Fig. 5). Surprisingly, these IS, which are obtained from a liquid configuration, are largely crystallized; see Figs. 7(a) and 7(b).42 To confirm this, we show in Fig. 8 the global orientational order parameter Q6 defined in Ref. 43 as a function of eIS, for selected volumes. Each point in the plot for Q6(eIS) corresponds to a single IS sampled by the liquid at the corresponding volume (T = 0.70). It follows that the anomalous IS with very low energy are also characterized by large orientational order, Q6 > 0.4. Since typical crystal configurations have Q6 > 0.20, these IS contain indeed crystalline domains (see Refs. 34 and 35). Crystallized IS are found at v ≥ 3.60 (T = 0.70), although they represent a very small subset of the IS sampled by the liquid. Interestingly, close to the Sastry volume, we also find a few crystallized IS that are also fractured, with very small cavities surrounded by crystalline domains; see Fig. 7(c).

FIG. 7.

Snapshots of the IS with unusually low energy eIS; see Fig. 8. (a) An IS that is homogeneous and fully crystallized (v = 3.90 < vS, eIS = −5.36, Q6 = 0.44). (b) An IS that is partially crystallized (v = 4.50 < vS, eIS = −5.21, Q6 = 0.21). (c) An IS that is crystallized and contains a small cavity (v = 5.16 < vS, eIS = −5.08, Q6 = 0.46) and (d) Instantaneous configuration corresponding to the IS shown in (c) [potential energy of −3.75 and Q6 = 0.06; similar instantaneous configurations are found for (a) and (b)].

FIG. 7.

Snapshots of the IS with unusually low energy eIS; see Fig. 8. (a) An IS that is homogeneous and fully crystallized (v = 3.90 < vS, eIS = −5.36, Q6 = 0.44). (b) An IS that is partially crystallized (v = 4.50 < vS, eIS = −5.21, Q6 = 0.21). (c) An IS that is crystallized and contains a small cavity (v = 5.16 < vS, eIS = −5.08, Q6 = 0.46) and (d) Instantaneous configuration corresponding to the IS shown in (c) [potential energy of −3.75 and Q6 = 0.06; similar instantaneous configurations are found for (a) and (b)].

Close modal
FIG. 8.

Global orientational order parameter Q6 as a function of the IS energy eIS for each of the 500 IS sampled by the liquid at T = 0.70 and selected volumes. Top panels: vvS; bottom panels: v > vS. Circles indicate IS that are homogeneous; fractured IS are indicated by red squares. IS (fractured or homogeneous) with large Q6 > 0.20–0.30 show crystallization. At v < vS, IS are homogeneous with a few exceptions at 5.10 < v < vS; at v > vS, the number of fractured IS increases with increasing v.

FIG. 8.

Global orientational order parameter Q6 as a function of the IS energy eIS for each of the 500 IS sampled by the liquid at T = 0.70 and selected volumes. Top panels: vvS; bottom panels: v > vS. Circles indicate IS that are homogeneous; fractured IS are indicated by red squares. IS (fractured or homogeneous) with large Q6 > 0.20–0.30 show crystallization. At v < vS, IS are homogeneous with a few exceptions at 5.10 < v < vS; at v > vS, the number of fractured IS increases with increasing v.

Close modal

The physical picture that results from the PEL formalism is very elegant and intuitive. At high temperatures and at approximately v < vS, the liquid samples only IS that are homogeneous. These IS are mostly amorphous, liquid-like, but they can also show crystallization. In other words, the liquid always has access to those regions of the PEL that belong to the crystal state. However, due to the large kinetic energy available, the system is able to escape the regions of the PEL associated with the crystal state. As shown below, only at low temperatures does the system get trapped in the regions of the PEL characteristic to the crystal state. At v > vS, the system can sample homogeneous IS, but it also starts to sample IS that are fractured.21,22 Most of the fractured IS correspond to an amorphous solid with cavities; however, a small number of IS show crystalline domains coexisting with cavities. As the volume of the system increases, above vS, (i) the number of homogeneous IS decreases, and (ii) the number of IS with crystalline domains (homogeneous or fractured) also decreases. At very large volumes, only fractured, amorphous IS are accessible to the system.

The picture described above is nicely visualized in Fig. 9, where we show eIS as a function of the configuration number n (T = 0.70). Here, n = 0, 1, …, 499 with n = 0 and n = 499 corresponding to the first and last configurations obtained during the MC simulation at a given v. At relatively small volumes, v = 3.50 < vS, the IS energy of the system exhibits small fluctuations during the MC simulation, since most of the IS are liquid-like (homogeneous); see upper panels in Fig. 9. With increasing volumes, but still below the Sastry volume, the system also samples IS that exhibit crystallization. These crystallized IS correspond to the sharp spikes in Fig. 9 for 3.90 ≤ v ≤ 5.10. The frequency of the exploration of crystalline IS varies with volume. At v = 3.50, no crystallized IS are accessible to the liquid. As the volume increases toward v = 4.30, more basins belonging to the crystal state are reached by the system, but these IS become less common at v = 5.10, as the liquid approaches the Sastry volume.

FIG. 9.

Energy of the IS sampled by the system during MC simulations performed at T = 0.70. n is the configuration number with n = 0 and n = 499 being the first and last configurations obtained during the MC simulations. Upper panels: eIS(n) at v < vs. At these volumes, the system mostly visits basins that are liquid-like and homogeneous, but it can also visit deep (low-eIS) basins of the PEL that correspond to IS with crystalline order (spikes). The frequency at which the liquid visits the basins belonging to the crystal region of the PEL varies with v. Lower panels: eIS(n) at v > vs. At these volumes, the system visits IS that are either homogeneous (mostly amorphous) or fractured. At v = 5.20–5.40, most of the IS are homogeneous, and the system visits the fractured IS sporadically (spikes). At large volumes, the situation is reversed; most IS are fractured and homogeneous IS are visited sporadically (spikes at v = 5.60, 5.80).

FIG. 9.

Energy of the IS sampled by the system during MC simulations performed at T = 0.70. n is the configuration number with n = 0 and n = 499 being the first and last configurations obtained during the MC simulations. Upper panels: eIS(n) at v < vs. At these volumes, the system mostly visits basins that are liquid-like and homogeneous, but it can also visit deep (low-eIS) basins of the PEL that correspond to IS with crystalline order (spikes). The frequency at which the liquid visits the basins belonging to the crystal region of the PEL varies with v. Lower panels: eIS(n) at v > vs. At these volumes, the system visits IS that are either homogeneous (mostly amorphous) or fractured. At v = 5.20–5.40, most of the IS are homogeneous, and the system visits the fractured IS sporadically (spikes). At large volumes, the situation is reversed; most IS are fractured and homogeneous IS are visited sporadically (spikes at v = 5.60, 5.80).

Close modal

At v > vS, the crystalline basins are rarely accessed, and the system spends most of the MC simulations sampling either homogeneous, amorphous IS or fractured IS. At v = 5.20, most of the basins visited by the liquid correspond to homogeneous IS; the spikes shown in Fig. 9 at v = 5.20 correspond to fractured IS. The number of fractured IS increases as v increases (see spikes in the panels of Fig. 9 for v = 5.30 and 5.40). Upon further increase in the volume, the situation reverses, and the system explores mostly fractured IS with infrequent visits to the homogeneous IS (spikes in Fig. 9 for v = 5.50–5.80).

3. The role of temperature

The results discussed in Sec. III A 2 are found to depend on temperature. In particular, the Sastry volume cannot be identified at low temperatures due to fast cavitation in the instantaneous configurations of the liquid. To show this, we include in Figs. 10(a)10(c) the P(v), PIS(v), and EIS(v) of the system at T = 0.70, 0.50, 0.35, and 0.20. Figure 10(a) shows that, at T < 0.50, the liquid-to-gas spinodal becomes inaccessible in our MC simulations, and the minimum in P(v) associated with the liquid-to-vapor spinodal is preempted by a sudden jump at a volume v*(T). Not surprisingly, the instantaneous configurations at v > v*(T) exhibit cavitation, while those at v < v*(T) remain homogeneous. The cavitation in the instantaneous configurations at v > v*(T) remains in the corresponding IS. Accordingly, a jump in PIS(v) also occurs at v = v*(T) (T = 0.50); see Fig. 10(b).

FIG. 10.

(a) Pressure, (b) IS pressure, and (c) IS energy of the FJ liquid as a function of volume across the LLPT at T = 0.70, 0.50, 0.35, 0.20 from MC simulations. For clarity, the curves of PIS(v) for T = 0.50, 0.35, 0.20 are shifted by ΔP = 0.35, 0.70, 1.05, respectively; curves of EIS(v) for T = 0.50, 0.35, 0.20 are shifted by ΔE = 0.2, 0.4, 0.6, respectively. For comparison, also included are the results obtained at T = 0.20 from MD simulations.

FIG. 10.

(a) Pressure, (b) IS pressure, and (c) IS energy of the FJ liquid as a function of volume across the LLPT at T = 0.70, 0.50, 0.35, 0.20 from MC simulations. For clarity, the curves of PIS(v) for T = 0.50, 0.35, 0.20 are shifted by ΔP = 0.35, 0.70, 1.05, respectively; curves of EIS(v) for T = 0.50, 0.35, 0.20 are shifted by ΔE = 0.2, 0.4, 0.6, respectively. For comparison, also included are the results obtained at T = 0.20 from MD simulations.

Close modal

As shown in Fig. 10(b), a Sastry volume can be identified only at T = 0.70 and 0.50. At T = 0.35, 0.20, there is no minimum in PIS(v) [v*(T) < vS(T)] and the Sastry volume is preempted by cavitation in the instantaneous configurations of the liquid. We confirmed these results by performing MD simulations of the FJ liquid at T = 0.20. The values of P(v), PIS(v), and EIS(v) obtained from our MD simulations practically overlap with the corresponding results from MC simulations (Fig. 10).

At T = 0.50, 0.35, 0.20, we also observe sporadic crystallization in the IS such as the case of T = 0.70 (while the system remains in the liquid state). Our analysis of crystallization at these temperatures presents similar results to those discussed in Sec. III A 2 for T = 0.70. At even lower temperatures, T = 0.15, the liquid becomes unstable, and it crystallizes at v > 3.60. It follows that there is no Sastry density at T = 0.15.

In this section, we explore the PEL-EOS of the FJ liquid during the LLPT at T = 0.15 < Tc = 0.18. In particular, we extend the concept of Sastry density/volume (originally, defined for the liquid–gas phase transition) to the phase transition between LDL and HDL. We also explore the nature of the IS sampled by the liquid during the LLPT and whether there is nucleation-and-growth of one liquid phase as the LLPT proceeds from one liquid state to the other. Briefly, we show that (i) Sastry volumes can be associated with both LDL and HDL. However, it is less evident whether these Sastry volumes mark the onset of IS sampled by the system that are heterogeneous (i.e., composed of LDL and HDL domains). Our results also show that (ii) the IS sampled across the LLPT show a large LDL cluster with embedded small HDL domains (at large volumes), or a large HDL cluster with small LDL domains (at small volumes). In addition, (iii) the small difference in energy between the IS sampled by the LDL and HDL leads to a smooth variation of EIS(v) during the LLPT, and a unimodal distribution Pb(eIS) at all volumes across the LLPT. Overall, results (ii) and (iii) indicate that the PEL description of the LLPT in the FJ liquid is reminiscent of the liquid–gas transition of binary mixtures where the attraction between particles is weak.22 

1. Sastry volumes for LDL and HDL

In order to identify Sastry volumes for LDL and HDL, we include in Fig. 11 the P(v) and PIS(v) of the FJ liquid at T = 0.15. At this temperature, crystallization occurs at volumes v > 3.6.33 The van der Waals loop in P(v), due to the LLPT, occurs at v < 3.6, i.e., just below the crystallization volumes. Fractured crystals, i.e., crystalline structures that exhibit cavitation are only found at v > 5.25. While the changes in P(v) and PIS(v) during crystallization are rather weak [Figs. 11(a) and 11(b)], crystallization leads to a sudden decrease in EIS(v) at v ≈ 3.60 [Fig. 11(c)] (crystallization is usually accompanied by the exploration of deep IS in the PEL).

FIG. 11.

(a) Pressure and (b) IS pressure as a function of volume across the LLPT at T = 0.15. The dashed-blue line in (a) is obtained by applying Maxwell’s construction from which the coexistence volumes vHDL = 2.51 and vLDL = 3.32 are obtained. A similar construction is included in (b) (dashed-magenta line; vHDLIS=2.20 and vLDLIS=3.55). Insets show P(v) and PIS(v) over the volumes at which crystallization occurs, approximately v ≥ 3.6; cavitation is observed at approximately v > 5.25. (c) Average IS energy EIS(v) of the liquid as a function of volume. The sudden drop in EIS at approximately v > 3.6 is due to the crystallization of the system.

FIG. 11.

(a) Pressure and (b) IS pressure as a function of volume across the LLPT at T = 0.15. The dashed-blue line in (a) is obtained by applying Maxwell’s construction from which the coexistence volumes vHDL = 2.51 and vLDL = 3.32 are obtained. A similar construction is included in (b) (dashed-magenta line; vHDLIS=2.20 and vLDLIS=3.55). Insets show P(v) and PIS(v) over the volumes at which crystallization occurs, approximately v ≥ 3.6; cavitation is observed at approximately v > 5.25. (c) Average IS energy EIS(v) of the liquid as a function of volume. The sudden drop in EIS at approximately v > 3.6 is due to the crystallization of the system.

Close modal

The maximum and minimum of P(v) define, respectively, the LDL-to-HDL and HDL-to-LDL spinodal volumes, vspl,LDL = 3.15 and vspl,HDL = 2.65. Analogous to Fig. 2(b), we associate the maximum and minimum in PIS(v) with the LDL-to-HDL and HDL-to-LDL Sastry volume, respectively. It follows from Fig. 11(b) that the Sastry volume associated with the HDL-to-LDL transition is vS,HDL ≈ 2.57, while the Sastry volume associated with the LDL-to-HDL transition is vS,LDL ≈ 3.27. Interestingly, vS,HDL < vspl,HDL, meaning that, as the HDL is decompressed and approaches the LDL state, a minimum develops, first, in PIS(v) and then, in P(v) (upon further decompression). Hence, as for the liquid-to-gas phase transition, the PEL-EOS anticipates the HDL-to-LDL transition that occurs in the system upon isothermal decompression. A similar situation occurs as LDL is compressed and approaches the HDL state, i.e., vS,LDL > vspl,HDL. In other words, as LDL is compressed, a maximum develops, first, in PIS(v) and then, in P(v). Again, the PEL-EOS anticipates the LDL-to-HDL transition that occurs in the system upon isothermal compression. We note that although one can identify Sastry volumes for LDL and HDL, the difference between the Sastry volume of LDL/HDL and the corresponding spinodal volume is not large, |vS,LDLvspl,LDL|, |vS,HDLvspl,HDL| ≈ 0.1. For comparison, in the case of the liquid-to-gas transition, we find |vSvspl| ≈ 0.3 (T = 0.70).

One may wonder whether the Sastry volumes for LDL and HDL have a physical meaning analogous to the Sastry volume associated with the liquid-to-gas phase transition. For example, if one decompresses HDL, are the IS sampled by the system homogeneous (with only HDL-like particles) at v < vS,HDL? and, in particular, are there heterogeneous IS (with LDL- and HDL-like particles) at v > vS,HDL? To address these questions, we classify the particles in the system as HDL and LDL and identify the volumes at which both LDL and HDL particles coexist in the IS sampled by the system. We classify the particles as HDL (LDL) if they have n ≥ 1 (n = 0) nearest-neighbors within a distance r < dc. As explained in the supplementary material, we find that a reasonable value for this cutoff distance is dc = 1.12. The average number of LDL and HDL particles at the IS sampled by the system during the LLPT, nLDLIS(v) and nHDLIS(v), is shown in Fig. 12(a); for comparison, the average number of LDL and HDL particles in the corresponding instantaneous configurations is shown in the supplementary material.

FIG. 12.

(a) Average number of LDL and HDL particles in the IS sampled by the system during the LLPT at T = 0.15, nLDLIS and nHDLIS. (b) Average number of LDL and HDL clusters in the IS sampled by the system at T = 0.15, mLDLIS and mHDLIS. Black dashed-lines indicate the Sastry volumes for HDL (left; vS,HDL = 2.57) and LDL (right; vS,LDL = 3.27). The red dashed-lines mark the volume range within which PIS(v) exhibits a van der Waals-like loop [see the Maxwell construction in Fig. 11(b)].

FIG. 12.

(a) Average number of LDL and HDL particles in the IS sampled by the system during the LLPT at T = 0.15, nLDLIS and nHDLIS. (b) Average number of LDL and HDL clusters in the IS sampled by the system at T = 0.15, mLDLIS and mHDLIS. Black dashed-lines indicate the Sastry volumes for HDL (left; vS,HDL = 2.57) and LDL (right; vS,LDL = 3.27). The red dashed-lines mark the volume range within which PIS(v) exhibits a van der Waals-like loop [see the Maxwell construction in Fig. 11(b)].

Close modal

Figure 12(a) shows that at small volumes, nLDLIS(v)=0 since, as expected, the liquid is in the stable HDL state. nLDLIS(v) increases sharply during the HDL-to-LDL transition and reached nLDLIS(v)=N=1000 at large volumes, as expected. In particular, we find that nLDLIS(v) starts to increase at approximately v > vS,HDL [see the left dashed-line in Fig. 12(a)]. Hence, the Sastry volume of HDL, indeed, marks the onset above which the system samples IS that are heterogeneous, containing both LDL and HDL particles. It follows that the Sastry volume of HDL, vS,HDL, plays a similar role to the Sastry volume vS defined in Fig. 2(b) in the case of the liquid-to-gas transition. The situation is somewhat different when considering LDL. While at large volumes nHDLIS(v)=0, as expected, nLDLIS(v) starts to increase at roughly v = 3.6. This volume is considerably larger than the Sastry volume of LDL (vS,LDL = 3.27; right dashed line). It follows that the Sastry volume of LDL does not mark the onset of volumes below which the IS become heterogeneous. We stress that our results are based on the specific classification of LDL and HDL particles given above. Our definitions of LDL and HDL particles are the simplest that one can imagine since they depend only on the interparticle separation, and they involve a single parameter dc. Yet, it is possible that alternative more complex definitions of LDL and HDL particles involving, e.g., orientational order parameters, may lead to slightly different results. In our LDL/HDL particle definitions, LDL particles have zero nearest-neighbors within a distance r < dc = 1.12. Yet, as shown in the supplementary material, the radial distribution function for the LDL liquid at the coexistence volume v = 3.3 ≈ vLDL (T = 0.15), while very small, is not strictly zero. This may lead to spurious HDL particles in the LDL liquid at v > vS,LDL (see the next paragraph).

We conclude this section with a brief description of the structure of the IS sampled by the system during the LLPT. To do so, we include in Fig. 12(b) the average number of clusters composed of LDL and HDL particles found in the IS sampled by the system, mLDLIS(v) and mHDLIS(v). Two HDL (LDL) particles are considered to belong to the same HDL (LDL) cluster if their separation is rdc=2.2. The values of dc correspond approximately to the second minimum of the radial distribution function of HDL and LDL at v = 2.5 ≈ vHDL and v = 3.3 ≈ vLDL, respectively; see Fig. S3(a) of the supplementary material. The main point of Fig. 12(b) is that there are coexisting LDL and HDL clusters at ∼2.4 ≤ v ≤ 3.7. Interestingly, as shown in Fig. 11(b), this is close to the volume interval that one obtains by applying the Maxwell construction to PIS(v) (vHDLIS=2.20v3.55=vLDLIS). This suggests a close connection between the nature of the IS and the PEL-EOS, and provides a potential meaning to the “Maxwell’s construction” applied to the PEL-EOS. It follows that the volume range at which the IS are composed of LDL and HDL clusters expands beyond the volumes between the Sastry volumes of LDL and HDL, vS,HDL < v < vS,LDL,

Figure 12(b) shows that mLDLIS(v)>1 only at 2.4 < v < 3.1. At these volumes, mHDLIS(v)=1. Similarly, mHDLIS(v)>1 only at 3.2 < v < 3.7, at which mLDLIS(v)=1. It follows that the heterogeneous IS are composed of a single cluster of particles belonging to a dominant species (LDL or HDL) that contains a small number (<25) of clusters of the minority species (HDL or LDL). Snapshots of the system at v = 2.6, 2.8, 3.0, and 3.2 are shown in Fig. 13.

FIG. 13.

Snapshots of IS sampled by the system at T = 0.15 and (a) v = 2.60, (b) v = 2.80, (c) v = 3.0, and (d) v = 3.20 across the LLPT (see also Fig. 12). Red and blue spheres correspond to HDL and LDL particles, respectively.

FIG. 13.

Snapshots of IS sampled by the system at T = 0.15 and (a) v = 2.60, (b) v = 2.80, (c) v = 3.0, and (d) v = 3.20 across the LLPT (see also Fig. 12). Red and blue spheres correspond to HDL and LDL particles, respectively.

Close modal

We note that Fig. 12(b) does not show a clear relationship between the Sastry volume of HDL and LDL (black dashed lines) and the volumes at which mLDLIS(v) and mHDLIS(v) deviate from zero. In other words, vs,HDL and vs,LDL are not related to the onset of LDL and HDL clusters in the IS visited by the system. We also looked at the average and maximum LDL and HDL cluster sizes during the LLPT and, again, we could not find a clear relationship between these quantities and the vs,HDL and vs,LDL.

2. IS energies of LDL and HDL

A comparison of Figs. 4 and 13 shows a very contrasting picture during the liquid–gas and LLPT. For example, upon expansion of HDL, there is no clear nucleus of LDL that develops across the LLPT. Instead, the LDL clusters are rather small and are dispersed throughout the sample. (i) Figure 4 is consistent with the snapshots shown in Ref. 21 during the liquid-to-gas transition for binary mixtures with strong attractive interactions. Instead, (ii) the snapshots of Fig. 13 are reminiscent of the snapshots shown in Ref. 22 during the liquid-to-gas transition for binary mixtures with weak attractive interactions.

In Ref. 22, it was shown that, the different nature of the IS sampled by binary mixtures [points (i) and (ii) above] was related to the IS energies available to the liquid and gas states. Specifically, the IS energies of the liquid and gas states of binary mixtures with strong attractive interactions were very different, which led to a bi-modal Pb(eIS) during the liquid-to-gas transition. Indeed, this is consistent with our findings in Fig. 5. Instead, in the case of binary mixtures with weak attractive interactions, the IS energies of the liquid and gas states are similar, and Pb(eIS) remains unimodal during the liquid-to-gas transition.

The distributions of IS energies sampled by the FJ liquid during the LLPT are shown in Fig. 14. Interestingly, Pb(eIS) remains unimodal during the LLPT, consistent with point (ii) above. Hence, it is apparent that the presence of IS with scattered small clusters (of one phase in another) during the corresponding phase transition is due to the similar IS energies accessible to the two phases involved in the corresponding (LDL-HDL or liquid–gas) phase transition.

FIG. 14.

Distribution probability of the IS energies sampled by the liquid, Pb(eIS), during the LLPT at T = 0.15. At all volumes, Pb(eIS) remains unimodal.

FIG. 14.

Distribution probability of the IS energies sampled by the liquid, Pb(eIS), during the LLPT at T = 0.15. At all volumes, Pb(eIS) remains unimodal.

Close modal

The liquid–liquid and liquid–gas first-order phase transitions are characterized by a range of volumes at which P/vT>0; see Figs. 2(a) and 11(a). At these states, the homogeneous system is thermodynamically unstable and hence it must phase-separate.44 Interestingly, the liquid–gas and liquid–liquid phase transitions are also accompanied by a range of volumes where PIS/vT>0; see, e.g., Figs. 2(b) and 11(b). This behavior in PIS(v) is not common; for example, in a system of soft-spheres, a model liquid with no liquid–liquid or liquid–gas phase transitions, one finds that at low temperatures PIS(v) decays monotonically with increasing v (see Ref. 5, Sec. 7.2). It is apparent that some connection exists between the thermodynamic conditions of stability [i.e., that P/vT>0 corresponds to an unstable system] and the anomalous region of the PEL-EOS where PIS/vT>0.

In previous studies, we argued that the anomalous behavior of the PEL-EOS across the LLPT can be understood in the context of thermodynamics.15,30,31 The argument is as follows: In the PEL approach, a liquid in equilibrium can be considered to be composed of two subsystems,4,5,45 a subsystem that depends solely on EIS and Pb(eIS) (IS-subsystem), and a subsystem that depends solely on the vibrational motion of the system within the basins of the PEL (vib-subsystem). The subsystems are characterized by Helmholtz free energies FIS and Fvib, and the total free energy of the liquid can be expressed as F = FIS + Fvib. Since, P=F/vT, it follows that P = PIS + Pvib, where

PIS=FIS/vT
(2)

and

Pvib=Fvib/vT.
(3)

Now, at low temperatures and constant volume, and within the harmonic approximation, PvibT4,5,46 and hence, PPIS. Accordingly, under these conditions, both P(v) andPIS(v) must exhibit anomalous behavior with positive slope during the LLPT [Figs. 11(a) and 11(b)] and liquid-to-gas transition [Figs. 2(a) and 2(b)]. Similar conclusions apply to the case of binary mixtures studied in Refs. 21 and 22, where PIS/vT>0 across the liquid-to-gas phase transition.

A similar argument may be used to explain why EIS(v) exhibits a concave region [2EIS/v2T<0] across the liquid–liquid and liquid–gas phase transitions [Figs. 6 and 11(c)]. This is also an anomalous property of the PEL; for example, in a system of soft spheres, where no liquid–liquid or liquid–gas phase transitions exist, 2EIS/v2T>0. During a first-order phase transition, volumes for which 2F/v2T<0 correspond to unstable states, where the system cannot remain homogeneous. Since F = ET S, where E and S are the internal energy and entropy of the system, respectively, one may expect that FE at low temperatures. If so, a concavity in F(v) should be accompanied with a concavity in E(v) (during a first-order phase transition).

Now, in the PEL approach, E = EIS + Evib, where EIS and Evib are the energy of the IS- and vib-subsystems, respectively. Moreover, in the harmonic approximation,5EvibT and hence, at low temperatures, FEEIS. It follows that, under these hypotheses, a liquid–liquid and liquid–gas first-order phase transition must be accompanied by an anomalous concave region in EIS(v). This may explain why, for the FJ liquid, EIS(v) exhibits a concave region during both the LLPT [Fig. 11(c)] and liquid-to-gas transition (Fig. 6). Similar conclusions apply to the case of the binary mixtures studied in Refs. 21–23, where 2EIS/v2T<0 across the liquid-to-gas phase transition.

In this work, we performed computer simulations of the FJ liquid to study the relationship between (i) the PEL-EOS PIS(v) at constant temperature and (ii) the nature of the corresponding IS sampled by the system. The focus of this study was on the liquid–gas and liquid–liquid first-order phase transitions.

The case of the liquid–gas phase transition was addressed in previous publications by Sastry, Stillinger, and Debenedetti.18,19,21,22 Consistent with those works, we find that at high temperatures (T = 0.50, 0.70); see Figs. 2 and 10, the PEL-EOS decreases with increasing volumes and reached a minimum at the Sastry volume. Importantly, the Sastry volume occurs at states where the liquid is stable and has not reached its spinodal limit. We also confirm that at high temperatures, the Sastry volume signals the onset volume above which the system starts to sample fractured (heterogeneous) IS where the liquid coexists with voids (e.g., Fig. 4). Hence, as found in the works by Sastry, Stillinger, and Debenedetti,18,19,21,22 the Sastry volume represents a spinodal-like limit at which PIS(v) has a minimum. Above the Sastry volume, the fractured IS exhibit strong voids akin to those found in Refs. 21 and 22 for the case of binary mixtures with strong attractive interactions. Consistent with those works, we find that the homogeneous (v < vS) and fractured IS (v > vS) have well separated IS energies (Fig. 5). Overall, our results at T = 0.70 obtained with the (monodisperse) FJ liquid are fully consistent with those shown in Refs. 21 and 22 for binary mixtures with strong attractive interactions.

By extending our studies to lower temperatures (T = 0.35, 0.20), we find that the Sastry volume cannot be defined. At low temperatures, the cavitation in the liquid becomes unavoidable and, as expected, cavitation in the liquid leads to fractured IS upon energy minimization. At the volumes accessible to the homogeneous liquid, PIS(v) is a monotonically decreasing function of v (Fig. 10).

An interesting finding of this study is the presence of crystallization in some of the IS sampled by the equilibrium liquid. Crystallization in the IS sampled by the liquid state has not been observed in the binary mixtures of Refs. 21 and 22, although in Ref. 47, it is reported that a few IS of a Lennard-Jones binary mixture, at specific conditions, exhibit small crystalline domains coexisting with amorphous regions. A comparison of our results at T = 0.20–0.70 gives an intuitive picture of crystallization within the PEL formalism in the case of a bad glass former. Specifically, at all temperatures, the equilibrium liquid is able to sample IS that belong to the crystal region of the PEL. While at high temperature the system always escapes from the crystal basins of the PEL (Fig. 9), the probability to escape from such basins decreases with decreasing temperatures. At low temperatures, T = 0.15, the system becomes trapped in the crystal basins of the PEL as the crystal becomes the stable state.

One property that makes the FJ liquid an interesting model system is that it presents a liquid–liquid phase transition at low temperatures as for the case of water and silicon, among others. One of the main aims of this work was to study the PEL-EOS during the LLPT at T = 0.15 as well as its relationship with the nature (homogeneous/heterogeneous) of the IS sampled by the system. As shown in Fig. 11, the PEL-EOS exhibits spinodal-like extrema during the LLPT. Such extrema are analogous to the Sastry volume vS (which was originally defined in the context of the liquid–gas phase) and hence, can be identified as the Sastry volumes for the LDL and HDL states. The HDL Sastry volume marks the onset above which IS are heterogeneous (composed of LDL and HDL particles) analogous to the liquid–gas phase transition case. However, the relationship between the LDL Sastry volume and the onset of heterogeneous IS structures is less evident and deserves further study. In particular, it would be important in the future to explore how our results (e.g., Fig. 12) are affected if alternative definitions are used to identify the LDL and HDL particles in the system. The role of system size on the PEL-EOS is also important and should be further studied in the future. Indeed, as shown in Refs. 21 and 22, the PEL-EOS can change with the system size if atoms in the liquid have strong attractive interactions.

The heterogeneous IS sampled by the system during the LLPT exhibit either (i) a single, large cluster of HDL particles that contains scattered, small clusters composed of LDL clusters; or (ii) a single, large LDL cluster with embedded small HDL clusters. Hence, the heterogeneous IS sampled during the LLPT are very different from those sampled in the case of the liquid–gas transition. Indeed, the IS sampled during the LLPT are reminiscent of the IS explored by binary mixtures with weak attractive interactions during the liquid–gas phase transition. In both cases, the interactions between the atoms belonging to the two phases (LDL and HDL in the case of the LLPT, and the liquid and gas phases in the case of the liquid-to-gas transition) are similar in strength. The different nature of the IS sampled by the FJ liquid during the liquid–gas and LLPT can also be related to the difference in the IS energies of the two phases involved in the phase transition (liquid and gas, or LDL and HDL). Specifically, during the liquid–gas transition, the IS energies of the gas and liquid are well separated (see bimodal distributions in Fig. 5); these IS exhibit large voids. Instead, during the LLPT, the IS energies of LDL and HDL are indistinguishable (see unimodal distributions in Fig. 14); these IS exhibit only small domains of LDL/HDL in an HDL/LDL background.

Overall, our study intends to unify the behavior of the PEL-EOS across (a) the liquid–gas and (b) liquid–liquid first-order phase transition in a polymorphic model liquid. In doing so, we extend the ideas and concepts, such as the Sastry volume, introduced previously to describe the liquid–gas phase transition, to the case of the LLPT. Accordingly, in the last section of this work, we also presented a simple thermodynamic argument that can explain the behavior of PIS(v) and EIS(v) during both the liquid–gas and liquid–liquid phase transitions. Specifically, we showed that, under simple assumptions, thermodynamic arguments may explain the anomalous behavior of PIS(v) and EIS(v) across both the liquid–gas transition and LLPT.

See the supplementary material for additional information on the definition of LDL and HDL particles.

This work was supported by the SCORE Program of the National Institutes of Health under Award No. 1SC3GM139673. We thank the NSF CREST Center for Interface Design and Engineered Assembly of Low-Dimensional systems (IDEALS), NSF Grant No. HRD-1547830, for additional support. This work was also supported, in part, by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grant Nos. CNS-0855217, CNS-0958379, ACI-1126113, and OEC-2215760.

The authors have no conflicts to disclose.

Yang Zhou: Data curation (lead); Software (lead); Visualization (lead); Writing – review & editing (equal). Gustavo E. Lopez: Conceptualization (lead); Funding acquisition (lead); Writing – review & editing (equal). Nicolas Giovambattista: Conceptualization (lead); Funding acquisition (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

1.
M.
Goldstein
,
J. Chem. Phys.
51
,
3728
(
1969
).
2.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. A
25
,
978
(
1982
).
3.
F. H.
Stillinger
and
T. A.
Weber
,
Science
225
,
983
(
1984
).
4.
F. H.
Stillinger
,
Energy Landscapes, Inherent Structures, and Condensed-Matter Phenomena
(
Princeton University Press
,
Princeton
,
2016
).
5.
F.
Sciortino
,
J. Stat. Mech.
2005
,
P05015
.
6.
P. G.
Debenedetti
and
F. H.
Stillinger
,
Nature
410
,
259
(
2001
).
7.
I.
Saika-Voivod
,
F.
Sciortino
, and
P. H.
Poole
,
Phys. Rev. E
69
,
041503
(
2004
).
8.
S.
Sastry
,
P. G.
Debenedetti
, and
F. H.
Stillinger
,
Nature
393
,
554
(
1998
).
9.
A.
Heuer
,
J. Phys.: Condens. Matter
20
,
373101
(
2008
).
10.
N.
Giovambattista
,
F. W.
Starr
, and
P. H.
Poole
,
J. Chem. Phys.
150
,
224502
(
2019
).
11.
A. J.
Liu
and
S. R.
Nagel
,
Annu. Rev. Condens. Matter Phys.
1
,
347
(
2010
).
12.
D. J.
Wales
,
J. Chem. Phys.
142
,
130901
(
2015
).
13.
S.
Bonfanti
and
W.
Kob
,
J. Chem. Phys.
147
,
204104
(
2017
).
14.
D.
Wales
,
Energy Landscapes: Applications to Clusters, Biomolecules and Glasses
(
Cambridge University Press
,
2004
).
15.
P. H.
Handle
,
F.
Sciortino
, and
N.
Giovambattista
,
J. Chem. Phys.
150
,
244506
(
2019
).
16.
N.
Giovambattista
and
G. E.
Lopez
,
Phys. Rev. Res.
2
,
043441
(
2020
).
17.
D. S.
Corti
 et al,
Phys. Rev. E
55
,
5522
(
1997
).
18.
S.
Sastry
,
P. G.
Debenedetti
, and
F. H.
Stillinger
,
Phys. Rev. E
56
,
5533
(
1997
).
20.
R. A.
LaViolette
,
Phys. Rev. B
40
,
9952
(
1989
).
21.
Y. E.
Altabet
,
F. H.
Stillinger
, and
P. G.
Debenedetti
,
J. Chem. Phys.
145
,
211905
(
2016
).
22.
Y. E.
Altabet
,
A. L.
Fenley
,
F. H.
Stillinger
, and
P. G.
Debenedetti
,
J. Chem. Phys.
148
,
114501
(
2018
).
23.
J.
Hernández-Rojas
and
D. J.
Wales
,
Phys. Rev. B
68
,
144202
(
2003
).
24.
K. H.
Kim
,
K.
Amann-Winkel
,
N.
Giovambattista
 et al,
Science
370
,
978
982
(
2020
).
25.
Y.
Suzuki
,
Proc. Natl. Acad. Sci. U. S. A.
119
,
e2113411119
(
2022
).
26.
J.
Bachler
,
V.
Fuentes-Landete
,
D. A.
Jahn
,
J.
Wong
,
N.
Giovambattista
, and
T.
Loerting
,
Phys. Chem. Chem. Phys.
18
,
11058
(
2016
).
27.
M.
Beye
,
F.
Sorgenfrei
,
W. F.
Schlotter
,
W.
Wurth
, and
A.
Föhlisch
,
Proc. Natl. Acad. Sci. U. S. A.
107
,
16772
16776
(
2010
).
28.
L.
Henry
,
M.
Mezouar
,
G.
Garbarino
 et al,
Nature
584
,
382
386
(
2020
).
29.
C.
Pierleoni
,
M. A.
Morales
,
G.
Rillo
,
M.
Holzmann
, and
D. M.
Ceperley
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
4953
4957
(
2016
).
30.
G.
Sun
,
L.
Xu
, and
N.
Giovambattista
,
Phys. Rev. Lett.
120
,
035701
(
2018
).
31.
N.
Giovambattista
,
F.
Sciortino
,
F. W.
Starr
, and
P. H.
Poole
,
J. Chem. Phys.
145
,
224501
(
2016
).
32.
Y.
Liu
,
G.
Sun
, and
L.
Xu
,
J. Chem. Phys.
154
,
134503
(
2021
).
33.
J. Y.
Abraham
,
S. V.
Buldyrev
, and
N.
Giovambattista
,
J. Phys. Chem. B
115
,
14229
(
2011
).
34.
S.
Reisman
and
N.
Giovambattista
,
J. Chem. Phys.
138
,
064509
(
2013
).
35.
A.
Gordon
and
N.
Giovambattista
,
Phys. Rev. Lett.
112
,
145701
(
2014
).
36.
E. A.
Jagla
,
J. Chem. Phys.
111
,
8980
(
1999
).
37.
E. A.
Jagla
,
Phys. Rev. E
63
,
061501
(
2001
).
38.
S. V.
Buldyrev
,
G.
Malescio
,
C. A.
Angell
,
N.
Giovambattista
,
S.
Prestipino
,
F.
Saija
,
H. E.
Stanley
, and
L.
Xu
,
J. Phys.: Condens. Matter
21
,
504106
(
2009
).
39.
G.
Sun
,
L.
Xu
, and
N.
Giovambattista
,
J. Chem. Phys.
146
,
014503
(
2017
).
40.
W. H.
Press
,
B. P.
Flannery
,
A. A.
Teukolsky
 et al,
Numerical Recipes: The Art of Scientific Computing
(
Cambridge University Press
,
Cambridge, UK
,
1986
).
41.
A. P.
Williard
and
D.
Chandler
,
J. Phys. Chem. B
114
,
1954
(
2010
).
42.

We confirmed that IS sampled by the liquid can indeed show crystallized structures by performing potential energy minimzations using two independent codes.

43.
P. J.
Steinhardt
,
D. R.
Nelson
, and
M.
Ronchetti
,
Phys. Rev. B
28
,
784
(
1983
).
44.
H. B.
Callen
,
Thermodynamics and an Introduction to Thermostatistics
(
John Wiley & Sons, Inc.
,
New York
,
1985
).
45.
P. G.
Debenedetti
 et al,
J. Phys. Chem. B
103
,
7390
(
1999
).
46.
47.
S.
Sarkar
and
B.
Bagchi
,
Phys. Rev. E
83
,
031506
(
2011
).

Supplementary Material