Melting and crystallization in Ni nanoclusters : The mesoscale regime

We studied melting and freezing of Ni nanoclusters with up to 8007 atoms ~5.7 nm! using molecular dynamics with the quantum-Sutten–Chen many-body force field. We find a transition from cluster or molecular behavior below;500 atoms to a mesoscale nanocrystal regime ~well-defined bulk and surface properties ! above;750 atoms~2.7 nm!. We find that the mesoscale nanocrystals melt via surface processes, leading to Tm,N5Tm,bulk2aN , dropping fromTm,bulk51760 K to Tm,336 5980 K. Cooling from the melt leads first to supercooled clusters with icosahedral local structure. For N.400 the supercooled clusters transform to FCC grains, but smaller values of N lead to a glassy structure with substantial icosahedral character. © 2001 American Institute of Physics. @DOI: 10.1063/1.1373664 #


I. INTRODUCTION
Particles with diameters of 1 nanometer to 10 nanometer exhibit properties that are often intermediate between molecular and crystalline.The properties of such mesoscale nanoparticles are expected to evolve gradually from molecular to solid like as the particle number increases, but it has been difficult to study the mesoscale region experimentally.It is well established that in the nanometer range the melting temperature decreases dramatically with the decreasing radius of the cluster, 1 but there is little quantitative data on the structures and energetics of the bulk and surface regions important to the properties of such system.Since transition metal clusters have a number of exciting potential applications in nanoscale electronic devices and catalysis, [2][3][4][5] it is important to develop quantitative understanding of the thermodynamic and structural properties of such metal clusters.Thus, in this paper we carry out molecular dynamics ͑MD͒ simulations using an accurate many-body force field to study melting and crystallization processes for Ni nanoclusters containing up to 8007 atoms ͑5.74 nm diameter͒.We find that at the upper of this range, the properties can be described accurately in terms of macroscopic concepts with welldefined bulk plus surface thermodynamic properties.
General considerations ͑see Sec.IV͒ suggest that when the crystal-vapor surface energy (␥ xv ) exceeds the sum of the liquid-vapor (␥ Lv ) and crystal-liquid (␥ xL ) surface energies ͓⌬␥ϭ␥ xv Ϫ(␥ Lv ϩ␥ xL )Ͼ0͔, then the formation of liquid skin on a planar or curved solid surface is energetically favorable.Using experimental bulk surface energies ͑Ref.11͒, and assuming that the driving force for melting on a surface is the reduction in the total interfacial energy, these considerations suggest that Au(⌬␥ϭ33 mJ/m 2 ), Sn(⌬␥ ϭ18 mJ/m 2 ), Pb(⌬␥ϭ22 mJ/m 2 ), and Pt(⌬␥ϭ29 mJ/m 2 ) would favor surface melting.Indeed experiments on Au 1 , Sn 7,8 , Pb 9 , Pt 10 , and Na 11 small particles, show melting at considerably lower temperature than bulk systems, suggesting that the melting process starts on the surface.By the surface energy criterion bulk Na(⌬␥ϭ8 mJ/m 2 ) and Ni(⌬␥ϭϪ2 mJ/m 2 ) are not expected to exhibit surface melting on planar surface.However small particles of Na ͑Ref.11͒ also show melting at much lower temperature than bulk systems, and the results have not yet been established for Ni.Experimental studies have not yet been reported on Ni, and hence we decided to study Ni using MD simulations to determine if surface melting is a general phenomenon in small particles.Previous melting simulations have been reported for small Ni clusters with 7-23 atoms 12 and 55 atoms, 13 but no systematic study of the size dependence of melting and crystallization behaviors, has been reported for Ni.
In this paper we apply molecular dynamics ͑MD͒ methods using the quantum Sutten-Chen ͑Q-SC͒ many-body force field to simulate melting and crystallization processes for Ni nanoclusters containing up to 8007 atoms ͑5.74 nm diameter͒.In particular, we examine how the melting temperature and heat of fusion depend on cluster size, including premelting prior to the transition.Then we examine the local structures formed upon solidification and the melting temperature during reheating of the newly formed crystals.
Section II summarizes various details of the calculations, while Sec.III describes the results for melting and crystallization of Ni clusters, with comparisons to results of bulk Ni.
The implications are discussed in Sec.IV and the Conclusion in Sec.V.

A. Dynamics
We use the quantum corrected Sutten-Chen ͑Q-SC͒ many-body force field 14 with parameters empirically fitted to data on density, cohesive energy, compressibility, and phonon dispersion.Q-SC leads to accurate values for surface energies, vacancy energy, and stacking fault energies and has previously been applied to studies of melting, glass transformation, and crystallization for bulk NiCu and CuAg alloys 15 and for studying plasticity and strain rate induced amorphization in one-dimensional nanorods. 16he MD calculations 17 were performed using the MPiSim software developed at Caltech. 18The calculations were carried out using two flavors of dynamics.
͑a͒ The nanoclusters were simulated using an MD Hamiltonian with constant temperature ͑T for Hoover 17 ͒, constant shape ͑h͒, and constant particle number ͑ThN͒ without periodic boundary conditions.͑b͒ The bulk systems were studied with periodic boundary conditions using an MD Hamiltonian with constant temperature, constant stress ͑t for Rahman-Parrinello 17 ͒ and constant particle number (TtN).
The MD equations of motion were integrated using a fifthorder predictor-corrector algorithm with a time step of 1.0 fs.This led to quite stable dynamics trajectories.

B. Clusters
All clusters started with geometries constructed from a large FCC block of Ni, using various spherical cutoff radii centered at a tetrahedral interstitial site.This leads to smooth but faceted surface structures.The ThN simulations for each cluster considered it to be in a large cubic unit cell with fixed size ͑no periodic boundary conditions͒.Here the cluster adjusts to keep a zero pressure.We considered clusters ranging from 336 atoms ͑2.144 nm diameter͒ to 8007 atoms ͑5.74 nm diameter͒.
For the simulations of the bulk system, we used a cubic periodic super cell with 500 atoms.Again, we started with the perfect bulk FCC structure.

C. Dynamics procedure
After equilibrated at 300 K for 25 ps, the clusters were subjected to a heating-cooling cycle consisting of a series of ThN/TtN MD simulations with temperature increments of ⌬Tϭ100 K and equilibration simulation times of dt ϭ25 ps.However, for the temperature region near the melting point, we used smaller temperature increment, ⌬T, but for shorter times dt so that the heating rate remains constant at ⌬T/dtϭ4 K/ps ͑e.g., 20 K for 5 ps͒.

D. Analysis of thermodynamic quantities
To calculate the heat capacity, we fitted smooth cubic splines to the average potential energy during the heating process as a function of temperature and then obtained the derivative from this fitted curve, We defined the melting point as the temperature with the maximum apparent heat capacity.
To obtain ⌬H melt ͑the heat of fusion at T melt ͒, we fitted the potential energy to a linear function of T in the solid and liquid phase ͑above and below the transition temperature͒ and determined the difference at T melt .To calculate the entropy change at the transition, we used ⌬S melt ϭ⌬H melt /T melt together with the computed values for the energy of fusion.
We defined the radius of the cluster as

͑2͒
where R g is the radius of gyration and the atomic radius R Ni is half the atomic distance in bulk system Comparing the surface area calculated from ͑3͒ with the solvent accessible surface ͑using R Ni as solvent radius͒ we find a difference of less than 5%.Thus, we have used R c to calculate the internal volume and surface area as a function of temperature.

E. Analysis of structural characteristics
To analyze how the atomic motions change near the melting point, we partitioned the cluster into six radial shells of equal dR.Within each shell we calculated the average root-mean-square thermal displacement ͑RMSD͒ as in ͑4͒, where To obtain RMSD we averaged over a 25 ps trajectory and used dt as 10 ps, which is sufficiently long for diffusion in the liquid phase.This provides a rough measure of the diffusion constant, which increases monotonically with the distance from the center.The RSMD changes dramatically upon melting, indicating a first order transition.
To investigate the local geometric structure, we searched for local FCC structural regions by finding neighboring atoms with an FCC packing sequence. 16We also used the Honeycutt-Andersen 19 index to distinguish between the local icosahedral and FCC structures.This has proved to be a useful tool for other studies of local structures. 19

A. Size effects on cluster melting
The analysis is illustrated in Fig. 1, which compares the potential energy for heating and cooling of the 1004 atom cluster with the corresponding results for bulk Ni.The melting transition is clearly identified by a rapid increase in energy, which is well defined but much broader than for the bulk system.We define the melting temperature, T melt , as the maximum in the heat capacity Cp (T), leading to T melt ϭ1160 K. Figure 1 shows that the potential energy begins to deviate from linearity at 800 to 900 K, but the Cp indicates that melting starts ͑becomes nonlinear͒ at ϳ1000 K and extends to over 1200 K.This indicates a heterogeneous melting process.
Figure 1 also shows the calculated heating curve for bulk Ni.These periodic boundary condition calculations provide no free surface for bulk Ni leading to an abrupt ͑ϳ40 K width͒ homogeneous melting transition at a temperature of T melt ϭ1760 K, which is 32 K ͑1.8%͒ higher than the experimental melting temperature of T melt ϭ1728 K for pure Ni.Such a higher melting temperature might be consistent with surface melting depressing the experimental value.
The total potential energy per atom of the cluster is higher than that of the bulk system due to the surface energy.For the 1004 atom system, the difference at 300 K is 33.6 kJ/mol.Using the calculated radius of 14.97 A ͑2816 A 2 surface area͒, we obtain a surface energy of 1990 mJ/m 2 , which compares well with the experimental surface energy for bulk Ni of 2104 mJ/m 2 and the calculated value for bulk Ni ͑same FF͒ of 2202 mJ/m 2 . 14The surface energy predicted for the various clusters at 300 K and 1500 K are shown in Table I.At 1500 K, we compared the liquid clusters with supercooled liquid bulk Ni to obtain the surface energy of liquid Ni.This leads to ␥ LV ϭ1537 mJ/m 2 for cluster with 1004 atoms, which can be compared to the calculated surface energy at 1500 K for bulk Ni of 2057 mJ/m 2 and the experiment value for bulk liquid Ni of 1750 mJ/m 2 .
Table I also shows how the various properties change with cluster size.The melting temperature drops dramatically with cluster size in the nanoscale region, leading to a T melt for 1004 atoms that is ϳ34% lower than for bulk Ni.
Figure 2͑a͒ shows that the melting temperature of Ni clusters scales linearly with N Ϫ1/3 .Extrapolating the results in Fig. 2 for finite clusters to Nϭϱ (N (Ϫ1/3) ϭ0), predicts a bulk melting point of 1590 K, significantly below the calculated value of 1760 K for bulk Ni.This is partly because the finite cluster calculations have a free surface, whereas the bulk calculations do not.In order to predict the dependence of T m c , on size, we write the free energies of solid and liquid clusters as the sum of central bulk region and surface.The superscript of b, c, and s represents bulk, cluster, and surface, respectively, 20 where is the free energy per atom of bulk Ni and A is the crystal surface area.Thus, at T m c , the melting temperature of the cluster, we have where we used and S b is per atom.Hence, Assuming a spherical cluster leads to a specific surface area of where V AT is the volume per atom.Thus, where a is a constant, ϭ9.83 J/mol K from our calculations and the average density calculated at the melting temperature from our calculation (ϭ1.09ϫ10 5 mol/m 3 ) leads to (␥ xv Ϫ␥ Lv )ϭ233 mJ/m 2 .Alternatively, we can estimate (␥ xv Ϫ␥ Lv ) from Eq. ͑7͒ for each size ͓derived from Eq. ͑5a͔͒ leading to a surface energy difference of ͑using the fitted value of T m b ϭ1590 K͒.Table I shows the (␥ xv Ϫ␥ Lv ) calculated for various clusters, which increases with cluster size.
Figure 2͑b͒ shows the differences of heat of fusion, entropy of fusion, and melting temperature for clusters compared with bulk as a function of cluster size.Considering that the cluster is composed of surface and bulk regions, and considering that the core part has the same heat of fusion and surface energy as the bulk phase, the enthalpy can be written as This leads to a heat of fusion for Ni clusters that scales linearly with From Fig. 2͑b͒, we see that this linear relation is obeyed above 736 particles, but that the smaller clusters are more b Calculated as the total PE of the cluster minus the bulk PE for the same number of atoms divided by the area using the cluster radius ͑second column͒.stable than what is suggested by ͑9͒.This will be discussed in Sec.IV.From ͑5͒ and ͑9͒, the entropy of fusion should have the form Figure 2͑b͒ shows that this relation ͑the dotted line͒ is obeyed above 736 particles while the smaller clusters have a lower entropy.

B. Surface melting and the Lindemann criterion
The concept of surface melting is that below the bulk T m a quasiliquid skin forms on the surface, which thickens as the temperature increases, leading finally to melting of the whole bulk solid.The bulk melting temperature is then the temperature at which the thickness of the ''skin'' diverges to infinity.
A useful way to relate the origins of melting to atomic phenomena is the Lindemann criterion, 21 which states that melting occurs when the root-mean-square ͑RMS͒ thermal displacement ͑RMSD͒ of the atoms in the lattice reaches a critical fraction ͑typically 10%-15%, depending on the crystal structure 22 ͒ of the equilibrium interatomic distance.The atoms on the surface have weaker restraining forces than the bulk atoms ͑since they have fewer near neighbors͒, allowing them to fulfill this criterion at a lower temperature, suggesting a much lower T melt .
To analyze the calculations in terms of such a model, we defined RMSD as in Eq. ͑3͒, and examined its dependence on the distance from the center of the cluster.Since the interatomic distance in pure Ni is about 2.49 A, we partitioned the atoms into six spherical bins, with dRϭ2.5 Å.The atoms were assigned to bins based on their average initial positions.The average RMSD was calculated for each spherical shell bin. Figure 3͑a͒ shows for various temperatures the average RMSD as a function of distance from the center ͑for the N ϭ1004 cluster͒.In all cases, the RMSD for the surface bin is larger than the central bins.Defining the critical RMSD to be molten as 0.6 Å ͑24% of the bulk interatomic distance of Ni͒, we find that at 1000 K only the surface bin is molten.As the temperature rises to 1100 K and 1150 K the average RMSD for the atoms in the fifth and fourth bins, respectively, reach the critical RMSD, while the atoms in the central three bins remain crystalline.Finally, at 1200 K ͑which is clearly above the melting temperature of the cluster͒ the RMSD for the whole cluster exceeds RMSD of 0.6 Å.Thus, using 0.6 Å as the melting criterion leads to Fig. 3͑a͒, where R crystal is the radius of the remaining solid crystal region.We can also estimate R crystal from the heat of fusion.The potential energy deviates from a linear function of temperature because of melting at the surface.We write this deviation as ⌬PEϭL*N L , where L is the latent heat at T m , and N L is the number of atoms in the liquid phase.This leads to a radius for the remaining solid crystal of the form Indeed, Fig. 4 shows that the radius of the remaining solid crystal from ͑8͒ leads to results very similar to the RMSD criteria when we use a melting criterion of 0.6 Å ͑24% of the interatomic distance in Ni͒.This is much larger than the typical Lindemann criterion of 10-15%.This could be because these small clusters have a large fraction of the atoms on the surface ͑leading to a larger diffusion region for the surface atoms͒, or it could be that 24% is more accurate.Following the changes in RMSD through the melting process, we find a discontinuity of the RMSD for the central atoms, indicating a first order phase transformation ͑with latent heat͒.Figures 3͑b͒ to 3͑d͒ show the snapshots of the cluster at 1000 K, 1140 K, and 1180 K, providing a clear picture of the increased thickness of liquid skin formed on the surface at 1000 K and 1140 K and the melted bulk cluster at 1180 K.

C. Crystallization
Starting with the liquid phase and cooling at a fixed rate, all clusters initially display supercooling below T melt and finally a large drop in the energy upon solidification.From the peak in the derivative of the energy ͑the specific heat͒, we can define the crystallization temperature, T x .
Figure 5 shows the potential energies for the heating and cooling, and reheating cycles of the 336 and 736 atom clusters.In both cases, after cooling to 300 K, the potential energy approaches a value close to that in the heating process ͑which started with a cluster having the FCC structure͒.But after reheating the 336 atoms cluster, the melting tempera-ture increased to 1050 K ͑the original FCC structure led to a melting temperature of 980 K͒, which indicates melting and quenching has led to a phase more stable than FCC.This contrasts with the 736 atom cluster, which leads to the same melting temperature, 1120 K, during reheating.
We developed an FCC searching algorithm to group atoms within one FCC grain and identify the stacking fault and grain boundary. 16In this method we first define an atom having 12 nearest neighbors with an ABC packing sequence as a part of a perfect FCC crystal, from it we check the atoms in far nearest neighbors, and stop at the atoms, which lost the FCC nearest neighbors configuration.Using this algorithm we found a perfect FCC crystal ͓Fig.5͑b͔͒ of about 500 atoms within the quenched 736-atom cluster.This FCC crystal had an orientation unrelated to the initial FCC structure, showing that memory or the original structure was lost.However, this new FCC structure still melts at 1120 K.
On the other hand, for the cluster with Nϭ336 atoms we could not find in the quenched cluster any FCC cluster having more than 13 atoms.Indeed, the 336-atom cluster has a slightly deformed icosahedral structure, as shown in Fig. 5͑c͒.This icosahedral structure is more stable than FCC by ϳ1 kJ/mol melting temperature 70 K ͑7%͒ higher than FCC.The Nϭ336 cluster is very close to the magic number of 309 atoms for a Mackay cluster with four shells, which could be the reason an icosahedral structure is more stable.On the other hand, it could be that the Nϭ736 atom cluster forms an FCC crystal easily because it is far from the nearby two magic numbers 561 and 923.Alternatively, icosahedral could be generally more favorable than FCC for sizes below 400 to 700.
The Honeycutt-Andersen (HA) analysis: To analyze the local ordering of the various clusters, we used the Honeycutt-Andersen ͑HA͒ algorithm 19 to characterize the environment for each pair of atoms.Table II shows the HA pair parameters for various structures ͑including FCC, HCP, and various icosahedral clusters͒ of Ni at 300 K. We used the minimum between the first and second nearest neighbor distance in FCC crystal as the cutoff for bonded atoms, 2.99 Å for Ni, and 3.55 for AgCu.We see that 1551 and 2331 HA pair types are characteristic of icosahedral ordering.Here 1551 corresponds to two neighbor atoms with five common neighbors forming a bonded pentagon, while 2331 corresponds to a pair of not bonded atoms that share three neighbor atoms forming a bonded triangle.Two other HA pair types, 1421 and 1422, characterize the FCC and HCP crystal structures.Normalizing by the total number of all nearest neighbor pairs, the perfect FCC structure contains 100% of 1421 pairs but no pairs of 1551 or 2331 type.
First we analyze the clusters with Nϭ336 and Nϭ736 atoms at 300 K formed before and after the heating-cooling cycle.Comparing the solidlike cluster of Nϭ336 atoms after quenching with the icosahedral cluster of 309 atoms, we see that they have almost the same ratio of every pair parameter.This suggests that the cluster with 336 atoms has a high content of fivefold symmetry, which is compatible with an icosahedral structure ͑it definitely does not have FCC-like local order͒.
Figure 6 shows the HA analysis for types 2331, 1551, and 1421 for the clusters with 336 and 736 atoms during the heating-cooling cycle.
Figure 6 also includes the analysis for the metallic glass formed in the AgCu system ͑studied in detail in Ref. 15͒.In the initial FCC structure, there were no 2331 or 1551 pairs but 70% of 1421 pairs ͑less than 100% because of the atoms on surface͒.As the temperature increases, the population of 1421 pairs decreases, dropping very quickly after surface melting begins and going close to zero after the cluster melts totally.Upon cooling AgCu, the 1421 pairs stay below 10% for the whole cooling process during which AgCu forms a metallic glass with T g ϭ550 K. 15 Below 550 K, the 2331 and 1551 pairs continue increasing with decreasing temperature, reaching values even larger than for the 55 atom icosahedral cluster ͑38% for 2331 pair and 10% for 1551 pair͒, but less than the 13 atom icosahedral cluster ͑71% for 2331 pair and 29% for 1551 pair͒.This indicates that the metallic glass has a great deal of icosahedral local structure.
For the 736-atom cluster the cooling process recovers 70% of the 1421 pairs after the new FCC crystal is formed.However, the 336 atoms cluster reaches only 33% of the number 1421 pairs characteristic of FCC ͑close to the ratio of 1421 pairs in the 309 icosahedral cluster͒.
The 2331 and 1551 pairs show quite interesting properties during the phase transition.During the heating process the fraction of both kinds of pairs increases at first and then drops upon melting of the whole structure.In the cooling process, the undercooled liquid shows an increasing fraction of 2331 and 1551 pairs ͑observed also in a Lennard-Jones system 12 ͒ until crystallization.For the Nϭ736 atom cluster, FIG. 5. ͑a͒ Potential energy U(T) as a function of heating from 300 K ͑starting from an FCC structure͒, 1600 K and then cooling back to 300 K. We then reheated the cluster to 1600 K.For 336 atoms, the heated and cooled cluster has an isocahedral structures, and this reheating shows that the isocsheadrla structure has a T melt ϳ1050 K compared to 980 K for the original cluster.For the clusters with 736 atoms, the cooled cluster has a large retion that is FCC and the reheating leads to T melt ϳ1120 K, just as for the original heating.͑b͒ ͓001͔ Projection of the new FCC crystal formed in the 736 Ni atoms cluster from heating and then cooling.Here the 494 dark balls form the new FCC crystal ͑orientation not related to the original crys-tal͒.The white balls are disordered.͑c͒ Snapshot of solid cluster with 336 atoms at 300 K.This is dominated by isosaheral character.we find that below the crystallization temperature of 600 K the amount of 2331 drops to 3% and the 1551 pairs drops to zero.However, for the 336 atom cluster, 2331 only drops to 18% while 1551 drops only to 2% ͑both close to the values of 19% 2331 pairs and 2% 1551 pairs in the icosahedral cluster with Nϭ309 atoms͒.

IV. DISCUSSION
We find that for Ni the thermodynamic properties relevant to melting and surfaces behave in a very regular way for clusters above 736 atoms, with the difference between the enthalpy and melting temperature for the clusters deviating from the bulk value by an amount proportional to N Ϫ1/3 as expected from the ratio of surface atoms to bulk atoms.We refer to this as the mesoscale regime.However, below ϳ500 atoms there are significant deviations from the trends observed in the mesoscale regime.We find that these deviations in the small cluster ͑molecular͒ regime are caused by a transition of the cluster to a structure different than for the bulk and mesoscale system.For Ni, the small clusters prefer a structure with significant icosahedral structure.For very small clusters ͑147 and 70 atoms͒, the icosahedral structure is obtained by annealing at moderate temperatures.For larger structures, it was necessary to melt the initial FCC structures in order to transform to the more stable icosahedral form.In any case the structure stabilized for small clusters leads to significant deviations from the thermodynamic properties extrapolated from the mesoscale regime.We will discuss below first the mesoscale regime and then the small cluster regime.

Decrease of melting point with decreasing particle size
The decrease of the melting point with particle size has been observed in experiments on supported clusters of Au, 1 Sn, 7 Pb, 9 and Pt, 10 for isolated Sn particles ͑with 500 atoms 8 ͒, and for gas phase Na particles ͑with only 139 atoms 11 ͒.][26][27][28][29][30] We find that the change in (T m b ϪT m c ) with size is quite proportional to N Ϫ1/3 for clusters with N above 736 atoms, which prefer the FCC structure.

Latent heat of fusion
The decrease of the latent heat of fusion with size, has been observed by calorimetric measurements on small Sn particles 7,8 and by numerical simulations. 28As shown in Fig. 2͑b͒, we find that the decrease scales as N Ϫ1/3 for clusters larger than 500 atoms.
We observed that for clusters with ϳ100 atoms, the FCC cluster transforms during heating to a lower energy icosahedral structure.This icosahedral cluster then melts at higher temperature and with a larger heat of fusion than the FCC cluster.A similar increased latent heat is exhibited for isolated tin clusters with about 500 atoms.It is not clear in our simulations that such a transition occurs for clusters with 484 and 363 atoms, but they might favor icosahedral structure leading to a higher heat of fusion, and a smaller ⌬H f b Ϫ⌬H f c , as found in Fig. 2͑b͒.

Critical size for zero latent heat of fusion
When the latent heat drops to zero, the cluster fluctuates between solid and liquid, leading to a phase transition that spreads over a finite temperature range.This results because of ͑i͒ formation of a size dependent liquid skin, and ͑ii͒ thermodynamic fluctuations, which become significant in small systems.We find that this skin grows with temperature.Based on the linear fit to the mesoscale FCC regime, we estimate that heat of fusion would go to zero for clusters of ϳ86 atoms, leading to a melting temperature of 633 K.According to Turnbull's glass forming ability criteria, a good glass former will have a ratio of glass transition temperature ͑Tg͒ to the melting point increased from values near 1/2 to 2/3, such that the homogeneous nucleation of crystals in the undercooled melt should become very sluggish.It's believed that Tg for pure metal is lower than 1/2 Tm, in Ni system it should be around 600 K. 31 For bulk Ni, this is only 1/3 of the melting point, making so it very difficult to form a bulk amorphous phase of Ni.However, if the melting temperature can be decreased to ϳ600 K, we expect that it would be easy to form glassy Ni.
For Ni clusters we did not find a clear formation of the amorphous phase, because FCC is not the stable structure for Ni clusters smaller than 500.Since the small FCC clusters ͑up to 147 atoms͒ transform into the icosahedral phase during heating, we obtain a larger melting temperature and latent heat for melting of the icosahedral phase.Thus the size for which the heat of fusion becomes zero for the icosahecral phase is likely much smaller than 86 atoms.

Surface melting
The role of surface melting in nanoclusters has also been studied for Pb particles on a Si substrate. 9These studies show that a ϳ0.5 nm liquid skin grows just below the melting temperature of a 50 nm Pb crystal, indicating size dependent melting.These results are consistent with our simulations.
The driving force for melting on a flat surface is the reduction in the total interfacial energy, 29,6 ⌬␥ϭ␥ xv Ϫ͑␥ LV ϩ␥ XL ͒Ͼ0, ͑12͒ where ␥ is the surface energy between the solid ͑X͒, liquid ͑L͒, and vapor (V) phases.For Ni, the experimental driving force ⌬␥ is close to zero 6 so that a surface melting on a flat surface is not expected.However, in our simulation the surface energy is decreasing with the radius of the particles, and since geometrical and capillary effects, might increase the thickness of quasiliquid skin in small particles 32 and decrease the interface energy of solid and quasiliquid.As a result, the small radius of these small clusters could favor surface melting, as we have seen in this study.Surfaces are expected to play an especially important role in crystallization processes, where optimizing the internal packing and minimizing surface energy are competitive driving forces.Thus, to understand the crystallization of a free nanocluster we must understand homogeneous nucleation.Below some critical radius ͑ϳ10-20 Å͒ the optimum structure might change from the bulklike structure ͑say, FCC͒ optimal for long range order to a new arrangement optimum for a cluster environment. 33,34 Small cluster "molecular… regime

Comparison with previous simulations
As shown in Fig. 5͑a͒, for 336 atoms the icosahedral structure of Ni melts at a temperature 70 K ͑7%͒ higher than the FCC structure.As a result, the N Ϫ1/3 rule does not apply at these small sizes.We also found that for very small clusters, such as Nϭ140, Nϭ87, the FCC structure transforms to icosahedral before melting, which also leads to a higher melting temperature than the Tm predicted from of the extrapolation of the N Ϫ1/3 rule.Since these icosahedral structures are formed during natural quenching or even heating processes, they are more stable than FCC.It is reasonable to propose that below a critical size, ϳ500 atoms in our study, the icosahedral phase actually is the stable phase favored by the small clusters.Thus these clusters will have a melting temperature higher than predicted by the N Ϫ1/3 extrapolation from FCC particles.This leads to melting temperature below the critical size that may depend on size differently than N Ϫ1/3 .
Cleveland and Landman reported in detailed studies of the size dependence of energetic and structure of small nickel clusters, 38 where they found icosahedral clusters to be favored for clusters smaller than 2300 atoms ͑EAM force field͒ and 1600 atoms ͑LJ force field͒.Since only the minimized energy was considered, the critical sizes are for 0 K.We quenched from a liquid drop and also observed the formation of icosahedral clusters.However, we find that for our FF at 300 K and a cooling rate of 4 K/ps, the critical size for forming icosahedral is ϳ500.

The Bachels model
Bachels et al. have combined their measurement of melting temperature on isolated tin nanoparticles 8 with experiments by Lai on supported tin particles 7 to propose a model in which below a critical radius ͑35 A͒, the melting temperature of the particle lies below the surface melting temperature, such that the surface premelting effect will be suppressed, leading to a constant melting temperature. 8But there is still a large difference between the predictions from their model and the experimental data on melting temperature for a cluster with 500 atoms, and they suggested that the difference is due to insufficient accuracy of the parameters used for evaluating their model.
However, we note that the actual melting temperature of tin particle with ϳ500 atoms is larger than that predicted from size dependence rule (N Ϫ1/3 ) based on nano particles ranging from 5-50 nm.This agrees with our finding that the critical size should be related to the icosahedral phase ͑or other stable structure͒ for small particles.It is very likely that the initial solid state cluster structure in the experiments with ϳ500 atoms were already in the icosahedral structure, and hence it is expected that the melting temperature will be higher than predicted by extrapolated from the larger nanoparticles of size 5-50 nm.

Simulations on Au particles
From simulations Cleveland et al. 30 showed that clusters smaller than 500 atoms lead to properties very different than larger particles.They studied gold clusters with 75, 146, and 459 atoms and found a solid ͑truncated-decahedra for Au75 and Au146, and FCC with truncated-octahedral for Au459͒to-icosahedral transformation.They did not find that a quasiliquid wetting layer formed on the surface, which thickens until the whole cluster melts in clusters Au75 and Au146, but there is solid-liquid coexistence in Au459.They calculated the melting temperature decreases with size ͑T melt ϳ550 K for Nϭ75, T melt ϳ625 K for Nϭ146, and T melt ϳ760 K for Nϭ459͒, all below that for bulk system (T melt ϭ1090 K).However, they did not report results for larger clusters, where the size dependence of melting temperature might be described by the N Ϫ1/3 rule.Based on our results, we expect that the melting temperatures for these cases are higher than the melting temperature of FCC structures having the same sizes.

V. CONCLUSION
Using MD simulations with the Q-SC FF on unsupported Ni nanoclusters we find that for the mesoscale regime above ϳ750 atoms, the melting temperature and heat of melting scales inversely as N 1/3 for FCC structures.We find that melting proceeds from the surface inwards and that the melting process for the central core is discontinuous ͑first order͒.Thus the mesoscale regime is characterized by surface melting.
Upon quenching from the melt, we find that larger clusters and bulk Ni crystallize easily to form FCC crystals.However, for smaller clusters, with less than ϳ500 atoms, we find that the stable structure is icosahedral, which has a higher melting temperature and larger heat of fusion than FCC structures of the same size.
We found that the supercooled liquid and glassy phase are highly icosahedral as is the premelted phase.

FIG. 1 .
FIG. 1. Potential energy U(T) and heat capacity Cp(T) for heating and cooling cycles of a cluster with N ϭ1004 Ni atoms and the bulk FCC phase of Ni.

Figure 2͑a͒ shows that
Figure 2͑a͒ shows that ͑6a͒ leads to a good fit to the calculations for the range of 336 to 8007 atoms, with a ϭ4220.2K. Using this value of a with S f bulk ϭ⌬H f bulk /T m b

c
FIG.2.͑a͒ Dependence of melting temperature on cluster size.The dashed line shows the best fit to a linear function of N (Ϫ1/3) .This leads to a predicted value for large N of 1590 K, well below the value 1760 K calculated for the bulk system ͑no surface͒.The melting temperature for N ϭ336͑ico͒, Nϭ141, and Nϭ87 are all higher than the Tm predicted from linear fitting.͑b͒ Tm b -Tm c , ⌬H f b -⌬H f c , and ⌬S f b -⌬S f c plotted as a function of N (Ϫ1/3) .The solid and long dashed lines show the fit of T and ⌬H, to a linear function of N (Ϫ1/3) .Here we see an excellent fit for N above 736 atoms.The dotted line is the fit to Eq. ͑10͒, which should be valid above 736 atoms.

FIG. 3 .
FIG. 3. ͑a͒ Average displacements for the atoms in various spherical shells, as a function of the distance from the center of the cluster.This shows surface melting starting at 1000 K ͑using RMSDϭ0.6A as the criterion for melting͒, and the cluster has already melted by 1200 K. ͑b͒ Snapshot of the 1004 Ni atom cluster at 1000 K, showing formation of the liquid ''skin'' at the outer layer.͑c͒ Snapshot of 1004 Ni atoms cluster at 1140 K, showing that the inner regions are still ordered.͑d͒ Snapshot of 1004 Ni atoms cluster at 1180 K, showing that the entire cluster is melted.

FIG. 4 .
FIG. 4. Radius of the remaining crystal in Nϭ1004 cluster as a function of temperature, melting process, calculated from two methods: The circles indicate the values obtained from assuming a Lindemann criterion of RMSDϾ0.6 Å as an indication of metling.The line is based on Eq. ͑11͒ based on analyzing the dependence of potential energy temperature.

FIG. 6 .
FIG. 6. Honeycutt-Andersen pair populations as a function of T. Shown are the 336 and 736 atom clusters, and the bulk CuAg system ͑which forms a glass upon cooling to 550 K͒. ͑a͒ 1421 HA pairs showing FCC character, ͑b͒ 2331 HA pairs showing icosadedral character, ͑c͒ 1551 HA pairs showing icosahedral character.

TABLE I .
Thermodynamic properties evaluated from the MD for various clusters as a function of size and the number of atoms within the cluster.

TABLE II .
The analysis of various structures in terms of Honeycutt-Andersen ͑HA͒ pairs ͑based on the final structure at 300 K͒.