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, *v*_{S}, 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, *P*_{IS}(*v*)] develops a spinodal-like minimum at *v*_{S}. 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 *P*_{IS}(*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 *P*_{IS}(*v*) across both the liquid–gas phase transition and LLPT.

## I. INTRODUCTION

The potential energy landscape (PEL) approach^{1–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 (3*N* + 1)-dimensional space defined by the potential energy of system V as a function of 3*N* degrees of freedom, $V(r\u20d71,r\u20d72,\u2026,r\u20d7N)$. The PEL approach allows one to characterize quantitatively the dynamic and thermodynamic properties of a system in terms of the topography of $V(r\u20d71,r\u20d72,\u2026,r\u20d7N)$ {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 liquid^{5} based on the topography of $V(r\u20d71,r\u20d72,\u2026,r\u20d7N)$. 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 *P*_{IS}(*v*) exhibits a spinodal-like behavior.^{18–23} Specifically, as the liquid is decompressed at constant temperature, *P*_{IS}(*v*) decreases and reaches a spinodal-like minimum at a volume *v*_{S} called the Sastry volume. Importantly, *v*_{S} 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* < *v*_{S} were homogeneous while, at *v* > *v*_{S}, 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 water^{24} and aqueous solutions^{25,26} and other polymorphic systems, such as silicon^{27} 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 *P*_{IS}(*v*) across the LLPT is analogous to the corresponding EOS, *P*(*v*). For example, *P*_{IS}(*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.

## II. SIMULATION DETAILS

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 *r* ≈ *a*, and an attractive minimum at *r* = *b* = 1.97 *a*; see Fig. 1. The functional form of the FJ pair potential is given by

The parameters *A*_{i} and *B*_{i} (*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,37}*U*_{J}(*r*), by (i) approximating the two linear ramps of *U*_{J}(*r*) with Fermi distributions and (ii) replacing the hard-core repulsion of *U*_{J}(*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 *T*_{c} = 0.18 and *P*_{c} = 0.35 (*v*_{c} = *V*_{c}/*N* = 2.9), in reduced units.^{33} Reduced units are defined by setting the atomic mass *m* = 1 and the Boltzmann constant *k*_{B} = 1; energies and distances are given in units of *ϵ*_{0} and *a*.

n
. | A_{0}
. | A_{1}
. | A_{2}
. | B_{0}
. | B_{1}
. | B_{2}
. |
---|---|---|---|---|---|---|

20 | 4.56 | 28.88 | 1.36 | 1.00 | 3.57 | 2.36 |

n
. | A_{0}
. | A_{1}
. | A_{2}
. | B_{0}
. | B_{1}
. | B_{2}
. |
---|---|---|---|---|---|---|

20 | 4.56 | 28.88 | 1.36 | 1.00 | 3.57 | 2.36 |

The FJ potential is truncated at the cutoff distance *r*_{c} = 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 > *T*_{c} and *T* = 0.15 < *T*_{c}. In a single MC step, all the atoms of the system are subjected to a move try. The system is equilibrated for 10^{6} MC steps followed by 10^{6}–5 × 10^{6} 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.

## III. RESULTS

### A. The liquid–gas phase transition in the PEL formalism

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), *P*_{IS}(*v*). *P*_{IS}(*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, *P*_{IS}(*v*) becomes negative and reaches a minimum at the so-called Sastry volume, *v*_{S}.^{18–20} In particular, it was shown that, at *v* < *v*_{S}, the liquid can only sample homogeneous IS, while at *v* > *v*_{S}, the liquid is able to sample IS that are fractured, i.e., they exhibit cavitation. Therefore, *v*_{S} 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 *v*_{spl}, 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, *v*_{spl} = 5.65, at which *P*(*v*) exhibits a minimum. As it is shown below, cavitation in the liquid occurs only at *v* > *v*_{spl}. Figure 2(b) shows the PEL-EOS at *T* = 0.70. As for the case of *P*(*v*), *P*_{IS}(*v*) also decreases monotonically with increasing volumes and reaches a minimum at the Sastry volume, *v*_{S} = 5.23 ± 0.03. In particular, it follows that for the FJ liquid at *T* = 0.70, *v*_{S} < *v*_{spl}. 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).

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 *f*_{LG} as the fraction of atoms in the system that belongs to the interface of a cavity. The behavior of *f*_{LG}(*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* < *v*_{S} while $fLGIS(v)$ increases rapidly as *v* increases above *v*_{S}. Hence, practically, only the IS at *v* > *v*_{S} may show cavitation. Instead, in the instantaneous configurations, *f*_{LG}(*v*) ≈ 0 for all volumes *v* < *v*_{sp},meaning that the system is always homogeneous in the stable/metastable liquid states; *f*_{LG}(*v*) > 0 at approximately *v* > *v*_{spl}. A snapshot of the system at an IS obtained at *v* = 5.40, slightly above *v*_{S}, is shown in Fig. 4.

It was also shown in Refs. 21 and 22 that at *v* < *v*_{S} the distribution of IS energies, *e*_{IS}, is unimodal, while at *v* > *v*_{S}, the distribution of *e*_{IS} splits into two separate subsets. The large-*e*_{IS} subset of IS was found to correspond to IS that were homogeneous or liquid-like, while the low-*e*_{IS} 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 *P*_{b}(*e*_{is}) at selected volumes for *v* < *v*_{S}. At first sight, *P*_{b}(*e*_{is}) seems unimodal, but there are a few apparently residual IS that exhibit very small values of *e*_{IS}. For example, at *v* = 3.90, *P*_{b}(*e*_{is}) shows a large maximum at *e*_{IS} ≈ −5.28 and a small subset of IS with energies about *e*_{IS} ≈ −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 *e*_{IS} ≈ −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 *P*_{b}(*e*_{is}) (e.g., at *e*_{IS} ≈ −5.28 for *v* = 3.90).

If the small set of IS that exhibits crystallization is removed, then *P*_{b}(*e*_{is}) is unimodal. This maximum in *P*_{b}(*e*_{is}) (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* > *v*_{S}, *P*_{b}(*e*_{is}) develops a second maximum at lower energies. For example, at *v* = 5.30, the large maximum in *P*_{b}(*e*_{is}) is located at *e*_{IS} ≈ −4.90, while a small maximum develops at *e*_{IS} ≈ −4.95. This extra subset of IS with *e*_{IS} ≈ −4.95 corresponds to fractured IS. The subset of IS with *e*_{IS} ≈ −4.90, corresponds to the liquid-like homogeneous IS. It follows from Fig. 5 that, as the volume increases (*v* > *v*_{S}), the number of homogeneous IS decreases while the number of fractured IS increases. Accordingly, at very large volumes, only those IS with low-*e*_{IS} remain, implying that all IS sampled by the liquid exhibit cavitation. The average IS energy of the FJ liquid as a function of volume, *E*_{IS}(*v*), is shown in Fig. 6. The fact that *P*_{b}(*e*_{is}) becomes bimodal at approximately *v* > *v*_{S} is reflected in the sudden increase in the error bars of *E*_{IS}(*v*) at *v* ≥ *v*_{S}.

#### 2. Crystallization in the IS

Next, we focus on that small subset of IS that exhibits very small IS energy (*e*_{IS} ≈ −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 *Q*_{6} defined in Ref. 43 as a function of *e*_{IS}, for selected volumes. Each point in the plot for *Q*_{6}(*e*_{IS}) 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, *Q*_{6} > 0.4. Since typical crystal configurations have *Q*_{6} > 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).

The physical picture that results from the PEL formalism is very elegant and intuitive. At *high* temperatures and at approximately *v* < *v*_{S}, 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* > *v*_{S}, 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 *v*_{S}, (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 *e*_{IS} 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 < *v*_{S}, 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.

At *v* > *v*_{S}, 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*), *P*_{IS}(*v*), and *E*_{IS}(*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 *P*_{IS}(*v*) also occurs at *v* = *v**(*T*) (*T* = 0.50); see Fig. 10(b).

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 *P*_{IS}(*v*) [*v**(*T*) < *v*_{S}(*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*), *P*_{IS}(*v*), and *E*_{IS}(*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.

### B. The liquid–liquid phase transition in the PEL formalism

In this section, we explore the PEL-EOS of the FJ liquid during the LLPT at *T* = 0.15 < *T*_{c} = 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 *E*_{IS}(*v*) during the LLPT, and a unimodal distribution *P*_{b}(*e*_{IS}) 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 *P*_{IS}(*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 *P*_{IS}(*v*) during crystallization are rather weak [Figs. 11(a) and 11(b)], crystallization leads to a sudden decrease in *E*_{IS}(*v*) at *v* ≈ 3.60 [Fig. 11(c)] (crystallization is usually accompanied by the exploration of deep IS in the PEL).

The maximum and minimum of *P*(*v*) define, respectively, the LDL-to-HDL and HDL-to-LDL spinodal volumes, *v*_{spl,LDL} = 3.15 and *v*_{spl,HDL} = 2.65. Analogous to Fig. 2(b), we associate the maximum and minimum in *P*_{IS}(*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 *v*_{S,HDL} ≈ 2.57, while the Sastry volume associated with the LDL-to-HDL transition is *v*_{S,LDL} ≈ 3.27. Interestingly, *v*_{S,HDL} < *v*_{spl,HDL}, meaning that, as the HDL is decompressed and approaches the LDL state, a minimum develops, first, in *P*_{IS}(*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., *v*_{S,LDL} > *v*_{spl,HDL}. In other words, as LDL is compressed, a maximum develops, first, in *P*_{IS}(*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, |*v*_{S,LDL} − *v*_{spl,LDL}|, |*v*_{S,HDL} − *v*_{spl,HDL}| ≈ 0.1. For comparison, in the case of the liquid-to-gas transition, we find |*v*_{S} − *v*_{spl}| ≈ 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* < *v*_{S,HDL}? and, in particular, are there heterogeneous IS (with LDL- and HDL-like particles) at *v* > *v*_{S,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* < *d*_{c}. As explained in the supplementary material, we find that a reasonable value for this cutoff distance is *d*_{c} = 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.

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* > *v*_{S,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, *v*_{S,HDL}, plays a similar role to the Sastry volume *v*_{S} 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 (*v*_{S,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 *d*_{c}. 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* < *d*_{c} = 1.12. Yet, as shown in the supplementary material, the radial distribution function for the LDL liquid at the coexistence volume *v* = 3.3 ≈ *v*_{LDL} (*T* = 0.15), while very small, is not strictly zero. This may lead to spurious HDL particles in the LDL liquid at *v* > *v*_{S,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 $r\u2264dc\u2032=2.2$. The values of $dc\u2032$ correspond approximately to the second minimum of the radial distribution function of HDL and LDL at *v* = 2.5 ≈ *v*_{HDL} and *v* = 3.3 ≈ *v*_{LDL}, 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 *P*_{IS}(*v*) $(vHDLIS=2.20\u2264v\u22643.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, *v*_{S,HDL} < *v* < *v*_{S,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.

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, *v*_{s,HDL} and *v*_{s,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 *v*_{s,HDL} and *v*_{s,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 *P*_{b}(*e*_{IS}) 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 *P*_{b}(*e*_{IS}) 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, *P*_{b}(*e*_{IS}) 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.

### C. Thermodynamics and the PEL-EOS during the liquid–liquid and liquid–gas first-order phase transitions

The liquid–liquid and liquid–gas first-order phase transitions are characterized by a range of volumes at which $\u2202P/\u2202vT>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 $\u2202PIS/\u2202vT>0$; see, e.g., Figs. 2(b) and 11(b). This behavior in *P*_{IS}(*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 *P*_{IS}(*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 $\u2202P/\u2202vT>0$ corresponds to an unstable system] and the anomalous region of the PEL-EOS where $\u2202PIS/\u2202vT>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 *E*_{IS} and *P*_{b}(*e*_{IS}) (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 *F*_{IS} and *F*_{vib}, and the total free energy of the liquid can be expressed as *F* = *F*_{IS} + *F*_{vib}. Since, $P=\u2212\u2202F/\u2202vT$, it follows that *P* = *P*_{IS} + *P*_{vib}, where

and

Now, at low temperatures and constant volume, and within the harmonic approximation, *P*_{vib} ∝ *T*^{4,5,46} and hence, *P* ≈ *P*_{IS}. Accordingly, under these conditions, both *P*(*v*) *and* *P*_{IS}(*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 $\u2202PIS/\u2202vT>0$ across the liquid-to-gas phase transition.

A similar argument may be used to explain why *E*_{IS}(*v*) exhibits a concave region $[\u22022EIS/\u2202v2T<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, $\u22022EIS/\u2202v2T>0$. During a first-order phase transition, volumes for which $\u22022F/\u2202v2T<0$ correspond to unstable states, where the system cannot remain homogeneous. Since *F* = *E* − *T S*, where *E* and *S* are the internal energy and entropy of the system, respectively, one may expect that *F* ≈ *E* 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* = *E*_{IS} + *E*_{vib,} where *E*_{IS} and *E*_{vib} are the energy of the IS- and vib-subsystems, respectively. Moreover, in the harmonic approximation,^{5} *E*_{vib} ∝ *T* and hence, at low temperatures, *F* ≈ *E* ≈ *E*_{IS}. 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 *E*_{IS}(*v*). This may explain why, for the FJ liquid, *E*_{IS}(*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 $\u22022EIS/\u2202v2T<0$ across the liquid-to-gas phase transition.

## IV. SUMMARY AND DISCUSSION

In this work, we performed computer simulations of the FJ liquid to study the relationship between (i) the PEL-EOS *P*_{IS}(*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 *P*_{IS}(*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* < *v*_{S}) and fractured IS (*v* > *v*_{S}) 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, *P*_{IS}(*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 *v*_{S} (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 *P*_{IS}(*v*) and *E*_{IS}(*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 *P*_{IS}(*v*) and *E*_{IS}(*v*) across both the liquid–gas transition and LLPT.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

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