Anomalous Size Dependence in the Melting Temperatures of Free Sodium Clusters: An Explanation for the Calorimetry Experiments

Anomalous Size Dependence in the Melting Temperatures of Free Sodium Clusters: An Explanation for the Calorimetry Experiments Andre ´s Aguado*and Jose ´M. Lo ´pez Departamento de Fı ´sica Teo ´rica, Universidad de Valladolid, Valladolid 47011, Spain (Received 11 December 2004; published 14 June 2005) The meltinglike transition in unsupported NaNclusters (N55, 92, 147, 181, 189, 215, 249, 271, 281 and 299) is studied by first-principles isokinetic molecular dynamics simulations. The irregular size dependence of the melting temperatures Tmobserved in the calorimetry experiments of Schmidt et al. [Nature (London) 393, 238 (1998)] is quantitatively reproduced. We demonstrate that structural effects alone can explain all broad features of experimental observations. Specifically, maxima in TmNcorrelate with high surface stability and with structural features such as a high compactness degree. DOI: 10.1103/PhysRevLett.94.233401 PACS numbers: 36.40.Ei, 31.15.Qg, 36.40.Sx, 64.70.Dv There is a fundamental interest in understanding the analog of the melting phase transition in small atomic clusters. Classical arguments predict that melting temperatures Tmof small particles should show a monotonic decrease from the bulk limit as their sizes shrink. However, Jarrold and co-workers [1] have demonstrated that small gallium and tin clusters melt at temperatures higher than Tbulk m. Also, calorimetry experiments by Haberland and co-workers [2] show that the size dependence of Tmis not monotonic for NaNclusters in the size range N50–350. The observed pattern of maxima and minima in the TmNcurve cannot be fully explained either by electronic or geometric (icosahedral) shell closings. Despite much theoretical effort [3–9], no definite explanation for the calorimetry experiments has emerged yet. In this Letter, we show that molecular dynamics (MD) simulations, with an accurate (first-principles) representation of interatomic forces, reproduce and provide an explanation for the calorimetry results. We employ an orbital-free (OF) version of density-functional-theory (DFT) [10], which expresses the energy as an explicit functional of the valence electron density, and pseudopotentials to represent the ionic field acting on the valence electrons. The technical details of our OF-DFT implementation are explained in previous work [5,11]. Here we describe just the main improvements incorporated into our energy functional. Compared to Kohn-Sham (KS) [12] calculations, two new ingredients appear in OF-DFT: (i) an explicit expression for the functional dependence of the noninteracting electron kinetic energy on density, Tsn, and (2) if pseudopotentials are employed, they must be local, that is, dependent only on spatial position and not on orbital angular momenta which are simply not defined in OF-DFT. In previous works [5], we employed the second-order gradient expansion for Tsn[13] and the local pseudopotential of Fiolhais et al. [14]. Those pseudopotentials were adjusted to reproduce bulk properties within linear-response-theory and might not be optimal for cluster DFT studies. In this work, the Tsnfunctional is given by [11] TsnTvWnTn;(1) T3 10 Zd~ rn~ r5=32~ k~ r2;(2) ~ k~ r2k0 f3Zd~ sk~ sw2k0 Fj~ r~ sj;(3) where TvW 1 8Rd~ rjrnj2 nis the von Weizsa ¨cker functional, k~ r321=3n~ r,k0 fis the Fermi wave vector corresponding to a mean electron density n0,0:51, and wxis a weight function, determined by requiring Tsn to recover the correct expression for both the uniform density and linear-response-theory limits. Further details may be found in Ref. [11]. This functional provides results of KS accuracy for bulk alkalis [15] and also at alkali surfaces [16]. Its application here to clusters is new. The mean electron density is defined as n0V=N, where Vis the volume of a sphere of radius equal to the mean gyration radius of the cluster. The appropriateness of this expression will be ultimately tested by our ability to reproduce ab initio KS results (see below). Regarding the pseudopotential vps, we retain the evanescent-core functional form of Fiolhais et al. [14], but the two parameters entering vps (named and Rin Ref. [14]) are reobtained by fitting to KS calculations. Specifically, we have performed KS static calculations on 20 representative cluster geometries obtained by preliminary MD simulations (including different sizes and structures—icosahedra, decahedra, disordered, etc.—at different temperatures). KS calculations are performed with the SIESTA code [17], employing nonlocal Troullier-Martins pseudopotentials [18] and the LDA approximation [19,20] to the exchange-correlation functional EXCn(this is the same EXCnemployed in the OFDFT calculations). We request that OF-DFT calculations, performed with our local vps, reproduce SIESTA results for the interatomic forces and energy differences between PRL 94, 233401 (2005) PHYSICAL REVIEW LETTERS week ending 17 JUNE 2005 0031-9007=05=94(23)=233401(4)$23.00 233401-1 2005 The American Physical Society cluster geometries. For the optimal vps thus obtained (opt 3:799;Ropt 0:445), the average deviation between OF and KS forces is smaller than 5%. The energy differences are reproduced with an average error smaller than 1%. Figure 1 provides an explicit demonstration that OF forces are of KS accuracy. This level of agreement is uniform with respect to cluster size and temperature and is the same for structures other than those included in the fitting set. As a first step, we have located candidate ground state isomers of NaNclusters, with N55, 92, 147, 181, 189, 215, 249, 271, 281, and 299. Each cluster is placed in a unit cell of a cubic superlattice of length L(L25 a:u:for Na55 and 48 a.u. for Na299). The set of plane waves periodic in that superlattice, up to an energy cutoff of 20 Ryd, is used as a basis set to expand the valence electron density. Forces on atoms are calculated after the ground state density has been found. Minimum-energy isomers were located by simulated annealing, performed at a cooling rate of 0:2K=ps. As starting geometries, we have employed high-Tliquidlike (disordered) and solidlike structures with icosahedral, decahedral, and cuboctahedral symmetries. The initial temperature Tof the solidlike structures was the highest one which does not lead to melting of the cluster. For all sizes, lowest-energy structures were found by annealing of the icosahedral structures. The final structures (excepting Na55,Na147, and Na299) do not resemble incomplete Mackay icosahedra, but rounded geometries which try to attain an optimal packing and be as spherical as possible. Two examples of these structures are shown in Fig. 2. A subsequent heating-cooling run of these structures did not locate lower-energy isomers. For clusters of this size, it is quite unlikely that we find the absolute minimumenergy isomer, irrespective of the sampling method. We have chosen simulated annealing because it provides a reasonable efficiency=cost ratio and is sufficiently inexpensive that it allows for the employment of first-principles forces at all stages of optimization. For each cluster size, we have performed isokinetic Born-Oppenheimer MD runs, in which the average kinetic energy is kept constant by velocity rescaling and the atomic forces are evaluated from the Hellmann-Feynman theorem. The time step is 3 fs, and the total simulation length for each size is 5 ns (for those temperatures close to the transition region, the simulations were 300 ps long, while shorter runs of 150 ps were employed at other temperatures). We employ multiple histogram techniques [21] in order to extract smooth caloric and specific heat curves. Statistical sampling is necessarily less complete than that achieved in Ref. [21] due to the computational expense of first-principles calculations, but we have checked that it is adequate to extract accurate caloric curves. A representative example of these curves is given in Fig. 3, while melting temperatures Tm(read from the specific heat maxima) and latent heats lm(estimated from the step height between liquid and solid branches of the caloric curve) are given in Table I and in Fig. 4. Both values are in quantitative agreement with the experimental determinations [2]. Although not shown explicitly, this means that the entropies of melting also agree with experiment. -0.002 0 0.002 0.004 0.006 Fx (a.u.) -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 Fy (a.u.) 010203040 50 60 Atom number -0.006 -0.004 -0.002 0 0.002 0.004 Fz (a.u.) FIG. 1 (color online). Comparison between SIESTA (full symbols) and OF-DFT (continuous line) interatomic forces, for a Na55 atomic configuration extracted from a run at 150 K, which was not directly employed in the pseudopotential fitting. FIG. 2 (color online). OF-DFT ground state isomers of Na181 and Na215 found by simulated annealing. Atoms at the surface are represented by light (golden) spheres and interior atoms by dark (blue) spheres. 160 180 200 220 240 260 280 T (K) -6,06 -6,04 -6,02 -6 Total Energy (eV/atom) 160 180 200 220 240 260 280 T (K) 4 6 8 10 12 14 16 Cv/NkB latent heat T* FIG. 3. Caloric curve (left) and specific heat per particle in units of kB(right) of Na299. The latent heat lmis obtained from the caloric curve as shown. Tsignals the beginning of melting. PRL 94, 233401 (2005) PHYSICAL REVIEW LETTERS week ending 17 JUNE 2005 233401-2 This level of agreement has not been achieved in previous simulations over a broad size range. The specific heat curves are similar for all sizes and contain a single peak which is slightly wider on the low-temperature side. Analysis of the diffusion constants of different sodium atoms indicates that melting initiates at surface atoms, providing an explanation for the peak asymmetry. Surface and homogeneous melting temperatures are so close that they merge into a single specific heat peak. The absence of premelting signatures in the experimental melting curves is thus reproduced by the OF-DFT calculations. In the following, we try to identify some energetic and/or structural trends that help to rationalize the oscillations in TmN. If, as we observe, melting nucleates at the cluster surface, it is sensible to think that the stability of surface atoms plays a role in determining Tm. This idea has already been put forward by Lee et al. [7] within the context of parametrized potential models. As a measure of surface stability, we employ the average surface evaporation energy, defined as EevapN 1 NsNs i1Ei evap, with Nsthe number of surface atoms and Ei evap EN EiN1ENa the energy required to remove atom ifrom the cluster surface. This is a zero temperature calculation, and EiN1is obtained from local structural relaxation to the nearest minimum after evaporation of atom i. As a measure of the stability of a surface atom, relative to an average atom in the cluster, we show in Table I the quantity NEevapNEbN, where EbNENNENa Nis the average binding energy per atom. Nprovides a measure of relative surface stability against evaporation. Although a corresponding measure involving diffusion instead of evaporation would be more appropriate in a discussion of melting, Nprovides at least a rough idea of the relative strength of surface bonds. Ntakes positive values precisely for the same sizes where a maximum in Tmis observed; moreover, the larger the value of N, the higher the melting temperature. This suggests that melting temperatures of sodium clusters are directly related to relative surface stability. It also agrees with the experimental observation that entropic effects play only a secondary role in determining melting temperatures, as compared to energetic effects [2]. Trying to find a structural origin for the enhanced surface stability of selected cluster sizes, we have evaluated the average volume per atom, defined as V=N 4r3 3N, with rthe gyration radius of the cluster. Let us note that this quantity decreases with decreasing cluster size as the average interatomic distance contracts. Figure 4 shows the variation of V=N with cluster size, together with the best linear fit to the data. It is immediately appreciated that clusters with an ‘‘anomalously’’ high melting temperature are more compact. The enhanced surface stability is thus shown to be induced by a higher compactness degree, which can be obtained only for some particular sizes. In order to proceed further, we have divided each cluster into surface and ‘‘bulk’’ atoms (as shown in Fig. 2) and evaluated the average distance between bulk atoms (dbb), between surface atoms (dss), and between surface and bulk atoms (dsb). We have found that dsb shows a clear correlation with V=N (see Table I): the higher the melting temperature of a given cluster, the shorter its dsb value. Table I also shows that clusters with shorter dsb have larger interatomic distances between surface atoms. As surface tension is driven by the tendency of dss to decrease (which contributes to the internal cluster pressure), those clusters with higher Tmhave also larger surface tensions. As the cluster is heated, we observe that dsb distances expand, while dss slightly contract (which leads to the expected decrease in surface tension as temperature increases). At the very initial stages of the meltinglike transition (identified from the initial increase in the specific heat curve; see Fig. 3), Table I shows that ddssTdsbTis very 0.005 0.01 0.015 1/N 13 13.5 14 14.5 15 V/N (A3) 50 100 150 200 250 300 N 200 220 240 260 280 Tm (K) Experiment OF-DFT FIG. 4 (color online). Size variation of the volume per atom (top) and melting temperature (bottom) in sodium clusters. The dashed line in the upper panel is the best linear fit to the data. TABLE I. Size variation of melting temperatures Tm, latent heats lm, and several energetic and structural properties described in the text. Size Tm(K) lmeV=atom(eV) dsb (A ˚)dss (A ˚)dT 55 280 0.011 0.122 3.552 3.683 0.009 92 206 0.004 0:182 3.612 3.670 0.008 147 256 0.014 0.057 3.604 3.708 0.005 181 252 0.013 0.037 3.603 3.709 0.006 189 217 0.006 0:095 3.640 3.698 0.002 215 249 0.011 0.012 3.612 3.716 0.004 249 225 0.006 0:203 3.645 3.708 0.003 271 245 0.007 0.009 3.638 3.710 0.002 281 222 0.008 0:067 3.650 3.704 0.001 299 251 0.013 0.014 3.615 3.720 0.002 PRL 94, 233401 (2005) PHYSICAL REVIEW LETTERS week ending 17 JUNE 2005 233401-3 close to zero. This result is independent of size and thus serves to identify a systematic behavior in the melting transition of Na clusters: all clusters show a reduction of dupon heating and all melt when dis close to zero, but different sizes approach this critical stage at different temperatures. Our calculations thus predict a simple systematics for the melting behavior of small sodium clusters in a broad size range. It is a well-known fact that surface metal atoms tend to undergo inwards relaxation (bond length contraction) in order to compensate for the reduced electronic density at the surface. At the same time, there is a strong energy penalty for dangling atoms, which would possess a very low coordination. Satisfaction of these rules leads to ground state isomers with smooth surfaces (without any surface steps or dangling atoms; see Fig. 2) which try to optimize ionic packing. The allocation of all surface atoms into a single, rounded surface shell induces oscillations in dss and dsb as a function of size. When dss is relatively large, dsb may be more efficiently reduced, resulting in an additional stabilization of the cluster surface which inhibits melting. For example, the ‘‘anomalous’’ very high melting temperature of Na55 is simply due to its higher compactness relative to other sizes. Garcı ´a-Gonza ´lez et al. [22] have shown that Tsnfunctionals with a higher degree of nonlocality than the one employed here would be needed in order to reproduce electronic shell closing effects. Our results thus show that electronic effects can play only a secondary role and are not needed in an explanation of the main Tmoscillations. We have calculated the energy eigenvalue gap at the Ne 58 electronic shell closing for the Na55 cluster, for both solid and liquid configurations at temperatures close to Tm, using SIESTA. These values are averages over 10 solidlike and 10 liquidlike independent configurations extracted from the MD runs. The results (Esolid gap 0:23 eV and Eliquid gap 0:20 eV) show that both solid and liquid phases are similarly ‘‘magic’’ from the electronic point of view, which explains why electronic shell effects do not affect the melting transition significantly. A recent analysis of calorimetric and photoelectron experiments by Haberland et al. [23] also supports the idea that the size variation of Tmis controlled by geometric effects, providing further, independent support for our conclusion. In summary, we have reported first-principles MD simulations that reproduce the main oscillation pattern observed in the calorimetric determination of melting temperatures, latent heats, and entropies of melting of small sodium clusters. The precise location of the maxima in the melting temperature [for example, TmNa142>T mNa147] continues to be an issue and will be the subject of future analysis. Our results provide just a plausible explanation of the experimental ‘‘anomalies’’ and find a systematic behavior in the melting of these clusters. Because of the practical impossibility of locating the exact ground state isomer of clusters with hundreds of atoms, alternative explanations might still be possible. This work was supported by Junta de Castilla y Leo ´n (Project No. VA073/02), by DGES (Project No. MAT200204393-C02-01), and by the ‘‘Ramo ´n y Cajal’’ program of the Spanish Ministry of Science and Technology. *Electronic address: [email protected] [1] A. A. Shvartsburg and M. F. Jarrold, Phys. Rev. Lett. 85, 2530 (2000); G. A. Breaux et al.,ibid. 91, 215508 (2003). [2] M. Schmidt et al., Phys. Rev. Lett. 79, 99 (1997); Nature (London) 393, 238 (1998); C.R. Physique 3, 327 (2002); Phys. Rev. Lett. 90, 103401 (2003); R. Kusche et al., Eur. Phys. J. D 9, 1 (1999). [3] P. Blaise et al., Phys. Rev. B 55, 15 856 (1997). [4] A. Rytko ¨nen et al., Phys. Rev. Lett. 80, 3940 (1998); K. Manninen et al., Phys. Rev. A 70, 023203 (2004). [5] A. Aguado et al., J. Chem. Phys. 111, 6026 (1999); J. Phys. Chem. B 105, 2386 (2001); Eur. Phys. J. D 15, 221 (2001). [6] F. Calvo and F. Spiegelmann, Phys. Rev. Lett. 82, 2270 (1999); J. Chem. Phys. 112, 2888 (2000); 120, 9684 (2004). [7] Y. J. Lee, E. K. Lee, S. Kim, and R. M. Nieminen, Phys. Rev. Lett. 86, 999 (2001). [8] J. A. Reyes-Nava, I. L. Garzo ´n, and K. Michaelian, Phys. Rev. B 67, 165401 (2003). [9] K. Manninen et al., Eur. Phys. J. D 29, 39 (2004). [10] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). [11] A. Aguado et al., in Progress in Chemical Physics Research (Nova Science Publishers, Inc., New York, to be published). [12] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). [13] J. P. Perdew, Phys. Lett. A 165, 79 (1992). [14] C. Fiolhais et al., Phys. Rev. B 51, 14 001 (1995); 53, 13 193(E) (1996). [15] D. J. Gonza ´lez et al., Phys. Rev. E 69, 031205 (2004); J. Blanco et al.,ibid. 67, 041204 (2003). [16] D. J. Gonza ´lez et al., Phys. Rev. Lett. 92, 085501 (2004). [17] J. M. Soler et al., J. Phys. Condens. Matter 14, 2745 (2002). [18] N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991). [19] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). [20] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). [21] P. Labastie and R. L. Whetten, Phys. Rev. Lett. 65, 1567 (1990). [22] P. Garcı ´a-Gonza ´lez et al., Phys. Rev. B 53, 9509 (1996). [23] H. Haberland et al., Phys. Rev. Lett. 94, 035701 (2005). PRL 94, 233401 (2005) PHYSICAL REVIEW LETTERS week ending 17 JUNE 2005 233401-4