Measuring kinetic coefficients by molecular dynamics simulation of zone melting ∗ Franck Celestini and Jean-Marc Debierre Laboratoire Mat´eriaux et Micro´electronique de Provence, Universit´e d’Aix-Marseille III and CNRS, Facult´e des Sciences et Techniques de Saint-J´erˆome, Case 151, 13397 Marseille Cedex 20, FRANCE 2 Moleculardynamicssimulationsareperformedtomeasurethekineticcoefficientatthesolid-liquid 0 interfaceinpuregold. Resultsareobtainedforthe(111),(100)and(110)orientations. BothAu(100) 0 and Au(110) are in reasonable agreement with the law proposed for collision-limited growth. For 2 Au(111), stackingfault domainsform, as firstreported byBurke,Broughton and Gilmer[J. Chem. n Phys. 89,1030(1988)]. Theconsequenceonthekineticsofthisinterfaceisdramatic: themeasured a kinetic coefficient is three times smaller than that predicted by collision-limited growth. Finally, J crystallization and melting are found to be always asymmetrical but here again the effect is much 1 more pronounced for the(111) orientation. 2 PACS numbers: 81.30.Fb, 68.45.-v, 02.70.Ns ] i c s - rl I. INTRODUCTION (EAM) [5], the glue model (GM) [4] and the effective t medium theory (EMT) [6]. In the near future, the in- m Solidification of pure elements is of technological in- crease of computer power should open the possibility to . t terest because the way a given material solidifies usu- addressthecaseofmorecomplicatedmaterialslikesemi- a m ally affects its structure and, as a consequence, its final conductors, molecular crystals and organic compounds, elastic and other macroscopic properties. From a fun- for which potentials do not simply reduce to pair inter- - d damental point of view, interest in free and directed so- actions. Very recently, the functions γℓmn and µℓmn(Ti) n lidication comes from the underlying nonlinear physics, have been determined and used in phase field simula- o morphological instabilities being at the origin of generic tions of dendritic growth for pure nickel [7]. The good c [ microstructures such as dendrites or cells. quantitative agreement found between experiments and Important theoretical and numerical contributions simulationsispromisingandshouldstimulateinthenear 1 have been made to solve this difficult physical problem futuretheconstructionofothermaterial-dedicatedphase v 1 [1]. Recently,a quantitative phasefield modelwasintro- field models. 7 duced[2]. Asubsequentrefinement,consistinginsolving New methods for the determination of the func- 3 thediffusionequationwiththehelpofBrownianwalkers, tions γℓmn and kℓmn(Vi) have been recently proposed 1 permitted to bridge the wide gap between the capillary [8,9]. In the present paper, we rather concentrate on 0 and diffusion lengths, allowing direct comparison with µℓmn(Ti). The kinetic response of a solid-liquid inter- 2 experiments [3]. As a consequence, there is currently face has been simulated quantitatively in the 80’s by 0 / an increasing need for accurate values of the interfac re- Broughton, Gilmer and Jackson (BGJ) for a Lennard- t a sponse functions that are used as input parameters for Jones (LJ) potential and a (100) orientation [10]. These m realistic phase field simulations. authors showed that growth is not diffusion-limited but - Inthe caseofapureelement, thesurfacetensionγℓmn rather that the interface velocity is related to the mean d must be known as a function of the interface orientation kinetic energy of the atoms. For this collision-limited n (ℓmn). In addition, the kinetic coefficient µ (T ) giv- growth regime, the growth rate should be directely pro- o ℓmn i c ing the relation between the interface velocity and the portional to the distance between two successive layers v: interface temperature Ti, should also be known for the dℓmn. Indeed, since the liquid atoms do not diffuse to i differentorientations. Forabinaryalloy,temperaturede- choose their adsorption sites but are almost instanta- X pendenceofthesolutediffusioncoefficient,D(T),aswell neously incorporated into the solid, the larger dℓmn the r asvelocityandorientationdependenceofthesegregation more effective and the faster the advance of the solid- a coefficient k (V ) are also necessary. liquid interface should be. The analytical expression for ℓmn i Both k and µ are hardly accessible in the experiments thegrowthvelocityofaroughsolid-liquidinterfacereads and convection effects often lead to overestimated val- [11]: uesofdiffusioncoefficients. Differentsimulationschemes ∆µ˜ have thus been proposed as an alternative. Such nu- V ∝dℓmnV˜h1−exp(cid:0)− kT (cid:1)i, (1) merical experiments have been rendered possible by the discovery of realistic interatomic potential models, such dℓmn being the interplane spacing, ∆µ˜ the chemical po- as, in the case of metals, the embedded atom model tential difference between solid and liquid phases, T the absolute temperature, k the Boltzmann constant, and V˜ 1 the thermal velocity. This law is confirmed by molecu- lar dynamics simulations for the (100) and (110) orien- tations: the expected √2 ratio between the correspond- ingkineticcoefficientsiswellrecoveredforseveralmetals cristalizing in a face centered cubic (fcc) structure (Ni, Ag and Au) [12,13]. Nevertheless, for these rough ma- terials, growth of the (111) interface does not obey this simple law: according to Eq. (1), the (111) orientation shouldbe muchfaster,andwhatis foundis preciselythe opposite. Burke, Broughton, and Gilmer [11] attribute this slowing-down to the growth of competing fcc and hcpdomainsinthesolidfyinglayer,followedbytheelim- ination of the defect lines between the two phases. Another question associated with solid-liquid inter- facesisthatofsymmetrybetweensolidificationandmelt- ing kinetics. Asymmetry has been already observed in different systems. It is not really surprising for faceted materialslikesiliciumwheresolidificationinvolvesnucle- ation while melting does not. The question is more del- icate when one considers rough materials with collision- limited growth. Indeed, available results are controver- sial: ifasymmetryhasbeenfoundforaNa(100)interface [14], it has not been observed for a LJ(100) [15]. More surprisingly, in the latter case an opposite asymmetry (crystal growing faster than the melt) can be found, de- pending on the way the solid germ is prepared. Inthis paper,weaddresstheabovequestionsconcern- ing the growth of a rough solid-liquid interface. We first presentourimplementationofanon-equilibriummolecu- FIG. 1. A typical simulation box with periodic bound- lardynamicsschemeforazonemelting experiment. The aryconditionsin allthreedirections ((111) solid-liquid inter- secondsectionis devotedtothe study of(100)and(110) faces). Atomsin dark grey are within thehot and cold slices orientations. The special case of (111) growth is exam- where temperature is fixed. ined in section III and asymmetry between melting and solidification in section IV. Finally, a summary of the different results and a discussion are given in the last U is the glue term associating an extra potential energy section. to atom i as a function of its coordination. This glue potentialhasdemonstateditsefficiencyinpredictingthe physical properties of gold as well as in describing sev- II. SIMULATION PROCEDURE eral experimentally observed phenomena such as surface melting and surface reconstructions [16]. For this study we use the Ercolessi glue potential for A distinctive feature of our method is to simulate a Au [4]. In this formalismthe total potentialenergy for a zone melting experiment in which both a solidification system of N atoms is given by: and a melting front are simultaneously advancing at a fixed velocity V. This velocity is that of the virtual fur- 1 N N nace which imposes two symmetric thermal gradients. U = 2 X Φ(rij)+XU(ni) (2) The particle coordinatesare defined in a reference frame i,j=1 i=1 moving at velocity V in the z direction, so that after equilibrationthe positions of the two interfaces are fixed Thefirsttermisaclassicalpairinteraction. Inthesecond in the simulation box. Heat transport from the furnace term, n is the coordination of atom i, i issimulatedbyimposingonetemperaturebelowandone N above the melting point inside two distant slices, 20˚A n = ρ(r ), (3) each in thickness (Fig. 1). Within each slice, tempera- i X ij j=1 tureiskeptconstantbyusingaclassicalvelocityrescaling procedure[17]. Periodicboundaryconditionsareapplied whereρ(rij)is afunctionofthe interatomicdistancerij, inthethreedirections. Moredetailsaboutthenumerical with a cut-off radius of 3.9˚A here. The energy function methodcanbefoundinarecentstudyofsolutetrapping 2 −3.4 1400 K) s) V=0 m ( o V=15 m/s T t a 1000 / V −3.5 e ( Y G R ) −3.4 E s N m E o −3.6 t a / −3.5 V e 1000 1200 1400 ( E T (K) −3.6 −80 −40 0 40 80 o FIG. 3. Caloric curves for V = 0 (circles) and V = 15 Z (A) ms−1 (diamonds). For non zero velocity, the kinetic effects split the curve in two parts. The dotted and full straight FIG.2. Temperature and energy profiles along the z axis lines represent respectively the functions Ei(T), EL(T) and perpendicular to theinterfaces. ES(T). inaLJbinaryalloy,whereasimilarsimulationtechnique run lasts 106 MD steps (3.5 103 ps), so that the atoms × was used [9]. in the system solidify and melt severaltimes. According First, the fcc solid and the liquid are equilibrated to a recent study by Tepper et al. [15], we know that separatly at zero pressure and at a temperature close the melting kinetics can be affected by the way the solid of the melting point. Our smallest system has a size is equilibrated. Thusmulti-cycling is necessaryto mimic S0 20 20˚A2 in cross-section, that is about 64 atoms themeltingofarealsolid,usuallyresultingfromprevious ≃ × per layer. After equilibration, the solid and liquid are solidification(s). broughtintocontactandplungedinthetemperaturegra- To conclude this section, the method used to estimate dient imposed by the two temperature-controled slices. the interface temperature T from the caloric curves is i The total system is about 220˚A in height. After a sec- described. We assume the energy of atoms lying at the ond equilibration period (during which the velocity of interace,E ,tobeaweightedaverageoftheperfectsolid i the furnace is zero), the two interfaces reach a station- and liquid energies at the same temperature T. arypositionandweroughlyhave50%ofsolidandliquid (see Fig. 1). Fig. 2 shows the temperature and energy Ei(T)=αES(Ti)+(1 α)EL(Ti). (4) − profiles along the z axis. Combining the two profiles to eliminate the z coordi- Linearrelations,ES(T)=aST+bS,andEL(T)=aLT+ nate, oneobtains a caloriccurve,i.e., a plotofenergyas bL, are fitted to the data points obtained onthe low and a function of temperature. In Fig. 3, the caloric curves high temperature side, respectively (Fig.3). The curve obtainedfortwodifferentvaluesofthepullingvelocityV Ei(T) is thus a line with a slope are displayed. For V =0, the data points corresponding p=αa +(1 α)a . (5) to the solidification and the melting fronts merge onto S − L the same curve: no kinetic effects are at play and the The value of coefficient α is then extracted from the interfacetemperatureisthe equilibriummeltingtemper- caloriccurveat zerovelocity,for which T must be equal ature T0 1330K. When a velocity is imposed, a dy- i ≃ toT0 (Fig. 3). Finally,theinterfacetemperatureisgiven namical hysteresis appears on the caloric curve. Kinetic by the intersection of the line E (T) with the caloric effects split the curve in two distinct parts: the interface i curve. An alternative method consists in building an or- temperature of the solidification front decreases while it derparameterthatdistinguishesbetweensolidandliquid increases on the melting front. We can deduce both in- atoms[8,18]: aplotofthisorderparameterasafunction terfaceundercoolingandinterfacesuperheatingfromthis of temperature also gives an interface temperature. We plot. An interest of this method is that, as said before, have checked that the two methods give equivalent re- the interface is fixed in the reference frame of the simu- sults. lationbox,sothatstatisticsareeasytorecord. Atipical 3 30 1.0 (100) s) (100) 0.8 / m Y ( 20 (110) N 0.6 T µ I C O 0.4 L 10 E (111) V 0.2 V=15 m/s 0.0 0 0 50 100 1 10 100 T − T (K) S/So o i FIG.4. Velocity of thesolid-liquid interfaceas afunction FIG. 5. Normalized kinetic coefficient as a function of ofundercoolingfor(100)and(110)orientations. Thestraight system size S for the (100) and (111) orientations. lines are best fits to a linear kinetic law. and we obtain roughly the same behavior for the (110) III. GROWTH OF (100) AND (110) INTERFACES interface. This size effect has never been reported in the past for (100) and (110)orientations: the fact that Hoyt In this section we compute the kinetic coeficient for andco-workersdonofindsizeeffectsforthesetwoorien- the Au(100) and Au(110) interfaces, using the method tations [8] is certainly due to the fact that their smaller described above. We concentrate here on pulling veloc- system is larger than ours. We can now propose extrap- ities ranging between V = 5 ms−1 and V = 30 ms−1, olated values for the kinetic coefficients: for which kinetics remain linear. We also perform a −1 −1 few simulations at higher velocities, where kinetics de- µ100 =18.8 1.0cms K (9) ± viates from linearity, but comments on nonlinear effects are postponed to the concluding section. In Fig. 4, we µ110 =12.6 1.0cms−1K−1 (10) ± plot the interface velocity as a function of the measured undercooling T0−Ti. Linear fits to the law Tgohoedcoargrreesepmoenndtinwgitrhatitoheµ1v0a0/luµe11√0 2=p1r.e4d9ic±ted0.1b5yisEiqn. V =µℓmn(T0 Ti) (6) (1). Hence, the assumption of collision-limited growth − for (100) and (110) orientations is confirmed to be the give the following estimates for the two kinetic coeffi- relevant one. At this point, we can compare our results cients: with those of Hoyt et al. for gold [13]. If they also find a √2 ratio between their two orientations, their µ val- µ⋆ =23.1 1.0cms−1K−1 (7) 100 ues are larger than ours by a factor 1.8. Linearizing the ± expression given by BGJ, we find and µ⋆ =15.5 1.0cms−1K−1. (8) V ∼T0−1Ti−1/2(T0−Ti) (11) 110 ± for the interface velocity. The potential used by Hoyt However, finite-size effects are expected to bias these et al. gives a melting point T0 of 1090K [19] much estimates because the system cross-section area, S0 = smallerthanthevalue1330KobtainedwithEroclessipo- 20 20˚A2, is rather small. tential. Introducing this temperature shift in Eq. (11) × Additionnalrunsarethusperformedinordertoquan- roughly accounts for the discrepancy between the values titativelyevaluatefinite-sizeeffects. Thepullingvelocity ofµ. SinceErcolessipotentialgivesameltingpointmuch is fixed to V =15ms−1, the system heightto H 222˚A, closerto the experimentalone,it shouldbe alsothe case ≃ and the cross-section area S is progressively increased. for our estimates of the kinetic coefficients. WedefinethenormalizedkineticcoefficientµN(S)asthe In order to understand the origin of the size effects on ratio of the kinetic coefficient obtained at size S to that the value of the kinetic coefficient, we take now a closer obtained at size S =S0 [Eqs (7-8)]. As shown in Fig. 5, look at the in-plane structure of gold layers in the vicin- the size effects are important and the kinetic coefficients ity of the solid-liquid interface. We compute a density appear to converge only for S 100 100˚A2. For the profile along the z axis from which we are able to sepa- ≃ × (100) direction, there is a decrease of about 20 percent rateatomsbelongingtodifferentlayers. Deepinthesolid 4 the in-plane square structure of the (100) orientation is effectively recovered without any significant amount of defaults and vacancies. For the two solid layers just be- low the interface the situation is more complex. To dis- tinguish between different symmetries, we first perform aVoronoiconstructionforallthe atomsinthe layer. We then collect the set of first neighbors for each atom. For a fcc solid with lattice parameter a, on the square lattice of the (100) orientation an atom has four nearest neighbors at a distance a/√2 and four second nearest neighbors at a distance a. On the other hand, for a tri- angular lattice (as the one of the (111) plane) the six neighbors all lie at the same distance a/√2. In Fig. 6, we showasnapshotofthe interfacesolidlayerwherethe Delaunay triangulationis only drawnfor the atoms that FIG. 6. Snapshot showing the atoms in the (100) solid have six first neighbors at comparable distances, in or- layer next to the interface. The Delaunay triangulation is der to reveal the local triangular structure. It is clear only drawn in regions with triangular underlyingsymmetry. that most of the atoms have reached their positions on the square lattice but several islands with a triangular 30 symmetry remain. Note that the number of atoms in the layerhasalreadyattainedthe value itwillhavedeep in the solid with a perfect square structure. To compen- 20 sateforthehigherdensityofthetriangularstructure,the ) s correspondingislands are surrounded by a border region m/ where the density is very low. This coexistence of two V ( symmetries is not observed in our smallest system: one 10 can imagine that for a small area the square structure is easily formed and hence triangular islands do not ap- pear. This phenomenon is very close to the well known 0 reconstruction of the (100) solid-vapor interface where 0 50 100 150 thefirstlayeradoptsatriangularstructure[21]. Turning T − T (K) back to the solid-liquid interface, the system apparently b i uses some of the solidification driving force to eliminate FIG. 7. Velocity of the (111) solid-liquid interface as a oneofthe twophasesandfinallyreachanalmostperfect function of undercooling. square symmetry. Hence, the interface velocity is lower for larger systems. orientation. The extrapolated value of the kinetic coeffi- Such an in-plane ordering is not taken into account cient, in the collision-limited model but in spite of this we re- coverthepredicted√2valuefortheratioµ100/µ110. This µ111 =7.0 1.0cms−1K−1, (12) suggests a similar effect, roughly of the same order, for ± the (110) orientation. We have not been able to visual- is now 60 percent below its value for S =S0. Relatively ize ordering at (110) interfaces but one could imagine a to the two other orientations, we find mechanismreminiscentofthemissingrowreconstruction observed for (110) solid-vapor interfaces. µ111 0.37µ100 0.56µ110. (13) ≃ ≃ These ratios largely differ from the values predicted by IV. THE SPECIAL CASE OF (111) INTERFACE Eq. (1), respectively 2/√3 1.15 and 2 2/3 1.63. ≃ p ≃ The (111) orientation, expected to grow faster because We now turn to the case of the (111) orientation. In of a larger interlayer spacing, is surprisingly found to be the same way as above for the (100) and (110) orienta- theslowestone. Thisdiscrepencytellsusthatthegrowth tions we calculate the interface temperature for different mechanismforthe(111)orientationisnot,oratleastnot velocities. As can be seen in Fig. 7, a linear kinetic only, a collision-limited one. law is also valid for the (111) orientation. Results of the Here again we look at the symmetries inside the lay- finite-size analysis, presented in Fig. 5, show that the ers close to the interface. For a (111) layer, there are sizeeffectsaremuchmorepronouncedthanforthe(100) threepossibleorderedphaseslyingonthreedifferentbut equivalentsub-latticesthatwewillcalla,bandc. Asthe 5 20 With our simulation scheme, this study is straightfor- ward,since both a melting anda solidificationfronts are simulated at once: no additional calculations are thus s) 0 required. Fig. 9 represents the velocities of both the m/ V ( meltingandsolidificationfrontsasfunctionsofT0−Tifor the(111)orientation(inourconventionsapositiveveloc- −20 ity corresponds to solidifation). The data are obtained in a system of size S = S0 and corrected according to thefinite-size analysisreportedabove. Itisimportantto −100 0 100 200 notethatnosizeeffectsareactuallyfoundforthemelting T − T (K) o i front: in contrast with the solidification front, the melt- inginterfacetemperatureremainsthesamewhateverthe FIG. 8. Snapshots of the three solid layers immediately belowtheinteface(toplayerincontactwiththeliquidphase). system size. This can be understood if one remembers Thegrey levelscorrespond tothethreedifferentsub-lattices: thatforsolidification,especiallyforthe(111)orientation, a (white), b (light grey), c (dark grey). growth is not only collision-limited but also requires in- plane ordering. This is no longer the case for melting, which justifies the absence of size effects. The asymme- stacking fault energy is weak for gold (it is actually zero try shown in Fig. 9 is larger for the (111) orientation. for the potential we use), once a perfect, say a, layer is The same analysis is also made for the two other orien- formed, the next layer to form is either b or c. In Fig. 8, tations and we find the following degrees of asymmetry: we show a snapshot of the three uppermost solid layers and we distinguish between atoms belonging to a, b or µm111 =25 4cms−1K−1 3.6µs111 (14) ± ≃ c phases. For the lowest solid layer, phase a is selected anditoccupiesthewholeplane. Forthelayerjustabove, µm100 =39 2cms−1K−1 2.1µs100 (15) ∗± ≃ thereiscoexistencebetweenbandcsub-lattices. Finally, inthehighestsolidlayer,allthreephasescoexist. Were- µm110 =20 2cms−1K−1 1.6µs110 (16) ± ≃ cover here the effect first observed by Broughton et al. wherethesuperscriptssandmreferrespectivlytosolid- [11]foraLJpotential. Fora(111)orientationthesystem ificationandmelting kinetics. An asymmetryis revealed hesitatesbetweenthe differentphasesitcanequivalently inthethreecasesbutitismorepronouncedforthe(111) form. Here again, the system dissipates a part of the orientation in the same way as size effects observed dur- available driving force to select one of the phases. As a ingsolidification. Weconcludeherethatthisasymmetry consequence, the velocity of the interface is reduced as is directly related to the ordering within the interface compared to the value expected for a purely collision- layers. The asymmetry is strong for (111) because of limited growth. The size effect is easily understood be- thepeculiargrowthmechanismdiscussedintheprevious cause in a small system, coexistence is strongly reduced. section. It would be of interest to determine the amount of driv- The melting front is found to be faster than the solid- ing force spent in this in-plane organisation in order to ification interface in agreement with the idea that disor- estimate the correspondingdecrease inV111. To perform dering is an easier task than ordering. Our results con- this, one could for instance use a 3-state Potts model firm the majority of experimental and numerical studies in three dimensions with ferromagnetic intra-plane and [14,22–25]. We also confirm the conclusions of a debate anti-ferromagnetic inter-plane interactions. To conclude betweenRichards[26]andOxtobyandco-workers[27,28] this section we have to point out that phase coexistence on the importance of density change on the asymmetry is related to the value of the stacking fault energy E . s betweenmeltingandsolidificationkinetics. Inagreement ForamaterialwithlargeE phasecoexistenceshouldbe s with the conclusions of Oxtoby, the gold density change less probable and the front velocity in better agreement at melting is small ( 2%) and can not be responsible with the prediction of Eq. (1). ≃ for such an important asymmetry. On another hand, Tepper [15] does not find asymmetry for the growth of a (100) LJ solid. Even if the materials differ, they both V. ASYMMETRY BETWEEN MELTING AND belong to the same class of rough materials and such a SOLIDIFICATION qualitative difference may be surprising. Nevertheless, one should remember the strong tendency to surface re- Asdiscussedintheintroduction,asymmetryisobvious constructionin Au, as observedfor the (100)orientation for faceted materials but is not as clear when consider- where triangular-like regions are formed. This tendency ing roughmaterials like metals. The question is to know isfuthermoreenhancedbytheuseofErcolessigluepoten- if, at equal absolute undercooling and superheating, the tial but is weaker for a LJ potential, what could explain solid-liquidandliquid-solidfrontshavethesamevelocity. the different behaviors observed. 6 20 extractthe amountofdrivingforcespentforphasesepa- rationin order to modify Eq. (1) and find an acceptable expression for the interface velocity. For melting we find the following order between the s) 0 differentkineticcoefficients: µm100 >µm111 >µm110. Toour m/ knoweldgethis hierarchydoes notobey anyexistinglaw. V ( Thisresultwillhopefullystimulatefurtherinvestigations to reach a clear understanding of the specifities of melt growth as compared to crystal growth. −20 The presentstudy is devotedto the linearrelationship between velocity and undercooling. For all the orienta- −100 0 100 200 tionsconsideredhere,nonlineareffectsappearatvelocity V 30 ms−1 and undercooling ∆T 200K. It is not T − T (K) ≃ ≃ o i possible to explain this deviation using either the diffu- sion limited [29,30] or the collision-limited growth law. FIG. 9. Velocity of the (111) solidification and melting Thissuggestsapossiblechangeintheinterfacestructure interfacesasafunctionofundercooling. Notethattheresults forsuchlargedeviationsfromequilibrium. Densitydiffer- for solidification incorporate finite-size corrections. encebetweentheliquidandsolidphasesshouldalsocon- tributetotriggernonlinearbehavior[26]. Understanding this cross-over would be of importance in the context of Finally, comparing the melting kinetic coefficients in very rapid solidification. the different orientations, we find µm > µm > µm . 100 111 110 Finally, we would like to stress that the kinetic effects We presently do not have a satisfactory explanation for can contribute to the anisotropy of the segregationcoef- this hierarchy in the melting kinetics. ficientk(V)forabinaryalloy. Atsufficientlylargeveloc- ity, one expects an important difference in the interface temperaturesfor (111)and(100)orientations. As a con- VI. DISCUSSION sequence, the diffusivity of solvent atoms and hence the segregationcoefficient, as predictedby the Aziz law [31], Our molecular dynamics simulations of zone melting should also differ. This effect may cause solute trapping experiments allow us to measure simultaneously the so- to appear at lower velocities for (111) than for (100) or lidification and melting kinetics for a pure element. (110) orientations. We are currently investigating such For (100)and(110)orientations,growthis apparently segregation effects induced by kinetic anisotropy in the well described by a collision-limited process. Neverthe- Al-Cu system. less, we observe small 2D islands with triangular sym- metry to form in the solid layer at the (100) solid-liquid interface. As a consequence, size effects and asymmetry ACKNOWLEDGMENTS betweenmeltingandsolidificationarefound. Wecannot decide whether this effect is solely due to the tendency Itisapleasureforus tothank M.Asta, B.Billia,J.J. of the glue potential to overestimate surface reconstruc- Hoyt and A. Karma for fruitfull discussions. tion, or if it is an intrinsic property of gold and/or other Electronic address: [email protected] metals. ∗ The case of the (111) orientation is rather special. Phasecoexistenceofthreetriangularsub-lattices,asfirst proposedby Broughtonet al. [11],is recovered. This pe- culiar behavior has a strong influence on the kinetics of the interface. Our finite-size analysis show that in order to measure a realistic value of the kinetic coefficient one [1] K. Kassner, Pattern Formation in Diffusion-Limited has tosimulate systems witha solid-liquidinterfacearea Crystal Growth (World Scientific,Singapore, 1996). largerthan100 100˚A2. Theconsequenceonasymmetry [2] A. Karma and W.J. Rappel, Phys. Rev. E 57, 4323 × between melting and solidification is also of importance: (1998). for a given driving force, the melting front is more than [3] M. Plapp and A. Karma, Phys. Rev. Lett. 84, 1740 three times faster than the solidificationone. Because of (2000). [4] F. Ercolessi, M.Parrinello, and E.Tossati, Phil.Mag. A this disagreementwith a purely collision-limited growth, 58, 213 (1988). no analytical model seems, at present, able to predict [5] M.S. Daw, and M.I. Baskes, Phys. Rev. B 29, 6443 the kinetic law of a (111) interface. As discussed previ- (1984). ously,itwouldbeinterestingtouseastatisticalmodelto 7 [6] K.W. Jacobsen, J.K. Nørskov, and M.J. Puska, Phys. Rev.B 35, 7423 (1987). [7] J. Bragard, A. Karma, Y. H. Lee and M. Plapp, private communication (2001). [8] J.J.Hoyt,M. AstaandA.Karma,Phys.Rev.Lett. 86, 5530 (2001). [9] F. Celestini and J.M. Debierre, Phys. Rev.B. 62, 14006 (2000). [10] J.Q.Broughton,G.H.GilmerandK.A.Jackson,Phys. Rev.Lett. 49, 1496 (1982). [11] E. Burke, J. Q. Broughton and G. H. Gilmer, J. Chem. Phys.89, 1030 (1988) [12] J. J. Hoyt, B. Sadigh, M. Asta and S. M. Foiles, Acta. Mater. 47, 3181 (1999). [13] J. J. Hoyth, M. Asta and A. Karma, private communi- cation (2001). [14] J. Tymczack and J. R. Ray, Phys. Rev. Lett. 64, 1278 (1990); J. Chem. Phys.92, 7520 (1990). [15] H. L. Tepper and W. J. Briels, J. Cryst. Growth 230, 270 (2001). [16] F. Ercolessi, A. Bartolini, M. Garofalo, M. Parrinello, and E. Tosatti, Surf. Sci. 189 55 (1989); F. Ercolessi, S. Iarlori, O. Tomagnini, E. Tosatti, and X. J. Chen, Surf. Sci. 251, 951 (1991). [17] M. P.Allen and D. J. Tildesley, Computer simulation of liquids.(Claredon Press, Oxford, 1987). [18] W.J.BrielsandH.L.Tepper,Phys.Rev.Lett.79,5074 (1997). [19] S. M. Foiles and J. B. Adams, Phys. Rev. B. 40, 5909 (1989). [20] F. F. Preparata and M. L. Shamos, Computational Ge- ometry : An Intoduction (Springer, New-York,1985). [21] D. Passerone, SISSA-ISASPhD Thesis (1998), available at http://www.sissa.it/cm/document/doc cm phd.htm. [22] M. D. Kluge and J. R. Ray, Phys. Rev. B. 39, 1738 (1989). [23] J.Y.Tsao,M.J.Aziz,M.O.ThompsonandP.S.Peercy, Phys.Rev.Lett. 56, 2712 (1986). [24] R. Moss and P. Harrowell, J. Chem. Phys. 100, 7630 (1994). [25] G. H. Rodway and J. D. Hunt, J. Cryst. Growth 112, 554 (1991). [26] P.M. Richards, Phys.Rev.B 38, 2727 (1988). [27] D. W. Oxtoby and P. R. Harrowell, J. Chem. Phys. 96, 3834 (1991). [28] Y. C. Shen and D. W. Oxtoby, J. Chem. Phys. 104, (1996). [29] H.A. Wilson, Phil. Mag. 50, 238 (1900). [30] J. Frenkel,Phys. Z. Sowjetunion. 1, 498 (1932). [31] M. J. Aziz, J. Appl.Phys. 53, 1158 (1981). 8