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

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
2][3][4][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. 6Together 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: where 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][14] 150 years after the pioneering work of Gibbs, 6 nucleation phenomena keep fascinating scientists.6][17][18][19][20][21][22][23][24][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. 1or pure fluids, CNT is qualitatively correct. 26However, 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. 27CNT 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. 28Its computational simplicity and robustness make DGT attractive.However, it is only the first order approximation to the full classical density functional theory (DFT). 29The accuracy and predictive ability of DGT has, therefore, been questioned, in particular at the highest metastabilities. 30Using full DFT has been touted as a promising alternative. 31,32In previous studies, DFT has been used with great success to estimate a range of important interfacial properties, such as surface tensions [33][34][35][36] and adsorption isotherms. 37The 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. 39or many applications, it is impractical and too timedemanding to use DGT or full DFT to calculate the work of formation. 5,40This 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.][43][44][45][46][47][48] Recently, it was shown how the Tolman length and rigidity constants can be calculated using DGT [49][50][51][52] or DFT. 53,54Incorporating the curvature corrections into nucleation theory has been shown to significantly improve agreement with experimental results for water. 49It has also solved some of the inconsistencies that arise when CNT is extended to multi-component systems. 52The 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. 55n 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 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
Consider a metastable vapor at temperature T and chemical potential μ.For a spherical critical cluster, the work of formation is given by the formally exact relation, 6,56 where γ is the surface tension of the cluster with the surface of tension as a dividing surface.The Laplace pressure is defined as i.e., the pressure difference between the (stable) bulk liquid and the (metastable) bulk vapor having the same intensive properties (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, Here, R is the radius of the critical cluster at the surface of tension. 51,57In this work, the quantities R, γ (with no subscript) are always associated with the surface of tension.
For nucleation in an open system, the formation of a critical droplet will increase the number of particles in the volume compared to the homogeneous metastable vapor.One definition of the excess number of particles is 58 which is independent of the dividing surface.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.

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
The classical nucleation theory (CNT) formula for the nucleation barrier is obtained from Eq. ( 2) by using the capillarity approximation, which assumes that the surface tension γ is the same as that of the bulk liquid with a planar interface γ 0 , Accordingly, the ratio of the exact work of formation of a spherical critical droplet to that of the CNT-approximation is Consequently, 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, where α = (∂γ/∂ΔP) T .Therefore, if the surface tension has a nonzero 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
Density gradient theory (DGT) models the grand free energy at the state defined by (T, μ, V) by the following functional of the density profile: The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
The influence parameter κ 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, δΩ[ρ(r)]/δρ(r) = 0. To solve the Euler-Lagrange equations for DGT, we used the approach described in Ref. 51.

D. Density functional theory
Heier et al. 59 derived a density functional theory (DFT) for the LJTS(2.5σ)fluid based on a weighted density average (WDA) approach. 34To 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. 64Density functional theory models the grand free energy at the state defined by (T, μ, V) by the following functional of the density profile: where A[ρ(r)] is the Helmholtz energy functional.For the LJTS(2.5σ)fluid, the Helmholtz energy has contributions from the hard-sphere reference (FMT), dispersion (disp), and an ideal gas where Φ is the reduced Helmholtz energy density, and n(r) is a weighted density.The weighted densities are calculated via convolutions with the density profiles.For the contribution from FMT, the weights are listed in Ref. 65.The weighted density for the dispersion term is given by where ω disp is a weight function of a sphere with diameter φd BH , Here, Θ is the Heaviside function, 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: where F is the specific reduced Helmholtz energy.
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,66In addition, convergence of the overall equation system can be accelerated by use of Anderson mixing. 67auer 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 Re 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
The Helfrich expansion is based on a quadratic expansion of the surface tension in terms of the curvature 1/R around the planar limit, where δ is the Tolman length, and ks = 2k + k is the spherical rigidity, with k and k being the bending and Gaussian rigidity, respectively.All these coefficients can be computed using either DGT or DFT. 50,51,54Such an expansion can be performed for an arbitrary choice of dividing surface, but in this work, we use the surface of tension.
The Helfrich expansion can be viewed as a second-order refinement of the zeroth-order capillarity approximation γ(R) = γ 0 .Incorporating this expansion into CNT, one obtains a curvature corrected expression for the work of formation of critical clusters, 52 where The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp is the CNT prediction for the radius of the critical cluster.The physical significance of this curvature corrected CNT is corroborated by its success in improving the predictions of classical nucleation theory for the experimental nucleation rates of water 50 and water-alcohol mixtures. 52In 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
0][71][72][73] Accordingly, many expressions for this scaling function f have been proposed in the literature, 69,72,[74][75][76][77][78][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 (ΔPsp) = 0, where ΔPsp represents the difference in bulk pressures between the stable and metastable phases at the spinodal.
The most straightforward scaling function fulfilling these two requirements is the linear function f (ΔP) = 1 − ΔP/ΔPsp, proposed (in terms of Δμ) by Talanquer. 695][76][77][78][79] In fact, the c-CNT expression, Eq. ( 18), can be rewritten as However, c-CNT is derived using only the Helfrich expansion, and although it predicts a vanishing of the nucleation barrier, this does not necessarily agree with the location of the thermodynamic spinodal (i.e., at ΔPsp, 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 (ΔPsp) = 0, and the value of Tolman length δ as conditions to evaluate the coefficients.
We can formulate a more accurate expression for the scaling function 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 ks, (iv) zero work of formation at the vapor spinodal f (ΔP v sp ) = 0, and (v) at the liquid spinodal f (ΔP ℓ sp ) = 0. Instead of using ks 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 namely, the same maximum as obtained from Eq. ( 20).An accurate scaling function that satisfies these five conditions is a fourth-order polynomial in the Laplace pressure, 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 ΔPm.
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 rc, corresponding to the first minimum in the radial distribution function of the liquid at the desired temperature.
The number of liquid-like molecules in the largest liquid cluster n ℓ was used as a reaction coordinate in the harmonic umbrella biasing potential, The Journal of Chemical Physics

ARTICLE scitation.org/journal/jcp
where kn 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 kn = 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(−βΔW b (n ℓ )), 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 ℓ (n ℓ ) 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 ℓ (n ℓ ), we used the formula where 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 kBT.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.The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpagainst 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. 30This 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).

A. Comparison of theories and simulation results
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 < 100k B T; for higher barriers, the nucleation rate is generally too low to be observable.For nearly all droplets for which W < 100k 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/W CNT = (γ/γ 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 δ, 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 γ DFT 0 − γ DGT 0 .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/Psat, 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 γ ∞ 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., 50 namely, δ = −0.074σ.Wilhelmsen et al. used another implementation of DGT with a less accurate equation of state. 50This 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 ks ∼ −0.45 for DGT (cf.Table I).This is in line with the findings by Wilhelmsen et al. 50How 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
For 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 γ ∞ 0 = 0.227 774, which is 3.2% lower than the value from DGT, γ DGT 0 = 0.235 23, and 4% lower than the value from DFT, γ DFT 0 = 0.237 23.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 where Tc is the critical temperature, and Tc = 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 γ bare /γ ∞ 0 = 1.044 at the critical temperature.For T = 0.900, this ratio equals 1.0364, which is in reasonable agreement with γ DGT 0 /γ ∞ 0 = 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

ARTICLE
scitation.org/journal/jcpaim 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.

Figures 2 -FIG. 2 .FIG. 3 .FIG. 4 .
Figures 2-4 compare the number of particles and the work of formation predicted from all theories considered in this work

The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp 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.

TABLE I .
Planar surface tension and Helfrich coefficients for the LJTS(2.5σ)fluid.