On the role of metastable intermediate states in the homogeneous nucleation of solids from solution^{1}^{1}1To appear in Advances in Chemical Physics
Abstract
The role of metastable liquid phases in vaporcrystal nucleation is studied using Density Functional Theory(DFT). The bulk free energy functional is modeled using a hybrid model consisting of a sum of a hardsphere functional (Fundamental Measure Theory) and a liquidlike contribution to account for the attractive part of the interactions. The model gives a semiquantitatively accurate description of both the vaporliquidsolid phase diagram for both simple fluids (LennardJones interactions) and of the lowdensity/highdensity/crystal phase diagram for model globular proteins (ten WoldeFrenkel interaction). The density profile is characterized by two local order parameters, the average density and the crystallinity. The bulk free energy model is supplemented by squaredderivative terms in these order parameters to account for inhomogeneities thus producing a model similar in spirit to phasefield theory. It is shown that for both interaction models, the vaporcrystal part of the phasediagram can be separated into regions for which metastable liquid phases are more or less stable than the vapor, but always less stable than the solid. The former case allows for the possibility of double nucleation whereby liquid droplets nucleate from the vapor followed by a separate nucleation of the solid phase within the liquid droplets. Whether or not this actually occurs depends on the relative free energy barriers for vaporsolid and vaporliquid nucleation and it is shown that for simple fluids, double nucleation is indeed favored at sufficiently large supersaturation. Finally, by studying the minimum free energy path from the vapor to the solid, the separate possibility of transient nucleation, where the vaporsolid transition involves a single nucleation event but the subcritical clusters tend to be liquidlike while small, is shown to also be possible when the metastable liquid exists, but the supersaturation is too small for double nucleation.
pacs:
64.60.Q, 64.70.Hz, 64.60.My, 82.60.Nh, 68.08.pContents
I Introduction
The goal of this work is to describe, at a semimicroscopic level, the process of homogeneous nucleation of solids from solution. This is studied in the approximation that the effects of the solvent can be accounted for by an effective interaction between the solvate molecules so that at low concentrations, the molecules in solution are treated as a singlecomponent gas with the crystalline phase being the corresponding solid phase in the effective single component system. In this picture, there is also the possibility of a dense liquidlike phase which simply corresponds to the liquid phase of the effective singlecomponent system. Thus, the following work is intended to be applicable both to nucleation in singlecomponent systems, such as simple fluids, and to nucleation of a solid phase from macromolecules in solution. For the sake of consistency, most of the discussion will be phrased in the language of nulceation of solids in a single component system. So, in the following, all statements concerning the “vapor phase” should be understood to be equally applicable to the “lowconcentration” or “low density” phase of molecules, such as proteins, in solution and statements concering the “liquid phase” are equally applicable to the “highconcentration” or “highdensity” phase of a solution.
In general, homogeneous nucleation involves at least two phases: the stable phase, S, and the metastable phase, M. Initially, the system consists of M in its homogeneous, bulk state. Fluctuations give rise to small volumes of S which will here always be taken to be spherical clusters. By definition, the free energy of bulk S is lower than that of bulk M but finite clusters also involve an interface which has higher free energy, on average, than either phase. Hence, sufficiently small clusters are typically unstable, having higher free energy even than the equivalent amount of M. Furthermore, for small clusters, there is no reason to believe that the material inside the cluster is in the bulk S state since the interface can have a volume comparable to that of the bulk region. Thus, small clusters will be unstable with respect to M: indeed, will be unstable with respect to smaller clusters thus leading to cluster dissolution. For sufficiently large clusters, the interior will consist mostly of material near the bulk S state, and the volume of the interface will be small compared to the volume of the bulk so that the cluster is stable with respect to M but unstable with respect to cluster growth thus driving the transition of the whole system from M to S (see, e.g. Kashchiev ()).
In recent years, it has become apparent, first for proteinstWF (); VekilovCGDReview2004 (); LN (); Gunton (); GuntonBook () but then even for simple fluidsLN (); FrenkelLJ (); Chen (), that homogeneous nucleation of solids from vapor can be more complicated than this simple picture. This is because in many cases, the conditions for solid nucleation also allow for a second metastable phase, namely that of the bulk liquid. In cases where the bulk liquid is more stable than the bulk vapor (i.e. has lower free energy) but less stable than the solid, it can be that, instead of directly nucleating a crystalline cluster from the vapor, it is energetically favorable to first nucleate a liquid cluster which then grows and, subsequently, to complete the transition to the crystalline state via a second homogeneous (liquidcrystal) nucleation event within the liquid cluster (thus creating a growing solid cluster within the growing liquid cluster). This process will be referred to as ”double nucleation”. Whether or not this actually occurs will depend on the free energy barrier encountered in the vaporcrystal and vaporliquid transition. Alternatively, even if direct vaporcrystal nucleation is favored, it is still possible that the growing cluster nevertheless passes through a liquidlike stage which, however, is always subcritical, as a transient on the way to creating a solidlike critical cluster. In this case, small clusters would be liquidlike in structure, becoming more crystalline as they grow larger. This scenario, involving only a single nucleation event, will be referred to as ”transient twostep nucleation”. The simulations of ten Wolde and FrenkeltWF () are clear illustrations of transient twostep nucleation: they exhibit a nucleation pathway that involves liquidlike small clusters followed by solidification as the cluster approaches the critical size. Vekilov discusses both scenarios and suggests that even when double nucleation is possible and energetically favored, it may be suppressed by kineticsVekilovCGDReview2004 (). Kaschiev and Vekilov have analyzed the effect of double nucleation on observed nucleation ratesKashchievVekilov (). van Meel et al report simulation results showing double nucleation for a LennardJones fluidFrenkelLJ ().
The key to understanding these processes is the construction of models for the free energy of inhomogeneous  multiphase systems. Indeed, Classical nucleation theory (CNT), from which we take our lead, is fundamentally based on a description of the free energy as a function of the size of the cluster. In CNT (which should be thought of in terms of the simpler vaporliquid transition) the width of the cluster interface is taken to be zero and the interior of the cluster is assumed to be in the bulk state, S. The only variable is then the size of the cluster and it is assumed that homogeneous nucleation consists of the growth of a cluster from size zero until it is arbitrarily largeLutskoEPL (); Lutsko_JCP_2008_2 (). This, in other words, is understood to be the ”nucleation pathway”. This idea can be developed more rigorously to include finite width of the interface and nonbulk properties of the interior in which case the nucleation pathway involves a description of the variation of all these properties  size, interfacial width and interior density  simultaneouslyLutsko_InPrep ().
The process of solid nucleation is more complex, as the phenomenology above would imply. First, the vaporsolid transition involves at least two order parameters: the density and the crystallinityOLW1 (); TalOxtoby (); Gunton (); LN (). Since the typical solid density is close to that of the metastable liquid (if it exists), the double nucleation scenario involves a separation of these order parameters: the vaporliquid transition involves two different densities but occurs at zero crystallinity while the liquidsolid transition occurs at (nearly constant) density and involves a change from zero to finite crystallinity. In the transient twostep scenario, a small (unstable) cluster begins to form at zero crystallinity but at some point, as the cluster becomes larger, the crystallinity increases so that the critical cluster has both finite density difference compared to the vapor and finite crystallinity. Both of these can be viewed as ”twostep” scenarios compared to the possibility that, from the beginning, the crystallinity and density change at the same time. Clearly, description of these processes must involve a model for the free energy of inhomogeneous systems which is sufficiently detailed so as to describe liquid and solid clusters embedded in a vapor bulk. In order to capture the transient scenario, this must be supplemented by some means of determining the nucleation pathway.
Here, the basic tool for calculating the free energies will be classical, finitetemperature Density Functional Theory (DFT). If DFT were sufficiently welldeveloped and technically simple, nothing more would need to be said about this part of the study: unfortunately, neither of these conditions holds for the solid phase. (In contrast, direct calculations for liquidvapor systems are possibleLutsko_JCP_2008 (); Lutsko_JCP_2008_2 (). DFT is only sufficiently well understood so as to give an a priori description of bulk fluid, bulk solid and inhomogeneous fluidsolid systems for hardsphere interactions. For any more realistic potential  particularly those with attractive interactions  some sort of modelling is necessary. Typically, this will involve the introduction of an effective hardsphere diameter and the representation of the free energy as the sum of the hardsphere functional and a second term accounting for the attractive part of the interactions. Futherrmore, even for the hardsphere functional, the calculation of the free energy for an inhomogeneous system (i.e. of a solid cluster embedded in a fluid background) is very complicated and computationally expensiveSong_FMT (). Fortunately, an easier alternative is available. It has been shown that the exact free energy for a solid can be expanded in terms of gradients of the order parametersLowen1 (); Lowen2 (); LutskoGL () thus providing a connection between DFT and the older gradient theories of Cahn et alCahnHilliard (); CahnHilliardNucleation (), as well as the phasefield theories commonly used to study solidsolid transitionsphase_field (). This is the justification for using a gradient model in the present work. A significant advantage is that the free energy functional is then only needed for the case of homogeneous (i.e. bulk) systems, thus placing less stress on the accuracry of the DFT model. On the other hand, expressions for the coefficients of the gradient terms must be found. In principle, the exact results express these in terms of derivatives of the bulk free energy but in practice, they are hard to calculate except for the case of a fluid. In this work, a semiempirical procedure will be used to fix these terms. The main difference between the present work and that of Gunton et alGunton (), which is similar in spirit, is that the latter made use of toy free energy models whereas here realistic models are used which give semiquantitatively accurate bulk phase diagrams as well as of the liquidvapor phase transitionLutsko_JCP_2008 ().
Given these approximations, the free energy for any configuration of order parameters can be calculated. The practical exploration of the models will make use of energysurface techniques commonly applied to the study of chemical reaction pathways and structural transitionsWales (). A recently developed method  involving the approximation of the spatiallyvarying order parameters as piecewisecontinuous functions, will be used to determine the critical clusters  i.e. saddle points of the free energy  for homogeneous nucleation. This will already give enough information to identify if and when double nucleation is possible. The nucleation pathway will be identified, as is commonly done for chemical reaction pathways, using steepest descent paths. These identify the most likely path for the transition given the free energy surface and are a natural generalization of the CNT pathway. They do not take into account kinetic factors, such as rates of mass transport, that could play a significant role particularly for small molecules.
Section II will describe the technical details of the free energy models used. The bulk thermodynamics is used in Section III to limit the regions of the phase diagram in which double nucleation is possible. A simple model for double nucleation is also used to illustrate the role of bulk free energy differences and of surface tension. The fourth Section describes the application of the model to planar interfaces and illustrates the role of wetting. Detailed calculations of the energy barriers for direct nucleation of the crystal and of double nucleation are presented in Section V. The possibility of transient double nucleation is also described. The paper ends with a summary and with a discussion of future directions.
Ii The free energy model
ii.1 Density
The central quantity defining the system state in DFT is the local density, . It is possible to formulate the theory with no a priori restrictions placed on the density, but this is computationally expensive. A commonly used alternative is to represent the density in terms of a set of basis functions. For bulk crystalline systems, this usually means a sum of Gaussians localized at the lattice sites,
(1) 
where is a lattice vector, controls the width, is the occupancy and is the lattice density. The average density is
(2) 
The density can also be written in Fourier representation as
(3) 
where is a reciprocal lattice vector. Notice that in this representation, it is clear that the limit gives a uniform density which describes a fluid. It is therefore natural to characterize the crystallinity by the size of the amplitude of a typical nonzero wavevector term such as
so that the density can also be written as
(4) 
thus showing that the density is parameterized entirely by , and . In the following, we neglect the variation of the lattice density and generalize to inhomogeneous systems by allowing the occupancy and the crystallinity to depend on position. It is also convenient then to replace the occupancy by the average density so that we have
(5) 
The order parameters are then the local average density, , and the crystallinity, . It will sometimes be more convenient to replace the latter by the amplitude of the smallest nonzero wavevector, .
ii.2 Gradient expansion
In order to determine the density, a model for the (grand canonical) free energy functional, is necessary. Good models exist for liquids and can be used to study, e.g., the liquidvapor transition in great detailLutsko_JCP_2008 (). However, the theory is less well developed for the solid phase and in any case calculations for inhomogeneous solids are very expensive. The present work therefore makes use of a gradient expansion of the free energy which focusses attention on the order parameters and only requires information about homogeneous solidsLowen1 (); Lowen2 (); LutskoGL (). The grandcanonical free energy is written as
(6) 
where is the chemical potential. In general, if the density can be written in terms of order parameters, as
(7) 
so that the density of the uniform, bulk system is
(8) 
and if in some sense ”slowly varying”, then the squaredgradient approximation (SGA) for the functional is
(9) 
where the summation convention is used. The first term of the free energy involves the bulk free energy density defined as
(10) 
For the order parameters used here, the free energy is explicitly
(11) 
ii.3 Bulk Free Energy
The bulk free energy model used here is based on the idea of separating the free energy into a hardsphere contribution, for which the DFT is well developed, and a second contribution that accounts for the longranged attractive interactions. A particularly simple model is based on the observation that the local structure of an FCC solid and a simple fluid are quite similar so that, as a simplest approximation, one can imagine that the correction to the hardsphere model is independent of the local structureCurtin (); LN () giving
(12) 
where is the effective hardsphere diameter and
(13) 
is the contribution of the attractive part of the interaction to the liquid free energy. The effective hardsphere diameter is determine using the WCA expression as modified by Ree et alRee1 (); Ree2 ().
ii.4 Bulk phase diagram
In DFT, one is always working in the grandcanonical ensemble so the external parameters are the temperature, the chemical potential, and the applied external field, . The latter includes any confining walls: if the walls are hard, then the volume is a fixed parameter (as will always be assumed here). The appropriate free energy is the grand potential,
and it should be noted that all information about the state is encoded in the local density function, . The equilibrium states (i.e. density distributions) are determined by minimization,
(14) 
which is to say
(15) 
For parameterized profiles, the requirement that the free energy be a minimum gives
(16) 
and if the parameters are constants, then
(17) 
For a given value of chemical potential, there may be multiple solutions for the density in which case the equilibrium state is the one corresponding to the absolute minimum of the grand potential. Two phase coexistence therefore requires that
(18)  
In particular, using the chain rule for functional differentiation,
(19)  
gives the usual thermodynamic relation
(20) 
Thus, for uniform densities, and e.g. , the conditions for coexistence are
(21)  
which is to say equality of chemical potentials and of pressures since in the bulk, .
For the crystalline system, the parameterization used here involves not just the average density but also the crystallinity and the lattice parameter. The conditions for an extremum of the free energy are then
(22)  
In principle, the chemical potential (an external field) are specified and these equations solved for and : in practice, it is simpler to choose a value of the lattice density and to use these conditions to determine and . Further technical details are discussed in Appendix A.
Iii Thermodynamics of two step nucleation
iii.1 Independent variables and Ensembles
The main calculational tool used here, DFT, is formulated in the grandcanonical ensemble in which the independent variables are the temperature, chemical potential and volume (or, more generally, applied field) and the free energy of interest is the grand potential. Experiments are typically performed at constant pressure, temperature and volume for which the Gibbs free energy is relevant. Fortunately, in discussing nucleation, one is always interested in free energy differences and these are independent of the ensembleOxtoby_Evans (). In the following, since it is based on DFT, the grandcanonical ensemble is always used, however in many cases the difference between ensembles can be suppressed by focussing on physical quantities. Thus, rather than specify the chemical potential (for the grand ensemble) or the pressure (for the PVT ensemble) one can specify the density of the final phase which is a meaningful variable in both formulations and in either case is uniquely related to the state variable.
iii.2 Interaction potentials and the fluid phase
In this study, the nucleation properties of two different systems will be considered: simple fluids as described by the LennardJones potential,
(23) 
and globular proteins as described by the ten WoldeFrenkel model interaction,
(24) 
which will be studied here for the value as is appropriate to model the phase behavior of globular proteins.
Both of these potentials are longranged in the sense of decaying as powerlaws. In simulation, infiniteranged potentials are difficult to deal with, so any longranged potential is typically cutoff at some distance, . For MonteCarlo, a simple shift of the potential to make the energy continuous at the cutoff is typically used so that the socalled truncated and shifted potential is
(25) 
In molecular dynamics simulations, the force is usually required to be continuous so that the forceshifted potential is commonly used,
(26) 
where . Other forms are also important, particularly the BroughtonGilmer modification of the LJ potential:
(27) 
where , , , and BG1 ().
All of these details significantly affect the liquidvapor phase diagram. The DFT model described above requires as input both the interaction potential and the equation of state for the fluid phase. For the LJ potential with no cutoff, essentially exact empirical equations of state are available for temperatures above the triple pointJZG (); Mecke (); Mecke1 (). For finite cutoffs, these must be modified with inexact meanfield corrections. Although useful for studying liquidvapor coexistenceLutsko_JCP_2008 (); Lutsko_JCP_2008_2 (), they are of limited utility for vaporsolid nucleation since the interesting affects occur below the triple point. Less accurate, but more broadly applicable, are meanfield equations of state based on thermodynamic perturbation theory such as the WeeksChandlerAndersen (WCA) theoryWCA1 (); WCA2 (); WCA3 (); Ree1 (); Ree2 () which will be used here.
iii.3 Solid phase diagrams and intermediate phases
Figure 1 shows the calculated vaporliquidFCC solid phase diagram for the infiniteranged LJ potential and for the truncated and shifted tWF potential compared to simulation in terms of the reduced temperature . The LJ phase diagram possesses both a liquidvapor critical point and a triple point whereas for the protein model, the liquidvapor transition is metastable. Since the fluid equation of state is being calculated from thermodynamic perturbation theory, which is a mean field theory, the critical point is, as expected, poorly described for both potentials. Apart from this expected inaccuracy, the model is in semiquantitative agreement with simulation for all three phases.
iii.3.1 Nucleation scenarios for globular proteins
In order to clarify the thermodynamics of metastable states, we consider in more detail the phase diagram derived from the tWF potential. Figure 2 shows a line crossing the phase diagram at constant temperature. The points A and B identify coexisting vapor and solid phases which, by definition, have the same free energy and the same control parameter (chemical potential or pressure). Starting at the vapor point A and moving along the isotherm in the direction of increasing density, i.e. to the right, corresponds to increasing the chemical potential and to each value of the density there will be a point on the isotherm to the right of the solidbranch of the binodal, i.e. to the right of point B, having the same chemical potential; however, under these conditions the solid phase will have lower free energy than the vapor phase so that the vapor points to the right of A are metastable with respect to the solid phase. Eventually, as one moves along the isotherm, the point in the vapor region reaches the vaporliquid coexistence curve so that there is a coexisting liquid phase. However, by definition, this particular liquid state has the same free energy as the vapor, one knows that it has higher free energy than the solid and so can only be metastable. The remains true as one moves along the isotherm to densities to the right of the vapor branch of the liquidvapor binodal. Eventually, one reaches the vaporliquid spinodal and above this density, the vapor phase does not exist. The set of liquid points with the same chemical potential as the vapor points on the spinodal will therefore divide region to the right of the liquid binodal into a lowerdensity region, for which a vapor with the same chemical potential can always be found, and a higherdensity region with no corresponding vapor.
What is the role of the metastable liquid phase in vaporsolid nucleation? We can only address this question here with respect to double nucleation. (Transient metastable states are a property of the nucleation pathway and its presence or absence cannot be answered solely from a consideration of the bulk phases.) First, consider again the vaporliquid coexistence curve. By a reversal of the reasoning above, vapor states to the left of the vapor branch correspond to liquid points to the left of the liquidbranch of the coexistence curve and are more stable than the corresponding liquid state. Moving to the left away from the vapor branch, the corresponding liquid point also moves left until it reaches the liquidvapor spinodal: beyond this point, there is no corresponding liquid. This therefore defines a line in the phase diagram with the property that vapor states to the left have no corresponding liquid states, shown as a red line in Fig. 2. In particular, points on the vapor branch of the vaporsolid coexistence curve have no corresponding liquid. This line therefore divides the metastable vapor region into two parts: a low density part in which there is no corresponding liquid state (i.e. no liquid with the same chemical potential) that can play a role in nucleation and a higher density region where the corresponding liquid exists as a metastable state. Hence, for systems prepared with the vapor density between the vaporbranch of the solidvapor coexistence curve and the new demarcation line, double nucleation is not possible as there is no metastable liquid. The region within which vapor and liquid states with the same chemical potential can be found is therefore an envelope around the binodal and will be referred to as the nodal curve.
Starting at the vapor branch of the vaporsolid coexistence curves, and moving to the right we therefore find that:

From coexistence to the nodal line, there is no liquid phase and so, double nucleation is not possible. In this region, the solid free energy is less than that of the vapor.

From the nodal line, to the vapor branch of the vaporliquid coexistence curves, there is a liquid state but it has higher free energy than the vapor (which has higher free energy than the solid). Nucleation via a liquid cluster is possible, with the solid cluster nucleating within the liquid but the liquid cluster would always be metastable. This is a form of transient twostep nucleation.

From the vapor branch of the liquidvapor coexistence curve to the vaporliquid spinodal, the liquid has lower free energy than the vapor but higher than the solid. True double nucleation is possible depending on the barriers for the formation of a critical liquid cluster within the vapor, a critical solid cluster in the liquid as compared to the barrier for directly forming a critical solid cluster in the vapor. Even if the latter is favored, a transient scenario is possible.
All of this is shown in a more direct way for the tWF potential in Fig. 3 where supersaturation, with the vapor pressure and the vaporsolid coexistence pressure, is used as the independent variable. By definition, for , the vapor is the stable phase and for the solid is the stable phase. The latter region is divided into three sections: that for which there is no corresponding liquid state, that in which there is a liquid but it has higher free energy than the vapor and that in which there is a liquid state with lower free energy than the vapor but higher free energy than the solid. Only in the latter case is liquid nucleation possible so it is only in this region of the phase diagram that double nucleation can occur. This therefore represents a set of necessary conditions for double nucleation. Nucleation pathways that do not involve double nucleation but which pass through liquidlike state, i.e. transient twostep nucleation, can in principle occur for any value of . The primary goal of detailed DFT calculations is to determine the necessary conditions for double nucleation and to assess for what conditions, if any, transient twostep nucleation occurs.
iii.3.2 Nucleation scenarios for simple fluids
In some ways, the phase diagram for simple liquids is more complex because the various coexistence curves cross. As shown in Fig.2, the calculations indicate that there is liquidvapor coexistence below the triple point. However, just as in the case of the globular proteins, the vapor branch of the vaporliquid coexistence curve is to the right of the vapor branch of the vaporsolid coexistence curves thus indicating that the coexisting liquid and vapor have higher free energy than does the solid phase (at the value of chemical potential or pressure corresponding to liquidvapor coexistence). Hence, below the triple point, the liquid phase is again metastable just as in the case of the globular proteins. Similarly, below the triple point the liquid branch of the liquidsolid coexistence curves lies to the left of the liquid branch of the liquidvapor coexistence curve: i.e. it is in the metastable (or even unstable) region of the liquidvapor phase diagram. So, these liquid states have higher free energy than the corresponding vapor phase and the coexisting liquid and solid phases are metastable with respect to the vapor. Taking all of this together, only the vaporsolid transition is usually drawn below the triple point, but there is a metastable liquid and a liquidvapor coexistence curve in this region. It is this metastable liquid phase which gives rise to the possibility of twostep nucleation even for simple liquids. However, the fact that the vapor branches of the vaporsolid coexistence curves and the vaporliquid coexistence curves are very close together makes it more difficult to display the various metastability boundaries. Figure 4 shows the phase information in terms of the supersaturation and the similarities and differences to the model protein behavior are evident. The main difference is that the intermediate liquid state exists for all values of the supersaturation leaving only the two regions distinguishing liquids which are less or more stable than the vapor phase.
iii.4 A simple picture of double nucleation
In this subsection, the goal is to anticipate the following, more technical, developments to give an idea of how doublenucleation can be described by something that looks like an extension of CNT. We begin with the free energy functional defined above specialized to spherical symmetry,
(28) 
where is the grand potential per unit volume and the coefficients of the gradient terms are taken to be constants. It is assumed that the temperature and chemical potential are such that there are three bulk states: a solid with order parameters , a liquid with order parameters and a vapor with parameters .We now introduce a very simple model for the order parameters for a (spherically symmetric) vaporsolid cluster which nevertheless captures the important physics of the real system:
(29)  
This piecewiselinear model can be extended to give arbitrarily complex approximations and will play an important role below. For now, it is substituted into the expression for the free energy which is then simplified to get
(30)  
with
(31)  
In the present discussion, it is assumed that the clusters are not small (in the sense that ) so that lower order terms in the expression for can be neglected.
The critical cluster is a stationary point of the free energy so that the width is determined from
(32) 
and the radius from
(33)  
giving
(34) 
The same model applied to the vaporliquid and liquidsolid critical clusters gives
(35)  
Finally, let us make the further approximation that the density of the liquid and solid states are almost the same so that
(36)  
Double nucleation will be possible if and will be the energetically preferred pathway provided which in turn is more likely if one or more of these conditions is fulfilled:

The barrier for a direct transition is larger than that for an indirect transition (It was this condition that was studied previously by Lutsko and NicolisLN ()).


are not small compared to .
Clearly, the factor occurring in the denominator of can be important in raising the vaporliquid barrier compared to the vaporsolid barrier. On the other hand, the factor involving the gradient coefficients is always going to be smaller for the vaporliquid cluster than for the vaporsolid cluster due to the terms related to crystallinity.
Iv Gradient coefficients and planar interfaces
In the standard CNT model, the excess free energy of a cluster is the sum of a bulk term and a surface term with the latter being proportional to the planar surface tension at coexistence. In the same way, for the model considered here, the planar surface tension will play a key role in determining the gradient coefficients.
iv.1 Gradient coefficients
Since the secondgradient approximation is the result of a formal expansion of the free energy, exact expressions for the gradient coefficients exist (see Appendix B). Their evaluation requires full knowledge of the direct correlation function in the bulk system for all densities and crystallinities which is of course not known. While reasonable models could be used to make approximate evaluations of the coefficients, the results would probably be disappointing when used to evaluate physical quantities such as the liquidsolid or vaporsolid planar surface tension simply because the squardgradient model is a crude approximation. In fact, the relevant interfaces tend to have widths of only a few times the typical atomic separation so that the assumption of slowlyvarying order parameters that underlies the SGA is unlikely to be very good, although for the particular case of the liquidvapor interface which involves only a single order parameter and no structural change, it is actually rather goodLutsko_InPrep (). Thus, while it would be interesting to make good evaluations of the exact expressions for the coefficients, the approach used here is more pragmatic. First, it is noted that a systematic expansion in crystallinity and density gives
(37) 
with
Second, there is evidence that the densitydensity coefficient can be well modeled in the liquid, i.e. for , by a densityindependent quantityLutsko_InPrep (). Thirdly, explicit calculations in the accessible limit of low density indicate that all three coefficients are relatively insensitive to the crystallinity (aside from the explicit factors shown above) at least up to crystallinities about half that of the bulk solid. Finally, Laird has notedLaird_solid_liquid_hs (); Davidchack_crystal_melt () that the structural properties tend to be dominated by the shortrange repulsion of the pair potential so that the liquidsolid surface tension can be approximated by that of the hardsphere solid evaluated with an effective hardsphere diameter giving the useful approximation that the excess surface free energy for a planar liquidsolid interface is
(38) 
. This suggests that and should be dominated by hardsphere contributions which would imply that they scale linearly with temperature. It is also consistent with the model free energy functional used here in which all of the dependence of the free energy on the crystallinity enters through the hardsphere part of the free energy. All of this suggests a simple approximation for the structural coefficients along the lines of
(39)  
In a final simplification, the lowdensity limit gives .
iv.2 Planar interfaces
Figure 5 shows the result of using this model to calculate the liquidsolid, vaporsolid and vaporliquid surface tensions at various temperatures as compared to the available simulation data. It is evident that, e.g., choosing to give, e.g., a value of the liquidsolid surface tension in agreement with a single point of either simulation or the DavidchackLaird approximation is enough to give a reasonable description of the liquidsolid and vaporsolid surface tensions over a range of temperatures. Furthermore, reasonable values for all of the physical quantities are found for a range of choices of the coefficient so that the model is relatively robust with respect to this choice.
Figure 6 shows the planar profiles calculate using for a temperature below the triple point, one near the triple point and one above the triple point. The profiles above and below the triple point share common features: the width and position of the transition regions for both the average density and the crystallinity are roughly the same and the overall width of the interfaces is about three atomic planes. The profile near the triple point is broader and is actually composed of two distinct regions: the first part in which there is a modest drop in density and the crystallinity goes to zero, and the second region in which the density drops to that of the vapor with the crystallinity equal to zero. The first region has the nature of a solidliquid transition while the second has the structure of a liquidvapor transition. In fact, this can be interpreted as a caricature of a wetted interface. True wetting, with a liquid region of finite width, is not possible in this model because the longranged van der Waals forces that give rise to the fluid are not being explicitly modeled. Hence, the results shown are as close as this model can come to representing wetting  basically, wetting with the bulk fluid region having zero width.
V VaporCrystal nucleation
v.1 General considerations
Given a model for the free energy of interfacial systems it is now possible to consider the process of nucleation. In this context, nucleation means the formation of clusters that eventually become supercritical. As mentioned above, the most important issue that arises in vaporsolid nucleation is the description of the nucleation pathway, for which there are at least three possibilities. The first is the conventional pathway in which the subcritical cluster has essentially bulk solid properties (density and crystallinity) except for an interfacial region. When the clusters are small, they are subcritical and the free energy increases with cluster size until a transition state  a critical cluster  is reached, after which further growth lowers the free energy. This pathway is therefore characterized by a single nucleation barrier and by nearbulk solid properties of the interior of the cluster. The second pathway involves double nucleation. First, a purely liquidlike cluster forms and grows until it reaches a critical size after which it is supercritical and stable with respect to the vapor. Within this liquid cluster, a second cluster forms, this time having the properties of the bulksolid, and goes through a similar process of subcritical growth reaching a transition state and then becoming stable with respect to both the vapor and the liquid. This pathway is therefore characterized by two nucleation events and two energy barriers. The final possibility is termed transient twostep nucleation and is in some sense intermediate between the classical and doublenucleation scenarios. Sufficiently small clusters, which are of course subcritical, are liquidlike but at a certain size, below the critical size for the liquid, the crystallinity begins to increase so that there is a single critical cluster, perhaps more solidlike than liquidlike, and a single energy barrier to be crossed on the path towards crystallization. In this case, the liquid need not even be stable with respect to the vapor  indeed, there is not even a priori reason why the liquid must exist as a thermodynamic state (e.g. minimum of the bulk free energy) at all.
The question therefore arises as to how one uses the free energy model to characterize the nucleation pathway? Nucleation is of course a rare, noisedriven event and a dynamical description would seem most appropriate. In fact, in this sense, characterizing nucleation pathways is similar to the problem of characterizing chemical reaction pathways, for which the same issues occur.
v.2 Double Nucleation
The first step in characterizing any structural transition or reaction is the identification of the saddle points of the free energy surface, i.e. the transition states. In general, the beginning state (the bulk vapor) and the end state (the bulk crystal) are local minima of the free energy and the transition states are the critical nuclei which define the energy barrier separating the minima. Here, it is assumed that the relevant density distributions are spherically symmetric and the order parameter fields, and are again modeled by piecewise continuous functions,
with a similar model for the crystallinity. In order to refer to different realizations of these models, the notation will be used indicating that there are links for the density profile (which means parameters since there is the initial radius and then one density and one width for each link) and links for the crystallinity profile, giving parameters for a total of parameters. The free energy is then a function of those parameters and the transition states are identified by searching the parameter space for stationary points of the free energy surface. This is done using standard eigenvectorfollowing techniquesWales ().
Three specific models will be studied here. First is the ”CNT” model in which there is a single link in both the density and crystallinity profiles together with the additional constraint that the radii and widths of the two profiles are the same and that the density and crystallinity in the bulk region are the same as for the bulk crystal for the given temperature and chemical potential. This model therefore depends on only two parameters (the radius and width of the profiles) and is the simplest possible model. The other models studied will be , a single link for each profile, and , in which there are two links in the density profile. As for the planar interfaces described above, this is necessary to allow for wetting of the surface and will be seen to lead to a substantial reduction of the free energy over the single link model. On the other hand, including additional links in the crystallinity has very little effect. (For example, at and , the change in free energy of the critical nucleus found using the model and the model is on the order of one percent.)
The interesting feature of this problem is that for sufficiently high supersaturations, both the liquid and the solid are more stable than the vapor so that there are three transition states: one corresponding to the vaporliquid transition, another to the vaporsolid transition and a third to the liquidsolid transition. All of this follows simply from the bulk free energies as discussed above and illustrated in Fig.7. The question addressed here is which of the possible paths will occur: direct transition from vapor to solid or a twostep transition via first a vaporliquid transition and then a liquidsolid transition. It is assumed that whichever transition involves the lowest free energy barriers will dominate.
Tables 1 and 2 give the relevant free energy barriers for transitions at different supersaturations for and respectively. It is found that for low supersaturation, corresponding to the case of large critical nuclei, the barrier for the direct vaporsolid transition is lower than that for the vaporliquid transition. However, at larger supersaturations, the vaporliquid transition becomes less costly thus implying that the doublenucleation scenario is favored.
Figures 8 show the structure of the critical nuclei for a few cases. In all cases, the crystallinity drops to zero before the density reaches that of the bulk vapor so, in some sense, the crystalline interior is separated from the bulk vapor by a liquid buffer. This is especially true if one considers that the structure is essentially liquidlike for crystallinities less than about 0.3. However, at the higher temperature, a more distinctive wetting of the cluster is seen whereby the crystallinity drops to zero over a region in which the density drops from a solidlike to a liquidlike value followed by a drop of the density from liquid to vapor values. This is analogous to the planar wetting illustrated above and shows that this is also a feature of the critical clusters near the triple point.
S  2.08  3.46  5.92  

1.00  1.005  1.01  
Model  
CNT  2063  1738  78  356  587  66  135  269  57 
m(1,1)  2053  1440  73  355  463  62  133  *  53 
m(2,1)  1899  1304  72  335  *  61  128  *  53 
0.32  0.74  0.41  0.78  1.24  0.46  1.27  1.78  0.51 
S  1.26  1.76  3.68  

0.975  0.98  0.99  
Model  
CNT  4045  7013  758  396  1142  407  59  189  184 
m(1,1)  4033  5171  708  392  783  379  58  *  171 
m(2,1)  3714  4176  704  368  *  377  57  *  170 
0.128  0.227  0.099  0.408  0.547  0.139  1.029  1.252  0.223 
v.3 Transient liquid state
Even in regions where double nucleation is not possible, the metastable liquid state could still play a role in nucleation. This becomes a question of the nucleation pathway which requires much more information than just the properties of the critical nucleus. A full description of the pathway would necessarily have to be dynamical in nature accounting for a variety of kinetic effects including heat and mass transport. It goes without saying that such a detailed description, while highly desirable, can be expected to be very difficult to implement.
Alternatively, one might imagine beginning at the transition state which is a saddle point of the free energy functional and is characterized by a size (in numbers of atoms), . It would be natural then to find other stationary points under the constraint of fixed number of atoms in the cluster and thus to trace a path from the critical cluster to the bulk vapor. This is in fact similar to the methods used in simulationstWF (). However, even in the case of the liquidvapor transition this gives discontinuous paths whereby the structure changes discontinuously for some value of . For the liquidvapor transition, this is a transition from a welldefined cluster with a liquidlike central density to a much larger structure with a central density slightly larger than the vapor. In the present case, the transition observed is from a welldefined crystalline cluster to one with zero crystallinity.
An alternative, widely used in quantum chemistry and in the study of clusters, is the construction of steepest descent pathways away from the transition state. These are expected to be approximations to the dynamical pathway, especially in the case of quasiequilibrium dynamics due to some sort of damping: an example might be the behavior of colloids (which can often be modeled as simple fluids) or macromolecules in solution. Many different methods are used to construct the steepest descent paths (also commonly called the Minimum Free Energy Pathway or MFEP) including heuristic methods such as the Nudged Elastic Band and the String Method. Both methods have been applied to nucleation problemsLutsko_JCP_2008 (); LutskoEPL () but here I expand on recent work which seeks to construct the exact MFEP for parameterized density profilesLutsko_InPrep ().
All methods of constructing steepest descent paths require the notion of a distance between two points in parameter space since ”steepest descent” literally means the path for which the energy varies most rapidly per unit distance in parameter space. The problem is that the various parameters used here  densities, widths, crystallinity  are incommensurate. However, since they exist only to specify a density profile, a natural solution is to define a metric in densityspace and then to use this to induce a metric in parameter space. In the study of fluids, the metric was taken to be the Euclidean distance between two density profiles,
(40) 
It would be natural to continue to use this definition however, it is unsuitable for two reasons. The first is simply that it is technically difficult to evaluate. The second, more importantly, is that it leads to the metric being sensitive to details of the lattice structure. For example, as the radius parameter varies there is little variation in the metric until the radius crosses an atomic shell at which point there is very rapid variation. This is unacceptable in the present context since the free energy surface constructed above is based on a separation of length scales according to which the order parameters vary slowly over atomic length scales.
In order to motivate a simple alternative more in keeping with the present approach, note that for the liquidvapor transition the crystallinity is identically zero, , so that the Euclidean distance function becomes
(41) 
which is the Euclidean distance between the amplitudes in the expansion of the density. In fact, the parameterization of the density used here is
(42) 
where the first sum is over the first (nearestneighbor) shell in wavevector space. Hence, the quantity is the spatiallyvarying amplitude of the smallest nonzero wavevector. It therefore seems reasonable to treat this on a par with the amplitude of the zerowavevector component and to define a distance function as
(43) 
which is what will be used henceforth. When the order parameters are expressed in terms of a collection of scalar parameters as and this then defines a distance between points in parameter space
(44) 
Assuming this function is sufficiently continuous, the distance function implies a metric
(45) 
The steepest descent paths are then determined by
(46) 
This equation is obviously similar to a dynamics consisting of simple relaxation driven by free energy gradients as discussed in Appendix C. However, since nucleation is, fundamentally, a noisedriven process, there is no reason to expect a priori that the correct path can be determined by running a deterministic dynamics backwards from the transition state. The idea of steepest descent is that it includes the idea of being a most probable path since it is in some sense the most efficient path over the barrier. (An analogy would be a multidimensional random walker in a potential field. In order for the walker to pass over a barrier, it must take some number of improbable steps in the right direction until it reaches the top of the barrier. The steepest descent path is the one involving the fewest number of steps. Other paths must cross the same barrier, but by including more steps, say in the “wrong” direction, the probability of falling backwards towards the local minimum increases.)
The nucleation pathways have been calculated using this model for values of the supersaturation such that the metastable liquid is more stable than the vapor, but below the doublenucleation threshold. Figure 10 shows the pathway for and plotted as a function of the the total number of atoms in the cluster. It is clear that both the central density and central crystallinity begin at the values of the bulk vapor. This is because the surface tension depends on the gradient of these quantities and for very small clusters, the dominance of the surface tension term over the bulk free energy contribution forces the gradients to be zero. As the cluster grows, the density increases very rapidly while the crystallinity increases more slowly. The core of the cluster therefore densifies more quickly than it crystallizes but the effect is minor. Figure 11 shows the same quantities for , close to the triple point, and for . In this case, transient twostep nucleation is clearly in evidence: the density increases while the crystallinity remains almost at zero until the cluster consists of well over 100 atoms. These contrasting behaviors are compared in Fig. 12 where the number of “crystalline” atoms, defined as
(47) 
is plotted as a function of . Near the triple point, the number of crystalline atoms does not increase appreciably until the cluster is over atoms in size whereas it increases almost immediately at the lower temperature. This behavior is, incidentally, quite similar to the observations of ten Wolde and Frenkel who noted that for their model of globular proteins, the favored nucleation path involved twostep nucleation near the triple point but that this was not seen further below the triple pointtWF ().
Vi Conclusions
The primary goal of this work has been they study of different pathways for the homogeneous nucleation of crystals from solution based on meanfield, DFT models. Simply from the bulk thermodynamics, it is possible to divide the phase diagram into regions for which double nucleation is and is not possible. When the bulk models are extended to inhomogeneous, multiphase systems  here via the squaredgradient approximation  it becomes possible to determine when double nucleation is more energetically favorable than singlestep nucleation. Finally, by studying steepestdescent pathways connecting the transition states to the bulk states it was possible to illustrate transient twostep nucleation for the LennardJones fluid. A summary of these investigations would be that:

there is a lower supersaturation limit for double nucleation: at the triple point, the limit goes to one and it increases rapidly as the temperature is lowered.

transient twostep nucleation seems to be closely tied to wetting and to therefore occur most distinctly near the triple point.
Almost all phases of this study could be improved. The underlying DFT model is highly simplified and in particular the inclusion of a longranged van der Waals term would give a more realistic description of wetting including a finite wettinglayer thickness. The modeling of the gradient coefficients is crude and quite empirical. One can estimate them directly from the DFT modelLowen1 (); Lowen2 (); LutskoGL (); OLW1 () and this might give a more realistic dependence on density and crystallinity. The piecewise continuous models could be used with more than the minimal number of links used here or some other, perhaps more physical, basis functions could be chosen (or the EulerLagrange equations could be solved directly without parameterizing the fields). Finally, instead of the Minimum Free Energy Pathways studied here, a more physically meaningful approach would be to search for the most probable pathways which requires the introduction of dynamics and noise, but which seems quite feasible based on the approach of Heymann and VandenEijndenMLP ().
Acknowledgements.
I am grateful to Gregoire Nicolis for reading an early draft of this manuscript and for a number of useful suggestions. This work was supported in part by the European Space Agency under contract number ESA AO2004070.Appendix A Bulk solid properties
The calculation of bulk thermodynamic properties for the hardsphere solid using the given parameterization is complex. The reason is that the hardsphere free energy functional diverges for slightly greater than  at typical solid densities, the divergence occurs for . It might seem that this is no problem as one simply takes however, upon closer inspection, this doesn’t work.
To explain the problem, let us consider the ”correct” way to fix these two densities. In principle, an equilibrium state must satisfy three conditions:
(48)  
Let the solution to these equations be . Solution of these three simultaneous equations is delicate due to the fact that (a) and (b) the surprising fact that
(49) 
In words, the free energy is a very rapidly varying function of for near .
This situation is somewhat better than it seems since in practical calculations, we are almost always performing some sort of search over, or tabulation in terms of, the chemical potential. Hence, rather than being given and having to solve three simultaneous equations, we can often take, e.g., as the independent parameter and search for satisfying
(50)  
We then calculate given these values. (Note however that even evaluating the various derivatives numerically is quite delicate when the range over which can be varied symmetrically might be or smaller.)
There is a heuristic which alleviates many of these problems. Intuitively, one tends to think in terms of the density and not to distinguish between the lattice density and the average density. In fact, in most calculations for the hardsphere solid, the approximation is made and the results appear quite reasonable. Since
(51) 
this suggests that it must be the case that
(52) 
In fact, this suspicion is borne out in practice. Combining these two points, and noting that the value of that makes the free energy stationary varies slowly as a function of the various densities, an efficient practical procedure is to choose a value of , to set , to determine the stationary value of and to use the left hand side of Eq.52 to evaluate . This involves a simple, controlled minimization and no need to evaluate the free energy near the divergence. A proof that Eq. 52 is exact, or at least a good approximation, is to my knowledge missing and would be useful.
Finally, an important point is that these technicalities play no role in the calculations presented here. Since we deal here with interfacial problems, the important thing is how the free energy varies as the average density varies from that of the solid down to a relatively low value (that of a liquid or vapor). Furthermore, in clusters, the interior density is always somewhat below that of the bulk solid (dwarfing the comparatively tiny difference between and ) . The only way these technical points would be important would be if we tried to do a free minimization of the free energy for a planar interface in which case the algorithm would have to find the correct values for the bulk system. This would essentially mean solving Eq. 50 which, because of the divergences very near the correct solution, leads to numerical challenges. However, using the piecewisecontinuous approximation, we always set the bulk values ”by hand” and thereby avoid this problem
Appendix B Formal results for gradient coefficients
They are given by
where is the direct correlation function for the bulk solid. Unfortunately, the DFT models used here do not give realistic expressions for this quantity although some useful results are possible. In particular, expanding in the crystallinity gives
(53)  
with
(54)  