Phase-field modeling of crystal nucleation in undercooled liquids -- A review

TL;DR

This review discusses how phase-field models, including classical and crystal-specific variants, have advanced understanding of crystal nucleation and microstructure evolution in undercooled liquids, covering various nucleation phenomena and morphologies.

Contribution

It provides a comprehensive overview of phase-field modeling approaches and their applications to diverse crystal nucleation processes in undercooled liquids.

Findings

Phase-field models capture homogeneous and heterogeneous nucleation.

Crystal growth morphologies depend on nucleation and growth dynamics.

Recent models operate on molecular scales using density functional theories.

Abstract

We review how phase-field models contributed to the understanding of various aspects of crystal nucleation including homogeneous and heterogeneous processes, and their role in microstructure evolution. We recall results obtained both by the conventional phase-field approaches that rely on spatially averaged (coarse grained) order parameters in capturing freezing, and by the recently developed phase-field crystal models that work on the molecular scale, while employing time averaged particle densities, and are regarded as simple dynamical density functional theories of classical particles. Besides simpler cases of homogeneous and heterogeneous nucleation, phenomena addressed by these techniques include precursor assisted nucleation, nucleation in eutectic and phase separating systems, phase selection via competing nucleation processes, growth front nucleation (a process, in which grains of new orientations form at the solidification front) yielding crystal sheaves and spherulites, and transition between the growth controlled cellular and the nucleation dominated equiaxial solidification morphologies.

Figures

Tables

Table 1. Table 1: A few g ​ ( ϕ ) 𝑔 italic-ϕ g(\phi) double-well and p ​ ( ϕ ) 𝑝 italic-ϕ p(\phi) interpolation functions employed in the literature.

Model	$g\left(\varphi \right)$	$p\left(\varphi \right)$	Ref.
standard PF	$\frac{1}{4}{\varphi }^2{\left(1-\varphi \right)}^2$	${\varphi }^3\left(10-15\varphi +6{\varphi }^2\right)$	ref85
GL bcc	$\frac{1}{4}{\varphi }^2{\left(1-\varphi \right)}^2$	${\varphi }^3\left(4-3\varphi \right)$	ref84
GL fcc	$\frac{1}{6}{\varphi }^2{\left(1-{\varphi }^2\right)}^2$	${\varphi }^4\left(3-2{\varphi }^2\right)$	ref84
GL bcc/diamond	$\frac{1}{4}{\varphi }^2{\left(1-\varphi \right)}^2$	${\varphi }^2$	ref86
dynamic $\tanh$	$\frac{1}{4}{\varphi }^2{\left(1-\varphi \right)}^2$	${\varphi }^2\left(3-2\varphi \right)$	ref85 ; ref87
double parabolic	$\min \{g_l,g_s\}$	$0$	ref88 ; ref88
ice	$\frac{1}{4}{\varphi }^2{\left(1-\varphi \right)}^2$	$\varphi$	ref90

Table 2. Table 2: Properties of Ni used in computing the nucleation barrier

Property	value	unit	Ref.
$m$	$58.6934\times {10}^{-3}$	kg mol-1	ref104
$\Delta H_f$	17.29	kJ mol-1	ref49 (a)
$T_f$	1728.3	K	ref105
$\Delta f$	$H_f\left(T-T_f\right)T_f^{-1}{v}_m^{-1}$	J m-3	ref106
$\rho$	8357	kg m-3	ref49 (a)
${\gamma }_{SL,eq}$	0.302	J m-2	ref107
$d_{10\%-90\%}$	1.0	nm	deduced from Ref. ref107

Equations

$$\dot{N_1}=a_2^{-}{N}_2-a_1^{+}{N}_1+\sum _{n>1}(a_n^{-}{N}_n-a_{n-1}^{+}{N}_{n-1}),$$

$$\dot{N_n}=a_{n-1}^{+}{N}_{n-1}+a_{n+1}^{-}{N}_{n+1}-(a_n^{+}+a_n^{-})N_n.$$

$$J_{SS,hom}=$$

$$J_{SS,het}=x_{a}q(\Psi )J_0\exp \{-W_{hom}^{\ast }f(\Psi )/kT\},$$

$$\displaystyle Y(t)=1-\exp\big{\{}-\tilde{Y}\big{\}},$$

$$\displaystyle Y(t)=1-\exp\big{\{}-[(t-t_{0})/\tau]^{p}\big{\}}.$$

$$\displaystyle F=\int d\mathbf{r}\bigg{\{}\frac{\epsilon_{\phi}^{2}}{2}\sum_{i,j=1}^{N}A_{ij}(\nabla\phi_{i}\cdot\nabla\phi_{j})+\frac{\epsilon_{c}^{2}}{2}\sum_{i,j=1}^{M}B_{ij}(\nabla c_{i}\cdot\nabla c_{j})+\cdots+f(\mathbf{\Phi},\mathbf{c},T,\cdots)\bigg{\}},$$

$$-\frac{\partial {\varphi }_i}{\partial t}=\sum _{j=1}^N{\kappa }_{ij}\left(\frac{\delta F}{\delta {\varphi }_i}-\frac{\delta F}{\delta {\varphi }_j}\right).$$

$$\frac{\partial c_i}{\partial t}=\nabla \cdot \left[\sum _{j\ne i}{\kappa }_{ij}\nabla \left(\frac{\delta F}{\delta c_i}-\frac{\delta F}{\delta c_j}\right)\right].$$

$$\frac{\delta F}{\delta {\varphi }_i}=\lambda (\mathbf{r}),i=1\dots N.$$

$$\frac{\delta F}{\delta {\varphi }_i}=\frac{\delta F}{\delta {\varphi }_j}\text{for}i,j=1\dots N.$$

$$\displaystyle F=\int d\mathbf{r}\bigg{\{}\frac{\epsilon^{2}s^{2}({\bf n})}{2}|\nabla\phi|^{2}+wg(\phi)+p(\phi)\Delta f\bigg{\}},$$

$$0=I-{\varphi }_x\frac{\partial I}{\partial {\varphi }_x},$$

$$\displaystyle\frac{\delta F}{\delta\phi}=\bigg{\{}\frac{\partial I}{\partial\phi}-\nabla\cdot\frac{\partial I}{\partial\nabla\phi}\bigg{\}}=0,$$

$$\frac{{\partial }^2\varphi }{\partial r^2}+\frac{2}{r}\frac{\partial \varphi }{\partial r}=\frac{1}{{ϵ}^2}\frac{\partial f(\varphi )}{\partial \varphi },$$

$$\displaystyle\dot{\phi}=-M_{\phi}\frac{\delta F}{\delta\phi}+\zeta=M_{\phi}\big{\{}\epsilon^{2}\nabla^{2}\phi-wg^{\prime}(\phi)-p^{\prime}(\phi)\Delta f\big{\}}+\zeta,$$

$$F=\int_{S}Z(\phi)dS+\int dV\bigg{\{}f(\phi)+\frac{\epsilon^{2}}{2}(\nabla\phi)^{2}\bigg{\}},$$

$$\delta\phi\big{\{}z(\phi)+\epsilon^{2}\nabla\phi\cdot\mathbf{n}\big{\}}=0\quad\mathrm{on}\ S.$$

$$\mathbf{n}\cdot \nabla \varphi =(2\Delta f/{ϵ}^2{)}^{1/2}\cos (\Psi )\propto \varphi (1-\varphi )\mathrm{on}S.$$

$$\displaystyle F=\int d\mathbf{r}\bigg{\{}\frac{\epsilon_{\phi}^{2}s^{2}(\theta)T}{2}|\nabla\phi|^{2}+w(c)Tg(\phi)+p(\phi)f_{s}(c,T)+[1-p(\phi)]f_{l}(c,T)\bigg{\}},$$

$${ϵ}^2{\nabla }^2\varphi$$

$${ϵ}_c^2{\nabla }^{2}c$$

$$(\mathbf{n}\cdot \nabla \varphi )=\sqrt{\frac{2\Delta f[\varphi ,c(\varphi )]}{{ϵ}^2}}\cos (\Psi ),$$

$$\displaystyle F=\int dV\bigg{\{}Z(\phi)|\nabla\phi_{W}|+\left[f(\phi,c)+\frac{\epsilon^{2}}{2}|\nabla\phi|^{2}\right](1-\phi_{W})\bigg{\}},$$

$$\displaystyle-\frac{1}{M_{\phi}}\frac{\partial\phi}{\partial t}=\frac{\delta F}{\delta\phi}=\bigg{\{}\frac{\partial f}{\partial\phi}-\epsilon^{2}\nabla^{2}\phi\bigg{\}}(1-\phi_{W})+\left[z(\phi)+\epsilon^{2}\nabla\phi\cdot\mathbf{n}\right]|\nabla\phi_{W}|.$$

$$F$$

$$\displaystyle\frac{\partial\theta}{\partial t}=-M_{\theta}\frac{\delta F}{\delta\theta}=M_{\theta}\nabla\bigg{\{}p(\phi)h\frac{\nabla\theta}{|\nabla\theta|}-\epsilon^{2}s\frac{\partial s}{\partial\theta}(\nabla\phi)^{2}\bigg{\}},$$

$$\displaystyle f_{ori}=2hp(\phi)\bigg{\{}\sum_{i=1}^{4}(\nabla q_{i})^{2}\bigg{\}}^{1/2},$$

$$\displaystyle\frac{\partial q_{i}}{\partial t}=-M_{q}\frac{\delta F}{\delta q_{i}}+\zeta_{i}=M_{q}\nabla\bigg{\{}\left[hp(\phi)\frac{\nabla q_{i}}{(\sum_{l}|\nabla q_{l}|^{2})^{1/2}}\right]-q_{i}\sum_{k}q_{k}\nabla\left[hp(\phi)\frac{\nabla q_{k}}{(\sum_{l}|\nabla q_{l}|^{2})^{1/2}}\right]\bigg{\}}+\zeta_{i},$$

$$\displaystyle\mathcal{F}=\int d\mathbf{r}\bigg{\{}\frac{\psi}{2}[-\epsilon+(1+\nabla^{2})^{2}]\psi+\frac{\psi^{4}}{4}\bigg{\}},$$

Full text

1. I Introduction

2. II Definitions, notions, classical theory

3. II.1 Homogeneous nucleation

4. II.2 Heterogeneous nucleation

5. II.2.1 Particle induced freezing: Free growth limited mechanism

6. II.3 Comparison to experiments and molecular simulations

7. II.3.1 Lennard-Jones (LJ) system

8. II.3.2 Hard-sphere (HS) system

9. II.3.3 Metals

10. II.3.4 Water-ice

11. II.3.5 Structural aspects

12. II.4 Formal theory of crystallization

13. III Coarse grained phase-field models

14. III.1 Application of the phase-field method to nucleation

15. III.1.1 Single-field models

16. III.1.2 Two-field models

17. III.1.3 Nucleation in PF models with three or more fields

18. III.2 Numerical methods

19. III.3 Nucleation vs. microstructure evolution

20. III.3.1 Transformation kinetics in alloys

21. III.3.2 Large-scale modeling of polycrystalline solidification morphology

22. III.3.3 Columnar to equiaxed transition

23. III.3.4 Polycrystalline growth

24. III.4 Discussion: Nucleation in coarse grained PF models

25. IV Molecular scale phase-field models of crystallization

26. IV.1 The Phase-Field Crystal approach

27. IV.1.1 The single-mode PFC model

28. IV.1.2 The two-mode PFC model

29. IV.1.3 The multi-mode PFC model

30. IV.1.4 Other advanced PFC models

31. IV.2 Application of the PFC models to crystal nucleation

32. IV.2.1 Euler-Lagrange equation and other methods to find free energy extrema

33. IV.2.2 PFC with diffusive dynamics (DPFC)

34. IV.2.3 Hydrodynamic PFC model of freezing (HPFC)

35. IV.2.4 Numerical methods

36. IV.3 Discussion

37. V Nucleation prefactor from phase-field/density functional theories

38. VI Summary and future directions

39. VII Disclaimer

I Introduction

The crystallization of ideally pure liquids cooled below their melting point starts with homogeneous nucleation, a process in which the internal fluctuations of the undercooled liquid lead to the formation of crystal-like seeds able to grow to macroscopic sizes. This process is normally assisted by the presence of surfaces (container walls, foreign particles, etc.) termed heterogeneous nucleation. These phenomena are of interest for various branches of science including physical chemistry, materials science, biophysics, geophysics, cryobiology, etc. and play important roles in a range of technologies.

Modeling of crystal nucleation has a long history covered in a number of reviews ref1 ; ref2 . The applied approaches range from discrete atomistic simulations relying on molecular dynamics (MD) ref3 ; ref4 ; ref5 ; ref6 , Monte Carlo (MC) ref7 ; ref8 ; ref9 , Brownian dynamics (BD) ref10 ; ref11 , and cluster dynamics techniques ref12 ; ref13 ; ref14 , to continuum models including the van der Waals/Cahn-Hilliard/Ginzburg-Landau/$\phi^{4}$ type models ref15 ; ref16 ; ref17 ; ref18 , based on the square gradient (SG) approximation, and the more complex phase-field ref19 ; ref20 and classical density functional methods ref21 ; ref22 . Among these, the phase-field (PF) approaches became the method of choice when the description of complex solidification structures (such as dendrites, eutectic patterns, spherulites, fractal-like structures, etc.) are required. The terminology changed with time: the SG approximation-based approaches were often considered as the simplest form of density functional theories (see e.g., Ref. ref21 ), whereas the phase-field theory originally meant an SG approach, in which a structural order parameter (the phase field) monitors the crystal-liquid transition ref23 . Recently, however, phase-field methods working on the molecular scale, termed phase-field crystal (PFC) models were introduced ref24 ; ref25 , which can be classified as simple dynamical density functional approaches.

Herein, we review the application of phase-field and phase-field crystal methods to nucleation problems. The models will be presented in a historical manner, illuminating their increasing predicitive power as the research progressed. The paper is structured as follows. In Section II, we briefly recall a few general ideas and notions that can be best introduced using the classical nucleation theory (CNT), then in Section III we review the work done using conventional PF models that rely on a coarse grained structural order parameter(s) and coupled fields in describing nucleation. The areas covered include homogeneous nucleation, phase selection via competing nucleation processes, heterogeneous nucleation, free growth limited particle induced freezing, growth front nucleation, techniques to implement nucleation into PF simulations, and microstructure evolution in the presence of nucleation. Section IV is devoted to the molecular scale phase-field studies: PFC results for homogeneous nucleation, amorphous precursor mediated crystal nucleation, heterogeneous nucleation on flat and modulated surfaces, particle induced crystallization, and growth front nucleation will be reviewed. In Section V, we cover developments concerning the nucleation prefactor. Finally, in Section VI a brief summary is given, and we outline directions in which further research appears promising.

II Definitions, notions, classical theory

II.1 Homogeneous nucleation

The classical approach views the crystal-like fluctuations appearing in the undercooled liquid as small spherical domains of the bulk crystalline phase bound by a mathematically sharp solid-liquid interface ref26 (known as the droplet model or capillarity approximation), while the formation, growth, and decay of these fluctuations is assumed to happen via a series of single-molecule attachment and detachment events. The work of formation of such crystallites is expressed as a sum of a volumetric and an interfacial term: $W_{hom}=(4\pi/3)R^{3}\Delta\omega+4\pi R^{2}\gamma_{SL}$, where $R$ is the radius of the surface on which the surface tension acts (the surface of tension), $\Delta\omega$ is the thermodynamic driving force of solidification (the volumetric grand potential difference between the solid and the liquid; $\Delta\omega<0$ for the undercooled liquid), and $\gamma_{SL}$ is the free energy of the solid-liquid interface. Since for small sizes the (positive) surface term dominates, while in the case of large sizes the (negative) volumetric term is the leading one, the work of formation shows a maximum of height $W_{hom}^{*}=(16\pi/3)(\gamma_{SL}^{3}/\Delta\omega^{2})$ as a function of the size at $R^{*}=-2\gamma_{SL}/\Delta\omega$. The critical size $R^{*}$ and the nucleation barrier $W^{*}$ are essential features of the critical fluctuation or nucleus, also called as the critical cluster.

The kinetic part of the classical approach relates the nucleation rate, i.e., the net formation rate of critical fluctuations, to the attachment/detachment rates of the molecules to/from the crystalline clusters. The master equations governing the time evolution of the cluster population can be formulated as follows:

[TABLE]

and for $n>1$

[TABLE]

In Eqs. (1) and (2), $N_{n}$ is the number $n$-molecule clusters, while $a_{n}^{+}=O_{n}\Gamma\exp\{-(W_{n+1}-W_{n})/2kT\}$ and $a_{n}^{-}=O_{n-1}\Gamma\exp\{-(W_{n-1}-W_{n})/2kT\}$ denote the frequencies for molecule attachment and detachment. $O_{n}=4n^{2/3}$ is the number of surface sites on an $n$-molecule cluster to which liquid molecules can be attached, whereas $\Gamma=6D/\lambda^{2}$ stands for the time scale of molecule attachment/detachment. $D$ is the self-diffusion coefficient (often related to the viscosity via the Stokes-Einstein relationship), and $\lambda$ is the molecular jump distance. $W_{n}$ stands for the free energy of formation of an $n$-molecule cluster, $k$ is Boltzmann’s constant, and $T$ the temperature. While Eq. (2) has been used in many works ref1 (a), ref2 (a), ref13 ; ref14 , Eq. (1) describes the time evolution of the number of monomers: the depletion of monomolecular clusters due to dimer formation and attachment to larger clusters, and its increase via dimer dissolution and detachment from larger clusters ref14 (b), ref27 .

Solving Eqs. (1) and (2) numerically, steady state nucleation occurs after a transient period of length $\tau\approx K\lambda^{2}kT(n^{*})^{2/3}/(Dv_{m}\Delta\omega)$, where $K$ is a geometrical factor, $n^{*}$ is the number of molecules in the critical cluster, and $v_{m}$ the molecular volume. The steady state nucleation rate can be expressed as follows [see e.g., Ref. ref13 (a)]:

[TABLE]

Here $N_{eq,n}=\rho_{0}\exp\{-W_{n}/kT\}$ is the equilibrium population of the $n$-molecule clusters, and $J_{0}=\rho_{0}O_{n^{*}}\Gamma Z$ is the pre-exponential factor of the nucleation rate, while $\rho_{0}$ denotes the number density of the single molecule clusters (the molecules of the liquid), whereas the Zeldovich factor $Z=\{|d^{2}W_{hom}/dn^{2}|_{n^{*}}/(2\pi kT)\}^{1/2}$ accounts for the decay of the critical size clusters.

Fitting Eq. (3) with unknown $D$ and $\gamma_{SL}$ to nucleation rate experiments on oxide glasses imply that the order of magnitude of the classical estimate for $J_{0}$ is reasonable ref1 (a), ref2 (a), ref14 (a). Molecular dynamics investigations suggest that the classical prefactor for crystal nucleation might be 2 orders of magnitude too low ref4 . For vapor condensation, molecular dynamics simulations found $J_{0}$ that is $\sim 1$ order of magnitude larger than the CNT prediction ref28 , whereas a dynamic extension of the classical density functional theory yielded a reasonable agreement with the CNT result for $J_{0}$ ref29 . In the case of Monte Carlo simulations for the 2D Ising model, a good agreement was observed between the CNT and simulations, if $W^{*}$ from the droplet model was replaced by the proper cluster free energies in the classical kinetic framework ref30 , implying that the kinetic part of the CNT is reasonably accurate. A field theoretical expression, J_{SS}=J_{0}^{\prime}\exp\big{\{}-W^{*}/kT\big{\}}, similar to Eq. (3) has been derived by Langer ref31 , which lead to comparable results to CNT in the few cases in which comparison was made ref32 ; ref33 .

Apparently, the large (several orders of magnitude) deviation between experimental and theoretical (CNT) nucleation rates reported for oxide glasses and other substances ref1 (a), ref2 (a), ref14 (a), ref34 is attributable primarily to the failure of the droplet model. This view is supported by a direct evaluation of the nucleation barrier via umbrella sampling (a biased Monte Carlo technique) that shows that the droplet model relying on a constant $\gamma_{SL}$ fails when predicting the nucleation barrier ref5 (a). The failure of the droplet model for small clusters was also borne out by Monte Carlo simulations for the Ising model ref35 . In analogy to the case of liquid droplets, for small clusters, corrections may be introduced for the surface tension (interfacial free energy) ref36 . A fairly frequently used correction is Tolman’s that introduces a size dependent interfacial free energy: $\gamma(R_{p})\approx\gamma_{\infty}/[1+2(\delta_{T}/R_{p})]$, where $R_{p}=[3W^{*}/(2\pi\Delta f)]^{(1/3)}$ is the radius of the surface of tension on which the surface tension acts, $\delta_{T}=\lim_{R_{e},R_{p}\rightarrow\infty}(R_{e}-R_{p})$ is the Tolman length, and $R_{e}$ is the equimolar surface (the position of the step function that has the same amplitude and radial integral as the density profile). For details see Ref. ref37 .

II.2 Heterogeneous nucleation

In real liquids, crystal nucleation takes place normally in a heterogeneous manner: formation of crystal-like fluctuations is assisted by heterogeneities, such as container walls, floating particles, molecular impurities, etc. In such heterogeneous processes, the nucleation barrier may be reduced significantly $(W_{het}^{*}<W_{hom}^{*})$, leading to higher nucleation rates. In the CNT, the spherical cap model is used to quantify the catalytic effect of foreign particles floating in the undercooled liquid. It is assumed that the foreign particles are distributed homogeneously in liquid, they are considerably larger than the nuclei, and are bound by flat walls. In equilibrium, the relevant interface free energies are related to each other by the Young-Laplace equation ref38 : $\gamma_{WL}=\gamma_{WS}+\gamma_{SL}\mathrm{cos}(\Psi),$ where $\gamma_{WL}$ and $\gamma_{WS}$ stand for the wall-liquid and wall-solid interfacial free energies, whereas $\Psi$ is the contact angle. Under such conditions the critical fluctuation is a spherical cap (a fraction of the homogeneous nucleus that realizes the contact angle; Fig. 1). The respective work of formation can be expressed as $W_{het}^{*}=W_{hom}^{*}f(\Psi)$, where $f(\Psi)=\frac{1}{4}[2+\mathrm{cos}(\Psi)]\{1-\mathrm{cos}(\Psi)\}^{2}$ is the catalytic potency factor. For small contact angles, $f(\Psi)$ can be small, reducing the nucleation barrier significantly. The number of sites on the spherical cap to which liquid molecules can be attached is $O_{n}=2\{1-\mathrm{cos}(\Psi)\}n^{2/3}$. Since only those molecules may participate in heterogeneous nucleation, which are effectively adsorbed on the surface of heterogeneities, the steady state nucleation rate is expressed as ref39

[TABLE]

where $x_{a}\ll 1$ denotes that fraction of all the molecules that is adsorbed on the surface of heterogeneities, and $q(\Psi)=\frac{1}{2}\{1-\mathrm{cos}(\Psi)\}/\sqrt{f(\Psi)}\in[1/\sqrt{3},1]$. The classical approach has been adapted to various geometries of the wall, including spherical particles and cavities, depressions, and rough surfaces ref40 .

A generalized form of Turnbull’s experimental test ref41 for the classical nucleation theory can be devised on the basis of Eq. (4) provided that there is a theoretical estimate for the temperature dependence of the interfacial free energy (incorporating the effect of surface curvature) $\chi(T)=\gamma_{SL}(T)/\gamma_{SL,eq}$, and $D(T)$, $\lambda$, and $\Delta\omega(T)$ are known: plotting the logarithm of the normalized experimental (steady state) nucleation rate, $\mathrm{log}(J_{SS,exp}/J_{0})$, vs. the temperature dependent part of the argument of the exponential function, $\chi(T)^{3}\Delta\omega(T)^{-2}T^{-1}$, one should obtain a straight line that intersects the ordinate at $\log(x_{a}q(\Psi))\approx\log(x_{a})\leq 0$ with a slope of $-16\pi\gamma_{SL,eq}^{3}f(\Psi)/3k$ ref34 ${(c),(d)}$. For example, the assumption of a curvature- and temperature-independent interfacial free energy ($\chi=1$) yields intersections many orders of magnitude too high ref1 (a), ref2 (a), ref14 (a), ref34 , indicating a failure of the original droplet model that relies on constant $\gamma_{SL}$.

II.2.1 Particle induced freezing: Free growth limited mechanism

A particularly interesting case, in which volumetric heterogeneities play a decisive role, is the free growth limited mechanism of particle induced freezing proposed by Greer and coworkers ref42 . In this approach, the foreign particles floating in the melt are viewed as cylinders of radius $R_{p}$ that have ideally wetting circular faces ($\Psi=0$), and non-wetting sides ($\Psi=\pi$). These idealized particles remain dormant at and below a critical undercooling, $\Delta T_{c}=2\gamma_{SL}/(\Delta s_{f}R_{p})$, at which the radius of the particles is equal to that of the critical radius for the homogeneous nucleus, $R_{p}=R^{*}$. (Here $\Delta s_{f}$ is the volumetric entropy of fusion.) For undercoolings $\Delta T>\Delta T_{c}$, free growth takes place. The model has been adapted for various geometries of the foreign particles, including triangular and hexagonal prisms, cubes, and regular octahedra ref43 . The free growth limited mechanism of particle induced freezing proved highly successful under a broad range of conditions in materials science, cryobiology, and other branches of science, where foreign particles initiate solidification ref1 (a).

II.3 Comparison to experiments and molecular simulations

The central problem preventing conclusive experimental tests of nucleation theory is the lack of means (other than nucleation theory) to determine the solid-liquid interfacial free energy with sufficient accuracy. Although there are experimental methods for evaluating $\gamma_{SL}$ in the vicinity of the melting point ref44 , the associated error is usually far too high (5$\%$ to 10 $\%$) for a conclusive test. Furthermore, the interfacial free energy may depend on temperature ref45 and curvature ref34 (c),(d), ref36 , both of which could influence the nucleation rate considerably. In addition, despite a wealth of highly accurate nucleation rate data available for oxide glasses ref34 ; ref47 , it is also rather difficult to determine whether homogeneous or heterogeneous nucleation took place in the individual experiments. These together with the uncertainties (or complete absence) of the experimental data for $\gamma_{SL}$ make a rigorous experimental test of nucleation theory practically impossible.

In the past decade or so, detailed information became available from molecular simulations (MD, MC, and BD) for model systems ref48 and pure metals ref49 , and experiments for crystallizing colloidal suspensions: ref50 ; ref51 some of these approaches ref48 ; ref49 ; ref51 deliver the trajectory of the particles of the crystallizing liquid, while the interfacial and thermodynamic properties might also be known. One may summarize the emerging results as follows:

II.3.1 Lennard-Jones (LJ) system

Direct investigations of the nucleation barrier in the Lennard-Jones (LJ) system at temperature $T=0.6$ and pressure $P=0.67$ performed using umbrella sampling indicate a reasonable agreement between the simulations and the droplet model ($W^{*}_{MC}/kT=19.4$ vs. $W^{*}_{CNT}/kT=17.4$) ref4 . In computing the latter, the orientation average of $\gamma_{SL,eq}=0.35$ evaluated at the triple point was used from a MD study performed with a modified LJ potential ref52 . It appears, however, that (much like in the hard sphere system) $\gamma_{SL,eq}\propto T$ in the LJ system ref45 , which means that the previous value of $W^{*}_{CNT}/kT$ needs to be multiplied by $(T/T_{3})^{3}$, yielding $W^{*}_{CNT}/kT=12.0$ that is considerably less than the simulation result. To recover $W^{*}_{MC}/kT=19.4$ at $T=0.6$, one needs to assume a positive curvature correction from $\gamma_{SL,eq}^{T=0.6}=0.309$ to $\gamma_{SL}^{nucl}=0.363=1.17\gamma_{SL,eq}^{T=0.6}$ of the nucleus.

II.3.2 Hard-sphere (HS) system

A significant difference was reported ref5 (a) for $W^{*}$ from MC simulations and the droplet model that relies on the equilibrium value $\gamma_{SL,eq}/(kT/\sigma^{2})=0.617$ ref53 ; ref54 from a numerical estimate. (The values from the droplet model are $30\%-50\%$ too low.) This difference in $W^{*}$ is attributed to the supersaturation dependence of $\gamma_{SL}$ that increases with increasing volume fraction: $\gamma_{SL}^{nucl}/(kT/\sigma^{2})=0.699,0.738$, and 0.748, for volume fractions 0.5207, 0.5277, and 0.5343, respectively ref5 (a). This again indicates a curvature effect that increases the effective interfacial free energy by factors of 1.13, 1.19, and 1.21 for the nuclei of free energy $W^{*}/kT\approx 18.8,27.7,$ and 41.3, in a reasonable agreement with the result for the LJ system.

It is worth noting that the nucleation rate data observed at small supersaturations in MD simulations are orders of magnitude smaller than the results from colloid experiments ref5 (a). (Given the extreme sensitivity of $W^{*}_{CNT}$ to $\gamma_{SL}$ and $\Delta\omega$, a possible explanation can be that the relevant properties of the colloids used in the experiments are not sufficiently close to the properties of the true hard-sphere system.) A more recent BD study by Kawasaki and Tanaka ref10 , however, indicates a reasonable agreement between simulations and colloid experiments, attributing the agreement to the use of diffusive dynamics in BD simulations, as opposed to the ballistic process in MD. In turn, in a subsequent paper ref55 , an extensive numerical study relying on three different techniques (BD, umbrella sampling, and forward flux sampling) yields similar nucleation rates for the three methods, which differ from those of Kawasaki and Tanaka, and indicate that a huge discrepancy between simulation and experiment still persists.

II.3.3 Metals

Recent excellent nucleation rate results obtained on single metal droplets by chip calorimetry ref56 are reported to be in fair agreement with molecular dynamics simulations based on embedded atom potentials and the CNT ref49 , provided that in the latter the HS relationship $\gamma_{SL}\propto T$ is adopted. While this might be a reasonable approximation for metals, apparently the oxide glass data fit better to a linear relationship $\gamma=a+bT$ ref1 (a), ref2 (a).

II.3.4 Water-ice

Owing to its importance for various branches of science (atmospheric sciences, cryobiology, climatology, etc.), nucleation of ice in undercooled water is among the best investigated nucleation problems. It has been studied by both experimental ref57 and molecular dynamics methods ref48 (b). While the experimental results from various sources show reasonable coherence for the nucleation rate, the MD results display orders of magnitude differences as a function of the chosen potential and other simulation details (see Fig. 11 in Ref. ref48 (b)). We are unaware of MD simulations that provide a complete set of data needed for a rigorous theoretical test of nucleation (specifically, barrier height from umbrella sampling, interfacial free energy, and thermodynamic driving force as a function of undercooling) for the same water potential and simulation details. As far as the experiments are concerned, there is a wealth of data for the nucleation rate of hexagonal ice $I_{h}$, and some for cubic ice $I_{c}$ at much deeper undercoolings ref57 , a range of estimated values for the interfacial free energy ref44 ; ref58 , and some accurate thermodynamic data for the undercooled state, at least for $I_{h}$ ref59 . Owing to the lack of data for the height of the nucleation barrier, only comparison of the theoretical and measured/simulated nucleation rates can be used as a test, into which uncertainties associated with the nucleation prefactor enter ref4 .

II.3.5 Structural aspects

Theoretical considerations, colloid experiments, and molecular simulations imply that structural ordering in the liquid plays an essential role in crystal nucleation. The concept of heterophase fluctuation has been enriched considerably by molecular scale studies in the past decade ref48 . The structure of the local molecular neighborhoods has been characterized in terms of bond order parameters ref60 . The importance of Medium Range Crystalline Order ref10 ; ref48 (MRCO, Fig. 2) and crystal-like precursor structures ref4 ; ref5 ; ref61 have been emphasized in crystal nucleation. Theoretical, experimental and simulation results were presented for the appearance of a disordered precursor preceding crystal nucleation in a range of systems, including 2D and 3D colloids ref62 , the LJ ref63 , and HS ref9 systems. Prediction of such complex structural phenomena is beyond the possibilities of the classical nucleation theory and represents an important challenge to more advanced theoretical approaches.

II.4 Formal theory of crystallization

The kinetics of crystallization taking place via nucleation and growth is often interpreted in terms of a simple mean field approach, the Johnson-Mehl-Avrami-Kolmogorov (JMAK) theory ref64 , which expresses the time evolution of the transformed fraction in a $d$-dimensional system as

[TABLE]

where $\tilde{Y}=K\int_{t_{0}}^{t}J(\xi)\{\int_{\xi}^{t}V(\zeta)d\zeta\}^{d}d\xi$ is the extended volume calculated by allowing a multiple overlapping of particles, $J$ is the nucleation rate, and $V$ the maximum of the anisotropic growth rate, while $K$ is a geometrical factor (volume of the $d$-dimensional particle that has unit radius in the direction of the maximum growth rate). Often transient nucleation cannot be neglected, and the formation of critical clusters starts only after an incubation time $t_{0}>0$. For constant nucleation and growth rates with an elliptical anisotropy, Eq. (5) transforms to ref40 (c)

[TABLE]

Here $\tau=[(\Omega_{d}/p)JV_{min}^{d}\Pi_{j=1}^{d}\beta_{j}]^{-1/p}$ is the characteristic time of the crystallization, and $p=1+d$ is the Avrami-Kolmogorov kinetic exponent, $\Omega_{d}$ is the volume of the $d$-dimensional unit sphere, $V_{min}$ the minimum growth rate, and $\beta_{j}=V_{j}/V_{min}\geq 1$, whereas $V_{j}$ $(j=1,\dots,d)$ are the growth rates in the direction of the principal axes of the $d$-dimensional ellipsoid. This relationship is exact if (i) the system is infinite, (ii) the nucleation rate is spatially homogeneous, and (iii) either a common time-dependent growth rate applies or anisotropically growing convex particles are aligned parallel {for derivation of Eq. (6) by Cahn’s the time cone method see Ref. ref65 }. The exponent $p$ is often evaluated from the slope of the “Avrami plot,” $\ln[-\ln(1-Y)]$ vs. $\ln(t)$, a method widely used to extract information from experiments on the transformation mechanism. Standard references ref40 (c) present $p$ values expected for different transformation mechanisms. (For example, in a 2D system, $p=3$ applies for constant nucleation and growth rates.) However, Monte Carlo simulations for randomly oriented anisotropic particles indicated a substantial deviation from the JMAK kinetics, and that under such conditions $p$ reduces with increasing transformed fraction ref66 . Apparently, $p$ alone cannot be taken as a reliable indicator of the crystallization mechanism. Advanced theoretical approaches were put forward to address these and more complex cases ref67 . The PF simulations proved useful in addressing such problems ref19 (a), ref68 , ref69 (c).

III Coarse grained phase-field models

The phase-field (PF) technique and its applications are described in a number of recent reviews ref69 . We recall only its main features that are needed for understanding the matter presented herein. The PF theory is a direct descendant of the Cahn–Hilliard/Ginzburg–Landau type classical field theoretic approaches to phase boundaries, and its origin can be traced back to a model of Langer from 1978 ref23 and developments presented later by others ref23 ; ref70 ; ref71 ; ref72 ; ref73 ; ref74 ; ref75 ; ref76 ; ref77 . In order to characterize the local state of matter, a non-conserved structural order parameter $\phi(\mathbf{r},t)$ is introduced that is termed the phase field. This structural order parameter is considered to be a measure of local crystallinity, and viewed as the Fourier amplitude of the dominant density wave representing the periodic number density in the crystalline phase. Another interpretation often used considers the phase field as the local volume fraction of the phase represented. While much depends on the details of the approach, the presence of $N$ phases can be monitored by $N$ phase fields $\{\phi_{i}(\mathbf{r},t)\}$ $(i=1,...,N)$. Some of the models, such as the multi-phase-field (MPF) theories by Chen et al. ref71 , Steinbach et al. ref72 , and others ref73 , introduce an independent phase field for every crystal grain, and work with tens of thousands of fields in describing multi-grain problems ref74 . Other more economic approaches employ orientation fields to monitor the local crystallographic orientation ref19 (a),(c), ref69 (c),(f), ref75 ; ref76 ; ref77 .

Expanding the free energy (or entropy) density of an inhomogeneous system consisting of the liquid and solid phases in terms of the structural order parameters $\mathbf{\Phi}=\{\phi_{i}\}$, and other slowly changing fields such as the chemical concentration field(s) $\mathbf{c}=\{c_{j}\}$, temperature, etc., while retaining only those spatial derivatives that are allowed by symmetry considerations, the free energy of the system can be written as a local functional of these fields, and their spatial derivatives:

[TABLE]

where spatial intergation is performed for the volume of the system. The first terms on the RHS emerge from the square gradient approximation, and penalize the spatial variation of the applied fields, giving rise to the excess free energy associated with the interfaces. The coefficient matrices $\mathbb{A},\mathbb{B},$ … denote general quadratic terms for the respective gradients: Choosing, for example, $\mathbb{A}=\mathbb{I}$ (here $\mathbb{I}$ is the identity matrix) yields a simple sum of the SG terms, $\mathbb{A}=\mathbb{I}-\mathbf{e}\otimes\mathbf{e}$ [where $\mathbf{e}=(1,1,...,1)$] yields a pure pairwise construction, whereas $A_{i,i}=\sum_{j\neq i}\phi_{j}^{2}$, $A_{ij\neq i}=-\phi_{i}\phi_{j}$ realizes the anti-symmetrized (Landau-type) gradient term. Coefficients $\epsilon_{\phi}$ and $\epsilon_{c}$ may vary with the temperature, the orientation of the interface, and the other field variables. The last term in the integrand of Eq. (7) is the free energy density, $f(\mathbf{\Phi},\mathbf{c},T,...)$, that has at least two minima: one for the bulk liquid phase, while the other(s) for the crystalline phase(s) or crystal grains. In some models, the local crystallographic orientation is represented by an orientation field $\mathbf{\Theta}$, which can be a scalar, vector, or tensor field depending on the dimensionality of the problem. To ensure the rotational invariance of the free energy, only differences of the orientation field and its spatial derivatives may appear in the free energy.

Although attempts have been made to derive the free energy functional of the crystal-liquid systems on statistical physical grounds, relying on different versions of the density functional theory of classical particles ref21 ; ref22 ; ref78 , these molecular scale approaches are usually too complicated to address complex solidification problems. As a result, phenomenological free energy (or entropy) functionals are used in the PF approaches, whose form owes much to the Ginzburg–Landau models used in describing magnetic phase transitions or phase separation ref79 . The PF models differ in the field variables considered, as well as the actual form chosen for their interaction. For example, there are models that prescribe restrictions for the sum of some of the applied fields: e.g., $\sum_{i=1}^{N}\phi_{i}(\mathbf{r},t)=1$ or $\sum_{i=1}^{M}c_{i}(\mathbf{r},t)=1$.

Once the free energy functional is defined, there are two ways to address nucleation: (i) one may solve the equations of motion (EOMs) in the presence of appropriate noise representing the thermal fluctuations and observing then nucleation, or (ii) evaluate nucleation barrier and other properties of the nucleus via solving the Euler-Lagrange equation (ELE) assuming 0 field gradient at the center of the nucleus and the undercooled liquid properties in the far-field. Since the addition of noise to the EOM may influence the thermodynamic properties (free energy minima, interfacial properties, etc.) ref80 , results from the two routes are expected to converge in the small noise limit, unless parameter renormalization is used to match the noisy and noiseless systems ref80 ; ref81 .

Making the assumption of relaxation dynamics together with the previously mentioned criteria for the sum of the fields, the time evolution of the system is described by the following equations of motion (EOMs) ref73 (g) for non-conserved and conserved fields ref82 , respectively:

Non-conserved dynamics:

[TABLE]

Conserved dynamics:

[TABLE]

Here $\delta F/\delta\phi_{i}$ and $\delta F/\delta c_{i}$ are the functional derivatives of the free energy with respect to $\phi_{i}$ and $c_{i}$, respectively. In different PF models different choices were made for the mobility matrix $\kappa_{ij}$ ref69 ; ref73 ; ref74 . {For a recent review and a physically consistent formulation see Ref. ref73 (g).} Since the free energy should decrease in any volume the mobility matrix must be positive semidefinite.

The equilibrium features, such as the phase diagram, the interfacial free energies, the nucleation barrier, etc. can be obtained by solving the Euler-Lagrange equation(s) [ELE(s)] under the appropriate boundary conditions. For a general $N$-phase field case, the multiphase Euler-Lagrange equations read as:

[TABLE]

Here $\lambda(\mathbf{r})$ is a Lagrange multiplier emerging from the local constraint of the sum of the phase fields. Eliminating the Lagrange multiplier yields

[TABLE]

Next, we illustrate the application of the PF approach to nucleation in a few simple cases.

III.1 Application of the phase-field method to nucleation

As specific examples that were used in addressing nucleation problems, we briefly outline here a few simpler PF models, where the local physical state is characterized by a single phase field ($\phi$) that is 0 in the bulk liquid and 1 in the crystal, and cases in which an additional field is used, such as concentration ($c$), number density ($\rho$), or volume fraction ($X$).

III.1.1 Single-field models

A realization of such an approach is given by the following free energy functional, EOM, and ELEs:

Free energy:

[TABLE]

where $s({\bf n})$ is an anisotropy function that depends on the components of $\nabla\phi=\{\phi_{x},\phi_{y}\}=|\nabla\phi|{\bf n}$ in two dimensions ref83 , n the surface normal, $w$ the free energy scale, $g(\phi)$ is a double well function, $p(\phi)$ is an interpolation function varying monotonically between 0 at $\phi=0$ and 1 at $\phi=1$, $\Delta f=(f_{s}-f_{l})<0$ is the thermodynamic driving force for solidification, whereas $f_{l}$ and $f_{s}$ are the grand free energy densities for the bulk liquid and solid states. In the Ginzburg-Landau (GL) approach, the double well and interpolation functions depend on the crystal structure ref84 . For various choices of $g(\phi)$ and $p(\phi)$ made in the literature, see Table I ref84 ; ref85 ; ref86 ; ref87 ; ref88 ; ref89 ; ref90 . We note that the $p(\phi)$ function in the first and fifth rows of Table I were deduced ref85 for non-isothermal problems following the thermodynamically consistent approach of Penrose and Fife ref91 .

Euler-Lagrange equation: The extremum of the free energy can be obtained as the solution of the ELE. In the case of the planar equilibrium interface ($\Delta f=0$) the integral form of the ELE applies ref16 (a):

[TABLE]

where $I$ stands for the integrand of Eq. (12). Inserting $I$, one finds that $\frac{\epsilon^{2}}{2}(d\phi/dx)^{2}=wg(\phi)$. This relationship can be used to relate the model parameters to measurable quantities: integrating the free energy density and $d\phi/dx$ from $\phi=0$ to $\phi=1$, and $(dx/d\phi)$ from $\phi=0.1$ to $\phi=0.9$, while assuming a quartic double well, $g(\phi)=\frac{1}{4}\phi^{2}(1-\phi)^{2}$, the solid-liquid interfacial free energy and the $10\%$ to $90\%$ interface thickness (along which the phase field varies between 0.1 and 0.9) can be expressed as $\gamma_{SL}=\frac{1}{6}\sqrt{\epsilon^{2}s^{2}w/2}$ and $d_{SL}=2\log(9)\sqrt{2\epsilon^{2}s^{2}/w}$, respectively, whereas the phase field profile \phi(x)=\frac{1}{2}\big{\{}1+\tanh[\log(9)x/d_{SL}]\big{\}} minimizes the excess free energy of the solid-liquid interface. Making $\epsilon^{2}\propto T$ and $w\propto T$, one finds $\gamma_{SL}\propto T$ and $d_{SL}=\mathrm{const.}$, a behavior consistent with MD simulations ref45 , and the results for the hard-sphere system.

Similarly to the classical theory, the nucleus represents here an unstable equilibrium (a saddle point in the function space), whose properties can be found by solving the ELE ref16 (b)

[TABLE]

while assuming 0 field gradients at the center for symmetry reasons, and bulk undercooled liquid properties in the far-field. Assuming spherical symmetry (a reasonable approximation for metallic systems), the ELE can be rewritten in the spherical coordinate system as

[TABLE]

where $f=wg(\phi)+p(\phi)\Delta f$, and the boundary conditions are $\partial\phi/\partial r=0$ at $r=0$, and $\phi=0$ for $r\rightarrow\infty$.

Equation of motion: As $\phi$ is a non-conserved field, an Allen-Cahn type EOM applies ref92 . For the isotropic case ($s=1$), one finds:

[TABLE]

where $\zeta$ stands for an additive noise of correlator $\langle\zeta(\mathbf{r},t),\zeta(\mathbf{r}^{\prime},t^{\prime})\rangle=2kTM_{\phi}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$ satisfying the fluctuation dissipation theorem. In the case of $p(\phi)=\phi^{2}(3-2\phi))$, for which the $\tanh$ profile is retained for steady state front propagation, the front velocity can be obtained analytically as $v_{PF}=M_{\phi}12\Delta f\sqrt{\epsilon^{2}/(2w)}$. We recall here that the steady state velocity seen in molecular dynamics simulations for the Lennard-Jones system could be fitted with a Wilson-Frenkel type expression ref93 , v=v_{0}\big{\{}1-\exp[-\Delta F_{m}/(RT)]\big{\}}\approx v_{0}\Delta F_{m}/(RT)\approx\beta\Delta T with different coefficients for the (111) and (100) interfaces: $v_{0,111}=(a\alpha D/\lambda^{2})e^{-\Delta S/k}$ and $v_{0,100}=(\alpha a/\lambda)\sqrt{3kT/m}$ ref94 . Here $\Delta F_{m}=-v_{m}\Delta f$ is the molar free energy difference, $\beta$ the kinetic coefficient, $a$ the spacing of the crystal planes, $\alpha$ the fraction of active sites at the crystal surface (taken to be 0.27), $D$ the self-diffusion coefficient, $\lambda$ the mean free path, $\Delta S$ the molecular entropy difference between the liquid and the crystal, and $v_{m}$ the molar volume. One approach to fixing the phase-field mobility is via equating $v_{PF}$ with $v_{0}\Delta F_{m}/(RT)$, yielding $M_{\phi}\approx\log(9)/(3\sqrt{2})[v_{0}v_{m}/(d_{SL}RT)]$. Another approach, which allows one to obtain an orientation dependent $M_{\phi}$, is to employ kinetic coefficients evaluated from molecular dynamics simulations ref69 (b), ref95 .

1.1 Applications for homogeneous nucleation

This simple model has been used with some success to describe crystal nucleation in single component and stoichiometric systems ref87 ; ref88 ; ref89 ; ref90 ; ref96 ; ref97 ; ref98 ; ref99 ; ref100 ; ref101 ; ref102 ; ref103 . Oxtoby and coworkers ref88 ; ref96 and Iwamatsu ref89 developed analytical solutions for the nucleation barrier on the basis of piecewise parabolic free energy and a variable-position step function $p(\phi)$. Interestingly, the density functional theory of condensation relying on Yukawa attraction can be transcribed exactly to a gradient theory using the chemical potential of the HS system as an order parameter ref97 , which together with a piecewise parabolic free energy leads to a simple but accurate analytic formulation for vapor condensation ref98 ; ref99 . A quartic double well combined with $p(\phi)=\phi$ has been adopted for describing crystal nucleation in undercooled H2O ref90 , Ar, and six stoichiometric oxide glass compositions ref14 (c). The standard PF model [for $g(\phi)$ and $p(\phi)$, see Table I] has been employed to address nucleation in the HS system ref100 , in undercooled Ar, and H2O ref19 (a) and methane hydrate ref101 . Following Stroud and coworkers ref102 , Iwamatsu and Horii put forward Ginzburg-Landau formulations for bcc and diamond cubic structures ref86 . Quantitative tests have been performed for crystal nucleation in the HS system using the SG approximation combined with Ginzburg-Landau free energies ref84 ; ref103 . It appears that in the cases where comparison has been made with simulation or experimental data, the PF approaches gave considerably more realistic results than the droplet model of the CNT ref14 (c), ref19 (a), ref84 ; ref90 ; ref96 ; ref97 ; ref98 ; ref99 ; ref103 . It also appears from studies for the HS system that better accuracy may be expected if physically motivated (Ginzburg-Landau) double well and interpolation functions are used ref84 ; ref103 .

In general, for large clusters ($R^{*}\gg d_{SL}$), the PF models predict essentially the same $W^{*}$ and $R^{*}$ as the droplet model. However, improved results are obtained for small clusters ($R^{*}\approx d_{SL}$), for which the PF result for $W^{*}$ is smaller than $W^{*}_{CNT}$ and tend to zero, when moving towards a spinodal temperature predicted by the PF theories. (Below the spinodal point the undercooled liquid becomes unstable against density fluctuations and crystallizes. Depending on the details of the PF model, the critical size may diverge or tend to zero at the spinodal. To illustrate these general features in a simple case, we employ the standard PF model for crystal nucleation in pure Ni, whose properties are summarized in Table II. ref49 (a), ref104 ; ref105 ; ref106 ; ref107 (The thermodynamic driving force of crystallization $\Delta f$ has been approximated using Turnbull’s linear relationship ref106 .) The results are presented in Fig. 3. The temperature dependence of the double well free energy is shown in Fig. 3(a), whereas the radial phase-field profiles are presented in Fig. 3(b). The nucleation barrier and the critical size are displayed in Figs. 3(c) and 3(d), together with the respective data from the CNT, in which assuming with $\gamma_{SL}=\gamma_{SL,eq}(T/T_{f})$, and from a phenomenological diffuse interface theory (DIT) ref34 (c),(d). For comparison, data from MC simulations and experiments ref49 (a) are also shown in Fig. 3(c). The PF, CNT, and DIT results appear to be in a reasonable agreement with each other, and the data from MC simulations and the experiments. For $\Delta T/T_{f}\rightarrow 1$ one observes $W^{*}/kT\rightarrow 0$. We note that the closeness of the PF and CNT results follows from the symmetric $g(\phi)$ and $p(\phi)$ functions of the standard PF model that yield a symmetric phase-field profile, known to result in a zero Tolman length ref108 , which in turn leads to an interfacial free energy independent of surface curvature. The standard PF model was also found to be in a good agreement with experimental data for crystal nucleation in undercooled liquid Ar and water ref19 (a).

Because of the importance of ice nucleation in undercooled water, we present here an updated comparison between available experiments ref57 and two models, the standard PF model and the CNT. In the latter, two cases are considered: (a) using the $\gamma_{SL}(T)=\gamma_{SL,eq}(T/T_{f})$ relationship that trivially occurs in the hard sphere system, and was close to the temperature dependence observed in MD simulations for the LJ system ref45 , and (b) a constant value $\gamma_{SL}=\gamma_{SL,eq}$. In calculating the nucleation rate, we use the classical nucleation prefactor, which was, however, multiplied with $10^{2}$ to agree with the MD results by ten Wolde et al. ref4 . The other properties were taken from Ref. ref90 , except that the polynomial fit for the specific heat difference was used down to $T_{1}=T_{f}-45$ K, where the enthalpy difference between the liquid and solid is $\Delta H_{1}=3658$ J/mol. Below this temperature transition to low density water is expected, which can be quenched into low density ice. To extend our treatment below $T_{1}$ in a thermodynamically consistent way, we follow the route described in Ref. ref90 , and use a simple model to estimate the specific heat difference between the undercooled liquid and the ice crystal on the basis of measured thermodynamic properties. Accordingly, an average excess specific heat difference of $\Delta C_{p,1}=80$ J/mol/K is assumed between temperatures $T_{1}$ and $T_{2}=T_{f}-\{\Delta H_{1}-\Delta H_{x}-\Delta C_{p,x}(T_{1}-T_{x})\}/\Delta C_{p,1}$, where $T_{x}=155$ K, and $\Delta H_{x}=1380$ J/mol are the temperature and heat of crystallization of the low density amorphous phase ref59 , whereas $\Delta C_{p,x}=3.6$ J/mol/K is the specific heat difference between liquid and the crystal above the glass transition ref59 . The predicted nucleation rates are presented as a function of undercooling in Fig. 4. Apparently, the predictions of the standard PF model and the CNT model with $\gamma_{SL}\propto T$ are in a good agreement with the experiments for small undercoolings, where the thermodynamic properties are accurate, and are yet in a reasonable agreement with them at large undercoolings, where the expected accuracy of the thermodynamic data is lower. In contrast, as seen in the case of Ni and the HS system, the CNT model with the choice $\gamma_{SL}=\gamma_{SL,eq}$ deviates considerably from the experiments.

Although the agreement between the standard PF theory and the experiments/computer simulations appear to be reasonable for both Ni and the ice-water system, we wish to note that the results are sensitive to the form of the double-well and interpolation functions, and further work is needed to sort out, which of these functions represent the best (physical) choice. A comparison of a few choices are presented for Ni in Fig. 5. Apparently, the height of the nucleation barrier is critically sensitive to the choice of the interpolation function $p(\phi)$, whereas it is far less sensitive to the form of $g(\phi)$. Remarkably, for Ni, the Ginzburg-Landau prediction for the fcc structure falls far away from the MC results. In turn, the symmetric $p(\phi)$ of the standard PF theory leads to results that are fairly close to the MC results, essentially independently of the double-well function. However, it is worth recalling in this respect that the nucleation pathway can be more complex than assumed here: the system may visit one or more nucleation precursors, in which case any observed agreement might turn out to be fortuitous. Such a phenomenon is expected to reduce the nucleation barrier, and could in principle account for the failure of the Ginzburg-Landau approach.

1.2 Treatment of heterogeneous nucleation

A foreign wall can be represented in a single-field PF theory as a boundary condition at a mathematical surface $S$ that defines the shape of the wall ref19 (b), (d), where either $\nabla\phi\cdot{\bf n}$ or $\phi$ can be specified, which in turn fixes the contact angle $\Psi$ that characterizes the wetting properties at $S$. To show how this can be done, the free energy of the system is written as a sum of a surface and a volumetric contribution, as done by Cahn ref109 :

[TABLE]

where $Z(\phi)$ is Cahn’s “surface function”, specifying the wetting properties of the surface. Minimizing $F$, yields Eq. 16 and the boundary condition

[TABLE]

Here $\bf n$ is the outward pointing normal to surface $S$, whereas $z(\phi)\equiv\partial Z/\partial\phi$. This boundary condition can be realized by setting either $-\epsilon^{2}\nabla\phi\cdot{\bf n}=z(\phi)$ (Model I) or $\phi=\mathrm{const.}\ \mathrm{e.g.}\ (\delta\phi=0)$ (Model II), at the boundary ref19 (b).

In Model I, Gránásy et al. ref19 (b) assumed that the presence of a flat wall does not perturb the structure of the equilibrium (planar) solid-liquid interface. $z(\phi)$ is then deduced as follows: Eq. (13) is employed at the melting point yielding $(\epsilon^{2}/2)(\nabla\phi)^{2}=f(\phi)-f(\phi_{\infty})=\Delta f(\phi)$, whereas the normal component of the phase-field gradient is expressed as ${\bf n}\cdot\nabla\phi=|\nabla\phi|\cos(\Psi)$, where $\Psi$ is the (contact) angle between the solid-liquid interface and the foreign surface. Combining these expressions one obtains the boundary condition ref19 (b)

[TABLE]

For $\Psi=90^{\circ}$, this relationship boils down to the no-flux boundary condition used in early PF studies ref110 ; ref111 for realizing this specific contact angle. A similar expression was proposed by Ding and Spelt ref110 (a) at essentially the same time, and an earlier work by Jacqmin implies without presenting the details that he probably developed a similar model ref110 (b). The solutions of the Euler-Lagrange equation in 2D using this boundary condition are presented in the upper row in Fig. 6, which shows that Model I can be regarded as a diffuse interface realization of the classical spherical cap model, however, with sharp wall-liquid and wall-solid interfaces.

In the case of Model II, the relationship between the phase field value $\phi_{0}=$ const. prescribed on $S$ can be related to the contact angle as $\mathrm{cos}[\Psi(\phi_{0})]=-1+6\phi_{0}^{2}-4\phi_{0}^{3}$ ref19 (b). Solving the Euler-Lagrange equation with this boundary condition, diffuse wall-liquid and wall-solid interfaces are realized that represent liquid ordering and solid disordering at the wall [see central row in Fig. 6]. A remarkable feature of Model II is the existence of a surface spinodal; i.e., there is a $\phi_{0}$ dependent critical undercooling, at which the nucleation barrier disappears [see bottom panel of Fig. 6] ref19 (b).

Simulations employing the boundary conditions of Models I and II were used to describe heterogeneous nucleation on complex surfaces (Figs. 7 and 8), including the crystallization of a liquid volume containing $N$ nanorods characterized by a contact angle of $\Psi=75^{\circ}$. With increasing $N$ and decreasing $\Psi$ the crystallization of the liquid volume accelerates (Fig. 8). Another straightforward application is modeling of particle induced crystallization (”athermal” nucleation) ref19 (b). Including a cylindrical particle of 20 nm diameter and 5 nm height, aligning its axis vertically, which is bound by wetting horizontal circular faces ($\Psi=45^{\circ}$) and non-wetting ($\Psi=175^{\circ}$) vertical sides, one can investigate the theoretical prediction ref42 that there exists a critical undercooling $\Delta T_{c}$ below which stable crystalline caps form, whereas beyond it free growth takes place. It has been found that indeed this happens qualitatively (see Fig. 9(a)) ref19 (b), however, owing to the presence of fluctuations in the equation of motion, $\Delta T_{c}$ in the simulations is about half of the theoretical prediction. Similar behavior has been reported recently in molecular dynamics simulations ref112 (Fig. 9(b)).

1.3 PF simulations of transformation kinetics

Next, we review results for crystallization kinetics obtained from simulations performed using single-field PF models in two dimensions by solving numerically the EOM [Eq. (16)], while incorporating an appropriate noise that represents the thermal fluctuations of the order parameter ref113 ; ref114 ; ref115 . Here nucleation happens automatically after an incubation period as a result of the accumulating effect of the phase-field noise. The time dependence of the transformed fraction appears to follow Johnson-Mehl-Avrami-Kolmogorov kinetics, while the kinetic exponent $p$ falls between 2.9 and 3.0 that compares well with the theoretical value $p=1+d=3$ expected for steady state nucleation and growth in two dimensions ($d=2$), (see Fig. 10). The results Jou and Lusk ref113 , Castro ref114 , and Iwamatsu ref115 reported are in a general agreement with each other and Monte Carlo simulations for steady state nucleation and isotropic growth ref66 , although Iwamatsu reported higher values for the kinetic exponent ref115 presumably because of a nucleation rate increasing with time. Different results were reported by Heo et al. ref116 , who found a time dependent growth rate for small sizes due to the Gibbs-Thomson effect. While this growth transient was observed in preceding works ref113 ; ref114 , owing to the much smaller nucleation rate in those works, it had a negligible effect on the transformation kinetics, whereas in Ref. ref116 the transient dominates. According to the study of Iwamatsu, JMAK kinetics is of limited relevance to phase transitions, in which heterogeneous nucleation plays a role ref115 , whereas Castro found that the correlation length of colored noise (of Gaussian correlator) can influence the kinetic exponent significantly ref114 . Although employing noise in the EOM to model the thermal fluctuations is an appealing approach in many respects (e.g., the transition happens automatically), the time and size scales of such simulations are rather limited.

Another way to model nucleation that circumvents this problem is to incorporate randomly distributed critical or supercritical particles into the simulation box so that they mimic the stochastic features of the nucleation process as proposed by Simmons et al. ref68 , Gránásy et al. ref19 (a), and Heo et al. ref116 , an approach that yields similar transformation kinetics as the ones based on modeling the fluctuations via adding noise to the EOM.

1.4 Nucleation in the presence of a metastable phase

A straightforward approach to modeling nucleation in the presence of metastable phases is to use a multi-well ($n\geq 2$) free energy-order parameter relationship ref96 ; ref111 ; ref116x , where the lowest well is the stable phase and the others are metastable. This approach gives rise to a splitting instability of the liquid-solid interface, leading to the formation of the metastable phase ref111 . A piecewise parabolic three-well free energy [Fig. 11(a)] was employed in a single-field PF theory (termed at that time a single order parameter Cahn-Hilliard theory) to address crystal nucleation in the presence of a metastable phase ref96 . It has been shown that besides the nucleus of the metastable phase, two types of composite nuclei exists: one with a thin interfacial layer of the metastable phase and another with a broad interfacial layer [see Figs. 11(b), 11(c), and 11(d)]. These composite nuclei converge in a bifurcation point at a critical undercooling [Fig. 11(e)].

III.1.2 Two-field models

The structural order parameter is often coupled to other slowly evolving fields. Some, like energy and mass, are subject to conservation laws and can then be further mapped through constitutive relationships to measureable quantities like temperature, density and concentration, while there are also couplings to non-conserved fields such as the orientation field, which evolves by entirely different principles. Here, we summarize the features of a simple two-field model akin to the one by Warren and Boettinger ref70 that describes isothermal solidification in a binary system; i.e., the phase field is coupled to a concentration field. This model offers a possibility to address crystal nucleation ref19 (a) and dendritic solidification ref70 .

2.1 Homogeneous nucleation in binary systems

2.1.1 Solving the Euler-Lagrange equation: Our starting point is a binary extension of the free energy functional (for constituents A and B) given by Eq. (12), in which the bulk free energy density now depends on not only the non-conserved phase-field, $\phi(\mathbf{r},t)$ and the (uniform) temperature, $T$, but also on a conserved field that specifies the local concentration of species B, $c(\mathbf{r},t)$. Following Warren and Boettinger ref70 , we incorporate here an SG term only for the phase field (strictly, this simplification is valid for an ideal solution ref16 ):

[TABLE]

where the $g(\phi)$ and $p(\phi)$ functions of the standard PF model are used, whereas $w(c)=(1-c)w_{A}+cw_{B}$, and the thermodynamic data, $f_{s}(c,T)$ and $f_{l}(c,T)$, can be taken from either databases or from the ideal or regular solution models ref19 (a). Owing to the lack of an SG term for the concentration, the respective ELE is a degenerate one ($\partial f/\partial c=0$), which defines a one-to-one relationship between the phase- and the concentration fields, $c=c(\phi)$. Assuming isotropic interfacial properties (spherical symmetry), the Euler-Lagrange equation for the phase field boils down to Eq. 15, however, now $f=f\{\phi,c(\phi)\}$. Under the conditions that define the binary nucleus, i.e., $\phi^{\prime}\rightarrow 0$ and $c^{\prime}\rightarrow 0$ for $r=|\mathbf{r}|\rightarrow 0$ (′ stands here for differentiation with respect to $r$), whereas $\phi\rightarrow 0$ and $c\rightarrow c_{\infty}$ for $r\rightarrow\infty$, $c(\phi)=c_{\infty}e^{-y}/(1-c_{\infty}+c_{\infty}e^{-y})$, where $y=v_{m}(w_{B}-w_{A})g(\phi)/R-v_{m}(\Delta f_{B}-\Delta f_{A})p(\phi)/RT$. The quantities $\Delta f_{A}$ and $\Delta f_{B}$ are the free energy density differences for the pure components at the actual temperature, which can be well approximated by Turnbull’s linear relationship ref106 for low undercoolings.

This approach has been used to address crystal nucleation in Cu-Ni, a system with nearly ideal solution behavior ref19 (a). Assuming homogeneous nucleation, $\alpha=0.6$ (close to 0.58 from MD simulations evaluated from the capillary wave spectrum at the crystal-liquid interface of Ni ref117 ), one obtains undercoolings for nucleation rates of $10^{-4}$ to $1$ drop*-1s-1* for droplets of 6 mm diameter that fall close to the experimentally observed values ref118 . We note, however, that the agreement might be fortuitous, as the accuracy of the MD result for the interfacial free energy critically depends on the accuracy of the applied potential. Another uncertainty is that amorphous precursor mediated two-step nucleation may be relevant here ref63 .

A more complex treatment is required when addressing eutectic systems ref119 , in which case two crystalline phases compete with each other. In such a case one needs to add an SG term acting on the concentration field $\frac{1}{2}\epsilon_{c}^{2}|\nabla c|^{2}$ to the free energy density in Eq. 20. Such an approach was used to investigate nucleation in the Ag-Cu system, whereas the $g(\phi)$ and $p(\phi)$ functions deduced for the fcc structure within the GL approach were used. Considering this form of the free energy functional, the respective ELEs read as follows:

[TABLE]

and

[TABLE]

where it has been utilized that $\phi\rightarrow 0$ and $c\rightarrow c_{\infty}$ for $r\rightarrow\infty$, and it was assumed that $\epsilon_{c}^{2}=\epsilon_{c}^{2}(\phi,c)$. The latter assumption was made to connect $\epsilon_{c}^{2}$ to the phase- and concentration dependent interaction coefficient of Ag-Cu in the spirit of Ref. ref16 (a). The thermodynamic data were taken from the database ThermoCalc ref105 .

The radial phase-field and concentration profiles and the free energy of formation for crystal nuclei in pure Ag, Cu, and Ag50Cu50 liquids are displayed in Fig. 12 together with experimental undercooling data ref120 . It was found that at $c=0.5$ a Cu rich nucleus forms. The map of the nucleation barrier as a function of temperature and liquid composition (Fig. 13) indicates two types of nuclei forming: Ag rich on the Ag side and Cu rich on the other side of the phase diagram. This is consistent with the results of time dependent phase-field simulations based on solving the respective equations of motion (see Fig. 9 in Ref. ref121 ). A far more complex behavior was reported for the metastable liquid miscibility gap occurring in the phase diagram at high undercoolings: various types of nuclei compete (see Fig. 14). These include types that have a solid core surrounded by a liquid layer of composition, which differs from that of the initial liquid, indicating that here crystal nucleation is coupled to liquid phase separation ref119 . In the neighborhood of the critical point an enhanced nucleation rate is observed ref119 , in agreement with experiment ref122 , density functional computations ref123 , and MD simulations ref124 .

2.1.2 Solving the equations of motion: The first two-field PF simulations of nucleation in eutectic systems (relying on phase- and concentration fields) were performed by Elder and coworkers ref125 for a model system of symmetric phase diagram. They added noise to the EOMs that satisfies the fluctuation-dissipation theorem. At the eutectic composition competing nucleation of the two solid solution phases were observed, whereas at asymmetric compositions nucleation was dominated by the solid solution of the majority component, with the minority phase nucleating on the surface of the growing majority phase (Fig. 15).

The interaction between nucleation (represented by supercritical particles of one of the solid phases, which were included into the simulation box ahead of the front of the other phase) and peritectic growth has been studied by Nestler et al. ref126 (a). A similar approach was used for solidification in monotectic systems, however, with two-phase liquid ahead of the front (Fig. 16) ref126 (b). An extension of Simmons’ treatment of nucleation to systems described by coupled non-conserved and conserved fields was put forward by Jokisaari et al. ref126 (c).

Patterns closely resembling experimental observations were reported in a generic peritectic system ref127 . Here a heterogeneous mechanism was assumed represented by particles placed spatially randomly at the solidification front when the undercooling surpassed a critical value. The mean separation between heterogeneous nuclei were treated as an adjustable parameter. Varying the separation of nuclei transition between bands and islands were observed (Fig. 17) ref127 .

A drawback of these models is that in them the two solid phases are distinguished only by their composition, and in principle no structural difference is incorporated (the same phase field applies to both), and the free energy of the interphase boundary is exclusively of chemical origin. One might also wish to address the solidification of multicomponent systems, or polycrystalline matter composed of differently oriented crystallites. Such problems can only be addressed in the framework of PF models by introducing additional fields, which, however, necessarily complicates the treatment of crystal nucleation.

2.2 Heterogeneous nucleation in binary systems

Binary generalizations of the boundary conditions representing a foreign wall is straightforward, when no SG term is incorporated for the concentration field in the free energy density. (See Pusztai et al. ref19 (c) and Warren et al. ref19 (d).) Under such conditions the general form of the boundary condition for the phase field that sets the contact angle to $\Psi$ at the surface of substrate is as follows ref19 (c):

[TABLE]

where $c=c(\phi)$ is the implicit solution from the degenerate ELE for the concentration field, whereas to realize a chemically inert foreign wall, one needs to prescribe a no-flux boundary condition at surface $S$ of the wall. These boundary conditions can be combined with both ELE and EOM with comparable effect. We note, however, that in off-equilibrium states (in undercooled/supersaturated liquid), this boundary condition sets a contact angle that only approximates $\Psi$. Examples are available in Refs. ref19 (c) and ref19 (d).

2.3 Nucleation with density change

Two-field models analogous to those used in binary systems have been proposed to address crystal nucleation in the hard-sphere system ref100 ; ref103 , except that the phase-field is now coupled to the particle density, which in turn controls the driving force of crystallization. This is an especially interesting case, as all the properties needed for fixing the model parameters are available from MD/MC simulations (see e.g., a recent collection of them in Ref. ref103 ), together with the height of the nucleation barrier ref5 . Relying on such data, one finds that the choice of the double well and interpolation functions is crucial to determining the height of the nucleation barrier. Ignoring, for the moment, that precursor structures also play a crucial role ref9 , and using the best available MD data for the hard-sphere system to fix all the model parameters, one finds that the PF calculations performed with $g(\phi)$ and $p(\phi)$ functions from the Ginzburg-Landau approach for the fcc structure (with or without coupling to the density field), fall relatively close to the results from umbrella sampling ref5 (see Fig. 18), together with the prediction of a simple phenomenological diffuse interface theory (DIT) ref34 (c), ref128 . In contrast, the single- and two-field PF models with the standard choice of the $g(\phi)$ and $p(\phi)$ functions envelope the predictions of the droplet model of the CNT and the self-consistent CNT (SCNT). In the latter $W^{*}_{SCNT}=W^{*}_{CNT}-W_{1,CNT}$ ref129 , where $W_{1,CNT}$ is the free energy of the monomer in the classical droplet model. Note that as pointed out by Auer and Frenkel $W^{*}_{CNT}$ falls far below the MC results ref5 .

It is worth noting in this respect that in the MC simulations of Schilling et al. ref9 , in which dense amorphous clusters were observed that acted as precursors for crystal nucleation, the volume fraction of the initial liquid was about $\phi_{L}=0.54$, i.e., somewhat larger than in the simulations of Auer and Frenkel ref5 . Since the appearance of the amorphous precursors may be strongly dependent on the supersaturation, one cannot be sure whether they influence the $W^{*}$ values from MC simulations in the volume fraction range shown in Fig. 18. Further MC investigations are required to settle this issue.

2.4 Viewing the wall as a second solid phase

Another possibility for a two-field representation of heterogeneous nucleation within the PF theory is via introducing an extra field for the foreign wall (substrate), $\phi_{W}$. Instead of using Eq. (17) with the surface function acting at surface $S$, one may extend the integrals over the volume of the inert wall bounded by $S$, incorporating thus both the solidifying substance and the wall material into the modeling space, a procedure that yields ref19 (d):

[TABLE]

where $|\nabla\phi_{W}|$ is a Dirac $\delta-$function that locates $Z(\phi)$ to surface $S$, while the new factor $1-\phi_{W}$ locates the free energy density for the liquid and solid phases to those regions, where $\phi_{W}=0$. Computing then $\delta F/\delta\phi$ yields the following EOM for $\phi$,

[TABLE]

where it has been utilized that $\mathbf{n}=\nabla\phi_{W}/|\nabla\phi_{W}|$, an expression that is in some sense “obvious”, as we added the Model I boundary condition multiplied by a $\delta$-function to the original variation over the volume bounded by the inert wall. Thus, introducing the auxiliary field $\phi_{W}$, the computation can be performed over all the space, and one does not need to impose the boundary conditions explicitly at the wall. Near equilibrium, contact angles obtained using this approach for a diffuse interface wall are shown in Fig. 19 ref130 .

III.1.3 Nucleation in PF models with three or more fields

Introducing additional fields, both the capabilities and complexities of PF models expand substantially: one is able to describe competing crystalline phases, multicomponent alloys, and polycrystalline structures. An especially interesting extension is the inclusion of an orientation field (a scalar field in 2D, the rotation tensor field or the quaternion field in 3D), which monitors the local crystallographic orientation, and offers another route to model polycrystalline structures.

3.1 Nucleation in multi-phase-field models

Polycrystalline solidification and grain coarsening is traditionally addressed in the framework of multi-order-parameter (MOP) ref71 or multi-phase-field (MPF) ref72 , in which an individual phase field is assigned to each orientation distinguished in the simulation. In such models, nuclei are incorporated “by hand”, i.e., supercritical seeds of random orientation and position are placed into the simulation window, either according to the computed nucleation rate {where in Eqs. 3 or 4, one may use either the solution of the respective ELE for $W_{hom}^{*}$ ref19 (a), ref121 or for $W_{het}^{*}$ ref19 (b)-(d), or take the nucleation barrier from other sources ref68 }. Alternatively, one may incorporate particles into the simulations and activate them in agreement with a criterion from Greer’s free growth limited model ref42 , as done in a number of simulations ref131 . The MPF approach was also used to address nucleation on curved surfaces ref132 .

3.2 Nucleation of competing crystalline phases

Computation of the properties by solving the respective ELE have been so far performed for modeling phase selection in the Fe-Ni system ref133 . A three-phase MPF model relying on three phase-fields $\phi_{i}$ (where $i=1,2,3$) coupled to a concentration field $c$ was developed on the basis of Ginzburg-Landau free energies of the liquid-fcc, liquid-bcc, and fcc-bcc subsystems. The bulk free energy density $\Delta f(\phi_{i},c)$ was supplemented with an antisymmetric differential operator term ref72 (a),(b).

[TABLE]

For the fcc-bcc transition the order parameter is related to the magnitude of Bain’s distortion ref134 and leads to the same $g(\phi)$ and $p(\phi)$ functions as those deduced for the liquid-bcc system. The thermodynamic data were taken from a CALPHAD (CALculation of PHAse Diagrams ) type assessment ref105 . Other required data were the interfacial energies and thicknesses in phase equilibria for the pure components. These were taken from molecular dynamics simulations (see Ref. ref134 ). The model was applied to map the properties of the nuclei as a function of composition, temperature, and structure. Typical radial field profiles, predicted with a realistic choice of model parameters, are shown in Fig. 20(a), (b) for three cases denoted by symbols in Fig. 20 (c). Remarkably, a significant bcc layer was predicted at the surface of fcc nuclei as reported for the Lennard-Jones system by MD simulations ref4 and density functional theory ref22 (a). Model predictions for the fcc-bcc phase-selection boundary also fall close to the experimental results ref135 for the Fe-Ni system (see Fig. 21).

3.3 Nucleation in the presence of orientation field

The coupling of the phase- and concentration fields to an orientation field that represents the local crystallographic orientation opens up new directions in modeling polycrystalline solidification. This idea has been explored in a number of works addressing polycrystalline solidification and grain coarsening ref69 (c),(f), ref75 ; ref76 ; ref136 ; ref137 . In this review we address only those aspects of the orientation fields that are related to crystal nucleation. To illustrate this approach in 2D, first we supplement the free energy density in Eq. 20 with an orientational contribution $f_{ori}=hp(\phi)|\nabla\theta|$, where the orientation energy scale is $h\propto T$, and $\theta(\mathbf{r},t)$ is a scalar field representing the local orientation angle relative to the laboratory frame. Accordingly, $\theta$ is a circular variable. As $\theta$ is a nonconserved field, the respective equation of motion reads as ref136 (c)

[TABLE]

where $M_{\theta}$ is the orientational mobility that sets the time scale of orientational ordering, and is proportional to the rotational diffusion coefficient of the molecules of the system ref136 (c), whereas the second term differs from 0 if the interfacial free energy is not isotropic $[s(\theta)\neq 1]$. Behavior of this type of EOMs has been addressed by Kobayashi and Giga in 1D and 2D, providing analytical solutions for testing ref138 . Eq. (III.1.3) tends to reduce the orientational differences in the system. A distinct possibility is to add randomly oriented supercritical seeds to the simulation that are distributed randomly in space and time ref76 (a), ref19 (a), yet a more physical approach that creates crystallites via adding noise to the equations of motion is also possible ref19 (a). Herein, we recapitulate the latter, which has been used to model a broad range of polycrystalline structures ref69 (c),(f), ref137 . For this, the orientation field is extended to the liquid, where it is made to fluctuate in time and space. To accomplish this, a noise of correlator $\langle\zeta_{\theta}(\mathbf{r},t),\zeta_{\theta}(\mathbf{r}^{\prime},t^{\prime})\rangle=\zeta_{\theta,0}^{2}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$ is added to the right hand side of Eq. (III.1.3), where $\zeta_{\theta,0}$ is the respective noise strength. To avoid double counting of the orientational contribution in the liquid, which is in principle incorporated into the bulk free energy of the liquid, $f_{ori}$ contains a multiplier $p(\phi)$, so it acts only in the solid (providing grain boundary energies) and the solid-liquid interface (where crystalline ordering takes place). Finally, since we are interested here in nucleation of grains in the undercooled liquid or at the solidification front, and not in grain coarsening of polycrystalline structures, where the latter happens on a much longer timescale than the former, $M_{\theta}=M_{\theta,0}[1-p(\phi)]$ is assumed, where $M_{\theta,0}$ is the orientational mobility of the liquid. This means that grain boundary motion is blocked after freezing. (Of course, this restriction can be relaxed, and similar orientation field models can be used to address grain coarsening ref137 .)

To illustrate how noise induced nucleation works in the orientation field model three simulations are shown. Snapshots of crystal nucleation in an undercooled liquid alloy of ideal solution (Cu-Ni) thermodynamics and faceted solid-liquid interfacial free energy of fourfold symmetry, are presented in Fig. 22, whereas in Fig. 23 the time evolution of the composition map in the case of two-step crystal nucleation in a phase separating liquid (a regular solution approximant of the peritectic Ag-Pt system) is displayed. Finally, a sequence of snapshots of eutectic solidification under conditions ref139 (a) that correspond to repeated laser melting (mimicking the thermal history during laser additive manufacturing) is shown in Fig. 24. Here the regular solution approximation was used for the Ag-Cu system, with a physical interface thickness ($\sim 1$ nm), in a constant temperature gradient. $f_{ori}$ of the model employed here ref139 (b) ensures a fix misorientation between the two solid phases. The morphology occurring in the simulation is characteristic to such processes ref139 (c), and includes alternating layers of large scale equiaxed eutectic grains (with radial lamellae) and small scale equiaxed grains (droplets of the nucleating phase surrounded by the matrix of the other phase).

The orientation field approach has been generalized to three dimensions ref9 (c), ref136 (d), ref140 . There are different mathematical representations of the crystallographic orientation in 3D, such as Euler angles, rotation matrices, Rodrigues vectors, quaternions, etc. ref141 . Here, we present an approach developed in the quaternion formalism ref19 (c), ref136 (d), ref142 . A 3D generalization of the 2D orientation field model described previously differs in $f_{ori}$, which now takes the form

[TABLE]

where $q_{i}$ are the components of the quaternion field. This definition of $f_{ori}$ boils down to a reasonably accurate approximation of the 2D model described above, provided that the orientation of the grains differs in only the angle of rotation around the normal of a plane. The EOMs obtained assuming an isotropic interfacial free energy read as ref19 (c), ref136 (d), ref142

[TABLE]

where Gaussian white noises of amplitude $\zeta_{i}=\zeta_{S,i}+(\zeta_{L,i}-\zeta_{S,i})[1-p(\phi)]$ were added to the EOMs, a form that ensured that the quaternion properties of the $q_{i}$ fields were retained. (Here $\zeta_{L,i}$ and $\zeta_{S,i}$ are the noise amplitudes in the liquid and solid phases, respectively.) While the main application of Eq. (29) has been the modeling of polycrystalline freezing ref19 (c) ref136 (d), ref142 , a similar quaternion based model was used for addressing grain coarsening ref143 .

III.2 Numerical methods

The Euler-Lagrange equations (ELEs) and the equations of motion (EOMs) employed in the PF models for crystal nucleation are fairly complex ordinary or partial differential equations, and analytical solution are available in only exceptional cases, such as the systems with piecewise parabolic free energies ref88 ; ref89 ; ref96 ; ref98 ; ref99 . Even there, implicit equations govern the matching of the analytical solutions corresponding to the individual parabolic segments of the free energy, which need to be solved numerically (by e.g. Newton-Raphson iteration ref144 ).

The boundary value problems for 1D or 2D ELE in single-field PF models (ordinary differential equations) are usually solved by combining the Runge-Kutta methods ref144 ; ref145 with a shooting algorithm ref144 , or employing a relaxation method ref144 . A few examples are given for the application of these methods in Refs. ref14 (c), ref19 ; ref100 ; ref101 ; ref103 . In solving the sets of ELEs for boundary value problems of multi-phase-field computations in 1D or 2D, the relaxation method was used ref133 . In $d\geq 2$, the use of an unstructured grid ref19 (b)-(d) is advantageous for avoiding lattice anisotropy, which emerges when using ordered grids (e.g., square or cubic).

Besides solving the ELE, the nucleation barrier can also be explored using surface walking methods that find the saddle point at the nucleation barrier starting from an initial state, and the path finding methods such as the nudged elastic band, string, and minimax methods that find the minimum energy path. A review of such methods can be found elsewhere ref146 . Path finding methods were used by several authors to address homogeneous and heterogeneous nucleation ref132 ; ref147 .

The EOMs used in PF simulations of solidification are usually sets of coupled nonlinear stochastic partial differential equations that are solved using standard numerical methods like explicit, semi-implicit, or implicit time stepping combined with finite difference ref144 ; ref145 ; ref148 , finite element ref149 , or spectral discretization ref144 ; ref150 ; ref151 . The importance of the proper choice of numerical methods (e.g., spectral spatial discretization) was nicely demonstrated in Ref. ref152 . Methods for producing uncorrelated/colored noise for non-conserved and conserved fields can be found in Ref. ref153 . Application of an adaptive grid can significantly reduce the computational cost and the lattice anisotropy ref154 , however, it does not help time stepping. It is also less efficient when nucleation is modeled via adding noise to the EOMs, as in this case high spatial resolution is needed everywhere.

III.3 Nucleation vs. microstructure evolution

Following Gránásy et al., ref69 (f), we present a few examples for phase-field modeling of polycrystalline morphologies, in which nucleation played a crucial role.

III.3.1 Transformation kinetics in alloys

Let us start by recalling that within the formal theory of crystallization (JMAK kinetics), in infinite systems, constant nucleation and growth rates are represented by the following values of the Avrami-Kolmogorov kinetic exponent: $p=3$ in 2D, and $p=4$ in 3D. A compilation of $p$ values expected for different transformation mechanisms are available in Ref. ref40 (c). The JMAK approach is widely used in different branches of sciences including materials science, chemistry, geophysics, biology, cosmology, etc. However the Avrami-Kolmogorov exponent has a limited validity as an indicator of the transformation mechanism. For example, theoretical ref67 (a),(d),(e),(f) and numerical ref66 studies conclude that in the case of anisotropic growth (e.g., elliptical crystals) $p$ decreases with increasing transformed fraction $Y$ due to a multi-level blocking of impingement events. Another essential class of transformations that deviate from the JMAK behavior, is the case of “soft impingement”, where crystal grains interact with each other indirectly via their diffusion fields. Although in the latter case, handbooks ref40 (c) assign a $d/2$ contribution to $p$ from $d$-dimensional growth, this often proves to be a crude approximation. Various approximate treatments were proposed for soft impingement ref67 (b),(c), ref155 . Numerical simulations based on the PF approach can incorporate both anisotropic growth and diffusion controlled front propagation in a natural way, and are expected to address even cases dominated by complex solidification morphologies such as dendrites. We explore a few examples, where PF modeling recovers the texbook result, while in more complex cases fitting of JMAK kinetics to the simulations is only possible by introducing a time- or transformed fraction dependent kinetic exponent.

(i) Noise induced nucleation and impinging polycrystalline spherulites that form under almost perfect solute trapping conditions (i.e., composition of the liquid and solid were close) yield an almost perfect fit to JMAK kinetics, with $p=3.04\pm 0.02$ falling close to the theoretical expectation ($p=3$) ref136 (c).

(ii) Nucleation and growth of elongated needle crystals in 3D led to an Avrami-Kolmogorov exponent $p$ that decreases with increasing $Y$ (Fig. 25) ref136 (d). Analogous phenomenon was reported in 2D MC simulations ref66 ; ref67 (d).

(iii) Although dendritic growth is controlled by diffusional instabilities, the dendrite tip is a steady-state solution of the diffusion equation. As a result, well developed equiaxed dendritic structures (with side arms filling the space between the main arms) may have steady state growth in a roughly self-similar way, so that the average of the concentrations of the solid and the interdendritic liquid domains is about that of the initial liquid, i.e., solidification without long range diffusion takes place as in the case of eutectic solidification. A 2D example of continuously nucleating dendritic structures is shown in Fig. 26(a), for which $p\approx 3$ was observed ref19 (a), ref69 (c), ref121 , consistent with the theoretically expected $p=1+d$ ($=3$ here). For a similar growth morphology, however, with a constant number of seeds, PF simulations tend to the theoretically expected value $p\approx 2(=d)$ ref156 . Apparently, however, one can produce dendrites, whose side arms are underdeveloped or missing, in which case exponent $p$ decreases with increasing $Y$ ref156 as in the case of elliptical needle crystals, ref66 indicating that various levels of blocking effects dominate the value of $p$. Another possibility for obtaining non-JMAK behavior is found in situations where nucleation dominates solidification, and the steady-state growth stage is not achieved for the majority of the dendrites, a situation that was studied in three dimensions ref19 (c), ref157 . If the particles interact via their diffusion fields, $p=1+d/2=2.5$ is expected in 3D ref40 (c). However, this value is valid only for compact growth morphologies, whereas with transition towards a fully developed dendritic structure growing in a nearly self-similar manner, $p\rightarrow 1+d=4$ is expected, i.e. the kinetic exponent shall be in the range $2.5\leq p\leq 4$. Indeed the values $p=2.99\pm 0.01$ ref19 (c) and $p=3.21\pm 0.01$ ref157 reported for large scale simulations of multi-dendritic solidification fall into this range (Fig. 27)(a). The larger the nucleation rate, the more compact the crystallites, and the closer is the interaction of the particles to the diffusion-controlled mechanism. Apparently, soft impingement is expected to reduce $p$ with an extent increasing with increasing $Y$; an expectation supported by theory ref158 and PF simulations ref69 (c). (Analogous phenomena were observed in PF simulations in 2D ref121 .) Kinetics of crystallization in thin films was also investigated using a 3D model ref157 . The Avrami–Kolmogorov exponent observed was, $p=2.37\pm 0.01$ ref157 . It falls between the 2D limits $p=1+d/2=2.0$ that corresponds to steady-state nucleation and fully diffusion controlled growth, and $p=1+d=3.0$, which applies for steady-state nucleation and growth rates. These limits are justifiable for thin films, as long as the thickness of the film is small relative to the size of the crystallites. Figs. 26(b) and (c), and Fig. 27)(b) show experimental images ref158a ; ref158b ; ref158c of microstructures that are qualitatively similar to the respective panels (a), however, the staistics is not sufficient to evalute the kinetic exponent $p$.

Summarizing, in the case of binary alloys, only under rather specific circumstances (e.g., fully developed dendrites, or close to full solute trapping) does one observe a constant kinetic exponent, whereas in cases of elongated crystals or diffusion controlled growth of compact particles $p$ decreases with the transformed fraction.

III.3.2 Large-scale modeling of polycrystalline solidification morphology

Recently extremely large scale ($4096\times 4104\times 4096$ voxels) binary phase-field simulations were performed on the TSUBAME 2.0 hybrid supercomputer (16,000 CPUs and 4,000 GPUs) at the Tokyo Institute of Technology, in which crystallization was started by placing randomly oriented seeds on a surface. The evolving polycrystalline dendritic growth morphology resembles closely the morphology observed in the experiments (c.f. two blocks of Fig. 28) ref159 ; ref159x .

III.3.3 Columnar to equiaxed transition

During the columnar to equiaxed transition (CET) the microstructure changes from elongated columnar growth morphology formed by directional solidification to equiaxed morphology controlled by nucleation (Fig. 29(a)), a phenomenon that influences the properties of the cast alloy . This phenomenon is described by the phenomenological model of Hunt ref160 in terms of parameters (nucleation undercooling, undercooling at the dendrite tip, and density of equiaxed grains) that are difficult to quantify. PF modeling was used to address CET more than a decade ago, however, with an approach that was unable to handle multiple orientations ref161 . The authors made a few simplifications, such as using the same crystallographic orientation for all the grains, and placing the nuclei on a crystal lattice by “hand”.

In a subsequent study that relied on an orientation field based PF model ref162 , these simplifications were removed, and quantitative computations were made for CET in Al-Ti alloys. For this purpose, the 2D and 3D versions of the orientation field model were combined with Kim’s quantitative PF model of multicomponent systems ref162 . In the simulations the matter traveled from right to left with a pulling velocity $v\in[4\times 10^{-4},32\times 10^{-4}]$ m/s in a linear temperature distribution of gradient $G\in[5\times 10^{3},4\times 10^{4}]$ K/m. At the upper and lower sides of the simulation box periodic boundary conditions were assumed. New $\Delta x$ thick layers of given composition and temperature entered the simulation box on the right hand side with timing consistent with $v$, while other columns were shifted left, and the leftmost column left the simulation box on the left hand side. Particle induced nucleation was assumed, which was approximated by the free-growth-limited model of Greer et al. ref42 . Foreign particles of a Gaussian size distribution with a mean diameter of $2R_{p}=(29\pm 3)$ nm were placed randomly into the new layers of thickness $\Delta x$, which were then activated (the phase field in the pixel or voxel was flipped to 1) once the local undercooling surpassed the critical undercooling $\Delta T_{c}=2\gamma_{SL}/(\Delta s_{f}R_{p})$. The thermodynamic properties of Al-Ti were taken from Ref. ref164 . Then, the parameters of Hunt’s model were evaluated from the simulations. The theoretical curve marking the CET are shown in Fig. 29(b). The microstructure evolution was mapped in two dimensions under conditions marked by the symbols in Fig. 29(b), whereas the respective results are displayed in Fig. 30. The predicted morphologies are consistent with Hunt’s analytical model: nucleation-controlled equiaxed structures appeared for the eight runs above the curve from Hunt’s model, whereas columnar dendritic structures were seen for the rest. Similar results were obtained in 3D (see Fig. 31) ref162 ; ref165 . A reasonable qualiative agreement can be seen between the predicted and experimental 3D columnar dendritic microstructures, thought the experimental result refers to an Al-Cu alloy ref165a .

III.3.4 Polycrystalline growth

A broad range of materials show “polycrystalline growth,” defined by solidification where new crystal grains nucleate at the perimeter of the growing crystal, a phenomenon termed growth front nucleation (GFN). This mechanism leads to the formation of spherulites, crystal sheaves, quadrites, disordered dendrites, and polycrystalline fractallike morphologies. Polycrystalline growth occurs in technical alloys, polymers, biopolymers, liquid crystals, minerals, crystallizing food products (e.g, fat, chocolate, etc.), and in biological systems including cholesterol crystals in the arteries, kidney stones, amyloid plaques associated with Alzheimer’s disease, etc. Orientation field based PF models successfully addressed a range of polycrystalline growth morphologies, such as “dizzy dendrites” ref136 (a), and a broad range of spherulitic structures and crystal sheaves, including quadrites ref136 (b)-(d). Recent reviews discuss these methods and their achievements ref69 (c),(f). Therefore, we present herein only a brief summary of various mechanisms proposed to lead to GFN.

(i) Inclusion of foreign particles into the growing crystal may lead to the formation of new grains, e.g., may yield “tip deflection” in the case of dendrites ref136 .

(ii) Decreasing the temperature in molecular liquids, the rotational diffusion coefficient decreases faster than the translational one, decreasing the ratio $M_{\theta}/M_{\phi}$. As a result, orientational ordering slows down relative to front propagation, leading to the formation of orientational defects (that can be interpreted as bundles of dislocations) freezing into the crystal ref136 (b)-(d).

(iii) Fixed angle branching of needle crystals in directions of low grain boundary energies ref136 (c).

(iv) MD simulations indicate that crystal nuclei often contain multiple orientations due to the formation of stacking faults, multiple twinning, etc. ref48 ; ref166 .

(v) Recent large scale MD simulations for Fe indicate that the probability of nucleating new crystals is higher in the vicinity of the crystallization front, leading to “satellite crystals,” a phenomenon attributed to a higher concentration of atoms of icosahedral neighborhood ref167 .

(vi) Amorphous cores can lead to the formation of polycrystalline growth forms (spherulites) ref168 .

Most of these mechanisms have been captured by orientation field based PF simulations ref136 . For example, they recovered the main features of the experimental orientation distribution in spherulites ref169 (see Fig. 32). To illustrate further the abilities of the orientation field based PF models in describing spherulitic solidification, we present a simulation in which the single crystal has a needle crystal morphology (see Fig. 33) ref169a . A continuous transition from spherulites to simpler polycrystalline objects is observed. We note, however, the micromechanisms realizing (i) to (iii) need to be further clarified, a task that requires a molecular scale approach. The progress of molecular scale PF models made in this direction is reviewed in Section IV.

III.4 Discussion: Nucleation in coarse grained PF models

After reviewing the application of phase-field models of various complexity for nucleation related problems, we attempt to answer the question whether or not the coarse grained-phase field models of crystal nucleation can be made quantitative. As pointed out earlier, information available for such tests is limited to only a very few systems. Below, we recall the relevant results, and try to use them to construct a coherent picture. Unfortunately, this proves to be a non-trivial task.

First, we address the single-component nucleation problems. We compare the results obtained for Ni and water in Section III.A.1.1 with those for the hard sphere system in Section III.B.2.3. It is remarkable that without adjustable parameters, for the HS system the PF model based on the Ginzburg-Landau free energy for fcc structure yields reasonable agreement with the MC simulations for the nucleation barrier, whereas the CNT and the standard PF model predict about half of the $W^{*}_{hom}$ from MC simulations. Unfortunately, in the case of Ni, the opposite relationship is seen: the results from the CNT (with $\gamma_{SL}\propto T$) and the standard PF model are in fair agreement with the MC results, whereas the Ginzburg-Landau (GL) free energy based PF model overestimates the nucleation barrier considerably. On one hand, the applicability of the standard PF model is further supported by the fact that it was found to be in fair agreement with experiments for crystal nucleation in undercooled liquid Ar and water. On the other hand, results for phase selection in the Fe-Ni system via heterogeneous nucleation seem to support that, at least from this viewpoint, an MPF model with GL free energy is in a qualitative agreement with the experiments. Summarizing, it appears that different versions of the PF models appear successful for different materials, while the free energy functional often has a plausible, though somewhat ad hoc, form.

It is not easy to resolve these apparent contradictions. We note that in the HS system, the density is the control parameter that determines the driving force of crystallization, whereas temperature is not a relevant parameter, e.g., the solid-liquid interface free energy depends on temperature in a trivial way, $\gamma_{SL,eq}\propto T$, just as assumed in the CNT version for Ni, that agrees reasonably with the respective MC simulations. As mentioned earlier, this temperature dependence is further supported by results for the Lennard-Jones system. In turn, in the HS system it was found that the interfacial free energy deduced from the nucleation barrier using the classical droplet model increases with density. In describing real liquids, one needs to consider both effects. In principle, both effects are incorporated into the PF model based on GL free energy: the density dependence appears to be correct, as this is tested in the case of the HS system. The temperature dependence of the equilibrium solid-liquid interfacial free energy is reproduced by making the coefficients of the SG term and the double well term proportional to $T$ ($\epsilon^{2}\propto T$ and $w\propto T$), however, the latter proves insufficient to recover the Ni data, raising the possibility of a different temperature dependence for these terms.

It is worth noting that there can be various explanations for the contradictions described above. For example, (a) one of the interfacial free energies used in computing $W^{*}_{hom}$ is significantly away from its true value, (b) the temperature dependence assumed for the interfacial free energy of Ni is inaccurate. Further uncertainty may come from (c) a possible presence of a two-/multi-step nucleation mechanism, in which crystal nucleation is assisted by one or more metastable precursor structures (either amorphous or crystalline) that may reduce the effective nucleation barrier significantly, a phenomenon that might be missed by umbrella sampling. Finally, it is possible that (d) the application of coarse grained models on the nanometer scale may be misleading, as well as the sharp-interface droplet model of the CNT.

Considering the contradictory results the PF models lead to, it is somewhat ironic that a simple analytical diffuse interface theory (DIT) ref34 (c),(d) ref128 seems to work fairly well for both the hard sphere system and Ni (and for a range of condensation/solidification problems ref128 ). Especially, that a density functional study for nucleation in vapor condensation suggests that the basic assumptions of the DIT are not satisfied, yet the predicted barriers appear to be quite accurate in the experimentally interesting regime ref170 . This happens in spite of the finding that the curvature dependence the DIT predicts for the interfacial free energy is quite away from the behavior MC simulation yield for the lattice gas model ref171 .

It is desirable to have a single approach that captures the essential features of nucleation in most cases. To handle some of these problems one would need atomic scale modeling of the nucleation phenomenon. Treatment of nucleation on the basis of ab initio (quantum mechanical) computations is infeasible, and remains so in the foreseeable future. The next best choice appears to emerge from statistical physics: the density functional theories of classical particles ref21 ; ref22 ; ref78 achieved some accuracy in predicting the properties of the HS system ref78 (e),(f). Yet, these methods are rather complex and their application to complex nucleation problems is not straightforward. A newcomer in this field is a relatively simple phase-field approach that works on the molecular scale, termed the Phase-Field Crystal model ref24 ; ref25 (PFC), which has been employed recently for exploring structural aspects of crystal nucleation ref172 . It may also serve as the basis for deriving a physically motivated coarse grained approach. Results achieved in this area are briefly reviewed in the following section.

IV Molecular scale phase-field models of crystallization

A sequence of approximations that lead to different levels of continuum models of crystalline freezing is visualized in Fig. 34 ref25 . In the classical density functional theories the local state is characterized by a time averaged particle or number density, whereas in the molecular scale phase-field models termed the phase-field crystal (PFC) theories, a gradient expansion emerging from a Taylor expansion of the direct correlation function in the Fourier space is employed to reduce the mathematical complexity of the free energy functional. For symmetry reasons only even powers of the wave number are allowed ref102 . Varying the number of the retained expansion terms, a range of PFC type models can be deduced ref24 ; ref25 . Remarkably, to achieve a reasonable reproduction of the Percus-Yevick analytical direct correlation function of the hard-sphere liquid up to its second peak the expansion needs to be extended beyond the twelfth order in $k$, corresponding to a term proportional to $\nabla^{12}\psi$ in the free energy density. With the increasing number of terms, one needs to define the value of respective coefficients. While all these models can be used to address crystal nucleation in undercooled liquids, herein we concentrate only on those that were in fact used for such purposes.

IV.1 The Phase-Field Crystal approach

Since the PFC models are newcomers to the field of nucleation modeling, we briefly review the basics of this type of approach. Here crystal structure plays a central role. The form of the free energy can be adjusted to realize different crystalline structures ref25 . We restrict this review only to those models that were in fact used in describing crystal nucleation. Dynamics is another essential issue; accordingly, we recall results obtained with models that rely on either diffusively or hydrodynamically controlled modes of density relaxation. In these models, as in the density functional approaches, the local state is characterized by time averaged particle densities that might be related to MD simulations under appropriate conditions ref173 .

IV.1.1 The single-mode PFC model

This version of the PFC model was derived from the perturbative density functional theory of Ramakrishnan and Yussouff ref78 (a) via expanding the two-particle direct correlation function, $c_{2}(k)$ up to the $4th$ order term, a procedure that yields a single peak in the direct correlation function, as detailed in Refs. ref24 (b), ref25 ; ref174 . The dimensionless free energy measured relative to a reference liquid can be expressed in terms of reduced temperature $\epsilon(>0)$ and reduced density $\psi$, leading to a Swift-Hohenberg/Brazowskii type expression ref175 :

[TABLE]

where all quantities (including r and the differential operators) are dimensionless. In 2D, the reduced temperature – reduced density ($\epsilon$ vs. $\psi$) phase diagram of the PFC model includes stability regions for the homogeneous liquid, for a crystalline phase of triangular structure, and for a striped phase (which is unphysical if modeling at the atomic level, and must be avoided during simulation), and appropriate coexistence regions in between {see Fig. 35(a)} ref24 ; ref25 . In contrast, in 3D, the phase diagram of this model contains stability domains for the homogeneous liquid, for the bcc, hcp, and fcc crystalline phases, and for the triangular rod and lamellar structures {the 3D analogues of the 2D triangular and striped phases; see Fig. 35(b), and similarly unphysical at the atomic scale} ref176 ; ref177 . The liquid phase becomes unstable against density fluctuations beyond a linear stability limit given by the relationship $\psi_{s}=-(\epsilon/3)^{1/2}$. The liquidus and solidus lines emerge from a critical point at $\epsilon_{0}=0$ and $\psi_{0}=0$, while the interfacial free energies and interface thicknesses scale with the respective mean field critical exponents ref177 . Formation of an amorphous solid via first order phase transition has also been reported ref178 . Note that the model contains a single model parameter that incorporates a combination of the expansion coefficients of the direct correlation function, which in turn can be related to the compressibility of the reference liquid, the bulk modulus of the crystal, and the interatomic distance ref24 ; ref25 .

Although there is another (equivalent) formulation of the single-mode PFC model ref24 (b), ref178 ; ref179 , it can be readily transformed into the Swift-Hohenberg/Brazowskii form ref25 ; ref177 .

Absolute density and temperature scales can be obtained for the single-mode PFC model via fitting its phase diagram to that of real matter in a limited temperature/density window as done by Elder et al. ref24 (a); whereas another route was proposed by Wu and Karma to fix the model parameters to data taken from molecular dynamics simulations ref180 .

IV.1.2 The two-mode PFC model

A two-mode extension of the PFC model was put forward by Wu et al. ref24 (c), which was designed to promote fcc solidification in the neighborhood of the critical point, where linear elasticity is valid. With a specific choice of model parameters this model reduces to the single-mode PFC model ref24 (c).

Re-casting the respective expression for the free energy in terms of $\lambda=R_{1}/(1+R_{1})\in[0,1]$ (where $R_{1}$ is the relative strength of the first- and second-mode contributions ref24 (c)), we obtain a form compatible with Eq. 30. This parameter can be used to interpolate between the two-mode ($\lambda=0$) and single-mode ($\lambda=1$) PFC formulations:

[TABLE]

Here $\psi\propto(\rho-\rho_{L}^{ref})/\rho_{L}^{ref}$ is the scaled density difference relative to the reference liquid of particle density $\rho_{L}^{ref}$. Again, the reduced temperature $\epsilon$ can be related to the bulk moduli of the fluid and the crystal, whereas $Q_{1}=q_{1}/q_{0}$ $(=2/3^{1/2}$ for fcc ref24 (c)) is the ratio of the wave numbers of the two modes.

For $\lambda=0$, the phase diagram of the three dimensional two-mode PFC model has the fcc structure down to the critical point [see Fig. 36(a)], whereas for $0<\lambda\leq 1$, a bcc domain extending with decreasing $\lambda$ appears [see Fig. 36(b)], possibly with an hcp domain between fcc and bcc domains, as in the single-mode PFC case.

IV.1.3 The multi-mode PFC model

An $N$-mode version of the PFC approach was proposed by Mkhonta et al. ref181 , assuming that the two-particle direct correlation has $N$ peaks and it is isotropic. Then the respective free energy was written as

[TABLE]

where $\Lambda$, $b_{i}$, $\tau_{x}$, and $\xi$ are phenomenological constants, whereas the $N$ peaks of the direct correlation function are located at wave numbers $Q_{i}$ $(i=0,1,...,N-1)$ emerging from a Landau-Brazowskii expansion of the direct correlation function ref24 (b), ref25 . To address the competition of bcc and fcc phases Tang et al. ref182 have chosen $N=3$, $Q_{0}=1$, $Q_{1}=\sqrt{2}$, and $Q_{2}=\sqrt{3/2}$. With appropriate values of the other parameters, a phase diagram, which resembles that of the single-mode PFC model was obtained ref182 .

IV.1.4 Other advanced PFC models

Another PFC model used in investigating nucleation is the eighth order fitting version, EOF-PFC, in which eighth order expansion of the two-particle direct correlation function around its maximum $(k=k_{m})$ was used to approximate the properties of bcc materials ref183 :

[TABLE]

Here the model parameters were fixed so as to recover the position, height, and the second derivative of $c_{2}(k)$ at the first peak, which was assured by the choices

[TABLE]

Using input from MD simulations for Fe as in Ref. ref180 , a fair agreement was achieved between the theoretical predictions and MD results for the volume change on melting, the bulk moduli of the liquid and solid phases, and the magnitude and anisotropy of the crystal-liquid interfacial free energy ref183 .

IV.2 Application of the PFC models to crystal nucleation

Once the free energy functional is defined, nucleation can be addressed in ways analogous to those seen in the case of coarse grained PF models: (i) one may compute the properties of nuclei via either solving the respective Euler-Lagrange equation, or finding the nucleation barrier on the free energy surface by methods like the elastic band approach, or (ii) one may perform dynamic simulations of nucleation, while relying on the equation of motion. As pointed out earlier, the two approaches lead to comparable results provided that dynamic effects do not overwhelm the thermodynamic preferences.

In this subsection, we review formulations of the Euler-Lagrange equation and the governing equations for diffusive and hydrodynamic modes of density relaxation, mentioning briefly the numerical methods that were used to solve these equations. Finally, we review the molecular scale results these models predict for homogeneous and heterogeneous crystal nucleation, and for the formation of new crystal grains at the solidification front (GFN).

In one approach to modeling a crystalline substrate, a term $V(\mathbf{r})\psi$ can be incorporated into the free energy density, which acts in the area/volume occupied by the substrate ref177 ; ref184 . The potential can be written in the form $V(\mathbf{r})=[V_{s,0}-V_{s,1}S(a_{s},\mathbf{r})]f(\mathbf{r})$, where $V_{s,0}$ governs the adsorption of crystal layers, and $V_{s,1}$ is the amplitude of the periodic part of the potential. The crystal structure of the substrate is set by the function $S(a_{s},\mathbf{r})$ ref177 , while $a_{s}$ is its lattice constant. The shape of the substrate is defined by an envelope function $f(\mathbf{r})\in[0,1]$.

IV.2.1 Euler-Lagrange equation and other methods to find free energy extrema

The Euler-Lagrange equation (ELE) for the single field PFC models can be formally written as

[TABLE]

where $\psi_{0}$ is the reduced particle density of the liquid, relative to which the free energy is to be calculated ref177 . For example, when determining the properties of the nuclei in an undercooled/supersaturated bulk liquid, it is the reduced particle density of the unperturbed liquid.

Methods based on solving the ELE have been developed to determine the equilibrium properties of the fcc-liquid ref185 , bcc-liquid ref186 , and the fcc-amorphous solid ref185 interfaces, including the anisotropy of the bcc-liquid ref186 interfacial free energies. It has been reported that the interface thickness diverges, whereas the interfacial free energy tends to 0 with the respective mean field exponents ref177 .

The solution of the ELE was used to map the free energy surface in 2D and 3D in the vicinity of the nucleation barrier. Unlike the case of coarse grained PF models, where the ELE has usually one solution (the unstable nucleus, represented by a saddle point in the function space), here multiple minima exist. Since in the PFC model the crystal structure of the molecular clusters evolves automatically, the free energy surface is rough on the molecular scale with many local minima corresponding to the individual cluster configurations, whose lower envelope outlines the shape of the nucleation barrier. Illustrative results are shown in Fig. 37 for homogeneous and heterogeneous nucleation of the triangular phase in 2D at a reduced temperature of $\epsilon=0.5$, where faceted crystal morphology develops ref177 .

The same approach has been used to evaluate the nucleation barrier for homogeneous nucleation ref177 ; ref185 , heterogeneous nucleation ref177 ; ref184 , and particle induced freezing ref177 ; ref184 within the single-mode PFC model both in 2D and 3D. It has been reported for the heterogeneous process that the nucleation barrier and the contact angle are non-monotonic functions of the lattice mismatch (Fig. 38) ref184 .

A simplified string method has been used to find crystal nuclei in the single-mode PFC model for homogeneous ref187 and heterogeneous ref188 nuclei forming on a structureless substrate at small $\epsilon$ values, where the interface appears to be anisotropic only for small clusters ref187 . The structure of the respective homogeneous crystal nuclei is shown in Fig. 39.

In the case of a structureless substrate the nucleation barrier for heterogeneous nucleation decreases monotonically, whereas the contact angle shows a non-monotonic variation with increasing interaction strength of the potential by which the substrate interacts with the solidifying matter (Fig. 40).

IV.2.2 PFC with diffusive dynamics (DPFC)

In the original PFC model an overdamped conservative dynamics was assumed, for which the equation of motion (EOM) has the dimensionless form

[TABLE]

where the fluctuations are represented by a colored Gaussian noise $\zeta$ of correlator $\langle\zeta(\mathbf{r},\tau)\zeta(\mathbf{r}^{\prime},\tau^{\prime})\rangle=-\alpha\nabla^{2}g(|\mathbf{r}-\mathbf{r}^{\prime}|,\sigma)\delta(\tau-\tau^{\prime})$. Here $\alpha$ is the noise strength, and $g(|\mathbf{r}-\mathbf{r}^{\prime}|,\sigma)$ a high frequency cutoff function ref153 ; ref184 ; ref185 for wavelengths shorter than the interatomic spacing ($\sigma=\pi\sqrt{6}$). This EOM can be deduced from the dynamical density functional theory of colloidal suspensions after making a number of approximations ref174 . The DPFC model is, therefore, expected to be applicable to colloidal crystal aggregation. In such systems colloidal particles float in a carrier fluid, and when the crystal grows into the liquid of lower density, a depletion zone forms ahead of the front, into which new particles can move in from the bulk liquid only via Brownian motion, leading to a diffusion controlled growth process, for which the front velocity scales as $v\propto t^{-1/2}$ ref190 . Beyond a critical undercooling or supersaturation diffusionless growth takes place, displaying steady state growth ref190 .

2.1 Homogeneous and heterogeneous nucleation

Besides describing crystallization kinetics, the long-time solutions of the EOM can also be used to find the equilibrium configurations (a method essentially equivalent to using relaxation methods in solving the ELE). This route was followed to determine the interfacial free energy for the fcc-liquid interface in the single-mode PFC model ref191 . Remarkably, in the single-mode PFC model the fcc-liquid and bcc-liquid interfacial free energies fall quite close to each other (Fig. 41) ref186 . The same model was employed to address homogeneous and heterogeneous crystal nucleation in 2D and 3D. Finally, a similar DPFC type approach was used to find the free energy for crystal nuclei in a binary version of the single-mode PFC model ref191 , where the size dependent interfacial free energy could be fitted well (see Fig. 42) by an expression

[TABLE]

derived ref192 from a phenomenological diffuse interface theory (DIT) ref128 . Here $r$ is the radius of the nucleus (considering the surface of tension), $\delta$ is the characteristic interface thickness, and $a$ a numerical constant.

2.2 Two-step homogeneous nucleation

As mentioned earlier, nucleation via an amorphous precursor appears to be a common mechanism of crystal nucleation in various systems including colloids. Since the EOM of the DPFC model can be viewed as an approximation of the EOM of the relevant dynamical density functional theory, it makes sense to expect that the DPFC model can recover this feature of the colloids. The first DPFC study that reported a competing nucleation of bcc and amorphous solids is by Berry et al. ref178 . It yielded comparable critical sizes for the two phases and a $10\%$ to $25\%$ smaller nucleation barrier for the amorphous structure (Fig. 43). Instantaneous and continuous quenching simulations by Tóth et al. ref177 ; ref185 have shown a two-step nucleation process: first, the nucleation of globules of amorphous structure formed by a disordered arrangement of the number density peaks, followed by heterogeneous nucleation of the bcc structure on the amorphous globules (Fig. 44), whose structure closely resembles that of the bulk amorphous structure reported in continuous cooling simulations ref178 .

Structural changes during the two-step nucleation process were analyzed in the single-mode DPFC model in terms of the bond-order parameters of Steinhardt et al. ref60 (a) and Lechner and Dellago ref60 (b). It was shown ref193 that on the $\bar{q_{4}}$ vs. $\bar{q_{6}}$ map, the domain of the amorphous precursor forming in the DPFC model essentially coincides with that of the liquid phase observed in molecular dynamics simulations for the Lennard-Jones system (Fig. 45). During the amorphous-precursor-assisted bcc nucleation process, parallel occurrence of neighborhoods displaying both amorphous, MRCO (as defined by Tan at al. ref62 ), and bcc-like structures were observed, but with strongly time dependent relative amounts (see Figs. 46 and 47) ref193 . An analysis within the framework of the two-mode PFC indicates that nucleation of the amorphous globules may precede the nucleation of the crystalline phase (fcc in this study), as the amorphous solid-liquid interfacial free energy is about 2/3 of the crystal-liquid interfacial free energy ref185 .

The EOF-PFC model provides a refined thermodynamic description for bcc systems ref183 . With parameters taken from MD simulations for Fe ref183 , the pair correlation function $g(r)$ for the amorphous nucleation precursor predicted by this model ref177 is in a remarkable agreement with the $g(r)$ from MD simulations ref194 for amorphous Fe (see Fig. 48). It is also interesting that the splitting of the second peak of $g(r)$ differs significantly for the single- and two-mode PFC models ref185 .

The possibility of two-step nucleation in 2D was investigated first in the framework of the single-mode DPFC model, in a study that reported the presence of amorphous nucleation precursor, while excluding the presence of the hexatic structure ref195 . In a recent study, two-mode DPFC simulations were used to investigate the role of nucleation precursors in 2D for competing square and triangular crystals ref196 . On the basis of structural analysis in terms of bond order parameters, the authors suggest that the bond-order fluctuations trigger the formation of intermediate precursors, while the packing density of these precursors determine the structural transformation pathways from the intermediate phases to the crystal ref196 .

2.3 Competing fcc and bcc nucleation

While seeds of bcc and fcc structures tend to grow in the single- and two-mode DPFC models, owing to an aggressive nucleation of the amorphous phase, competing nucleation of the bcc and fcc phases was not observed in the regimes accessible for dynamic simulations ref185 . Study of the latter phenomenon, however, turned out to be tractable ref182 in the three-mode DPFC model ref181 that relies on three sets of reciprocal lattice vectors, when expanding the free energy of the crystal. Starting from the undercooled liquid, clusters of short-range order (SRO) and medium range order (MRO) and finally crystalline long range order (LRO) evolve. (The authors defined these regimes in terms of the magnitude of the bond-order parameter $\bar{q_{6}}$ as follows. SRO: $\bar{q_{6}}<0.28$, MRO: $0.40>\bar{q_{6}}>0.28$, and LRO: $\bar{q_{6}}>0.4$). Crystal nucleation begins with the formation of MRO clusters structurally similar to the crystal that will eventually nucleate. The bcc and fcc nuclei form from MRO clusters in several steps: (i) a thin platelet of MRO appears first; (ii) it grows to form three dimensional MRO clusters; (iii) crystal embryos form in these clusters; finally (iv) crystal embryos transform into nuclei of the stable crystalline phase. The authors also identified the mechanisms by which bcc precursors assist in the formation of fcc nuclei: it was found that (112) surfaces and steps on (110) surfaces of the bcc precursors are energetically favorable for initiating the nucleation of the fcc phase, a phenomenon reflected in the specific orientational relationships observed between the fcc and bcc structures (Pitsch, Nishiyama-Wassermann ref197 and Kurdjumov-Sachs ref198 relationships). With increasing reduced temperature the following structural sequences were observed during nucleation (Fig. 49): close to the critical point $(\epsilon=0.05)$ the bcc phase forms from the MRO domains, at $(\epsilon=0.15)$ formation of the fcc phase from bcc precursor is reported, whereas at $(\epsilon=0.4)$ solid amorphous domains consisting of SRO and MRO precede the formation of the fcc structure from the MRO regions. Direct transition from SRO to crystalline structure has not been observed.

2.4 Crystallization kinetics

The binary version of the single-mode DPFC model ref24 (b) was used to model crystallization kinetics in an undercooled alloy in 2D ref19 (c). The dependence of transformation kinetics on the number of initial seeds has been investigated (Fig. 50). It was found that, probably because of the complex geometry, freezing does not follow a JMAK-type transformation kinetics, the Avrami-plot is not a straight line, and the Avrami-Kolmogorov exponent depends strongly on the transformed fraction with a functional form that differs for the seed numbers $N_{s}=$ 5, 50, and 500. The origin of the discrepancy from the theoretically expected exponent $p=d/2=1$ for fixed number of seeds and diffusion controlled growth in two dimensions ref40 (c) is not clear. We can speculate that this phenomenon can be associated with the complex growth geometry: early stage dendritic growth forms with only a few side arms (that are far less developed than in Fig. 26), which does not evolve self-similarly, were observed for $N_{s}=$ 5 and 50, whereas more rounded crystal grains developed in the case of $N_{s}=$ 500. Other sources of deviations from theoretical expectations include the small sample size and that a fraction of the liquid phase could not crystallize on the accessible time scale. Further work is needed to clarify the origin of this complex behavior.

2.5 Free growth limited mechanism of particle-induced freezing

The single-mode DPFC model was used to address particle-induced freezing (athermal nucleation) in undercooled liquids, representing the foreign particles/substrate by periodic potentials confined to their respective part of the simulation ref184 . Solving the ELE it was shown that for low reduced temperatures (small anisotropy of the interfacial free energy) the growth limited mechanism proposed by Greer and coworkers works quite well ref184 . In turn, for high anisotropies of the interfacial free energy (that occur at large reduced temperatures), free growth takes place at a critical supercooling/supersaturation, where the critical thickness of the adsorbed crystalline layer cannot be related to the size/shape of the homogeneous crystal nucleus. It was also shown that the shape and amount of the crystalline matter adsorbed on the substrate just before the onset of free growth depends strongly, and in a non-monotonic way, on the misfit between the lattice constants of the foreign particle and the adsorbed crystal ref184 . The shape of the adsorbed crystal volume just before free growth depends strongly on the size of the foreign nanoparticle ref199 (Fig. 51): in the case of cube shaped foreign particles, with increasing linear size $L$ of the cube, the adsorbed crystalline matter forms pyramid shapes, then spherical caps on the faces, finally converging towards a four-cornered curved surface that emerges from macroscopic simulations ref42 (b).

2.6 Nucleation on patterned surfaces

The DPFC strategies developed to handle surfaces interacting with crystallizing liquids allow the modeling of crystal nucleation on surfaces modulated on the nanoscale. Examples are nucleation in 2D rectangular corners, where owing to the mismatch of the triangular lattice with the rectangular corner, the latter is not a preferred nucleation site ref172 . This approach allows the investigation of nucleation on substrates of more complex geometries (see e.g. Fig. 52) ref195 ; ref199 . Another way to model the effect of a substrate is to use a periodic potential to represent the substrate layer underlying the first layer of the crystal. This approach has been used widely ref200 to model various pattern formation problems ref201 (involving nucleation and growth) on specific substrates (Fig. 53) ref195 ; ref200 .

Evidently, lattice mismatch plays an important role in phase-selection in heterogeneous nucleation on patterned surfaces: appropriate choice of the lattice constant of the substrate initiates heteroepitaxial growth of the fcc and bcc structures, and the formation of an amorphous structure at large mismatch (Fig. 54) ref184 .

IV.2.3 Hydrodynamic PFC model of freezing (HPFC)

In contrast to colloidal suspensions, in simple liquids there is no carrier fluid, and the liquid flows to low density domains, rapidly homogenizing the density. Accordingly, in this case crystal growth occurs with a constant velocity proportional to the inverse of viscosity $v\propto\mu^{-1}$. To achieve this, hydrodynamic modes of density relaxation need to be considered in the dynamic equations. A hydrodynamic PFC (HPFC) model, proposed recently, realizes this ref202 . It relies on the momentum transport and continuity equations of fluctuating nonlinear hydrodynamics ref203 valid down to the nanometer scale ref204 :

[TABLE]

[TABLE]

where $\mathbf{p}(\mathbf{r},t)$ is the momentum field, $\rho(\mathbf{r},t)$ the mass density, $\mathbf{v}=\mathbf{p}/\rho$ the velocity, $\nabla\cdot\mathbf{R}=-\rho\nabla\frac{\delta F[\rho]}{\delta\rho}\approx-\rho_{0}\nabla\frac{\delta F[\rho]}{\delta\rho}$ the divergence of the reversible stress tensor, $\frac{\delta F[\rho]}{\delta\rho}$ the functional derivative of the free energy with respect to density. $\rho_{0}$ stands for the reference density, and $\mathbf{D}=\mu_{S}\{(\nabla\otimes\mathbf{v})+(\nabla\otimes\mathbf{v})^{T}\}+[\mu_{B}-\frac{2}{3}\mu_{S}]\mathbf{I}(\nabla\cdot\mathbf{v})$ the dissipative stress tensor, while $\mu_{S}$ and $\mu_{B}$ are the shear and bulk viscosities. A momentum noise satisfying the fluctuation-dissipation theorem is added to Eq. 38, which is characterized by the covariance tensor:

[TABLE]

To avoid interatomic flow due to the extreme density gradients, when solving the Navier-Stokes equation in the crystalline phase, the momentum and density fields were coarse-grained, and these coarse grained fields were used in computing the velocity field, $v=\hat{\mathbf{p}}/\hat{\mathbf{v}}$, both in the convection and dissipation terms ref200 . It has been shown that with model parameters corresponding to a hypothetical 2D Fe the model recovers steady state growth $v\propto\mu^{-1}$, a linear dispersion relation for the long wave acoustic phonons, and a reasonable stress relaxation ref202 .

This model was used to address homogeneous, heterogeneous, particle induced, and growth front nucleation using model with parameters from Refs. ref202 ; ref205 , corresponding to a hypothetical 2D Fe ref202 .

Other PFC-based models of freezing were put forward very recently, which rely on either simplified or full hydrodynamics ref206 ; ref207 , on a modified Chapman-Enskog procedure to derive the hydrodynamic equations ref208 , or on amplitude expansion when coupling the PFC model to the Navier-Stokes equation with additional terms in the continuity equation ref209 . However, to our knowledge, none of these models have been used to address crystal nucleation so far.

We stress that the HPFC models differ from the DPFC model in only the dynamic equations; accordingly, solutions for the ELE are the same whether the HPFC or DPFC model is considered. Yet, the time scales can differ considerably, as well as the outcome of dynamic simulations.

3.1 Two-step homogeneous nucleation

The open question, whether in simple liquids amorphous precursors may assist crystal nucleation, is of practical importance. Although a study of Lutsko and Nicolis ref63 indicates that even in simple systems as the Lennard-Jones fluid, formation of a dense liquid precursor is expected to precede crystal nucleation, it is unclear whether it is so in the case of metal liquids. Although metals are expected to behave similarly to the HS and LJ systems, a direct investigation of the problem is yet unavailable.

3.2 Crystallization kinetics

In recent studies ref210 ; ref211 , homogeneous nucleation and growth was investigated in a deeply undercooled metastable liquid that was close to its linear stability limit (the simulation was performed at $\epsilon$ = 0.1158 and $\psi_{0}=-0.1982$, on a $2048^{2}$ rectangular grid). The results are displayed in Fig. 55. Nucleation started after a dimensionless incubation time of $t_{0}=3490$. After further growth the crystallites impinged upon each other, forming grain boundaries (coherent dislocation walls; see 55(a)-(f)). Remarkably, even at this highly undercooled state, crystallization follows the Johnson-Mehl-Avrami-Kolmogorov kinetics [see Eq. (5)], with a time constant of $\tau=1763$, and an Avrami-Kolmogorov exponent of $p=3.31\pm 0.03$, obtained from the slope of the Avarami plot, $\ln[-\ln(1-Y)]$ vs. $\ln(t-t_{0})$; a value which indicates a slightly increasing nucleation rate combined with two-dimensional steady state growth ref40 (c).

3.3 Free growth limited particle induced freezing

Crystallization of undercooled melts is normally initiated by foreign particles floating in the liquid volume. In Sections III.A.1.2 and IV.B.2.5, we reviewed simulation results for this process obtained in the framework of a coarse grained PF theory and the DPFC model (colloidal system). These models, in agreement with the macroscopic description ref42 , predict that when the foreign particles of a given size are wet by the crystalline phase, there occurs a critical undercooling, beyond which the particles initiate free growth of the crystal. Apparently, this behavior remains valid also for simple liquids as predicted by the HPFC model ref210 . Analogously to the approach used in the DPFC model, a term $V(\mathbf{r})\rho(\mathbf{r},t)$ was introduced to the free energy density to model the foreign particle. Here $V(\mathbf{r})=0$ outside of the volume of the foreign particle, whereas a periodic potential, which determines the crystal lattice in the foreign particle was prescribed inside. A square lattice of the same lattice constant $a$ as the growing triangular crystal was used, which yielded a nearly ideal wetting at the foreign particle-crystal interface. The linear dimension of the square particle was $L=32a$. Varying the undercooling by multiplying the model parameter $C_{0}$ by factors in the range $0.990<\xi<0.99$, it was found that the same behavior was observed, indicating that the free growth limited mechanism applies generally for all these systems. The results are summarized in Fig. 56, which displays the long-time solutions of the dynamic equation as a function of $\xi$, and shows a transition from stable configurations to free growth.

3.4 Growth front nucleation near the linear stability

The formation of new grains at the growth front has been studied using the HPFC model in the neighborhood of the stability limit of the liquid phase (Fig. 57) ref193 . The simulation was performed in the presence of noise in a regime, where the liquid is metastable with respect to the fluctuations. A crystal seed was used to initiate crystal growth. The reduced temperature was $\epsilon=0.1158$, and the initial density $\psi_{0}=-0.1982$, while the reduced temperature at the stability limit was $\epsilon_{c}=0.1178$. The HPFC simulation was performed on a $2048^{2}$ rectangular grid, using periodic boundary conditions on all sides. A complex order parameter $g_{6}=\sum_{j}\\ exp\{6i\theta_{j}\}$, was used to monitor the local structure, where the summation was for the nearest neighbors, while $\theta_{j}$ is the angle of the $j^{th}$ neighbor in the laboratory frame. $|g_{6}|$ then represents the degree of order, whereas the phase of $g_{6}$ is the local crystallographic orientation. A Voronoi polyhedral analysis has also been performed for the growing crystal employing the same coloring scheme as in Fig. 56. As the crystallite grows, gradually new orientations appear via two mechanisms of GFN. (a) Dislocations enter first near the corners and the center of the edges of the crystallite, and later at cusps of the front. These were interpreted as dislocations, emerging presumably due to the stress emerging from the non-equilibrium lattice constant forming at the growth front ref212 ; ref213 , however, further analysis of this process is needed. (b) Small crystallites nucleate in the vicinity of the solidification front, which was associated with the interference of the density waves emanating from a rough solid-liquid interface. The latter mechanism resembles “satellite” crystal formation, as reported in large scale MD simulations ref167 ; ref214 . A more detailed analysis is available in Ref. ref210 . These modes of growth front nucleation lead to the formation of nanoscale spherulite-like polycrystalline morphologies ref210 , raising the possibility that similar processes may play a role in spherulitic solidification on a larger scale.

IV.2.4 Numerical methods

Equilibrium methods: A specific mathematical difficulty associated with the PFC models and their derivatives is that the free energy density contains high order spatial derivatives. Such derivatives occur both in the Euler-Lagrange and the dynamic equations.

For example, the properties of nuclei in PFC type models can be determined by solving either the Euler-Lagrange equation or by using other saddle point finding methods. A pseudo-spectral successive approximation scheme combined with the operator-splitting technique proved used for solving ELE in a number of studies ref172 ; ref177 ; ref184 ; ref185 ; ref186 , whereas the string/elastic band methods ref146 ; ref147 ; ref213 were employed to find the nucleation barrier and the minimum free energy path for homogeneous and heterogeneous nucleation ref187 ; ref188 . Other possible methods are reviewed in Ref. ref146 . A specific feature of the nucleation barrier in the PFC models, especially at large reduced temperatures, is that due to the atomic structure of the cluster, the free energy surface at the nucleation barrier is rough with many local minima corresponding to different cluster configurations, the lower envelope of which minima outlines the barrier itself both in 2D and 3D ref172 ; ref177 .

In solving the dynamic equations for the DPFC and HPFC models various numerical methods were used, including the finite difference ref24 (b), finite element ref216 and spectral approaches combined with various operator splitting strategies ref217 . A numerical method that proved fairly successful in solving the equation of motion of the DPFC model is a spectral semi-implicit scheme that relies on parallel fast Fourier transform, while assuming periodic boundary conditions at the perimeter ref217 .

Numerical problems associated with combining PFC type models with hydrodynamics were addressed in Refs. ref207 ; ref218 .

IV.3 Discussion

Summarizing, the PFC type models are able to address on the molecular scale a broad range of phenomena associated with homogeneous and heterogeneous crystal nucleation ref172 such as

• •

the role of precursor structures ref177 ; ref185 ; ref195

• •

phase selection between competing phases ref182

• •

mismatch vs. wetting and nucleation barrier ref184

• •

effect of patterned and non-patterned surfaces ref184 ; ref199

• •

particle induced freezing ref184 ; ref210

• •

growth front nucleation ref193 ; ref210

• •

formation of satellite crystals ref193 ; ref210

• •

surface alloying ref200 (a),(b)

• •

crystallization kinetics ref19 (c), ref210

• •

hydrodynamic effects on solidification ref193 ; ref199 ; ref210 .

Apparently, the family of PFC models offers a flexible framework to address microscopic details of crystallization in a qualitative way. Although a large variety of physical phenomena are automatically incorporated into such models (elasticity, plasticity, grain boundary dynamics, atomic scale structure, anisotropies, etc.), the adaptation of early PFC models to real materials does not seem to be straightforward. A possible solution is offered by the XPFC model ref24 (d), ref219 . Introducing appropriate peaks in the two-particle correlation function, this model is able to capture a variety of solid structures including crystalline (fcc, bcc, hcp, sc in 3D, and hexagonal, kagome, honeycomb, square, rectangular, etc. lattices in 2D ref220 ), the quasi-crystalline, and the amorphous phases, in multiphase/multicomponent systems. In a few cases, quantitative modeling was made possible by specific PFC approaches ref200 ${(a),(b)}$ ref221 . Revisiting previous work done using simpler PFC models in the framework of the XPFC theory raises the possibility of addressing crystal nucleation in real materials with improved accuracy.

A recent step in this direction is a study of two-step crystal nucleation in a phase separating liquid that has been addressed within a combination of a binary version of XPFC model with diffusive dynamics ref222 . Precipitation from the solution happens here via crystal nucleation in the solute-rich droplets, which formed during liquid phase separation (Fig. 58). The predictions accord qualitatively with experiments performed on various systems.

V Nucleation prefactor from phase-field/density functional theories

The previous two sections mostly dealt with the magnitude of the nucleation barrier and the phenomena influencing it. To compute the nucleation rate, however, one needs to evaluate the nucleation prefactor $J_{0}$ as well (see Eqs. (3) and (4)). As mentioned in Section II, in the case of crystal nucleation, the nucleation prefactor from the classical kinetic approach appears to be acceptable for orders of magnitude estimates: it is $\sim 2$ orders of magnitude too low, when compared to the result from MD simulations. The thermal transport limited nucleation prefactor derived by Grant and Gunton ref223 is probably less relevant here, as the rate limiting process in crystal nucleation is usually molecular mobility (self-diffusion/chemical diffusion) with the exception of perhaps pure metals, where barrierless crystal growth was observed ref224 . Besides evaluating the nucleation barrier for a GL type phase-field model with a short-range nonlocal interaction, Roy et al. determined the nucleation rate with the exception of a dynamical prefactor ref225 .

In a recent study, the single-mode DPFC model was used to study the kinetics of crystal nucleation in the undercooled state ref226 . Various aspects of crystal nucleation were explored, including transient nucleation, and the steady state nucleation rate was also evaluated as a function of the reduced temperature. The results are in qualitative agreement with the scaling behavior predicted by the CNT, however, it is suggested that a multivariable approach ref226 would be necessary to develop a quantitatively correct nucleation theory.

Lutsko and coworkers ref22 (a),(b), ref227 put forward a nucleation theory for colloidal systems, which is based on nonlinear fluctuating hydrodynamics. They pointed out that the single molecule attachment/detachment assumption of the classical kinetic approach is invalid at large supersaturations, where cluster-cluster reactions have relatively higher probability. Such problems can be automatically treated within DPFC and HPFC simulations in which thermal fluctuations are represented by noise added to the dynamic equations. However, further work is needed to develop improved analytical expressions for the nucleation prefactor accurate in simple and in highly undercooled liquids.

VI Summary and future directions

We have reviewed the advances made by coarse grained and molecular scale phase-field models of crystal nucleation during the past decades.

It is evident that coarse-grained PF models based on the Ginzburg-Landau free energy are of limited use at present, although this is the only PF approach that works fairly well for the hard sphere system, when comparison is made between the nucleation barriers predicted without adjustable parameters and the nucleation barrier evaluated from atomistic simulations using umbrella sampling. However, in the case of bcc and fcc metals, the parameter-free GL predictions are far less successful. A possible reason for this failure might be that the trivial temperature dependence that applies for the hard sphere system is not valid for real liquids. Owing to the simplifications made during their derivation the PF models can be viewed as advanced phenomenological models, whose input parameters can be evaluated from microscopic modeling such as MD ref228 . Improvements of the coarse grained models might be expected from e.g., the amplitude expansion of molecular scale theories.

While the simple dynamical density functional theories, such as the DPFC and HPFC approaches, are able to recover a broad range of interesting nucleation phenomena (precursors in two- and multistep nucleation, competing phases, interaction with foreign walls, heteroepitaxy, flow effects, dispersion of the acoustic phonons, etc.) and can address nucleation on the molecular scale, their application to real materials remains limited. Recent activities aimed at making the PFC models quantitative for single- and multicomponent systems ref180 ; ref183 ; ref219 ; ref220 ; ref221 ; ref222 ; ref229 open up the road to addressing nucleation under close to realistic conditions. An attractive feature of these models is that being molecular scale approaches they connect structure and the physical properties. Furthermore, the PFC–type models may serve as a starting point for developing physically reasonably well established coarse grained nucleation theories that are able to consider the effect of different crystal structures and the related anisotropies. Despite promising results achieved in this area in the recent years, a few essential problems remain open. For example, it is yet unclear how much of the knowledge learned from experiments on colloidal systems remains transferable to simple liquids ref230 . For example, it is not known whether amorphous precursors indeed assist crystal nucleation in simpler liquids, such as molten metals. Another interesting question is, how far these findings might apply for biological systems. Owing to the obvious experimental difficulties associated with following the trajectories of individual atoms in a crystallizing monatomic liquid or in liquids of small molecules, atomic scale modeling relying on the MD, MC, HPFC, and other classical dynamical density functional techniques are expected to play a decisive role in clarifying these issues.

A further interesting development is the combination of the atomic scale studies (MD, PFC) with coarse grained phase-field modeling to form a multi-scale approach that bridges the atomic and macroscopic scales ref231 .

Finally, a practical issue: One often needs a quick “educated guess” for the nucleation rate. In the case of homogeneous nucleation, this can be best done on the basis of the classical expression $J_{SS}=J_{0}\exp\{-W^{*}/kT\}$. However, one needs to use the nucleation prefactor from the CNT multiplied by a factor of $\sim 100$ ref4 , while taking the nucleation barrier $W^{*}$ from

(i) the droplet model of the CNT with $\gamma=\gamma_{SL,eq}(T/T_{f})$, or

(ii) numerically solving the Euler-Lagrange equation of the standard PF model, or

(iii) a simple analytical diffuse interface theory (DIT) that considers the thickness of the solid-liquid interface ref128 .

Note, however, that the accuracy of such predictions hinges on the validity of a number of assumptions, such as the lack of precursors, the absence of cluster-cluster interaction, diffusion controlled dynamics, large distance from metastable critical points, and validity of the hard-sphere-like temperature scaling of the interfacial free energy.

In the heterogeneous case, one may rely on the free growth limited mechanism of particle induced freezing by Greer and coworkers ref42 , for which the size- and perhaps the shape distribution of the foreign particles are needed. Handling of other, more complex situations ref1 (a), ref232 , may require an extended knowledge on the heterogeneities.

Summarizing, since establishment of classical nucleation theory more than 90 years ago ref233 and its adaptation to crystal nucleation in undercooled liquids about a quarter of a century later ref234 , an incredible amount of information accumulated on crystal nucleation experimentally, theoretically, and via computer simulations; and a broad range of theoretical approaches were proposed that consider various aspects of the phenomenon. Yet, clearly, there are substantial opportunities for theoretical progress.

VII Disclaimer

The mention of commercial products, their source, or their use in connection with the material reported herein is not to be construed as either an actual or implied endorsement by the National Institute of Standards and Technology.

Acknowledgements.

We thank G. C. Sosso and A. Michaelides for the tabulated data on ice nucleation. This work was supported by the National Agency for Research, Development, and Innovation, Hungary (NKFIH) under Contract Nos. K-115959, KKP-126749, and NN-125832.