Ginzburg-Landau theory of crystalline anisotropy for bcc-liquid interfaces Kuo-An Wu1, Alain Karma1, Jeffrey J. Hoyt2 and Mark Asta3 1Department of Physics, Northeastern University, Boston, Massachusetts 02115 2Sandia National Laboratories, Albuquerque, New Mexico 87185 3Department of Chemical Engineering and Materials Science, Center for Computational Science and Engineering, 6 University of California Davis, Davis, CA 95616 0 0 The weak anisotropy of the interfacial free-energy γ is a crucial parameter influencing dendritic 2 crystal growth morphologies in systems with atomically rough solid-liquid interfaces. The physical n origin and quantitative prediction of this anisotropy are investigated for body-centered-cubic(bcc) a formingsystemsusingaGinzburg-Landautheorywheretheorderparametersaretheamplitudesof J densitywavescorrespondingtoprincipalreciprocallatticevectors. Wefindthatthistheorypredicts 1 the correct sign, γ100 > γ110, and magnitude, (γ100−γ110)/(γ100+γ110) ≈ 1%, of this anisotropy 1 ingood agreement with theresultsofMD simulationsfor Fe. Theresultsshowthatthedirectional dependenceoftherateofspatialdecayofsoliddensitywavesintotheliquid,imposedbythecrystal ] structure, is a main determinant of anisotropy. This directional dependence is validated by MD i c computations of density wave profiles for different reciprocal lattice vectors for {110} crystal faces. s Our results are contrasted with the prediction of the reverse ordering γ100 < γ110 from an earlier - l formulation of Ginzburg-Landau theory [Shih et al.,Phys. Rev. A 35, 2611 (1987)]. r t m PACSnumbers: 64.70.Dv,68.08.De,68.70.+w,81.30.Fb . t a I. INTRODUCTION magnitude of the crystalline anisotropy has been tradi- m tionallycharacterizedbycomparingthevaluesofγcorre- - d The advent of microscopic solvability theory [1, 2, 3] sponding to {100} and {110} crystal faces. MD calcula- n in the 1980s lead to the prediction that the anisotropy tionshaveyieldedanisotropyvalues(γ100−γ110)/(γ100+ o of the excess free-energy of the crystal-melt interface is γ110) = 0.5−2.5%, and experimental values extracted c from equilibrium shape measurements fall generally [ a crucialparameterthatdetermines the growthrate and morphology of dendrites, which shape the microstruc- within this range. What determines physically the pos- 2 itive sign and the magnitude of this anisotropy, how- tures of many commercial metallic alloys. This predic- v ever, remains unclear. One interesting clue is that tionwaslargelyvalidatedbyphase-fieldsimulations[4,5] 7 MD-calculated anisotropies generally depend more on of dendritic solidification during the 1990s. More recent 5 the crystal structure than on the microscopic details of 0 work in the present decade has focused on the quantita- inter-molecularforcesforthesamecrystalstructure,and 1 tivepredictionofboththemagnitudeandtheanisotropy anisotropies tend to be consistently smaller for bcc than 0 ofthe interfacialfree-energyγ using moleculardynamics 6 (MD) simulations [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, for fcc [20]. A striking example of the former is the fact 0 that MD studies of Fe [16, 17] and succninonitrile [19], 17, 18, 19, 20, 21, 22]. In parallel, experimental progress / with the same bcc structure but entirely different inter- t has been made to determine this anisotropy in metallic a molecular forces, have yielded comparable anisotropies systems from accurate equilibrium shape measurements m around half a percent. The weaker anisotropy of γ for [23, 24, 25], which extend pioneering measurements of - bcc compared to fcc is also consistent with experimen- d this anisotropy in transparent organic crystals [26, 27]. tal measurements of anisotropy values in the range of n MD-based methods, including the cleaving technique 0.5−0.7% and 2.5−5% for the bcc and fcc transparent o [6, 7, 8, 9] and the capillary fluctuation method (CFM) c organic crystals succinonitrile and pivalic acid, respec- [10, 14], have been successfully developed to compute γ : tively [27, 28]. v andtoaccuratelyresolveitsnotoriouslysmallanisotropy Xi of the order of 1%. These methods have been applied to The fact that anisotropy appears to depend more awiderangeofsystems,includingseveralelementalface- strongly on crystal structure than inter-molecular forces r a centered-cubic (fcc) [10, 11, 14, 16, 17], body-centered- suggests that it may be possible to predict this critical cubic (bcc) [16, 17, 20], and hexagonal-close-packed[18] parameterfrom a continuum density wavedescriptionof metals,aswellasonefccmetallicalloy[13]modeledwith the solid-liquid interface, which naturally incorporates interatomicpotentialsderivedfromthe embedded-atom- anisotropy because of the broken symmetry of the solid. method (EAM), the Lennard-Jonessystem [8, 15], hard- ThesimplestofsuchdescriptionsistheGinzburg-Landau sphere [7, 21, 22] and repulsive power-law potentials [9], (GL)theoryofthebcc-liquidinterfacedevelopedbyShih and,mostrecently,abccmolecularorganicsuccinonitrile et al. [29]. The order parameters of this theory are the [19] used extensively in experimental studies of crystal amplitudes of density waves corresponding to the set of growth patterns. principal reciprocal lattice vectors {K~ }, and the free- i In systems with an underlying cubic symmetry, the energy functional is derived from density functional the- 2 ory (DFT) [30, 31, 32] with certain simplifying assump- largerK~ modedenotedby“...”isneglected. Expanding tions. the free-energy as a power series of the u ’s around its i Thistheoryhasyieldedpredictionsofγ forvariousbcc liquid value Fl ≡F[n0] yields elementsthatareinreasonablygoodagreementwithex- periments. Moreover,it has providedan elegantanalyti- n0kBT calderivationoftheproportionalitybetweenγandthela- ∆F = 2  d~ra2 cijuiuj δ0,K~i+K~j tentheatoffusion,inagreementwiththescalingrelation Z Xi,j rγonb−o2r/a3te=dαbLyMprDopsoismeudlabtyioTnusr[n9b,u1l6l][,3w3]haenrednreicsetnhtelyscoolird- −a3 cijkuiujuk δ0,K~i+K~j+K~k (2) i,j,k atomic density and L is the latent hear per atom. This X theory,however,predictsthewrongorderingγ100 <γ110. +a4 cijkluiujukul δ0,K~i+K~j+K~k+K~l This apparent failure could be due to the truncation of i,j,k,l X larger |K~| modes, i.e. density waves corresponding to du 2 i preacriapdrooxcaiclallalyt,tiacestvreocntgordseKp~enwdietnhce|K~o|f>the|K~ain|i.soHtroowpeyveorn, +b Xi ci (cid:12)(cid:12) dz (cid:12)(cid:12) !, largerK~ modeswouldbehardtoreconcilewiththeweak where we have defin(cid:12)(cid:12)ed ∆(cid:12)(cid:12)F ≡ F −Fl and the gradient dependence of this quantity on details of inter-molecular square terms arise from the spatial variation of the or- forces seen in MD studies, since these forces dictate the der parameters along the direction normal to the inter- amplitudes of these modes in the crystal where the den- face parameterized by the coordinate z. The Kronecker sity is sharply peaked around atomic positions. delta δm,n, which equals 0 or 1 for m 6= n or m = n, To shed light on this paradox, we revisit the simplest respectively, enforces that only combinations of princi- GL theory of the bcc-liquid interface based on the min- pal reciprocal lattice vectors that form closed polygons imal set of principal reciprocal lattice vectors. Our cal- K~ + K~ + ··· = 0 contribute to the free-energy func- i j culation differs principally from the one of Shih et al. tional. Closed triangles generate a non-vanishing cubic [29] in the derivation of the coefficients of the gradient term that makes the bcc-liquid freezing transition first squaretermsintheGLfree-energyfunctional. Eachterm order. The multiplicativefactorsa andbareintroduced i measuresthe free-energycostassociatedwiththe spatial since it is convenient to normalize the sums of the c’s variation, in the direction nˆ normal to the solid-liquid to unity (i.e. c = 1, c δ = 1, etc). To i i i,j ij 0,K~i+K~j interface, of a subset of equivalent density waves with completethe theory,one needsto determine allthe coef- the same magnitude of the direction cosine Kˆi ·nˆ and ficients appearPing in the GPL free energy functional. hence the same amplitude. Shih et al. choose these co- The coefficients of the quadratic terms are obtained efficients to be proportional to the number of principal from the standard expression for the free-energy func- reciprocallattice vectorsin eachsubset. This procedure, tional that describes small density fluctuations of an in- however, turns out to yield an incorrect directional de- homogeneous liquid pendence (i.e. dependence onKˆ ·nˆ) ofthe rateofdecay i ′ k T δ(~r−~r ) of density waves into the liquid. Here we find that the ∆F = B d~rd~r′δn(~r) −C(|~r−~r′|) δn(~r′)(3) inclusion of the correct directional dependence, as pre- 2 n0 Z Z (cid:20) (cid:21) scribed by DFT, yields the correct ordering γ100 > γ110 where δn(~r) ≡ n(~r) − n0 and C(|~r −~r′|) is the direct andareasonableestimateofthemagnitudeofanisotropy. correlationfunction whose fourier transform Furthermore, we validate this directional dependence by MD computations of density wave profiles. C(K)=n0 d~rC(|~r|)e−iK~·~r (4) Z −1 is related to the structure factor S(K)=[1−C(K)] . II. GINZBURG-LANDAU THEORY The two expressions for ∆F, Eqs. (2) and (3), can now be related by assuming that the amplitudes of den- We review the derivation of the GL theory and com- sity waves vary slowly across the interface on a scale pare our results to MD simulations in the next section. ∼ 1/K where K is the value of K corresponding max max The theory is derived in the DFT framework where the to the peak of the structure factor. Accordingly, u (z′) i free energy of an inhomogeneous system is expressed as can be expanded in a Taylor series about z a functional F =F[n(~r)] of its density distribution n(~r), du (z) which can be expanded in the form ′ i ′ δn(~r ) ≈ n0 ui(z)+ (z −z) dz i (cid:20) X n(~r)=n0 1+ ui(~r) eiK~i·~r+... (1) +1d2ui(z)(z′−z)2+... eiK~i·~r′, (5) i ! 2 dz2 X (cid:21) where the order parameters u ’s are the amplitudes of where the contribution “...” involving higher-order i density waves corresponding to principle lattice vectors derivatives can be neglected. Namely, terms propor- h110iofthereciprocalfcclattice,andthecontributionof tional to (z′ − z)ndnu (z)/dzn ∼ 1/(K w)n, where i max 3 w is the characteristicwidth of the solid-liquidinterface, the corresponding value of (Kˆ ·nˆ)2 for the correspond- i i.e. the scale over which order parameters vary from a ingsubsetI,II,orIII,respectively. Thesecoefficientsare constant value in the solid to zero in the liquid. Hence, c =1/16 for subset I, c =1/4 for subset II, and c =0 i i i these terms vanish at large n under the assumption that forsubsetIII.Thesecoefficientsyieldthegradientsquare w ≫ 1/Kmax. Substituting this expression in Eq. (3) terms −C′′(|K~110|)|du/dz|2 and −C′′(|K~110|)|dv/dz|2 in ′ and carrying out the integral over~r , we obtain theGLfreeenergyfunctional(2)sincethereare8equiv- alent reciprocal lattice vectors in subset I and 2 in sub- n0kBT 1 set II, respectively. The coefficient of |dw/dz|2 vanishes ∆F ≈ d~r u u δ 2  S(|K~ |) i j 0,K~i+K~j since principal reciprocallattice vectors in subset III are Z Xi,j i orthogonalto nˆ and ci =0.  2 1 du − C′′(|K~ |)(Kˆ ·nˆ)2 i , (6) In contrast, Shih et al. choose the c ’s to be equal 2 i i dz i i (cid:12) (cid:12) # for all subsets with a non-vanishing direction cosine X (cid:12) (cid:12) where C′′(K)≡d2C(K)/dK2. Compa(cid:12)(cid:12)ring(cid:12)(cid:12)Eqs. (6) and (subsets I and II), and ci = 0 for subsets with prin- cipal reciprocal lattice vectors orthogonal to nˆ (subset (2), we obtain at once III). Since there is a total of 10 reciprocal lattice vec- 1 tors in subsets I and II, the normalization condition a2cij = S(|K~110|), (7) theicgira=die1ntysiqeuldasrectier=ms1−/1(08./5)TCh′′e(s|eK~1c1o0e|ffi)|dcuie/ndtzs|2yaienldd P −(2/5)C′′(|K~110|)|dv/dz|2 in the GL free energy func- 1 bci =− C′′(|K~110|)(Kˆi·nˆ)2, (8) tional(2),whichareweightedproportionallytothenum- 2 ber ofreciprocallattice vectorsineachsubset,anddiffer where we have used the fact that all reciprocal lattice from the correct terms derived above. vectors have the same magnitude |K~i| = |K~110|. Sum- For the {100} and {111} crystal faces, the weighting ming both sides of Eq. (7) and using the normalization procedure ofShih et al. andEq. (11) give coincidentally c δ =1 gives i,j ij 0,K~i+K~j the same coefficients of the gradient square terms. Thus P a2 = δ0,K~i+K~j = 12 , (9) tthheesdeiffcaerseenstncereydstnaoltfabceesraepreeastuemdmhaerriez.edTihneTraebsuleltIs.for i,j S(|K~110|) S(|K~110|) X The determination of all the other coefficients in the and cij = 1/12. Similarly, summing over i both sides of GL free energy functional is identical to the calculation Eq. (8), and using the normalization ici =1, yields of Shih et al.. The coefficients of the cubic and quar- tic terms, c and c , respectively, are determined by b=−1 C′′(|K~110|)(Kˆi·nˆ)2 =−P2C′′(|K~110|), (10) the ansatz itjhkat all pijokllygons with the same number of 2 sides have the same weight, which yields c = 1/8 and i ijk X c =1/27;for quadratic terms, this ansatz reproduces ijkl where the second equality can be shown to be indepen- theresultc =1/12derivedabovesincetherearetwelve ij dent of the direction of nˆ, and two-sidedpolygonsformedbytheprinciplereciprocallat- 1 ticevectors. Usingthesecoefficientsandidentifyingeach ci = 4(Kˆi·nˆ)2 (11) ui with the order parameter u, v, or w, depending on whether the corresponding K~ on one side of a polygon i To make the difference between the above derivation belongs to subset I, II, or III, respectively, Eq. (2) re- of the c ’s and the one of Shih et al. explicit, con- i sider one of the {110} crystal faces with nˆ pointing in the [110] direction. The set of 12 principal recip- rocal lattice vectors {K~ } corresponding to the h110i i directions can be separated into three subsets with TABLEI: Comparisonofcoefficientsofsquaregradientterms the same value of (Kˆ · nˆ)2: subset I with 8 vec- tors ([011],[0¯11],[01¯1],[1i01],[¯101],[10¯1],[0¯1¯1],[¯10¯1]) and ci predicted by Eq. (11) (DFT) and Shih et al. [29] for the {100}, {110}, and {111} crystal faces. For each orientation, (Kˆi · nˆ)2 = 1/4, subset II with 2 vectors ([110],[¯1¯10]) the12principalreciprocallatticevectoraregroupedintosub- and (Kˆ ·nˆ)2 = 1, and subset III with 2 vectors (¯110, sets where Kˆi·nˆ havethesame magnitudein each subset. i [1¯10]) and(Kˆi·nˆ)2 =0. Density wavesin a givensubset 100 110 111 havethesameamplitudedenotedherebyu,v,andw for (Kˆi·nˆ)2 0 1/2 1/4 1 0 0 2/3 subsets I, II and III, respectively. Numberof K~i’s 4 8 8 2 2 6 6 It follows that the correct coefficient of the gradient ci (Eq. 11) 0 1/8 1/16 1/4 0 0 1/6 squaretermsforagivenorderparameteru,v,orw,isob- tained using the expressionfor c givenby Eq. (11) with ci (Ref. [29]) 0 1/8 1/10 1/10 0 0 1/6 i 4 duces for {110} crystal faces to 30 MD MH(SA)2 potential 25 Ginzburg Landau (present calculation) ∆F = n0k2BT Z d~r(cid:20)a2(cid:18)32u2+ 16v2+ 61w2(cid:19) -3m)20 1 1 12 1 1 c 15 −a3 (cid:18)2u2v+ 2u2w(cid:19)+a4(cid:18)27u4+ 27v4+ 27w4 -220 10 1 + 4 u2v2+ 4 u2w2+ 1 w2v2+ 4 u2vw n ( 5 27 27 27 27 (cid:19) 2 2 0 du dv −C′′(|K~110|) −C′′(|K~110|) , (12) (cid:12)dz(cid:12) (cid:12)dz(cid:12) # 100 110 120 130 140 150 (cid:12) (cid:12) (cid:12) (cid:12) z (Å) (cid:12) (cid:12) (cid:12) (cid:12) The corresponding expr(cid:12)essi(cid:12)on of Shih et a(cid:12)l. d(cid:12)iffers by the coefficients of |du/dz|2 and |dv/dz|2 that have an FIG. 1: Comparison of planar densityprofiles n(z) from MD extra multiplicative factor of 8/5 and 2/5, respectively, simulations(solidline)andthepresentGinzburg-Landauthe- as discussed above. Their expressions for ∆F for the ory (dashed line) for {110} crystal faces. {100} and {111} crystal faces are identical to ours since Eq. (11)andtheequalweightansatzyieldcoincidentally the same c ’s for these faces. [17], the amplitudes were computed by averaging over i Finally, the coefficients a3 and a4 are determined by many configurations the instantaneous value of a planar the constraints that the equilibrium state of the solid is structure function (i.e., the magnitude of the complex aismthineimvaulumeoofffarleletehneeorgrdy,er∂∆paFr/a∂mueit|euris=uins =the0,swolhide,reanuds Fthoeuraimerpclioteuffidecsieonftsdeonfstithyewdaevnessitsya)t.uWratitehtothaissmapapllronaocnh- that solid and liquid have equal free energy at the melt- vanishingvalueintheliquid. Theseamplitudes,however, ing point, ∆F(u ) = 0. These two constraints yield the aregenerallyexpectedtovanishintheliquidthathasno s relations a3 =2a2/us and a4 =a2/u2s that determine a3 long range order, consistent with the GL theory. and a4 in terms of a2 given by Eq. (9), which completes To calculate amplitudes that vanish in the liquid, the the determination of all the coefficients. followingprocedureisfollowed. Wefirstcomputetheav- erage number density n(~r)=n(x,y,z), with z measured fromafixedreferenceplaneofatomsinthesolid. During the MD simulations,only thoseconfigurationswhere the III. RESULTS AND COMPARISON WITH MD solid-liquidinterfacehasthesameaveragepositionalong SIMULATIONS z were considered. As described in detail by Davidchack and Laird, [35] this procedure avoids an artifical broad- The order parameter profiles for {110} crystal faces ening of the density profiles due to either the natural were calculated by minimizing ∆F given by Eq. (12) fluctuations in the average position of the interface or with respect to the order parameters u, v and w, and by Brownianmotionofthecrystal. Theinterfacepositionis solving numerically the resulting set of coupled ordinary foundbyfirstassigningtoeachatomanorderparameter differential equations with the boundary condition u = proportional to the mean square displacement of atoms v = w = u in solid and u = v = w = 0 in liquid. The s from their positions on a perfect bcc lattice. (The or- value of γ110 was computed using Eq. (12) with these derparametercalculationisthesameasthatusedinthe profiles. The same procedurewas repeatedfor the {100} capillaryfluctuationmethodandisdescribedinmorede- and {111} crystal faces. tail in reference [16]). Then an order parameter profile WeusedinputparametersfortheGLtheorycomputed as a function of z is computed by averaging within the directly from the MD simulations in order to make the x−y plane and the interface position is that value of z comparison with these simulations as quantitative and where the averaged order parameter is midway between precise as possible. These parameters include the peak the bulk liquid and bulk solid values. of the liquid structure factor ≈ S(|K~110|), which yields Next, we compute the x−y averageddensity a2 = 3.99 using Eq. (9), C′′(|K~110|) = −10.40 ˚A2, and theamplitudeofdensitywavescorrespondingtoprincipal 1 Lx Ly n(z)= dxdyn(~r), (13) reciprocal lattice vectors in the solid us =0.72. LxLy Z0 Z0 The MD simulations were carried out using the EAM potentialforFefromMendelev,Han,Srolovitz,Ackland, which is illustrated in Fig. 1. Lastly, we calculate the Sun and Asta (MH(SA)2) [34] and the same thermody- amplitude of density waves from the fourier transform namic ensemble and geometries as in Ref. [17], which 1 Lx Ly zj+1 need not be repeated here. The main difference of the u = dxdydzn(~r)exp(iK~ ·~r), i i present simulations is the way in which the MD results LxLy∆z Z0 Z0 Zzj were used to calculate density wave profiles. In Ref. (14) 5 TABLE II: Comparison of interfacial free energies for different crystal faces (in erg/cm2) and anisotropy parameter ǫ4 ≡ (γ100−γ110)/(γ100+γ110) predicted by Ginzburg-Landau theory with input parameters from MD simulations for Fewith the EAM potential of MH(SA)2 [34] (Table III), and obtained from MD with theMH(SA)2 potential and two other potentials. 100 110 111 ǫ4(%) MD (ABCH) (Ref. [16]) 207.3 (10.1) 205.7 (10.0) 205.0 (10.0) 0.4(0.4) MD (Pair) (Ref. [16]) 222.5 (14.1) 220.2 (14.0) 220.8 (14.0) 0.5(0.5) MD (MH(SA)2) (Ref. [16]) 177.0 (10.8) 173.5 (10.6) 173.4 (10.6) 1.0(0.6) GL (present calculation) 144.26 141.35 137.57 1.02 GL (Shihet al. [29]) 144.26 145.59 137.57 −0.46 6 where zj and zj+1 correspond to sequential minima of 0.4 n(z) and ∆z ≡zj+1−zj. In addition, ui is evaluated at MGiDnz MbuHrg( SLAan)2d apuo t(epnrteisaelnt calculation) the midpoint of this interval, (zj +zj+1)/2. The order 0.3 Ginzburg Landau (Shih et al.) parametersuandvwerecomputedfor{110}crystalfaces using K~110 and K~101, respectively. The results of the present GL theory are compared u0.2 to those of Shih et al. and MD simulations in Fig. 2 and Table II. Using Eq. (12) with the c ’s given by i 0.1 Eq. (11), we obtain the correct ordering of interfacial free energiesγ100 >γ110 andaweakcapillaryanisotropy (γ100 − γ110)/(γ100 + γ110) ≈ 1%, consistent with the 0 110 120 130 resultsofMD simulationsforbcc elements [9,16, 17,19, 20], while the ansatz of equally weighted c ’s of Shih et z (Å) i al. (with the values listed in Table I) gives the reversed ordering γ100 < γ110. Note that the predictions of GL theory are to be compared to the MD results with the 0.4 MH(SA)2 potentialinTableIIsincethispotentialisused MD MH(SA)2 potential here to compute input parameters for this theory given Ginzburg Landau (present calculation) in Table III. MD results for the other potentials are 0.3 Ginzburg Landau (Shih et al.) mainly includedto illustratethe dependence ofγ andits v0.2 0.8 MD MH(SA)2 potential Ginzburg Landau (present calculation) 0.1 Ginzburg Landau (Shih et al.) 0.6 0 100 110 120 130 u0.4 z (Å) FIG. 3: Comparison of analytically calculated linearized or- 0.2 der parameter profiles u and v for (110) crystal faces near theliquidfromthepresentGLtheory(solidline)andtheGL theory of Shih et al. [29] (dashed line), and computed from 0 100 110 120 130 140 150 MDsimulationsusingEq. (14)withK~101 andK~110 foruand z (Å) v, respectively (solid circles). anisotropy on details of interatomic forces. 0.8 Fig. 1 shows that the planar density profile predicted MD MH(SA)2 potential Ginzburg Landau (present calculation) by GL theory, obtained by substituting Eq. (1) with Ginzburg Landau (Shih et al.) the numerically calculated order parameter profiles into 0.6 Eq. (13), is in remarkably good agreement with MD simulations on the liquid side. The discrepancy on the v0.4 solid side is due to the fact that GL theory neglects the contribution of larger |K~| reciprocal lattice vectors that contributetothelocalizationofdensitypeaksaroundbcc 0.2 lattice positions in MD simulation. Fig. 2 shows that the amplitude profiles and the in- 0 terfacewidthspredictedbyGLtheoryareingoodagree- 100 110 120 130 140 150 ment with MD simulations. The MD results clearly val- z (Å) idate the directional dependence of the rate of spatial decay of density waves in the liquid that is the main de- FIG. 2: Comparison of numerically calculated nonlinear or- terminantoftheanisotropyofγ. Thisdirectionaldepen- derparameterprofilesuandvfor(110)crystalfacesobtained dence is most clearly seen by examining the amplitude from thepresentGL theory(solidline) andtheGL theoryof profileson the liquidside ofthe interface. In this region, Shih et al. [29] (dashed line) and computed form MD sim- theamplitudesofdensitywavesaresufficientlysmallthat ulations using Eq. (14) with K~101 and K~110 for u and v, one can neglect the cubic and quartic terms in the GL respectively (solid circles). free energy functional. The resulting linear second order 7 differentialequationsforuandv obtainedbyminimizing III, Eq. (15) yields a value of γ =147.4 erg/cm2 in rea- this functional can be solved analytically, and have ex- sonably good agreement with the averagevalues of γ for ponentially decaying solutions that are compared to the the different crystal faces in Table II obtained from MD MD results in Fig. 3. The coefficients of the gradient simulationsandthefullyanisotropicGLcalculationwith squaretermsinthefreeenergyfunctionalcontrolthede- different order parameter profiles. With the same coef- cay rates. The u and v profiles (i.e. the amplitude of ficients, Eq. (16) yields a latent heat value L = 0.114 density waves corresponding to K~101 and K~110) calcu- eV/atom about 30% lower than the MD value in Table lated with coefficients that depend on the angle between III, where the difference can be attributed to the con- the principal reciprocal lattice vectors and the interface tribution of larger K~ modes that are neglected in GL normalthrough Eq. (11), which is consistent with DFT, theory. Eq. (17) in turn predicts a value of the Turnbull have different decay rates that are in good quantitative coefficient α = 0.45 that is about 25% larger than the agreement with the MD results. In contrast, u and v MD value α≈0.36,owing to the underestimation of the profiles calculated based on the ansatz of equal weights latent heat of melting in GL theory with input param- for the c ’s [29] have the same spatial decay rate, which eters of Table III from the present MD simulations. In i does not agree with the MD results. thefuture,itwouldbeinterestingtotesthowtheTurbull It is interesting to note that Mikheev and Chernov coefficient predicted by GL theory (Eq. 17) varies with (MC) [36, 37], in a formulation of the anisotropy of the input parameters computed from MD simulations using solid-liquidinterfacemobility,alsostresstheimportance different interatomic potentials. of the decayrate of the amplitude of density waves. The MCmodelpredictscrystalgrowthratesandanisotropies that are in qualitative agreement with MD simulations IV. CONCLUSIONS of FCC systems. The theory, however, is linear in the sensethatonlytheeffectivewidthsofthedensityprofiles, We have revisited the simplest GL theory of the bcc- whichareallowedtovarywithK~ andnˆ,arerequiredand liquid interface whose order parameters are the ampli- the authors make no attempt to compute, as was done tudes of density waves corresponding to principle recip- here, the full amplitude profile as a function of z. rocallattice vectors. We find that, despite its simplicity, Finally, even though we have focused primarily in this this theory is able to predict the density wave structure paperoncrystallineanisotropy,itisusefultore-examine of the interface and the anisotropy of the interfacial en- the prediction of GL theory for the magnitude of γ and ergy,inreasonablygoodquantitativeagreementwiththe for the Turnbull coefficient using input parameters from results of MD simulations. the present MD simulations. Shih et al. [29] derived Amaindeterminantoftheanisotropyoftheinterfacial an analytical expression for the magnitude of γ in the energy in this theory is the rate of spatial decay of den- isotropic approximation where all the order parameters sitywavesinthe liquid. Thisdecayratemustdependon are assumed to have the same profile through the inter- theanglebetweenprincipalreciprocallatticevectorsand face, i.e. u (z) = u for all i. In this approximation, thedirectionnormaltotheinterfaceforthistheorytobe i the free energy density reduces to the sum of the gradi- consistentwithDFT.Thisdirectionaldependence,which entsquaretermb|du/dz|2 andaquarticpolynomialinu. we validated quantitatively by MD simulations, is a di- The stationary profile u(z) that minimizes the free en- rectreflectionoftheunderlyingcrystalstructure. There- ergy is then an exact hyperbolic tangent profile and the fore,thepresentresultsprovideasimplephysicalpicture analytical expression for the interfacial energy is of the strong relationship between crystal structure and crystalline anisotropy, consistent with the findings of a γ = n0k6BTmu2s(a2b)1/2. (15) growing body of MD-based and experimental studies of crystalline anisotropy. Furthermore, Shih et al. related the latent heat (per An interesting future prospect is to extend the GL atom) to the temperature variation of the inverse of the theory to other crystal structures, and in particular fcc- peak of structure factor proportional to a2 (Eq. 9), liquid interfaces. This requires, however, to consider the coupling of density waves corresponding to the principal L= Tm ∂∆F = kBTm2 u2 da2 (16) reciprocallatticevectorstolargerK~ modes,whichmakes N ∂T (cid:12)T=Tm 2 s dT (cid:12)T=Tm the theory intrinsically more complicated. (cid:12) (cid:12) where N is the nu(cid:12)mber of atoms in the(cid:12)system. This (cid:12) (cid:12) yields the expression for the Turbull coefficient Acknowledgments α= γn0−2/3 = n10/3(a2b)1/2 , (17) This research is supported by U.S. DOE through L 3Tm da2/dT|T=Tm Grants No. DE-FG02-92ER45471 (KW and AK) and which Shih et al. evaluate using parameters for the hard No. DE-FG02-01ER45910 (JJH and MA) as well as sphere system [29]. Using the values of the various co- theDOEComputationalMaterialsScienceNetworkpro- efficients obtained from MD simulations listed in table gram. Sandia is a multiprogram laboratory operated by 8 TABLEIII:Valuesof inputcoefficients forGinzburg-Landautheorycomputedfrom MDsimulations usingtheEAMpotential for Fe from MH(SA)2 [34] and average value of γ and latent heat of melting from these simulations. a2 b (˚A2) da2/dT (K−1) us γ (erg/cm2) L (eV/atom) MD (MH(SA)2) 3.99 20.81 0.00163 0.72 175(11) 0.162 9 Sandia Corporation, a Lockheed Martin Company, for der contract DE-AC04-94AL85000. theDOE’sNationalNuclearSecurityAdministrationun- [1] J.S.Langer,inChanceandMatter,LecturesontheThe- [19] X. Feng and B. B. Laird, J. Chem. Phys. (submitted). ory of Pattern Formation, Les Houches, Session XLVI, [20] J.J.Hoyt,M.AstaandD.Y.Sun,Phil.Mag.(inpress). edited by J. Souletie, J. Vannimenus, and R. Stora [21] Y.Mu, A.Houkand x.Y.Song,J. Phys.Chem. B 109, (North-Holland,Amsterdam, 1987), pp. 629-711. 6500 (2005). [2] D.Kessler,J.Koplik,andH.Levine,Adv.Phys.37,255 [22] R.L.Davidchack,J.R.MorrisandB.B.Laird,J.Chem. (1988). Phys. (submitted). [3] M. Ben Amar and E. Brener, Phys. Rev. Lett. 71, 589 [23] R. E. Napolitano, S. Liu and R. Trivedi, Interf. Sci. 10, (1993). 217 (2002). [4] A. Karma and W. J. Rappel, Phys. Rev. Lett. 77, 4050 [24] S.Liu,R.E.NapolitanoandE.Trivedi,ActaMater. 49 (1996); Phys.Rev. E57, 4323 (1998). 4271 (2001). [5] N. Provatas, N. Goldenfeld, and J. Dantzig, Phys. Rev. [25] R. E. Napolitano and S. Liu, Phys. Rev. B 70, 214103 Lett.80, 3308 (1998). (2004). [6] J. Q. Broughton and G. H. Gilmer, J. Chem. Phys. 84, [26] S.-C. Huang and M. E. Glicksman, Acta Metall 29, 701 5759 (1986). (1981). [7] R. L. Davidchack and B. B. Laird, Phys. Rev. Lett. 85, [27] M. E.Glicksman and N.B. Singh,J. Cryst.Growth 98, 4751 (2000). 277 (1989). [8] R. L. Davidchack and B. B. Laird, J. Chem. Phys. 118, [28] M. Muschol, D. Liu and H. Z. Cummins, Phys. Rev. A 7651 (2003). 46, 1038 (1992). [9] R. L. Davidchack and B. B. Laird, Phys. Rev. Lett. 94, [29] W. H. Shih, Z. Q. Wang, X. C. Zeng and D. Stroud, 086102 (2005). Phys. Rev.A 35, 2611 (1987). [10] J.J.Hoyt,M. AstaandA.Karma,Phys.Rev.Lett. 86, [30] T. V.Ramakrishnan and M. Yussouff,Phys.Rev.B 19, 5530 (2001). 2775 (1979). [11] J.J.HoytandM.Asta,Phys.Rev.B65,214106(2002). [31] A.D.J.HaymetandD.Oxtoby,J.Chem.Phys.74,2559 [12] J. J. Hoyt, M. Asta and A. Karma, Mat. Sci. Engin. R (1981). 41, 121 (2003). [32] D.W.OxtobyandA.D.J.Haymet,J.Chem.Phys.76, [13] M. Asta, J. J. Hoyt and A. Karma, Phys. Rev. B 66, 6262 (1982). 100101(R) (2002). [33] D. Turnbull,J. Appl.Phys. 24, 1022 (1950). [14] J. R. Morris, Phys.Rev. B 66, 144104 (2002). [34] M. I. Mendelev, S. Han, D. J. Srolovitz, G. J. Ackland, [15] J. R. Morris and X. Y. Song, J. Chem. Phys. 119, 3920 D. Y. Sunand M. Asta, Philos. Mag. 83, 3977 (2003). (2003). [35] R.L.DavidchackandB.B.Laird,J.Chem.Phys., 108, [16] D. Y. Sun, M. Asta, J. J. Hoyt, M. I. Mendelev and D. 9452 (1998). J. Srolovitz, Phys.Rev.B 69, 020102(R) (2004). [36] L.V.MikheevandA.A.Chernov,Sov.Phys.JETP 65, [17] D. Y. Sun, M. Asta and J. J. Hoyt, Phys. Rev. B 69, 971 (1987). 174103 (2004). [37] L.V.MikheevandA.A.Chernov,J.Cryst.Growth112, [18] D.Y.Sun,M.I.Mendelev,C.Becker,M.Asta,K.Kudin, 591 (1991). D.J.Srolovitz, J. J.Hoyt,T. Haxhimaliand A.Karma, Phys.Rev.B (submitted).