Arguably, the main challenge of nucleation theory is to accurately evaluate the work of formation of a critical embryo in the new phase, which governs the nucleation rate. In Classical Nucleation Theory (CNT), this work of formation is estimated using the capillarity approximation, which relies on the value of the planar surface tension. This approximation has been blamed for the large discrepancies between predictions from CNT and experiments. In this work, we present a study of the free energy of formation of critical clusters of the Lennard-Jones fluid truncated and shifted at 2.5*σ* using Monte Carlo simulations, density gradient theory, and density functional theory. We find that density gradient theory and density functional theory accurately reproduce molecular simulation results for critical droplet sizes and their free energies. The capillarity approximation grossly overestimates the free energy of small droplets. The incorporation of curvature corrections up to the second order with the Helfrich expansion greatly remedies this and performs very well for most of the experimentally accessible regions. However, it is imprecise for the smallest droplets and largest metastabilities since it does not account for a vanishing nucleation barrier at the spinodal. To remedy this, we propose a scaling function that uses all relevant ingredients without adding fitting parameters. The scaling function reproduces accurately the free energy of the formation of critical droplets for the entire metastability range and all temperatures examined and deviates from density gradient theory by less than one *k*_{B}*T*.

## I. INTRODUCTION

^{1–5}They are initiated by thermal fluctuations that lead to the formation of a nano-sized embryo of the new phase. In order to grow and for the phase transition to proceed, the size of the embryo must be sufficiently large to overcome the free energy cost of forming a surface. The thermodynamic principles of such processes were explained in the seminal work of Gibbs.

^{6}Together with the theory to describe the kinetics developed by Volmer and Weber,

^{7}Farkas,

^{8}Becker and Döring,

^{9}and Zeldovich,

^{10}among others, established the basis of Classical Nucleation Theory (CNT), which still today leads our current understanding of nucleation phenomena. According to CNT, the rate of formation of embryos in the new phase per unit volume (

*J*), called the nucleation rate, is described by the following formula:

*W*is the work of formation of the critical embryo,

*k*

_{B}is Boltzmann’s constant,

*T*is the temperature, and

*K*is a kinetic pre-factor.

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.^{11} Upon approaching the spinodal, the nucleation barrier and the critical size are expected to go to zero, and the nucleation time becomes extremely short.^{1} Therefore, these states are mostly studied indirectly.^{12–14}

150 years after the pioneering work of Gibbs,^{6} 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.^{1}

For pure fluids, CNT is qualitatively correct.^{26} 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.^{27} 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.^{5}

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.^{28} Its computational simplicity and robustness make DGT attractive. However, it is only the first order approximation to the full classical density functional theory (DFT).^{29} The accuracy and predictive ability of DGT has, therefore, been questioned, in particular at the highest metastabilities.^{30} 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.^{37} The performance of DFT is limited by a lack of knowledge about the exact functional to extremize^{34} and knowledge about the real intermolecular potential, even for simple substances like argon^{38} or water.^{39}

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.^{49} It has also solved some of the inconsistencies that arise when CNT is extended to multi-component systems.^{52} 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.^{55}

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*_{B}*T* for the entire metastability range and all temperatures examined.

## II. THEORY

*T*and chemical potential

*μ*. For a spherical critical cluster, the work of formation is given by the formally exact relation,

^{6,56}

*γ*is the surface tension of the cluster with the surface of tension as a dividing surface. The Laplace pressure is defined as

*T*,

*μ*). The calculation of the work of formation thus requires the evaluation of two main ingredients: the pressure difference Δ

*P*and the surface tension

*γ*.

*P*is the difference in

*bulk*thermodynamic properties. It can be related to the surface tension through the Laplace relation,

*R*is the radius of the critical cluster at the surface of tension.

^{51,57}In this work, the quantities

*R*,

*γ*(with no subscript) are always associated with the surface of tension.

^{58}

### A. Equation of state

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.*^{59} satisfies this criterion. (Note that an EoS of similar accuracy was also developed for the even-shorter-ranged Lennard-Jones spline potential,^{60} 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).

### B. The capillarity approximation and classical nucleation theory

*γ*is the same as that of the bulk liquid with a planar interface

*γ*

_{0},

*W*/

*W*

_{CNT}→ 1 in the planar limit (Δ

*P*→ 0), and this is often interpreted as the work of formation from CNT being exact in the limit of large droplets. At the same time,

*α*= (

*∂γ*/

*∂*Δ

*P*)

_{T}. Therefore, if the surface tension has a non-zero derivative with respect to Δ

*P*in the planar limit, the absolute error of CNT for the work of formation can be substantial for large droplets. More generally, the difference in work of formation from two theories will diverge upon approaching the binodal if the two theories predict differing values of

*γ*

_{0}or

*∂γ*/

*∂*(Δ

*P*) at the binodal (or both). This must be kept in mind when comparing theories close to the binodal. But obviously, the relative error (

*W*−

*W*

_{CNT})/

*W*

_{CNT}goes to zero at the binodal, as expected for the thermodynamic limit.

### C. Density gradient theory

*T*,

*μ*,

*V*) by the following functional of the density profile:

*κ*is a positive parameter, usually fitted to reproduce the planar surface tension, and

*a*

^{EoS}(

*ρ*) is the Helmholtz energy density (per volume) of a fluid with uniform density

*ρ*at the specified temperature. The function

*a*

^{EoS}(

*ρ*) = −

*P*

^{EoS}(

*ρ*) +

*μ*

^{EoS}(

*ρ*)

*ρ*is calculated from the PeTS EoS. We used the constant influence parameter

*κ*= 2.7334 (in reduced, LJ-units) suggested by Heier

*et al.*,

^{59}which was shown to reproduce planar surface tension simulations almost within their error bars.

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, $\delta \Omega \rho r/\delta \rho r=0$. To solve the Euler–Lagrange equations for DGT, we used the approach described in Ref. 51.

### D. Density functional theory

*et al.*

^{59}derived a density functional theory (DFT) for the LJTS(2.5

*σ*) fluid based on a weighted density average (WDA) approach.

^{34}To arrive at a formulation that is consistent with Barker–Henderson perturbation theory

^{63}as represented by the PeTS EoS, the functional used the White Bear fundamental measure theory (FMT) for the hard sphere contribution.

^{64}Density functional theory models the grand free energy at the state defined by (

*T*,

*μ*,

*V*) by the following functional of the density profile:

*σ*) fluid, the Helmholtz energy has contributions from the hard-sphere reference (FMT), dispersion (disp), and an ideal gas (id) term,

*ω*

^{disp}is a weight function of a sphere with diameter

*φd*

_{BH},

*d*

_{BH}is the temperature dependent Barker–Henderson diameter,

^{63}and

*φ*is a parameter that depends on the potential range. Using the WDA approach, the dispersion contribution is a function of the dispersion term from the EoS as follows:

The profiles from DFT can be obtained by solving the Euler–Lagrange equations, $\delta A\rho r/\delta \rho r=\mu $. 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.^{67}

Sauer and Gross^{34} 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;^{68} 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,^{55} they have been found to yield similar free energy predictions^{54} 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^{30} 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.

### E. The Helfrich expansion and c-CNT

*R*around the planar limit,

*δ*is the Tolman length, and $ks=2k+k\u0304$ is the spherical rigidity, with

*k*and $k\u0304$ being the bending and Gaussian rigidity, respectively. All these coefficients can be computed using either DGT or DFT.

^{50,51,54}Such an expansion can be performed for an arbitrary choice of dividing surface, but in this work, we use the surface of tension.

*γ*(

*R*) =

*γ*

_{0}. Incorporating this expansion into CNT, one obtains a curvature corrected expression for the work of formation of critical clusters,

^{52}

^{,}

^{50}and water-alcohol mixtures.

^{52}In this work, we shall perform a more stringent test of the accuracy of this formula for the LJTS(2.5

*σ*) using the values of the coefficients calculated from DGT according to the procedure of Ref. 51.

### F. Scaling function for the work of formation

Since the pioneering work of McGraw and Laaksonen on scaling relations,^{26} 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.

*f*(Δ

*P*) = 1 − Δ

*P*/Δ

*P*

_{sp}, proposed (in terms of Δ

*μ*) by Talanquer.

^{69}Along the same lines, different linear and parabolic scaling functions have been proposed in the literature.

^{72,74–79}In fact, the c-CNT expression, Eq. (18), can be rewritten as

*P*

_{sp}, where the isothermal compressibility diverges).

A very similar parabolic scaling function was recently proposed by Kashchiev^{79} 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.

*f*(Δ

*P*), using all the ingredients that we know, namely, (i)

*f*(0) = 1, (ii) the slope at the binodal given by the Tolman length

*δ*, (iii) the spherical rigidity

*k*

_{s}, (iv) zero work of formation at the vapor spinodal $f(\Delta Pspv)=0$, and (v) at the liquid spinodal $f(\Delta Psp\u2113)=0$. Instead of using

*k*

_{s}to specify the second derivative at the binodal, it is more convenient to use it to enforce the location of the maximum of the work of formation. We have enforced that the maximum of

*f*occurs at

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.

## III. SIMULATION DETAILS

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,^{81} 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.

*n*

_{ℓ}was used as a reaction coordinate in the harmonic umbrella biasing potential,

*k*

_{n}and

*n*

_{0}set the width and location of the window. A total of 15–40 windows located at intervals of 20–150 molecules with typical values of

*k*

_{n}= 0.005 − 0.02 were simulated at each temperature and pressure. For each window, 500k MC steps were used for equilibration, followed by 500k MC steps used for production. To speed up the simulations, a staging scheme was implemented, where a short trajectory of 25 MC steps was performed without bias and then accepted with a probability given by $exp\u2212\beta \Delta Wb(n\u2113)$, where Δ

*W*

_{b}(

*n*

_{ℓ}) is the difference in biasing potential before and after the trajectory. At each temperature and pressure, we also performed an unbiased simulation to obtain the cluster size distribution $P(n)$ of all clusters. In each window, the probability distribution function of the

*largest*cluster size $P\u2113(n\u2113)$ was recorded. The weighted histogram analysis method described in Ref. 82 was then used to reconstruct the free energy from the umbrella sampling simulations, joining the different windows. In order to obtain the free energy of

*any*cluster size,

*β*Δ

*G*(

*n*), from the probability distribution of the

*largest*cluster $P\u2113(n\u2113)$, we used the formula

*n*

_{min}is the size that we used to match the two distributions (typically

*n*

_{min}= 5). Figure 5 in the supplementary material shows a typical example of the reconstructed free energy landscape for

*T*= 0.900 and

*p*= 0.037. The location and value of the maximum of the reconstructed free energy

*β*Δ

*G*(

*n*) provide the size Δ

*N*and energy of formation of the critical cluster

*W*.

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*_{B}*T*. 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,^{81} 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.

## IV. RESULTS AND DISCUSSION

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.

### A. Comparison of theories and simulation results

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.^{30} 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.*^{55} 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*_{B}*T* 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* < 100*k*_{B}*T*; for higher barriers, the nucleation rate is generally too low to be observable. For nearly all droplets for which *W* < 100*k*_{B}*T*, 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=(\gamma /\gamma 0)3$, where *γ*_{0} 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 *γ*_{0} from the two theories (cf. Table I), as well as the different Helfrich coefficients $\delta ,k,k\u0304$. 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 $\gamma 0DFT\u2212\gamma 0DGT$.

. | γ_{0}
. | δ
. | k
. | $k\u0304$ . | $2k+k\u0304$ . |
---|---|---|---|---|---|

T = 0.625 | |||||

DGT | 0.732 11 | −0.081 26 | −0.322 32 | 0.192 49 | −0.452 16 |

DFT | 0.715 81 | −0.072 81 | −0.323 49 | 0.200 50 | −0.446 47 |

T = 0.700 | |||||

DGT | 0.585 99 | −0.078 29 | −0.323 00 | 0.189 84 | −0.456 16 |

DFT | 0.580 60 | −0.070 32 | −0.336 53 | 0.204 46 | −0.468 59 |

T = 0.741 | |||||

DGT | 0.509 38 | −0.076 35 | −0.323 10 | 0.188 45 | −0.457 75 |

DFT | 0.507 58 | −0.068 40 | −0.341 72 | 0.205 52 | −0.477 92 |

T = 0.800 | |||||

DGT | 0.403 18 | −0.073 18 | −0.322 17 | 0.186 09 | −0.458 24 |

DFT | 0.404 30 | −0.065 01 | −0.345 78 | 0.205 16 | −0.486 39 |

T = 0.900 | |||||

DGT | 0.235 23 | −0.067 36 | −0.313 89 | 0.178 93 | −0.448 85 |

DFT | 0.237 23 | −0.058 28 | −0.338 83 | 0.197 00 | −0.480 67 |

T = 1.000 | |||||

DGT | 0.088 27 | −0.064 58 | −0.281 11 | 0.158 62 | −0.403 59 |

DFT | 0.088 89 | −0.054 40 | −0.296 58 | 0.169 26 | −0.423 90 |

. | γ_{0}
. | δ
. | k
. | $k\u0304$ . | $2k+k\u0304$ . |
---|---|---|---|---|---|

T = 0.625 | |||||

DGT | 0.732 11 | −0.081 26 | −0.322 32 | 0.192 49 | −0.452 16 |

DFT | 0.715 81 | −0.072 81 | −0.323 49 | 0.200 50 | −0.446 47 |

T = 0.700 | |||||

DGT | 0.585 99 | −0.078 29 | −0.323 00 | 0.189 84 | −0.456 16 |

DFT | 0.580 60 | −0.070 32 | −0.336 53 | 0.204 46 | −0.468 59 |

T = 0.741 | |||||

DGT | 0.509 38 | −0.076 35 | −0.323 10 | 0.188 45 | −0.457 75 |

DFT | 0.507 58 | −0.068 40 | −0.341 72 | 0.205 52 | −0.477 92 |

T = 0.800 | |||||

DGT | 0.403 18 | −0.073 18 | −0.322 17 | 0.186 09 | −0.458 24 |

DFT | 0.404 30 | −0.065 01 | −0.345 78 | 0.205 16 | −0.486 39 |

T = 0.900 | |||||

DGT | 0.235 23 | −0.067 36 | −0.313 89 | 0.178 93 | −0.448 85 |

DFT | 0.237 23 | −0.058 28 | −0.338 83 | 0.197 00 | −0.480 67 |

T = 1.000 | |||||

DGT | 0.088 27 | −0.064 58 | −0.281 11 | 0.158 62 | −0.403 59 |

DFT | 0.088 89 | −0.054 40 | −0.296 58 | 0.169 26 | −0.423 90 |

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.

### B. The Tolman length and rigidity constants

We will next compare our estimates of the Tolman length to those of van Giessen and Blokhuis,^{84} 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 $\gamma 0\u221e$ 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.*,^{50} namely, *δ* = −0.074*σ*. Wilhelmsen *et al.* used another implementation of DGT with a less accurate equation of state.^{50} 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\u0304$ 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.*^{50} How the *W*/*W*_{CNT} curve varies with temperatures is therefore strongly connected to the value of the Tolman length.

### C. Discussion of a possible influence of capillary wave contributions to the planar surface tension

*T*= 0.900, van Giessen and Blokhuis

^{84}calculated the planar surface tension of the LJTS(2.5

*σ*) fluid for increasingly large systems (up to box sizes of 50

*σ*× 50

*σ*× 100

*σ*) and extrapolated to the thermodynamic limit. This yielded the value $\gamma 0\u221e=0.227774$, which is 3.2% lower than the value from DGT, $\gamma 0DGT=0.23523$, and 4% lower than the value from DFT, $\gamma 0DFT=0.23723$. We have attempted to refit the influence parameter of DGT to the presumably more accurate value from van Giessen and Blokhuis

^{84}and found that this

*removes*the good agreement between DGT and the simulations shown in Fig. 4 (cf. Fig. 4 in the supplementary material). One possible reason for this is the larger influence of capillary waves on large surfaces. Using the formulation of Planková

*et al.*,

^{85}capillary waves are expected to decrease the surface tension

*γ*according to

*T*

_{c}is the critical temperature, and

*T*

_{c}= 1.085 for LJTS(2.5

*σ*). The capillary wave contribution to the surface tension is the largest close to the critical point, with a limiting value of $\gamma bare/\gamma 0\u221e=1.044$ at the critical temperature. For

*T*= 0.900, this ratio equals 1.0364, which is in reasonable agreement with $\gamma 0DGT/\gamma 0\u221e=1.0327$. This may support the conclusions of Planková

*et al.*,

^{85}namely, that the small length scales of critical droplets tend to inhibit capillary waves, and it is therefore the bare surface tension that is relevant. It also seems to support the findings of Checa

*et al.*,

^{86}which indicate that DFT only captures fluctuations at length scales of (10 ± 2)

*σ*and that large capillary waves must be separately accounted for.

## V. CONCLUSIONS

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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

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

## DATA AVAILABILITY

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

## REFERENCES

*Metastable Liquids: Concepts and Principles*

*Nucleation: Basic Theory with Applications*

*Nucleation in Condensed Matter*

*Classical Nucleation Theory in Multicomponent Systems*

*Fundamentals of Inhomogeneous Fluids*

*Molecular Theory of Capillarity*

*Thermodynamics of Small Systems*

*Computer Simulation of Liquids*

*Understanding Molecular Simulation: From Algorithms to Applications*