On the role of metastable intermediate states in the homogeneous nucleation of solids from solution1footnote 11footnote 1To appear in Advances in Chemical Physics

On the role of metastable intermediate states in the homogeneous nucleation of solids from solution111To appear in Advances in Chemical Physics

James F. Lutsko Center for Nonlinear Phenomena and Complex Systems CP 231, Université Libre de Bruxelles, Blvd. du Triomphe, 1050 Brussels, Belgium jlutsko@ulb.ac.be http://www.lutsko.com
July 14, 2019

The role of metastable liquid phases in vapor-crystal 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 hard-sphere functional (Fundamental Measure Theory) and a liquid-like contribution to account for the attractive part of the interactions. The model gives a semi-quantitatively accurate description of both the vapor-liquid-solid phase diagram for both simple fluids (Lennard-Jones interactions) and of the low-density/high-density/crystal phase diagram for model globular proteins (ten Wolde-Frenkel 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 squared-derivative terms in these order parameters to account for inhomogeneities thus producing a model similar in spirit to phase-field theory. It is shown that for both interaction models, the vapor-crystal part of the phase-diagram 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 vapor-solid and vapor-liquid 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 vapor-solid transition involves a single nucleation event but the sub-critical clusters tend to be liquid-like while small, is shown to also be possible when the metastable liquid exists, but the supersaturation is too small for double nucleation.

64.60.Q-, 64.70.Hz, 64.60.My, 82.60.Nh, 68.08.-p

I Introduction

The goal of this work is to describe, at a semi-microscopic 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 single-component 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 liquid-like phase which simply corresponds to the liquid phase of the effective single-component system. Thus, the following work is intended to be applicable both to nucleation in single-component 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 “low-concentration” or “low density” phase of molecules, such as proteins, in solution and statements concering the “liquid phase” are equally applicable to the “high-concentration” or “high-density” 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 (liquid-crystal) 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 vapor-crystal and vapor-liquid transition. Alternatively, even if direct vapor-crystal nucleation is favored, it is still possible that the growing cluster nevertheless passes through a liquid-like stage which, however, is always sub-critical, as a transient on the way to creating a solid-like critical cluster. In this case, small clusters would be liquid-like in structure, becoming more crystalline as they grow larger. This scenario, involving only a single nucleation event, will be referred to as ”transient two-step nucleation”. The simulations of ten Wolde and FrenkeltWF () are clear illustrations of transient two-step nucleation: they exhibit a nucleation pathway that involves liquid-like 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 Lennard-Jones 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 vapor-liquid 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 non-bulk 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 vapor-solid 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 vapor-liquid transition involves two different densities but occurs at zero crystallinity while the liquid-solid transition occurs at (nearly constant) density and involves a change from zero to finite crystallinity. In the transient two-step 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 ”two-step” 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, finite-temperature Density Functional Theory (DFT). If DFT were sufficiently well-developed 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 liquid-vapor 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 fluid-solid systems for hard-sphere 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 hard-sphere diameter and the representation of the free energy as the sum of the hard-sphere functional and a second term accounting for the attractive part of the interactions. Futherrmore, even for the hard-sphere 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 phase-field theories commonly used to study solid-solid 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 semi-empirical 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 semi-quantitatively accurate bulk phase diagrams as well as of the liquid-vapor 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 energy-surface techniques commonly applied to the study of chemical reaction pathways and structural transitionsWales (). A recently developed method - involving the approximation of the spatially-varying order parameters as piecewise-continuous 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,


where is a lattice vector, controls the width, is the occupancy and is the lattice density. The average density is


The density can also be written in Fourier representation as


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


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


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 non-zero 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 liquid-vapor 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 grand-canonical free energy is written as


where is the chemical potential. In general, if the density can be written in terms of order parameters, as


so that the density of the uniform, bulk system is


and if in some sense ”slowly varying”, then the squared-gradient approximation (SGA) for the functional is


where the summation convention is used. The first term of the free energy involves the bulk free energy density defined as


For the order parameters used here, the free energy is explicitly


ii.3 Bulk Free Energy

The bulk free energy model used here is based on the idea of separating the free energy into a hard-sphere contribution, for which the DFT is well developed, and a second contribution that accounts for the long-ranged 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 hard-sphere model is independent of the local structureCurtin (); LN () giving


where is the effective hard-sphere diameter and


is the contribution of the attractive part of the interaction to the liquid free energy. The effective hard-sphere 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 grand-canonical 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,


which is to say


For parameterized profiles, the requirement that the free energy be a minimum gives


and if the parameters are constants, then


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


In particular, using the chain rule for functional differentiation,


gives the usual thermodynamic relation


Thus, for uniform densities, and e.g. , the conditions for coexistence are


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


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 grand-canonical 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 grand-canonical 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 Lennard-Jones potential,


and globular proteins as described by the ten Wolde-Frenkel model interaction,


which will be studied here for the value as is appropriate to model the phase behavior of globular proteins.

Both of these potentials are long-ranged in the sense of decaying as power-laws. In simulation, infinite-ranged potentials are difficult to deal with, so any long-ranged potential is typically cutoff at some distance, . For Monte-Carlo, a simple shift of the potential to make the energy continuous at the cutoff is typically used so that the so-called truncated and shifted potential is


In molecular dynamics simulations, the force is usually required to be continuous so that the force-shifted potential is commonly used,


where . Other forms are also important, particularly the Broughton-Gilmer modification of the LJ potential:


where , , , and BG1 ().

All of these details significantly affect the liquid-vapor 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 mean-field corrections. Although useful for studying liquid-vapor coexistenceLutsko_JCP_2008 (); Lutsko_JCP_2008_2 (), they are of limited utility for vapor-solid nucleation since the interesting affects occur below the triple point. Less accurate, but more broadly applicable, are mean-field equations of state based on thermodynamic perturbation theory such as the Weeks-Chandler-Andersen (WCA) theoryWCA1 (); WCA2 (); WCA3 (); Ree1 (); Ree2 () which will be used here.

iii.3 Solid phase diagrams and intermediate phases

Figure 1: Phase diagrams for the Lennard-Jones potential, left panel, and the tWF potential, right panel, together with simulation data. The tWF potential was cutoff and shifted at .

Figure 1 shows the calculated vapor-liquid-FCC solid phase diagram for the infinite-ranged 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 liquid-vapor critical point and a triple point whereas for the protein model, the liquid-vapor 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 semi-quantitative agreement with simulation for all three phases.

iii.3.1 Nucleation scenarios for globular proteins

Figure 2: The phase diagram for the tWF potential as in Fig. 1 but without the simulation data and showing the spinodal for the vapor-liquid transition (broken line). The figure also shows, in red, the boundaries for the existence of the intermediate liquid phase: there is no liquid phase for chemical potentials corresponding to vapor densities to the left of the red line. A similar line for the solid phase is nearly indistinguishable from the solid binodal. Similar boundaries also exist on the Lennard-Jones phase diagram but are so close to the binodals as to be almost indistinguishable. The horizontal line, an isotherm, picks out coexisting states with the coexisting vapor and solid states being labeled A and B, respectively.

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 solid-branch 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 vapor-liquid 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 liquid-vapor binodal. Eventually, one reaches the vapor-liquid 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 lower-density region, for which a vapor with the same chemical potential can always be found, and a higher-density region with no corresponding vapor.

What is the role of the metastable liquid phase in vapor-solid 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 vapor-liquid 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 liquid-branch 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 liquid-vapor 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 vapor-solid 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 vapor-branch of the solid-vapor 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 vapor-solid coexistence curves, and moving to the right we therefore find that:

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

  2. From the -nodal line, to the vapor branch of the vapor-liquid 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 two-step nucleation.

  3. From the vapor branch of the liquid-vapor coexistence curve to the vapor-liquid 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.

Figure 3: Phase diagrams for the tWF potential with (vapor-solid) supersaturation as the independent variable.

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 vapor-solid 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 liquid-like state, i.e. transient two-step 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 two-step nucleation occurs.

iii.3.2 Nucleation scenarios for simple fluids

Figure 4: Phase diagram for the LJ potential with (vapor-solid) supersaturation as the independent variable.

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 liquid-vapor coexistence below the triple point. However, just as in the case of the globular proteins, the vapor branch of the vapor-liquid coexistence curve is to the right of the vapor branch of the vapor-solid 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 liquid-vapor 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 liquid-solid coexistence curves lies to the left of the liquid branch of the liquid-vapor coexistence curve: i.e. it is in the metastable (or even unstable) region of the liquid-vapor 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 vapor-solid transition is usually drawn below the triple point, but there is a metastable liquid and a liquid-vapor coexistence curve in this region. It is this metastable liquid phase which gives rise to the possibility of two-step nucleation even for simple liquids. However, the fact that the vapor branches of the vapor-solid coexistence curves and the vapor-liquid 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 double-nucleation 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,


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) vapor-solid cluster which nevertheless captures the important physics of the real system:


This piecewise-linear 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




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


and the radius from




The same model applied to the vapor-liquid and liquid-solid critical clusters gives


Finally, let us make the further approximation that the density of the liquid and solid states are almost the same so that


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:

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

  2. are not small compared to .

Clearly, the factor occurring in the denominator of can be important in raising the vapor-liquid barrier compared to the vapor-solid barrier. On the other hand, the factor involving the gradient coefficients is always going to be smaller for the vapor-liquid cluster than for the vapor-solid 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 second-gradient 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 liquid-solid or vapor-solid planar surface tension simply because the squard-gradient 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 slowly-varying order parameters that underlies the SGA is unlikely to be very good, although for the particular case of the liquid-vapor 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



Second, there is evidence that the density-density coefficient can be well modeled in the liquid, i.e. for , by a density-independent 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 short-range repulsion of the pair potential so that the liquid-solid surface tension can be approximated by that of the hard-sphere solid evaluated with an effective hard-sphere diameter giving the useful approximation that the excess surface free energy for a planar liquid-solid interface is


. This suggests that and should be dominated by hard-sphere 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 hard-sphere part of the free energy. All of this suggests a simple approximation for the structural coefficients along the lines of


In a final simplification, the low-density limit gives .

Figure 5: Excess surface free energy for the solid-vapor, liquid-vapor and liquid-solid interfaces for a system interacting via the LJBG pair potential. The small symbols are simulation data taken from Broughton and GilmerBG4 () (solid-vapor and liquid-vapor) and from Davidchack and LairdDavidchack_crystal_melt ()(solid-liquid). The diamonds show the Laird approximation for the solid-liquid planar surface tensionLaird_solid_liquid_hs (). The simulation results at each temperature are, from highest to lowest, for the 111, 100 and 110 planes respectively.

iv.2 Planar interfaces

Figure 5 shows the result of using this model to calculate the liquid-solid, vapor-solid and vapor-liquid 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 liquid-solid surface tension in agreement with a single point of either simulation or the Davidchack-Laird approximation is enough to give a reasonable description of the liquid-solid and vapor-solid 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: The order parameters (upper panels) and the 100 planar-averaged density (lower panels) for the solid-fluid interface at three temperatures which are, respectively, below, near and above the triple point. In the upper panels, the higher curve is and the lower curve is . All quantities are in reduced units.

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 solid-liquid transition while the second has the structure of a liquid-vapor 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 long-ranged 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 Vapor-Crystal 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 super-critical. As mentioned above, the most important issue that arises in vapor-solid 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 near-bulk solid properties of the interior of the cluster. The second pathway involves double nucleation. First, a purely liquid-like 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 bulk-solid, 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 two-step nucleation and is in some sense intermediate between the classical and double-nucleation scenarios. Sufficiently small clusters, which are of course sub-critical, are liquid-like 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 solid-like than liquid-like, 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, noise-driven 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 eigenvector-following 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 vapor-liquid transition, another to the vapor-solid transition and a third to the liquid-solid 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 two-step transition via first a vapor-liquid transition and then a liquid-solid transition. It is assumed that whichever transition involves the lowest free energy barriers will dominate.

Figure 7: Phase diagram for the LJ potential with (vapor-solid) supersaturation as the independent variable. Double nucleation occurs in the region to the right and above the broken line: between the broken line and red line, the liquid is metastable but double nucleation has a higher energy barrier than single-step nucleation of the solid.

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 vapor-solid transition is lower than that for the vapor-liquid transition. However, at larger supersaturations, the vapor-liquid transition becomes less costly thus implying that the double-nucleation scenario is favored.

Figure 8: The structure of the critical nucleus for and as determined using the CNT, m(1,1) and m(2,1) parameterizations of the profiles. In each figure, the upper curve is the average density and the lower curve is the crystallinity. In the figure for the m(2,1) model, a dashed line is drawn at the radius where the crystallinity becomes zero: the fact that the second link of the density profile begins at nearly this point is a clear illustration of the role of the intermediate liquid state in wetting the cluster.
Figure 9: The same as Fig. 8 for and . Desipite the existence of a metastable intermediate liquid, there is no significant wetting behaviour in this case.

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 liquid-like 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 solid-like to a liquid-like 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
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
Table 1: Properties of critical nuclei for the LJBG interaction as a function of supersaturation at . The supersaturation and lattice density are given in the first two rows followed by the excess free energy for the nuclei for vapor-liquid (VL), vapor-solid (VS) and liquid-solid(LS) nucleation. The last row gives the difference in bulk free energies for the various transitions. Blank entries occur for cases where it was not possible to find a transition state.
S 1.26 1.76 3.68
0.975 0.98 0.99
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
Table 2: Same as Table 1 for .

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 liquid-vapor transition this gives discontinuous paths whereby the structure changes discontinuously for some value of . For the liquid-vapor transition, this is a transition from a well-defined cluster with a liquid-like 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 well-defined 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 quasi-equilibrium 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 density-space 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,


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 liquid-vapor transition the crystallinity is identically zero, , so that the Euclidean distance function becomes


which is the Euclidean distance between the amplitudes in the expansion of the density. In fact, the parameterization of the density used here is


where the first sum is over the first (nearest-neighbor) shell in wave-vector space. Hence, the quantity is the spatially-varying amplitude of the smallest non-zero wavevector. It therefore seems reasonable to treat this on a par with the amplitude of the zero-wavevector component and to define a distance function as


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


Assuming this function is sufficiently continuous, the distance function implies a metric


The steepest descent paths are then determined by


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 noise-driven 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 double-nucleation 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 two-step 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


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 two-step nucleation near the triple point but that this was not seen further below the triple pointtWF ().

Figure 10: The average density (upper curve) and crystallinity (lower curve) at the center of the cluster as a function of cluster size for . The inset shows an expanded view of the early stages of nucleation.
Figure 11: Same as Fig. 10 for .
Figure 12: The number of crystalline atoms as a function of the total number of atoms in the cluster showing transient two-step nucleation near the triple point, , and almost completely one-step nucleation away from the triple point.

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 mean-field, 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 squared-gradient approximation - it becomes possible to determine when double nucleation is more energetically favorable than single-step nucleation. Finally, by studying steepest-descent pathways connecting the transition states to the bulk states it was possible to illustrate transient two-step nucleation for the Lennard-Jones fluid. A summary of these investigations would be that:

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

  2. transient two-step 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 long-ranged van der Waals term would give a more realistic description of wetting including a finite wetting-layer 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 Euler-Lagrange 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 Vanden-EijndenMLP ().

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 AO-2004-070.

Appendix A Bulk solid properties

The calculation of bulk thermodynamic properties for the hard-sphere solid using the given parameterization is complex. The reason is that the hard-sphere 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:


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


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


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 hard-sphere solid, the approximation is made and the results appear quite reasonable. Since


this suggests that it must be the case that


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 piecewise-continuous 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