J. Phys. Chem. B 2001, 105, 11605-11614
11605
Melting and Freezing of Gold Nanoclusters† Yaroslav G. Chushak and Lawrence S. Bartell* Department of Chemistry, UniVersity of Michigan, Ann Arbor, Michigan 48109 ReceiVed: March 12, 2001; In Final Form: July 3, 2001
Molecular dynamics simulations were performed on series of gold nanoclusters comprised of 459, 1157, and 3943 atoms, to study their structures and properties during heating and cooling. The increased depression of melting point as particle size decreases has been interpreted in terms of Pawlow’s triple point theory, the liquid shell model, and extensions of the two. The solid-liquid interfacial free energy σsl of ∼0.15 J/m2 inferred from the liquid shell model was close to the values predicted by several empirical relations. When molten particles in a given series were frozen, several different final structures were obtained even though conditions had been identical. Most of the frozen clusters attained an icosahedral structure, a structure found to be stable to mild annealing. Other structures observed were imperfect truncated decahedra, twinned truncated octahedra, and hexagonal close packed structures. The rate of crystal nucleation in liquid clusters was analyzed in terms of classical nucleation theory and the diffuse interface theory. The interfacial free energy of 0.084 J/m2 inferred from the kinetics of freezing was appreciably lower than that inferred from the melting points. Since existing theories of the melting of small particles and the nucleation of crystals in them are far from rigorous, a close agreement between the quasi-thermodynamic and kinetically based free energies cannot be taken for granted. Sources of error in the models of melting are discussed.
Introduction Nanoclusters, or nanoparticles containing tens to thousands of atoms or molecules, play important roles in science and technology. Applications of clusters in innovative technology seem to be endless, as they range from industrial catalysis to the miniaturization of electronic devices. They have also been attractive objects for study from the academic viewpoint because they occupy a niche between single molecules and bulk condensed phases. Properties of nanoparticles strongly depend on their internal structure as well as their size. Experiments and computer simulations on clusters of polyatomic molecules such as CCl4,1,2 NH3,3 CO2,4 and hexafluorides5-9 have shown that nanoparticles containing as few as 100 molecules exhibit bulklike crystalline structures. On the other hand, atomic metal or rare gas clusters up to several thousand atoms in size exhibit noncrystalline icosahedral (Ih) or other multiply twinned structures (MTP) even though the bulk materials are crystalline. Such a dramatic difference in structure between atomic and molecular nanoparticles may depend on the effective range of the intermolecular forces relative to the distances between centers of mass.10,11 Alternatively, this difference has been attributed to the anisotropy of the repulsive part of the intermolecular potential.4 A fully satisfactory resolution of the problem has not been achieved. Considerable effort, both experimental and theoretical, has been devoted to the study of the structure of atomic clusters. One difficulty in experiments is that by most experimental techniques it is practically impossible to prepare a large ensemble of nanoparticles with a narrow size distribution. Furthermore, even for a given particle size, a distribution of structural isomers occurs. Therefore, the interpretation of experimental data can be very complicated.12 That is why computer simulations are considered as not only complementary †
Part of the special issue “Howard Reiss Festschrift”.
but also essential components of structural studies. The majority of simulations have focused on exhaustive searches for the lowest energy configuration (global minimum on the energy surface) at T ) 0 K and possible changes in the optimal configuration related to particle size. Inasmuch as the formation and growth of clusters are usually nonequilibrium processes, a distribution of structural isomers is generally encountered even for a given cluster size. Final atomic arrangements strongly depend both on thermodynamic factors and on growth kinetics. Customarily, crystalline clusters in computer simulations are generated by construction of idealized structures which may be subjected to thermal agitation, or by the freezing of molten clusters. In a recent molecular dynamics (MD) study, however, silver particles were generated by growth from the gas phase.13 In this paper we present the results of computer simulations of free gold nanoparticles. Small gold particles have been a particularly popular subject for investigation for well over a century,14,15 partly because of their stability in colloidal sols, which makes them convenient objects for testing theories of light scattering,16 partly because of their use in medicine,17 and partly because they are interesting in their own right. Various methods have been used to prepare gold clusters in the laboratory. These include condensation of a metallic vapor on different substrates,18-20 nucleation in an inert-gas atmosphere,21 and, more recently, as cores of particles passivated by selfassembled monolayers.22 Different experimental methods under different conditions gave a rich diversity of structures. The main types identified have been cuboctahedra (CO) with facecentered-cubic (FCC) internal structures, twinned FCC (containing one or several parallel twin planes),20 twinned hexagonal close packed (HCP),23 icosahedral and truncated icosahedral (Ih),24,25 and truncated decahedral (Dh). Computational studies have also ranged over a variety of approaches. Recently, Cleveland et al.26,27 have used an embedded atom (EAM) potential to determine the optimal structure
10.1021/jp0109426 CCC: $20.00 © 2001 American Chemical Society Published on Web 09/01/2001
11606 J. Phys. Chem. B, Vol. 105, No. 47, 2001 of nanometer-size Au clusters comprised of 100-1000 atoms (1.4-3.0 nm in diameter). A close competition among several structural types was found across the entire range of sizes. The predominant types were truncated decahedral, truncated octahedral (TO+), and its symmetrically twin-faulted variant (tTO+). All other configurations, including Ih, which is the lowest energy configuration for Ni or Cu clusters in the same size range,28 are much less stable. A similar result was obtained by Baletto et al.29 using quenched molecular dynamics simulations with a potential function developed on the basis of the secondmoment approximation to the tight-binding model. On the other hand, it was found30 that during heating gold nanoparticles initially possessing their stable, low-temperature structures underwent a solid-state transformation to an icosahedral structure before melting. In the present paper we report the results of molecular dynamics simulations of melting and freezing of gold nanoclusters of three different sizes. To better understand the process of crystallization, while at the same time acquiring sets to enable a determination of nucleation rates, we performed simulations on sets of molten clusters of the same size. Even though the individual clusters in each of the sets were generated and frozen under identical conditions, several different final structures were obtained. These will be discussed in the following.
Chushak and Bartell TABLE 1: Bond-Order Parameters for Different Local Structures structure
q4
q6
q8
FCC HCP Ih Dh
0.191 0.097 0.0 0.053
0.575 0.485 0.663 0.430
0.404 0.317 0.0 0.139
Structure Analysis Several criteria can be used to identify the local atomic arrangement in molecular dynamics simulations. We employ the bond-order parameter method34,35 to analyze the clusters’ structures as well as to distinguish between atoms with a liquidlike environment and atoms with one that is solidlike. To avoid tedious circumlocution we will refer to such atoms, respectively, as “liquid atoms” and “solid atoms.” We define “bonds” based merely on a geometrical relationship in terms of unit vectors rij connecting an atom i with its neighboring atoms j that are within a given radius rcut of i. The orientation of a bond rij with respect to some reference system is specified by the spherical harmonics functions Ylm(rij) ) Ylm(θij;φij), where θij and φij are the polar and azimuthal angles of vector rij in this reference frame. Only even-l spherical harmonics, which are invariant under inversion, were considered. The local order around an atom i is determined by averaging over all bonds with its neighbors Nnb(i)
Details of Computer Simulation Molecular dynamics simulations were performed on nanoclusters comprised of 459, 1157, and 3943 gold atoms represented by an EAM interaction potential.31 This is an n-body potential proven to reproduce both bulk and surface properties of transition metals.32 In the present work, simulations were carried out at constant temperature by rescaling velocities to bring the system to the desired temperature with a tolerance of (5 K. Equations of motion were integrated by a fifth-order Gear predictor-corrector algorithm with a time step ∆t ) 2.8 fs for 459 and 1157 atoms clusters and ∆t ) 4.3 fs for the largest one.33 For each of the cluster sizes, an initial, approximately spherical cluster was constructed to possess an FCC structure. Each cluster was heated from 300 K with temperature steps of 50 K and equilibrated at each temperature for ∼200-250 ps. In the vicinity of the melting transition the temperature step was reduced to 10 K. Clusters were heated 200 K above their melting temperatures to ensure complete melting. The final configuration was additionally equilibrated, generating along the way 20 saved configurations for 459 and 1157 atom clusters and 12 saved configurations for a 3943-atom cluster to serve as starting configurations for cooling runs. These clusters were cooled to 750 K in decrements of 50 K, spending ∼5 ps at each temperature, corresponding to a cooling rate of 15 × 1012 K/s. All configurations generated for production runs remained liquid, even though deeply supercooled, as determined by structure analyses. Production runs of 1 ns duration were carried out at 700 and 720 K for 459-atom clusters, at 700, 720, and 740 K for 1157atom clusters, and at 720 and 740 K for 3943-atom clusters. The frozen clusters which materialized were cooled to 300 K at a cooling rate of 3 × 1011 K/s for further mild annealing. Phase transitions were recognized by the evolution of configuration energy and by the structure analyses described in the next section.
qlm(i) )
Nnb(i)
1
Ylm(rij) ∑ j)1
Nnb(i)
(1)
To avoid a dependence of the local order parameter on the choice of reference system, a second-order invariant can be constructed:34
ql(i) )
(
∑ |qlm(i)|2 2l + 1 m) -l l
4π
)
1/2
(2)
The global order parameter Ql can be obtained by averaging qlm(i) over all N atoms in a cluster:
Ql )
(
4π
l
∑
2l + 1 m)-l
| |)
1/2
Q h lm
2
(3)
where N
Q h lm )
Nnb(i) qlm(i) ∑ i)1 N
(4)
Nnb(i) ∑ i)1 In an isotropic liquid state, the global bond-order parameter goes to zero in the thermodynamical limit for all values of l > 0. In a solid cluster where the bond orientations are correlated throughout the cluster, the value of the order parameter Ql is nonzero and depends on the crystalline structure. For symmetry reasons, the first nonzero value of ql(i) (other than the constant value for l ) 0) occurs at l ) 4 for atoms with cubic or decahedral local symmetry and at l ) 6 in aggregates with icosahedral symmetry. The values of bond-order parameters for different structures are listed in Table 1. At finite temperatures the crystal structure is distorted by thermal vibrations of atoms around their equilibrium positions. As a result, a given structure
Melting and Freezing of Gold Nanoclusters
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11607
Figure 1. Temperature dependence of the global order parameter Q6 and configurational energy per atom U for gold nanoparticles during heating stages. Open circles, 459-atom clusters; solid circles, 1157atom clusters; open triangles, 3943-atom clusters. The spike in Q6 and plateau in U for the smallest particle just before melting are real. They correspond to a change in structure.
will be characterized by a distribution of bond-order parameters rather than a single idealized value. We have used the distributions of q4(i), q6(i), and q8(i) signatures for different configurations as the criterion to identify the solidlike particles and the global structure of clusters produced in simulations. Once atoms in supercooled liquids have been identified as solidlike, they are considered to be members of embryos, the precursors of nuclei. If an atom so identified is a neighbor of another solidlike atom, it is assigned to the same embryo. The cutoff distance for neighbors was chosen to be rcut ) 3.5 Å, the position of the first minimum in the pair correlation function of FCC gold. The structure of a cluster is also reflected in its diffraction pattern.36 We have used the Debye scattering equation37 to calculate the X-ray diffraction elastically scattered intensity I(s) corresponding to a rotationally averaged ensemble as a function of the scattering variable s ) 2 sin(Θ)/λ:
I(s) )f 2(s)
∑ i,j
sin(2πsrij) 2πsrij
(5)
where rij is the distance between atoms i and j in the cluster, f(s) is the atomic scattering factor, and Θ is the Bragg angle. Results Melting. The temperature of a melting transition can be identified in several ways. We employ the variations of potential energy U and global order parameter Q6 with temperature. They are shown in Figure 1 for three clusters studied. For large clusters, containing 1157 and 3943 atoms, both potential energy and order parameter vary smoothly with temperature and exhibit a sharp jump near the transition. On the other hand, clusters composed of 459 atoms undergo structural reorganizations upon heating. The global order parameter significantly drops well
Figure 2. Computed melting point as a function of N-1/3. The asterisk shows the bulk melting temperature for the EAM potential function. Solid line, Pawlow’s triple point theory (t ) 0); crosses, second-order corrections included for t ) 0; heavy dashed curve, Sambles’s liquid shell theory, t ) 5 Å; light dashed curve, second-order corrections included, t ) 5 Å.
before the melting transition. Furthermore, it executes a sudden spike indicative of a reorganization close to the cluster’s melting point. Small particles melt at temperatures Tm ) To + ∆T, depressed from the bulk melting point To by surface effects. Several theories have been developed to account for this behavior. In a thermodynamic approach to the problem, Pawlow38 considered the triple point equilibrium between a solid and a liquid particle of equal mass in contact with their vapor. Neglecting secondorder terms, he derived the following relation for quasi-spherical solid particles of radius rs:
∆T ) -
[ ()]
Fs 2V h sTo σ - σl Lrs s Fl
2/3
(6)
where |∆T| is the melting point depression, L is the molar latent heat of fusion, V h s is the molar volume of the solid, Fs and Fl are the mass densities of the solid and liquid phases, respectively, and σs and σl are the surface tensions of the corresponding condensed phases. In another phenomenological approach based on a treatment by Reiss and Wilson39 and developed by others,40,41 the solid particle is supposed to be surrounded by a thin liquid shell of thickness t at the still lower equilibrium melting point Tlsm. According to Sambles,41 the change in melting point, ∆T ) Tlsm - To for a particle of radius r is
∆T ) -
[
( )]
σl Fs 2V h sTo σsl - 1L r-t r Fl
(7)
where σsl represents the solid-liquid interfacial tension. Equation 7 is often used to estimate the solid-liquid interfacial tension, a quantity which is very difficult to measure directly, from the size dependence of melting temperature.41-43 The melting points obtained in MD simulations are plotted in Figure 2 as a function of N-1/3. In this graph are also shown the results expected for several theoretical models. For the sake of consistency, model calculations were performed with parameters appropriate to EAM potential for gold, rather than to experimental data for the real material (see Table 2). Pawlow’s first-order triple point theory, eq 6, implies a linear dependence of the melting temperature in a plot of Tm vs N-1/3.
11608 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Chushak and Bartell
TABLE 2: Thermodynamic and Physical Properties of Bulk Gold from EMA Simulations and Experiment quantity melting temp mass density liquid solid molar latent heat surf. tension liquid solid solid-liquid diffusion coeff (T ) 1500 K)
symbol
EAM potential
exptl
Tm (K)
109049
133659
Fl (kg/m3) Fs (kg/m3) L (J/mol)
17 28048 19 00048 10 59748
17 28059 18 40060 12 36259
σl (J/m2) -dσ/dT (mJ/m2 K) σs (J/m2) -dσ/dT (mJ/m2 K) σsl (J/m2) D (10-9 m2/s)
0.7448
1.1361 0.1461 1.4062 0.1462 0.2741
0.9048 0.13-0.15 4.0 (this work)
Such a dependence agrees well with the MD data only for large clusters. In adjusting the liquid shell model (eq 7) to agree with the plot, the solid-liquid interfacial tension and shell thickness were considered as unknown parameters. This procedure yielded the results σsl ) 0.145 ( 0.01 J/m2 and t ) 5.0 ( 0.5 Å. The thickness of the liquid shell derived is approximately twice the gold-gold internuclear distances in the metal and, hence, is not implausible in magnitude. The value of σsl obtained agrees reasonably well with the value of 0.16 J/m2 implied by Antonow’s rule44 σsl ) σs - σl, and with those inferred from several empirical relationships. Such empirical expressions to estimate the solid-liquid interfacial free energy include one developed by Tyson and Miller,45 who proposed that
σsl )R σs
or
σsl R ) σl 1 - R
(8)
where R ) 0.15 ( 0.03 for metallic systems. By employing these relationships, we obtained σsl in the range 0.11-0.16 J/m2. Another empirical estimation of σsl, proposed by Turnbull,46 can be written as
σsl ≈
kTL (V h 2NA)1/3
(9)
〈r2(t)〉 6t
TABLE 3: Distributions of Final Configurations Materializing Spontaneously during the Freezing of Gold Nanoclusters 459-atom clusters
1157-atom clusters
3943-atom clusters
structure
700 K
720 K
700 K
720 K
740 K
720 K
740 K
Ih Dh TO HCP ???
18 1 1
19 1
12 2 6
13 1 4 2
12 4 3
7 1 2 1 1
7 4 1
total
20
20
20
20
19
12
12
This coefficient is needed for the prefactor of the nucleation rate equation. Although the surface atoms have a higher diffusivity than the internal ones, we did not discriminate between surface and interior atoms in computing D. The results of our diffusion investigations are presented in Figure 3. Since the experimental value of the diffusion constant for bulk gold is not available, we compared our results for clusters with the data from MD simulations of bulk liquid gold.47 A good agreement between the two sets of data was observed. The calculated temperature dependence of the diffusion coefficient was fitted by the empirical Arrhenius functional form
D ) D0 exp(-Ea/RT)
where V h is the molar volume (not originally specified as of solid or liquid) and NA is Avogadro’s number. Using the reported value of the proportionality constant kT, namely 0.45 for metal systems, we obtained a value for σsl of 0.12 J/m2. However suggestive the above results are, they should be regarded with some skepticism for two reasons. For one thing, as will be discussed subsequently, the Pawlow and Sambles models themselves are based on equivocal assumptions. For another, the data upon which the empirical relations were calibrated are difficult to obtain accurately and, at least in some cases, depended upon inferences incorporating the equivocal assumptions just referred to. Freezing. The final liquid configurations were additionally equilibrated and, as mentioned above, 12-20 configurations were generated for the cooling runs. These clusters were cooled to 750 K and served as starting configurations for freezing. During the cooling runs a self-diffusion coefficient was calculated from the slope of the mean-square displacement curve r2(t) according to the Einstein equation
D)
Figure 3. Arrhenius plot of the diffusion coefficient for normal and supercooled liquid gold nanoclusters calculated from MD simulations. Solid circles, solid triangles, and open squares represent runs for 1157atom clusters, 3943-atom clusters, and bulk gold particles, respectively.
(10)
(11)
which yields an activation energy Ea ) 23.88 kJ/mol and a prefactor D0 ) 27.9 × 10-9 m2/s. Despite the fact that individual clusters for each size were frozen under identical conditions, they froze to different final structures. One of the clusters with 1157 atoms failed to freeze at 740 K during the 1 ns simulation run. In Table 3 we summarize the distributions of final configurations for gold clusters frozen at different temperatures. Clearly the icosahedral structure formed preferentially for all cluster sizes, but others also formed spontaneously. Figure 4 shows typical structures obtained during the freezing of supercooled clusters. The dark gray spheres represent gold atoms with an FCC local structure, and light gray spheres are atoms with an HCP environment. In the case of Ih or Dh structures HCP atoms are located in the twin boundary of FCC tetrahedra with (111) faces. Stacking faults in the FCC planes of TO clusters cause the planes identified as HCP. Correspondingly, the stacking faults of planes in the HCP structure lead to the appearance of planes with FCC local symmetry. Figure 5 presents the time evolutions of the number of solid atoms with different local environments during the freezing of 459-atom clusters at 700 K. Only core atoms (atoms that have 12 or more neighbors) with well-defined local
Melting and Freezing of Gold Nanoclusters
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11609
Figure 6. Time evolution of a global order parameter Q6 for 1157atom clusters of different final structures.
Figure 4. Typical structures materializing spontaneously during the solidification of supercooled gold clusters. Dark gray spheres represent gold atoms with an FCC local structure, light gray spheres are atoms in an HCP environment, and black atoms indicate sites with 5-fold (Dh) symmetry.
Figure 5. Time evolution of the number of solid atoms for a particular 459-atom cluster (a) that froze to an icosahedral structure and (b) for another that froze to a truncated octahedral structure. Solid lines represent the number of atoms in an environment with 5-fold local symmetry and dotted curves one with FCC local symmetry.
symmetry (FCC, Ih, or Dh) are presented on these plots. The evolution of a typical cluster that froze to the icosahedral structure is shown in Figure 5a. In this configuration only a central atom is in site with a truly icosahedral local symmetry. Other atoms in the Ih cluster are located on three different sites: internal atoms in the tetrahedral units that form an icosahedron have an FCC environment, atoms on the internal
twinning plane have an HCP local symmetry, and atoms that lie on the line connecting the central atom and one of its surface vertex atoms have a decahedral symmetry. Before nucleation occurred, supercooled clusters contained dozens of solid atoms, mostly with a 5-fold symmetry (either Ih or Dh). At ∼90 ps an embryo with FCC and HCP atoms appeared, but it quickly dissolved back. In the cluster under discussion another embryo appeared at ∼120 ps and started to grow. The number of atoms with a 5-fold symmetry rapidly decreased to ∼20 and then increased to ∼30-35 atoms. The completed nth icosahedral shell contains 12 atoms on vertexes, 30(n - 1) atoms on edges, and 10(n2 - 3n + 2) atoms in the faces of the shell.50 The cluster with three completed Ih shells then contains 37 atoms with 5-fold symmetry, 20 FCC atoms, and 90 HCP atoms. The final configuration of a cluster, presented in Figure 5a, contains nearly the same number of different types of atoms, although the number of HCP atoms is smaller. The freezing picture is completely different in Figure 5b, where the time evolutions of different solid atoms in a cluster which developed into a truncated octahedral final structure are presented. When the crystal started to grow, the number of atoms with a 5-fold symmetry rapidly decreased and finally disappeared, while the majority of atoms acquired an FCC symmetry. Small numbers of atoms with an HCP local symmetry are mostly located in the stacking-fault planes of the FCC cluster. Different cluster structures have different degrees of crystal order or crystallinity. Figure 6 shows the time evolutions of the global order parameter Q6 for 1157-atom clusters with different final structures. In the TO cluster with an FCC structure all atoms have the same local bond orientations. As a result, the global order parameter Q6 (eq 3) tends to its ideal FCC value. A truncated decahedron is composed of five slightly distorted truncated tetrahedra that share a common edge. In the Dh structure atoms have the same orientations of nearest bonds within the tetrahedral units but different orientations between the tetrahedra. Therefore, Q6 is lower in such a cluster compared with the FCC structure. An icosahedron may be considered as being composed of 20 distorted tetrahedra joined at the center of the cluster with different orientations of units. It results in a very low degree of crystallinity in Ih clusters, as one can see in Figure 6. The X-ray diffraction (XRD) powder patterns for various structures as a function of cluster size are presented in Figure 7. It is apparent that sizes and structures of nanoparticles have a significant impact on the diffraction pattern. The diffraction pattern of TO 3943-atom clusters has well-defined peaks, which can be indexed as Bragg reflections from planes in an FCC lattice. Decreasing particle size causes a broadening of peaks
11610 J. Phys. Chem. B, Vol. 105, No. 47, 2001
Chushak and Bartell
Figure 8. Plots of ln(Nn/N0) vs time of nucleation for the sets of data given in Table 4.
ation rate J can be derived from the fraction of unfrozen particles Nn/N0 at time tn by comparing MD results with the first-order rate law
Nn ) exp[-JVc(tn - t0)] N0
Figure 7. Calculated X-ray diffraction powder patterns for clusters of different sizes and structures at T ) 300 K.
TABLE 4: Nucleation Times (in ps), and Final Configurations for 1157-Atom Clusters Frozen at Different Temperatures run no.
700 K
720 K
740 K
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
88 (TO) 50 (Ih) 29 (TO) 51 (TO) 109 (Ih) 48 (Ih) 45 (TO) 139 (Ih) 37 (Ih) 38 (Dh) 110 (Ih) 39 (TO) 105 (Ih) 72 (TO) 85 (Ih) 25 (Ih) 68 (Ih) 57 (Ih) 38 (Ih) 94 (Ih)
81 (TO) 225 (Ih) 80 (Ih) 108 (HCP) 106 (Ih) 75 (Ih) 87 (Ih) 158 (Ih) 137 (Ih) 30 (TO) 72 (Ih) 256 (Ih) 166 (Ih) 163 (Dh) 149 (Ih) 131 (TO) 70 (Ih) 38 (Ih) 94 (HCP) 180 (TO)
10 (Dh) liquid 611 (Ih) 145 (Dh) 275 (Dh) 60 (Ih) 150 (Ih) 50 (Ih) 190 (Dh) 850 (Ih) 85 (Ih) 325 (TO) 65 (Ih) 155 (Ih) 277 (TO) 62 (Ih) 115 (Ih) 162 (Ih) 288 (TO) 320 (Ih)
in the diffraction profiles, leading to a distinct asymmetry of the main peak. XRD reflections from Dh and especially Ih clusters are weaker and much more diffuse because domains with ordered structures are smaller in these structures. It is known that the freezing of supercooled liquids is an activated process which occurs by crystal nucleation and growth. Nucleation times for freezing tn can be readily recognized from plots such as those in Figure 6. In Table 4 we list tn for 1157atom clusters at three different temperatures, together with the indication of the final configuration. As can be seen from the table, there is no discernible correlation between the nucleation time and final cluster structure. Therefore, freezing was counted whenever a cluster froze, regardless of its final structure, a criterion usually adopted in nucleation experiments. The nucle-
(12)
where N0 is the number of clusters, Vc is the volume of the cluster, and t0 is the nucleation time lag (here assumed to be constant at a given temperature). Since the time tn is that at which the nth nucleation event has occurred, Nn is taken to be51 N0 - n + 1. Figure 8 presents plots of ln(Nn/N0) vs the time of nucleation for the sets of data given in Table 4. From the slope of these plots, the quantity JVc was obtained, and from the intercept, the time lag t0 was estimated. In the theory of homogeneous nucleation, the nucleation rate is expressed as
J ) A exp(-∆G*/kBT)
(13)
where A is a prefactor and ∆G* is the free energy of formation of a critical nucleus. Two different approaches to the crystal nucleation were considered in present analysis. One is the classical (capillary) nucleation theory (CNT) developed by Turnbull and colleagues,52,53 and another is Granasy’s diffuse interface theory (DIT)54 that explicitly takes into account a thickness of interface between the solid and liquid phases. The unknown quantity in the CNT is the kinetic parameter σsl, supposed to represent the liquid-solid interfacial free energy per unit area, while in DIT the unknown parameter is δ, considered to express the distance between the dividing surfaces for enthalpy and for entropy in the interface between the liquid and solid phases. In Figure 9 we present the temperature dependence of the nucleation rate for different sizes of gold nanoclusters derived from MD simulations in comparison with the results of the CNT and DIT theoretical expressions. We took the CNT interfacial free energy σsl and the DIT distance δ to be independent of temperature for the purposes of the figure. These two parameters were estimated by adjusting them to reproduce the nucleation rate for the 1157-atom cluster at 720 K obtained from MD simulations. Values so determined were σsl ) 0.084 J/m2 and δ ) 1.2 Å. The value of σsl obtained from the nucleation experiment is appreciably lower than the corresponding parameter obtained by applying the Sambles equation for the size dependence of the melting temperature.
Melting and Freezing of Gold Nanoclusters
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11611
Figure 10. Stages of freezing in the model of Reiss, Mirabel, and Whetten. The figure also identifies the symbols used in the treatment of Sambles liquid shell model.
and the heats of vaporization and sublimation and surface tensions to be linear functions of the temperature, so that
V h i(To+∆T, p°+∆p) ) V h i°[1 + 3Ri(∆T) - κi(∆p)] Figure 9. Temperature dependence of the nucleation rate for different sizes of gold nanoclusters derived from MD simulations in comparison with the results of the CNT (solid line) and DIT (dashed line) theoretical calculations. Open circles, diamonds, and solid circles represent runs for 459-, 1157-, and 3943-atom clusters, respectively.
Discussion Melting. It is worthwhile to examine several popular models devised to treat the size dependence of the melting points of nanoparticles. These models include one presented in the classic paper by Buffat and Borel42 and those proposed by Pawlow (eq 6) and by Sambles (eq 7) which were applied in the foregoing. In all of these the underlying idea was that, at the melting temperature of a nanoparticle, the chemical potential of the matter in the solid particle is the same as that in a melted particle of the same mass. It was also assumed that the particle is spherical, or nearly so. Several years ago Reiss, Mirabel, and Whetten55 (RMW) showed that this condition of equality of chemical potential is incorrect, and a later paper2 added additional caveats. Nevertheless, it is instructive to develop the popular relations which apply if the equalization condition is adopted. It is far simpler to derive these relations by taking advantage of the intimate connection between vapor pressure and chemical potential than to use the approach of Buffat and Borel. If the chemical potentials of the solid and liquid particles are the same, so are the vapor pressures. Therefore, the Clapeyron and Kelvin equations suffice to yield the general results. We take it for granted that the vapor pressures of the metal clusters involved are so low that the molar volumes of the vapors are enormous compared with those of the condensed phases, and that the vapors are adequately treated as ideal. Therefore, for the bulk condensed phases, the Clapeyron equation reduces to
ln[pi (T)] ) ln[pi (To)] + ∆
∆
∫T [Li(T)/RT ] dT T
2
(14)
o
where the vapor pressure pi∆ refers the vapor pressure that the condensed phase i experiences at the standard pressure p°, and Li(T) is the corresponding temperature-dependent molar heat of vaporization or sublimation. Vapor pressures of the particles are augmented by the Laplace pressures brought to bear on them, as expressed by the Kelvin equation:
ln[pi(T,ri)] ) ln[pi∆(T)] +
∫p°p°+p
L
(V h i/RT) dP
(15)
where ri is the radius of the particle and pL is the applicable Laplace pressure. Following the prior treatments, we take the molar volumes to be linear functions of temperature and pressure
(16)
Li(To+∆T) ) Li° + Cpi∆T
(17)
σi(To+∆T) ) σi° + (dσi/dT)∆T
(18)
where Ri is the linear coefficient of expansion, κi is the compressibility of phase i, and V h i°, Li°, and σi° represent values at To. The melting point depression |∆T| can be calculated from the above equations by setting the vapor pressures of the liquid and solid particles at the melting point
Tm ) To + ∆T
(19)
to be equal, and solving for ∆T. The procedure for doing this is described in the Appendix. If terms higher than linear in ∆T and ∆p are neglected, the solution of Buffat and Borel is obtained. If the change in heat capacity from solid to liquid is neglected as well as the coefficients of thermal expansion and compressibility and temperature dependences of the surface tensions, the Pawlow result, eq 6, is obtained. If these simplifications are retained but it is assumed, in addition, that a liquid shell of thickness t surrounds the solid particle at the (still lower) melting temperature, Tlsm, then the liquid shell result of Sambles, eq 7, is obtained, very nearly. In all of these models, bulklike properties are attributed to the very small particles, which is one reason, but not the only reason, for their quantitative deficiencies. A schematic sketch of the systems in the liquid shell model is given in Figure 10, where the various radii are defined. In the liquid shell model, the Laplace pressure exerted upon the solid core is then
pL ) 2σl/rlsm + 2σsl/rs
(20)
making the melting point change
∆T ) -
(
[
)]
Fs rlsm 2V h sTo σsl σl 1+ L h fus rlsm - t rlsm Fl r l
(21)
the same expression published by Sambles provided that rlsm ) rl, an equality that is approximately but not completely accurate. Actually, Sambles was somewhat vague about the meaning of the radius r in eq 7 that corresponds to rlsm in eq 21, and therefore it is not possible to compare expressions rigorously. The ratio rlsm/rl in eq 21 is
() ( )
rlsm Fl ≈ rl Fs
1/3
+
Fs t 1rl Fl
Therefore, for small t, eq 21 becomes
(22)
11612 J. Phys. Chem. B, Vol. 105, No. 47, 2001
∆T ) -
[
{ ( ) }]
2V h sTo σsl Fs σl 1+ L h fus rlsm - t rlsm Fl
Chushak and Bartell
2/3
(23)
When the interfacial free energy in this expression, σsl, is replaced by the Antonow relation,44 which relation was shown by Reiss and Wilson39 to hold for cases corresponding to the present treatment, it can be seen that the corrected Sambles expression reduces identically to the Pawlow expression when thickness t goes to zero. In Figure 2 the melting points of the gold particles are compared with those calculated not only by the approximations of Pawlow and of Sambles but also by those including the second-order corrections of Buffat and Borel, while taking the thickness parameter t to be 5 Å. In the original papers by Buffat and Borel and by Castro, Reifenberger, Choi, and Andres,43 in which the liquid film thickness was taken to be zero, the secondorder corrections began to make observed melting points depart considerably from the Pawlow curve as the particle size is decreased, particularly because of the value adopted for the slope -dσs/dT for the solid. This value now seems unrealistically high. When the parameters (Table 2) selected to be as consistent as possible with the potential function of the present simulations were used, the second-order corrections became much smaller than in the papers by Buffat and Borel and by Castro et al. It is apparent that the fit with MD simulations is improved if the thickness of the liquid shell at melting is assumed to be about twice the gold-gold interatomic distance. Such a value is not implausible, but as suggested in the foregoing, it is not based on a sound foundation. None of the above treatments of melting should be considered to be rigorous. To support this assertion, we refer to a treatment by Reiss, Mirabel, and Whetten (RMW). These authors introduced a model (Figure 10) in which the free energy G(rs) of a spherical particle with a solid core of radius rs and liquid shell of thickness, t, is calculated as a function of the core radius. It is well-known that, in MD simulations of clusters, as solid clusters are heated, they melt on the outside well before the core melts. Therefore, a treatment in which a cluster exists partially as solid and partially as liquid is quite appropriate. In the RMW approach, the so-called thermodynamic condition for melting is that the chemical potentials the contents of the liquid and the solid particles are equal at Tm ) To′, in the notation of RMW. This condition corresponds to a maximum in the G(rs) curve. The maximum is at t ) 0 and the temperature To′ is essentially the Pawlow result (in view of simplifications in the RMW model). This result implies that the equilibrium at To′ is an unstable equilibrium so that, if the molten system were cooled to To′ and freezing began, any disturbance of the system would melt the system. This can be seen in Figure 11 illustrating the behavior of G(rs) at various temperatures from above To′ to far below To′. It is true that the RMW model is simplified to a case where the densities of the solid and liquid are the same, and the temperature dependences of volumes, surface tensions, and the heat of fusion are neglected, but these simplifications should not affect the essential implications of the results. It should be noted that the liquid shell model of melting also yields a maximum in G(rs) at thickness t, again implying an unstable equilibrium. It is clear that the liquid shell model is better at reproducing the MD results than the plain “thermodynamic model”, because it introduces an additional parameter helping to account for the accelerating decrease in melting point with decreasing particle size. Still, even though it is well-known that small particles begin to melt at the surface before the core melts, the model is an ad hoc model, not one based on fundamental
Figure 11. Variation of Gibbs free energy of a nanoparticle (in arbitrary units) with isothermal freezing, according to the Reiss, Mirabel, and Whetten model as the radius of the solid core grows. The curves going from top to bottom correspond to increasing degrees of supercooling. The solid curve is for the temperature To′, the so-called thermodynamic melting point at which the chemical potentials for the contents of the completely liquid and fully frozen nanoparticles (clusters) are the same. A more realistic melting temperature would be that estimated as follows.2 If the curves were weighted by the fraction of frozen cluster and the Boltzmann distributions of clusters of different G(rs) were calculated, the temperature at which the population of mostly liquid clusters equaled that of the mostly solid clusters could be considered the melting point. This temperature would be lower than To′ but not as low as the RMW temperature To′′ at which G(rs) for the clusters at rs ) 0 equals that for the clusters at rs ) rmax.
principles. Moreover, even though the RMW model is qualitatively correct in implying a lower melting point than the Pawlow model and in forecasting a range of temperatures over which melting takes place for a given particle size, it is quantitatively very poor in reproducing results of MD simulations, particularly for small molecular clusters. The reasons for its deficiencies include the assumption of bulklike thermodynamic properties for the liquid and solid portions of the particles even though those regions may contain comparatively few molecules. In addition, the model totally neglects possible effects of interfacial thicknesses. The liquid shell model crudely simulates a surface transition layer to the vapor phase, but it totally neglects the interfacial thickness between liquid and solid phases of the melting particle. What can be concluded from this discussion? First, that interfacial properties inferred by applying the model equations to experimental results or results of MD simulations can, at best, be accepted with considerable reservation. Next, interfacial properties derived from the kinetics of freezing need not be close to those derived from the thermodynamics of melting, even though both are based on a simplified capillary model. A somewhat similar discussion of weaknesses of the popular nucleation theories is given in ref 35. Currently, the best theoretical accounts of melting and freezing are given by density functional computations. Unfortunately, these calculations are complex and depend on specific properties of particular systems and hence are far more difficult to apply to a general case than the capillary-based models discussed above. A realistic model readily applicable to general systems would be a very helpful advance. Freezing. Intensive computer investigations26,27,56 which sought the optimal structures for gold clusters have indicated that the icosahedral structure is energetically noncompetitive
Melting and Freezing of Gold Nanoclusters
J. Phys. Chem. B, Vol. 105, No. 47, 2001 11613
even in the very small size range. As the size of a cluster increased, this structure was calculated to become less and less favorable due to accumulated strain energy, especially at the center of a cluster. On the other hand, icosahedral gold nanoparticles generated in the laboratory by various methods have been observed by high resolution electron microscopy.19,20,24,25 Similar results were obtained for free copper and silver clusters produced in an inert-gas aggregation source and flowing in a molecular beam.57,58 In our own computer experiments on gold clusters, most of the clusters of three different sizes studied froze to an icosahedral structure. One could argue that it is not an equilibrium structure due to the limited time of computer simulations. We performed further annealing of clusters by relatively slow cooling to 300 K. This treatment did not induce icosahedral clusters to transform to any of the other possible low-energy structures. Another striking result of our simulations was the formation of nanoparticles with a predominantly HCP structure. Such a structure was also observed in experimental investigations.23 This is of interest because the HCP structure had not been considered to be stable. We have observed a transition from HCP to the FCC structure in 3949-atom clusters. Nevertheless, two clusters with 1157 atoms remained in a principally HCP arrangement even during cooling, indicating the presence of a significant energy barrier for the HCP-FCC transformation. Clearly, the kinetics of nucleation and growth in clusters play a key role in determining the structure of clusters. Different structures were obtained under identical conditions due to the stochastic nature of nucleation. Acknowledgment. This research was supported by a grant from the National Science Foundation. Appendix To make the present treatment for calculating the melting points of small particles fully explicit, we enumerate the steps involved. The second-order terms of Buffat and Borel and higher-order terms are generated in the following way. Equations 14 and 16-18 are incorporated into eq 15, and then the vapor pressures pl(To+∆T, p°+∆pl, rl) and ps(To+∆T, p°+∆ps, rs) are equated. Upon rearrangement is obtained
∫TT +∆T o
o
(
)
Lvap - Lsub 2
RT
dT )
∫p°p°+∆p
s
( )
V hs dp RT
∫p°p°+∆p
l
( )
V hl dp (24) RT
where the Laplace pressures are
∆pl ) 2[σl° + (dσl/dT)∆T]/[rl°(1 + Rl∆T)]
(25)
for a liquid particle, and
∆ps ) 2[σl° + (dσl/dT)∆T]/rlsm + 2[σsl + (dσsl/dT)]/[rs°(1 + Rs∆T)] (26a) with rlsm ) rs°(1 + Rs∆T) + t(1 + Rl∆T) in the case of the liquid shell model, or
∆ps ) 2[σs° + (dσs/dT)∆T]/[rs°(1 + Rs∆T)]
(26b)
for a fully frozen particle. When applying the liquid shell model for a given rlsm and shell thickness t, in order to obtain the corresponding rs, it is necessary to solve a cubic equation derived by equating the number of atoms in the liquid particle of radius
rl to the number of atoms in the partially frozen particle. The left-hand side (LHS) of eq 24 becomes
∫T
To+∆T o
(
)
-[Lfus° + ∆Cpfus(T - To)
-
RT2 Lfus° - ∆CpfusTo
(
RTo(To + ∆T)
)
∆T -
dT )
(
)
∆Cpfus ∆T ln 1 + (27) R To
and the right-hand side (RHS) becomes
{V h s°[(1 + 3Rs∆T)∆ps - (1/2)κs∆ps2] V h l°[(1 + 3Rl∆T)∆pl - (1/2)κl∆pl2]}/[R(T + ∆T)] (28) Solving for the ∆T that makes the LHS equal the RHS gives the desired solution. Although the above solution is not readily expressed in a closed form, it can be obtained by a simple computer program. References and Notes (1) Bartell, L. S.; Dibble, T. S. J. Phys. Chem. 1991, 95, 1159. (2) Bartell, L. S.; Chen, J. J. Phys. Chem. 1992, 96, 8801. (3) Huang, J.; Bartell, L. S. J. Phys. Chem. 1994, 98, 4543. (4) Torchet, G.; de Feraudy, M.-F.; Boutin, A.; Fuchs, A. H. J. Chem. Phys. 1996, 105, 367. (5) Valente, E. J.; Bartell, L. S. J. Chem. Phys. 1983, 76, 2683. (6) Dibble, T. S.; Bartell, L. S. J. Phys. Chem. 1992, 96, 8603. (7) Chushak, Y. G.; Bartell, L. S. J. Phys. Chem. B 1999, 103, 11196. (8) Torchet, G.; de Feraudy, M.-F.; Raoult, B.; Farges, J.; Fuchs, A. H.; Pawley, G. S. J. Chem. Phys. 1990, 92, 6768. (9) Torchet, G.; de Feraudy, M.-F.; Raoult, B. J. Chem. Phys. 1995, 103, 3671. (10) Doye, J. P. K.; Wales, D. J.; Berry, R. S. J. Chem. Phys. 1995, 103, 4234. (11) Bartell, L. S. Unpublished. (12) Marks, L. D. Rep. Prog. Phys. 1994, 57, 603. (13) Baletto, F.; Mottet, C.; Ferrando, R. Phys. ReV. Lett. 2000, 84, 5544. (14) Faraday, M. Philos. Trans. 1857, 147, 152. (15) Zsigmondy, R. Lieb. Ann. 1898, 301, 30. (16) Steubing, W. Ann. Phys. 1907, 24, 1; Ann. Phys. 1908, 26, 329. (17) Weiser, H. B. Colloid Chemistry; Wiley: New York, 1939; p 251. (18) Ino, S. J. Phys. Soc. Jpn. 1966, 21, 346. (19) Allpress, J. G.; Sanders, J. V. Surf. Sci. 1967, 7, 1. (20) Buffat, P.-A.; Flueli, M.; Spycher, R.; Stadelmann, P.; Borel, J. P. Faraday Discuss. 1991, 92, 173. (21) Kimoto, K.; Nishida, I. Jpn. J. Appl. Phys. 1967, 6, 1047. (22) Whetten, R. L.; Khoury, J. T.; Alvarez, M. M.; Murthy, S.; Vezmar, I.; Wang, Z. L.; Stephens, P. W.; Cleveland, C. L.; Luedtke, W. D.; Landman, U. AdV. Mater. 1996, 5, 428. (23) Jose-Yacaman, M.; Herrera, R.; Gomez, A. G.; Tehuacanero, S.; Schabes-Retchkiman, P. Surf. Sci. 1990, 237, 248. (24) Komoda, T. J. Appl. Phys. 1968, 7, 27. (25) Ascencio, J. A.; Perez, M.; Jose-Yacaman, M. Surf. Sci. 2000, 447, 73. (26) Cleveland, C. L.; Landman, U.; Schaaff, T. S.; Shafigullin, M. N.; Stephens, P. W.; Whetten, R. L. Phys. ReV. Lett. 1997, 79, 1873. (27) Cleveland, C. L.; Landman, U.; Shafigullin, M. N.; Stephens, P. W.; Whetten, R. L. Z. Phys. D 1997, 40, 503. (28) Cleveland, C. L.; Landman, U. J. Chem. Phys. 1991, 94, 7376. (29) Baletto, F.; Mottet, C.; Ferrando, R. Phys. ReV. Lett. 2000, 84, 5544. (30) Cleveland, C. L.; Luedtke, W. D.; Landman, U. Phys. ReV. B 1999, 60, 5065. (31) Foiles, S. M.; Baskes, M. I.; Daw, M. S. Phys. ReV. B 1986, 33, 7983. (32) Daw, M. S.; Foiles, S. M.; Baskes, M. I. Mater. Sci. Rep. 1993, 9, 251. (33) The Classical Molecular Dynamics code (ALCMD) developed at Ames Laboratory of the U.S. Department of Energy has been used in the present simulations. (34) Steinhard, P. J.; Nelson, D. R.; Ronchetti, M. Phys. ReV. B 1983, 28, 784. (35) Chushak, Y.; Bartell, L. S. J. Phys. Chem. A 2000, 104, 9328. (36) Hall, B. D. J. Appl. Phys. 2000, 87, 1666. (37) Warren, B. E. X-ray Diffraction; Dover: New York, 1990. (38) Pawlow, P. Z. Phys. Chem. 1909, 65, 1; Z. Phys. Chem. 1909, 65, 545.
11614 J. Phys. Chem. B, Vol. 105, No. 47, 2001 (39) Reiss, H.; Wilson, I. B. J. Colloid Sci. 1948, 3, 551. (40) Hanszen, K.-J. Z. Phys. 1960, 157, 523. (41) Sambles, J. R. Proc. R. Soc. London, A 1971, 324, 339. (42) Buffat, P.; Borel, J.-P. Phys. ReV. A 1976, 13, 2287. (43) Castro, T.; Reifenberger, R.; Choi, E.; Andres, R. P. Phys. ReV. B 1990, 42, 8548. (44) Antonow, G. N. J. Chim. Phys. 1907, 5, 372. (45) Tyson, W. R.; Miller, W. A. Surf. Sci. 1977, 62, 267. (46) Turnbull D. J. Appl. Phys. 1950, 21, 1022. (47) Bogicevic, A.; Hansen, L. B.; Lundqvist, B. I. Phys. ReV. E 1997, 55, 5535. (48) Lewis, L. J.; Jensen, P.; Barrat, J.-L. Phys. ReV. B 1997, 56, 2248. (49) Foiles, S. M.; Adams, J. B. Phys. ReV. B 1989, 40, 5909. (50) Mackay, A. L. Acta Crystallogr. 1962, 15, 916. (51) Chushak, Y.; Santikary, P.; Bartell, L. S. J. Phys. Chem. A 1999, 103, 5636. (52) Turnbull, D.; Fisher, J. C. J. Chem. Phys 1949, 17, 71.
Chushak and Bartell (53) Turnbull, D. J. Chem. Phys. 1952, 20, 411. (54) Granasy, L. J. Non-Cryst. Solids 1994, 162, 301; Mater. Sci. Eng. 1994, A178, 121; J. Phys. Chem. 1995, 99, 14183. (55) Reiss, H.; Mirabel, P.; Whetten, R. L. J. Phys. Chem. 1988, 92, 7241. (56) Michaelian, K.; Rendon, N.; Garzon, I. L. Phys. ReV. B 1999, 60, 2000. (57) Reinhard, D.; Hall, B. D.; Ugarte, D.; Monot, R. Phys. ReV. B 1997, 57, 7868. (58) Reinhard, D.; Hall, B. D.; Berthoud, P.; Valkealahti, S.; Monot, R. Phys. ReV. B 1998, 58, 4917. (59) Hultgren, R.; et al. Selected Values of the Thermodynamic Properties of the Elements; American Society for Metals: Metals Park, OH, 1973. (60) CRC Handbook of Chemistry and Physics, 64th ed.; Chemical Rubber Co.: Cleveland, 1983. (61) Miedema, A. R.; Boom, R. Z. Metallkd. 1978, 69, 183. (62) Miedema, A. R. Z. Metallkd. 1978, 69, 287.
