Minimum in the pressure dependence of the interfacial free energy between ice Ih and water

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.³ 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.⁷ 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.⁸ Nevertheless, both notably increase at high pressure.⁷ 

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²⁵ and also in water transpiration inside plants²⁶ where negative pressure occurs. Creating and maintaining negative pressure over a sample is non-trivial in experiments.²⁷ 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,³² centrifugation,³³ 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.⁷ 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,³⁵ 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⁴¹ to measure the interfacial free energy at coexistence.

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,⁴⁴ both with a relaxation time of 0.5 ps. Electrostatic interactions are accounted for via the particle-mesh-Ewald summation algorithm⁴⁵ 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,⁵⁸ $q̄6(T,p)$⁠, is used in combination with the mislabeling criterion⁹ to identify ice-like and water-like molecules. We obtain $q̄6(T,p)$ for each molecule. The molecules with $q̄6(T,p)$ above a certain threshold $q̄6,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̄6,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,⁵⁹ 

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,⁴¹ 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 $γ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 $γ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.

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.

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⁷ 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₁₀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.⁷ 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.

For N_c, we obtain

and for ΔG_c (in kJ/mol), we obtain

where T₀ 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², we obtain

and, finally, for J in m⁻³ s⁻¹, one should use Eq. (4) along with Eq. (8) (after converting into in k_BT units), and 10³⁶ m⁻³ s⁻¹ as the prefactor $(ρw|Δμ|/(6π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.⁵⁶ 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²³ and for the mW model⁶³ 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₁₀J/(m⁻³ s⁻¹) = 15, and we compare it to the experimental HNL.⁴ In addition, the coexistence lines of the model²³ and the experimental one²² 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₁₀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.⁷ 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.⁶⁷ 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⁴¹ and also in salty water where salt plays a similar role to pressure in terms of disruption of the hydrogen bond network.⁶⁸ 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 − Δμ₁)/Δμ₁ 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 Δμ₁, 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.⁶⁹ 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.

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, $ηγ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 $ηγe$ must be finite. In fact, Bianco et al.²³ 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=−8.5⋅10−3$ K/bar, and $aTm=−2⋅10−6$ K/bar². In this way, we assume that $ηγ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 $ηγe$ and quadratic fits with the constraint of having the vertex at the same p. Therefore, assuming that $ηγe$ is constant seems to be a reasonable approximation.

At this level of approximation, $ηγe$ is found to be 0.32 mJ/m² K. Note that $ηγ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². Therefore, dγ_m/dp_m > 0 and dT_m/dp_m < 0, which mean that $ηγ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². Hence, dγ_m/dp_m < 0 and dT_m/dp_m > 0 so that the same sign in $ηγ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⁷² and from 1.29 to ∼1 kcal/mol⁷³ 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 $ηγe$ is the excess internal energy $eγe$⁠,

The excess internal $eγ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 $ηγe$ is constant, then $eγe$ must be constant as well. If we approximate $ηγe$ as constant with the value of 0.32 mJ/m² K, we find $eγe=$ 115 mJ/m².

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² 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² K and 115 mJ/m² is enough to provide a good description of the thermodynamics of the ice Ih–water interface within the studied pressure range.

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.

The authors have no conflicts to disclose.

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

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

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

As can be seen, at the equimolar dividing surface where Γ = 0, one recovers Eqs. (11) and (12). These expressions are relevant when nucleation data are extrapolated to coexistence because the relevant dividing surface in nucleation is usually the surface of tension for which Γ ≠ 0.