# Fileset

[manuscript-20240216 R3.docx](https://mdr.nims.go.jp/filesets/2bc6ce7e-5542-408a-be63-12e8008311eb/download)

## Creator

Kongping Wu, Leng Zhang, Fangzhen Li, [Liwen Sang](https://orcid.org/0000-0003-0946-1025), [Meiyong Liao](https://orcid.org/0000-0003-1361-4266), Kun Tang, Jiandong Ye, Shulin Gu

## Rights

[Creative Commons BY-NC-ND Attribution-NonCommercial-NoDerivs 4.0 International](https://creativecommons.org/licenses/by-nc-nd/4.0/)

## Other metadata

[Enhancement of interfacial thermal conductance by introducing carbon vacancy at the Cu/diamond interface](https://mdr.nims.go.jp/datasets/225ea149-a6e9-4936-bbfb-1eb6f5c088d4)

## Fulltext

Enhancement of interfacial thermal conductance by introducing carbon vacancy at the Cu/diamond interfaceKongping Wu[footnoteRef:1],†, Leng Zhang1, Fangzhen Li1, Liwen Sang2, Meiyong Liao[footnoteRef:2],‡, Kun Tang3, Jiandong Ye3, Shulin Gu3† Corresponding author. Tel: +86(025) 8618-8572. E-mail address: kpwu@jit.edu.cn (K. P. Wu).‡ Corresponding author. Tel: +81(29) 860-4508. E-mail: Meiyong.Liao@nims.go.jp (M. Y. Liao).1School of Electronics and Information Engineering, Jinling Institute of Technology, Nanjing, Jiangsu, 211169, China2 National Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0044, Japan3School of Electronic Science and Engineering, Nanjing University, Nanjing Jiangsu, 210093, ChinaAbstractEffective heat dissipation of semiconductor devices is crucial for their extended lifespan and operational stability, and the interface of semiconductors provides an effective window for thermal design and management. In this work, we have systematically investigated the effect of the carbon vacancy on the thermal conductivity of diamond and the interface thermal conductance (ITC) of Cu/diamond by using both first- principles calculation and molecular dynamics methods. Although the carbon vacancy leads to a decrease in the thermal conductivity of diamond, a marked increase in ITC from 37.98 MWm-2K-1 to about 177 MWm-2K-1 for diamond (1 1 1) plane and from 78.8 MWm-2K-1 to about 241 MWm-2K-1 for diamond (0 0 1) plane is observed between Cu and diamond with carbon vacancy. The increase of the ITC is mainly due to the anharmonic phonon scattering, revealed by the phonon density of states and phonon participation ratio. Keywords: Diamond, Thermal conductivity, Interfacial thermal conductance, First- principles simulation, Molecular dynamics calculationI. INTRODUCTION Miniaturization of size and rapid increase in power density of microelectronic devices have made heat dissipation become the most important challenge to stability and reliability of device performance, for example, the minimization of microelectronic devices [1] and the III-nitrides based electronic devices [12]. It is an important route for solving heat dissipation problem of semiconductor devices to enhance interfacial thermal conductance (ITC) [2-3]. Diamond has the highest thermal conductivity (more than 2000 Wm-1K-1) [4] and helps to improve heat dissipation. Combined with its high breakdown field and high electron mobility, it is considered as the extreme semiconductor for next-generation power devices [5-6]. If metal electrodes in electronic devices can serve as heat dissipation channels, they undoubtedly provide aa new facitilial  solution for the thermal management of diamond-based power electronics. A deep understanding of thermal transport at metal/diamond interfaces is crucial not only for reliability of the diamond-based power devices performance, but also for fundamental understanding and engineering application of thermal science, like the thermal management of power electronics.The metal/semiconductor interfaces not only appear in electrodes, but also widely present in junction field-effect transistors and Schottky barrier diodes, which provide opportunities and windows for the structure design and performance adjustment of semiconductor devices. According to the current existing research results, besides introducing intermediate layers in the interface region to enhance the ITC [7-9], the interfacial details also have large influence on the ITC, such as interface mismatch [10-11], bonds coupling strength [12-13], surface roughness [14-15], chemical components [16], and so on. Vacancy defects, one of the important point defects, widely exists in semiconductor bulks and their interfaces. Although it is also difficult to precisely determine the effect of vacancy defects on the thermal transport of semiconductors or semiconductor interfaces in engineering technology, simulation calculations at the atomic scale have a prominent advantage in predicting the vital role of vacancy defects in thermal transport [17].Recently, the interaction relationship between phonons and vacancy defect in semiconductors has attracted extensive attention [18-20]. Especially, due to the scattering effect of vacancy defects on phonon, the thermal conductivity of semiconductors is further limited or lowered [21-24]. However, few calculations on the vacancy at semiconductor interfaces are available to date [17, 25-26], especially, there is still no relevant report on the effect of vacancies on the thermal transport at the diamond-based interfaces. For Cu and diamond, they both have high thermal conductivity, and their lattices match well with relatively small thermal expansion. As composite materials [27-28], the ITC of the Cu/diamond (Cu/C) interface is the determining factor of the overall thermal conductivity. But the experimental results of ITC at the Cu/C interface varied significantly according their preparation methods. In recent publication work by Fan et al. [29], they reported an experimental values of 55 MW/(m2·K) for the (100) diamond/Cu and 48 MW/(m2·K) for the (111) diamond/Cu was reported. Besides, the experimental values of 61.5 MW/(m2·K) [12], 66 MW/(m2·K) [8], 20 ~100 MW/(m2·K) [30], 35 MW/(m2·K) [31], 40~80 MW/(m2·K) [32] and 132 MW/(m2·K) [33] were also shown with large difference. The thermal transport across the interface at the nanoscale greatly depends on the interface details. Although vacancies are the most common point defects in these details, the effect of vacancy defects on thermal transport is less well-known. Particularly, a systematic study on the influence of vacancy defects on the interfacial thermal conductance between diamond and Cu is still sparse.In this work, we systematically examine how carbon vacancies (Vc) affect the thermal conductivity of diamond and the ITC of the Cu/C interface in order to realize the thermal transport in diamond-based devices and composites on the atomic scale. First, the thermal conductivity of diamond was calculated based on the first principles and molecular dynamics. Then, based on an abrupt interface model and non-equilibrium molecular dynamics (NEMD) simulation, the ITC of Cu/C interface was calculated and predicted. Finally, by comparing the calculated results with experimental data, we elucidated the physical mechanism behind the influence of the vacancy defect on thermal transport. The findings can serve as a theoretical reference for improving and enhancing the ITC of the Cu/C interface by engineering the Cu/C interface details.Fig. 1. Relationship between total energy and lattice constant for diamond without (a) and with Vc (b). (c) Relationship between total energy and lattice constant for Cu. (d) The interface structured Cu/C interface with fixed boundary conditions in the z direction in our NEMD simulations. Lattice structures of diamond without and with Vc and Cu are also shown in the inset. II. MODELS AND METHODSA. Models For diamond bulk, itDiamond crystal  has 8 C atoms in the conventional cell. One Vc can be created by removing a C atom from bulk diamond. The equilibrium lattice constant values of 3.572 and 3.631 Å for diamond and Cu bulks were reported in our recent work are selected [34]. The (0 0 1) Cu/C heterostructure was constructed using the supercell method as follows. A supercell of diamond (0 0 1) surface (5×5) with the thickness of 42 diamond conventional cell (consisting of 16800 C atoms) was built and a similar Cu supercell was b placed on the top of the diamond (0 0 1) surface supercell to minimize the effect of interface strain resulted from lattice mismatch on the ITC (a compression strain of the in-plane Cu lattice by 1.65%). Consequently, a rectangular box of 25.68×25.68×304.53 Å3 to simulate the (0 0 1) Cu/C heterostructure was successfully constructed. In order to obtain a rectangular box for the (1 1 1) Cu/C heterostructure, 1×√3 diamond (1 1 1) surface was firstly used, then a supercell of diamond (1 1 1) surface (5×5√3) with the thickness of 74 diamond conventional cell (consisting of 29600 C atoms) was built, and a similar Cu supercell was placed on the top of the diamond (1 1 1) surface supercell. B. Computational Methods   All first-principles calculations in this article are performed based on VASP (Vienna Ab-initio Simulation Package) software with the PBE (Perdew-Burke-Enzerhof) potential function [35-36]. To have a convergence, a truncated energy of 500 eV was set and a 9×9×9 Monkhorst-Pack mesh was is used in the first-principles calculations. In the process of structural optimization, the convergence condition is set to be that the forces on atoms is less than 1×10-6 eV/Å. The density functional perturbation theory and the finite displacement method are used to calculate the second-order and third-order force constants, respectively [37-38]. The third-order force constant is calculated using the 2×2×2 supercell with a cutoff radius of 6.0 Å, which is sufficient for obtaining well-converged phonon dispersion relationship. Based on the second-order and third-order force constants, the phonon properties and thermal conductivity of the diamond were are obtained by solving the phonon Boltzmann transport equation (BTE) using the ShengBTE code [39-40]. It should be noted that when exploring low concentration of the Vc in bulk diamonds bulk, molecular dynamics methods was is used, mainly because first-principles calculation for the thermal conductivity of larger systems cannot be supported by the computational power of the supercomputer.Besides, the nonequilibrium molecular dynamics (NEMD) approach is also used to calculate the ITC of the Cu/C interfaces. The simulations of thermal transport along the z axis perpendicular to the interface are performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software package [41]. In the NEMD simulation, the Tersoff [42-43] potential function is adopted to describe the atomic interactions in diamond bulk, and the modified embedded atom Method (MEAM) [44-45] interatomic potential is used to describe the atomic interactions between Cu and C at the interfaces. The metal/semiconductor interface thermal transport can be described by the dual-temperature model proposed by Majumdar-Reddy [46]. According to the model, the electron-phonon coupling has little contribution to the interface thermal conduction below 500 K [47-48]. For example, for the interfacial thermal conductance at the Al/Si interface, electron-phonon coupling contributes less than 9% of the total interfacial thermal conductance [49]. Therefore, phonon-phonon interaction is considered to play a decisive role in interface thermal conductance at the Cu/diamond interfaces.The three directions simulated in this article all adopt periodic boundary conditions. So the phonon boundary scattering effect is not considered. The atoms at the two ends of the heterostructure are fixed along the heat conduction direction. The heat source (320 K) is on the left and the  heat sink (280 K) is on the right, and the time step is 0.001 ps. Before simulation, static equilibrium is firstly performed at 300 K and the conjugate gradient method is used to minimize the energy of the system. Next, the system is dynamically optimized under a microcanonical ensemble (NVE) at 300 K, and h Langevin heat bath method for 4 million steps of relaxation is used to replace the NVE. It was found that in NEMD simulations, Langevin heat bath method had more advantages in controlling local temperatures than Nosé-Hoover method [50]. After relaxation, 2 million steps in Langevin heat bath method run continuously, and the system reaches a steady state with a linear temperature distribution along the direction of thermal transport. The system was is divided into 22 segments in the direction of thermal transport to calculate the temperature of every segment. The Ttemperature of each segment is recorded every 10000 steps and the data sample collects the average temperature once every 1000 steps. At the same time, the energy and temperature gradient were are obtained to calculate the ITC. The same calculation is performed by many times with changing the random number in the simulation to ensure the reliability of the simulation results, and the ITC is obtained according to the average value.III. RESULTS AND DISCUSSIONA. Thermal conductivity of diamond without and with VcAccording to the equilibrium molecular dynamics (EMD) method, the thermal conductivity of diamond can be calculated by the Green-Kubo expression [51-52].                  (1)Where where kB, T and V are Boltzmann constant, temperature 300K and the system volume, respectively. And <Jz(0) Jz(t)> presents the heat flux autocorrelation function (HFACF) along the thermal transport direction. The  changes over time and is usually called running thermal conductivity [53-54]. Fig. 2. (a) The relationship between the heat flux autocorrelation function (HFACF) and running thermal conductivity and correlation time for diamond at 300 K. (b) The relationship between energy of the heat source and sink region and time in steady state, (c) temperature profile with the diamond length along z direction, and (d) the relationship between the reciprocal of thermal conductivity and the reciprocal of the bulk diamond model length (1/Lz) in the MD simulations. The corresponding linear fit is also shown in these figures.In EMD simulation, the thermal conductivity of bulk crystal is heavily dependent on the size of the model. An extrapolation method was proposed by Schelling [55] to obtain the thermal conductivity of infinite systems, as follow.                    (2)where  is the value of  when Lz equals infinity. The simulation was is carried out for diamond with different sizes through the supercell method, and the 5×5×22, 5×5×28, 5×5×45 and 5×5×50 supercells of the diamond conventional cell contain 4400, 5600, 9000 and 10000 C atoms, respectively. In EMD simulation for 5×5×50 diamond supercell, its cross-sectional area and length along the thermal transport direction are about 318.98 Å2 and 178.6 Å, respectively. The HFACF and running thermal conductivity (κrunning) changing with the correlation time were are shown on the left of the panel in Fig. 2(a). The HFACF eventually converges with increasing the correlation time. After calculating many times, the average value of HFACF is represented by a black line in Fig. 2(a). The corresponding values of were shown on the right of the panel using blue line in Fig. 2(a). For diamond, the calculated thermal conductivity value of 1589±12 Wm-1K-1 is lower than the experimental values of 2000 Wm-1K-1, which is mainly because the size of the 5×5×50 diamond supercell is not large enough. Similarly, the calculated thermal conductivity values of the other three diamond supercells 5×5×22, 5×5×28 and 5×5×45 are 1255±7 Wm-1K-1, 1341±10 Wm-1K-1 and 1577±17 Wm-1K-1, respectively.In the NEMD simulation, a larger diamond supercell 5×5×80 was is used in order to eliminate the effect of finite supercell size,. the The thermal conductivity can beis calculated by Fourier’s law as [56]. Here,  is the heat flux, which is defined as the energy transfer rate through cross-sectional area (S) in the form of  . And the ▽T is the temperature gradient regarding to geometric size along the direction of thermal transport in the process of energy transfer. According to the Langevin heat bath method [57], the temperature values of heat source and heat sink are set to be 325 K and 275 K, respectively. The energy varies over time at the position of heat source and heat sink as shown in the Fig. 2(b). It can also be seen in the Fig. 2(b) that this functional relationship conforms to linear fitting. The energy transfer rate is the average of the absolute values of the slopes of two fitted straight lines. Besides, The temperature profile along the direction of thermal transport is shown in the Fig. 2(c). Similarly, it can also be linearly fitted, which indicates the simulated system is in a steady-state during the imposed heat flux process. The temperature gradient ▽T is achieved by taking the slope of the fitted line. According to Fourier’s law, the thermal conductivity of this system is about 1916.41 Wm-1K-1, which is very close to the experimental value of 2000 Wm-1K-1. Similarly, the calculated thermal conductivity values of the other three diamond supercells 5×5×100, 5×5×120 and 5×5×150 are 1989.57 Wm-1K-1, 2021.45 Wm-1K-1 and 2088.11 Wm-1K-1, respectively.All values calculated for thermal conductivity are shown in the Fig. 2(d). According to extrapolation method (Lz→∞), we obtain the thermal conductivity values of 2010.25 Wm-1K-1 for diamond by EMD method and 2339.45 Wm-1K-1 for diamond by NEMD, respectively. The calculated values are in good agreement with the experimental values of 2000 Wm-1K-1. For the thermal conductivity of diamond, there were some other calculated values reported recently, like 1950 Wm-1K-1 by EMD method using the Tersoff potential [58], 840 Wm-1K-1 by EMD method using the anharmonic Brenner potential [51], and 2740 Wm-1K-1 by NEMD method using the Tersoff potential [59]. These computational results coincide with our calculating results very well. The diamond supercells 5×5×80, 5×5×100, 5×5×120 and 5×5×150 contain 16000, 20000, 24000 and 30000 C atoms. The Vc can be generated by removing C atoms randomly according to the certain percentage of the total number of atoms from these supercells, and the concentration (nVc) of Vc is the ratio of the number of removed atoms to the total number of atoms in the supercell. For the Vc concentration of 10-4 in 5×5×100 diamond supercells at 300 K, the energy varies over time at the position of heat source and heat sink as shown in the Fig. 3(a). The energy increases at the heat sink region, while decreases at the heat source region. And its temperature profile along the direction of thermal transport is also shown in the Fig. 3(b). According to thermal conductivity , the calculated value is listed in Fig. 3(c). Similarly, the values of thermal conductivity for the other three supercells 5×5×80, 5×5×100, 5×5×120 and 5×5×150 were are also calculated, as shown in Fig. 3(c). Finally, based on the extrapolation method, the thermal conductivity of diamond with the Vc concentration of 10-4 is about 1337 Wm-1K-1. The value is slightly bigg largerer than the experimental value of 1206 Wm-1K-1 [60] and the theoretical value of 1100 Wm-1K-1 [21, 22]. The values of thermal conductivity for diamond with the other defect Vc concentrations were are also calculated using a similar method, as shown in the Fig. 3(d). The highest Vc concentration in the diamond is 12.5%, and itwhich can be generated by removing one C atom randomly from diamond conventional cell with eight C atoms. Thermal conductivity of the diamond with 12.5% Vc concentration can be obtained through ShengBTE code based on a full iterative solution to the Boltzmann transport equation [39-40]. The calculated value is about 14.15 Wm-1K-1 and also shown in the the Fig. 3(d).FIG. 3. (a) The relationship between energy of the heat source and sink region and time in steady state, (b) temperature profile with the diamond length along z direction, (c) the relationship between the reciprocal of thermal conductivity and the reciprocal of the length (1/Lz) for 5×5×100 diamond supercell with the Vc concentration of 10-4 in the NEMD simulations at temperature of 300 K. The corresponding linear fit is also shown in the figure. And (d) the relationship between the thermal conductivity of diamond and the Vc concentration (nVc) at temperature of 300 K. The first-principles calculation provides fundamental understanding of the effect of the Vc on the thermal conductivity of diamond. Based on the 2nd-order and 3rd-order force constants, the phonon properties and thermal conductivity of diamond with and without Vc are obtained using ShengBTE code. Convergence testing was is firstly performed on scalebroad, which has a significant impact on computational speed [40]. According to our testing results shown in the inset of Fig. 4(a), the value of scalebroad is taken as 0.4, which is large enough and does not to affect the calculated result. The thermal conductivity of diamond with and without the Vc has been calculated from 50 K up to about 1350 K and is shown in the Fig. 4(a). The values of thermal conductivity for diamond with and without the Vc decrease with increasing temperature as are shown in the Fig. 4(a). Besides, the thermal conductivity of diamond with the Vc is significantly lower than that of diamond without the Vc, which indicates that the Vc in diamond has a strong suppression of thermal conductivity of diamond. In the Fig. 4(b), two transverse acoustic modes (TA1 and TA2) and a longitudinal acoustic mode (LA) in the C7-Vc have moved to a lower frequency position compared to the pristine diamond (C8). Besides, the effective anharmonic scattering rate are is estimated that isto be several hundred times larger than that in the pristine diamond. This is expected to a large suppression ofsuppress the thermal conductivity of diamond [22]. It is mainly because that the phonon-Vc anharmonic scattering rates in diamond with the Vc are much larger than the phonon-phonon scattering rates in diamond without the Vc at the lower frequency range. FIG. 4. (a) The relationship between the thermal conductivity of diamond without (C8) and with Vc (C7-Vc) and temperature. (b) Anharmonic scattering rates in diamond without (C8) and with (C7-Vc) the VC at temperature of 300 K. B. Interfacial thermal conductance It can be understood that the Vc intrinsic defect induces a decrease in the thermal conductivity of diamond due to the enhanced local the phonon-VC anharmonic scattering by the Vc defects. However, it was reported recently that the presence of the intrinsic vacancy defect enhanced interfacial thermal transport [25,26]. The Cu/C heterostructure has important application value in the field of power devices, and the effect of the intrinsic Vc defect on interfacial thermal conductance at the Cu/C interface is still unclear. The values of ITC at the (0 0 1) Cu/C and (1 1 1) Cu/C interface were are calculated by NEMD method using the MEAM interatomic potential, which is also appropriate for calculating ITC at Cu/C interface with the intrinsic Vc defect [17].  The model systems, including two types of interfaces: the (0 0 1) Cu/C and (1 1 1) Cu/C interface, are constructed and calculated. In the process of minimizing system energy, the interface energy can be calculated using the formula as the follow [61]:                         (3)where is the total potential energy of the model systems, and are the cohesive energy of C and Cu, which can be obtained from the relaxed diamond and copper bulks. The calculated values of and are −7.148 eV/atom and −3.537 eV/atom, which is in good agreement with the experimental values −7.371 eV/atom for C [62] and −3.915 eV/atom for Cu [63]. andare the number of Cu atoms and C atoms in the model system, respectively. Sint is the cross-section area of the model system. Besides, the interface space is also obtained from the relaxed model systems. The calculated interface energy and space were listed in Table 1. The interface energy of 4.13 J/m2 at the (0 0 1) Cu/C interface is larger than the interface energy of 3.26 J/m2 at the (1 1 1) Cu/C interface, and the interface space (1.55 Å) of the former was alsoalso is  smaller than that (1.46 Å ) of the latter based on the analysis and calculation of the position coordinates of interface atoms. Table 1. Interface energy, interface space and ITC at the calculated model systems. Interface types Interface energy (J/m2) Interface space (Å) ITC(MWm-2K-1) ITC (MWm-2K-1) in Literatures (0 0 1) Cu/C 4.13 1.55 78.88 (Exp.) 55 in Ref. [29]     (Exp.) 66 in Ref. [8]     (Theo.) 72 in Ref. [29] (1 1 1) Cu/C 3.26 2.28 37.98 (Exp.) 48 in Ref. [29]     (Theo.) 78 in Ref. [29]“Exp.” is the abbreviation of “the experimental value”“Theo.” is the abbreviation of “the theoretical calculation value”For The ITC at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces, it can be , is calculated by NEMD method using the formula as the follow [34,64]:                   (4)where J is defined as the average heat flux along the thermal transport, which can beis obtained from the average of the absolute values of the fitted slope. And ΔT is the temperature jump at the interface. For the (0 0 1) Cu/C interface, Fig. 5(a) shows the energy of the the (0 0 1) Cu/C interface varies varying over time at the heat source region and heat sink region. And temperature jump (ΔT) at the interface is shown in the Fig. 5(b), where the interface is denoted with a dashed line. Similar calculations also apply to the (1 1 1) Cu/C interface. The relationship between the energy and time at the heat source region and heat sink region and the temperature jump (ΔT) at the interface are shown in the Figs. 5(c) and 5(d). respectively. The calculated values of ITC for the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces are 78.88 MWm-2K-1 and 37.98 MWm-2K-1, respectively. Both values of ITC are also included in the table 1.FIG. 5. (a) The relationship between energy of the heat source and sink region and time in steady state, and (b) temperature profile with the system length along z direction for the (0 0 1) Cu/C interface. (c) The relationship between energy of the heat source and sink region and time in steady state, and (b) temperature profile with the system length along z direction for the (1 1 1) Cu/C interface at temperature of 300 K. The corresponding linear fit is also shown in the figure. The calculated ITC of 78.88 MWm-2K-1 at the (0 0 1) Cu/C interface is slightly higher than those of reported experimental values of 55 MWm-2K-1 [29] and 66 MWm-2K-1 [8], as well as . It is also very close to the theoretical calculation value of 72 MWm-2K-1 [29] calculated by Lennard-Jones potential. For the (1 1 1) Cu/C interface, its the ITC of 37.98 MWm-2K-1 is only about one half of the calculation value of 78 MWm-2K-1[29], while it is in good agreement with the experimental values of 35 MWm-2K-1 reported by Monachon et al. [7, 31], as well as slightly smaller than the experimental values of 48 MWm-2K-1 reported by Yang [29]. The theoretical calculated results demonstrate that the ITCs of Cu/C interface is are dependent on the diamond orientation to a certain extent. One of the main reasons is that compared to the (1 1 1) Cu/C interface, there is a stronger intensity of interfacial coupling at the (0 0 1) Cu/C interface with a bigger higher interface energy and a shorter interface space. Therefore, the transmission of interface phonons and further promote the interfacial thermal transport is enhancedced [35, 65,66]. Fig. 6 (a) The dependence of ITC on the distance from the position of the Vc to diamond surface (dz), and (b) the dependence of ITC on the concentration of the Vc (nVc) at the (0 0 1) Cu/C interface and the (1 1 1) Cu/C interface.The intrinsic defect Vc is widely present in diamond-based semiconductor heterostructures, and clarifying its impact on the ITC is of great important for thermal management of diamond-based semiconductor devices. Firstly, it is necessary to explain localization properties of Vc in diamond. The atom projected phonon density of states (PPDOS) in C8 and C7-Vc are is shown in Fig. S1 in the Supplementary Material. And the first (1st-C) and second nearest (2nd-C) C atoms around Vc site are distinguished by gray and blue as shown in Fig. 1(b). Compared to the PPDOS of C atoms in C8, the PPDOS of C atoms in C7-Vc are is significantly enhanced in the low-frequency region (less than 15 THz). Additionally, the PPDOS of the 1st-C atom is noticeably stronger than that of the 2nd-C atoms. This also indicates that the influence of the Vc on phonon properties is localized. At the same time, in order to precisely regulate Vc in the interface region, we first need to determine which parts of the diamond in interface region. Since we only examine the intrinsic defect Vc here, only one Vc is moved from the surface of the diamond towards the interior. The distance from the position of the Vc to the surface of the diamond is denote by dz as shown in Fig. 1(d). Then the dependence of ITC on the dz is examined in the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces as shown in Fig. 6(a) to further determine the thickness of the interface region on the diamond side.It can be found in Fig. 6(a) that the ITC at the interface gradually tends to a stable value as the intrinsic defect Vc in diamond moves away from the interface. The stable values of the ITC are about 92.24 MWm-2K-1 and 40.15 MWm-2K-1 at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces, respectively. They are  slightly larger than 78.88 MWm-2K-1 and 37.98 MWm-2K-1 at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces without the Vc, as shown in the Fig. 6(a). Besides, the critical values of distance (dz) corresponding to the stable ITC are about 17.88 Å and 8.86 Å in the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces, respectively. Based on the distance (dz), 1800 and 800 C atoms are contained in the side of diamond within the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces. Among these atoms, randomly deleting C atoms according to a certain percentage produces the defect Vc as shown in Fig. 6(b). Following the utilization of four distinct random numbers to establish a 0.6% vacancy ratio, the system undergoes relaxation to reach its minimum energy state. The testing procedure is shown in Figs. S2-S5 in the Supplementary Material. Consequently, the discrepancy between the maximum and minimum values of ITCs is about 5.3 MWm-2K-1. All the values of ITC at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces with concentration of the Vc (nVc) are calculated and shown in the Fig. 6(b). The Vc far from the interface causes a slight increase in the value of ITC from 78.88 MWm-2K-1 to 92.24 MWm-2K-1 for the (0 0 1) Cu/C interface and from 37.98 MWm-2K-1 to 40.15 MWm-2K-1 for the (1 1 1) Cu/C interface, respectively. Compared to the ITC at Cu/C interface with Vc far from the interface, the Vc in interface region has a more significant effect in enhancing the ITC. Therefore,  the effect of Vc on ITC has strong localization characteristics at the Cu/C interface. Besides, at lower concentrations nVc, as the percentage of Vc increases, the ITC increases more significantly as shown in the Fig. 6(b). The maximum values of the ITC at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces with the Vc are about three times and over four times higher than those at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces without the Vc, respectively.To reveal the physics mechanism of the enhanced ITC resulted from the Vc, it is necessary to understand the transmittance of interface phonons and the phonon density of states at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces [67,68,69,70]. Firstly, the phonon transmittance at the interface can be is calculated using two classical physical models [68]. According to the interface phonon transport models (diffuse mismatch model and acoustic mismatch model), the phonon transmittance can beis calculated using two formulas as the follow [68]:                   (5)and ,                         (6)where the parameterrepresents the transmission of phonons from the diamond to the metal Cu. V is the velocity of the phonon group and the subscript represents the corresponding phonon modes. DOS is the density of states of each vibration mode. ρ1 and ρ2 are the density of diamond and Cu, respectively [71]. The phonon group velocity (V) and density of states (DOS) of each vibration mode for diamond and Cu can baree obtained using the first-pricipleprinciples calculation combined with Phonopy code [72].Fig. 7. The phonon dispersions (left-half plane) and phonon density of states (PDOS) (right-half plane) of diamond (C8) (a), diamond with Vc (C7-Vc) (b) and Cu (c), and their corresponding group speed (d). The relationship between transmission of phonons from the diamond to the Cu and the phonon frequency from diffuse mismatch model (e) and acoustic mismatch model (f).The phonon dispersions and phonon density of states (PDOS) of diamond without Vc and with Vc and Cu are firstly calculated and shown in Figs. 7(a), 7(b) and 7(c), respectively. Our calculated results are in good agreement with the reported experimental results [73,74]. Then the corresponding group velocity are is obtained, as  and shown Fig. 7(d). The interaction between the same phonon modes at the interface is only considered in the theoretical model. The contribution of acoustic phonon to the ITC is also considered for the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces. Two transverse phonon modes (TA1 and TA2) and a longitudinal phonon mode (LA) are shown in the in Figs. 7(a), 7 (b) and 7(c), respectively. In addition, it is clearly observed that the Vc in diamond causes the cutoff frequency of acoustic phonon modes to decrease from about 20 THz to 12 THz, as shown in the Fig. 7(d). The cutoff frequency of acoustic phonon modes for Cu is less than 7.5 THz, and the PDOSs of diamond without Vc and with Vc under 7.5 THz are also shown in the inset of the right plane in the Fig. 7(a) and 7(b), respectively. It is revealed that the Vc in diamond causes an enhancement of PDOSs below the phonon frequency of 7.5 THz. The phonon transmittance coefficients are calculated using the DMM and AMM and shown in the Figs. 7(e) and 7(f), respectively. By comparison comparing with perfect diamond, the phonon transmittance coefficient below the phonon frequency of 7.5 THz has been a certain increase in the diamond with Vc. The low phonon frequency is considered as  the main reason for the enhancement of ITC from 78.88 to 92.24 MWm-2K-1 at the (0 0 1) Cu/C interface and from 37.98 to 40.15 MWm-2K-1 at the (1 1 1) Cu/C interface, respectively.C. Phonon density of states analysisAs mentioned above, the ITCs at the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces strongly depend on the Vc in interface region ant its concentration, which is also considered to be correlated with the overlap of PDOSs for diamond and Cu in the interface region [26,75]. On the Cu side, a thickness of 20 Å near the interface is considered as a sufficiently long to accommodate all Cu atoms in the interface region.  Correspondingly, there are 88 and 324 Cu atoms in the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces region, respectively. The PDOS is an effective way to understand phonon activity in materials, which can be calculated based on the Fourier transform of the velocity autocorrelation function VACF(t) for all atoms as the follow[76]:,                      (7)where VACF(t) can be defined as . N is the number of atoms, and  and are the initial velocity of the i-th atom and its velocity at time of t, respectively. < > means taking the average of the system. At the Cu/C interfaces, all phonons that collide with the interface are elastically scattered only one time and then enter the other side with a probability proportional to the PDOSs. The degree of mismatch in the PDOSs of Cu and C atoms can be evaluated by the overlap factor S, which can beis described by the equation[26],                 (8)where PDOSC(ω) and PDOSCu(ω) are the PDOSs of C and Cu atoms in the Cu/C interfaces, respectively. The overlap factor S of the PDOS can be used to estimate phonon transmission across the Cu/C interfaces through elastic scattering.In addition, calculating phonon participation ratio (PPR) is another effective method to measure the proportion of atoms participating in a certain intrinsic vibrational mode. It is also reliable way to describe the phonon localization quantitatively. When phonon localization occurs, it is no longer a traveling wave mode in which all atoms participate. And the stronger the localization, the smaller the phonon participation rate. Besides, the PPR can be directly computed from the trajectory and velocity of atoms using NEMD method without lattice dynamic calculation and  can contain the contribution of anharmonic scatterings [76-77]. The phonon participating ratio PPR (ω) can be written as [78] :,             (9)where PDOS(ω)n is the PDOS of the n-th atom with a frequency of ω. And N is the number of C or Cu atoms within the interface region.The velocity autocorrelation function VACF(t) and the PDOSs for the Cu and C atoms in the (0 0 1) Cu/C interface region are calculated , and shownas shown in Fig. 8(a) and 8(b), respectively. Similarly, The VACF(t) and the PDOSs of the Cu and C atoms in the (1 1 1) Cu/C interface region are shown in Figs. 8(c) and 8(d), respectively. Through integrating the PDOSs, we found that the overlap factor S of the (0 0 1) Cu/C is nearly twice that of the (1 1 1) Cu/C interface. And the ITC listed in the tTable 1 at the (0 0 1) Cu/C is a bit more than around twice larger than that of the (1 1 1) Cu/C interface. It showsThis is consistency consistent with the statement that the ITC correlates strongly with the overlap of the PDOSs of the two interface materials [79].Fig. 8 The relationship between the velocity autocorrelation function and correlation time (a) and the PDOS and the frequency (b) for the (0 0 1) Cu/C interface, and the relationship between the velocity autocorrelation function and correlation time (c) and the PDOS and the frequency (d) for the (1 1 1) Cu/C interface.To further explore the effect of nCv on the ITC, the PDOSs of these Cu/C interfaces with the different nCv are calculated. The PDOSs of the (0 0 1) Cu/C and (1 1 1) Cu/C interfaces with the concentration of 0.1%, 0.75 and 1.5% are shown in Figs. 9(a)-9(c) and 9(d)-9(f), respectively. The ITCs corresponding to these interfaces have already been listed in the tTable 2. According to the Eq. (8), the overlap factor S for PDOSs of these Cu/C interfaces are calculated. The values of S for the (0 0 1) Cu/C interfaces with the concentration of 0.125%, 0.75 and 1.5% are 0.00327, 0.00307 and 0.0058, respectively. And the values of S for the (1 1 1) Cu/C interfaces with the concentration of 0.125%, 0.75 and 1.5% are 0.00267, 0.00308 and 0.00314, respectively. All these values are also listed in table Table 2. For the (0 0 1) Cu/C interface, the values of S don’t increase monotonously with the nCv increasing, and they are monotonically increasing, but their proportion relationship is completely inconsistent with that of the ITCs. Besides, their values of the interface space are calculated after optimization of these interface models. As the nCv increases, the interface space shows a monotonically decreasing trend, which means that the interaction force at these interfaces are increasing. It This is in good agreement with the statement that the stronger interaction at the interface will enhance the ITC. [34,65,66].Fig. 9 The relationship between PDOSs and phonon frequency at the (0 0 1) Cu/C interface with the concentration of 0.1% (a), 0.75 (b) and 1.5% (c). The relationship between PDOSs and phonon frequency at the (1 1 1) Cu/C interface with the concentration of 0.1% (d), 0.75 (e) and 1.5% (f).Table 2. Interface C vacancy concentration (nCv), interface space and ITC at the calculated model systems. the overlap factor S Interface types nCv Interface space (Å) ITC(MWm-2K-1) the overlap factor S (0 0 1) Cu/C 0.125 1.51 158.81 0.00327  0.75 1.49 219.56 0.00307  1.5 1.46 240.82 0.0058 (1 1 1) Cu/C 0.125 2.15 91.58 0.00267  0.75 2.03 144.69 0.00308  1.5 1.95 176.67 0.00314According to Eq. (9), the PPR is calculated and shown in Fig. 10 to evaluate the contribution of anharmonic scattering to the ITCs. The Figs. 10(a)-10(c) shows the PPR corresponding to the (0 0 1) Cu/C interfaces with the concentration of 0.1%, 0.75 and 1.5%, respectively. And the Figs. 10(d)-10(f) show the PPR corresponding to the (1 1 1) Cu/C interfaces with the concentration of 0.1%, 0.75 and 1.5%, respectively. The PPR value of 0.4 is regarded as the critical line represented by dashed line in these figures. The PPR value above the critical line is regarded as delocalized mode, while below the critical line is regarded as localized mode. And the delocalized mode is the main contributor to thermal transport [76, 78].Fig. 10 The relationship between phonon participation ratio (PPR) and phonon frequency at the (0 0 1) Cu/C interface with the concentration of 0.1% (a), 0.75 (b) and 1.5% (c). The relationship between the PPR and phonon frequency at the (1 1 1) Cu/C interface with the concentration of 0.1% (d), 0.75 (e) and 1.5% (f).Recently, it is reported that the optical modes in diamond with the higher frequency (lager than the cut-off frequency of Cu) is also the important heat-carrying phonons for contributing to the interfacial ITC through anharmonic scattering [75]. For the (0 0 1) Cu/C interfaces with the concentration of 0.1%, 0.75 and 1.5%, as the concentration of nVc increases, the localized model characteristics gradually strengthens at the lower frequency region (less than the cut-off frequency of Cu) as shown in the lower left of the Figs. 10(a)-10(c). From the first-principles calculation, although abundant generation of non-harmonic phonons in low-frequency regions occurs , shown in Fig. 4(b), these phonons exhibit strong localized characteristics due to the scattering interaction between non-harmonic phonons and Vc. It is well in line with the reported by Loh et. al [78] that the acoustic modes localized at locations nearest the vacancy, suppressing the thermal transport at the (0 0 1) Cu/C interface. In contrast, the PPR above frequency of 20 THz (corresponding to the optical modes of diamond) has a significant increase with the increasing of the nVc. These phonons are the main contributor to interface thermal transport through the anharmonic scattering across these Cu/C interfaces. The localization enhancement of phonon at low-frequency region and the delocalization enhancement of phonon at high-frequency region also occur in the (1 1 1) Cu/C interfaces with the concentration of 0.1%, 0.75 and 1.5% as shown in the Figs. 10(d)-10(f), respectively. The contribution of anharmonic scattering to the ITCs is also from the scattering channels including the phonon-Vc scattering and phonon-boundary scattering at these interfaces. These anharmonic scattering plays a major role in enhancing interface thermal transport [49, 80]. The contribution increased to the ITC through anharmonic phonon scattering not only makes up for the decrease due to localization at lower-frequency region, but also leads to the increase of ITC with the increase of nVc. It is expected that the introduction of the Vc on diamond surface would help improve the electronics devices power handling with higher reliability [81,82]. Finally, we note that our calculation may be adopted for different metals/diamond junctions. However, one need to consider the specific Debye temperature of each metal, the acoustic mismatch at the junction, interface bonds and lattice mismatch, etc.  IV. CONCLUSIONThe effect of the Vc on the thermal conductivity of diamond and the ITC at its interface with Cu were systematically investigated using both first-principles calculation and molecular dynamics methods. According to the first-principles calculation, a substantial decrease was observed in the thermal conductivity of diamond with Vc compared to the pristine diamond. This suggests a strong suppressive effect of Vc on the thermal conductivity of diamond. It is mainly because that the presence of Vc significantly boosts anharmonic scattering rates in diamond, particularly at lower frequencies, far exceeding the phonon-phonon scattering in pristine diamond. For the Cu/C interface, the calculated ITC of 78.88 MWm-2K-1 at the (0 0 1) Cu/C interface and 37.98 MWm-2K-1 at the (1 1 1) Cu/C interface are in good agreement with the experimental values of 66 MWm-2K-1 and 35 MWm-2K-1, respectively. The introduction of Vc defect at interfaces leads to a significant increase in the ITC. Based on the PDOS and PPR calculation, it was verified that anharmonic scattering played a major role in enhancing the interface thermal transport. This work helps understand the fundamentals of the effect of the intrinsic defects (Vc) on the thermal conductivity of diamond and the ITC at metal/diamond interface and provides a guideline for thermal management of the diamond-based electronic devices and composite.ACKNOWLEDGMENTSThis work is financially supported by the Science Foundation of Jinling Institute of Technology (jit-rcyj-202001) and the scholarship from China Scholarship Council (CSC) under the Grant No.202308320197. We thank the technology support from the National Supercomputer Center (TianHe-2) in Lvliang and Dawning Intelligent Computing AC Platform (KunShan) for providing computer resources.References [1] Y. R. Wang, D. Gao, T. Zhang, H. Zhang, Y. Zhang, Q. Fu, H. Zhao, Z. B. Ma, Functional Diamond, 3, 2183097 (2023).[2] J. L. Smoyer and P. M. Norris, Heat Transfer Eng., 40, 269 (2019).[3] X. D. Zhang, G. Yang, and B. Y. Cao, Adv. Mater. Interfaces, 9, 2200078 (2022).[4] D. G. Onn, A. Witek, Y. Z. Qiu, T. R. Anthony, and W. F. Banholzer, Phys. Rev. Lett., 68, 2806 (1992).  [5] H. Umezawa, M. Nagase, Y. Kato, S. Shikata, Diam. Relat. Mater., 24, 201 (2012).[6] N. Donato, N. Rouger, J. Pernot, G. Longobardi, F. Udrea, J. Phys. D: Appl. Phys., 53, 093001 (2020).[7] G. Chang, F. Sun, L. Wang, Z. Che, X. Wang, J. Wang, M. J. Kim, H. Zhang, ACS Appl. Mater. Interfaces, 11, 26507 (2019).[8] G. Chang, F. Sun, J. Duan, Z. Che, X. Wang, J. Wang, M. J. Kim, H. Zhang, Acta Mater., 160, 235 (2018).[9] M. Blank, L. Weber, J. Appl. Phys., 124, 105304 (2018).[10] R. J. Stevens, L. V. Zhigilei, P. M. Norris, Int. J. Heat Mass Transf., 50, 3977 (2007). [11] Q. S. Li, F. Liu, S. Hu, H. F. Song, S. S. Yang, H. L. Jiang, T. Wang, Y. K. Koh, C. Y. Zhao, F. Y. Kang, J. Q. Wu, X. K. Gu, B. Sun, X. Q. Wang, Nat. Commun., 13, 4901 (2022).[12] B. Xu, S. Q. Hu, S. W. Hung, C. Shao, H. Chandra, F. R. Chen, T. Kodama, J. Shiomi, Sci. Adv., 7, eabf8197 (2021).[13] N. Q. Chen, K. M. Yang, Z. Y. Wang, B. A. Zhong, J. J. Wang, J. Song, Q. Li, J. M. Ni, F. Y. Sun, Y. Liu, T. X. Fan, ACS Appl. Mater. Interfaces, 15, 34132 (2023)[14] Y. Zhang, D. K. Ma, Y. Zang, X. Wang and N. Yang, Front. Energy Res., 6, 48 (2018).[15] J. W. Zhao, R. Zhao, Y. K. Huo, W. L. Cheng, Int. J. Heat Mass Transf., 140, 705 (2019).[16] L. N. Yang, X. Wan, D. K. Ma, Y. Jiang and N. Yang, Phys. Rev. B, 103, 155305 (2021).[17] Z. X. Lu, A. M. Chaka, and P. V. Sushko, Phys. Rev. B, 102, 075449 (2020).[18] N. Mingo, K. Esfarjani, D. A. Broido, and D. A. Stewart, Phys. Rev. B, 81, 045408 (2010).[19] T. Wang, J. Carrete, A. van Roekeghem, N. Mingo, and G. K. H. Madsen, Phys. Rev. B, 95, 245304 (2017).[20] C. A. Polanco and L. Lindsay, Phys. Rev. B, 98, 014306 (2018).[21] N. A. Katcho, J. Carrete, W. Li, and N. Mingo, Phys. Rev. B, 90, 094117 (2014).[22] N. H. Protik, J. Carrete, N. A. Katcho, N. Mingo, and D. Broido, Phys. Rev. B, 94, 045207 (2016). [23] A. Katre, J. Carrete, T. Wang, G. K. H. Madsen, and N. Mingo, Phys. Rev. Mater., 2, 50602 (2018).[24] Z. Q. Yan and S. Kumar, Phys. Chem. Chem. Phys., 20, 29236 (2018).[25] Y. Y. Zhang, X. Qian, Z. Peng, N. Yang, MRS Online Proceedings Library, 1753, 7 (2015).[26] X. J. Liu, G. Zhang, and Y. W. Zhang, Nano Lett., 16, 4954 (2016).[27] S. G. Daia, J. W. Li, N. X. Lu, Diam. Relat. Mater., 108, 107993 (2020).[28] L. Lei, L. Bolzoni, F. Yang, Carbon, 168, 553 (2020)[29] K. M. Yang, Z. Y. Zhang, H. H. Zhao , B. H. Yang , B. Zhong , N. Q. Chen , J. Song , C. Chen , D. W. Tang, J. Zhu , Y. Liu , T. X. Fan, Acta Mater., 220, 117283 (2021).[30] C. Monachon, L. Weber, C. Dames, Annu. Rev. Mater. Res., 46, 433 (2016).[31] C. Monachon, L. Weber, Emerging Mater. Res., 1, 89 (2012).[32] S.-W. Hung, S. Hu, J. Shiomi, ACS Appl. Electron. Mater., 1, 2594 (2019).[33] G. Chang, F. Sun, L. Wang, Y. Zhang, X. Wang, J. Wang, M. J. Kim, H. Zhang, Compos. Part A, 135, 105921 (2020).[34] K. P. Wu, L. Zhang, D. B. Wang, F. Z. Li, P. Z. Zhang, L. W. Sang, M. Y. Liao, K. Tang, J. D. Ye, S. L. Gu. Sci. Rep., 12, 19907 (2022).[354] G. Kresse and J. Furthm¨uller, Phys. Rev. B, 54, 11169 (1996).[36] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett., 77, 3865 (1996).[37] D.A. Broido, M. Malorny, G. Birner, N. Mingo, D.A. Stewart, Appl. Phys. Lett. 91, 231922 (2007).[38] K. Esfarjani, H.T. Stokes, Phys. Rev. B, 77, 144112 (2008).[39] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, Comput. Phys. Commun., 185, 1747 (2014).[40] Z. R. Han, X. L. Yang, W. Li, T. L. Feng and X. L. Ruan, Comput. Phys. Commun., 270, 108179 (2021).[41] S. Plimpton, J. Comput. Phys., 117, 1 (1995).[42 J. Tersoff, Phys. Rev. B, 39, 5566 (1989).[43] J. Tersoff, Phys. Rev. B, 37, 6991 (1988).[44] A. Agrawal, R. Mirzaeifar, Comput. Mater. Sci., 188, 110204 (2021).[45] B. Jelinek, S. Groh, M. F. Horstemeyer, J. Houze, S. G. Kim, G. J. Wagner, A. Moitra, and M. I. Baskes, Phys. Rev. B, 85, 245102 (2012).[46] A. Majumdar and P. Reddy, Appl. Phys. Lett. 84, 4768 (2004).[47] P. Singh, M. Seong and S. Sinha, Appl. Phys. Lett. 102, 181906 (2013).[48] S. Sadasivam, N. Ye, J. P. Feser, et al., Phys. Rev. B 95, 085310 (2017). [49] Y. X. Xu, H. Z. Fan, Z. G. Li, Y. G. Zhou, Int. J. Heat Mass Transf., 201, 123628 (2023).[50] Z. Li, S. Y. Xiong, C. Sievers, Y. Hu, Z. Y. Fan, N. Wei, H. Bao, S. D. Chen, D. Donadio, T. Ala-Nissila, J. Chem. Phys., 151, 234105 (2019).[51] J. Che, T. Cagin, W. Deng, and W. A. Goddard, J. Chem. Phys., 113, 6888 (2000).[52] S. G. Volz and G. Chen, Phys. Rev. B, 61, 2651 (2000).[53] H. K. Dong, Z. Y. Fan, L. B. Shi, A. Harju and T. Ala-Nissila, Phys. Rev. B, 97, 094305 (2018).[54] H. K. Dong, S. Y. Xiong, Z. Y. Fan, P. Qian, Y. J. Su, and T. Ala-Nissila, Phys. Rev. B, 103, 035417 (2021).[55] P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B, 65, 144306 (2002).[56] A. Dhar, Adv. Phys., 57, 457 (2008). [57] P. Jund and R. Jullien, Phys. Rev. B, 59, 13707 (1999).[58] Z. Y. Fan, L. F. C. Pereira, H. Q. Wang, J. C. Zheng, D. Donadio and A. Harju, Phys. Rev. B, 92, 094301 (2015).[59] B. Ni, T. Watanabe and S. R. Phillpot. J. Phys.: Condens. Mat., 21, 084219 (2009).[60] E. A. Burgemeister and C. A. J. Ammerlaan, Phys. Rev. B, 21, 2499 (1980).[61] W. S. Yu, L. P. Wu, S. P. Shen. Mater. Sci. Eng. A, 695, 239 (2017).[62] R. Q. Hood, P. R. C. Kent, R. J. Needs, and P. R. Briddon, Phys. Rev. Lett., 91, 076403 (2003).[63] K. Sankaran, S. Clima, M. Mees, and G. Pourtois, ECS J. Solid State Sci. Technol., 4, N3127 (2015).[64] X. B. Li, R. G. Yang, J. Phys. Condens. Mat., 24, 155302 (2012).[65] D. Wu, X. H. Cao, P. Z. Jia, Y. J. Zeng, Y. X. Feng, L. M. Tang, W. X. Zhou, K. Q. Chen, Sci. China Phys. Mech., 63, 276811 (2020).[66] P. Z. Jia, Y. J. Zeng, D. Wu, H. Pan, X. H. Cao, W. X. Zhou, Z. X. Xie, J. X. Zhang and K. Q. Chen, J. Phys. Condens. Mat., 32, 055302 (2019).[67] J. Chen, X. F. Xu, J. Zhou, and B. W. Li, Rev. Mod. Phys., 94, 025002 (2022).[68] E. T. Swartz, R. O. Pohl, Rev. Mod. Phys., 61, 605 (1989).[69] M. Kazana, Appl. Phys. Lett., 95, 141904 (2009).[70] Z. X. Lu, N. P. Smith, M. P. Prange, R. A. Bunker, J. L. Orrell, A. M. Chaka, Phys. Rev. Mater., 5, 086002 (2021).[71] K. P. Wu, L. Zhang, D. B. Wang, G. C. Chen, F. Z. Li, P. Z. Zhang, L. W. Sang, M. Y. Liao, Diam. Relat. Mater., 129, 109390 (2022).[72] A. Togo, I. Tanaka, Scr. Mater., 108, 1 (2015).[73] J. L. Warren, J. L. Yarnell, G. Dolling, R. A. Cowley, Phys. Rev., 158, 805(1967). [74] G. Nilsson, S. Rolandson, Lattice dynamics of copper at 80 K, Phys. Rev. B, 7, 2393 (1973).[75] T. S. English, J. C. Duda, J. L. Smoyer, D. A. Jordan, P. M. Norris, and L. V. Zhigilei, Phys. Rev. B, 85, 035438 (2012).[76] T. Liang, M. Zhou, P. Zhang, P. Yuan, D. G. Yang, Int. J. Heat Mass Transfer, 151, 119395 (2020).[77] H. Bao, J. Chen, X. K. Gu, B. Y. Cao, ES Energy Environ., 1, 16 (2018).[78] G. C. Loh, E. H. T. Teo and B. K. Tay, Diamond Relat. Mater., 23, 88 (2012).[79] B. w. Li, J. H. Lan, and L. Wang, Phys. Rev. Lett., 95, 104302 (2005).[80] Y. X. Xu, L. N. Yang, Y. G. Zhou, Phys. Chem. Chem. Phys., 24, 24503 (2022).[81] C. Du, R Ye, X. Cai, X. Duan, H. Liu, Y. Zhang, G. Qiu, M. Mi, J. Semicond., 44, 121801 (2023). [82] A. Kobayashi, H. Tomiyama, Y. Ohno, Y. Shimizu, Y. Nagai, N. Shigekawa, J. Liang, Functional Diamond, 2, 142 (2022).2oleObject1.binimage3.wmf()ktoleObject2.binimage4.jpegimage5.wmf011(1)()zzLLlkk=+oleObject3.binimage6.wmf0koleObject4.binimage7.wmf()zLkoleObject5.binoleObject6.binimage8.wmf()||zQLTk=ÑoleObject7.binimage9.wmfQoleObject8.binimage10.wmf1dEQSdt=oleObject9.binoleObject10.binimage11.jpegimage12.jpegimage13.wmfintintE()/SCTotalCdiamondCuCuENENE=--oleObject11.binimage14.wmfTotalEoleObject12.binimage15.wmfCdiamondEoleObject13.binimage16.wmfCuEoleObject14.binoleObject15.binoleObject16.binimage17.wmfNCuoleObject17.binimage18.wmfNColeObject18.binimage19.wmfintITC/(S)JT=D·oleObject19.binimage20.jpgimage21.jpegimage22.wmfmod2,2,mod112modmod1,1,2,2,mod1mod1ejnjjejeknejnkkjjekejDOSVDOSVDOSVa--®----×=×+×åååoleObject20.binimage23.wmf112212211224()VVVVrrarr®=+oleObject21.binimage24.wmf12a®oleObject22.binimage25.jpegimage26.wmfPDOS()()iteVACFtdtww+¥--¥=òoleObject23.binimage27.wmf()11(0)()NiiiNVVVACtFt-==åoleObject24.binimage28.wmfi(0)VoleObject25.binimage29.wmf()iVtoleObject26.binimage30.wmfCCu0CCu00PDOS()PDOS()S=PDOS()PDOS()dddwwwwwww¥¥¥òòòoleObject27.binimage31.wmf224(PDOS())1PPR()=NPDOS()nnnnwwwååoleObject28.binimage32.jpegimage33.jpgimage34.jpegimage1.jpegimage2.wmf21z0()()(0)()BzkTVJJtdttkt-=ò