2074
J. Phys. Chem. C 2008, 112, 2074-2078
Long-Time Scale Molecular Dynamics Study of Diffusion Dynamics of Small Cu Clusters on Cu(111) Surface Jianyu Yang,†,‡ Wangyu Hu,*,† Jianfeng Tang,† and Maichang Xu‡ Department of Applied Physics, Hunan UniVersity, Changsha 410082, China, and Department of Maths and Physics, Hunan Institute of Engineering, Xiangtan 411104, China ReceiVed: June 19, 2007; In Final Form: NoVember 6, 2007
Diffusion dynamics of small two-dimensional atomic clusters Cun (with n ranging from 1 to 8) on Cu(111) surface have been studied using molecular dynamics in the temperature range from 300 to 800 K. The diffusion coefficients of these clusters are derived from mean square displacement of cluster’s mass-center, which is obtained by long simulation times (∼0.1 µs) and tracing of surface interstitial atoms. A drop in the diffusion migration energy for tetramer compared to the trimer or pentamer is found, and the migration energy for octamer is larger than that for the compact heptamer cluster. The diffusion prefactor of a heptamer is about 30 times larger than that for single atom diffusion. The characteristic temperature, derived from the Arrhenius plots, indicates that the heptamer starts to move at very low temperature (167 K).
1. Introduction A comprehensive understanding of the mechanisms and energetics by which small atomic clusters migrate across the single-crystal terraces is of both fundamental and technological importance. In addition to a providing physical insight into the details of adatom-surface and adatom-adatom interactions on surfaces, information on the mechanism and energetics of cluster motion can contribute to a better understanding for crystal and thin-films growth. For some time now, it has been recognized that clusters are able to migrate over crystal surfaces1-11 and many different mechanisms have been proposed to account for this mobility. For example, a leapfrog diffusion mechanism on FCC(110) surfaces,8 a rotation of a cluster, on FCC(100) surfaces,9 concerted rotation and translation for compact-cluster on FCC(111), dislocation mechanism of a cluster, on Fcc(111)4 have been found. Despite that, rather little quantitative information is available from theoretical studies about the diffusion properties of clusters, especially for the diffusion prefactor. A considerable amount of static quantitative calculation has been devoted to clusters self-diffusion.11,12 The diffusion coefficient is intrinsically a dynamic quantity. It is the prefactor which contains the true dynamical information, and the theoretical challenge is to understand this quantity properly within a dynamical theory.13 Therefore, the molecular dynamics (MD) simulation, which can involve all diffusion mechanisms, is necessary to describe the cluster diffusion. Also dynamic simulations at finite temperatures both reveal the diffusion mechanism and give the corresponding diffusion migration energy, and should also be the first choice.14 To our knowledge, only few dynamical studies were performed.15 In the present work, the long-time scale molecular dynamics simulations are used to investigate the diffusion of twodimensional (2D) clusters with Nc (Nc from 1 to 8) Cu atoms on Cu(111) surfaces in the temperature range from 300 to 800 * Corresponding author. E-mail:
[email protected]. † Hunan Unviersity. ‡ Hunan Institute of Engineering.
K. The diffusion prefactors and migration energies for these clusters are calculated. It is of interest to study cluster diffusion on the Cu(111) surface, because it proceeds in a particularly simple manner. In addition, FCC(111) surface is extremely stable in temperature, and it is therefore possible to study the phenomenon over a relatively wide range of temperature and to study the cluster size and temperature dependence of the diffusion prefactors and migration energies. This allows excellent quality diffusion data to be accumulated for a description of film growth. 2. Simulation Approach In the present MD simulation, the modified analytic embedded atom potential (MAEAM) are used.16,17 The potential provided a good description for many-body interaction, so the potential can distinguish the difference between atomic interaction in the bulk and that near the surface without any additional parameter. Our group constructed a MAEAM model for facecentered cubic (fcc) metals which has been successfully used for the study of both clusters and surface systems.17,18 In our MD calculations, the simulation box contains 21 atomic layers with 192 atoms each, arranged in an fcc lattice. Periodic boundary conditions are applied in the lateral directions, that is, parallel to the surface. Two free surfaces are obtained by fixing the dimensions of the supercell size to a value twice as large as the thickness of the crystal in the direction normal to the surface. In the present MD simulations, the six-value Gear predictor-corrector algorithm and Nose´ constant-temperature technique19 are employed, except for a series of bulk calculations in the constant temperature, constant pressure (NPT ensemble) controlled by Andersen method20 in order to determine the lattice constant at each simulated temperature, used for properly constructing a simulation box. A 2D atomic cluster with 1 to 8 atoms was deposited on the top of the surface layer of the slab. The stable configurations for all clusters are obtained by quenching molecular dynamics.21 Simulations are performed of 40-100 ns (depending on the temperature and the cluster size) in the temperature range from 300 to 800 K (depending on the cluster size), with the time step of 3 fs.
10.1021/jp074754c CCC: $40.75 © 2008 American Chemical Society Published on Web 01/24/2008
Small Cu Clusters on Cu(111) Surface
J. Phys. Chem. C, Vol. 112, No. 6, 2008 2075
In experiment, the diffusion coefficient for cluster can generally be obtained by measuring the mean-square displacement of the cluster’s mass center.22 In the present simulation, once the system is in thermo-dynamical equilibrium, that is, the system total energy should remain constant, the positions of interstitial atoms on surface for every time step are recorded for further analyses. In the scenario of two-dimensional randomwalk motion, the diffusion coefficient D is calculated using a long-time Einstein relation:
〈R2(t)〉 tf∞ 4t
D ) lim
(1)
where 〈R2(t)〉 is the mean square displacement of the mass-center of cluster and is given by
〈R2(t)〉 ) |∆x2(t) + ∆y2(t)|
(2)
In order to have reasonable statistics, the quantities |∆x2(t) + ∆y2(t)| are averaged over a number of time origins. Then the prefactor D0 and migration energy Em are obtained from the Arrhenius relation for the diffusivity:
( )
D ) D0 exp -
Em kT
Figure 1. Schematic picture of the structure of stable small Cu clusters on the Cu(111) surface. The gray balls show the surface atoms. All adatoms stay FCC sites. The lines joining the adatoms are drawn to show which atoms belong to the same cluster.
(3)
In order to ascertain that the potential used describes well the system under study, the diffusion of single Cu adatom on the Cu(111) surface is first considered. The migration energy obtained from the mean square displacement is 0.06 eV, which is in agreement with other theoretical calculations. In these computations the migration energy is found between 0.04 and 0.087 eV.12,15,27,28 And the present result is close to the existing experimental values of the Cu adatom diffusion on (111) surface, 0.10 eV.30 The calculated prefactor is 3.01 × 10-4cm2/s, which is consistent with other dynamical calculations. For example, the prefactor determined from the Arrhenius diagram of MD simulation has been found to be 5.2 × 10-4cm2/s when a Morse potential was used and 6.4 × 10-4cm2/s with a FS (FinnisSinclair) potential.28 With the EAM or EMT model, the calculated prefator was 3.0 × 10-4cm2/s or 4.6 × 10-4cm2/s,29 respectively. With the tight binding potential, Kallinteris et al.27 found that the self-diffusion for Cu adatom on (111) surface exhibited Arrhenius behavior at two distinct temperature regions, from which the diffusion prefactor was 2.0 × 10-4cm2/s and 7.0 × 10-4cm2/s, respectively. These comparisons indicate that the present EAM potential provides a rather good description of the dynamics properties of the Cu surface. 3. Results and Discussion On an FCC(111) surface, there exist two nonequivalent binding sites, which are called FCC and HCP sites. The FCC sites correspond to the ABC stacking of an fcc crystal. In contrast, adatoms located in HCP sites correspond to an ABAB stacking. In the present work, simulations using molecular dynamics show that the binding energies of Cun clusters are slightly higher (average by 0.007 eV) at an HCP site than at an FCC one. This small difference is less than the expected error of EAM and therefore cannot be taken seriously. This phenomenon has also been observed for Ir clusters on Ir(111) by a fieldion microscopy.23 In addition, the theoretical calculation with the effective medium theory (EMT) for Pt on Pt(111) has also obtained similar result.24 The initial configurations of Cun clusters are shown in Figure 1.
Figure 2. Binding energy as a function of cluster size.
The dissociation of a cluster is expected to depend on the size of the cluster as well as on the shape. According to nucleation theory, the density of islands on the surface during growth is dependent on the size of the critical nucleus, that is, the smallest adatoms configuration, which is stable. Therefore, knowledge of the stability of adatom clusters, as a function of size and shape is important for a detailed understanding of epitaxial growth and of the structure of the overlayers formed during growth.26 The binding energy of each n atom cluster is calculated and plotted in Figure 2. The binding energy is the energy reduction when the cluster breaks up into two parts with one being an adatom.25 The binding energies change clearly indicates that a cluster consisting of seven atoms is the most stable and followed by a trimer. Because of the structure stability, a heptamer is likely the critical nuclei of threedimensional growth of thin films on the Cu(111) surface. Whether it does act as the nuclear depends on its mobility, to be presented in the follow paragraphs. In the movement of particles on solid surfaces, for very short time the particles move ballistically and 〈R2(t)〉 ∼ t2. The longtime behavior of MSD gives us the diffusion coefficient. For the simulation of the Cu heptamer at 700 K, the MSD as a function of time is plotted in Figure 3. As shown in the figure, the plot is nearly linear and gives a reliable value for the diffusion coefficient. The Arrhenius plots of the diffusion coefficients of Cun are shown in Figure 4. It can be seen that the diffusion coefficients for all clusters follow well the
2076 J. Phys. Chem. C, Vol. 112, No. 6, 2008
Yang et al.
Figure 3. MSD of the mass-center of Cu heptamer vs time at 700 K.
Figure 4. Arrhenius plots of the diffusion coefficients for Cun clusters. The solid lines are Arrhenius fits to the molecular dynamics data.
TABLE 1: Diffusion Barriers (Em) and Corresponding Prefactor (D0) for Pt1 and Pt7 Clusters Pt1 Em (eV)
Pt7 D0 (cm2/2)
Em (eV)
D0 (cm2/2)
present 0.19 5.1 × 1.08 6.1 × 10-1 expt36 0.260 ( 0.003 2.0( × 1.4(1) × 10-3 1.17 ( 0.04 5.1( × 3.8(1) × 10-1 10-3
Arrhenius law. From these plots, one obtains the diffusion prefactors D0 and migration energies Em, associated to each cluster, which are plotted in Figures 5 and 6. In general, the migration energies derived from the static calculations are lower than that obtained from the MD simulation, because the static calculations may be based on a diffusion path with the lowest energy. The actual diffusion involves a number of diffusion paths and mechanisms, so the MD simulation should be more reasonable. Several results stand out from the present calculations. First, one can see clearly that the diffusion migration energies increase from adatom, dimer, to trimer, but drop at tetramer, and then rise again for pentamer, up to the highest cluster checked, octamer. The energy barrier drop from trimer to tetramer is expected because only two atoms are needed to move at each step for the tetramer diffusion, namely, dislocation mechanism.4 Three atoms are needed to move together in the concerted motion of a trimer diffusion,31 because atom-by-atom movement is unlikely and the energy to break a bond is large (also see Figure 2). Second, the migration energy for a Cu8 is
larger than that for a Cu7 cluster. The result is inconsistent with field ion microscopy observations of 2-13 atom 2D fcc metal clusters, showing that close-packed heptamer have the highest diffusion migration energies.6 Our results show that octamer has a very low diffusivity compared with the other clusters, which is consistent with Chirita’s conclusion.32,33 The detailed MD trajectory analysis shows that the periphery motion of a single atom leads to Cu8 cluster shape change, but the motion cannot contribute to the diffusion of Cu8 cluster. The mainly diffusion mechanism for Cu8 is similar to heptamer, namely, the concerted rotation and translation,3 which have a high migration energy. Third, the prefactor for Cu7 is about 30 times larger than that for single atom diffusion, whereas the experiments in Pt and Ir show that the prefactor for compacted surface cluster diffusion is almost 3-4 orders of magnitude larger than the usual prefactors.34-36 The reason for the present result is as follow. Ku¨rpick et al. and Marinica et al.12,37 have calculated the prefactors of a number of hopping mechanisms for self-diffusion of Ir and Cu clusters and found that, for all hopping mechanisms, the prefactors are of the same order of magnitude as in a single atom diffusion. There is no hint of the large prefactor values found in experiments. So Ku¨rpick et al. suggested that the high prefactors found in experiments came about because a large number of nonequivalent processes, for instance the cluster deformed, can participate in diffusion of compact clusters.37 Silva’s MD simulation for Pd7 on Pd(111) indicated that the distortion events do occur in a short time interval. During this time interval, more frequent diffusion events occur, with subsequent recomposition of the original shape of the cluster.38 On the basis of the above analyses, we suggest that the high prefactors found in experiments for the diffusion of compacted Pt7 and Ir7 clusters could arise from different surface character. The fcc metal surfaces can be distinguished in terms of hard, middle, and soft.39,40 The softer the surface is, and the closer the adatom is to the surface, then the stronger the interaction between the adatoms and the surface atoms is. On the hard surfaces, the dominant diffusion mechanism is hopping, while on the soft surfaces, the exchange mechanism is preferred. Cu surface is a hard surface, and the Pt and Ir surfaces are soft surfaces.39 The present simulations indicate that the distance between the Pt adatoms and the surface (0.48a) is about 10% lower than Cu adatoms (0.53a). So the nonequivalent processes on Pt and Ir surfaces, which can lead to a high prefactor, are easier to occur than on Cu surface. Our calculated diffusion prefactor for Cu7 is only 30 times larger than that for a single adatom. Using the local vibrational density of states, Marinica et al.12 calculated the diffusion prefactor for Cu7 cluster on Cu(111) surface is also 15 times larger than that for a single adatom. To further explain this, we also calculate Em and D0 for Pt1 and Pt7 clusters diffusion on Pt(111) surface (Table 1) using the same molecular dynamics simulations. The results are in good agreement with the experimental data. The height prefactors for heptamer is also obtained and in the present MD simulation, we also observe the distortions of the Pt7 atomic structure. In conclusion, because of the different surface character, the prefactor for Cu7 cluster diffusion is only the 1 order of magnitude larger than the single adatom. The mobility of cluster can be measured by characteristic temperature (TD), defined here as the temperature, derived from the Arrhenius plots, at which a cluster moves one nearestneighbor ((x6/6)a) per second.36 The calculated characteristic temperature (TD) as a function of cluster size is shown in Figure
Small Cu Clusters on Cu(111) Surface
Figure 5. Diffusion barrier Em as a function of cluster size for Cun clusters on (111). The fill circles and open circles correspond to the present and static calculated31 values.
Figure 6. Diffusion prefactor D0 as a function of cluster size for Cun clusters on (111).
Figure 7. Characteristic temperature TD for diffusion as a function of cluster size.
7. The trend as the cluster size increases is quite similar with diffusion migration energy. A sharp drop in the diffusion temperature for Cu4 compared with the Cu3 or Cu5 is observed. After the Cu4, the temperature for cluster motion rises monotonically up to the largest cluster investigated, Cu8. Compared to the experiment for Re clusters diffusion on Re(0001),41 our results should be reasonable. As shown in Figure 7, the TD is very low for Cu3, only 93 K, so the Cu3 do not act as nuclear
J. Phys. Chem. C, Vol. 112, No. 6, 2008 2077
Figure 8. Temperature dependence of the order parameter for the Cu7 cluster on the Cu(111) surface. In addition, the temperature dependence of R and h is also plotted.
of three-dimensional growth of thin films on the Cu(111) surface, whereas the heptamer starts to move at 167 K. Using the diffusion prefactors and migration energies, we can estimate the efficiency of the heptamer as a critical nucleus in threedimensional growth of thin films.25 The distance that the heptamer diffuses in time ∆t is x4D∆t, where D is the diffusion coefficient of the cluster. The typical deposition rate is taken to be 17 nm/s, the time of growing one monolayer under typical deposition condition is 0.012 s, which is the time of the heptamer diffusing 60 nm (the dimension for most facets25) at 200 K. Because of the high mobility, the Cu heptamer is ineffective as the critical nucleus at 200 K and above. In order to get deeper insight into the film growth behavior, the thermal stability of the Cu7 cluster on the surface is also investigated. The order parameter η is calculated using the following relation:42,43
η ) R/h
(4)
where R is the mean cluster radius and h is the height of the center of mass from the substrate. The mean radius R was obtained by averaging the distance of atoms from a vertical axis passing through the center of mass. The melting transition of cluster would be signaled by an abrupt increase in η. The temperature dependence of the order parameter for Cu7 cluster is shown in Figure 8. It is clear that the dissociating temperature of Cu7 cluster is about 1150 K, which is about 15% lower than the bulk melting temperature in MAEAM. The result is in good agreement with Gallego’s calculation.42 In addition, in order to illuminate the dissociating process of Cu7 cluster, the temperature dependence of mean cluster radius (R) and height of the center of mass from the substrate (h) is also plotted in Figure 8. As shown, the dissociating for Cu7 cluster is mainly caused by the increasing of R, it indicates that the cluster dissociates into the smaller clusters and loses their hexagonal configuration at the temperature. 4. Conclusion In the present work, the diffusion coefficients for Cu clusters are calculated by the mean square displacement of the masscenter of Cu cluster. The binding energies indicate that heptamer is the most stable, in accordance with the atomic structure of the underlying substrate. The diffusion migration energies increase from adatom, dimer, to trimer, but drop at tetramer. The calculated diffusion prefactor for compact Cu7 cluster is
2078 J. Phys. Chem. C, Vol. 112, No. 6, 2008 only 1 order of magnitude larger than that for single atoms diffusion. The result is different with the experimental results for other fcc metals (Pt and Ir). The reason could be the different surface character for these metals. The calculated characteristic temperature indicates that the heptamer start to move at as low as 167 K. Using the diffusion factor and activation energy, we find that the Cu heptamer is ineffective as critical nucleus in film growth at 200 K and above. The thermal stability for Cu7 cluster is investigated by calculating the order parameter. The computed dissociating temperature of the close-packed cluster is 1150 K, which is about 15% lower than the bulk melting temperature in the present EAM potential. Acknowledgment. We thank Shifang Xiao for useful discussions. This work is financially supported by NSFC (No. 50671035) and High Performance Computing Center of Hunan University. References and Notes (1) Shi, Z. P.; Swan, A. K.; Wendelken, J. F. Phys. ReV. Lett. 1996, 76, 4927. (2) Chirita, V.; Mu¨nger, E. P.; Greene, J. E.; Sundgren, J. E. Surf. Sci. 1999, 436, L641. (3) Hamilton, J. C.; SφÄ rensen, M. R.; Voter, A. F. Phys. ReV. B 2000, 61, R5125. (4) Hamilton, J. C.; Daw, M. S.; Foiles, S. M. Phys. ReV. Lett. 1995, 74, 2760. (5) Voter, A. F. Phys. ReV. B 1986, 34, 6819. (6) Wang, S. C.; Ehrlich, G. Surf. Sci. 1990, 239, 301. (7) Ehrlich, G. Surf. Sci. 1991, 246, 1. (8) Montalenti, F.; Ferrando, R. Phys. ReV. Lett. 1999, 82, 1498. (9) Trushin, O. S.; Salo, P.; Ala-Nissila, T. Phys. ReV. B 2000, 62, 1611. (10) Salo, P.; Hirvonen, J.; Koponen, I. T.; Trushin, O. S.; Heinonen, J.; Ala-Nissila, T. Phys. ReV. B 2001, 64, 161405(R). (11) Yildirim, H.; Kara, A.; Durukanoglu, S.; Rahman, T. S. Surf. Sci. 2006, 600, 484. (12) Marinica, M. C.; Barreteau, C.; Desjonque´res, M. C.; Spanjaard, D. Phys. ReV. B 2004, 70, 75415.
Yang et al. (13) Ala-Nissila, T.; Ferrando, R.; Ying, S. C. AdV. Phys. 2002, 51, 949. (14) Wang, J.; Huang, H.; Cale, T. S. Model. Mater. Sci. Eng. 2004, 12, 1209. (15) Goyhenex. C. Surf. Sci. 2006, 600, 15. (16) Hu, W. New Frontiers of Process Science and Engineering in AdVanced Materials. 2004, 7. (17) Yang, J.; Hu, W.; Deng, H.; Zhao, D. Surf. Sci. 2004, 572, 439. (18) Zhang, Z.; Hu, W.; Xiao, S. J. Chem. Phys. 2005, 122, 214501. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (20) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (21) Ferrando, R.; Tre´glia, G. Surf. Sci. 1995, 331-333, 920. (22) Kellogg, G. L. Prog. Surf. Sci. 1996, 53, 217. (23) Wang, S. C.; Ehrlich, G. Phys. ReV. Lett. 1989, 62, 2297. (24) Liu, S.; Zhang, Z.; Nørskov, J.; Metiu, H. Surf. Sci. 1994, 321, 161. (25) Liu, W. C.; Woo, C. H.; Huang, H. J. Comput. Aid. Mol. Des. 1999, 6, 311. (26) Breemanm, M.; Barkema, G. T.; Boerma, D. O. Surf. Sci. 1995, 323, 71. (27) Kallinteris, G. C.; Evangelakis’, G. A.; Papanicolaou, N. I. Surf. Sci. 1996, 369, 185. (28) Shiang, K. D. Phys. Lett. A 1993, 180, 444. (29) Liu, C. L.; Cohen, J. M.; Adams, J. B.; Voter, A. F. Surf. Sci. 1991, 253, 334. (30) Wynblatt, P.; Gjostein, N. A. Surf. Sci. 1968, 12, 109. (31) Chang, C. M.; Wei, C. M.; Chen, S. P. Phys. ReV. Lett. 2000, 85, 1044. (32) Chirita, V.; Mu¨nger, E. P.; Sundgren, J.-E.; Greene, J. E. Appl. Phys. Lett. 1998, 72, 127. (33) Mu¨nger, E. P.; Chirita, V.; Sundgren, J.-E.; Greene, J. E. Thin Solid Films 1998, 318, 57. (34) Wang, S. C.; Ehrlich, G. Phys. ReV. Lett. 1997, 79, 4234. (35) Wang, S. C.; Ku¨rpick, U.; Ehrlich, G. Phys. ReV. Lett. 1998, 81, 4923. (36) Kyuno, K.; Ehrlich, G. Surf. Sci. 1999, 437, 29. (37) Ku¨rpick, U.; Fricke, B.; Ehrlich, G. Surf. Sci. 2000, 470, L45. (38) Silva, E. Z. da; Antonelli, A. Surf. Sci. 2000, 452, 239. (39) Zhuang, J.; Liu, L. Phys. ReV. B 1999, 59, 13278. (40) Zhuang, J.; Liu, L.; Ning, X.; Li, Y. Surf. Sci. 2000, 465, 243. (41) Goldstein, J. T.; Ehrlich, G. Surf. Sci. 1999, 443, 105. (42) Longo, C. R.; Rey, C.; Gallego, L. J. Surf. Sci. 2000, 459, L441. (43) Antoneli, A.; Khanna, S. N.; Jena, P. Phys. ReV. B 1993, 48, 8263.
