Temperature dependent structure of low index copper surfaces studied by molecular dynamics simulation

Temperature dependent structure of low index copper surfaces studied by molecular dynamics simulation

F.J. Resende^I; V.E. Carvalho^I; B.V. Costa^I; C.M.C. de Castilho^II

^IDepartamento de Física ICEX, Universidade Federal de Minas Gerais, Caixa Postal 702, 30123-970, Belo Horizonte, MG, Brazil

^IIInstituto de Física, Universidade Federal da Bahia, Rua Barão de Geremoabo sn, Campus Universitário da Federação, 40210-340, Salvador, BA, Brazil

ABSTRACT

The thermal behavior of the (010), (110) and (111) copper surfaces is studied by molecular dynamics simulation. We have used a many-body potential based on the tight-binding model in order to describe the Cu-Cu interaction. The calculations we have performed correspond to simulations in the temperature range between 600 and 1800 K. The observed order in the stability follows the same order as in the packing density, i. e., (110), (010) and (111). The (110) disorder results from anharmonic effects and by vacancy-adatom formation. On the other end, the (111) surface is very stable, and remains so up to temperatures of the order of the bulk melting point. The melting proceeds by a layer-by-layer mechanism.

1 Introduction

The thermal behavior of solid surfaces has raised a lot of attention in recent years as a result, in a large extent, of refinements in either experimental or theoretical techniques. Low-energy electron diffraction (LEED) [1, 2], medium-energy ion scattering (MEIS) [3, 4, 5, 6] and X-ray diffraction [7] are some of the experimental techniques well suited to the study of the structure of metallic surfaces and their behaviour as the temperature is varied [8]. These experiments indicate that quantities such as surface thermal expansion, mean-square displacements and structural parameters often exhibit a behaviour which are anomalous when compared to the bulk one. These behaviours have been attributed to strong anharmonic vibrations of the atoms at the surface. This anharmonicity brings disorder to the surface, mainly at high temperatures [8, 9, 10]. These experimental techniques have been applied to study the copper (010), (110) and (111) surfaces. The results for the (110) face have shown an enhancement of the disorder for temperatures above T = 550K and a roughening transition at T_R = 870K [98] and T_R = 1070K[10]. For the (010) face, LEED experiments show an enhancement of the atoms vibrational amplitude at the surface [11]. In the case of the (111) face, experimental studies by MEIS show an anomalous behaviour of the thermal expansion coefficient for temperatures above 1000K and an anisotropy of the vibrations at the surface, with the atoms vibrational amplitudes in the direction normal to the surface higher than those in the plane of the surface [12]. Numerical simulational methods, like Monte Carlo (Lattice-gas model and solid-on-solid model) and Molecular Dynamics (MD), have been used, with a relative success, to study the structure and dynamics of the surfaces. The main ingredient of the MD approach is the interatomic interaction potential. Realistic potentials, for metallic systems, have been available only during the past few years. Up to now, MD simulations, implemented with these potential models, have been made for Cu(110), using the effective-medium theory (EMT) [13, 14] and the Embedded-atom method (EAM) [15]. For Cu(010), the Effective-medium theory (EMT) [16] was applied and for Cu(111) a tight-binding (TB) method in the second moment approximation was used [23]. A premelting phenomenon, driven by the vacancy-adatom mechanism, has been observed in simulations of the Cu(110) surface [13, 15]. The (111) face of copper, on the other hand, seems to be quite stable, and maintains the crystalline state until high temperatures. In these studies, strong effects of the anharmonic vibrations of the atoms at the surface have been also observed for Cu(110) [15]. The anharmonicity manifests itself as a larger thermal-expansion coefficient of the surface area and also as an anomalous increase in the vibrational amplitudes of the surface atoms. In these previous simulations, the anharmonicity effects are less pronounced at (010) and (111) faces [16].

In this work we have studied the structure of the (110), (010) and (111) faces of copper, in the temperature range from 600 up to 1800K, by classical Molecular Dynamics simulation. With the purpose of directly comparing the behaviour of the different crystal faces, a task which is generally difficult if the comparison involves distinct simulation approaches, we have employed comparable sample sizes and the same many-body potential, based on the tight-binding model, to describe the Cu-Cu interaction.

This work is organized as follows. In section 2 we describe the simulation scheme, presenting the results in section 3. The conclusions and some final remarks are presented in section 4.

2 Simulation Method

Our simulation was carried out using molecular dynamics (MD). We have used a semi-empirical potential [17, 18, 19], based on the tight-binding method for the second moment approximation, in order to describe the attractive interatomic interaction, while a Born-Mayer-type interatomic repulsion was assumed. This potential model has given good results in previous MD calculations of bulk [18, 19] and surface [20, 21, 22, 23, 24, 25, 26] properties of transition metals. For a system with just one kind of atomic species, the total potential energy is given by

In this equation, indexes i and j stands for particles at position _i and _j, respectively, N is the total number of particles, r_i,j = |_i – _j| is the distance between atoms i and j and r₀ is the nearest-neighbour equilibrium distance in the bulk metal. The parameters A, x, p and q are adjusted in order to reproduce the cohesive energy, the equilibrium condition (zero pressure at T= 0 K) and the elastic constants of the metal under consideration. For Cu-Cu we have used A = 0.0855eV, x = 1.224eV, p = 10.96, q = 2.278 and R₀ = 2.5562, according to ref.[42]. A cutoff is introduced in the potential, so that it is possible to account for interactions up to the fifth nearest neighbours, as a route to accelerate the simulation process. Particles in the simulation move according Newton's law and this generates a set of 3N coupled equations of motion, which are simultaneously solved by increasing forward in time the physical state of the system and that is done in small time-steps of size dt = 10^–14s. To do that we have used Beeman's method of integration [26, 27, 28, 29, 30, 31], which is a fourth order predictor-corrector scheme. To improve the integration procedure we have also used a Verlet table [26, 27, 28, 29, 30, 31].

The copper substrate is composed by distributing N particles over a frozen layer in order to simulate the presence of the semi-rigid solid. Particles are then distributed around their equilibrium position, filling 4 layers in the y direction, which is the one perpendicular to the substrate. We also performed some simulations in a system with seven layers as a control procedure.

Periodic boundary conditions are used for the directions parallel to the substrate while an open condition is set for the y direction. Once the temperature (T) of the system is fixed, a process of velocity renormalization is used to equilibrate the system by letting it to evolve for 2 × 10³dt time steps. The system is assumed to be in equilibrium provided the deviations in the mean kinetic energy becomes less than 1%. With the aim of estimating the melting temperature of a slab, we have chosen, as a reference, the (010) surface. In Fig. 1 is presented the total energy (per molecule) as a function of temperature for systems with 4 (N » 400 atoms) and 7 layers (N » 700 atoms) in the y direction and surface dimensions of about (25Å × 25Å). It is easy to see that the melting temperature (T_m) does not change very much as the number of layers is varied; in both cases T_m > 1400K. Thus, we have chosen to work with systems of 4 layers and N » 400 atoms, restricting our attention to temperatures in the range 600K to 1100K, well bellow the estimated melting point. In order to have reasonable statistics, our results correspond to averages over simulations of 200ps and 6 independent initial configurations.

3 Results

To evaluate the degree of disorder in the structure and how it is affected as a result of temperature variation, we have calculated the structural parameter l, which is defined by [32]

where k^m correspond to the reciprocal-lattice vectors, () are the coordinates of the atom i in the plane parallel to the surface and N is the number of atoms in the plane. The value of l is usually considered as an order parameter, being 1.0 in the case of a perfectly ordered layer at low temperature, approaching zero for a completely disordered layer.

Figure 2 shows, as a function of temperature, the values of l for the first (free surface) layer and also for the second one, in case of the Cu(010) face. For low temperatures l » 1 and, as the temperature arises, l tends to zero. As can be seen from this figure the values of l for the first layer are lower than those for the second layer and also go to zero faster. This behaviour is associated with a higher disorder in the surface layer, as compared with the second one, and it is a strong evidence of the surface premelting process. In this process, the melting of the material starts from the surface, with the formation of a thin liquid-like film at a temperature bellow the thermodynamical bulk melting point T_m. This film propagates inward the solid as the temperature increases until the process is concluded at T_m.

Figure 3 shows for comparison the parameter l as calculated for Cu(010), Cu(110) and Cu(111) surfaces in the case of the first layer, whereas, in figure 4, the same is shown for the second layer. It is possible then to observe that the surface stability follows the packing density order, i. e., the disorder is most pronounced on the (110) face, much weaker on the (100), and the face (111) behaves as a very well ordered surface up to T_m.

3.1 Vacancy-adatom Formation

It has been observed in both theoretical and experimental studies that at elevated temperatures the thermal generation of vacancy-adatom pairs, at the outermost surface layer, might be responsible for some of the disorder observed at the surface [13, 14, 15, 16, 23]. Vacancies at the surface have then an indirect effect on the order parameter, by inducing distortions in the surface region. In our simulations we have studied the vacancy-adatom formation by evaluating the layer occupation over the crystal surface at each time step in the temperature range from 600 to 1100K. We have observed that on the (110) face, the formation of the adatom layer begins at around 700K, whereas for the (010) and (111) faces the onset of the vacancy-adatom formation takes place at higher temperatures, 800K and 900K, respectively. For the three explored faces, the adatom related vacancies are mainly found in the first crystal layer. These observations are in agreement with other simulational studies [13, 14, 15, 16, 23]. Figure 5 shows an Arrhenius plot for the adatom concentration. In this figure it has been included just the high-temperature points, where good statistics have been achieved. From the slope of the curves it is possible to derive the activation energy for adatom formation as 0.43(6)eV, 0.9(1)eV and 1.4(1)eV for the (110), (010) and (111) faces, respectively. In Table I is presented the obtained results from this work as well as those from studies using molecular-dynamics (MD), the tight-binding method (TB) [23, 35], effective medium theory (EMT) [16] and embedded atom method (EAM) [36] and experimental studies of X-ray diffraction [37]. As can be seen from the table our results agree well with those from other works and Cu(111) is the face presenting the largest discrepancies.

3.2 Surface Relaxation

Thermal expansion is a direct manifestation of anharmonicity. For a surface, the broken inversion symmetry at the surface causes an asymmetry in the potential, as seen by surface atoms. The potential is quite generally expected to be softer in the region near the surface, thereby causing atoms to be displaced further away, as compared with the situation deeper in the bulk. This motion, of course, depends strongly on temperature and, in general, it is also affected by the particular nature of the surface under consideration. For layered systems, thermal expansion is directly observed in the dependence with temperature of the interlayer distance d. While the thermal expansion coefficient for the bulk shows only a weak temperature dependence, even at values approaching melting, the situation is quite different for the surface. The reported results have shown that, in general, the thermal-expansion coefficient at the surface (a_surface) is larger than that of the bulk (a_bulk)[38, 39, 16]. We have calculated the surface relaxation using the difference between the mean distance of the first and second layers of the crystal. As it is shown in Fig. 6, the calculated surface relaxation dd₁₂ for the faces (010), (110) and (111) as a function of the temperature, defined as the percent change of the inter-layer spacing with respect to the value of the bulk lattice parameter at each temperature. It should be expected that, in this kind of plot, if the thermal expansion coefficient for the surface is the same as that of the bulk, these data will fall on a horizontal line parallel to that corresponding to the bulk. In fact, the curves seem to be nearly linear for the three faces although showing different slopes. This is an indication of an occurrence of different thermal expansion coefficient for the three surfaces, which can be calculated from a dd₁₂ × T plot. We have carried out this calculation and found a₀₁₀ = 87(6) × 10^–6K^–1 ( about five times a_bulk), a₁₁₀ = 108(9) × 10^–6K^–1 (six times a_bulk) and a₁₁₁ = 32(6) × 10^–6K^–1 (two times a_bulk). Analysing figure 6, the first point that may be mentioned is that the expansion of the three surfaces shows a behaviour similar to that found in other theoretical and experimental works, that is, an occurrence of a contraction at low temperature and an expansion at higher temperatures [38, 39]. Furthermore, for the Cu(010) surface in the temperature range T ~ 500 – 700K, the results from MEIS experiments, reported by Fowler [40], show a contraction (dd₁₂ = –2.3 %) that is very close to our result. On the other hand, results from MD using EAM, show a contraction of 0.8% for T = 600 K[41]. Using EMT, the results of Hakkinen and Manninen [16] show contraction between 0.3% and 0.6% for a temperature below 400K and an expansion between 0.3% and 1% for a temperature between 600K and 900K. In this same study, a linear behaviour of the relaxation curves with a thermal-expansion coefficient a = 38(8) × 10^–6K^–1, was observed i.e., twice a_bulk. In another MD simulation, a contraction of 0.6 % for T = 1000K was found while above T = 1000K a strong expansion was observed [23]. The thermal-expansion coefficient, in the low temperature range, is 37(5) × 10^–6K^–1, i. e., twice a_bulk. For higher temperatures, above T = 1000K, a grows up to five times a_bulk. In our simulations we have calculated a thermal-expansion coefficient of a₀₁₀ = 87(6) × 10^–6K^–1, i. e., five times a_bulk. For the (110) copper face the results obtained by Hakkinen and Manninen show a strong contraction of 4 %, for low temperatures, and for high temperatures this contraction goes to 1 % at T = 900K [16]. The relaxation curve, in these studies, is approximately linear with a = 32(5) × 10^–6, i. e. , twice a_bulk. A similar behaviour is observed in MD simulations using TB [35], with a contraction of 2 % at T = 400K and of 0.4 % at T = 1000K. Experimental results by LEED and X-rays, at low temperatures, show contractions up to 6 %. The MD study of Ditlevsen et al [14] using EMT shows a significant expansion of the surface interlayer distance. For low temperatures the curve is linear, with a thermal-expansion coefficient tree times a_bulk. At high temperatures the surface inter-layer distance rises dramatically. This is attributed to anharmonic effects. The value of the thermal-expansion coefficient for (110) face, obtained in the present work, is a₁₁₀ = 108(9) × 10^–6K^–1, that is six times the bulk value a_bulk.

MEIS experimental studies on the (111) copper surface [12] show a contraction of the first inter-layer distance of 1% in the temperature range of 300K up to 800 K. For temperatures between 800K and 900K, the contraction vanishes and for temperatures above 900K it is observed a strong expansion of 4% at T = 1180K with a thermal-expansion coefficient six times the bulk value. The results of MD using TB show a similar contraction, in the low temperature region, but an expansion of only 0.2 % for T = 1200K [23] and a linear behavior with a thermal-expansion coefficient of 10(5) × 10^–6 , very close to the bulk value. This behaviour is also observed in the results obtained by Hakkinen and Manninen and the contraction at low temperatures being nearly 1 %. Again, the thermal-expansion coefficient is very close to the bulk value. In our study we have observed a linear behavior and the thermal-expansion coefficient obtained is a₁₁₁ = 32(6) × 10^–6K^–1, i. e., approximately twice a_bulk. So, the results for the thermal behaviour of the interlayer distance of the three surfaces investigated in this work show an overall qualitative agreement with the reported results from other techniques.

3.3 Conclusions

In the study of the thermal behavior of the (010), (110) and (111) copper surfaces by molecular dynamics (MD) simulations, we have observed that the order in the stability follows the same order as in the packing density. The (110) face presents a considerable disorder at relatively low temperatures (Fig. 3), probably due to anharmonic effects, what is evidenced by the surface thermal-expansion coefficient higher than the bulk value (Fig. 6). Above T = 700K the onset of the vacancy-adatom formation takes place. This observation, although do not represent results of direct measurements, shows that roughening and premelting transitions, as observed in other works [16], are possible in the Cu(110) surface. The (010) copper face presents disorder only for temperatures above T = 800K and, at this temperature, it is observed the beginning of the vacancy-adatom formation. In this surface, the anharmonic effects are also present, with a surface thermal-expansion coefficient higher than the bulk value. The (111) surface is ordered up to high temperatures and presents no significant vacancy-adatom pairs formation, with a surface thermal-expansion coefficient only two times that bulk value, i. e., the presence of only weak anharmonic effects. We have not observed an anomalous behaviour of the thermal-expansion coefficient on these faces. These results show qualitative agreement with previous theoretical and experimental studies.

Financial support from the Brazilian agencies CNPq, CAPES, FAPEMIG and CIAM-02 49.0101/03-8 (CNPq) are gratefully acknowledged.

Received on 6 January, 2004

• [1] E. A. Soares, V. B. Nascimento, V. E. Carvalho, C. M. C. de Castilho, A. V. de Carvalho, R. Toomes, and D. P. Woodruff, Surf. Sci. 419, 89 (1999).
• [2] E. A. Soares, G. S. Leatherman, R. D. Diehel, and M. A. Van Hove, Surf. Sci. 468, 129 (2000).
• [3] B.W. Busch and T. Gustafsson, Surf. Sci. 407, 7 (1998).
• [4] G. Bracco, M. Canepa, P. Cantini, F. Fossa, L. Mattera, S. Terreni, and D. Truffelli, Surf. Sci. 269/270, 61 (1992).
• [5] G. Bracco, A. Pizzorno, and R. Tatarek, Surf. Sci. 377/379, 94 (1997).
• [6] L. Pedemonte, G. Bracco, R. Tatarek, M. Aschoff, K. Brüning, and W. Heiland, Surf. Sci. 482-485, 1457 (2001).
• [7] G. Helgesen, D. Gibbs, A. P. Badddorf, D. M. Zehner, and S. G. J. Mochrie, Phys. Rev. B48, 15320 (1993).
• [8] D. Gorse, J. Lapuloulade, Surf. Sci. 162, 847 (1985).
• [9] S. G. J. Mocrie, Phys. Rev. Lett. 59, 304 (1987).
• [10] P. Zeppenfeld, et al. Phys. Rev. Lett. 62, 63 (1989).
• [11] G. Armand, et al., Europhys. Letters. 3, 1113 (1987).
• [12] K. H. Chae, H. C. Lu, and T. Gustafsson Phys. Rev. B 54, 14082 (1996).
• [13] H. Hakkinen, J. Merikoski, Phys. Rev. lett. 70, 2451 (1993).
• [14] P. D. Ditlevsen, et al. Phys. Rev. B 44, 13002 (1993).
• [15] M. S. Daw and M. I. Baskes. Phys. Rev. B 29, 6443 (1984); M. Foiles, M. I. Baskes and M. S. Daw, Phys. Rev. B 33, 7983 (1986).
• [16] H. Hakkinen and M. Manninen Phys. Rev. B 46, 1725 (1992).
• [17] D. Tomanek, A. A. Aligia, and C. A. Balseiro, Phys. Rev. B 32, 5051 (1985).
• [18] Carlo Massobrio, V. Pontikis, and G. Martin, Phys. Rev. Lett. 62, 1142 (1989).
• [19] F. Willaime and Carlo Massobrio, Phys. Rev. Lett. 63, 2244 1989).
• [20] K. -D. Shiang et al, Surf. Sci. 301, 136 (1994).
• [21] N. I. Papanicolaou et al, Surf. Sci. 337, L819 (1995).
• [22] G. A. Evangelakis, N. I. Papanicolaou , Surf. Sci. 347, 376 (1996).
• [23] G. C. Kallinteris et al, Surf. Sci. 369, 185 (1996).
• [24] G. A. Evangelakis et al, Surf. Sci. 394, 185 (1997).
• [25] N. Levanov et al, Surf. Sci. 400, 54 (1998).
• [26] F. J. Resende and B.V. Costa, Surf. Sci. 481, 54 (2001).
• [27] M. P. Allen and D. J. Tildesley,Computer Simulation of Liquids, Oxford Science Publications, New York (1992).
• [28] D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, New York (1997).
• [29] P. Z. Coura and B.V. Costa, Int. J. Mod. Phys. C 9, 857 (1998).
• [30] P. Z. Coura, O. N. Mesquita, and B.V.Costa, Phys. Rev. B 59, 3408 (1999).
• [31] F. J. Resende and B.V.Costa, Phys. Rev. B 61, 12697 (2000).
• [32] C. P. Toh, C. K. Ong, and F. Ercolessi, Phys. Rev. B 50, 17507 (1994).
• [33] D. Passerone, F. Ercolessi, F. Celestini, and E. Tosatti Surf. Rev. Lett. 6, 5, 663 (1999).
• [34] Ngoc-Thanh Vu and D. B. Jack, Surf. Rev. Lett. 6,5, 683(1999).
• [35] G. A. Evangelakis, D. G. Papageorgiou, G. C. Kallinteris, Ch. E. Lekka, and N. I. Papanicolaou, Vacuum 50, 165 (1998).
• [36] M. Karimi et al, Phys. Rev. B 52, 5364 (1995).
• [37] H. Durr, R. Shneider, Phys. Rev. B 43, 12187 (1991).
• [38] B. W. Bush and T. Gustafsson, Surf. Science. 6, 5, 683 (1999).
• [39] Q. T. Jiang, P. Fenter, and T. Gustafsson, Phys. Rev. B 44, 5773 (1991).
• [40] D. E. Fowler, J. V. Barth, and T. Gustafsson Phys. Rev. B 52, 2117 (1995).
• [41] Liqiu Yang, Talat S. Rahman, and Murray S. Daw, Phys. Rev. B 44, 13725 (1991).
• [42] F. Cleri and V. Rosato, Phys. Rev. B 48, 22 (1993).

Publication Dates

• Publication in this collection
  01 Sept 2004
• Date of issue
  June 2004

History

• Accepted
  06 Jan 2004
• Received
  06 Jan 2004