Free energy of critical droplets—from the binodal to the spinodal

Since nucleation is an activated process, the pressure of a vapor can be increased far above the saturation pressure, and the temperature of a liquid can be decreased well below the melting temperature before condensation or crystallization occurs. In thermodynamic theory, there is a limit to the degree of metastability that can be achieved called the spinodal (Fig. 1), beyond which the system becomes intrinsically unstable and phase transition is a spontaneous process. For a pure fluid, the thermodynamic spinodal can be defined as the state of diverging isothermal compressibility.¹¹ Upon approaching the spinodal, the nucleation barrier and the critical size are expected to go to zero, and the nucleation time becomes extremely short.¹ Therefore, these states are mostly studied indirectly.^12–14 

150 years after the pioneering work of Gibbs,⁶ nucleation phenomena keep fascinating scientists. The theoretical developments struggle to keep up with the increasingly accurate results from simulations and experiments that measure the nucleation rates of different substances and types of phase transitions.^15–25 Arguably, the main challenge of nucleation theory is to accurately estimate the work of formation for the critical embryo. Since the nucleation rate depends exponentially on this quantity, a small error will cause the predicted nucleation rate to deviate by orders of magnitude from experimental results or simulations.¹ 

For pure fluids, CNT is qualitatively correct.²⁶ However, it does not account for the expected vanishing of the nucleation barrier at the spinodal, and the predicted rates show systematic deviations from experiments, with errors reaching 20 orders of magnitude for the simple case of argon condensation.²⁷ CNT assumes that the incipient embryo of the new phase is spherical and that its properties are identical to those of the bulk macroscopic stable phase, including the value of the surface tension of the planar vapor–liquid surface. This is referred to as the capillarity approximation, and it has been the main suspect in the problems of CNT.⁵ 

To overcome the limitations of CNT, several routes have been pursued to provide more accurate estimates of the work of formation. Density gradient theory (DGT) assumes that thermodynamic properties depend on density gradients in addition to the local thermodynamic state. It obtains density profiles that minimize a constrained functional for the total Helmholtz energy.²⁸ Its computational simplicity and robustness make DGT attractive. However, it is only the first order approximation to the full classical density functional theory (DFT).²⁹ The accuracy and predictive ability of DGT has, therefore, been questioned, in particular at the highest metastabilities.³⁰ Using full DFT has been touted as a promising alternative.^31,32 In previous studies, DFT has been used with great success to estimate a range of important interfacial properties, such as surface tensions^33–36 and adsorption isotherms.³⁷ The performance of DFT is limited by a lack of knowledge about the exact functional to extremize³⁴ and knowledge about the real intermolecular potential, even for simple substances like argon³⁸ or water.³⁹ 

For many applications, it is impractical and too time-demanding to use DGT or full DFT to calculate the work of formation.^5,40 This motivates the development of simple corrections to CNT that maintain the inherent simplicity of the theory. Among the most successful are the curvature corrections of the planar surface tension. The second order curvature expansion of the surface tension is called the Helfrich expansion, where the coefficients of the first- and second-order corrections are given by the Tolman length and rigidity constants, respectively. The challenge of using this approach has been to evaluate the magnitude of the corrections. There have even been controversies about the sign of the leading-order correction, the Tolman length.^41–48 

Recently, it was shown how the Tolman length and rigidity constants can be calculated using DGT^49–52 or DFT.^53,54 Incorporating the curvature corrections into nucleation theory has been shown to significantly improve agreement with experimental results for water.⁴⁹ It has also solved some of the inconsistencies that arise when CNT is extended to multi-component systems.⁵² The ultimate test of the accuracy of this curvature-corrected CNT (c-CNT) should be performed using a model system where all the relevant ingredients and properties are as much as possible under control. In this context, the Lennard-Jones (LJ) fluid stands out because it is a simple potential that can be implemented in simulations and has been extensively studied in the literature.⁵⁵ 

In this work, we have performed a comprehensive study of droplets in the Lennard-Jones fluid truncated and shifted at 2.5σ, combining umbrella sampling Monte Carlo simulations, DGT, and DFT to test the accuracy of the curvature-corrected CNT. The focus has been to evaluate the work of formation of liquid droplets for a wide range of conditions from the binodal to the spinodal and at different temperatures. We show that DFT and DGT yield similar results that are in excellent agreement with the molecular simulation results. Remarkably, c-CNT reproduces very accurately the results from simulations, DGT, and DFT for all temperatures and moderate degrees of metastability.

Since the c-CNT is based on the Helfrich expansion, which is a Taylor expansion around the planar surface, it does not give a vanishing nucleation barrier at the spinodal, unless by chance. This hampers its precision at high metastabilities and for sufficiently small embryos. We mend this shortcoming by formulating a scaling expression for the nucleation barrier that incorporates both the coefficients of the Helfrich expansion and the requirement of vanishing nucleation barriers at the spinodals as corrections to CNT. We demonstrate that the modified theory gives an excellent representation of the work of formation that differs from the predictions from DGT by less than one k_BT for the entire metastability range and all temperatures examined.

For a given metastable phase (T, μ), Eqs. (2)–(4) constitute three equations for four properties. The closure of this system of equations is obtained by specifying a relation between γ and R. The closures considered in this work are detailed in Secs. II B–II F.

To calculate ΔP in Eq. (2) for condensation from a metastable vapor, one needs an equation of state (EoS) that is accurate for the stable liquid and the metastable vapor. For the Lennard-Jones fluid truncated and shifted at 2.5σ, LJTS(2.5σ), the PeTS equation of state by Heier et al.⁵⁹ satisfies this criterion. (Note that an EoS of similar accuracy was also developed for the even-shorter-ranged Lennard-Jones spline potential,⁶⁰ but for this potential simulations of droplets are missing.) We have verified the accuracy of PeTS in predicting the binodal by performing Gibbs ensemble Monte Carlo simulations^61,62 for the three temperatures studied here (see the supplementary material).

The grand canonical functional is stationary for a system at a fixed temperature T, volume, and chemical potential μ. The necessary conditions for a stationary point of Eq. (11) are provided by the Euler–Lagrange equations, $δΩρr/δρr=0$⁠. To solve the Euler–Lagrange equations for DGT, we used the approach described in Ref. 51.

The profiles from DFT can be obtained by solving the Euler–Lagrange equations, $δAρr/δρr=μ$⁠. This can be efficiently accomplished by using the convolution theorem to solve the convolution integrals in Fourier space with cosine and sine transforms, for both the planar and spherical geometries.^65,66 In addition, convergence of the overall equation system can be accelerated by use of Anderson mixing.⁶⁷ 

Sauer and Gross³⁴ regressed the φ parameter to surface tension data of n-alkanes and found that φ = 1.3862 gave good results for most fluids. For the short ranged LJTS(2.5σ) fluid, Heier et al. had to use a smaller value of φ = 1.21 to achieve agreement with the simulation results for the surface tension;⁶⁸ this is also the value used in this work.

For both DFT and DGT, the bulk behavior is governed by the underlying EoS, i.e., PeTS. Since the work of formation of a droplet (bubble) in DGT will be zero at the vapor (liquid) spinodal of the EoS, the location of this spinodal may severely affect the prediction of free energies of small clusters. Although non-local DFT yields more realistic density profiles than DGT,⁵⁵ they have been found to yield similar free energy predictions⁵⁴ in the planar limit. In this work, we will also compare them for small droplets.

Of interest to the current work is the behavior of mean field theories, such as DGT and DFT, close to the spinodal. Wilemski and Li³⁰ studied in detail the scaling behavior of density gradient theory and mean field density functional theory. Both of these theories always predict that the work of formation vanishes at the spinodal, whereas the excess number of particles [Eq. (5)] of the critical cluster diverges to infinity. The mean field theories predict that the equimolar radius R_e and the equimolar surface tension γ_e both diverge to infinity upon approaching the spinodal. For the surface of tension, however, γ → 0, R → 0 at the spinodal, whereas ΔP goes to a finite value.

Since the pioneering work of McGraw and Laaksonen on scaling relations,²⁶ there is increasing evidence suggesting that the work of formation of critical clusters normalized by the CNT prediction, f ≡ W/W_CNT, may be a universal function of the degree of metastability, depending only weakly on the temperature and on the details of the intermolecular potential.^69–73 Accordingly, many expressions for this scaling function f have been proposed in the literature,^{69,72,74–79} where the degree of metastability is typically expressed either in terms of the difference in chemical potential Δμ or the Laplace pressure ΔP. This scaling function should fulfill at least two requirements. First, CNT is expected to be valid at the binodal when the size of the critical cluster goes to infinity, such that f(0) = 1. Second, the work of formation of the critical cluster should be zero at the spinodal, such that f(ΔP_sp) = 0, where ΔP_sp represents the difference in bulk pressures between the stable and metastable phases at the spinodal.

A very similar parabolic scaling function was recently proposed by Kashchiev⁷⁹ based on the nucleation theorem and using the assumption that the equimolar critical cluster radius only differs from the CNT prediction by the Tolman length. The main difference is that the second order term is obtained by imposing that the free energy of the critical cluster must vanish at the vapor spinodal, instead of using the value of the spherical rigidity. In fact, Kaschiev’s expression can be recovered assuming that the scaling function is parabolic and using f(0) = 1, f(ΔP_sp)= 0, and the value of Tolman length δ as conditions to evaluate the coefficients.

The first two factors in the numerator of Eq. (22) warrant that the work of formation vanishes at both spinodals. The coefficients of the parabolic third factor are obtained from the slope at the binodal given by the Tolman length and the location of the maximum of the work of formation at ΔP_m.

In Sec. IV, we will show that Eq. (22) accurately reproduces results from simulations, DGT, and DFT. The proposed scaling function in Eq. (22) is a modification of the Helfrich expansion so that the surface tension goes to zero at the right location of the spinodals while maintaining the correct behavior in the binodal limit described by the Tolman length and the spherical rigidity.

In order to test the accuracy of c-CNT and of the proposed scaling relation, we have performed umbrella sampling Monte Carlo simulations to calculate the free energy of critical droplets of the LJTS(2.5σ) fluid at two different temperatures: T = 0.625 and T = 0.900. Monte Carlo simulations in the NPT ensemble were implemented in a cubic box with periodic boundary conditions, using a self-developed code based on the algorithms described in Refs. 80 and 81. Simulations were performed with typically N = 4000 particles for the largest supersaturations and N = 10 000 for small supersaturations. In the MC simulations, each trial move consisted either of an attempted particle displacement or a trial volume change, chosen at random. The maximum displacement for both types of MC moves was adjusted to give a 50% acceptance ratio. A MC step is defined as N + 1 trial moves.

Liquid clusters were identified using the ten Wolde–Frenkel cluster criterion,⁸¹ where a molecule is classified as liquid-like if it has at least five neighbors within a distance smaller than the Stillinger radius r_c, corresponding to the first minimum in the radial distribution function of the liquid at the desired temperature.

For T = 0.900, simulations were performed in the pressure range of 0.034 25 ≤ P ≤ 0.039 and for T = 0.625 in the range of 0.0028 ≤ P ≤ 0.078. For pressures smaller than that, the critical cluster size was too large to be simulated. For pressures larger than those values, the system was too close to the spinodal, and many clusters were formed spontaneously, hampering a reliable reconstruction of the barrier. Therefore, the range of nucleation barriers accessible through our umbrella sampling simulations is in the range of 20–150 k_BT. Details from the simulations, such as the work of formation and the number of particles, are included in the supplementary material.

Simulation results for a third intermediate temperature T = 0.741 were also included in our comparison. For that temperature, we took advantage of the results of ten Wolde and Frenkel,⁸¹ who also used umbrella sampling to calculate the free energy at P = 0.007 83, and then used thermodynamic integration to obtain the nucleation barrier at other pressures.

We next discuss the results. The numerical values of the simulation results are listed in the supplementary material. All properties are reported in reduced Lennard-Jones units. Conversion formulas for real fluids are listed in Table I of Ref. 83.

Figures 2–4 compare the number of particles and the work of formation predicted from all theories considered in this work against simulations for the temperatures T = 0.625, T = 0.741, and T = 0.900. The x axes start from the binodal and end at the vapor spinodal for each temperature.

All theories, namely, DGT, DFT, CNT, and c-CNT (labeled “Helfrich expansion” in Figs. 2–4), and the scaling function, reproduce the excess particle number in the critical cluster ΔN obtained in the simulations, mostly within the simulation uncertainty. Although the excess number of particles is closely related to the work of formation through Eq. (5), in the simulations, it is calculated directly from the number of molecules in the critical cluster identified by the ten Wolde-Frenkel cluster criterion. Therefore, the fact that a theory can reproduce both the work of formation and the excess number of particles represents a stringent validation as well as an excellent consistency check of the simulations. Mean field theories, such as DFT and DGT, predict a diverging excess number of particles upon approaching the spinodal.³⁰ This is not shown in Figs. 2–4 [panel (a)] since this divergence only occurs very close to the spinodal. However, note that none of the theories predict zero excess particle number at the spinodal; this is most evident from Fig. 3, panel (a).

For the work of formation, DGT and the proposed spinodal scaling relation are in closest agreement with the molecular simulation results. This can be seen by comparing the work of formation ratio W/W_CNT and the value of W − W_DGT [Figs. 2–4, panels (c) and (d)] to the simulations.

The excellent match of the simulations and DGT with a temperature-independent influence parameter is remarkable, considering that the uncertainty of the planar surface tension values for the LJTS(2.5σ) fluid is on the order of several percent (cf. the supplementary material). Intuitively, one might expect DFT to give results that are closer to simulations than DGT. After all, DGT is only the first approximation to DFT. It should be kept in mind, however, that the exact functional to minimize in DFT is unknown, and the functional presented in Sec. II D represents only an approximation. A fitting parameter is needed in both theories to reproduce surface tensions from molecular simulations.

Stephan et al.⁵⁵ showed that although the surface tension of the planar vapor-interface is accurately represented by both DFT and DGT, the theories give density and pressure tensor profiles that deviate from molecular simulation results. One theory did not show any clear advantage over the other (see Figs. 2–6 in Ref. 55). These discrepancies do not seem to influence the curvature dependence of the surface tension to a large degree, since both theories reproduce the work of formation from umbrella sampling simulations with impressive accuracy.

The c-CNT that uses the Helfrich expansion is also remarkably accurate when measuring the deviation in units of k_BT for the droplet sizes studied in this work [Figs. 2–4, panel (d)]. Still, the deviation from simulations is only sligthly outside the simulation error bars at the largest metastabilities, as seen in Fig. 2. Only the CNT approximation results in large deviations from the simulations.

With the exception of CNT, all theories predict that the surface tension increases with metastability at low, positive values of ΔP, reaches a maximum, and then decreases monotonically toward the spinodal. In practice, homogeneous nucleation only occurs for metastabilities where W < 100k_BT; for higher barriers, the nucleation rate is generally too low to be observable. For nearly all droplets for which W < 100k_BT, the theories predict that the surface tension is lower than the planar surface tension for the LJTS(2.5σ) fluid.

The theories that only use the properties of the surface tension associated with the binodal, namely, CNT and the Helfrich expansion, differ markedly at high supersaturations from the theories that incorporate the spinodal. CNT does not predict a spinodal at all, whereas the Helfrich expansion predicts spinodals at significantly higher metastabilities than the PeTS EoS. Both theories, therefore, predict a positive work of formation W at the PeTS spinodal, which is largest for the lowest temperatures [compare Figs. 2–4 (panel d) for the three temperatures]. For CNT, this overprediction persists to the metastabilities of the simulations, showing clearly that CNT overpredicts the free energy of small droplets and that this is most severe at low temperatures and/or high metastabilities.

Figures 2–4 [panel (c)] illustrate that DFT and DGT yield very similar shapes of the ratio $W/WCNT=(γ/γ0)3$⁠, where γ₀ was calculated from the corresponding theory. For the value of W [Figs. 2–4 panel (d)], DGT performs slightly better than DFT in reproducing the simulation results. One possible reason is the different values of γ₀ from the two theories (cf. Table I), as well as the different Helfrich coefficients $δ,k,k̄$⁠. That the curves W_CNT − W_DGT in panel (d) diverge to infinity at the binodal follows from Eqs. (8)–(10), since the Tolman length of DGT is negative. We also see that W_DFT − W_DGT diverges when approaching the binodal, where the divergence is to ∞ or −∞ depending on the sign of $γ0DFT−γ0DGT$⁠.

Although we have not compared to molecular simulation results for bubbles, we have included a comparison between the theories and methods used in this work for bubbles in the supplementary material. This comparison shows that the performance of the theories discussed above applies to the bubble branch as well. Remarkably, the proposed scaling function Eq. (22) reproduces almost exactly the predictions of DGT both for bubbles and droplets in the entire range of metastabilities. The supplementary material also includes versions of Figs. 2–4, where the quantities are plotted as functions of the supersaturation ratio P/P_sat, which is commonly used as a measure of metastability in nucleation theory.

We will next compare our estimates of the Tolman length to those of van Giessen and Blokhuis,⁸⁴ who determined the Tolman length δ of the LJTS(2.5σ) fluid at T = 0.900 from bulk pressures of large liquid drops via molecular dynamics simulations. Using the planar surface tension $γ0∞$ as a boundary condition in their analysis, they found δ = −0.100σ or δ = −0.0845σ depending on how they calculated the bulk pressures. This deviates somewhat from the value determined by DGT and DFT (cf. Table I). In fact, the values of the Tolman length we find, −0.067 36 from DGT and −0.058 28 from DFT, are closer to those estimated by Wilhelmsen et al.,⁵⁰ namely, δ = −0.074σ. Wilhelmsen et al. used another implementation of DGT with a less accurate equation of state.⁵⁰ This shows that the coefficients in the Helfrich expansion are sensitive toward the underlying thermodynamic description of the bulk fluid.

For T = 0.625, the simulation results for the largest droplets yield values of the surface tension that are higher than the planar value (Fig. 2), which is consistent with a negative Tolman length. To our knowledge, this is the first time this has been shown by umbrella sampling simulations. The regime with surface tensions that are higher than the surface tension of the planar surface was only reachable at the lowest temperatures. This is because the critical cluster size corresponding to the maximum and the required number of particles in the simulation volume increase quickly with temperature.

Interestingly, the value of the spherical rigidity $ks=2k+k̄$ appears to be approximately independent of temperature for T ≤ 0.9, with k_s ∼ −0.45 for DGT (cf. Table I). This is in line with the findings by Wilhelmsen et al.⁵⁰ How the W/W_CNT curve varies with temperatures is therefore strongly connected to the value of the Tolman length.

In this work, we have performed a careful analysis of the work of formation of critical clusters of the truncated and shifted Lennard-Jones fluid LJTS(2.5σ) at three different temperatures and a wide range of pressures from the binodal to the spinodal. We have performed Monte Carlo umbrella sampling simulations and used density gradient theory (DGT) and density functional theory (DFT) to study the work of formation and the size of critical clusters. The aim was to test the accuracy of classical nucleation theory (CNT) and curvature corrected CNT (c-CNT). c-CNT relies on the Helfrich expansion, which is a second order expansion of the surface tension in curvature around the planar limit, where the first order correction is the Tolman length and the second order correction is the rigidity constants.

We have found that both DGT and DFT reproduce the simulation results with remarkable accuracy. As expected, since CNT uses the capillarity approximation, it overpredicts the nucleation barrier at high supersaturations and fails to predict a vanishing nucleation barrier at the spinodal. Remarkably, c-CNT reproduces accurately the results from molecular simulations at all temperatures and for most of the experimentally accessible region. This provides a solid confirmation of the validity of c-CNT, using a model system like the LJTS(2.5σ), with a simple intermolecular potential, and for which the equation of state and relevant properties have been characterized accurately by simulations.

The Tolman lengths from both DGT and DFT were found to be negative. Moreover, the second-order correction of the surface tension, the spherical rigidity, was found to be nearly independent of temperature. At the lowest temperature investigated, the work of formation from umbrella sampling simulations revealed a region where the surface tension is higher than the value of the planar vapor–liquid surface. This strongly suggests that the Tolman length is negative, in line with the predictions from DGT and DFT.

The molecular simulations for T = 0.900 were shown to be in agreement with capillary wave theory, suggesting that the capillary wave contributions to the planar surface tension must be removed before using it to tune DGT/DFT when modeling droplets. This warrants further investigation of the size-dependence of the surface tension for other temperatures, e.g., T = 0.625 and T = 0.741, and for other fluids to investigate whether this finding is general.

The only regime where c-CNT deviates from simulations is for the smallest droplets at the highest metastabilities in the vicinity of the spinodal. In fact, c-CNT mispredicts the location of the spinodal, identified as the vanishing of the nucleation barrier. This is expected, since c-CNT is based on a curvature expansion around the planar limit, i.e., at the binodal.

In order to solve this limitation, we have formulated a new scaling function that, in addition to the leading order curvature corrections (i.e., Tolman length and rigidity constants), also uses the location of the vapor and liquid spinodals. The resulting fourth order polynomial reproduces with excellent accuracy the results from simulations, DGT and DFT, in the whole range of studied temperatures and degrees of metastability, from the binodal to the spinodal. It also reproduces, with remarkable accuracy, the DGT and DFT predictions for the bubble branch.

Our results suggest that both c-CNT and the new scaling function constitute accurate yet simple models to accurately predict nucleation rates and solve the limitations of CNT. Here, we have confirmed this for the work of formation of droplets. However, the new scaling function is expected to be equally accurate for bubbles and could also possibly be extended to crystallization. Therefore, the proper incorporation of curvature corrections and the location of the spinodals might open the door to quantitatively accurate predictions of nucleation rates for other substances and a wide range of problems of scientific and technological interest.

The supplementary material contains a comparison of the predictions of the different theories on the bubble branch; a comparison between simulations and theories for critical droplets for T = 0.900 when refitting the influence parameter of DGT to reproduce the planar surface tension from van Giessen and Blokhuis;84 a table with all the simulation results; a validation using Gibbs ensemble Monte Carlo simulations of the PeTS EoS for the three temperatures studied in this work; a comparison of the calculated surface tensions from DFT and DGT against the simulation results by Vrabec et al.;68 and an alternative representation of Figs. 2–4 in terms of the supersaturation.

This work received funding from the Norwegian Research Council, Project No. 328679. A.A., Ø.W., and M.H. acknowledge funding from the Research Council of Norway (RCN), the Center of Excellence Funding Scheme, Project No. 262644, PoreLab. D.R. acknowledges funding from the Spanish government through Grant Nos. PGC2018-098373-B-I00 and PID2021-126570NB-I00 (MINECO/FEDER, UE).

The authors have no conflicts to disclose.

Ailo Aasen: Conceptualization (equal); Data curation (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Øivind Wilhelmsen: Conceptualization (equal); Investigation (equal); Supervision (equal); Writing – review & editing (equal). Morten Hammer: Data curation (equal); Investigation (equal); Writing – review & editing (equal). David Reguera: Conceptualization (equal); Data curation (equal); Investigation (equal); Supervision (equal); Writing – review & editing (equal).

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