Beyond the Ginzburg-Landau theory of freezing:Anisotropy of the interfacial free energy in the Phase-Field Crystal model

Beyond the Ginzburg-Landau theory of freezing: Anisotropy of the interfacial free energy in the Phase-Field Crystal model


This paper re-visits the weakly fourth order anisotropic Ginzburg-Landau (GL) theory of freezing. First we determine the anisotropy of the interfacial free energy in the Phase-Field Crystal (PFC) model analytically, and prove that it remains finite at the critical point as a direct consequence of the one-mode dominance of the model. Next, we derive the leading order PFC amplitude model and show the formal analogy to traditional weakly 4th order anisotropic GL theories. We conclude that the material-independent anisotropy appearing in emergent GL theory coincides with the remnant anisotropy of the generating PFC model. As a result, we show that the reduced temperature does not enter into the interfacial free energy anisotropy for metallic materials in both the Phase-Field Crystal model and the emerging Ginzburg-Landau theories. Finally, we investigate the possible pathways of calibrating anisotropic Ginzburg-Landau theories.

61.50.Ah, 68.35.Md, 68.08.-p, 64.60.F-

I Introduction

The anisotropy of the crystal-liquid interfacial free energy is regarded as the key factor of dendritic solidification, since it determines the microstructure of the crystallizing material, including many commercial metallic alloys. Many attempts have been made to determine the shape and the value of the anisotropy of the interfacial free energy, including equilibrium shape measurements Napolitano et al. (2002); Liu et al. (2001); Napolitano and Liu (2004) and molecular dynamics simulations. Molecular dynamics-based methods, such as the cleaving technique Broughton and Gilmer (1986); Davidchack and Laird (2000); Baidakov et al. (2013); Davidchack and Laird (2005) and the capillary fluctuation method Hoyt et al. (2001); Morris (2002) predict the anisotropy in the order of for several metallic systems. (For bcc systems, see References Sun et al. (2004a, b).) Since it has been revealed that the anisotropy critically depends on the crystal symmetry, and its magnitude depends mostly on the ratio of the crystal-liquid interface thickness and the interatomic distance, continuum descriptions also can be relevant tools for describing the anisotropic properties.

The first order parameter theory that captures anisotropy was developed by Haymet and Oxtoby Haymet and Oxtoby (1981); Oxtoby and Haymet (1982). The description is based on the classical Density Functional Theory (DFT) of freezing of the Ramakrishnan-Yussouff type Ramakrishnan and Yussouff (1979), which chareacterizes the system by the time-averaged local one-particle density. Since the theory works on the molecular scale in space, it inherently contains the crystalline symmetries of the system. Later a more convenient description, the Ginzburg-Landau (GL) theory of bcc-liquid interfaces was developed by Shih et al. Shih et al. (1987). In the GL theory the free energy of nonuniform phases is expressed in terms of space-dependent reciprocal lattice vector amplitudes, which are constant in the bulk phases and vary on the scale of the crystal-liquid interface thickness. The revised theory of Shih et al by Wu et al. Wu et al. (2006) predicts for iron. (In Reference Wu et al. (2006) the anisotropy parameter is defined as , where and are the interfacial free energies for the and crystal-liquid equilibrium planar interfaces, respectively.) This value is also supported by the simpler, DFT motivated Phase-Field Crystal (PFC) model Elder et al. (2002); Emmerich et al. (2012) and its amplitude theory Wu and Karma (2007), while two versions of the PFC model of Jaatinen et al Jaatinen et al. (2009) yielded (GL-PFC) and (Eight-order fit PFC), respectively.

Although the results of continuum theories are fair agreement with the experimental results and the results of atomistic simulations, both the -order GL and PFC amplitude theories of pure materials have a quite worrisome common property pointed out by Majaniemi and Provatas Majaniemi and Provatas (2009): they are ”weak” in a manner that all material parameters (except the crystal structure) scale out from the free energy functional. Consequently, the anisotropy parameters in these models depend exclusively on the crystal structure but not on the temperature, which results in a limited applicability of these models, and necessitates proper modifications. Such modifications may be including further reciprocal lattice vector sets and/or applying higher order polynomials in the free energy density Shih et al. (1987).

The starting point of developing consistent anisotropic Ginzburg-Landau theories is classical Density Functional Theory. The classical DFT inherently contains the crystal symmetries, and its amplitude expansions lead to particular Ginzburg-Landau theories. Since the PFC is a -order density functional theory with relatively simple spatial operators, it is a good candidate to employ for showing the relationship between the mathematical form of the anisotropy in a GL theory and how it emerges from an underlying classical DFT. In addition, the PFC amplitude theories show formal analogy to the anisotropic GL theories in the sense of the ”weak” nature, which seems to be more than just a coincidence.

The paper is organized as follows: In Section II we discuss the invariant formulations of the Phase-Field Crystal free energy functional. In Section III we calculate the equilibrium properties of the bulk (liquid and crystal) phases, and determine the properties (exponents and coefficients) of the equilibrium crystal amplitude and equilibrium density. Using the results, in Section IV we calculate the interfacial free energy, and prove that the anisotropy remains finite at the critical point, which is a direct consequence of the one-mode dominant behavior of the PFC. Finally, we derive the free energy functional of the anisotropic amplitude expansion of the PFC model in the leading order, and show that it is equivalent to a weakly fourth order Ginzburg-Landau theory. In section V we discuss the results.

Ii The Phase-Field Crystal model

In the first part we investigate the crystal-liquid equilibrium in the Phase-Field Crystal model introduced in Ref. Elder et al. (2002). After defining the free energy functional, we investigate the behavior of the PFC model close to the critical point, and prove that the first reciprocal lattice vector (RLV) set dominance of the model is related to the critical exponents of the RLV set amplitudes.

Figure 1: Direct correlation functions: (a) Schematic correlation function of a real system (gray) and typical Phase-Field Crystal correlation function (black). (b) Scaled PFC correlation function , where is the position of the maximum of the PFC (indicated by the horizontal dashed gray line in panel a). Note that the zero-valued minimum of at is independent from the particular form of the PFC .

ii.1 Minimal form of the free energy functional

In the single-component Phase-Field Crystal model the free energy of the system relative to a reference homogeneous state of density reads as Elder et al. (2002):


where is the scaled density field, and is a single-peaked direct correlation function in the wavelength space with peak position (see Fig 1.a). As a first step, we scale the model in order to identify the important parameters: Scaling the length as , the order parameter as and the free energy as results in a simplified form of Eq. (1):


The choice of and


results in the scales and , and the parameters

and . Here is an arbitrary scaling parameter: for example, choosing generates . Note that is a non-negative function with a single minimum at with (see Fig 1.b). This transformation of the direct correlation function will play a crucial role in our derivation. Taking into account that is an even function, it can be written as , which corresponds to in real space. (For the sake of simplicity, we won’t use from this point). Consequently, the term in Eq. (2) is equivalent to in the variational sense (note that both formulae results in the same functional derivative with respect to ). Using this equivalence, the cubic term can be eliminated: Substituting into Eq. (2) simply results in , while the terms up to the first order in can be neglected (since such terms vanish in both the Euler-Lagrange equation and the equation of motion). The ”minimal” form of the original free energy functional then reads as


where . This is a fairly simple form compared to Eq. (1) and shows that the important parameters of the model are only and .

ii.2 Periodic solutions

Eq. (4) generates a first order phase transition between homogeneous (liquid) and lattice periodic (crystal) solutions. These phases represent extrema of the free energy functional, therefore, they can be found by solving the Euler-Lagrange equation: by definition, where is the functional derivative of with respect to , and , i.e. the chemical potential of a homogeneous background liquid of density . Since the ELE is a nonlinear, higher order PDE, usually it is solved numerically. Instead, however, we can parametrize the lattice periodic solution in the following general form:


where is the average density, the amplitude of the RLV set, and the RLV in the RLV set. The bulk free energy density is defined as the volumetric average of the free energy in a unit cell:


where is the integrand of Eq. (4). For practical reasons we define the free energy density difference as:


Using the definitions (6) and (7), and substituting Eq. (5) into Eq. (4) together with results in (see Appendix A):


where we introduced the shorthand notation


where denoted here as the Kronecker-delta function , which gives 1 if the sum of the reciprocal lattice vectors in the argument is zero, otherwise it is 0. Therefore, is just the total number of N-term vector sums resulting in zero in which the first vector is from the RLV set , the second is from and so on. Consequently, is just the number of RLVs in the RLV set. Note that is invariant for the permutation of the indices.

ii.3 Equilibrium conditions

Eq. (8) realizes a parametrization of the free energy functional, which has to be minimized with respect to the set amplitudes and the selected wavelength at a constant average density . Introducing , where , the minimization equations read as:


From Eq. (10) two qualitatively different types of solutions emerge: (i) the trivial solution: for (homogeneous solution, the liquid phase), and (ii) a nontrivial lattice periodic solution (crystalline phase), where . Neglecting the crystal-liquid density jump for the crystal-liquid equilibrium is simply defined by equal free energy densities of the phases at the same average density, i.e.


where is the nontrivial solution. Eq. (11) together with Eq. (10) defines the atomic distance , the equilibrium solid amplitudes and the equilibrium density as a function of and .

ii.4 Critical behavior

In this section we show that the general PFC model described by Eq. (1) generates a mean-field Brazowskii/Swift-Hohenberg critical point at . We determine the critical exponents of the equilibrium density () and crystal RLV set amplitudes () and show that for any , implying the the one-mode dominance of the model.

Wavelength selection

For the particular choice Eq. (4) reduces to the well-known Brazowskii/Swift-Hohenberg form, which has a critical point at Elder and Grant (2004). It is reasonable to assume that this behavior doesn’t depend on the particular form of , and the model has a critical point as long as is a positive semidefinite function with a single, zero-value minimum at , i.e. . Indeed, it is relatively easy to see that the only solution of Eqns. (10) and (11) for is and . Therefore, we can write and for in general. In order to determine the critical exponents first we assume that there are more than one dominant RLV sets, meaning that , where and for all . Using this, the leading order term of Eq. (11) reads as:


where is the selected wavelength satisfying . Since and , Eq. (12) can be satisfied only if for all dominant RLV sets. Since has only one minimum at for which , only one RLV set can be dominant. In addition, this must be the first RLV set (thus ), since we’re searching for a crystal structure (in other words, the only dominant RLV set cannot be a harmonic). Moreover, since , the term has no effect on the phase diagram. This is in accordance with the original assumption, that the existence of the critical point doesn’t depend on the particular choice of . The critical point exists as long as and has a single minimum at with .

Critical exponents

Taking into account that for and using , the equilibrium condition reads as:


where we used the shorthand notations , and (details are shown in Appendix A). From Eq. (5) it is trivial that , otherwise, there is no first order transition for . In addition, in order to find nontrivial solution for and , the term in the free energy functional must contribute to the leading order. Taking these facts into account, the first row of Eq. (13) together with implies


therefore, the leading order of Eq. (8) is . In the next order of Eq. (13) (the second and the third lines) the minimization equations for are decoupled:


resulting in


on the same basis, therefore, the next order of Eq. (13) is proportional to . In addition, from it can be shown that , therefore, the first correction from this in Eq. (13) is in the order of . This means that our calculation is self-consistent.

Finally, one can determine the coefficients and by substituting , , and into Eq. (8) then taking the leading order of Eqns. (10) and (11). The equations then can be solved analytically for and :


showing that the leading order equilibrium density and crystal amplitude depend exclusively on the crystal structure (apart from , naturally).

It is noteworthy that our results stay valid when the equilibrium density jump is considered in the calculations (for details, see Appendix B).

Iii Interfacial free energy

In this section first we define the crystal-liquid interfacial free energy in the Phase-Field Crystal model, then we will approximate it analytically by using the results of the previous section. Considering the isotropic case first, we determine the interface thickness(es) and the interfacial free energy, and their critical exponents. As a key contribution of this work, we prove that the one-mode dominance of the PFC model, shown in the previous section, results in a remnant equilibrium crystal-liquid interfacial free energy anisotropy at the critical point. In the final part of this section we will determine the remnant anisotropy for the bcc structure and verify the result by comparing it to the results of numerical solutions of the Euler-Lagrange equation.

iii.1 Definition of the anisotropic crystal-liquid interfacial free energy

When the density jump between the equilibrium crystal and liquid is neglected, the anisotropic interfacial free energy reads as


where is the normal of the planar crystal-liquid interface, the orthogonal distance from the interface, while denotes an average calculated for a plane parallel to the interface at a constant value of . The integrand of Eq. (19) reads as

Here represents the equilibrium crystal-liquid density distribution, where is approximated as


where (the equilibrium crystal amplitudes) and we have used the following shorthand notations:


where is the characteristic interface width of the plane wave in the RLV set Majaniemi and Provatas (2009). Note that far from the interface Eq. (20) recovers the density distribution of the equilibrium bulk phases: and . Using Eqns. (20) and (5) in Eq. (75), after a straightforward but lengthy algebra one can come to a reasonably simple parametrized form of the leading order anisotropic crystal-liquid interfacial free energy (for details, see Appendix C):


where the sums for run for all different pairs in and , respectively, while we used the shorthand notation . The function


is responsible for the anisotropic contribution [here and are constants emerging from the particular form of ]. For example, for the theory (Brazowskii/Swift-Hohenberg), . Note, that the appearance of the anisotropic contribution to the leading order of is the consequence of the one-mode dominance of the theory, i.e. for any .

iii.2 Critical exponent of the interface thickness

Close to the critical point the interface thickness (correlation length) diverge as , where . Note that all interface thicknesses diverge with the unique critical exponent (for details, see Appendix D). In case of the isotropic limit (), Eq. (22) reads as:




is constant for geometrical reasons, and we used that and . Using (where is a constant specific to the isotropic case) in the minimization equation yields




Using these in Eq. (24) the isotropic interfacial free energy reads as


Note that the particular form of appear exclusively in the constant . Moreover, scales as with from Eq. (3), and as (since ), which results in the simple scaling relation


where and denote the isotropic interfacial free energy at and an arbitrary , respectively. Eq. (29) shows that is not a relevant parameter of the theory, and only helps to choose a convenient form of .

iii.3 Critical behavior of the anisotropy

Using the critical exponents and the facts that and in Eq. (22) yields


where the indices denote isotropic and anisotropic contributions, respectively. The anisotropy parameter is defined as


Applying Eq. (30) in Eq. (31) results in


where are defined by and , respectively. From Eq. (22) one can see that . However, for in case of in , therefore, the anisotropy is seen to be finite at the critical point:


which apparently contradicts to former expectations of Podmaniczky et al. Podmaniczky et al. (2014). Note that the remnant anisotropy () is a direct consequence of the one-mode dominant nature of the free energy functional: together with may yield a non-vanishing anisotropic contribution to the leading order of to the interfacial free energy.

Figure 2: Crystal-liquid interfacial free energy anisotropy for the (Brazowskii/Swift-Hohenberg) model. For , .

We have to mention at this point that it would also be useful to investigate the critical behavior of the anisotropy in the presence of fluctuations. Fluctuations destroy the mean-field behavior, and we know that the anisotropy vanishes at the critical point in the triangular Ising system Shneidman and Zia (2001), however, we also know that the Brazowskii system has its own universality class Brazowskii (1975); Fredrickson and Binder (1989).

iii.4 Verification of the remnant anisotropy

To determine in Eq. (33), first we calculate Eq. (22) divided by Eq. (28) for the general anisotropic case using a reasonable approximation of the envelope function integrals, which are defined as and in Eq. (22) and approximated in detail Appendix D. It yields the coupled minimization equations for interface thickness constants relative to the isotropic one, i.e. :


which have to be solved numerically for for the model. For the bcc structure preferred by the Phase-Field Crystal model close to the critical point in 3 dimensions , and . We started the numerical calculations from the isotropic solution defined by Eq. (27) for the [111] and [100] crystal planes, which give the minimal and the maximal interfacial free energies, respectively. Our calculation resulted in a significant remnant anisotropy


For comparison, following the method of Podmaniczky et al. Podmaniczky et al. (2014) we evaluated the interfacial free energy by solving the Euler-Lagrange equation numerically for bcc-liquid equilibrium interfaces at and . We have found , a nearly constant anisotropy parameter, which is in a fair agreement with the analytical result, and moreover, is in a perfect agreement with the results from the GL theory of Wu et al. Wu et al. (2006) or the PFC amplitude equations of Wu and Karma Wu and Karma (2007). This unexpected coincidence, however, suggests a deeper relationship between the weakly order Ginzburg-Landau / amplitude models of classical density functional theories having a critical point.

Iv Connection to Ginzburg-Landau theories

In this section we will investigate the connection between the critical behavior of the Phase-Field Crystal model and Ginzburg-Landau theories. First we derive an isotropic amplitude model from the critical PFC model (i.e. in the leading order in case of ), then - following the recent work of Provatas and Majaniemi - extend it for the anisotropic case. A key result of this paper is to formally show that the leading-order amplitude model of the PFC close to the critical point is analogous to a weakly order anisotropic one-mode Ginzburg-Landau theory, and the material parameter independent interfacial free energy anisotropy appearing in the GL theory is precisely the critical point remnant anisotropy inherited from the generating density functional theory.

iv.1 Isotropic limit

Ginzburg-Landau polynomial

In equilibrium one can define the normalized amplitudes , where and denotes the equilibrium amplitudes: and . Note that for a planar equilibrium interface , respectively. With this re-scaling, the equilibrium bulk liquid and solid phases are described by and , respectively. Considering only the leading order terms of Eq. (13) and substituting and yields


where and are defined by Eq. (17) and (18). Substituting these into Eq. (36) yields




Note that Eq. (37) is exactly the well-known order Ginzburg-Landau polynomial for triangular and bcc structures. (For the fcc structure , therefore, there is no fcc-liquid first-order phase transition in the Swift-Hohenberg formalism in leading order, i.e. close to the critical point.)

Amplitude equation

The isotropic single order parameter amplitude equation in equilibrium can be written as:


where is defined by Eq. (37) and . The equilibrium solution of the Euler-Lagrange equation is the kink-function , where . The interfacial free energy can be obtained by using the integral Euler-Lagrange equation: . The model parameters and can be then related to the interfacial free energy and interface thickness as:


Substituting Eqns. (28) and (27) into the above equation yields


Note that Eq. (41) consistently recovers Eq. (37), showing that our calculation is self-consistent. Also note that Eq. (39) consistently verifies the divergence of the interface thickness found in Eq. (24).

iv.2 Anisotropic extension

Introducing the order parameters gives the anisotropic extension of Eq. (36):


Following Majaniemi and Provatas Majaniemi and Provatas (2009), for a planar interface the anisotropic interfacial free energy can be written as:


where is defined by Eq. (43). Here and are defined by Eqns. (41) and (42) again. The elements of the coefficient matrix can be determined by substituting into Eq. (44), and comparing the result with Eq. (22) after substituting , with Eqns. (17) and (18), and considering Eqns. (41) and (42) (for details, see Appendix E). The calculation then yields a diagonal matrix with elements . Finally, the corresponding anisotropic Ginzburg-Landau free energy functional then reads as:




while , and are defined by Eqns. (43), (41) and (42), respectively. Note that and scale out from the free energy functional [and also from Eq. (91)], therefore, the anisotropy of the interfacial free energy is constant and depends exclusively on the crystal structure [the form of ], while its magnitude is just the magnitude of the critical point remnant anisotropy of the Phase-Field Crystal model [Eq. (91) shows the leading order of the interfacial free energy close to the critical point].

V Conclusions

Our calculations show that a weakly -order anisotropic one-mode Ginzburg-Landau theory inherits the properties of a leading-order amplitude model of a one-mode dominant -order classical Density Functional Theory close to its critical point. The constant anisotropy appearing in weakly order Ginzburg-Landau theories originates from the fact that all material parameters (except the crystal structure) scale out from the free energy functional in the determination of the crystalline anisotropy. We have shown that the magnitude of the GL anisotropy coincides with the critical point remnant anisotropy of the generating density functional theory. We have to emphasize that the non-vanishing behavior of the anisotropy doesn’t contradict to the mean-field theory, since the anisotropy is a secondary quantity, i.e. its critical exponent is not related to the important exponents.

Our results have consequences on the quantitative applicability of both the Phase-Field Crystal model and Ginzburg-Landau theories emerging from it. In the case of the PFC model the numerical calculations resulted in a remnant () anisotropy in the range . In this range , where is the usual interface thickness and the bcc lattice constant. Since this is true for simple metals, is not a relevant parameter in quantifying the anisotropy for metallic materials. In contrast, it has been found that inherited by the GL theory exclusively depends on the form of the scaled direct correlation function . Since the symmetry breaking of the GL coefficient matrix is trivially related to properties of the direct correlation function, one can calibrate the anisotropy in the Ginzburg-Landau theory by investigating the critical behavior of the generating PFC.

A possible pathway of deriving consistent GL theories, in accordance with the original idea of Shih et al. Shih et al. (1987), is to choose such a PFC description, in which more than one RLV set is dominant, i.e. we at least two peaks of the direct correlation function are considered. The best candidate is the so-called structural PFC (or XPFC) model Greenwood et al. (2010), in which the peak peak heights are weighted by the Debye-Waller factor. Since the peak heights are not equal, the critical point vanishes, meaning that the dependence appears in the amplitude theory. Nevertheless, combining the XPFC model with the recently published fluctuating hydrodynamic theory of freezing Tóth et al. (2014) might result in a continuum description of crystallization of simple liquids on the (classical) fundamental length scale of the material. Moreover, comparing the results of the model with molecular dynamics data will hopefully anchor to the physical temperature, making the model fully quantitative.


The authors wish to thank professor L. Gránásy and F. Podmaniczky from the Wigner Research Centre for Physics, Hungary, for the valuable discussions which significantly contributed to the quality of the work. This work has been supported by the Postdoctoral Programme of the Hungarian Academy of Sciences and the Natural Sciences and Engineering Research Council of Canada.

Appendix A: Evaluation of the bulk free energy density

In order to evaluate Eq. (7) for , where , first we re-formulate Eq. (4) as follows:


where we used that the functional derivative


results in the same for both and . The spatial derivatives of read as:


where . Introducing the shorthand notation for the lattice cell average the following terms emerge from in the free energy density:




is the (Kronecker) delta-function giving for , and otherwise. Therefore,