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

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.


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

A. 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 [26] (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π/3)R 3 ∆ω + 4πR 2 γ SL , where R is the radius of the surface on which the surface tension acts (the surface of tension), ∆ω is the thermodynamic driving force of solidification (the volumetric grand potential difference between the solid and the liquid; ∆ω < 0 for the undercooled liquid), and γ 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π/3)(γ 3 SL /∆ω 2 ) as a function of the size at R * = −2γ SL /∆ω. 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: and for n > 1Ṅ n = a + n−1 N n−1 + a − n+1 N n+1 − (a + n + a − n )N n .
In Eqs.
(1) and (2), N n is the number n-molecule clusters, while a + n = O n Γ exp{−(W n+1 − W n )/2kT } and a − n = O n−1 Γ 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 Γ = 6D/λ 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 λ 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 [1](a), [2](a), [13,14], 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 [14](b), [27].
(1) and (2) numerically, steady state nucleation occurs after a transient period of length τ ≈ Kλ 2 kT (n * ) 2/3 /(Dv m ∆ω), 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. [13](a)]: Here N eq,n = ρ 0 exp{−W n /kT } is the equilibrium population of the n-molecule clusters, and J 0 = ρ 0 O n * ΓZ is the pre-exponential factor of the nucleation rate, while ρ 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πkT )} 1/2 accounts for the decay of the critical size clusters. Fitting Eq. (3) with unknown D and γ SL to nucleation rate experiments on oxide glasses imply that the order of magnitude of the classical estimate for J 0 is reasonable [1](a), [2](a), [14](a). Molecular dynamics investigations suggest that the classical prefactor for crystal nucleation might be 2 orders of magnitude too low [4]. For vapor condensation, molecular dynamics simulations found J 0 that is ∼ 1 order of magnitude larger than the CNT prediction [28], whereas a dynamic extension of the classical density functional theory yielded a reasonable agreement with the CNT result for J 0 [29]. 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 [30], implying that the kinetic part of the CNT is reasonably accurate. A field theoretical expression, J SS = J 0 exp − W * /kT , similar to Eq. (3) has been derived by Langer [31], which lead to comparable results to CNT in the few cases in which comparison was made [32,33].
Apparently, the large (several orders of magnitude) deviation between experimental and theoretical (CNT) nucleation rates reported for oxide glasses and other substances [1](a), [2](a), [14](a), [34] 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 γ SL fails when predicting the nucleation barrier [5](a). The failure of the droplet model for small clusters was also borne out by Monte Carlo simulations for the Ising model [35]. In analogy to the case of liquid droplets, for small clusters, corrections may be introduced for the surface tension (interfacial free energy) [36]. A fairly frequently used correction is Tolman's that introduces a size dependent interfacial free energy: γ(R p ) ≈ γ ∞ /[1 + 2(δ T /R p )], where R p = [3W * /(2π∆f )] (1/3) is the radius of the surface of tension on which the surface tension acts, δ T = lim Re,Rp→∞ (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. [37].

B. 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 [38]: γ W L = γ W S + γ SL cos(Ψ), where γ W L and γ W S stand for the wall-liquid and wall-solid interfacial free energies, whereas Ψ 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 (Ψ), where f (Ψ) = 1 4 [2 + cos(Ψ)]{1 − cos(Ψ)} 2 is the catalytic potency factor. For small contact angles, f (Ψ) 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 − cos(Ψ)}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 [39] J SS,het = x a q(Ψ)J 0 exp{−W * hom f (Ψ)/kT }, where x a 1 denotes that fraction of all the molecules that is adsorbed on the surface of heterogeneities, and q(Ψ) = 1 2 {1 − cos(Ψ)}/ f (Ψ) ∈ [1/ √ 3, 1]. The classical approach has been adapted to various geometries of the wall, including spherical particles and cavities, depressions, and rough surfaces [40].
A generalized form of Turnbull's experimental test [41] 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) χ(T ) = γ SL (T )/γ SL,eq , and D(T ), λ, and ∆ω(T ) are known: plotting the logarithm of the normalized experimental (steady state) nucleation rate, log(J SS,exp /J 0 ), vs. the temperature dependent part of the argument of the exponential function, χ(T ) 3 ∆ω(T ) −2 T −1 , one should obtain a straight line that intersects the ordinate at log(x a q(Ψ)) ≈ log(x a ) ≤ 0 with a slope of −16πγ 3 SL,eq f (Ψ)/3k [34](c), (d). For example, the assumption of a curvature-and temperature-independent interfacial free energy (χ = 1) yields intersections many orders of magnitude too high [1](a), [2](a), [14](a), [34], indicating a failure of the original droplet model that relies on constant γ SL .

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 [42]. In this approach, the foreign particles floating in the melt are viewed as cylinders of radius R p that have ideally wetting circular faces (Ψ = 0), and non-wetting sides (Ψ = π). These idealized particles remain dormant at and below a critical undercooling, ∆T c = 2γ SL /(∆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 ∆s f is the volumetric entropy of fusion.) For undercoolings ∆T > ∆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 [43]. 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 [1](a).

C. 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 γ SL in the vicinity of the melting point [44], the associated error is usually far too high (5% to 10 %) for a conclusive test. Furthermore, the interfacial free energy may depend on temperature [45] and curvature [34](c), (d), [36], both of which could influence the nucleation rate considerably. In addition, despite a wealth of highly accurate nucleation rate data available for oxide glasses [34,46], 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 γ 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 [47] and pure metals [48], and experiments for crystallizing colloidal suspensions: [49,50] some of these approaches [47,48,50] 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: 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 * M C /kT = 19.4 vs. W * CN T /kT = 17.4) [4]. In computing the latter, the orientation average of γ SL,eq = 0.35 evaluated at the triple point was used from a MD study performed with a modified LJ potential [51]. It appears, however, that (much like in the hard sphere system) γ SL,eq ∝ T in the LJ system [45], which means that the previous value of W * CN T /kT needs to be multiplied by (T /T 3 ) 3 , yielding W * CN T /kT = 12.0 that is considerably less than the simulation result. To recover W * M C /kT = 19.4 at T = 0.6, one needs to assume a positive curvature correction from γ T =0. 6 SL,eq = 0.309 to γ nucl SL = 0.363 = 1.17γ T =0.6 SL,eq of the nucleus.

Hard-sphere (HS) system
A significant difference was reported [5](a) for W * from MC simulations and the droplet model that relies on the equilibrium value γ SL,eq /(kT /σ 2 ) = 0.617 [52,53]  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 [5](a). (Given the extreme sensitivity of W * CN T to γ SL and ∆ω, 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 [10], 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 [54], 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.

Metals
Recent excellent nucleation rate results obtained on single metal droplets by chip calorimetry [55] are reported to be in fair agreement with molecular dynamics simulations based on embedded atom potentials and the CNT [48], provided that in the latter the HS relationship γ SL ∝ T is adopted. While this might be a reasonable approximation for metals, apparently the oxide glass data fit better to a linear relationship γ = a + bT [1](a), [2](a).

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 [56] and molecular dynamics methods [47](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. [47](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 [56], a range of estimated values for the interfacial free energy [44,57], and some accurate thermodynamic data for the undercooled state, at least for I h [58]. 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 [4].
FIG. 2. Formation of the bcc crystalline phase (green particles) from clusters showing Medium Range Crystalline Ordering (red particles) in Brownian Dynamics simulation for the hard sphere system [10]. The coloring scheme is shown below, where Q6 =q6 is an average bond-order parameter of the Lechner-Dellago-type (see Ref. [59]). Reproduced with permission from Ref. [10] c 2010 National Academy of Sciences of USA.

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 [47]. The structure of the local molecular neighborhoods has been characterized in terms of bond order parameters [59]. The importance of Medium Range Crystalline Order [10,47] (MRCO, Fig. 2) and crystal-like precursor structures [4,5,60] 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 [61], the LJ [62], and HS [9] 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.

D. 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 [63], which expresses the time evolution of the transformed fraction in a d-dimensional system as 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 [40](c) Here τ = [(Ω d /p)JV d min Π d j=1 β j ] −1/p is the characteristic time of the crystallization, and p = 1 + d is the Avrami-Kolmogorov kinetic exponent, Ω d is the volume of the d-dimensional unit sphere, V min the minimum growth rate, and β j = V j /V min ≥ 1, whereas V j (j = 1, . . . , 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. [64]}. 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 [40](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 [65]. 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 [66]. The PF simulations proved useful in addressing such problems [19](a), [67], [68](c).

III. COARSE GRAINED PHASE-FIELD MODELS
The phase-field (PF) technique and its applications are described in a number of recent reviews [68]. 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 [23] and developments presented later by others [23,[69][70][71][72][73][74][75][76]. In order to characterize the local state of matter, a non-conserved structural order parameter φ(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 {φ i (r, t)} (i = 1, ..., N ). Some of the models, such as the multi-phase-field (MPF) theories by Chen et al. [70], Steinbach et al. [71], and others [72], introduce an independent phase field for every crystal grain, and work with tens of thousands of fields in describing multi-grain problems [73]. Other more economic approaches employ orientation fields to monitor the local crystallographic orientation [19](a),(c), [68](c),(f), [74][75][76].
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 Φ = {φ i }, and other slowly changing fields such as the chemical concentration field(s) 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: 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 A, B, ... denote general quadratic terms for the respective gradients: Choosing, for example, A = I (here I is the identity matrix) yields a simple sum of the SG terms, A = I − e ⊗ e [where e = (1, 1, ..., 1)] yields a pure pairwise construction, whereas A i,i = j =i φ 2 j , A ij =i = −φ i φ j realizes the anti-symmetrized (Landau-type) gradient term. Coefficients φ and 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 (Φ, 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 Θ, 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 [21,22,77], 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 [78]. 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., 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.) [79], 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 [79,80].
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) [72](g) for non-conserved and conserved fields [81], respectively: Non-conserved dynamics: Conserved dynamics: Here δF/δφ i and δF/δc i are the functional derivatives of the free energy with respect to φ i and c i , respectively. In different PF models different choices were made for the mobility matrix κ ij [68,72,73]. {For a recent review and a physically consistent formulation see Ref.
[72](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: Here λ(r) is a Lagrange multiplier emerging from the local constraint of the sum of the phase fields. Eliminating the Lagrange multiplier yields Next, we illustrate the application of the PF approach to nucleation in a few simple cases.
A. 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 (φ) 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 (ρ), or volume fraction (X).

Single-field models
A realization of such an approach is given by the following free energy functional, EOM, and ELEs: Free energy: where s(n) is an anisotropy function that depends on the components of ∇φ = {φ x , φ y } = |∇φ|n in two dimensions [82], n the surface normal, w the free energy scale, g(φ) is a double well function, p(φ) is an interpolation function varying monotonically between 0 at φ = 0 and 1 at φ = 1, ∆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 [83]. For various choices of g(φ) and p(φ) made in the literature, see Table I [ [83][84][85][86][87][88][89]. We note that the p(φ) function in the first and fifth rows of Table I were deduced [84] for non-isothermal problems following the thermodynamically consistent approach of Penrose and Fife [90].
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 (∆f = 0) the integral form of the ELE applies [16](a): Model Here GL = Ginzburg-Landau, whereas g l (φ) = λ l φ 2 ,and gs(φ) = λs(1 − φ) 2 + ∆f .
where I stands for the integrand of Eq. (12). Inserting I, one finds that 2 2 (dφ/dx) 2 = wg(φ). This relationship can be used to relate the model parameters to measurable quantities: integrating the free energy density and dφ/dx from φ = 0 to φ = 1, and (dx/dφ) from φ = 0.1 to φ = 0.9, while assuming a quartic double well, g(φ) = 1 4 φ 2 (1 − φ) 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 γ SL = 1 6 2 s 2 w/2 and d SL = 2 log(9) 2 2 s 2 /w, respectively, whereas the phase field profile φ(x) = 1 2 1 + tanh[log (9)x/d SL ] minimizes the excess free energy of the solid-liquid interface. Making 2 ∝ T and w ∝ T , one finds γ SL ∝ T and d SL = const., a behavior consistent with MD simulations [45], 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 [16](b) 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 where f = wg(φ) + p(φ)∆f , and the boundary conditions are ∂φ/∂r = 0 at r = 0, and φ = 0 for r → ∞.

Equation of motion:
As φ is a non-conserved field, an Allen-Cahn type EOM applies [91]. For the isotropic case (s = 1), one finds:φ where ζ stands for an additive noise of correlator ζ(r, t), ζ(r , t ) = 2kT M φ δ(r − r )δ(t − t ) satisfying the fluctuation dissipation theorem. In the case of p(φ) = φ 2 (3 − 2φ)), for which the tanh profile is retained for steady state front propagation, the front velocity can be obtained analytically as v P F = M φ 12∆f 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 [92], (111) and (100) interfaces: v 0,111 = (aαD/λ 2 )e −∆S/k and v 0,100 = (αa/λ) 3kT /m [93]. Here ∆F m = −v m ∆f is the molar free energy difference, β the kinetic coefficient, a the spacing of the crystal planes, α the fraction of active sites at the crystal surface (taken to be 0.27), D the self-diffusion coefficient, λ the mean free path, ∆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 Another approach, which allows one to obtain an orientation dependent M φ , is to employ kinetic coefficients evaluated from molecular dynamics simulations [68](b), [94].
for the nucleation barrier on the basis of piecewise parabolic free energy and a variable-position step function p(φ). 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 [96], which together with a piecewise parabolic free energy leads to a simple but accurate analytic formulation for vapor condensation [97,98]. A quartic double well combined with p(φ) = φ has been adopted for describing crystal nucleation in undercooled H 2 O [89], Ar, and six stoichiometric oxide glass compositions [14](c). The standard PF model [for g(φ) and p(φ), see Table  I] has been employed to address nucleation in the HS system [99], in undercooled Ar, and H 2 O [19](a) and methane hydrate [100]. Following Stroud and coworkers [101], Iwamatsu and Horii put forward Ginzburg-Landau formulations for bcc and diamond cubic structures [85]. Quantitative tests have been performed for crystal nucleation in the HS system using the SG approximation combined with Ginzburg-Landau free energies [83,102]. 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 [14](c), [19](a), [83,89,[95][96][97][98]102]. 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 [83,102].  Table 1 of Ref. [89] were used.
In general, for large clusters (R * 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 * ≈ d SL ), for which the PF result for W * is smaller than W * CN T 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. [48](a), [103][104][105][106] (The thermodynamic driving force of crystallization ∆f has been approximated using Turnbull's linear relationship [105].) 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 γ SL = γ SL,eq (T /T f ), and from a phenomenological diffuse interface theory (DIT) [34](c), (d). For comparison, data from MC simulations and experiments [48](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 ∆T /T f → 1 one observes W * /kT → 0. We note that the closeness of the PF and CNT results follows from the symmetric g(φ) and p(φ) functions of the standard PF model that yield a symmetric phase-field profile, known to result in a zero Tolman length [107], 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 [19](a).
Because of the importance of ice nucleation in undercooled water, we present here an updated comparison between available experiments [56] and two models, the standard PF model and the CNT. In the latter, two cases are considered: (a) using the γ SL (T ) = γ 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 [45], and (b) a constant value γ SL = γ 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. [4]. The other properties were taken from Ref. [89], 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 ∆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. [89], 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 ∆C p,1 = 80 J/mol/K is assumed between temperatures T 1 and where T x = 155 K, and ∆H x = 1380 J/mol are the temperature and heat of crystallization of the low density amorphous phase [58], whereas ∆C p,x = 3.6 J/mol/K is the specific heat difference between liquid and the crystal above the glass transition [58]. 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 γ SL ∝ 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 γ SL = γ 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(φ), whereas it is far less sensitive to the form of g(φ). Remarkably, for Ni, the Ginzburg-Landau prediction for the fcc structure falls far away from the MC results. In turn, the symmetric p(φ) 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.

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 [19](b), (d), where either ∇φ · n or φ can be specified, which in turn fixes the contact angle Ψ 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 [108]: where Z(φ) is Cahn's "surface function", specifying the wetting properties of the surface. Minimizing F , yields Eq. 16 and the boundary condition Here n is the outward pointing normal to surface S, whereas z(φ) ≡ ∂Z/∂φ. This boundary condition can be realized by setting either − 2 ∇φ · n = z(φ) (Model I) or φ = const. e.g. (δφ = 0) (Model II), at the boundary [19](b). In Model I, Gránásy et al.
[19](b) assumed that the presence of a flat wall does not perturb the structure of the equilibrium (planar) solid-liquid interface. z(φ) is then deduced as follows: Eq. (13) is employed at the melting point yielding , whereas the normal component of the phase-field gradient is expressed as n · ∇φ = |∇φ| cos(Ψ), where Ψ is the (contact) angle between the solid-liquid interface and the foreign surface. Combining these expressions one obtains the boundary condition [19](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 Ψ = 75 • . With increasing N and decreasing Ψ the crystallization of the liquid volume accelerates (Fig. 8). Another straightforward application is modeling of particle induced crystallization ("athermal" nucleation) [19](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 (Ψ = 45 • ) and non-wetting (Ψ = 175 • ) vertical sides, one can investigate the theoretical prediction [42] that there exists a critical undercooling ∆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)) [19](b), however, owing to the presence of fluctuations in the equation of motion, ∆T c in the simulations is about half of the theoretical prediction. Similar behavior has been reported recently in molecular dynamics simulations [111] (Fig. 9(b)).

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 [112][113][114]. 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 [112], Castro [113], and Iwamatsu [114] reported are in a general agreement with each other and Monte Carlo simulations for steady state nucleation and isotropic growth [65], although Iwamatsu reported higher values for the kinetic exponent [114] presumably because of a nucleation rate increasing with time. Different results were reported by Heo et al. [115], 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 [112,113], owing to the much smaller nucleation rate in those works, it had a negligible effect on the transformation kinetics, whereas in Ref. [115] 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 [114], whereas Castro found that the correlation length of colored noise (of Gaussian correlator) can influence the kinetic exponent significantly [113]. 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. [67], Gránásy et al.
[19](a), and Heo et al. [115], an approach that yields similar transformation kinetics as the ones based on modeling the fluctuations via adding noise to the EOM.

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 ≥ 2) free energy-order parameter relationship [95,110,116], 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 [110]. 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 [95]. 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)].

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 [69] 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 [19](a) and dendritic solidification [69].

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, φ(r, t) and the (uniform) temperature, T , but also on a conserved field that specifies the local concentration of species B, c(r, t). Following Warren and Boettinger [69], we incorporate here an SG term only for the phase field (strictly, this simplification is valid for an ideal solution [16]): where the g(φ) and p(φ) 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 [19](a). Owing to the lack of an SG term for the concentration, the respective ELE is a degenerate one (∂f /∂c = 0), which defines a one-to-one relationship between the phase-and the concentration fields, c = c(φ).
Assuming isotropic interfacial properties (spherical symmetry), the Euler-Lagrange equation for the phase field boils The quantities ∆f A and ∆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 [105] for low undercoolings. This approach has been used to address crystal nucleation in Cu-Ni, a system with nearly ideal solution behavior [19](a). Assuming homogeneous nucleation, α = 0.6 (close to 0.58 from MD simulations evaluated from the capillary wave spectrum at the crystal-liquid interface of Ni [117]), one obtains undercoolings for nucleation rates of 10 −4 to 1 drop −1 s −1 for droplets of 6 mm diameter that fall close to the experimentally observed values [118]. 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 [62].
A more complex treatment is required when addressing eutectic systems [119], 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 1 2 2 c |∇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(φ) and p(φ) 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: where it has been utilized that φ → 0 and c → c ∞ for r → ∞, and it was assumed that 2 c = 2 c (φ, c). The latter assumption was made to connect 2 c to the phase-and concentration dependent interaction coefficient of Ag-Cu in the spirit of Ref.
The radial phase-field and concentration profiles and the free energy of formation for crystal nuclei in pure Ag, Cu, and Ag 50 Cu 50 liquids are displayed in Fig. 12 together with experimental undercooling data [120]. 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. [121]). 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 [119]. In the neighborhood of the critical point an enhanced nucleation rate is observed [119], in agreement with experiment [122], density functional computations [123], and MD simulations [124].
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 [125] 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.
[126](a). A similar approach was used for solidification in monotectic systems, however, with two-phase liquid ahead of the front (Fig. 16) [126](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. [126](c).
Patterns closely resembling experimental observations were reported in a generic peritectic system [127]. 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) [127].
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.

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. angle to Ψ at the surface of substrate is as follows [19](c): where c = c(φ) 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 Ψ. Examples are available in Refs.

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 [99,102], 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. [102]), together with the height of the nucleation barrier [5]. 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 [9], 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(φ) and p(φ) 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 [5] (see Fig. 18), together with the prediction of a simple phenomenological diffuse interface theory (DIT) [34](c), [128]. In contrast, the single-and two-field PF models with the standard choice of the g(φ) and p(φ) functions envelope the predictions of the droplet model of the CNT and the self-consistent CNT (SCNT). In the latter W * SCN T = W * CN T − W 1,CN T [129], where W 1,CN T is the free energy of the monomer in the classical droplet model. Note that as pointed out by Auer and Frenkel W * CN T falls far below the MC results [5]. It is worth noting in this respect that in the MC simulations of Schilling et al. [9], in which dense amorphous clusters were observed that acted as precursors for crystal nucleation, the volume fraction of the initial liquid was FIG. 18. W * /kT as a function of the initial volume fraction of the liquid phase (φL) as predicted by various models [102]: CNT -droplet model of the classical nucleation theory; SCCT -self-consistent classical theory [129]; PFT/S1 -single order parameter PF theory with the standard double-well and interpolation functions; PFT/GL1 -single order parameter phase field theory with Ginzburg-Landau free energy; PFT/S2 -two-field (phase-and density fields) PF model with standard double-well and interpolation functions; PFT/GL2 -two-field PF model with double-well and interpolation functions from Ginzburg-Landau expansion; and DIT -a phenomenological diffuse interface theory [128]. For comparison, the nucleation barrier height from MC simulations by Auer and Frenkel [5]  about φ L = 0.54, i.e., somewhat larger than in the simulations of Auer and Frenkel [5]. 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.

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), φ 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 [19](d): where |∇φ W | is a Dirac δ−function that locates Z(φ) to surface S, while the new factor 1 − φ W locates the free energy density for the liquid and solid phases to those regions, where φ W = 0. Computing then δF/δφ yields the following EOM for φ, where it has been utilized that n = ∇φ W /|∇φ W |, an expression that is in some sense "obvious", as we added the Model I boundary condition multiplied by a δ-function to the original variation over the volume bounded by the inert wall. Thus, introducing the auxiliary field φ 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 [130].
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.

Nucleation in multi-phase-field models
Polycrystalline solidification and grain coarsening is traditionally addressed in the framework of multi-orderparameter (MOP) [70] or multi-phase-field (MPF) [71], 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 [19](a), [121] or for W * het [19](b)-(d), or take the nucleation barrier from other sources [67]}. Alternatively, one may incorporate particles into the simulations and activate them in agreement with a criterion from Greer's free growth limited model [42], as done in a number of simulations [131]. The MPF approach was also used to address nucleation on curved surfaces [132].

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 [133]. A three-phase MPF model relying on three phase-fields φ 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 ∆f (φ i , c) was supplemented with an antisymmetric differential operator term [71](a),(b).
For the fcc-bcc transition the order parameter is related to the magnitude of Bain's distortion [134] and leads to the same g(φ) and p(φ) functions as those deduced for the liquid-bcc system. The thermodynamic data were taken from a CALPHAD (CALculation of PHAse Diagrams ) type assessment [104]. Other required data were the Phase selection in the Fe-Ni alloy system as predicted by the multi-phase-field theory for heterogeneous nucleation [133]. The gray solid and dashed lines stand for the liquidus and solidus curves, whereas the symbols indicate the structure nucleated in the experiment: squares-bcc; circles-fcc [135]. The fcc-bcc phase-selection boundary predicted with a reasonable estimate of the fcc-bcc interface energy is denoted by a heavy solid line. (Reproduced with permission from Ref. [133] c 2011 American Physical Society.) interfacial energies and thicknesses in phase equilibria for the pure components. These were taken from molecular dynamics simulations (see Ref. [134]). 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 [4] and density functional theory [22](a). Model predictions for the fcc-bcc phase-selection boundary also fall close to the experimental results [135] for the Fe-Ni system (see Fig. 21).

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 [68](c),(f), [74,75,136,137]. 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(φ)|∇θ|, where the orientation energy scale is h ∝ T , and θ(r, t) is a scalar field representing the local orientation angle relative to the laboratory frame. Accordingly, θ is a circular variable. As θ is a nonconserved field, the respective equation of motion reads as [136](c) where M θ 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 [136](c), whereas the second term differs from 0 if the interfacial free energy is not isotropic [s(θ) = 1]. Behavior of this type of EOMs has been addressed by Kobayashi and Giga in 1D and 2D, providing analytical solutions for testing [138]. Eq. (27) 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 [75](a), [19](a), yet a more physical approach that creates crystallites via adding noise to the equations of motion is also possible [19](a). Herein, we recapitulate the latter, which has been used to model a broad range of polycrystalline structures [68](c),(f), [137]. 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 ζ θ (r, t), ζ θ (r , t ) = ζ 2 θ,0 δ(r − r )δ(t − t ) is added to the right hand side of Eq. (27), where ζ θ,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(φ), 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 θ = M θ,0 [1 − p(φ)] is assumed, where M θ,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 [137].) 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 [139](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 (∼ 1 nm), in a constant temperature gradient. f ori of the model employed here [139](b) ensures a fix misorientation between the two solid phases. The morphology occurring in the simulation is characteristic to such processes [139](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 [9](c), [136](d), [140]. There are different mathematical representations of the crystallographic orientation in 3D, such as Euler angles, rotation matrices, Rodrigues vectors, quaternions, etc. [141]. Here, we present an approach developed in the quaternion formalism [19](c), [136](d), [142]. A 3D generalization of the 2D orientation field model described previously differs in f ori , which now takes the form 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.
where Gaussian white noises of amplitude ζ i = ζ S,i + (ζ L,i − ζ S,i )[1 − p(φ)] were added to the EOMs, a form that ensured that the quaternion properties of the q i fields were retained. (Here ζ L,i and ζ 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 [19](c) [136](d), [142], a similar quaternion based model was used for addressing grain coarsening [143].

B. 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 [87,88,95,97,98]. 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 [144]).
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 [144,145] with a shooting algorithm [144], or employing a relaxation method [144]. A few examples are given for the application of these methods in Refs.
[14](c), [19,99,100,102]. In solving the sets of ELEs for boundary value problems of multi-phase-field computations in 1D or 2D, the relaxation method was used [133]. In d ≥ 2, the use of an unstructured grid [19](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 [146]. Path finding methods were used by several authors to address homogeneous and heterogeneous nucleation [132,147].
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 [144,145,148], finite element [149], or spectral discretization [144,150,151]. The importance of the proper choice of numerical methods (e.g., spectral spatial discretization) was nicely demonstrated in Ref. [152]. Methods for producing uncorrelated/colored noise for non-conserved and conserved fields can be found in Ref. [153]. Application of an adaptive grid can significantly reduce the computational cost and the lattice anisotropy [154], 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.

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.
[40](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 [66](a), (d),(e),(f) and numerical [65] 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 [40](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 [66](b),(c), [155]. 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 ± 0.02 falling close to the theoretical expectation (p = 3) [136](c).
(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 ≈ 3 was observed [19](a), [68](c), [121], 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 ≈ 2(= d) [156]. Apparently, however, one can produce dendrites, whose side arms are underdeveloped or missing, in which case exponent p decreases with increasing Y [156] as in the case of elliptical needle crystals, [65] 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 [19](c), [157]. If the particles interact via their diffusion fields, p = 1 + d/2 = 2.5 is expected in 3D [40](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 → 1 + d = 4 is expected, i.e. the kinetic exponent shall be in the range 2.5 ≤ p ≤ 4. Indeed the values p = 2.99 ± 0.01 [19](c) and p = 3.21 ± 0.01 [157] 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 diffusioncontrolled mechanism. Apparently, soft impingement is expected to reduce p with an extent increasing with increasing Y ; an expectation supported by theory [158]  simulations in 2D [121].) Kinetics of crystallization in thin films was also investigated using a 3D model [157]. The Avrami-Kolmogorov exponent observed was, p = 2.37 ± 0.01 [157]. 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 [159-161] of microstructures that are qualitatively similar to the respective panels (a), however, the staistics is not sufficient to evalute the kinetic exponent p.

and PF simulations [68](c). (Analogous phenomena were observed in PF
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.

Large-scale modeling of polycrystalline solidification morphology
Recently extremely large scale (4096 × 4104 × 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) [162,163].

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 [164] 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 [165]. 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 [166], 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 [166]. In the simulations the matter traveled from right to left with a pulling velocity v ∈ [4 × 10 −4 , 32 × 10 −4 ] m/s in a linear temperature distribution of gradient G ∈ [5 × 10 3 , 4 × 10 4 ] K/m. At the upper and lower sides of the simulation box periodic boundary conditions were assumed. New ∆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. [42]. Foreign particles of a Gaussian size distribution with a mean diameter of 2R p = (29 ± 3) nm were placed randomly into the new layers of thickness ∆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 ∆T c = 2γ SL /(∆s f R p ). The thermodynamic properties of Al-Ti were taken from Ref. [168]. 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) [166,169]. 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 [170].

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" [136](a), and a broad (ii) Decreasing the temperature in molecular liquids, the rotational diffusion coefficient decreases faster than the translational one, decreasing the ratio M θ /M φ . 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 [136](b)- (d).
(iii) Fixed angle branching of needle crystals in directions of low grain boundary energies [136](c).
(iv) MD simulations indicate that crystal nuclei often contain multiple orientations due to the formation of stacking faults, multiple twinning, etc. [47,172].
(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 [173].
Most of these mechanisms have been captured by orientation field based PF simulations [136]. For example, they recovered the main features of the experimental orientation distribution in spherulites [175] (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) [176]. 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.

D. 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 γ SL ∝ 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, γ SL,eq ∝ 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 ( 2 ∝ T and w ∝ 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) [34](c),(d) [128] seems to work fairly well for both the hard sphere system and Ni (and for a range of condensation/solidification problems [128]). 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 [177]. 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 [178].
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 [21,22,77] achieved some accuracy in predicting the properties of the HS system [77](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 [24,25] (PFC), which has been employed recently for exploring structural aspects of crystal nucleation [179]. 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 [25]. 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 [101]. Varying the number of the retained expansion terms, a range of PFC type models can be deduced [24,25]. 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 ∇ 12 ψ 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.

A. 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 [25]. 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 [180].

The single-mode PFC model
This version of the PFC model was derived from the perturbative density functional theory of Ramakrishnan and Yussouff [77](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.
[24](b), [25,181]. The dimensionless free energy measured relative to a reference liquid can be expressed in terms of reduced temperature (> 0) and reduced density ψ, leading to a Swift-Hohenberg/Brazowskii type expression [182]: where all quantities (including r and the differential operators) are dimensionless. In 2D, the reduced temperaturereduced density ( vs. ψ) 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 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} [183,184]. The liquid phase becomes unstable against density fluctuations beyond a linear stability limit given by the relationship ψ s = −( /3) 1/2 . The liquidus and solidus lines emerge from a critical point at 0 = 0 and ψ 0 = 0, while the interfacial free energies and interface thicknesses scale with the respective mean field critical exponents [184]. Formation of an amorphous solid via first order phase transition has also been reported [185]. 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 [24,25]. Although there is another (equivalent) formulation of the single-mode PFC model [24](b), [185,186], it can be readily transformed into the Swift-Hohenberg/Brazowskii form [25,184].
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.
[24](a); whereas another route was proposed by Wu and Karma to fix the model parameters to data taken from molecular dynamics simulations [187]. Re-casting the respective expression for the free energy in terms of λ = R 1 /(1 + R 1 ) ∈ [0, 1] (where R 1 is the relative strength of the first-and second-mode contributions [24](c)), we obtain a form compatible with Eq. 30. This parameter can be used to interpolate between the two-mode (λ = 0) and single-mode (λ = 1) PFC formulations: Here is the scaled density difference relative to the reference liquid of particle density ρ ref L . Again, the reduced temperature 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 [24](c)) is the ratio of the wave numbers of the two modes.
For λ = 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 < λ ≤ 1, a bcc domain extending with decreasing λ appears [see Fig.  36(b)], possibly with an hcp domain between fcc and bcc domains, as in the single-mode PFC case.

The multi-mode PFC model
An N -mode version of the PFC approach was proposed by Mkhonta et al. [188], assuming that the two-particle direct correlation has N peaks and it is isotropic. Then the respective free energy was written as where Λ, b i , τ x , and ξ 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 [24](b), [25]. To address the competition of bcc and fcc phases Tang et al. [189] have chosen N = 3, Q 0 = 1, Q 1 = √ 2, and Q 2 = 3/2. With appropriate values of the other parameters, a phase diagram, which resembles that of the single-mode PFC model was obtained [189].

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 [190]: 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 Using input from MD simulations for Fe as in Ref. [187], 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 [190].

B. 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 (r)ψ can be incorporated into the free energy density, which acts in the area/volume occupied by the substrate [184,191]. The potential can be written in the form V (r) = [V s,0 − V s,1 S(a s , r)]f (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 , r) [184], while a s is its lattice constant. The shape of the substrate is defined by an envelope function f (r) ∈ [0, 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 where ψ 0 is the reduced particle density of the liquid, relative to which the free energy is to be calculated [184]. 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 [192], bcc-liquid [193], and the fcc-amorphous solid [192] interfaces, including the anisotropy of the bcc-liquid [193] 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 [184]. 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 = 0.5, where faceted crystal morphology develops [184].
The same approach has been used to evaluate the nucleation barrier for homogeneous nucleation [184,192], heterogeneous nucleation [184,191], and particle induced freezing [184,191] 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) [191].
A simplified string method has been used to find crystal nuclei in the single-mode PFC model for homogeneous [194] and heterogeneous [195] nuclei forming on a structureless substrate at small values, where the interface appears to be anisotropic only for small clusters [194]. 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).  (112) and (011) are parallel with the (10) surface of a square lattice substrate [191]: Misorientation dependence of (a) the nucleation barrier and (b) the contact angle. Data were evaluated from solutions of ELE obtained by the relaxation method [184,191]. Note the minima in the relative nucleation barrier height at special misfit values, for which the surface structures of the crystal and the substrate match. The red symbols indicate results for the simulations shown in the upper panels.  [194,195]. In panel (a), the spatial variation of the reduced density is shown at supersaturations decreasing from curve a, to f; whereas in panel (b) the reduced density map is shown for half of the nucleus formed on a structureless flat substrate at the bottom of the image (by the courtesy of R. Backofen).

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 where the fluctuations are represented by a colored Gaussian noise ζ of correlator ζ(r, τ )ζ(r , τ ) = −α∇ 2 g(|r − r |, σ)δ(τ − τ ). Here α is the noise strength, and g(|r − r |, σ) a high frequency cutoff function [153,191,192] for wavelengths shorter than the interatomic spacing (σ = π √ 6). This EOM can be deduced from the dynamical density functional theory of colloidal suspensions after making a number of approximations [181]. 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 ∝ t −1/2 [197]. Beyond a critical undercooling or supersaturation diffusionless growth takes place, displaying steady state growth [197].  [193] (filled symbols) with those for the fcc-liquid interfaces from solving the dynamic equation of DPFC [196] (empty symbols).

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 [198]. Remarkably, in the single-mode PFC model the fcc-liquid and bcc-liquid interfacial free energies fall quite close to each other (Fig. 41) [193]. 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 [198], where the size dependent interfacial free energy could be fitted well (see Fig. 42) by an expression derived [199] from a phenomenological diffuse interface theory (DIT) [128]. Here r is the radius of the nucleus (considering the surface of tension), δ is the characteristic interface thickness, and a a numerical constant.

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. [185]. 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. [184,192] 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 [185]. 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. [59](a) and Lechner and Dellago [59](b). It was shown [200] that on theq 4 vs.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. [61]), and bcc-like structures were observed, but with strongly time FIG. 47. Structural analysis of the two-step nucleation process shown in Fig. 44 [198]. Particles with amorphous, MRCO, and bcc-like neighborhoods are colored white, red, and green, respectively. (Reproduced with permission from Ref. [200] c 2017 Elsevier.)  [200]. 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 [192].
The EOF-PFC model provides a refined thermodynamic description for bcc systems [190]. With parameters taken from MD simulations for Fe [190], the pair correlation function g(r) for the amorphous nucleation precursor predicted by this model [184] is in a remarkable agreement with the g(r) from MD simulations [201] 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 [192].
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 [202]. 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 [203]. 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 [203].

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 [192]. Study of the latter phenomenon, however, turned out to be tractable [189] in the three-mode DPFC model [188] 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 parameterq 6 as follows. SRO:q 6 < 0.28, MRO: 0.40 >q 6 > 0.28, and LRO: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 [204] and Kurdjumov-Sachs [205] relationships). With increasing reduced temperature the following structural sequences were observed during nucleation (Fig. 49): close to the critical point ( = 0.05) the bcc phase forms from the MRO domains, at ( = 0.15) formation of the fcc phase from bcc precursor is reported, whereas at ( = 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.

Crystallization kinetics
The binary version of the single-mode DPFC model [24](b) was used to model crystallization kinetics in an undercooled alloy in 2D [19](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 [40](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 Shape of the adsorbed crystal on the surface of a cube shaped substrate of simple cubic structure, as predicted by the ELE, at the critical undercooling, beyond which free growth takes place [95]. The size increases from left to right: (a) L = 16a f cc and (b) L = 32a f cc . (c) For the sake of comparison the theoretical shape corresponding to L/a f cc → ∞ obtained by a numerical surface solver [43] is also shown.
(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.

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 [191]. 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 [191]. 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 [191]. The shape of the adsorbed crystal volume just before free growth depends strongly on the size of the foreign nanoparticle [206] (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 [42](b).

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 [179]. This approach allows the investigation of nucleation on substrates of more complex geometries (see e.g. Fig.  52) [202,206]. 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 [207] to model various pattern FIG. 52. Heterogeneous crystal nucleation and growth on crystalline substrates of complex shape in the single-mode DPFC model [206]. From top to bottom: the substrate has atomic scale ledges, zigzagging shape, and rectangular grooves. Heteroepitaxy on an fcc substrate with a rectangular pit within the DPFC model [191]. Cross-sectional views are presented. Order parameters q4 and q6 were used to determine the local structure. Hues from dark to light denote the substrate, the (metastable) fcc, bcc, and (metastable) amorphous structures, respectively. From left to right the ratio of the lattice constant of the substrate and the stable fcc phase are as/a f cc = 1.0, 1.098, and 1.42, respectively.
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) [191].

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 ∝ µ −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 [209]. It relies on the momentum transport and continuity equations of fluctuating nonlinear hydrodynamics [210] valid down to the nanometer scale [211]: where p(r, t) is the momentum field, ρ(r, t) the mass density, v = p/ρ the velocity, δρ the divergence of the reversible stress tensor, δF [ρ] δρ the functional derivative of the free energy with respect to density. ρ 0 stands for the reference density, and D = µ S {(∇ ⊗ v) + (∇ ⊗ v) T } + [µ B − 2 3 µ S ]I(∇ · v) the dissipative stress tensor, while µ S and µ 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: 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 =p/v, both in the convection and dissipation terms [207]. It has been shown that with model parameters corresponding to a hypothetical 2D Fe the model recovers steady state growth v ∝ µ −1 , a linear dispersion relation for the long wave acoustic phonons, and a reasonable stress relaxation [209]. This model was used to address homogeneous, heterogeneous, particle induced, and growth front nucleation using model with parameters from Refs. [209,212], corresponding to a hypothetical 2D Fe [209].
Other PFC-based models of freezing were put forward very recently, which rely on either simplified or full hydrodynamics [213,214], on a modified Chapman-Enskog procedure to derive the hydrodynamic equations [215], or on amplitude expansion when coupling the PFC model to the Navier-Stokes equation with additional terms in the continuity equation [216]. 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.

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 [62] 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.

Crystallization kinetics
In recent studies [217,218], 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 = 0.1158 and ψ 0 = −0.1982, on

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 [42], 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 [217]. Analogously to the approach used in the DPFC model, a term V (r)ρ(r, t) was introduced to the free energy density to model the foreign particle. Here V (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 < ξ < 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 ξ, and shows a transition from stable configurations to free growth.

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) [200]. 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 = 0.1158, and the initial density ψ 0 = −0.1982, while the reduced temperature at the stability limit was 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 = j exp{6iθ j }, was used to monitor the local structure, where the summation was for the nearest neighbors, while θ 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 [219,220], 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 [173,221]. A more detailed analysis is available in Ref. [217]. These modes of growth front nucleation lead to the formation of nanoscale spherulite-like polycrystalline morphologies [217], raising the possibility that similar processes may play a role in spherulitic solidification on a larger scale.

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 [179,184,[191][192][193], whereas the string/elastic band methods [146,147,220] were employed to find the nucleation barrier and the minimum free energy path for homogeneous and heterogeneous nucleation [194,195]. Other possible methods are reviewed in Ref. [146]. 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 [179,184].
In solving the dynamic equations for the DPFC and HPFC models various numerical methods were used, including the finite difference [24](b), finite element [223] and spectral approaches combined with various operator splitting strategies [224]. 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 [224].
Numerical problems associated with combining PFC type models with hydrodynamics were addressed in Refs. [214,225].
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 [24](d), [226]. 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 [227]), the quasi-crystalline, and the amorphous phases, in multiphase/multicomponent systems. In a few cases, quantitative modeling was made possible by specific PFC approaches [207](a), (b) [228]. 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 [229]. 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.
FIG. 58. Two-step crystal nucleation in phase separating liquid solution as predicted by the binary XPFC model. Note that crystal nucleation happens in the solute-rich droplets [see panels (d) and (e)]. (Reproduced with permission from Ref. [229] c 2017 American Physical Society.)

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 ∼ 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 [230] 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 [231]. 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 [232].
In a recent study, the single-mode DPFC model was used to study the kinetics of crystal nucleation in the undercooled state [233]. 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 [233] would be necessary to develop a quantitatively correct nucleation theory.
Lutsko and coworkers [22](a),(b), [234] 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 [235]. 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 [187,190,[226][227][228][229]236] 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 [237]. 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 [238].
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 ∼ 100 [4], while taking the nucleation barrier W * from (i) the droplet model of the CNT with γ = γ 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 [128]. 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 [42], for which the size-and perhaps the shape distribution of the foreign particles are needed. Handling of other, more complex situations [1](a), [239], may require an extended knowledge on the heterogeneities.
Summarizing, since establishment of classical nucleation theory more than 90 years ago [240] and its adaptation to crystal nucleation in undercooled liquids about a quarter of a century later [241], 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.