Free Energy Barriers for Crystal Nucleation from Fluid Phases

TL;DR

This paper introduces a simulation-based method to measure crystal nucleation energy barriers from fluid phases without assuming nucleus shape, demonstrated on colloid-polymer mixtures, revealing significant underestimation by spherical approximations.

Contribution

The study develops a novel finite-volume approach to determine surface free energies of nuclei directly from simulations, avoiding shape assumptions and enabling accurate barrier estimations.

Findings

Nucleation barriers scale with nucleus volume as predicted by classical theory.

Spherical nucleus approximation often underestimates nucleation barriers.

Effective surface tension can be precisely estimated, accounting for nonspherical shapes.

Abstract

Monte Carlo simulations of crystal nuclei coexisting with the fluid phase in thermal equilibrium in finite volumes are presented and analyzed, for fluid densities from dense melts to the vapor. Generalizing the lever-rule for two-phase coexistence in the canonical ensemble to finite volume, "measurements" of the nucleus volume together with the pressure and chemical potential of the surrounding fluid allows to extract the surface free energy of the nucleus. Neither the knowledge of the (in general non-spherical) nucleus shape nor of the angle-dependent interface tension is required for this task. The feasibility of the approach is demonstrated for a variant of the Asakura-Oosawa model for colloid-polymer mixtures, which form face-centered cubic colloidal crystals. For a polymer to colloid size ratio of $0.15$, the colloid packing fraction in the fluid phase can be varied from melt values to zero by the variation of an effective attractive potential between the colloids. It is found that the approximation of spherical crystal nuclei often underestimates actual nucleation barriers significantly. Nucleation barriers are found to scale as $\Delta F^*=(4\pi/3)^{1/3}\bar{\gamma}(V^*)^{2/3}+const.$ with the nucleus volume $V^*$, and the effective surface tension $\bar{\gamma}$ that accounts implicitly for the nonspherical shape can be precisely estimated.

Equations

$$\Delta F^{\ast }(V_n)=-(p_c-p_l)V_n+F_{surf}(V_n).$$

$$F_{surf}(V_n)=4\pi R_n^2\bar{\gamma }=A_{iso}\bar{\gamma }{V}_n^{2/3},A_{iso}=(36\pi {)}^{1/3}.$$

$$F_{surf}(V_n)=V_n^{2/3}{\int }_{A_W}\gamma (\vec{n})d\vec{s}\equiv A_W\bar{\gamma }{V}_n^{2/3},$$

$$\bar{\gamma }=(A_W{)}^{-1}{\int }_{A_W}\gamma (\vec{n})d\vec{s}$$

$$V_n^{{\ast }^{1/3}}=\frac{2A_W\bar{\gamma }}{3(p_c-p_l)}$$

$$\Delta F^{\ast }=\frac{1}{3}{A}_W\bar{\gamma }{V}_n^{{\ast }^{2/3}}=\frac{1}{2}(p_c-p_l)V_n^{\ast }.$$

$$J\propto exp(-\Delta F^{\ast }/k_{B}T),$$

$$\rho V={\rho }_m{V}_n+{\rho }_f(V-V_n).$$

$$p(\rho ,T)=p_{coex}=p_l({\rho }_f,T)=p_c({\rho }_m,T),$$

$$\mu (\rho ,T)={\mu }_{coex}={\mu }_l({\rho }_f,T)={\mu }_c({\rho }_m,T).$$

$$\mu (\rho ,T)={\mu }_l({\rho }_l,T)={\mu }_c({\rho }_c,T)>{\mu }_{coex},$$

$$p_c({\rho }_c,T)>p_l({\rho }_l,T)>p_{coex}.$$

$$\rho V={\rho }_c{V}_n+{\rho }_l(V-V_n).$$

$$U_{\text{eff}}({\sigma }_c<r$$

$$\frac{\pi }{6}{k}_{B}T{\sigma }_p^3\frac{z_p}{(1+\frac{1}{q}{)}^3}$$

$$U_{\text{rep}}($$

$$-$$

$${\eta }_p^r=\frac{\pi {\sigma }_p^3}{6}{e}^{\frac{{\mu }_p}{k_{B}T}},$$

$$\eta =\frac{\pi }{6}{\sigma }_c^3\frac{N}{V}.$$

$${\mu }_c(p_c)={\mu }_{coex}+{\int }_{p_{coex}}^{p_c}d{p}^{\prime }\frac{\pi }{6{\eta }_c(p^{\prime })},$$

$${\mu }_c(p_c)={\mu }_{coex}+\frac{\pi }{6{\eta }_m}(p_c-p_{coex}),$$

$$p_l\cong p_{coex}+\frac{6}{\pi }{\eta }_f({\mu }_l-{\mu }_{coex}),$$

$${\mu }_l={\mu }_{coex}+\frac{\pi }{9}\frac{A_W\bar{\gamma }}{{\eta }_m-{\eta }_f}(V_n^{\ast }{)}^{-1/3},V_n^{\ast }\to \infty .$$

$$q_{lm}(i)=\frac{1}{N(i)}\sum _{j=1}^{N(i)}{Y}_{lm}\left(\theta ({\vec{r}}_{ij}),\varphi ({\vec{r}}_{ij})\right),$$

$$\overline{q_{lm}}(i)=\frac{1}{N(i)+1}\sum _{k_i=0}^{N(i)}{q}_{lm}(k_i),$$

$$\overline{q_l}(i)={\left(\frac{4\pi }{(2l+1)}\sum _{m=-l}^{+l}|\overline{q_{lm}}(i){|}^2\right)}^{1/2},$$

$$U_{\text{exc}}=-\frac{\sqrt{3}}{r_{min}^2}U(r_{min})$$

Full text

I Introduction and Overview

Understanding nucleation phenomena in condensed matter is a fundamental problem of statistical mechanics and has important applications, for atmospheric sciences (formation of raindrops or snow crystals), geosciences (crystallization of silicate melts), and material science (e.g. microstructure of metallic alloys obtained by cooling their melts), for instance 1 ; 2 ; 3 ; 4 ; 5 ; 6 . The classical theory assumes that nanoscopic nuclei of the crystalline solid form spontaneously by thermal fluctuations when one suddenly has brought the system from the region where the fluid phase is the thermodynamically stable phase beyond the fluid-solid coexistence curve, so that the fluid phase is only metastable (“homogeneous nucleation”). Although in the metastable region the bulk free energy density of the crystal is lower than the free energy density of the metastable fluid, large-scale phase transformation from fluid to solid is often prevented by the cost in surface excess free energy caused by the fluid-solid interface. This means a nucleation event amounts to overcoming a free energy barrier. The classical theory of homogeneous nucleation estimates this barrier assuming that the formation free energy of the nucleus can be written in terms of macroscopic nucleus volume ($V_{n}$) and surface terms,

[TABLE]

Here $p_{c}$ is the pressure inside of the crystal nucleus, and $p_{l}$ the pressure in the metastable liquid phase surrounding it. The traditional basic assumption of the theory then neglects the fact that for a crystal-fluid interface the interfacial tension $\gamma(\vec{n})$ depends on the orientation of the interface normal ($\vec{n}$) relative to the crystal axes. If one could replace $\gamma(\vec{n})$ by a constant $\bar{\gamma}$, assuming that the nucleus has a spherical shape, we could write ($R_{n}$ is the radius of the nucleus)

[TABLE]

When the anisotropy of the interfacial tension is taken into account, however, for large $V_{n}$ the equilibrium shape of a crystal nucleus clearly is not a sphere, but rather given by the Wulff construction 7 ; 8 ; 9 ; 10 ; 11 . In fact, at low temperatures (T lower than the appropriate roughening transition temperatures $T_{R}$ 12 ; 13 ) crystal shapes with strictly planar facets occur. As an extreme example, Fig. 1 shows crystal shapes of face-centered cubic crystals that result 14 when only three interface orientations (111), (100), and (110) would be permitted (recall that an (111) interface in the ideal lattice structure is a closed packed plane with triangular lattice structure, and it is natural to expect that the interface tension $\gamma_{111}$ of this surface is smaller than $\gamma_{100}$ or $\gamma_{110}$).

While shapes as shown in Fig. 1 are expected as $T\rightarrow 0$, in general one expects to have facets with rounded corners and edges (for nonzero $T<T_{R}$) and nonspherical but curved shapes for $T>T_{R}$. An example for the gradual change of a crystal shape with increasing temperature has been presented for the nearest neighbor simple cubic Ising model recently 15 : there occurs a gradual change from a cube (such as shown in the rightmost figure in the top row of Fig. 1) at $T=0$ to a sphere which only occurs for temperatures very close to the critical point of this lattice model. Of course, this model completely lacks an off-lattice fluid phase, and hence this model cannot capture the essential features of the crystal/fluid transition.

When the interface tension is anisotropic, one can in principle generalize Eq. (2) by writing the surface free energy $F_{surf}(V_{n})$ as an integral over the crystal surface,

[TABLE]

where $A_{W}$ is the surface area of a crystal having the shape resulting from the unrestricted Wulff construction 7 ; 8 ; 9 ; 10 ; 11 and unit volume. Note that the geometric factor $A_{W}$ must always exceed $A_{iso}$ (e.g., for the cube resulting in the Ising model $A_{W}=6$, $A_{W}/A_{iso}=(6/\pi)^{1/3}$), and the average interfacial tension

[TABLE]

in general will exceed the smallest of the interfacial tensions, when more than one type of interface contributes.

Now, according to the classical theory of homogeneous nucleation, one simply combines Eqs. (1) and (2) or (3) to identify the nucleation barrier from the maximum of $\Delta F(V_{n})$. This yields the nucleus linear dimension at the barrier

[TABLE]

and the associated free energy barrier

[TABLE]

The nucleation rate then is assumed to be proportional to

[TABLE]

where nonexponential prefactors occur describing both kinetic effects and corrections due to various fluctuation effects. The nucleation rate is defined as the number of critical nuclei (with volume $V_{n}^{*}$) formed per unit volume and unit time in the limting case where the metastable state has had enough time to reach at least local equilibrium, and nucleation events are seldom enough so the change of this state due to the forming nuclei can be neglected, and nucleation can be approximated as a steady-state process 1 ; 2 ; 3 ; 4 ; 5 ; 16 ; 17 ; 18 ; 19 .

Fluctuation effects are important, since for the case of barriers that are of physical interest ($\Delta F^{*}<10^{2}k_{B}T$) critical nuclei are nanoscopically small objects. Eqs.(1-6) describe the nuclei in terms of a simple variable $V_{n}$ and then $\Delta F(V_{n})$ has just a simple maximum at $V_{n}^{*}$, however, in reality nuclei have many degrees of freedom, and instead of a maximum of a simple function one deals with a saddle point of the coarse-grained free energy density in a high-dimensional configuration space 18 ; 20 ; 21 . Some of the resulting fluctuation corrections are commonly attributed to capillary wave excitations of the interface 22 . The latter occur only for rough fluid-crystal interfaces, of course, but have been observed in computer simulations of various models 23 ; 24 ; 25 ; 26 ; 27 and their analysis is useful for the estimation of the interfacial stiffness $\tilde{\gamma}$ of these interfaces. While such fluctuation corrections are typically of order $ln(V_{n})$, a curvature correction of order $V_{n}^{1/3}$ is also controversially discussed (mostly for the case of vapor-liquid interfaces, e.g. 28 ; 29 ; 30 ; 31 ; 32 ; 33 ), sign and magnitude of this correction are uncertain.

While for nucleation of liquid droplets from the vapor the two phases differ by the value of a scalar order parameter, the density, it is rather natural to use a description of the nucleus in terms of a simple variable as a starting point, namely the nucleus radius or volume. However, for the liquid to solid transition the distinction of the phases in terms of order parameters is very different: with crystallization both the translational and orientational symmetry of the liquid is broken 34 , and there is an accompanying density change (see e.g. 35 for a simulation study of synchronized ordering and densification). Clearly, such considerations must enter a first-principles discussion of crystal nucleation, but this is beyond the scope of the present work, which focusses on the effect that the nucleus shape is in general nonspherical.

As a consequence, quantitatively reliable predictions of nucleation barriers are difficult. For the vapor-liquid transition the nucleus on average has a spherical shape, and the interface tension $\gamma$ can be reliably measured (in the case of experiments) or predicted from models for the intermolecular forces (in the case of computer simulations), and thus the limiting case of the classical theory for the nucleation barrier (Eq. (6) with $A_{W}=A_{iso}$, cf. Eq. (2)) of large enough nuclei can be easily worked out. For the liquid-solid transition, however, $\gamma(\vec{n})$ and hence $\tilde{\gamma}$ cannot be easily measured in experiments, and finding it in simulations is a major effort and only possible approximately (e.g. 27 ; 36a ). Thus, working out $\Delta F^{*}\{$Eq. (6)$\}$ according to the classical nucleation theory is problematic, and hence tests of the theory by experiment can hardly yield conclusive results (see e.g. 36 ).

In this situation, a computer simulation approach to the problem appears as the most promising choice, but it also encounters conceptual problems: many techniques for obtaining $\Delta F^{*}$ of a nucleus require a distinction which particles are part of the crystal and which particles are part of the surrounding fluid 37 ; 38 . There is some ambiguity concerning the particles in the interfacial region at the surface of the nucleus, however. Moreover, a single critical nucleus surrounded by “supersaturated” (or supercooled, respectively) fluid at constant pressure in a macroscopically large volume is intrinsically unstable. As an object just at the top of the nucleation barrier it is bound to either grow or decay: thus performing a sampling of its properties one usually needs to apply a biasing potential 30 ; 37 to stabilize the nucleus at a pre-chosen size.

In the present paper, we present an alternative simulation approach, where the stable equilibrium between crystal nucleus and surrounding liquid in a finite volume is studied in the NVT ensemble. This stable equilibrium is a consequence of the lever rule 39 which describes two-phase coexistence in the canonical ensemble in the thermodynamic limit: when the density $\rho$ of the particles in the volume V exceeds the density $\rho_{f}$ describing the onset of freezing, but is still much smaller than the density $\rho_{m}$ describing the onset of melting of a crystal in equilibrium, we must have a crystal of volume $V_{n}$ such that ($N=\rho V$ while $N_{n}=\rho_{m}V_{n}$)

[TABLE]

Eq. (8) holds in the limit $V\rightarrow\infty$, $N\rightarrow\infty$ taken such that $\rho_{f}<\rho<\rho_{m}$, and both pressure $p$ and chemical potential $\mu$ are constants independent of $\rho$,

[TABLE]

[TABLE]

Here we have denoted the branches of the pressure $p(\rho,T)$ in the pure phases as $p_{l}(\rho,T)$ in the liquid and $p_{c}(\rho,T)$ in the crystal, and in the same way for $\mu$ ($\mu_{l}(\rho,T)$, $\mu_{c}(\rho,T)$). When one generalizes Eq. (8) to a finite volume, we need to specify the boundary conditions (we shall consider periodic boundary conditions (PBC) throughout) and the convention used for dividing up the particles between crystal nucleus and surrounding fluid: we use here the Gibbs convention of an equimolar dividing surface, so that $N=N_{f}+N_{n}$, $N_{f}$ being the particle number in the fluid and $N_{n}$ the particle number in the nucleus, no particle number excess associated with the interface being permitted. Then Eq. (8) still holds, but the applicable density range is restricted to $\rho_{1}<\rho<\rho_{2}$, where $\rho_{1}$ is the density where the “droplet evaporation/condensation transition” 40 ; 41 ; 42 ; 43 ; 44 ; 45 ; 46 occurs, while $\rho_{2}$ is the density where a transition from a compact droplet shape to a cylindrical droplet shape (stabilized by the PBC) occurs (Fig. 2).

Now due to nonnegligible interfacial contributions to the thermodynamic potential Eqs. (9),(10) no longer hold, both chemical potential and pressure are enhanced in comparison with their coexistence values,

[TABLE]

While in the presence of density inhomogenities in a system and corresponding interfacial contributions to its total free energy the chemical potential in equilibrium must still be spatially homogeneous, the pressure is not (Laplace pressure effect).

Due to the shift of equilibrium conditions from $\mu_{coex}$ to the chemical potential $\mu$ as specified in Eq. (11), the coexisting domains in the finite system no longer have the densities $\rho_{f}$ and $\rho_{m}$, but also have enhanced values $\rho_{l}$ and $\rho_{c}$, as anticipated in Eq. (11),(12), and hence Eq. (8) in the finite system gets replaced by

[TABLE]

The idea of the simulation approach described in this paper hence is a numerical analysis of the equilibrium described by Eqs. (11)-(13): choosing $\rho$ and $V$ suitably, a crystal with volume $V_{n}$ is equilibrated. The densitiy $\rho_{l}$ of the fluid in its surroundings as well as the associated pressure $p_{l}(\rho_{l},T)$ and chemical potential $\mu_{l}(\rho_{l},T)$ can then be calculated from the simulation. Since $\mu_{c}(\rho_{c},T)=\mu_{l}(\rho_{l},T)$, the density $\rho_{c}$ can also be extracted (from a suitable thermodynamic integration procedure). Then Eq. (13) yields the associated value of $V_{n}$ readily. If $V_{n}$ is large enough such that corrections to Eqs. (1),(2) are negligible, one can readily apply Eq. (6) to obtain the nucleation barrier from $\Delta F^{*}=\frac{1}{2}(p_{c}-p_{l})V_{n}^{*}$, and no consideration of $\gamma(\vec{n})$ and the resulting crystal shape is required. There will be two important consistency checks for the validity of this approach: (i) for large $V_{n}$ we must have $\Delta F^{*}\propto V_{n}^{2/3}$, and (ii) the results for $\Delta F^{*}$ must be independent from the total volume $V$. Both conditions are not expected to be true if $V$ is too small, when finite size effects related to fluctuations (which have nonzero correlation lengths $\xi_{l}$, $\xi_{c}$ in the pure fluid and crystal phases) come into play. Obviously this approach is not well suited, and not intended for, the study of ultrasmall nuclei containing a few particles only; we rather address the prediction that the classical theory of nucleation makes for crystal nucleation barriers. Although for very small nuclei corrections to this classical prediction are to be expected, it nevertheless is a very useful reference result, also for extensions to heterogeneous nucleation (e.g. 47 ; 48 ; 49 ; 50 ; 51 ; 52 ; 53 ; 54 ; 55 ; 56 ). The classical approach there considers sphere-cap shaped droplets on planar substrates under conditions of partial wetting 57a , and clearly for nucleation of crystals interfacial tension anisotropies are again expected to matter 55 .

In section II of the present paper we shall describe the model which we use to test the practical applicability of the approach outlined above, namely a version of the well-known Asakura-Oosawa model 57 ; 58 ; 59 ; 60 ; 61 of colloidal particles immersed in a fluid where added polymers provide a tunable effective attractive potential among the colloids 62 . Unlike the original model where the colloid-colloid interaction is a hard-core potential, we use a variant with a strongly repulsive continuous potential 63 , the so-called softEffAO model 63 ; 64 . The advantage of this model is that the pressure in the liquid phase, which is of crucial importance in our analysis, can be straightforwardly computed from the virial formula 65 , which cannot be applied for hardcore potentials. As far as a possible comparison of our results to real colloid-polymer mixtures [66] would be concerned, we note that real colloids anyway never interact strictly with hard-core potentials 67 ; 68 . Using Monte Carlo simulations 69 ; 70 ; 71 in the $NpT$ ensemble with PBC, the equation of state of homogeneous fluid and crystal phases is readily obtained very accurately, varying the temperature-like control parameter over a wide range.

Sec. III then discusses the precise estimation of coexistence conditions. We compare results from a method based on kinetics of a “slab” configuration (if $p\neq p_{coex}$ the slab grows or shrinks) 25 ; 26 to finite size extrapolations of the equilibrium described by Eqs. (11)-(13) to the thermodynamic limit, and find both approaches in very satisfactory agreement.

Sec. IV describes the preparation of the systems containing a nucleus surrounded by fluid in thermal equilibrium and discussses their analysis. It is shown that $\Delta F^{*}(V_{n}^{*})$ indeed is nicely consistent with Eq. (6), and does allow a rather accurate estimation of the effective surface tension $\gamma_{\text{eff}}=(A_{W}\bar{\gamma})/A_{iso}$ of the fluid-crystal interface. Evidence is presented that $\gamma_{\text{eff}}$ exceeds the surface tension $\gamma_{111}$ of the close-packed (111)-interface of the face-centered cubic crystal distinctly. Finally Sec. V contains some conclusions and gives an outlook on open problems.

II The soft effective Asakura-Oosawa (SoftEffAO) model and its bulk phase behavior

The Asakura-Oosawa (AO) model 57 ; 58 ; 59 ; 60 ; 61 describes colloids as hard spheres of diameter $\sigma_{c}$ and polymers as soft spheres (which may overlap each other with zero energy cost) of diameter $\sigma_{p}$, while overlap of colloids and polymers is strictly forbidden. Due to the ideal gas-like behavior of the polymers, one can integrate out their degrees of freedom exactly, thereby generating an effective colloid-colloid attraction (depletion potential). For a size ratio of $q=\sigma_{p}/\sigma_{c}<q^{*}=2/\sqrt{3}-1\cong 0.1547$ 62 this effective potential is a simple pairwise interaction,

[TABLE]

while $U_{\text{eff}}(r<\sigma_{c})=\infty$ and $U_{\text{eff}}(r\geq\sigma_{c}+\sigma_{p})=0$. Here $r$ is the distance between the centers of mass of the colloids, and $z_{p}=\exp(\mu_{p}/k_{B}T)$. To avoid the hard core singularity at $r=\sigma_{c}$, we rather add another repulsive potential to Eq. (14),

[TABLE]

which is a smooth function, diverging at $r=\epsilon\sigma_{c}$ and vanishing at $r=\sigma_{c}+q$. Useful choices for the constants $b$ and $\epsilon$ are $b=0.01$, $\epsilon=0.98857$ (of course, other choices of these constants or other choices of a smooth function replacing the hard core potential would equally be possible). Fig. 3 shows the resulting potential $U(r)=U_{\text{eff}}(r)+U_{\text{rep}}(r)$ of this SoftEffAO model for three choices of the parameter $\eta_{p}^{r}$, commonly denoted as “polymer reservoir packing fraction”,

[TABLE]

for the size ratio $q=0.15$.

Note that a solution of the polymers at concentration $\rho_{p}$ with no colloids present would have a polymer packing fraction of $\eta_{p}=(\pi\sigma_{p}^{3}/6)\rho_{p}=\eta_{p}^{r}$, of course; so variation of the polymer concentration in the colloidal suspension allows to control the strength of the colloid-colloid attraction. The fact that the range of this potenial, Eqs. (14), (15) is extremely short (and unlike the Lennard-Jones potential no artificial truncation is required) ensures a fast performance of the Monte Carlo simulation code.

To study the pure fluid and crystal phases of the system, it is most convenient to carry out standard $NpT$ Monte Carlo simulations, where apart from particle displacements in the cubic simulation box also volume changes are carried out 69 ; 70 ; 71 . It was found that using $N=4000$ finite size effects are completely negligible, and hence the observed volume $\left<V\right>$ at a given pressure $p$ and inverse temperature-like variable yields the same equation of state as for a calculation where the pressure $\left<p\right>$ is sampled in the canonical $NVT$ ensemble. It is common to convert the density $\rho=N/V=1/v$, where $v$ is the volume per particle, to the colloid packing fraction $\eta$,

[TABLE]

Fig. 4 shows the equation of state found by these $NpT$ simulations for four choices of $\eta_{p}^{r}$, in the $\eta-p$ plane.

Ideally the fluid branch $p_{l}(\eta)$ should end at $p_{coex}=p_{l}(\eta_{f})$ and exactly at the same pressure $p_{coex}=p_{c}(\eta_{m})$ the crystal branch should start, so both phases together should exist only for the single pressure $p_{coex}$. However, the data reveal considerable hysteresis: there is a rather extended pressure interval $\Delta p$ aroud $p_{coex}$, such that for pressures from this interval both phases are found. For a pressure $p\neq p_{coex}$ from this interval only one phase is stable, the other being metastable. In fact, the estimation of the equilibrium coexistence conditions ($p_{coex}$, $\eta_{f}$, $\eta_{m}$) is a delicate and nontrivial task, and will be discussed in the next section. When this task is solved, the phase diagram can be drawn as a diagram of variables $\eta_{p}^{r}$ and $\eta$ that clearly reveals the fluid + fcc crystal two-phase coexistence region (Fig. 5).

From Fig. 5 it is evident that the character of two-phase coexistence gradually changes from liquid/solid at small $\eta_{p}^{r}$, where $\eta_{f}$ is only a bit smaller than $\eta_{m}$, to gas/solid for large $\eta_{p}^{r}$, when $\eta_{f}\rightarrow 0$ while $\eta_{m}$ almost corresponds to close packing of hard spheres. The possibility that $\eta_{f}$ changes continuously from large values (typical for crystallization of melts) to very small values (typical for crystallization from the vapor phase) is a particular advantage of the present model which lacks a vapor to liquid transition and hence no triple point occurs; for the Lennard-Jones model, on the other hand, $\eta_{f}^{lc}$ for the liquid to crystal transition can decrease only a little with increasing temperature $T$, until the value $\eta_{f}^{lc}(T_{t})$ at the temperature of the triple point is reached, and then vapor to crystal nucleation starts at $\eta_{f}^{vc}(T_{t})\ll\eta_{f}^{lc}(T_{t})$; thus there is a wide regime of values for $\eta_{f}$ for which no crystallization/melting transition can be studied.

In Sec. I it was already pointed out that a particularly useful quantity to interpret simulations of phase coexistence is the chemical potential $\mu$. For the choice $\eta_{p}^{r}=0.28$ the packing fractions at the liquid branch are small and then the Widom particle insertion method 72 can be straightforwardly applied and yields accurate results (Fig. 6).

However, for $\eta_{p}^{r}=0.1$ and $0.2$ the Widom particle insertion method cannot be applied straightforwardly. For $\eta_{p}^{r}=0.1$ the packing fractions of interest are near $\eta_{f}=0.494$ and there the success rate of particle insertions simply is too small. Therefore Statt et al. 63 ; 64 ; 73 proposed to use a system in a $L\times L\times L_{z}$ geometry with $L_{z}\gg L$, where near $z=0$ a soft repulsive wall potential $V_{\text{rep}}(z)$ is applied, which has a range $\sigma_{w}$ of about two colloid diameters. For the opposite wall a short-range strongly repulsive wall potential is chosen, while in $x$ and $y$ directions, PBC are still applied. Due to these walls at $z=0$ and $z=L_{z}$, a nonuniform packing fraction profile $\eta(z)$ develops in the system, but for large enough $L_{z}$ (e.g. $L_{z}=30$) this profile develops a strictly constant part $\eta(z)=\eta_{bulk}$ in the center of the system, e.g. for $z_{min}\leq z\leq z_{max}$ we can compute $\eta_{bulk}$ from an average over the profile accurately, $\eta_{bulk}=(z_{max}-z_{min})^{-1}\int_{z_{min}}^{z_{max}}\eta(z)dz$. Note that $\eta_{bulk}$ differs from the average packing fraction $\eta$ in the system due to wall excess contributions, of course 63 ; 73 . On the other hand, general principles of statistical mechanics 39 imply that in equilibrium the chemical potential of such a system with inhomogeneous density is strictly constant throughout the system, $\mu(z)\equiv\mu$, even for $z$ close to the repulsive wall where $\eta(z)$ is much smaller than $\eta_{bulk}$. Measuring $\mu(z)$ by the particle insertion method for $z<\sigma_{w}$ this fact has been verified in the simulation, and although one uses only a small part $\Delta z\ll L_{z}$ of the total simulated system for the estimation of $\mu$, a very good accuracy can practically be established 63 ; 64 ; 73 . We also have verified that for $\eta_{p}^{r}=0.2$ in the range $p\leq 2.5$, which corresponds to $\eta_{c}\leq 0.39$, the agreement between both methods to estimate $\mu(p)$ is very good. For larger pressures, however, the Widom particle insertion method leads to a systematic overestimation of $\mu$; we attribute this problem to the fact that for large $\eta$ (and hence large $p$) the distribution of the chemical potential that is sampled by the insertion method gets very asymmetric around its maximum. Analogous data for $\mu$ vs $p$ in the case $\eta_{p}^{r}=0.1$ can be found in 63 ; 64 . For the estimation of the nucleation barrier $\Delta F^{*}=\frac{1}{2}(p_{c}-p_{l})V_{n}^{*}$ $\{$Eq. (6)$\}$ we shall need the pressure $p_{l}$ of the fluid surrounding the crystal nucleus (of volume $V_{n}^{*}$) as well as the associated crystal pressure $p_{c}$. However, this is not easily accessible directly, and hence a thermodynamic integration procedure has been applied. We use the relation for the chemical potential of the crystal

[TABLE]

making use of the fact that $\mu_{c}(p_{c})=\mu_{l}(p_{l})$, which follows from Eqs. (11), (12). For large enough nuclei it turns out that $\eta_{c}(p_{c})$ differs only very little from $\eta_{c}(p_{coex})=\eta_{m}$, and then Eq. (18) gets simplified as

[TABLE]

and hence the knowledge of $p_{coex}$, $\mu_{coex}$ together with $V_{n}^{*}$, $p_{l}$, $\mu_{l}(p_{l})$ will suffice to allow the computation of the desired free energy barrier $\Delta F^{*}$. Hence we proceed to discuss the estimation of phase coexistence conditions in the bulk.

III Characterizing bulk two-phase coexistence

In the discussion of Fig. 4 we have already pointed out that in the $NpT$ simulations of bulk systems we find an extended range of pressures near $p_{coex}$ where either the fluid phase or the crystal phase is metastable, and during affordable lengths of the simulation run no transition to the stable phase occurs. Thus these data cannot be used to estimate $p_{coex}$ and hence $\eta_{f}$ and $\eta_{m}$.

To overcome this problem, we apply the approach of Zykova-Timan et al. 24 ; 25 setting up a two-phase configuration of coexisting fcc crystal and fluid separated by two interfaces in an elongated box of linear dimensions $L\times L\times L_{z}$ with $L_{z}=5L$ with PBC throughout, simulated in the $Np_{z}T$ ensemble. The orientation of the crystal axes is chosen such that 2 interfaces of (100)-type occur, consistent with the geometry, and $L$ is chosen such that an integer number of lattice planes fits into the simulation box in $x$- and $y$- directions without misfit at the chosen pressure; the data for $\eta_{c}(p)$ at the crystal branches in Fig. 4 are used to choose the lattice parameter and hence $L$ (ranging from $n=4$ to $n=12$ lattice cells in $x$- and $y$- directions) appropriately. Fig. 7 shows an example for $n=6$ and $\eta_{p}^{r}=0.28$, choosing both phases in one half of the box volume. Note that in this case the fluid rather is a vapor-like phase, while for $\eta_{p}^{r}=0.1$ it rather is like a dense melt, and the quantity that is actually measured is the total box volume $V(t)$, or equivalently the linear dimension $L_{z}(t)$, as a function of simulation time t; note that after initialization the system is first equilibrated in the $NVT$ ensemble, choosing particle displacement Monte Carlo moves only, to allow that local equilibrium is established in the interfacial regions. So the origin $t=0$ of simulation time $t$ means the end of this $NVT$ equilibration period, when one switches to the $Np_{z}T$ simulation, now including volume changes by moves $L_{z}\rightarrow L_{z}\pm\Delta L_{z}$, as is standard 69 ; 70 ; 71 . If $p=p_{z}<p_{coex}$, the crystal should be unstable and hence particles “evaporate” from the interfacial region into the fluid, and thus the volume of the system must increase, to maintain the proper (low) density of the fluid at the chosen pressure. On the other hand, if $p=p_{z}>p_{coex}$, the fluid is unstable, particles “condense” in the interfacial region, and the total volume shrinks. The actual observations of $V(t)$ or $L_{z}(t)$ are consistent with these expectations; of course, there are huge random fluctuations from run to run, and thus 50 to 120 runs were averaged in order to obtain meaningful accuracy (Fig. 8). Of course, the inspection of the “raw data” for these time evolutions of crystal growth or shrinking only allows for a rather rough estimation of $p_{coex}$. Thus we tried to fit an interface velocity by fitting $V(t)$ or $L_{z}$ as a linear function over some time interval, but typically there is a systematic dependence of the resulting interface velocity on the size $\Delta t$ of this interval.

Fortunately, the resulting estimate of $p_{coex}$ as the pressure where this interface velocity changes sign shows almost no dependance on this interval. However, there seems to be a significant finite size effect: thus our final estimates were obtained by carrying out an extrapolation versus $1/n^{2}$ towards the thermodynamic limit (Fig. 9). For $\eta_{p}^{r}=0.1$ (data shown in 63 ; 64 and not repeated here) the result was $p_{coex}=8.45\pm 0.04$, excluding the smallest value of n from the fit. Clearly the fluctuations in the data of Fig. 9 (and of 63 ; 64 ) are large, and do not prove that the asymptotic regime of the extrapolation versus the inverse interfacial area really has been reached; thus a possible systematic error in these estimates for $p_{coex}$ cannot be excluded.

Clearly, it would be nice to have alternative methods to obtain $p_{coex}$. Deb et al. 55 suggested to use a $L\times L\times L_{z}$ geometry in the $NVT$ ensemble, with two planar repulsive walls at $z=0$ and $z=L_{z}$, and choose a packing fraction such that $\eta_{f}<\eta<\eta_{m}$. The system is initialized using an initial guess for $\eta_{m}$ (note that only when $p_{coex}$ has been determined are $\eta_{f}$ and $\eta_{m}$ accurately known!) to choose $L$ compatible with an integer number of lattice planes parallel to the $xy$-plane, and one chooses a few crystal layers parallel to the walls at $z=0$ and $z=L_{z}$, such that the remaining volume is compatible with the fluid phase at packing fraction $\eta_{f}$. The idea is that the system is self-regulating: if the packing fraction in the volume region available for the fluid is too high, particles will condense in the interfacial region, to make the crystalline layers a bit thicker, until $\eta_{f}$ in the fluid region would be established. However, if the initial guess for $\eta_{m}$ is not accurate enough, the crystal at $p_{coex}$ would be under strain, and hence this method would lead to a systematic error in the estimation of $p_{coex}$ and $\eta_{f}$ as well. Thus one needs to test a range of estimates for $\eta_{m}$ to find the choice where a $p_{coex}$ results for which the crystalline layers are undistorted. Of course, since not too few lattice planes parallel to the planar walls must be used at both walls, $L_{z}$ has to be large, and hence equilibration is rather slow. As a consequence, the accuracy reached by this method is not better than the method described in Figs. 7-9, and yields significantly lower values for $p_{coex}$ (see Fig. 10). Still another variant of phase coexistence simulation method due to Pedersen 74 studies the interface between crystal and fluid under the action of a force that prevents its motion if $p\neq p_{coex}$, and the magnitude of the force is varied until one finds the pressure for which the interface does not move for zero force. However, we have not yet tried out this approach for the present model, but rather we have used an approach based on the phase coexistence of compact droplets with surrounding fluid, as will be described in the next section. Combining Eqs. (5), (19) with an analogous expansion for $p_{l}$

[TABLE]

we find

[TABLE]

This relation suggests to plot the chemical potential $\mu_{l}$ of the fluid phase surrounding the crystal nucleus linearly vs. $(V_{n}^{*})^{-1/3}$: if the data fall on a straight line that is independent of the total volume $V$ in Eq. (13), we can estimate $\mu_{coex}$ as an intercept of the straight line for $(V_{n}^{*})^{-1/3}\rightarrow 0$. Fig. 11 shows that this recipe indeed works, and the accuracy seems to be competitive with the accuracy of the interface velocity method 25 ; 26 . While the two methods do not agree within the error bars, resulting from the respective fits, we stress that these errors do not include systematic errors, e.g. due to higher order terms in Eq. (21) or the problem that in Fig. 9 not all data points fall in the regime of $n^{-2}$ variation.

In principle, one can try a linear expansion similar to Eq. (21) also for the pressure $p_{l}$ to find $p_{coex}$, and the packing fraction $\eta_{l}(p_{l})$ to find $\eta_{f}$. When one does this, one finds for $\eta_{p}^{r}=0.1$ that $p_{coex}=8.47\pm 0.03$ and $\eta_{f}=0.498\pm 0.001$, close to the estimates extracted from the interface velocity method. For $\eta_{p}^{r}=0.2$, the corresponding results are $p_{coex}=1.59\pm 0.02$ and $\eta_{f}=0.376\pm 0.002$, while for $\eta_{p}^{r}=0.28$ the functions $p_{l}$ vs. $(V_{n}^{*})^{-1/3}$ and $\eta_{l}$ vs. $(V_{n}^{*})^{-1/3}$ exhibit pronounced curvature, and hence a linear extrapolation then does not work, while a nonlinear extrapolation certainly would be compatible with the results found from the equation of state when $\mu_{coex}$ is used.

As a result, we conclude that inaccurate knowledge of coexistence conditions is an important factor limiting the possibility to perform a high precision test of the predictions of homogeneous nucleation theory.

IV Analysis of the equilibrium between the crystalline nucleus and surrounding fluid

While in the case of vapor-liquid phase coexistence in a finite simulation box in the $NVT$ ensemble one can sample the free energy of the system (and the chemical potential as function of density, see Fig. 2, via a suitable numerical derivative) by umbrella sampling methods 43 ; 44 ; 45 , this is not yet feasible for fluid-crystal phase coexistence, at least not for the desired large simulation box volumes and particle numbers. Thus we have opted for an approach where a crystal nucleus of roughly the desired size is put into a volume that is suitably cut out from a configuration of a uniform fluid phase, choosing the density of that phase (as well as the packing fraction of the crystal) roughly in the region where we expect that the system will find its equilibrium.

Of course, even if the chosen packing fraction of the crystal (close to $\eta_{m}$) and of the surrounding fluid (somewhat larger than $\eta_{f}$) are in the right range, so that the system can settle down to an equilibrium of the type of Eq. (13), one must make sure that this equilibrium does not depend on irrelevant details of the inital state preparation. To test this issue, we have prepared different crystalline seeds for the same value of $\eta$ for a system of $N=10000$ particles (Fig. 12).

It is seen that the equilibrium that is established indeed is uniquely defined by the chosen parameters $N$, $V$, and $\eta_{p}^{r}$, and independent of the particular choice of initialization.

For the precise “measurement” of the density $\rho_{l}$ and the pressure $p_{l}$ of the fluid surrounding the nucleus it is important to exclude the volume region taken by the nucleus, as well as the adjacent interfacial region, from the sampling. For this purpose, we have to identify which particles belong to the nucleus, and which are definitely in the fluid region, no longer affected by the proximity of the nucleus. Following standard practice, we use the Steinhardt bond orientational parameters $q_{lm}(i)$, defined in terms of the spherical harmonics as 75 ; 76

[TABLE]

where $\vec{r}_{ij}$ is the distance between the considered particle $(i)$ to one $(j)$ of its $N(i)$ neighbors within a cutoff radius which defines what the “nearest neighbors” are. The polar angles $\theta(\vec{r}_{ij})$ and $\phi(\vec{r}_{ij})$ of this distance vector are chosen with respect to a fixed reference frame (which is arbitrary). The cutoff radius is determined from the first minimum of the radial distribution function $g(r)$ of the bulk fcc crystal. Following 76 $q_{lm}(i)$ is averaged over a particle and all its nearest neighbors,

[TABLE]

with $k_{i}=0$ meaning that one takes $q_{lm}(i)$ from Eq. (22). Averaging over $m$ we then obtain the root mean square order parameter $\bar{q_{l}}(i)$

[TABLE]

relevant choices being $l=4$ and $l=6$, allowing the distinction of fcc, hcp, bcc crystals from fluids (in our study, hcp and bcc order did not play any role). From sampling the histograms $P(\bar{q_{4}}$, $\bar{q_{6}})$ both in the pure fluid and crystal phases, as well as for systems where a crystal coexists with surrounding fluid, one can easily identify suitable criteria 64 ; 73 ; 77 to define regions in the $(\bar{q_{4}}$, $\bar{q_{6}})$-plane such that entries in these regions can be uniquely attributed to fluid-like or part of the fcc crystal, respectively.

Of course, one can also compute approximately the crystal volume from a radial density distribution of the nucleus around its center of mass 64 ; while $\bar{V_{n}}$ defined in this way is roughly compatible with the crystal volume $\bar{V_{n}}$ extracted from the lever rule, Eq. (13), we emphasize that the latter method does have the clear advantage that it does not assume anything about the shape of the crystal nucleus, and thus the results that will follow are based on Eq. (13) exclusively. The analysis in terms of Eqs. (24) is only used for identifying the region in which both $\eta_{l}$ and $p_{l}$ for the fluid coexisting with the nucleus are “measured”. The result of these “measurements” are now shown in Fig. 14.

One immediately recognizes that $p_{l}$ decreases with $\eta$, as expected from Fig. 2 where it was anticipated that in the region where actual coexistence of crystal nuclei with surrounding fluid occurs the chemical potential $\mu$ decreases with $\eta$, and $p_{l}$ is monotonously increasing with $\mu$ in equilibrium. It also is seen that the values $p_{l}(\eta)$ that we find are enhanced in comparison with $p_{coex}$ (Fig. 14 c); in cases a), b) $p_{coex}$ is even beyond the scale), but this enhancement systematically decreases with increasing $N$, and our analysis implies that these enhancements must converge to zero as $N\rightarrow\infty$, and then $p_{l}=p_{coex}$ irrespective of $\eta$.

It is very instructive to study the variation of the fluid pressure $p_{l}$ observed in our simulations with the density (or packing fraction, respectively) observed in this fluid. These do yield the left ascending curves in the plots, and it clearly is gratifying that these data no longer depend on the choice of $N$, and simply coincide with the equation of state that has been determined from $NpT$ simulations of homogeneous fluid phases, see Sec. II. This fact implies that our choices of $N$ are large enough that finite size effects that are related to fluctuations of the bulk phases are already negligibly small. Each choice of $N$ yields data points only for a limited range of $\eta$: if $\eta$ becomes too small, the data become unreliable due to the proximity of the droplet evaporation/condensation transition; if $\eta$ becomes too large, the data become unreliable due to the proximity of the transition where the droplet shape changes from compact to cylindrical. Of course, it is an implicit assumption in Eq. (13) that the effect of such fluctuations still is negligible: if e.g. the system would exist part of the observed simulation time in a state when the droplet either has evaporated or has made a transition to the cylindrical shape or even the slab configuration, the density and pressure of the fluid would significantly change in this part of the simulation time, leading to significant deviations in the curves on the left hand side of Fig. 14. In fact, setting up systems for each $N$ far outside the range shown in Fig. 14 one can verify that compact nuclei are unstable 77 , and either homogeneous fluid (at small $\eta$) and cylindrical nuclei or slab configurations (at large $\eta$) actually result 77 .

It is also evident from Fig. 14 that the absolute values of the slopes of the decreasing part clearly decrease with increasing $N$. This property, of course, must result as the thermodynamic limit is approached: for $N\rightarrow\infty$, $p_{l}(\eta)=p_{coex}$ irrespective of $\eta$, and thus the slope then is zero.

So the qualitative features of Fig. 14 are straightforward to understand from general principles. Using the fact that in Eq. (13) both $\rho$ and $V$ are fixed when $\eta$ and $N$ are chosen, $\rho_{l}=(6/\pi)\eta_{l}$ $\{$Eq. (17)$\}$ and the crystal density $\rho_{c}$ follows when we use the fact that the chemical potential $\mu_{c}(p_{c})$ of the crystal nucleus is equal to the chemical potential of the surrounding fluid (at density $\rho_{l}$, pressure $p_{l}$, which we have separately determined, see Sec. II above). For this purpose, in general one uses Eq. (18), once $p_{coex}$ and $\mu_{coex}$ have been determined, to compute $p_{c}$ and extract $\eta_{c}(p_{c})$ from the bulk equation of state (Sec. II). Note that in practice we have found that it suffices to linearize the equation of state around coexistence conditions (e.g. Eq. (19), see also Fig. 6 b)). But in the case of $\eta_{p}^{r}=0.28$ it is important to take into account the nonlinear variation of $\mu$ vs $p$; note also that in Fig. 14 c) the coexistence pressure $p_{coex}$ clearly falls in the region where the ideal gas law $p=(6/\pi)\eta$ holds, while most of the actual data for $p_{l}(\eta_{l})$ already deviate from it. So we conclude that Eq. (13) contains a single unknown parameter, namely the nucleus volume which hence is straightforwardly extracted. We have already anticipated that knowledge of $V_{n}^{*}$ provides an alternative method, Eq. (21), based on assuming the validity of Eq. (5), and indeed Fig. 11 has demonstrated the selfconsistency of this procedure, although slight systematic discrepancies between the different methods to precisely locate bulk coexistence conditions remain.

Since we have mentioned in Sec. I, that Eq. (2) is strictly correct only in the limit $V_{n}\rightarrow\infty$, and so one should discuss the effect of possible corrections that are expected to occur in Eq. (2) on our analysis. One correction that has been often identified for other systems (Ising model 15 ; 78 ; 79 , vapor-liquid transition 32 ; 78 etc.) is an additive constant as a correction term in Eq. (2); alternatively, this can be interpreted as a curvature correction to $\bar{\gamma}$ inversely proportional to the surface area. However, it trivially follows that such a constant correction to Eqs. (2, 3) does not affect the validity of Eq. (5).

This is no longer true for other corrections, such as fluctuation corrections leading to a $ln(V_{n})$ term in Eqs. (2, 3). One then expects that for large $V_{n}$ Eq. (3) is modified by a correction term of relative order $(V_{n}^{*})^{-2/3}$ and similarly in Eq. (21) a correction of order $(V_{n}^{*})^{-1}$ appears. Unfortunately, our data are not accurate enough to identify such a term. A Tolman-type 28 correction would be asymptotically larger, but since the prefactor of such corrections typically is very small 28 ; 29 ; 30 ; 31 ; 32 ; 33 , it cannot be identified either.

Disregarding all these corrections, we hence can use our results for $V_{n}$ and the associated liquid and crystal pressures $p_{l}$, $p_{c}$ to test Eq. (6). This is done in Fig. 15. It is seen that the data are well compatible with the relation $\Delta F^{*}=(4\pi/3)^{1/3}\gamma_{\text{eff}}V_{n}^{*^{2/3}}+const.$, as it should be for $V_{n}^{*}\rightarrow\infty$. We note that for $\eta_{p}^{r}=0.1$ the effective interfacial tension $\gamma_{\text{eff}}$ turns out to be only weakly enhanced ($\gamma_{\text{eff}}\approx 1.082$) in comparison with the actual interfacial tension $\gamma_{111}=1.013$ for the close packed (111) surface; also the surface tensions $\gamma_{110}=1.044$, $\gamma_{100}=1.039$ are not very different, and hence it is clear that the deviations from spherical nucleus shape can only be very minor, since $\gamma_{\text{eff}}=(A_{W}/A_{iso})\bar{\gamma}$ and $\bar{\gamma}\geq\gamma_{111}$, $A_{W}/A_{iso}>1$. However, for $\eta_{p}^{r}=0.2$ the anisotropy effects already are larger, $\gamma_{111}=2.078$, $\gamma_{110}=2.224$ and $\gamma_{100}=2.256$, while $\bar{\gamma}=2.406$. The surface tensions $\gamma_{111}$, $\gamma_{110}$, and $\gamma_{100}$ have been obtained elsewhere 73 ; 80 ; 81 with the “ensemble switch method”, a thermodynamic integration method constructing a path in phase space from two $L\times L\times L_{z}/2$ systems with PBC throughout without interfaces, one system being fluid and the other a bulk crystal, to a system of size $L\times L\times L_{z}$ where the two previously separate systems were joined together, such that two interfaces are formed. Clearly, we expect that when $\eta_{p}^{r}$ increases and the difference $\eta_{m}-\eta_{f}$ also increases strongly (Fig. 5) that also the differences between $\gamma_{111}$, $\gamma_{110}$, and $\gamma_{100}$ will increase. In the limit $\eta_{p}^{r}\rightarrow\infty$ (corresponding to the limit $T\rightarrow 0$ for a molecular system), the interface tension becomes equal to the excess of the potential energy due to the interface, since then entropic contribution to the interfacial tension becomes negligible. As discussed in Sec. II, the attractive potential is strictly short range, and in the fcc crystal it acts only between nearest neighbors. For a (111) plane of the fcc crystal structure, there are 6 nearest neighbors in the plane, and 3 nearest neighbors per particle in the adjacent planes. Hence for a (111) surface, 3 nearest neighbor bonds are missing. In the pairwise energy, each bond per particle enters with a factor $1/2$, and hence the surface excess energy per particle in the (111) surface energy is $-\frac{3}{2}U(r_{min})$, in the ground state, where $r_{min}$ is the distance where the potential $U(r)=U_{\text{eff}}(r)+U_{\text{rep}}(r)$ $\{$Eqs. (14), (15)$\}$ has its minimum. Converting from the surface excess energy per particle to the surface excess energy per unit area, one concludes

[TABLE]

for the potential energy part of the surface tension (which is an upper bound for the surface tension, since thermal fluctuations will reduce $U_{\text{exc}}$ in comparison with Eq. (25), and a further reduction occurs due to entropic contributions to $\gamma$).

But one sees that for $\eta_{p}^{r}=0.28$ the upper bound for $\gamma_{111}$ is already rather close to $\gamma_{\text{eff}}$, suggesting that for still larger $\eta_{p}^{r}$ this upper bound for $\gamma_{111}$ can already be used to obtain order of magnitude estimates for the actual surface tension. Of course, one must keep in mind that $\gamma_{\text{eff}}$ deviates from $U_{\text{exc}}$ due to two competing effects: $\gamma_{\text{eff}}$ is higher than $\gamma_{111}$ due to anisotropy effects; but there is a reduction of the surface tension (in comparison to $U_{\text{exc}}$) due to entropy 57a .

Finally, Fig. 17 shows a typical snapshot of a crystal nucleus for $\eta_{p}^{r}=0.28$.

Although it clearly does not a have a perfectly regular facetted shape, it is not a rough spherical shape either.

V Conclusions

In this paper a method has been described by which nucleation barriers for homogeneous nucleation of crystals from the fluid phase can be computed, emphasizing the limting case where the barrier is due to the competition between the difference in bulk thermodynamic potential between the fluid and solid phases with the surface free energy of a crystal nucleus which has the Wulff equilibrium shape when the crystal nucleus is large enough. It has been shown that the problem of finding this shape can be bypassed, analyzing the equilibrium properties (pressure, density, chemical potential) of the fluid coexisting with the crystal nucleus in a finite volume already contains all the information that is required to compute the nucleation barrier. This concept is exemplified for a simple model of colloidal particles with a short range effective attraction (caused physically by the depletion effect due to polymers, for instance). Choosing the size ratio of polymers and colloids $q=0.15$, the effective colloid-colloid attraction is strictly short-range and pairwise. In this model the fluid density where crystallization starts (or associate packing fraction $\eta_{f}$, respectively) is close to its counterpart in the crystal where melting starts ($\eta_{m}$) if the attraction is small, but the difference $\eta_{m}-\eta_{f}$ can be continously increased up to the limit where $\eta_{f}\rightarrow 0$ and $\eta_{m}$ approaches the packing fraction of close-packed spheres. The chosen softEffAO model can be studied by Monte Carlo simulation both in the $NVT$ and in the $NpT$ ensemble very efficiently; but the fluid-solid transition typically suffers from hysteresis, and very accurate location of coexistence conditions (pressure $p_{coex}$, chemical potential $\mu_{coex}$, etc.) is still one of the main factors limiting the accuracy of the present work. However, we feel that the situation is far more favorable than in experimental studies of colloidal systems, where neither the interaction potential nor the properties of the colloids themselves are known with sufficient accuracy to perform an experimental counterpart of our study (e.g., polydispersity of the colloids matters, ions in the solution adsorbing on the particle surface change the interactions slightly, etc.). Thus the accuracy with which experimental phase diagrams are measured does not match data such as shown in Figs. 4, 5; it is hoped that experimentalists will take up this challenge.

In our work, we have found that the conventional theory of homogeneous nucleation works nicely, the nucleation barriers $\Delta F^{*}(V_{n}^{*})$ scale with the nucleus volume $V_{n}^{*}$ as $\Delta F_{n}^{*}=(\frac{4\pi}{3})^{1/3}\gamma_{\text{eff}}V_{n}^{*^{2/3}}+const.$, $\gamma_{\text{eff}}$ being slightly but systematically enhanced in comparison with the surface tension $\gamma_{111}$ of close-packed (111) surfaces in the fcc crystal lattice. Experimental work often uses a similar expression, but $\gamma_{\text{eff}}$ is taken as a fit parameter that depends distinctly on the fluid density. We do not need to make such unjustified ad hoc assumptions here.

For the range of the free energy barriers for which we confirm the validity of the classical theory of homogeneous nucleation in general, namely $10^{2}k_{B}T\leq\Delta F^{*}\leq 10^{3}k_{B}T$, homogeneous nucleation is not directly observable, of course. However, we note that for many instances of heterogeneous nucleation the corresponding barriers are much smaller, and can only be predicted if the barriers for homogeneous nucleation are known. E.g., when one considers heterogeneous nucleation on planar substrates for conditions of partial wetting, the corresponding barriers are $\Delta F_{het}^{*}=\Delta F_{homo}^{*}~{}f_{VT}(\theta)$, where the reduction factor $f_{VT}(\theta)$ predicted by Vollmer and Turnbull is of order $f_{VT}\approx 10^{-1}$ if the contact angle is small. In any case, the knowledge of what the classical theory would predict also is a useful point of reference under conditions where the classical theory does not hold.

The present approach should be of general applicability, irrespective of the type of material that is studied. Thus, we hope to report on further applications, such as nucleation of ice crystals from supercooled water, in the future.

Acknowledgements.

Partial support from the Deutsche Forschungsgemeinschaft (grant VI 237/4-3) and the Graduate School of Excellence Materials Science in Mainz is acknowledged. P. Koß is a recipient of a DFG-fellowship/DFG-funded position through the Excellence Initiative by the Graduate School Materials Science in Mainz (GSC 266). The authors gratefully acknowledge the computing time granted on the supercomputer Mogon at Johannes Gutenberg University Mainz (hpc.uni-mainz.de) and Hornet at the Höchstleistungsrechenzentrum Stuttgart (HLRS).