J. Phys. Chem. C 2010, 114, 8717–8720
8717
Melting of Cu Nanowires: A Study Using Molecular Dynamics Simulation W. X. Zhang*,† and C. He‡ School of Materials Science and Engineering, Chang’an UniVersity, Xi’an 710064, Shannxi, China, and State Key Laboratory for Mechanical BehaVior of Materials, Xi’an Jiaotong UniVersity, Xi’an 710049, China ReceiVed: October 27, 2009; ReVised Manuscript ReceiVed: April 12, 2010
The melting behaviors and transport properties of Cu nanowires (NWs) are studied for the future application as interconnects in microelectronics by using first-principles molecular dynamics (FP-MD) and classical molecular dynamics (MD) methods. The results of the potential energy, Cp(T), the quantum conduction, and the density of states are used to monitor the phase transition. It is found that the melting temperature using the MD method is 10% higher than that using the FP-MD method, and the corresponding FP-MD results of Cu NWs are predicted with the diameter of 2-6 nm. Meanwhile, the size effects of the nanowires on the melting temperature are validated by the theoretical model. 1. Introduction The physical properties of metallic nanowires (NW) have been investigated in both industrial and academic fields because of their dramatically different behavior from bulk materials, such as their electronic, magnetic, optic, catalytic, and thermodynamic properties.1-3 Among these properties, much attention has been paid to the studies of the thermal stability and melting behavior of metallic NWs.4-7 Many researchers have reported the anomalous melting behavior of NWs due to their large surface-to-volume ratio by means of thermodynamic models and molecular dynamics (MD) simulations.6-10 For instance, Wang et al. have reported that the melting temperatures of ultrathin Ti NWs show dependence on wire sizes and structures.8 Ouyang et al. have discussed a thermodynamic model to address the surface energy of Zn, Ni, and Pd NWs and deduced the analytic expression of the sizedependent surface energy of NWs.11 Liu et al.12 have found that Pt NWs with a 3 nm diameter are unstable at 673 K, whereas the melting point of the bulk is 2045 K. As the surface fraction of atoms increases with decreasing size, the melting temperature of a free-standing clear nanowire is expected to decrease with decreasing size; it is because the atomic coordination number imperfection lowers the atomic cohesive energy that dominates the thermal stability, such as melting13 and phase transition.14,15 Therefore, the melting proceeds from the surface.16,17 Cu as interconnects in microelectronics is the most useful metal while the electronic transport of nanosized interconnects is one of the important characteristics for future microelectronic applications.9 Most studies had been dedicated to structures, quantum behaviors9,18 under electric fields, and optical properties2 of Cu NWs. However, there are a few reports involved in the theories and experiments about the thermodynamic properties of Cu NWs in the present literature.2,9,18 The thermodynamic properties of metallic NWs differ from that of the bulk counterpart, which play significant roles in the melting behavior.2,9,18 It is of fundamental interest to understand how the properties’ dependence on size and, in particular, on the bulklike limit. Furthermore, most studies are very useful for * To whom correspondence should be addressed. E-mail: wxzhang@ chd.edu.cn. Fax: 86-29-82337340. † Chang’an University. ‡ Xi’an Jiaotong University.
understanding the general melting features of the ultrathin NWs. The diameter of the metallic NWs being studied in the literature ranges from 0.1 to 2 nm. Most efforts have been focused on sizes below 1 nm for both Ti and other metals.8,18 To facilitate predicting simulation data, we investigate the melting behavior of the Cu NWs with a diameter of 2-6 nm. In general, it is difficult to observe the melting process of NWs directly by experiment. Computer simulation methods, in particular, firstprinciples molecular dynamics (FP-MD), are ideal tools for investigating the structural, mechanical, and thermodynamic properties of such nanostructures and improving our understanding in various applications. MD simulations have been employed to study the structures and properties of free-standing metal Ti6,8 and Zr7 NWs by Wang et al. In this paper, melting characteristics of Cu NWs are described by ab initio density functional theory (DFT) with the generalized gradient approximation (GGA). The simulation study of Cu NWs provides an opportunity for further understanding its unique role in experimental phenomena. The atomic structures of Cu NWs with diameters from 2 to 6 nm have been optimized, and their size-dependent melting properties have been discussed. The quantum conduction G(T) function and the density of states (DOS) are performed to determine changes of atomic and electronic structures of Cu NWs. Meanwhile, an analytic thermodynamic model with available experimental data is applied to verify the melting that is dependent on Cu NW size on the basis of the thermodynamics mechanics. The findings in this work may provide new insight into the fundamental understanding of the thermal stability and structural evolution of the materials with the size from bulk to diameter. 2. Simulation Details For bulk Cu and Cu NWs (D ) 2 nm), simulations were performed based on density functional theory (DFT) with the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) exchange-correlation functionals.19 Similar functionals have been successfully used to study the structural and electronic properties of water.20-22 In our work, the electron-ion interaction is described by ultrasoft pseudopotentials23 with a cutoff energy of 380 eV, where sufficient numbers of wave functions are included so as to get precise information about the electronic structure of the crystals.
10.1021/jp910246t 2010 American Chemical Society Published on Web 04/26/2010
8718
J. Phys. Chem. C, Vol. 114, No. 19, 2010
Zhang and He
Figure 1. Potential energy of bulk Cu for heating and cooling cycles with different simulation methods.
In the first-principles molecular dynamics methods (FP-MD), the Cp(T) function plays an important role in understanding thermodynamic functions of Cu. It is determined by24
Cp )
∂H 〈∆H2〉 ) ∂T kT2
(1)
where k is the Boltzmann’s constant and 〈∆H2〉 is the meansquared fluctuation in enthalpy. Classical MD simulations are performed to simulate Cu in the NPT statistical ensemble with the COMPASS (condensedphase optimized molecular potentials for atomistic simulation studies) force field based on the earlier class II CFF9x and PCFF force fields, which is the first ab initio base force field and has been parametrized using extensive data for molecules in the condensed phase. P and T are constants and N is the atom number. T is imposed by the Nose´-Hoover algorithm.25 The integration step is chosen as 1 fs using the Verlet leapfrog algorithm.26 Configurations saved every 1 ps with a 2 K step per state at the range of 300 K < T < 2100 K are kept with 0.1 ns. The constant-temperature MD simulations start from a low temperature (300 K). The initial configuration for any given T is taken to be the final one from the previous T. To confirm the equilibrium time t being adequate, the energy fluctuation in the energy evolution E(t) is analyzed. When t ) 100 ps, the error range is 1%, which is allowed in this simulation. This is also confirmed from our earlier simulation works and other literature.10,27,28 Therefore, in the later simulation, t ) 100 ps is taken. After the MD simulation, fluctuations in the NPT ensemble are analyzed. Thus, Cp at a given T can be calculated by
Cp(T) )
1 〈δ(κ + p + PV)2〉 RT2
(2)
where κ and p denote the instantaneous values of the kinetic and potential energy and P, V, and T show the familiar thermodynamic state variables. In addition, the notation δX means X - , where denotes the equilibrium ensemble average value of quantity X. In the simulation, δ(κ + p + PV)2 is directly given by analyzing results. 3. Results and Discussion Figure 1 shows the potential energy of bulk Cu versus temperature using FP-MD, COMPASS, and the embedded-atom method (EAM)17 during heating and cooling processes. The melting and cooling transition occurs at the rapidly increasing
Figure 2. Cp vs T plots of bulk Cu for heating cycles with different simulation methods. The experimental melting temperature of Tmelt point ) 1356 K is for bulk Cu.
and decreasing regions in energy, respectively. Considering the potential energy of heating versus T plot for bulk Cu in the heating process, this exhibits three temperature regions. At the beginning, the potential energy has no big changes, and the bulk Cu exhibits a solid-state behavior. In a solid crystal, the atoms are located in well-defined lattice equilibrium positions. Upon heating, energy undergoes a strong increase, and bulk Cu undergoes a transition and exhibits a behavior intermediate between the solid and liquid states. Upon increasing the temperature, energy is relatively high compared with the case of the solid state. As for the potential energy of the cooling process, the plot also exhibits three regions. The trend corresponds with other simulations.17 The melting and solidifying transition occurs at the rapidly increasing and decreasing regions in energy, respectively. Here, we obtain the Tmelt and Tsolid from the derivative of potential energy, which is shown in Figure 2, leading to Tmelt(FP-MD) ) 1460 K, Tsolid(FP-MD) ) 1240 K, Tmelt(COMPASS) ) 1620 K, and Tsolid(COMPASS) ) 1060 K. The experimental melting temperature of Tmelt point ) 1356 K for pure Cu is between them, and the difference is less than 100 K for FP-MD and about 200 K for the COMPASS method, which are more exact than EAM: Tmelt(EAM) ) 1800 K, and Tsolid(EAM) ) 824 K.17 According to the results above, FP-MD based on density functional theory is used to calculate the melting point of bulk Cu, and the difference between COMPASS and FP-MD is about 10%, which is very useful for us to predict the corresponding FP-MD melting point of Cu NWs. Meanwhile, we conclude that the FP-MD method is the most suitable method of simulating the melting point of Cu, then is the COMPASS method, and both of them are more suitable for simulating the melting of Cu than the EAM method. Starting from the optimized structures at low temperature, the melting behavior of the bulk Cu has been simulated by using different methods. It also can be proved by the plot of Cp versus T in Figure 2, which is a characteristic for the first-order solidto-liquid phase transition. The rapid jump with increasing temperature demonstrates the occurrence of melting. We define Tm as the max value of heating Cp. The Cp versus T plots indicate that the melting points are 1460, 1620, and 1771 K, respectively, according to FP-MD, COMPASS, and EAM methods.17 Comparing the results from these methods with the experimental results, we can easily obtain the same sequence of suitability of the three methods. The periodic boundary condition calculation provides no surface for the bulk system, leading to a homogeneous melting transition with a high superheating temperature. However, the existence of the free surface for NWs leads to a heterogeneous melting transition temperature lower than the experimental value.
Study of the Melting of Cu NWs Using MD Simulation
J. Phys. Chem. C, Vol. 114, No. 19, 2010 8719
Figure 3. Potential energy vs T plots when the diameters of Cu NWs are 2, 3, 4, 5, and 6 nm.
TABLE 1: G(T) Function of Cu Nanowire with D ) 2 nm Obtained from FP-MD Calculations Cu NWs 300 K 400 K 420 K 440 K 460 K 480 K 500 K 800 K G(T)
24
24
24
20
18
16
16
16
To explore the melting temperature of Cu NWs with different diameters, we plot the potential energy versus heating temperature in Figure 3 by the COMPASS method. Although we admit the best method to simulate the Tm of Cu NWs is the FP-MD method, it cannot be achieved when D > 2 nm for the limitation of computational power at present. We just simulate the potential energy of a 2 nm Cu NW by the FP-MD method and have compromised to choose the COMPASS method to simulate the potential energy above 2 nm. For comparison, the potential energy versus T plot of bulk Cu is also put into Figure 3. Tm(FP-MD) (2 nm) ) 440 K, and Tm(COMPASS) (2 nm) ) 480 K. It is found the melting temperature using the MD method is 10% higher than that using the FP-MD method, which corresponds to the result of bulk Cu. Therefore, the difference between COMPASS and FP-MD is about 10% is suitable for both bulk and Cu NWs. In the similar way as above, Tm(COMPASS)(D) ) 480, 900, 1280, and 1360 K when D ) 2, 3, 4, 5, and 6 nm are also obtained according to Figure 3. With the decreasing of diameters, the change of the energy exhibits the opposite trend. That the energy per atom of the NWs is higher than that of the bulk at the same temperature indicates the existence of the surface energy. Besides the melting behavior, the transport properties are very important for Cu NWs as interconnects. It is found that the changes of atomic structure and related structure decide the value of the quantum conduction G. According to the Landauer formula, the number of bands crossing the Fermi level Ef attributes to the number of conductional channels or the value of G.29 Therefore, FP-MD calculations are performed to determine the G function of Cu nanowire with D ) 2 nm, which is limited by the present simulation power. The determined G(T) values of Cu nanowire are shown in Table 1. The details can be found in our previous paper.9 Meanwhile, the density of states (DOS) is performed to determine changes of atomic and electronic structures of Cu nanowire with D ) 2 nm, which are present in Figure 4. The largest peak of DOS below Ef is located at -1.67 eV under 300 K and -1.54 eV under 800 K. As T increases, the largest peak obviously shifts right about 0.13 eV, which implies the energy increase. Meanwhile, the value of DOS at Ef under 800 K is 0.67 times that under 300 K, which is in great agreement with G (800 K)/G (300 K) ) 0.66 in Table 1. According to Table 1, we find that the melting process is between 420 and 480 K. Once T is above the melting point,
Figure 4. DOS of Cu nanowire with D ) 2 nm: T ) 300 K (solid line) and T ) 800 K (dashed line). Ef ) 0 (vertical dotted line) is taken.
the value of G decreases obviously. When T is higher than 480 K, the value of G is constant. Therefore, the value of G will not change with temperature until the liquefaction is completed. Hence, with Cu as interconnects in microelectronics, the ambient temperature should not be too high. Otherwise, once the melting process between the solid and liquid states appears, the conductivity of interconnects will drop dramatically and will lead to the failure. In particular, as the size of interconnects further reduces, the melting temperature would decrease and then the impact of ambient temperature on the conductivity would be more apparent. The size-dependent melting temperature of metallic nanowires Tm(D) can be described by Jiang’s model30-34
Tm(D)/Tm(∞) ) exp[2∆Svib(∞)/(3R)/(1 - D/4h)]
(3) where Tm(∞) is the melting temperature for the corresponding bulk crystal, D denotes the corresponding diameter, R is the ideal gas constant, and h is the atomic diameter. ∆Svib shows the vibrational part of the overall melting entropy ∆Sm.35-37 For metallic crystals, ∆Svib(∞) ) ∆Sm(∞) + R(xA ln xA + xB ln xB). For the melting transition, xA ) 1/(1 + ∆Vm/Vs) and xB ) 1 xA, where ∆Vm ) Vl - Vs, with V being the molar volume and the subscripts s and l denoting the crystal and the liquid, respectively. According to eq 3, the melting temperature of Cu metallic NWs can be quantitatively determined without the free parameter, which is shown in Figure 5. The related parameters are listed in the caption of Figure 5. For convenient comparison, our COMPASS simulation results (solid symbols), available simulation results (open symbols),3 and the liquid drop model with different β values (dashed lines) are also listed, where β is a material parameter.38-41 Meanwhile, according to the difference of 10% between COMPASS and FP-MD methods, the corresponding FP-MD results are predicted, which is also shown in Figure 5. It is shown that the melting temperature decreases with the decreasing diameter of Cu NWs, and the variation tendency of the results and these of the liquid drop model is consistent with the simulation results in the full size range. Note that, although the simplest function of Tm(D) is linearly proportional to 1/D that denotes the surface/volume ratio,19,42-44 the Tm(D) value in terms of eq 3 is apart from this linear relationship as D decreases. This fact implies that the size dependence of melting does not solely correspond to the surface/volume ratio, but to the energetic change of interior atoms of the NWs, which becomes evident when D is small. To explore the size dependence of the melting temperature for
8720
J. Phys. Chem. C, Vol. 114, No. 19, 2010
Figure 5. Tm(D) function of Cu NWs (solid line) in terms of eq 3. The necessary parameters are h ) 0.2556 nm,31 Tm(∞) ) 1356 K,24 Hm(∞) ) 13.05 KJ mol-1,31 ∆Vm/Vs ) 4.9,31 and ∆Svib(∞) ) 8.06 J mol-1 K-1. The dashed lines show other theoretical predictions with β ) 1.0240 and 0.9.41 The symbol 2 shows our COMPASS simulation results, and ∆3 shows other simulation results. The symbols 0 show our predicted FP-MD results.
Cu NWs, we plotted the depression of the melting temperature Tm versus the reciprocal diameter 1/D of Cu NWs in Figure 5b. The agreement shown in Figure 5 indicates that the COMPASS method is suitable to simulate the Cp(T) plots and the potential energy with different diameters, but it needs to be modified to correspond to the FP-MD results. Equation 3 can satisfactorily describe the size-dependent melting temperature of metallic Cu NWs. As shown above, the Tm(D) function for Cu NWs can be described as long as the related thermodynamic parameters are already known. The physical nature of the suppression of the physical property is related to both the surface/volume ratio and the elastic energetic change.40 4. Conclusion In summary, we have studied the melting process and transport properties of Cu NWs by different simulation methods, FP-MD and classical MD methods with the COMPASS potential. By analyzing the results of the potential energy, Cp(T), and the present simulation power, COMPASS is a compromise method between FP-MD and EAM. According to the difference of 10% between COMPASS and FP-MD for both bulk and Cu NWs with a diameter of 2 nm, the corresponding FP-MD melting points Tm(D) of Cu NWs with a diameter of 2-6 nm are predicted indirectly. The quantum conduction G of Cu NWs with D ) 2 nm between 420 and 480 K decreases obviously, which corresponds to our calculated melting point. Meanwhile, our COMPASS simulations and predicted FP-MD results of Cu NWs are validated by the size dependence of the Tm(D) function and in accordance with other simulation results. References and Notes (1) Pullini, D.; Innocenti, G.; Busquets, D.; Ruotolo, A. Appl. Phys. Lett. 2007, 90, 133106.
Zhang and He (2) Locharoenrat, K.; Sano, H.; Mizutani, G. Sci. Technol. AdV. Mater. 2007, 8, 277. (3) Hwang, H. J.; Kang, J. W. Surf. Sci. 2003, 532, 536–535. (4) Li, H.; Pederiva, F.; Wang, B. L.; Wang, J. L.; Wang, G. H. Appl. Phys. Lett. 2005, 86, 011913. (5) Wen, Y. H.; Zhu, Z. Z.; Zhu, R. Z.; Shao, G. F. Physica E 2004, 25, 47. (6) Li, H.; Wang, B. L.; Wang, J. L.; Wang, G. H. Chem. Phys. Lett. 2004, 399, 20. (7) Li, H.; Wang, B. L.; Wang, J. L.; Wang, G. H. J. Chem. Phys. 2004, 121, 8990. (8) Wang, B. L.; Wang, G. H.; Chen, X. S.; Zhao, J. J. Phys. ReV. B 2003, 67, 193403. (9) He, C.; Zhang, P.; Zhu, Y. F.; Jiang, Q. J. Phys. Chem. C 2008, 112, 9045. (10) Ao, Z. M.; Jiang, Q. Langmuir 2006, 22, 1241. (11) Ouyang, G.; Li, X. L.; Tan, X.; Yang, G. W. Nanotechnology 2008, 19, 045709. (12) Liu, Z.; Sakamoto, Y.; Ohsuna, T.; Hiraga, K.; Terasaki, O.; Ko, C. H.; Shin, H. J.; Ryoo, R. Angew. Chem., Int. Ed. 2000, 39, 3107. (13) Sun, C. Q.; Wang, Y.; Tay, B. K.; Li, S.; Huang, H.; Zhang, Y. B. J. Phys. Chem. B 2002, 106, 10701. (14) Li, S.; White, T.; Sun, C. Q.; Fu, Y. Q.; Plevert, J.; Lauren, K. J. Phys. Chem. B 2004, 108, 16415. (15) Sun, C. Q.; Bai, H. L.; Li, S.; Tay, B. K.; Jiang, E. Y. Acta Mater. 2004, 52, 501. (16) Hendy, S. C. Nanotechnology 2007, 18, 175703. (17) Wang, L.; Zhang, Y. N.; Bian, X. F.; Chen, Y. Phys. Lett. A. 2003, 310, 197. (18) Wang, B. L.; Zhao, J. J.; Chen, X. S.; Shi, D.; Wang, G. H. Nanotechnology 2006, 17, 3178. (19) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (20) Schwegler, E.; Galli, G.; Gygi, F.; Hood, R. Q. Phys. ReV. Lett. 2001, 87, 265501. (21) Hahn, P. H.; Schmidt, W. G.; Seino, K.; Preuss, M.; Bechstedt, F.; Bernholc, J. Phys. ReV. Lett. 2005, 94, 037404. (22) Wang, S. W.; Cao, Y. Z.; Rikvold, P. A. Phys. ReV. B 2004, 70, 205410. (23) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892. (24) Sharp, K. A.; Madan, B. J. Phys. Chem. B 1997, 101, 4343. (25) Nose, S. J. Chem. Phys. 1984, 81, 511. (26) Haile, J. M. Molecular Dynamics Simulation; Wiley: New York, 1992. (27) Karayiannis, N. C.; Mavrantzans, V. G.; Theodorou, D. N. Macromolecules 2004, 37, 2978. (28) Soldera, A. Polymer 2002, 43, 4269. (29) Mares, A. I.; van Ruitenbeek, J. M. Phys. ReV. B 2005, 72, 205402. (30) Jiang, Q.; Shi, H. X.; Zhao, M. J. Chem. Phys. 1999, 111, 2176. (31) Zhao, M.; Jiang, Q. Solid State Commun. 2004, 130, 37. (32) Yang, C. C.; Xiao, M. X.; Li, W.; Jiang, Q. Solid State Commun. 2006, 139, 148. (33) Lang, X. Y.; Han, L. P. J. Phys. Chem. C 2009, 113, 16036. (34) Jiang, Q.; Li, S. J. Nanosci. Nanotechnol. 2008, 5, 2346. (35) Regel’, A. R.; Glazov, V. M. Semiconductors 1995, 29, 405. (36) Jiang, Q.; Zhou, X. H.; Zhao, M. J. Chem. Phys. 2002, 117, 10269. (37) Zhang, Z.; Li, J. C.; Jiang, Q. J. Phys. D: Appl. Phys. 2000, 33, 2653. (38) Qi, W. H.; Wang, M. P. Mater. Chem. Phys. 2004, 88, 280. (39) Lu, H. M.; Han, F. Q.; Meng, X. K. J. Phys. Chem. B 2008, 112, 9444. (40) Wautelet, M. J. Phys. D: Appl. Phys. 1991, 24, 343. (41) Nanda, K. K.; Sahu, S. N.; Behera, S. N. Phys. ReV. A 2002, 66, 013208. (42) Prince, K. C.; Breuer, U.; Bonzel, H. P. Phys. ReV. Lett. 1988, 60, 1146. (43) Allen, G. L.; Bayles, R. A.; Giles, W. W.; Hesser, W. A. Thin Solid Films 1986, 144, 297. (44) Skripov, V. P.; Koverda, V. P.; Skokov, V. N. Phys. Status Solidi A 1981, 66, 109.
JP910246T
