Despite the importance of ice nucleation, this process has been barely explored at negative pressures. Here, we study homogeneous ice nucleation in stretched water by means of molecular dynamics seeding simulations using the TIP4P/Ice model. We observe that the critical nucleus size, interfacial free energy, free energy barrier, and nucleation rate barely change between isobars from −2600 to 500 bars when they are represented as a function of supercooling. This allows us to identify universal empirical expressions for homogeneous ice nucleation in the pressure range from −2600 to 500 bars. We show that this universal behavior arises from the pressure dependence of the interfacial free energy, which we compute by means of the mold integration technique, finding a shallow minimum around −2000 bars. Likewise, we show that the change in the interfacial free energy with pressure is proportional to the excess entropy and the slope of the melting line, exhibiting in the latter a reentrant behavior also at the same negative pressure. Finally, we estimate the excess internal energy and the excess entropy of the ice Ih–water interface.

## I. INTRODUCTION

Water crystallization is an essential phase transition in nature and technology. However, in the cryopreservation of biological samples,^{1,2} ice formation can be disastrous. The low temperature preserves the biological material but causes the water within the sample to be in a metastable state subject to crystallization.^{3} Interestingly, by keeping the sample under high pressure, ice nuclei are less likely to form, keeping water liquid for a longer time.^{4–6} The frequency of the nucleation process is mainly determined by the thermodynamic driving force and the cost of creating the interface between the emerging nucleus and the metastable liquid. The reason why high pressure slows down the nucleation process is that the difference in the chemical potential between ice and water, Δ*μ*, which represents the thermodynamic driving force, barely changes between isobars as a function of supercooling, whereas the cost of creating the interface notably increases.^{7} The interfacial free energy is the variable that quantifies this cost. At coexistence, through a planar interface, the interfacial free energy *γ*_{m} differs from the value for a critical nucleus *γ* due to the curvature of the surface.^{8} Nevertheless, both notably increase at high pressure.^{7}

Homogeneous ice nucleation at standard and high pressure has been extensively explored.^{7,9–21} However, ice nucleation in water under negative pressure, i.e., stretched water, has caught less attention.^{22–24} This process is relevant in porous media containing water solutions^{25} and also in water transpiration inside plants^{26} where negative pressure occurs. Creating and maintaining negative pressure over a sample is non-trivial in experiments.^{27} This is because a liquid at negative pressure is metastable.^{28,29} In general, this metastability is considered with respect to the vapor phase,^{30,31} although at certain conditions, it also can be metastable with respect to ice. Some ingenious approaches to create negative pressure in metastable water include the use of a Berthelot tube,^{32} centrifugation,^{33} and, more recently, the use of acoustic waves.^{27,31,34} In contrast, in computer simulations, it is straightforward to work under negative pressures.

In this work, we investigate how ice nucleation properties are affected by negative pressure at different degrees of supercooling. In fact, we find little effect when pressure changes from strongly negative to moderately positive. We investigate the role of the interfacial free energy since it is a key property in determining the phase behavior of water at high pressure.^{7} We find that the slope of the melting line is crucial to describe the change with pressure of the interfacial free energy, which displays a shallow minimum at negative pressure. Our study is based on molecular dynamics simulations with the TIP4P/Ice model,^{35} which has been extensively used to describe ice nucleation^{7,9,12,23} and growth^{36,37} as well as in supercooled water.^{38,39} In particular, we employ the seeding technique^{9,40} to study nucleation and the mold integration technique^{41} to measure the interfacial free energy at coexistence.

## II. SIMULATION METHODS

All simulations have been done with the GROMACS package (4.6.7-version in double precision) with the TIP4P/Ice water model. The simulations are performed in the isothermal–isobaric (NpT) ensemble with a time step of 2 fs using the Noose–Hoover thermostat^{42,43} and the Parrinello–Rahman barostat,^{44} both with a relaxation time of 0.5 ps. Electrostatic interactions are accounted for via the particle-mesh-Ewald summation algorithm^{45} with order 4 and a Fourier spacing of 0.1 nm. The cutoff for the Lennard-Jones and the Coulombic interactions is set to 0.9 nm, and the long-range corrections to the Lennard-Jones part of the potential are included in energy and pressure.

To study nucleation, we use the seeding technique,^{9,46,47} which involves the combination of molecular dynamics simulations and classical nucleation theory (CNT).^{48,49} This technique is based on the behavior of a critical nucleus, which has equal probability of growing and melting when surrounded by the metastable phase at the critical pressure and temperature. In practice, one inserts a spherical ice Ih seed in metastable water and then keeps track of the time evolution of the size of the cluster. One can vary *T*, *p*, and the seed size in order to find at which conditions a certain nucleus size is critical (*N*_{c}). Once *N*_{c} is known, CNT is used to find the interfacial free energy *γ*, the barrier height Δ*G*_{c}, and the nucleation rate *J*. Our system sizes ranged between 80 000 and 250 000 water molecules in total. The duration of the trajectories is between 40 and 115 ns.

It is important to note that in the Gibbsian description of interfaces, one has two bulk phases separated by a dividing surface. However, there is some arbitrariness in the location of the dividing surface, which also affects the interfacial free energy *γ* when the interface has a curvature.^{50–53} Within the CNT framework, the relevant dividing surface is the surface of tension.^{54,55} Besides, the accuracy of the seeding technique relies on the right selection of the order parameter given that *N*_{c} is very sensitive to that choice. In order to approximate the number of molecules in the critical nucleus, we employ an empirical approach that has been successfully applied in crystal nucleation for a large variety of systems.^{7,23,40,41,56,57} In this approach, the averaged Steinhardt bond order parameter,^{58} $q\u03046(T,p)$, is used in combination with the mislabeling criterion^{9} to identify ice-like and water-like molecules. We obtain $q\u03046(T,p)$ for each molecule. The molecules with $q\u03046(T,p)$ above a certain threshold $q\u03046,t(T,p)$ are labeled as ice, whereas those below are labeled as liquid. This threshold depends weakly on the considered thermodynamic range covering pressures from −2600 to −1000 bars and temperatures from 250 to 270 K [see the supplementary material in Ref. 23 for the isothermal change in $q\u03046,t(T,p)$ with pressure]. In this work, the value of the threshold determined by the mislabeling criterion changes between 0.365 for the highest temperature and pressure to 0.385 for the lowest temperature and pressure. Two molecules labeled as solid belong to the same solid cluster if their distance is smaller than 3.5 Å. However, due to the intrinsic arbitrariness of the crystallinity definition and the statistical limitations of the seeding method, it remains some uncertainty in *N*_{c}, which is difficult to estimate.

Once *N*_{c} is known, we employ the CNT equations^{48,49} to determine other important parameters. The interfacial free energy *γ* is given as

where *N*_{c} is the size of the critical nucleus, *ρ*_{ice} is the number density of ice Ih in the bulk at the metastable conditions at which the nucleus is critical, and |Δ*μ*| is known as the driving force to nucleation, i.e., the difference in chemical potential between the liquid and ice phases in the bulk at the conditions that cause the nucleus to be critical. This property can be obtained by thermodynamic integration along an isobar,^{59}

where *k*_{B} is the Boltzmann constant, *T*_{m} is the melting temperature, and *H* is the enthalpy, which can be obtained from simulations of bulk ice Ih and bulk water along the isobar of interest.

Then, the free energy barrier is given as

which allows us to obtain the nucleation rate *J*, the number of critical nuclei forming per unit of time and volume. According to CNT, *J* is given as

where *f*^{+} is the attachment rate that can be approximated through the expression^{7,23}

where *D* is the diffusion coefficient of the metastable liquid and *λ* is a characteristic length, the typical distance that a water molecule covers in order to attach into the nucleus, whose value is ∼3.8 Å for water.^{7,23}

To find the ice Ih–water interfacial free energy at coexistence for a planar interface, *γ*_{m}, we use the mold integration technique,^{41} which consists in computing the reversible work *W* that is necessary to form a crystal slab within a liquid at coexistence. This work is related to the interfacial free energy at coexistence, *γ*_{m}, by *W* = 2 *Aγ*_{m}, where *A* is the interfacial area, and the number 2 accounts for the two interfaces of the slab. The slab formation is induced by switching on an attractive interaction between the mold of potential energy wells and the particles of the initial liquid. The wells are arranged in the equilibrium positions of the oxygen atoms in the ice facet under investigation at coexistence conditions, i.e., for temperatures and pressures located along the ice Ih–water equilibrium line for the TIP4P/Ice water model. First, one has to obtain $\gamma rw$, which is given as

where *r*_{w} indicates the radius of the potential wells and *ϵ* is their energy (with maximum depth equal to *ϵ*_{w}). *N*_{w} is the number of wells in the mold, and ⟨*N*(*ϵ*)⟩ is the average number of occupied wells at a given potential depth *ϵ*. The integration needs to be reversible. To ensure this, thermodynamic integration is performed for wells whose radius is larger than a certain value $rw0$. At $rw0$, the slab is fully formed and the stability no longer depends on the mold–liquid interactions, hence leading to potentially irreversible ice growth. However, since this is the radius that recovers the actual value of *γ*_{m}, thermodynamic integration is repeated for several values of $rw>rw0$ and then $\gamma rw$ is extrapolated to its value at $rw0$ giving *γ*_{m}. For further details on the implementation of this technique we refer the reader to Refs. 41 and 60.

## III. RESULTS

### A. Universality in ice nucleation variables at negative and moderate pressure

First, we study nucleation along the isobars of −2600, −2000, and −1000 bars by means of the seeding approach. For pressures below −3000 bars, we observed spontaneous cavitation occurring within the timescale of the trajectories needed in the seeding method. We obtain the critical nucleus size *N*_{c}, the driving force to nucleation |Δ*μ*|, the interfacial free energy *γ*, the free energy barrier to nucleation Δ*G*_{c}, and the nucleation rate *J*. These results are presented in Table I. As can be seen, even though the pressure significantly differs, the results are surprisingly similar for nuclei of similar size for equivalent supercoolings. This behavior is considerably different from what has been found when comparing the nucleation scenario of normal vs high pressure (i.e., 1 vs 2000 bars; Ref. 7), where the increase in pressure brings down the ice nucleation rate.

N_{c}
. | T (K)
. | ΔT (K)
. | p (bars)
. | ρ_{w} (g/cm^{3})
. | ρ_{ice} (g/cm^{3})
. | |Δμ| (kJ/mol)
. | γ (mJ/m^{2})
. | ΔG_{c} (kJ/mol)
. | log_{10} (J[m^{−3} s^{−1}])
. |
---|---|---|---|---|---|---|---|---|---|

1650 | 255 | 23 | −1000 | 0.9208 | 0.8999 | 0.367 | 21.62 | 303 | −24 |

7450 | 264 | 14 | −1000 | 0.9285 | 0.8985 | 0.237 | 23.03 | 883 | −136 |

1750 | 255 | 25 | −2000 | 0.8855 | 0.8916 | 0.367 | 21.87 | 321 | −28 |

7600 | 266 | 14 | −2000 | 0.8876 | 0.8894 | 0.224 | 21.74 | 851 | −128 |

1950 | 255 | 24 | −2600 | 0.8674 | 0.8866 | 0.340 | 20.89 | 332 | −30 |

N_{c}
. | T (K)
. | ΔT (K)
. | p (bars)
. | ρ_{w} (g/cm^{3})
. | ρ_{ice} (g/cm^{3})
. | |Δμ| (kJ/mol)
. | γ (mJ/m^{2})
. | ΔG_{c} (kJ/mol)
. | log_{10} (J[m^{−3} s^{−1}])
. |
---|---|---|---|---|---|---|---|---|---|

1650 | 255 | 23 | −1000 | 0.9208 | 0.8999 | 0.367 | 21.62 | 303 | −24 |

7450 | 264 | 14 | −1000 | 0.9285 | 0.8985 | 0.237 | 23.03 | 883 | −136 |

1750 | 255 | 25 | −2000 | 0.8855 | 0.8916 | 0.367 | 21.87 | 321 | −28 |

7600 | 266 | 14 | −2000 | 0.8876 | 0.8894 | 0.224 | 21.74 | 851 | −128 |

1950 | 255 | 24 | −2600 | 0.8674 | 0.8866 | 0.340 | 20.89 | 332 | −30 |

To further understand this behavior, we connect our results with those from previous studies where nucleation had been studied for the TIP4P/Ice model at different pressures including negative, moderate, and high pressure states.^{7,23} In Fig. 1(a), we show the critical nucleus size as a function of supercooling for several isobars. We provide results at moderate supercoolings at −2600, −2000, and −1000 bars. For these same isobars and for the 1 bar isobar, we show the values reported in Ref. 23. For the 1 bar isobar, we also show the values given in Ref. 7, which also provides the values at the 2000 bars isobar. As can be seen, only the points corresponding to the 2000 bars isobar^{7} exhibit a different trend. The isobars at −2600, −2000, and −1000 bars from this work as well as from Ref. 23, and the 1 bar isobar from both Refs. 7 and 23 follow approximately the same curve. Notice that even a point at 450 bars reported in Ref. 23 was included being in agreement with this group of isobars. In fact, as shown in Fig. 1(b), pressure hardly affects the nucleation free energy barrier as a function of supercooling from −2600 to 450 bars.

Our results from Fig. 1 suggest that a similar nucleation behavior as a function of supercooling may take place from −2600 to 450 bars. That is a strikingly different behavior to the one observed when increasing pressure to 2000 bars. Thus, we propose universal empirical expressions for the variation of different homogeneous ice nucleation properties with the supercooling independently of the pressure as long as it lies within this regime. Nevertheless, we first need to confirm that what was observed for *N*_{c} and Δ*G*_{c} also applies to *J*. In Fig. 2, we show again (a) *N*_{c} and (b) Δ*G*_{c} as well as (c) *γ* and (d) log_{10}*J*. This time, for each magnitude, we include a common fit to data from moderately positive to deeply negative pressure (including our own data and those from Refs. 7 and 23) along a separate fit at high pressure.^{7} In Fig. 2(c), we show *γ* which exhibits higher variance. Finally, in Fig. 2(d), we show how very different pressures (from largely negative to moderately positive) lead to approximately the same nucleation rate *J* as a function of supercooling, Δ*T* = *T*_{m} − *T*. The values of *T*_{m} are given in Table II. Hence, we can use the respective common fit as universal empirical expressions to describe the change with supercooling within this broad range of pressures.

p_{m} (bars)
. | T_{m} (K)
. | γ_{m} (mJ/m^{2})
. |
---|---|---|

−2600 | 279.0 | 27.1(1.5) |

−2000 | 280.0 | 26.5(1.5) |

−1000 | 278.0 | 25.6(1.5) |

1 | 270.0 | 27.2 (0.8) |

1000 | 260.0 | 29.0(1.5) |

2000 | 246.5 | 37.2(1.5) |

p_{m} (bars)
. | T_{m} (K)
. | γ_{m} (mJ/m^{2})
. |
---|---|---|

−2600 | 279.0 | 27.1(1.5) |

−2000 | 280.0 | 26.5(1.5) |

−1000 | 278.0 | 25.6(1.5) |

1 | 270.0 | 27.2 (0.8) |

1000 | 260.0 | 29.0(1.5) |

2000 | 246.5 | 37.2(1.5) |

For *N*_{c}, we obtain

and for Δ*G*_{c} (in kJ/mol), we obtain

where *T*_{0} equals 1 K for correctness of units. The uncertainty in the exponents of both Eqs. (7) and (8) after combining all data in the pressure regime below 500 bars is roughly 5%. This is consistent with the significant uncertainty in the critical nucleus size and with the exponents obtained from power law fits to each isobar alone, which are shown in the caption of Fig. 1 only with the case of *N*_{c} at 1 bar differing slightly.

For *γ* in mJ/m^{2}, we obtain

and, finally, for *J* in m^{−3} s^{−1}, one should use Eq. (4) along with Eq. (8) (after converting into in *k*_{B}*T* units), and 10^{36} m^{−3} s^{−1} as the prefactor $(\rho w|\Delta \mu |/(6\pi kBTNc)f+)$.

The results shown in Fig. 2 have interesting consequences. First, taking into account that *N*_{c} and *γ* [panels (a) and (c), respectively] are roughly independent of *p* when it goes from largely negative to moderately positive pressures, the isobaric Tolman length, which determines the change in *γ* with the inverse of the radius of curvature of the cluster along an isobar,^{8,61,62} is roughly constant too and equal to 0.24(5) nm, where the parentheses indicates uncertainty in the last digit. This result is in agreement with previous work.^{56} Second, in panel (d), one can see that from strongly negative to moderately positive pressure we obtain the same nucleation rate with respect to the supercooling, which means that the homogeneous nucleation line (HNL) should be at a constant distance to the melting line in this regime as predicted recently for this water model^{23} and for the mW model^{63} in Ref. 24. In Fig. 3, we show the estimates for the model^{7,23} assuming that the HNL corresponds to an iso-nucleation rate of log_{10}*J*/(m^{−3} s^{−1}) = 15, and we compare it to the experimental HNL.^{4} In addition, the coexistence lines of the model^{23} and the experimental one^{22} are presented showing how the distance between the coexistence line and the HNL is roughly constant until pressure increases enough such that the required supercooling to reach log_{10}*J* = 15 becomes larger. However, even though this result might be useful, a physical explanation is still missing. In order to answer this question, we look at the pressure-induced deceleration of ice nucleation. In 2016, Espinosa *et al.*^{7} showed that the origin of this phenomenon arises from the increase with pressure of the interfacial free energy both at coexistence *γ*_{m} and for nucleation (*γ* at a given supercooling Δ*T*), while the difference in chemical potential Δ*μ* does not change so much with Δ*T*. Thus, one needs a larger Δ*T* to obtain the same *J* at high pressure. In this work, we observe approximately the same *J* as a function of Δ*T* from strongly negative to moderately positive pressure.

Since we obtain roughly the same *γ* as a function of Δ*T* at different negative pressures, we also expect *γ*_{m} to barely change with *p*. The term *γ*_{m} refers to a planar interface between ice and water at certain conditions along the coexistence line, whereas the term *γ* refers to a curved interface between a critical nucleus of ice and water at a certain supercooling Δ*T* along an isobar. In both cases, thermodynamic equilibrium holds. However, when the interface is planar, then the pressure is equal in both phases, while in a spherical interface, the pressure changes between phases following the Young–Laplace equation. Then, we compute *γ*_{m} for several points. In addition to the negative pressure isobars, we compute two points at 1000 and 2000 bars. We study only the basal plane as we do not expect severe anisotropy (as much as ∼10%) with the prismatic ones.^{41,64–66} For example, the change along coexistence of this anisotropy has been studied in detail for a Lennard-Jones system by Laird *et al.*^{67} There, it was shown that the relative difference in the interfacial free energy among planes is almost constant with pressure. Furthermore, the anisotropy at 1 bar is comparable between water models^{41} and also in salty water where salt plays a similar role to pressure in terms of disruption of the hydrogen bond network.^{68} The results of *γ*_{m} are presented in Table II and in Fig. 4. As shown, *γ*_{m} barely changes along the coexistence line when *p* varies from strongly negative to moderately positive. Interestingly, *γ*_{m} displays a shallow minimum. Thus, as long as Δ*μ* does not change significantly with Δ*T* at negative *p*, one can explain why in Fig. 2, *N*_{c}, Δ*G*_{c}, *γ*, and *J* seem to be independent of *p* against the supercooling when *p* is negative or moderate. In order to confirm this, we evaluate the effect of *p* on Δ*μ* as a function of supercooling Δ*T* by comparing with the value at 1 bar. To do so, we compute (Δ*μ*_{p} − Δ*μ*_{1})/Δ*μ*_{1} for the different isobars *p* = −2600, −2000, −1000, 1, and 2000 bars as a function of Δ*T* (for 1 and 2000 bars, we use the data from Ref. 7). As can be seen in Fig. 5(a), the 2000 bars isobar is very similar to the −1000 bars one in terms of Δ*μ* with respect to Δ*μ*_{1}, and the −2600 bars is the one that deviates the most with up to 18%. This deviation is, however, compensated by *γ*, which is rather dispersed and, in the end, *N*_{c}, Δ*G*_{c}, and *J* are very well-described by universal empirical expressions. Moreover, in Fig. 5(b), we show Δ*G* obtained as *N*_{c}|Δ*μ*|/2 by setting *N*_{c} to the common fit of Eq. (7) and changing Δ*μ* to that of the different isobars. As can be seen, from strongly negative to moderately positive pressure, the change in Δ*μ* does not significantly affect the free energy barrier for isobars between −2600 and 450 bars. Thus, we confirm that the universality in nucleation properties presented in Figs. 1 and 2 is the consequence of the small variation with *p* of the difference in chemical potential Δ*μ* as well as in the interfacial free energy both at coexistence *γ*_{m} and for the nucleation *γ* at a given Δ*T*.

Regarding the nucleation rate *J*, one may also consider the effect of water kinetics in order to confirm previous conclusions. It is true that the thermodynamic contribution of Δ*G*_{c} has a larger effect than the kinetic prefactor, as can be seen in Eq. (4), since the first is in the exponential, while the latter affects the nucleation rate linearly. In fact, the diffusion coefficient *D* of supercooled water can change by orders of magnitude along isotherms.^{69} In Fig. 6, we show *D* as a function of Δ*T* along different isobars from −2600 up to 1400 bars. As can be seen between −1000 and 1400 bars, the change is little for supercoolings below 30 K. Only when the supercooling is larger than 30 K or the pressure is strongly negative, one starts to appreciate the difference. Furthermore, these differences are only up to one order of magnitude, which is not much in terms of *D* considering how different the thermodynamic states are in pressure and temperature. In fact, this variation is not significant in terms of the nucleation rate because *J* can vary tens of orders of magnitude within a range of a few kelvin of supercooling. Therefore, this result is compatible with the universality found in *J*.

### B. Interfacial free energy and melting line of the ice Ih–water interface

We now understand the small variability with pressure of the nucleation properties as a function of supercooling at negative and moderate pressure. In order to understand why *γ*_{m} displays a shallow minimum, we use the thermodynamic formalism of Gibbs for interfaces.^{51,70} The interfacial Gibbs–Duhem relation is given by

where Γ = *N*_{γ}/*A* is the surface excess density, also called adsorption, and *η*_{γ} = *S*_{γ}/*A* is the excess contribution to the entropy. Since the location of the dividing surface is arbitrary, excess functions also depend on this choice with the exception of *γ*_{m}. For a planar interface, *γ*_{m} does not change with the location of the dividing surface unlike in the case of curved interfaces, where *γ* does change with its location.^{51,53,55} The choice that most simplifies the thermodynamic treatment in our case is the equimolar dividing surface, usually denoted as the Gibbs dividing surface, where the excess components *N*_{γ} is zero, and so is Γ (see the Appendix for a general dividing surface treatment). Hence, we can write

where the superscript *e* denotes the equimolar dividing surface. Equation (11) provides us with the temperature dependence of the interfacial free energy. It is crucial to note that this derivative must be taken along the coexistence line so that *p* is not constant. In fact, we can change Eq. (11) to describe the change of *γ*_{m} with pressure along the melting line *p*_{m} as

In our case, Eq. (12) is more convenient due to the reentrant behavior of the melting curve, i.e., for each *T*_{m}, one has two values of *p*_{m}, whereas for each *p*_{m}, there is only one value of *T*_{m} (see the solid red curve in Fig. 3). From Eq. (12), one can see that the change in *γ*_{m} with *p*_{m} is determined by the slope of the melting line and the value of the excess entropy per area at the equimolar dividing surface, $\eta \gamma e$. This means that if there is a reentrant behavior for the melting point, there must also be a reentrant behavior for *γ*_{m} as a function of pressure exactly at the same *p*_{m} because $\eta \gamma e$ must be finite. In fact, Bianco *et al.*^{23} reported a reentrant behavior in the ice Ih-liquid coexistence line of TIP4P/Ice, whose turning point occurred at −2000 bars.

Next, we want to confirm that the maximum in the melting line *T*_{m}(*p*_{m}) is consistent with the minimum in *γ*_{m}(*p*_{m}) that we have obtained from the mold integration technique. Thus, we fit the data for *γ*_{m}(*p*_{m}) from mold integration with a quadratic fit with the constraint of having the vertex at the same *p*_{m} (−2000 bars) as the quadratic fit for *T*_{m}(*p*_{m}). The latter, $Tm(p)=aTmp2+bTmp+cTm$ has the parameters $cTm=$ 271 K, $bTm=\u22128.5\u22c510\u22123$ K/bar, and $aTm=\u22122\u22c510\u22126$ K/bar^{2}. In this way, we assume that $\eta \gamma e$ is constant. The result is shown in Fig. 7. In the left panel, we show the melting line with points from the direct coexistence simulations of Ref. 23 and the quadratic fit. An inset is included zooming in the region of the maximum of the melting line. On the right panel, we show the points of *γ*_{m} from mold integration from this work and Ref. 41 as well as the quadratic fit. As can be seen in the right panel, the fit is fairly good even though we impose constant $\eta \gamma e$ and quadratic fits with the constraint of having the vertex at the same *p*. Therefore, assuming that $\eta \gamma e$ is constant seems to be a reasonable approximation.

At this level of approximation, $\eta \gamma e$ is found to be 0.32 mJ/m^{2} K. Note that $\eta \gamma e>0$ as expected from Eq. (12). For instance, from 1 to 2000 bars, *T*_{m} decreases from 270 to 246.5 K, and *γ*_{m} increases from 27.2 to 37.2 mJ/m^{2}. Therefore, *dγ*_{m}/*dp*_{m} > 0 and *dT*_{m}/*dp*_{m} < 0, which mean that $\eta \gamma e$ is positive. On the other side of the vertex, from −2600 to −2000 bars, *T*_{m} increases from 279 to 280 K, while *γ*_{m} decreases from 27.1 to 26.5 mJ/m^{2}. Hence, *dγ*_{m}/*dp*_{m} < 0 and *dT*_{m}/*dp*_{m} > 0 so that the same sign in $\eta \gamma e$ holds. Note that Eq. (11) is only valid for planar interfaces along the melting line. If one tries to apply this equation away from of this line as was done in previous studies,^{7,14,71} one should probably incorporate terms that account for the change in *γ* due to curvature. Note also that the empirical relation proposed by Turnbull, which states that *γ*_{m} is proportional to the change in melting enthalpy Δ*H*_{m}, does not describe *γ*_{m} well at high pressure. From 1 to 2000 bars, Δ*H*_{m} decreases from 1.44 to almost 1 kcal/mol in experiments^{72} and from 1.29 to ∼1 kcal/mol^{73} for the TIP4P/Ice model. Thus, the Turnbull relation predicts a decreasing *γ*_{m}, which is not supported by our direct calculations via the mold integration technique.

As can be seen in Fig. 7, the knowledge of the equilibrium melting curve and the assumption of a constant value for the interfacial excess entropy are sufficient to understand the complex variation of *γ*_{m} along the melting line. The fact that the melting line does not depend on the exposed plane or in an order parameter shows that the shallow minimum found in *γ* should be robust. Another relevant excess variable that depends on *γ*_{m}, *T*_{m}, and $\eta \gamma e$ is the excess internal energy $e\gamma e$,

The excess internal $e\gamma e$ is the difference in energy between the actual system having an interface and a virtual system where the two phases remain unchanged up to the dividing surface (the equimolar one in this case). As a result of Eq. (12), the following relation holds:

so that if $\eta \gamma e$ is constant, then $e\gamma e$ must be constant as well. If we approximate $\eta \gamma e$ as constant with the value of 0.32 mJ/m^{2} K, we find $e\gamma e=$ 115 mJ/m^{2}.

## IV. CONCLUSIONS

In conclusion, we perform seeding simulations to study ice nucleation at negative pressures. Such conditions can be relevant in porous media and water transport in plants, where supercooled water can be at negative pressure. By comparing with the previous results, we show that universal empirical expressions describe *N*_{c}, Δ*G*_{c}, *γ*, and *J* as a function of supercooling for isobars in the regime from strongly negative (−2600 bars) to moderately positive pressures (500 bars). Only when pressure is high (2000 bars) do these relations break down. In the regime where pressure hardly plays any role, the isobaric Tolman length is predicted to be positive and roughly constant with a value of 0.24 nm. In addition, our results suggest that the homogeneous nucleation line should be parallel to the coexistence line when pressure is below ∼500 bars (while at higher pressure they are not). We explain this result by inspecting how the interfacial free energy at coexistence changes with pressure. We evaluate the interfacial free energy at coexistence at different states from strongly negative to high pressure by means of the mold integration technique. We show that the interfacial free energy at coexistence barely changes with pressure as long as the system is below 500 bars. In fact, a shallow minimum is reported at negative pressure, suggesting that the minimum interfacial free energy between ice Ih and water is around 26 ± 1 mJ/m^{2} for the basal plane expanding for a broad range of pressure centered around −2000 bars. Then, we use the Gibbsian formalism to explain that this minimum in the interfacial free energy is connected to a maximum in the melting temperature as a function of pressure. In particular, we show that the change in the interfacial free energy with pressure is proportional to the excess entropy and to the slope of the melting line. Thus, the reentrance in the interfacial free energy occurs because of the reentrance in the melting line, which happens due to the crossover in density between ice and water. Finally, we estimate the excess entropy and the excess energy of the ice Ih–water interface. We suggest that a constant value of 0.32 mJ/m^{2} K and 115 mJ/m^{2} is enough to provide a good description of the thermodynamics of the ice Ih–water interface within the studied pressure range.

## ACKNOWLEDGMENTS

The authors thank Eduardo Sanz and Salvatore Romano for fruitful discussions. P.M.d.H. acknowledges support from the SFB TACO (Project No. F81-N) funded by the Austrian Science Fund. J.R.E. acknowledges funding from the Oppenheimer Fellowship, the Roger Ekins Fellowship from Emmanuel College, and a Ramon y Cajal Fellowship (Grant No. RYC2021-030937-I). C.V. acknowledges support from Project No. PID2019-105898GB-C21 of the Ministerio de Educacion y Cultura. This work was performed using resources provided by the Spanish Supercomputing Network (RES), the Vienna Scientific Cluster (VSC), and the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (http://www.hpc.cam.ac.uk) funded by EPSRC Tier-2 Capital Grant No. EP/P020259/1.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**P. Montero de Hijes**: Conceptualization (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **J. R. Espinosa**: Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **C. Vega**: Conceptualization (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **C. Dellago**: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: INTERFACIAL FREE ENERGY ALONG THE MELTING LINE FOR A GENERAL DIVIDING SURFACE

In this work, we used the equimolar dividing surface for simplicity. However, Eqs. (11) and (12) can be generalized for any choice of the dividing surface. To do so, it is necessary to involve not only the interfacial Gibbs–Duhem relation (10) but also the ice and liquid Gibbs–Duhem relations. Respectively, these are

where *v* is the volume per molecule (the inverse of the number density) and *s* is the entropy per molecule. Since the phase equilibrium holds, *dμ*_{m}, *dp*_{m}, and *dT*_{m} are common in all phases. Note that from Eqs. (A1) and (A2), one can obtain the Clausius–Clapeyron relation that explains the slope of the melting line,

By also including Eq. (10) in the relation, one can obtain the temperature and pressure dependence of the interfacial free energy without imposing a specific dividing surface. For the temperature, one obtains

whereas for the pressure, one finds