Classical nucleation theory

The central result of classical nucleation theory is a prediction for the rate of nucleation ${\displaystyle R}$, in units of (number of events)/(volume·time). For instance, a rate ${\displaystyle R=1000\ {\text{m}}^{-3}{\text{s}}^{-1}}$ in a supersaturated vapor would correspond to an average of 1000 droplets nucleating in a volume of 1 cubic meter in 1 second.

The CNT prediction for ${\displaystyle R}$ is[3]

${\displaystyle R\ =\ N_{S}Zj\exp \left(-{\frac {\Delta G^{*}}{k_{B}T}}\right)}$

where

• ${\displaystyle \Delta G^{*}}$ is the free energy cost of the nucleus at the top of the nucleation barrier, and ${\displaystyle k_{B}T}$ is the average thermal energy with ${\displaystyle T}$ the absolute temperature and ${\displaystyle k_{B}}$ the Boltzmann constant.
• ${\displaystyle N_{S}}$ is the number of nucleation sites.
• ${\displaystyle j}$ is the rate at which molecules attach to the nucleus.
• ${\displaystyle Z}$ is the Zeldovich factor, (named after Yakov Zeldovich[6]) which gives the probability that a nucleus at the top of the barrier will go on to form the new phase, rather than dissolve.

This expression for the rate can be thought of as a product of two factors: the first, ${\displaystyle N_{S}\exp \left(-\Delta G^{*}/k_{B}T\right)}$, is the number of nucleation sites multiplied by the probability that a nucleus of critical size has grown around it. It can be interpreted as the average, instantaneous number of nuclei at the top of the nucleation barrier. Free energies and probabilities are closely related by definition.[7] The probability of a nucleus forming at a site is proportional to ${\displaystyle \exp[-\Delta G^{*}/kT]}$. So if ${\displaystyle \Delta G^{*}}$ is large and positive the probability of forming a nucleus is very low and nucleation will be slow. Then the average number will be much less than one, i.e., it is likely that at any given time none of the sites has a nucleus.

The second factor in the expression for the rate is the dynamic part, ${\displaystyle Zj}$. Here, ${\displaystyle j}$ expresses the rate of incoming matter and ${\displaystyle Z}$ is the probability that a nucleus of critical size (at the maximum of the energy barrier) will continue to grow and not dissolve. The Zeldovich factor is derived by assuming that the nuclei near the top of the barrier are effectively diffusing along the radial axis. By statistical fluctuations, a nucleus at the top of the barrier can grow diffusively into a larger nucleus that will grow into a new phase, or it can lose molecules and shrink back to nothing. The probability that a given nucleus goes forward is ${\displaystyle Z}$.

Taking into consideration kinetic theory and assuming that there is the same transition probability in each direction, it is known that ${\displaystyle x^{2}=2Dt}$. As ${\displaystyle Zj}$ determines the hopping rate, the previous formula can be rewritten in terms of the mean free path and the mean free time ${\displaystyle \lambda ^{2}=2D\tau }$. Consequently, a relation of ${\displaystyle Zj}$ in terms of the diffusion coefficient is obtained

${\displaystyle Zj={\frac {1}{\tau }}={\frac {2D}{\lambda ^{2}}}}$.

Further considerations can be made in order to study temperature dependence. Therefore, Einstein-Stokes relation is introduced under the consideration of a spherical shape

${\displaystyle D={\frac {k_{B}T}{6\pi \eta \lambda }}}$, where ${\displaystyle \eta }$ is the material's viscosity.

Considering the last two expressions, it is seen that ${\displaystyle Zj}$ ${\displaystyle \propto T}$. If ${\displaystyle T\approx T_{m}}$, being ${\displaystyle T_{m}}$ the melting temperature, the ensemble gains high velocity and makes ${\displaystyle Zj}$ and ${\displaystyle \Delta G^{*}}$ to increase and hence, ${\displaystyle R}$ decreases. If ${\displaystyle T\ll T_{m}}$, the ensemble has a low mobility, which makes ${\displaystyle R}$ to decrease as well.

To see how this works in practice we can look at an example. Sanz and coworkers[8] have used computer simulation to estimate all the quantities in the above equation, for the nucleation of ice in liquid water. They did this for a simple but approximate model of water called TIP4P/2005. At a supercooling of 19.5 °C, i.e., 19.5 °C below the freezing point of water in their model, they estimate a free energy barrier to nucleation of ice of ${\displaystyle \Delta G^{*}=275k_{B}T}$. They also estimate a rate of addition of water molecules to an ice nucleus near the top of the barrier of ${\displaystyle j=10^{11}\ {\text{s}}^{-1}}$ and a Zeldovich factor ${\displaystyle Z=10^{-3}}$. The number of water molecules in 1 m³ of water is approximately 10²⁸. These leads to the prediction ${\displaystyle R=10^{-83}\ {\text{s}}^{-1}}$, which means that on average one would have to wait 10⁸³s (10⁷⁶ years) to see a single ice nucleus forming in 1 m³ of water at -20 °C!

This is a rate of homogeneous nucleation estimated for a model of water, not real water — in experiments one cannot grow nuclei of water and so cannot directly determine the values of the barrier ${\displaystyle \Delta G^{*}}$, or the dynamic parameters such as ${\displaystyle j}$, for real water. However, it may be that indeed the homogeneous nucleation of ice at temperatures near -20 °C and above is extremely slow and so that whenever water freezes at temperatures of -20 °C and above this is due to heterogeneous nucleation, i.e., the ice nucleates in contact with a surface.

The classical nucleation theory hypothesis for the form of ${\displaystyle \Delta G}$ can be examined more rigorously using the tools of statistical mechanics.[13] Specifically, the system is modeled as a gas of non-interacting clusters in the grand canonical ensemble. A state of metastable equilibrium is assumed, such that the methods of statistical mechanics hold at least approximately.[14] The grand partition function is[15]

${\displaystyle Q=\sum _{N=0}^{\infty }z^{N}\sum _{\mu _{S}|N}e^{-\beta {\mathcal {H}}_{N}(\mu _{S})},\qquad z\equiv e^{\beta \mu }.}$

Here the inner summation is over all microstates ${\displaystyle \mu _{S}}$ which contain exactly ${\displaystyle N}$ particles. It can be decomposed into contributions from each possible combination of clusters which results in ${\displaystyle N}$ total particles.[16] For instance,

${\displaystyle \sum _{\mu _{S}|N=3}e^{-\beta {\mathcal {H}}_{3}(\mu _{S})}=q_{3}+q_{2}q_{1}+{\frac {1}{3!}}q_{1}^{3}}$

where ${\displaystyle q_{n}}$ is the configuration integral of a cluster with ${\displaystyle n}$ particles and potential energy ${\displaystyle U_{n}\left(\{\mathbf {r_{n}} \}\right)}$:

${\displaystyle {\begin{aligned}q_{n}={\frac {1}{n!}}{\frac {1}{\Lambda ^{3n}}}\int d\mathbf {r} ^{n}e^{-\beta U_{n}\left(\{\mathbf {r_{n}} \}\right)}.\end{aligned}}}$

The quantity ${\displaystyle \Lambda }$ is the thermal de Broglie wavelength of the particle, which enters due to the integration over the ${\displaystyle 3n}$ momentum degrees of freedom. The inverse factorials are included to compensate for overcounting, since particles and clusters alike are assumed indistinguishable.

More compactly,

${\displaystyle Q=\exp \left[\sum _{n=1}^{\infty }q_{n}z^{n}\right]}$.

Then, by expanding ${\displaystyle Q}$ in powers of ${\displaystyle q_{m}z^{m}}$, the probability ${\displaystyle P(m,\ell )}$ of finding exactly ${\displaystyle \ell }$ clusters which each has ${\displaystyle m}$ particles is

${\displaystyle {\begin{aligned}P(m,\ell )&={\frac {z^{\ell m}q_{m}^{\ell }}{\ell$ !}}\cdot {\frac {\exp \left(\sum _{n\neq m}^{\infty }q_{n}z^{n}\right)}{\exp \left(\sum _{n=1}^{\infty }q_{n}z^{n}\right)}}\\&={\frac {z^{\ell m}q_{m}^{\ell }}{\ell !}}\exp \left(-q_{m}z^{m}\right).\end{aligned}}}

The number density ${\displaystyle \rho _{m}}$ of ${\displaystyle m}$-clusters can therefore be calculated as

${\displaystyle {\begin{aligned}\rho _{m}&={\frac {1}{V}}\sum _{\ell =0}^{\infty }\ell {\frac {z^{\ell m}q_{m}^{\ell }}{\ell$ !}}\exp \left(-q_{m}z^{m}\right)\\&={\frac {z^{m}q_{m}}{V}}.\end{aligned}}}

This is also called the cluster size distribution.

The grand potential ${\displaystyle \Omega }$ is equal to ${\displaystyle -k_{B}T\ln Q}$, which, using the thermodynamic relationship ${\displaystyle \Omega =-pV}$, leads to the following expansion for the pressure:

${\displaystyle {\frac {p}{k_{B}T}}=\sum _{n=1}^{\infty }{\frac {q_{n}}{V}}z^{n}.}$

If one defines the right hand side of the above equation as the function ${\displaystyle \Pi (T,z)}$, then various other thermodynamic quantities can be calculated in terms of derivatives of ${\displaystyle \Pi }$ with respect to ${\displaystyle z}$.[17]

The connection with the simple version of the theory is made by assuming perfectly spherical clusters, in which case ${\displaystyle U_{n}\left(\{\mathbf {r_{n}} \}\right)}$ depends only on ${\displaystyle n}$, in the form

${\displaystyle U_{n}=-E_{0}n+wn^{1/2}}$

where ${\displaystyle E_{0}}$ is the binding energy of a single particle in the interior of a cluster, and ${\displaystyle w>0}$ is the excess energy per unit area of the cluster surface. Then, ${\displaystyle q_{n}\propto \exp(-\beta \cdot (-E_{0}n+wn^{1/2}))}$, and the cluster size distribution is

${\displaystyle \rho _{n}\propto e^{-\beta [-(E_{0}+\mu )n+wn^{1/2}]},}$

which implies an effective free energy landscape ${\displaystyle \Delta G(n)=-(E_{0}+\mu )n+wn^{1/2}}$, in agreement with the form proposed by the simple theory.

On the other hand, this derivation reveals the significant approximation in assuming spherical clusters with ${\displaystyle U_{n}\left(\{\mathbf {r_{n}} \}\right)=-E_{0}n+wn^{1/2}}$. In reality, the configuration integral ${\displaystyle q_{n}}$ contains contributions from the full set of particle coordinates ${\displaystyle \{\mathbf {r_{n}} \}}$, thus including deviations from spherical shape as well as cluster degrees of freedom such as translation, vibration, and rotation. Various attempts have been made to include these effects in the calculation of ${\displaystyle q_{n}}$, although the interpretation and application of these extended theories has been debated.[5][18][19] A common feature is the addition of a logarithmic correction ${\displaystyle \sim \ln n}$ to ${\displaystyle \Delta G(n)}$, which plays an important role near the critical point of the fluid.[20]

Classical nucleation theory makes a number of assumptions which limit its applicability. Most fundamentally, in the so-called capillarity approximation it treats the nucleus interior as a bulk, incompressible fluid and ascribes to the nucleus surface the macroscopic interfacial tension ${\displaystyle \sigma }$, even though it is not obvious that such macroscopic equilibrium properties apply to a typical nucleus of, say, 50 molecules across.[21][22] In fact, it has been shown that the effective surface tension of small droplets is smaller than that of the bulk liquid.[23]

In addition, the classical theory places restrictions on the kinetic pathways by which nucleation occurs, assuming clusters grow or shrink only by single particle adsorption/emission. In reality, merging and fragmentation of entire clusters cannot be excluded as important kinetic pathways in some systems. Particularly in dense systems or near the critical point – where clusters acquire an extended and ramified structure – such kinetic pathways are expected to contribute significantly.[23] The behavior near the critical point also suggests the inadequacy, at least in some cases, of treating clusters as purely spherical.[24]

Various attempts have been made to remedy these limitations and others by explicitly accounting for the microscopic properties of clusters. However, the validity of such extended models is debated. One difficulty is the exquisite sensitivity of the nucleation rate ${\displaystyle R}$ to the free energy ${\displaystyle \Delta G^{*}}$: even small discrepancies in the microscopic parameters can lead to enormous changes in the predicted nucleation rate. This fact makes first-principles predictions nearly impossible. Instead, models must be fit directly to experimental data, which limits the ability to test their fundamental validity.[25]

For simple model systems, modern computers are powerful enough to calculate exact nucleation rates numerically. An example is the nucleation of the crystal phase in a system of hard spheres, which is a simple model of colloids consisting of perfectly hard spheres in thermal motion. The agreement of CNT with the simulated rates for this system confirms that the classical theory is a reasonable approximation.[26] For simple models CNT works quite well; however it is unclear if it describes complex (e.g., molecular) systems equally well. For example, in the context of vapor to liquid nucleation, the CNT predictions for the nucleation rate are incorrect by several orders of magnitude on an absolute scale — that is, without renormalizing with respect to experimental data. Nevertheless, certain variations on the classical theory have been claimed to represent the temperature dependence adequately, even if the absolute magnitude is inaccurate.[27] Jones et al. computationally explored the nucleation of small water clusters using a classical water model. It was found that CNT could describe the nucleation of clusters of 8-50 water molecules well, but failed to describe smaller clusters.[28] Corrections to CNT, obtained from higher accuracy methods such as quantum chemical calculations, may improve the agreement with experiment.[29]