Thermal Transport Across Metal Silicide-Silicon Interfaces: First-Principles Calculations and Green’s Function Transport Simulations

Thermal Transport Across Metal Silicide-Silicon Interfaces: First-Principles Calculations and Green’s Function Transport Simulations


Heat transfer across metal-semiconductor interfaces involves multiple fundamental transport mechanisms such as elastic and inelastic phonon scattering, and electron-phonon coupling within the metal and across the interface. The relative contributions of these different transport mechanisms to interface conductance remains unclear in the current literature. In this work, we use a combination of first-principles calculations under the density functional theory framework and heat transport simulations using the atomistic Green’s function (AGF) method to quantitatively predict the contribution of the different scattering mechanisms to the thermal interface conductance of epitaxial CoSi-Si interfaces. An important development in the present work is the direct computation of interfacial bonding from density functional perturbation theory (DFPT) and hence the avoidance of commonly used ‘mixing rules’ to obtain the cross-interface force constants from bulk material force constants. Another important algorithmic development is the integration of the recursive Green’s function (RGF) method with Büttiker probe scattering that enables computationally efficient simulations of inelastic phonon scattering and its contribution to the thermal interface conductance. First-principles calculations of electron-phonon coupling reveal that cross-interface energy transfer between metal electrons and atomic vibrations in the semiconductor is mediated by delocalized acoustic phonon modes that extend on both sides of the interface, and phonon modes that are localized inside the semiconductor region of the interface exhibit negligible coupling with electrons in the metal. We also provide a direct comparison between simulation predictions and experimental measurements of thermal interface conductance of epitaxial CoSi-Si interfaces using the time-domain thermoreflectance technique. Importantly, the experimental results, performed across a wide temperature range, only agree well with predictions that include all transport processes: elastic and inelastic phonon scattering, electron-phonon coupling in the metal, and electron-phonon coupling across the interface.

I Introduction

Interfaces between heterogeneous materials provide a plethora of possibilities for the design of devices with engineered electronic and optical properties. This work concerns the study of heat transport across metal-semiconductor heterojunctions that form a technologically important class of interfaces used in electronic devices. The understanding of charge and heat transport through metal contacts to semiconductor channels is critical to ensure reliable operation of field effect transistors that form the basic building block of high-power electronic devices. Understanding of thermal transport through metal-semiconductor interfaces is also important in the design of modern memory storage devices such as heat assisted magnetic recording Kryder et al. (2008) and phase change memory Reifenberg et al. (2008). Apart from their technological relevance, metal-semiconductor interfaces also provide a material system in which various physical mechanisms of heat transport such as elastic interfacial phonon scattering, inelastic phonon scattering, and electron-phonon coupling co-exist. In this work, we provide a rigorous modeling framework to understand the contribution of various interfacial scattering mechanisms to thermal transport across cobalt silicide (CoSi) - silicon interfaces that are extensively used in microelectronic devices Murarka (1995).

Elastic scattering of phonons at an interface is the most widely studied framework to understand and predict thermal interface conductance at heterojunctions. Under the elastic transport framework, a phonon of energy incident from one side of an interface is either transmitted across the interface or reflected back into the same material. For elastic interfacial transport, the primary quantity of interest is the phonon transmission function that represents the probability that a phonon incident from one side of the interface transmits to the other side. Anharmonic or three-phonon scattering processes typically become important at room temperature and above, in which a phonon of energy incident on the interface could transmit or reflect multiple phonons with appropriate energies to ensure energy conservation. This mechanism has been postulated to be important in acoustically mismatched interfaces such as Pb-diamond Kosevich (1995); Hohensee et al. (2015).

Electrons are the primary energy carriers in metals while phonons are dominant in intrinsic semiconductors. Hence, electron-phonon coupling can be another important energy transfer mechanism that affects thermal interface conductance in metal-semiconductor interfaces. Electron-phonon coupling within the metal provides an additional resistance to heat transfer Majumdar and Reddy (2004). However, electron-phonon coupling across an interface, i.e., coupling between metal electrons and semiconductor phonons, provides a parallel heat flow path in addition to phonon-phonon heat transfer across the interface. Time domain thermoreflectance (TDTR) experiments in the literature Hopkins et al. (2009); Guo et al. (2012) suggest that direct electron-phonon coupling can contribute significantly to heat transport across metal-semiconductor interfaces, and models Huberman and Overhauser (1994); Sergeev (1999); Mahan (2009); Ren and Zhu (2013); Lombard et al. (2014); Lu et al. (2016a) have been developed to quantify its contribution. The different mechanisms of heat transport at a metal-semiconductor interface are summarized in Figure 1.

Simplified empirical models are commonly used to interpret experimental thermal conductance data for metal-semiconductor interfaces. Elastic interfacial phonon scattering is commonly modeled using the acoustic Little (1959) and diffuse Swartz and Pohl (1989) mismatch models (AMM, DMM) which are heuristic approaches applicable for smooth and rough interfaces respectively. Also, simplifying assumptions such as the Debye approximation to phonon dispersion can compromise the quantitative accuracy of such models. Even atomistic simulation approaches for elastic interfacial thermal transport such as the atomistic Green’s function (AGF) method often involve empirical force constant models that can produce significant discrepancies when compared to calculations that employ harmonic force constants obtained from ab initio approaches Tian et al. (2012). The contribution of inelastic phonon scattering to thermal interface conductance has also been modeled in a simplified manner using heuristic extensions to the elastic diffuse mismatch model Hopkins (2009). The strength of electron-phonon coupling is typically modeled using idealized approximations such as bulk metal deformation potentials Huberman and Overhauser (1994); Sergeev (1998), and such approximations are expected to be inaccurate for the direct coupling of metal electrons with joint or interface phonon modes. Little work exists on rigorous first-principles determination of the strength of coupling between electrons and joint interface phonon modes at a metal-semiconductor interface.

Apart from the complexity of various thermal transport mechanisms described above, the uncertainty in interfacial atomic structure has historically made direct comparisons between simulations and experiments difficult. Much of the existing experimental data Stoner and Maris (1993); Stevens et al. (2005); Guo et al. (2012) on thermal conductance of metal-semiconductor interfaces involves materials with mismatched lattice constants, for which the interface atomic structure is likely to be at least partially amorphous. Experimental studies that simultaneously characterize interfacial atomic structure along with interface conductance are scarce Costescu et al. (2003); Wilson et al. (2015). However, predictive atomistic transport simulations that involve first-principles approaches are typically limited to crystalline epitaxial interfaces because of the associated computational tractability. This disparity between simulations and experimental studies often makes quantitative comparisons challenging. To overcome this difficulty, we choose to work with CoSi (metal) - Si (semiconductor) interfaces in the present work. Both CoSi and Si have FCC lattice structures with similar lattice constants of 5.36 Å and 5.43 Å respectively. Measurements of thermal interface conductance on CoSi (111)/ Si (111) interface using the TDTR technique have been reported in our recent work Ye et al. (2016), and the interface has been verified to be epitaxial and smooth using TEM imaging (see ref. Ye et al., 2016 for a TEM image of the interface). We use the same experimental data to compare with the present simulation predictions on a lattice-matched CoSi (111)/Si (111) interface; the interfacial atomic configuration was also chosen to match with the atomic configuration of samples used in the experiment (see Section II.2 for details of the various interfacial atomic configurations). The close correspondence between the atomic structures used in the present work and the experimental data reported in ref. Ye et al., 2016 enables a direct comparison between simulations and experiments.

Although the primary focus of the present work is the study of thermal transport across metal-semiconductor interfaces, the methods developed and reported here are also expected to be useful for a broad class of problems that use the non-equilibrium Green’s function (NEGF) method for atomistic transport simulations. From a methodology standpoint, we report a framework that combines first-principles calculations of interatomic force constants with the atomistic Green’s function method and evaluate the validity of the ‘mixing rule’ that is commonly used to approximate interfacial bonding at a heterojunction. The conventional AGF method that is suitable for elastic phonon transport is extended to include anharmonic phonon scattering using a Büttiker probe approach Miao et al. (2016). Since the Büttiker probe approach is not directly compatible with the conventional recursive Green’s function (RGF) method (see Section II.3.1), we develop a modification that enables the use of the RGF method in simulations that involve Büttiker probe scattering. The new RGF algorithm enables computationally efficient simulations of phonon-phonon scattering using the Büttiker probe approach and is expected to be applicable for a wide range of problems that require efficient representation of dephasing processes under the NEGF framework. Ab initio calculations of electron-phonon coupling are also integrated into the AGF transport simulations.

Apart from the development of new methods, the present work also provides useful insights into the physics of thermal transport across metal-semiconductor junctions. Rigorous first-principles calculations indicate that elastic phonon transport under-predicts the experimental data over a wide temperature range. Analysis of the cross-interface heat flux accumulation function provides useful insights on the microscopic mechanisms responsible for increased interface conductance due to anharmonic phonon scattering in the bulk materials forming the interface. First-principles calculations of electron-phonon coupling on an interface supercell along with a detailed analysis of the contribution from different kinds of phonon modes to the Eliashberg function reveal that delocalized phonon modes mediate cross-interface energy transfer between metal electrons and the semiconductor lattice. We also obtain an effective length scale of electron-phonon interaction in the semiconductor by comparing simulation predictions with experimental data and evaluate the accuracy of prior approximations to the length scale of joint or interface phonon modes.

Figure 1: Schematic of various mechanisms involved in heat transfer between the dominant energy carriers, i.e., electrons in the metal and phonons in the semiconductor. Phonon-phonon energy transfer across the interface could involve elastic and inelastic interfacial scattering processes. Electron-phonon coupling could involve coupling between electrons in the metal with phonons in the metal and with phonons in the semiconductor.

Ii Methods

ii.1 First-principles Calculations

All first-principles calculations in this paper were performed under the framework of density functional theory (DFT) using a planewave basis set as implemented in the Quantum Espresso suite of codes Giannozzi et al. (2009). Rappe-Rabe-Kaxiras-Joannopoulos (RRKJ) ultrasoft pseudopotentials were used for both Co and Si atoms, and the exchange correlation energy was approximated under the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) functional form. Three sets of first-principles calculations are performed for the results reported in this paper; these involve calculations on bulk Si (6 atom non-primitive unit cell along [111] direction), bulk strained CoSi (9 atom unit cell along [111] direction) where a tensile strain is applied along the (111) plane, and a Si (111)-CoSi (111) interface supercell. The relaxed lattice constants of bulk Si and bulk CoSi are 5.44 Å and 5.36 Å respectively. For all simulations considered in this paper, a tensile strain of 1.5% is applied on CoSi along the (111) plane to match the lattice constants of Si and CoSi. Table 1 shows the cutoff energies and k-point grids used for DFT calculations on all three systems. Structural relaxation is carried out to reduce the Hellmann-Feynman forces on every atom below 10 eV/Å. A full stress relaxation is carried out for bulk Si while the stresses on bulk strained CoSi and the interface supercell are relaxed only along the transport direction. CoSi is stretched along the in-plane direction to match its lattice with Si, and hence the in-plane stresses are not relaxed.

Phonons are analyzed using density functional perturbation theory (DFPT) where the dynamical matrices are obtained on a q-point grid for bulk Si (6 atom unit cell along the [111] direction), bulk strained CoSi (9 atom unit cell along the [111] direction) and a q-point grid for the interface supercell. The real-space inter-atomic force constants (IFCs) needed for Green’s function transport calculations are obtained from a Fourier transform of the dynamical matrices.

Parameter Bulk Si Bulk strained CoSi Interface supercell
Kinetic energy cutoff (eV) 680 820 820
Charge density cutoff (eV) 6800 8200 8200
Electron k-point grid
Phonon q-point grid
Table 1: Parameters used for DFT, DFPT calculations on bulk Si, bulk strained CoSi, and the Si-CoSi interface supercell.

ii.2 Coherent Phonon Transport using the Atomistic Green’s Function Method

The terms ‘coherent’ and ‘ballistic’ phonon transport are used interchangeably in the present manuscript and refer to simulations performed under the conventional AGF framework. These simulations do not model phonon dephasing (cf., coherent) and inelastic phonon scattering (cf., ballistic). However, elastic interfacial scattering, i.e., reflection and transmission of phonon waves at the interface is included in this framework. The next section describes modifications to the conventional AGF approach to model phonon dephasing and inelastic phonon scattering. The details of the AGF method are available in prior reports Zhang et al. (2007); Sadasivam et al. (2014), and a brief description is provided in this section. The basic conceptual framework for the AGF method involves a ‘device’ region that is connected to semi-infinite ‘leads’. The device Green’s function is given by:


where is the force constant matrix corresponding to the device region, and , are the contact self-energies. The contact self-energies are obtained from the surface Green’s functions , as follows:


where , represent the force constant matrices for interaction between atoms in the device region and the semi-infinite contacts. and are the surface Green’s functions of the contacts that are obtained using the Sancho-Rubio method Sancho et al. (1985); Guinea et al. (1983). The phonon transmission function across the device is obtained from the Caroli formula:


where denotes the imaginary part of the contact self-energies and physically represents the phonon ‘escape rate’ Sadasivam et al. (2014) from the device into the respective contacts. In all the above expressions, the dependence of the Green’s function, contact self-energy, and the transmission function on the transverse phonon wavevector (for structures with periodicity in the transverse or in-plane direction) and frequency is implicitly assumed. After obtaining the transmission function, the thermal interface conductance can be obtained using the Landauer formula:


ii.3 Inelastic Scattering Using Büttiker Probe Approach

The AGF formulation presented in the previous section is applicable only for elastic phonon transport, i.e., anharmonic scattering mechanisms such as Umklapp scattering are not considered. The formulation for extension of AGF to include anharmonic phonon scattering has been developed in ref. Mingo, 2006; however, the approach is computationally expensive, and we are not aware of its application to study phonon transport through realistic three-dimensional crystals. The authors recently proposed a phenomenological Büttiker probe approach to model anharmonic phonon scattering within the AGF method Miao et al. (2016). The approach is an extension to phonons of the widely used Büttiker probe method to model inelastic electron scattering processes in the NEGF framework Venugopal et al. (2003); Maassen et al. (2009); Afzalian (2011). Although the method is heuristic, it provides a computationally efficient alternative to the self-consistent Born approximation (SCBA) Luisier (2012) that is not phenomenological but is computationally intensive in both memory and time. The essence of the Büttiker probe approach involves attaching fictitious contact probes to every atom in the device, and the temperatures of these fictitious contacts are then iteratively solved to ensure energy conservation in the device region. The Büttiker probes contribute an additional self-energy to the device Green’s function (in addition to the self-energies due to the real contacts):


In the present formulation, the Büttiker probe self-energy is assumed to be a diagonal matrix whose diagonal elements are of the form:


where denotes the Büttiker probe self-energy at atom and vibrational direction (, , ). Similar to the matrices , , we also define that represents the imaginary part of the Büttiker probe self-energy. denotes the frequency dependent scattering time due to Umklapp scattering and is assumed to be of the form for both Si and CoSi. Quadratic frequency dependence of the Umklapp scattering rate has been used in prior studies involving the BTE Singh et al. (2011) and Landauer approach Mingo (2003); Jeong et al. (2012). The parameter is chosen by fitting (see Supplemental Material) to the experimental thermal conductivity of bulk Si (at different temperatures) and bulk CoSi (at room temperature). Due to lack of experimental data on the temperature dependence of the lattice thermal conductivity, the scattering parameter is assumed to be independent of temperature for CoSi. This assumption is acceptable because the thermal conductivity of CoSi is dominated by electrons, and the interface conductance is found to exhibit a weak dependence on the lattice thermal conductivity of CoSi (see Supplemental Material).

Recursive Green’s Function Method for Efficient Solution of Büttiker Probe Temperatures

Büttiker probes offer a heuristic but efficient method to implement scattering in NEGF simulations. However, the popular recursive Green’s function (RGF) method Anantram et al. (2008) that avoids full inversion in the calculation of the retarded Green’s function and the lesser Green’s function is not compatible with Büttiker probes. This incompatibility can be understood from the following equation used to enforce heat current conservation in each Büttiker probe .


where the summation in the variable runs over all other Büttiker probes () and the contacts. The computation of the transmission function matrix between every pair of Büttiker probes requires calculation of the full Green’s function matrix . The equation for charge current conservation in electronic transport is similar to the foregoing equation for heat current conservation Venugopal et al. (2003); Afzalian (2011). Hence, prior implementations Venugopal et al. (2003); Afzalian (2011); Miao et al. (2016) of the Büttiker probe formalism have employed direct matrix inversion instead of the RGF method to calculate the full device Green’s function matrix.

Eq. (7) enforces the condition that the total integrated energy flux in each Büttiker probe is zero, i.e., inelastic scattering between different energy levels is allowed. Alternative implementations of the Büttiker probe approach invoke energy flux conservation at each phonon frequency instead of the total integrated energy flux over all phonon frequencies Venugopal et al. (2003); Maassen et al. (2009). Such an approach is suitable only for elastic dephasing and is hence not adopted in this work.

Eq. (7) can be cast in a slightly different form as:


where and denotes the spectral function. In the above equation, the matrices and are block-diagonal. Hence, only the block-diagonals of and need to be computed, and this can be done using the RGF algorithm. However, the computation of and require knowledge of the Büttiker probe temperatures. Hence the use of RGF with Büttiker probes requires that and are recalculated during every Newton iteration of the solution for Büttiker probe temperatures. Although this step is not needed in a conventional Büttiker probe implementation with full matrix inversion, the computational advantage of RGF over full inversion far outweighs the computational expense of repeated RGF calculations in every Newton iteration. Also, the memory required to store and invert the full Green’s function matrix can become prohibitively large with increasing device length.

Anantram et al. Anantram et al. (2008) provide a detailed discussion of the RGF methodology for computation of the block diagonals of and . Here, we provide only the modifications needed to combine RGF with the Büttiker probe formalism. The RGF algorithm for computation of the block-diagonal elements of the retarded Green’s function remains unchanged. However, the RGF algorithm for computation of requires modification to also calculate the derivative of the diagonal elements of with respect to the Büttiker probe temperatures. We also assume that all Büttiker probes within a RGF ‘block’ have the same temperature, i.e., the number of Büttiker probe temperatures that need to be solved is equal to the number of blocks in the device region.

The equation for left-connected is given by Anantram et al. (2008):


where , . Our terminology follows that of ref. Anantram et al., 2008 where denotes the left-connected retarded Green’s function, and denotes the lesser self-energy. The derivative of with respect to the Büttiker probe temperature is given by:


The equation for is given by:


The derivative of with respect to Büttiker probe temperatures can be computed using the derivatives of the left connected Green’s function computed in Eq. (10).


Overall, the RGF algorithm for needs to be modified to compute the derivatives of with respect to the Büttiker probe temperatures. The new RGF algorithm’s commutation involves the following steps:

  1. , ,

  2. For and , compute Eqs. (9) & (10)

  3. , for

  4. For and , compute Eqs. (11) & (12).

The algorithm proceeds as follows.

  1. Start with an initial guess for the Büttiker probe temperatures.

  2. For each phonon frequency, compute , , and using the RGF algorithm described above.

  3. Compute energy current densities in each Büttiker probe using Eq. (8).

  4. Compute the Jacobian matrix whose element is given by the following equation:


    where is the Kronecker Delta function.

  5. Update the temperature of Büttiker probes using the Newton equation:

  6. If , go back to Step 1 with the new guess for Büttiker probe temperatures as .

An alternative to the Newton-Raphson method is the secant method in which the exact Jacobian needs to be computed only in the first iteration. For further iterations, the Jacobian could be updated using the Broyden’s update formula Broyden (1965), and this method was also found to give satisfactory convergence. With the secant method, the derivative of the lesser Green’s function with respect to Büttiker probe temperatures need only be computed in the first iteration. For the remaining iterations, the traditional RGF function is sufficient. For large device lengths, the computation and storage of the full Jacobian matrix can become prohibitively expensive; an alternative approach involves approximation of the Jacobian by a sparse block diagonal matrix. Such an approximate Jacobian was also found to lead to convergence, however with an increased number of iterations compared to the exact Jacobian. The approximate Jacobian provides a memory-time tradeoff as storage of the sparse block diagonal matrix involves lesser memory but the computational time increases relative to the Newton Raphson scheme with exact Jacobian. All the results presented in this paper involve the secant method in which the exact Jacobian is computed in the first iteration and the Broyden’s update formula is used for further iterations.

Figure 2 shows a comparison of the computational times for AGF simulations of bulk silicon with Büttiker probe scattering using direct inversion and the RGF algorithm described above. The computational times were obtained using MATLAB, and both direct inversion and the RGF methods were parallelized over transverse wavevectors. The computational time for full matrix inversion increases rapidly with device length (matrix inversion scales as where , denote the number of atoms per slab and the number of slabs in the transport direction respectively) while that for the RGF algorithm proposed above shows a more gradual scaling with device length (RGF scales as ).

Apart from the computational time improvement, another important advantage of the RGF method over full inversion is the reduced memory needed to store and invert full Green’s function matrices. The device sizes considered in the present work (see Section IV) are computationally intractable with direct matrix inversion. Hence, the extension of the RGF algorithm to Büttiker probes is expected to be critical for application of the Büttiker probe method to realistic device sizes. Also, the results in Figure 2 confirm that the computational expense of repeated AGF calculations of for each RGF iteration is far less than the computational expense for a single calculation of the full retarded Green’s function of the device through direct inversion.

Figure 2: Comparison of computational times for different device lengths obtained using direct inversion and the RGF algorithm.

ii.4 Fourier Diffusion of Electrons Coupled with Phonons

Electrons are the primary heat carriers in metals, and they transfer energy to phonons near the interface between metal and semiconductor. Intrinsic Si is the semiconductor of interest in this paper, and hence the contribution of electrons in Si to thermal transport is neglected. We neglect any cross-interface electron tranport through the CoSi-Si interface and consider diffusive transport of electrons in CoSi. Electrons are included in the AGF simulation within the framework of a two-temperature model that is commonly used to interpret ultra-fast laser experimental data and is also used to model energy transfer between electron and phonon subsystems within the Eliashberg function framework. The primary assumption involved in the definition of a local electron and phonon temperature is the existence of electron-electron and phonon-phonon collisions that enable local equilibrium separately within electron and phonon subsystems. The electron and phonon subsystems exchange energy through electron-phonon coupling that is expressed in terms of the Eliashberg function.The steady state Fourier diffusion equation for electrons in the metal is given by:


where denotes the volumetric heat source term due to coupling between electrons and phonons and is expressed in terms of the Eliashberg function as derived by Allen in ref. Allen, 1987:


In the foregoing equation , denote the local electron and lattice temperatures respectively at location , and is the electronic density of states at the Fermi energy. The above term can be included as a source term in the Büttiker probe at location , i.e., the energy current conservation equation for the Büttiker probe in the metal is given by:


where is the length that the Büttiker probe occupies along the transport direction. In a phonon-only simulation, the total energy current in each Büttiker probe is set to zero to ensure energy flux conservation (see Eq. (8)), i.e., Büttiker probes redistribute the energy, but with no net transfer of energy between electrons and phonons via the fictitious Büttiker contacts. When electrons are included in the transport calculation, the total energy current in each Büttiker probe is set to the electron-phonon energy exchange given by Eq. (16).

The local lattice temperature in Eq. (16) is obtained by equating the local phonon energy density to the product of a local Bose-Einstein distribution at temperature and the local phonon density of states:


The above equation makes use of the following expressions for the local phonon number density and the local phonon DOS in terms of the lesser Green’s function and the spectral function respectively:


Eqs. (15), (17), (18) constitute a set of coupled non-linear equations that are solved iteratively to obtain the electron temperature, the Büttiker probe temperature, and the local device temperatures. Similar to the methodology for Büttiker probe temperatures in Section II.3.1, the Newton-Raphson method is used for the solution of the above equation and details of the algorithm are provided in the Supplemental Material.

Iii Coherent Phonon Transport

This section contains results for the phonon transmission function and thermal interface conductance of Si-CoSi interface from ballistic phonon transport calculations (i.e., Büttiker probe scattering turned off). Different interfacial atomic configurations are possible for the Si-CoSi interface depending on the coordination number of the Co atom closest to the interface (possible values of 5, 7, 8) and the relative crystal orientation between the (111) surfaces of Si and CoSi. The ‘A’ type orientation occurs when the Si-CoSi stacking is continuous while the ‘B’ orientation occurs when the CoSi crystal is rotated by 180 about the [111] direction. Previous first-principles calculations in the literature Stadler et al. (1999); Wardle et al. (2005) indicate that the 8A and 8B configurations have the lowest interfacial energies and are hence the most probable interfacial atomic structures. Both configurations are considered for the present phonon transport calculations using AGF.

While bulk IFCs are quite commonly obtained from first-principles calculations, little work exists on the use of rigorous DFPT calculations to determine the force constants between atoms belonging to different materials across a heterogeneous interface. Force constants (in AGF) and inter-atomic potentials (in molecular dynamics) between atoms belonging to different bulk materials are commonly represented using simplifying approximations without rigorous calculations of the actual strength of interfacial bonding Huang et al. (2011); Ong and Pop (2010). In the present work, both bulk and cross-interface IFCs needed for AGF transport simulations are determined entirely from DFPT calculations.

Figure 3a,b shows supercells of the 8A and 8B interfacial atomic configurations respectively. Each supercell contains two interfaces because of periodic boundary conditions. Although DFPT calculations are performed on a finite interface supercell, transport simulations are performed on a single Si-CoSi interface formed by semi-infinite Si and CoSi crystals. The red dotted boxes enclose atoms around the interface for which the force constants are obtained from the interface supercell DFPT calculation. In the atomic structure considered for transport calculations, the IFCs for atoms outside the red dotted box are assumed to equal the bulk IFCs of Si (left of the box) and CoSi (right of the box). Results that illustrate the convergence of cross-interface force constants with respect to the size of interface supercell are provided in the Supplemental Material.

Enforcement of the acoustic sum rules is an important consideration when IFCs obtained from DFPT calculations are used in thermal transport simulations. Acoustic sum rules constitute a set of translational invariance conditions on the IFCs to ensure that long wavelength acoustic modes of a crystal have zero vibrational frequency:


where , denote atom indices while , denote the directions. The spatial range of inter-atomic interactions is artificially truncated in DFPT by the finite q-point grid used in the calculations. Although the neglected long-range interactions are typically small, this procedure results in small violations of the translational invariance conditions. Hence, the raw force constants obtained from DFPT do not satisfy the acoustic sum rules exactly and need to be enforced as a post-processing step on the IFCs Mingo et al. (2008). Common DFT codes such as Quantum Espresso automatically enforce translational invariance for the IFCs of the crystal on which DFPT calculations are performed. However in the present calculations, the IFCs obtained for bulk Si and bulk CoSi are combined with that obtained for the interface supercell. Hence the IFCs for Si and CoSi unit cells nearest to the interface will require modifications to ensure that the acoustic sum rules are satisfied. In the present work, the diagonal blocks of the force constant matrix are modified to ensure that Eq. (20) is satisfied.

Figure 3: Atomic structures of Si-CoSi interface supercells used in DFPT calculations (two unit cells shown along the in-plane direction for clarity). The red dotted boxes indicate the region around the interface for which IFCs are extracted from the interface supercell calculation. a) 8A configuration. b) 8B configuration.

Cross-interface force constants between heterogeneous materials are commonly approximated using simplifying assumptions due to the computational complexity of performing direct DFPT calculations on an interface supercell with a large number of atoms. If the materials on both sides of the interface have the same lattice structure such as Si-Ge interfaces, a common approximation is to assume the same force constants for both materials with the assumption that interfacial scattering is primarily affected by the change in atomic mass across the interface Tian et al. (2012). Other approximations include the use of empirical corrections to obtain the cross-interface force constants from the bulk force constants Huang et al. (2011). Another common approximation involves the use of mixing rules to obtain the strength of cross-interface interactions from an average of the bulk material parameters Luckyanova et al. (2012).

In order to evaluate the validity of a simple averaging approximation for the cross-interface force constants, Figure 4a compares the phonon transmission function at normal incidence () when the cross interface IFCs are assumed to be a simple arithmetic average of the bulk IFCs and when the cross-interface IFCs are obtained from DFPT on the interface supercell shown in Figure 3a. In the averaging approach, the cross-interface Co-Si and Si-Si IFCs are obtained by averaging the interactions in bulk Si and bulk CoSi. The averaging approximation is found to over-estimate the transmission function for most of the frequency range except at very low frequencies or long wavelengths (see inset in Figure 4a) where the predictions from both the average and DFPT IFCs converge to the acoustic mismatch limit. Long wavelength phonons are insensitive to the local details of interfacial bonding, and hence the transmission function at low phonon frequencies is expected to be insensitive to the exact interfacial force constants. However, rigorous predictions of cross-interface force constants are necessary for accurate prediction of transmission at higher frequencies that are expected to dominate phonon transport at room temperature and beyond. The thermal interface conductance computed directly from Eq. (4) includes contributions from the ballistic contact conductances at the contact-device interfaces in addition to the conductance of the Si-CoSi interface in the middle of the device region. To obtain the conductance of the Si-CoSi interface alone, we use the following expression to subtract the ballistic contact resistances from the total resistance obtained from Eq. (4) Tian et al. (2012):


where is the interface conductance computed from Eq. (4), and is the thermal interface conductance of a single Si-CoSi interface and plotted in Figure 4b. and are the ballistic conductances of homogeneous Si and CoSi slabs respectively.

The use of a simple arithmetic average for cross-interface IFCs over-estimates the thermal interface conductance (see Figure 4b) by almost 70% at room temperature, and the errors increase at higher temperatures. Hence, the prediction of phonon thermal interface conductance for temperatures beyond a few tens of K requires the rigorous prediction of interfacial bonding strength, and simple averaging approximations are not expected to be quantitatively accurate. Similar conclusions on the over-estimation of interface conductance due to simple approximations that neglect local changes in the force field near a heterogeneous interface were found in ref. Gu et al., 2015.

Figure 4: Results from ballistic phonon transport calculations for a Si-CoSi interface. a) Phonon transmission function at normal incidence computed using average and DFPT force constants for the 8A interface. The inset shows the same graph for small phonon frequencies or long wavelengths. b) Thermal interface conductance for 8A and 8B Si-CoSi interfaces.

Iv Effect of Anharmonic Scattering on Thermal Interface Conductance

In this section, the effect of anharmonic phonon scattering in Si and CoSi on the thermal interface conductance is presented. This section contrasts with results in the previous section for which phonon transport in Si and CoSi were assumed to be ballistic. The Büttiker probe scattering rates for both Si and CoSi were assumed to be of the form , and the parameter was fitted to obtain the bulk thermal conductivity of Si and CoSi (see Supplemental Material). Since Si has a relatively high phonon thermal conductivity compared to CoSi, this circumstance corresponds to a low scattering rate on the Si side of the interface and a high scattering rate in CoSi. To understand the effect of bulk scattering rates on the interface conductance, we also performed simulations in which the scattering rate in CoSi is reduced by a factor of 100 while the Si scattering rate is maintained the same (Case A) and the scattering rate in Si is increased by a factor of 100 while that in CoSi is maintained the same (Case C). Case B corresponds to the nominal scattering rate in both Si and CoSi, i.e., the scattering rate that reproduces the bulk thermal conductivity of Si and CoSi. The present simulations assume that all the Büttiker probes on the Si and CoSi sides of the interface have scattering rates of bulk Si and bulk CoSi respectively. However, the local anharmonicity near the Si-CoSi interface is likely to differ from the bulk anharmonicities of Si and CoSi. Future work on determining the change in Umklapp scattering rates near a heterogeneous interface is needed to improve the present model.

Figures 5a,b,c show the local device temperature profile corresponding to all three cases and the associated thermal interface conductance. A temperature difference of 10 K is applied across the leads in all cases. The conductance is enhanced with increase in the bulk scattering rates of the materials comprising the interface. As expected from conventional scattering theory, the bulk material conductances however decrease with increased scattering rates (observe the progressive rise in temperature drops within Si and CoSi in Figures 5a,b,c). The foregoing results suggest that the inclusion of inelastic phonon scattering in the AGF simulations produces contrasting effects on the interface and bulk material conductances.

To elucidate the microscopic mechanisms responsible for the enhancement in interface conductance with inelastic scattering, Figures 5d,e,f show the spectral variation of heat flux from the Si and CoSi contacts. For Case A with low scattering on both sides of the interface, the spectral heat fluxes from the two contacts follow each other, suggesting that ‘vertical transport’, i.e., mixing between different energy levels is insignificant. However, higher scattering rates result in a shift of the frequencies at which the spectral heat flux is a maximum. Also, the maximum allowed phonon frequency in strained CoSi is about rad/s while that in Si is close to rad/s. Hence the phonons in Si between rad/s and rad/s do not contribute to cross-interface heat transport in a ballistic simulation (see Figure 4a). Inelastic scattering enables phonon scattering into the high energy optical modes of Si whose contribution is enhanced with a rise in the scattering rates of Si and CoSi. The elevated participation of high-energy optical phonons in Si can also be observed in Figures 5g,h,i where phonons with frequencies larger than rad/s contribute 8%, 10%, and 18% of the total energy flux in Si for cases A, B, and C respectively. Although the spectral heat flux shows spatial variation, the total energy flux integrated over all phonon frequencies is independent of position. Similar conclusions on the enhancement of thermal interface conductance due to inelastic scattering have been reported in prior work Kosevich (1995); Landry and McGaughey (2009); Hohensee et al. (2015).

Figure 5: a,b,c) Device temperature profile for Cases A, B, C where Case B corresponds to nominal scattering rates in Si and CoSi while Cases A and C correspond to artificially decreased and increased scattering rates respectively. The magenta lines correspond to linear fits of the temperature profiles on either side of the interface. d,e,f) Spectral variation of the energy flux from Si and CoSi contacts for Cases A, B, C respectively. g,h,i) Accumulation of energy flux in the Si and CoSi contacts with respect to phonon frequency for Cases A, B, C respectively.

V Effect of Electron-Phonon Coupling on Thermal Interface Conductance

v.1 First-Principles Calculations

The results from first-principles calculations of electron-phonon coupling, both in bulk strained CoSi and Si-CoSi interface supercells, are reported in this section. The phonon linewidth due to electron-phonon scattering is given by Allen (1987); Bauer et al. (1998):


where is the electron-phonon coupling matrix element for scattering of an electron with energy in band by a phonon of energy into a state with energy in band . The above expression is valid at low temperatures when electron-phonon scattering is restricted to a narrow energy window around the Fermi surface. The phonon linewidth can be used to compute the spectral Eliashberg function which quantifies the strength of electron-phonon coupling:


The spectral Eliashberg function can be used to obtain an effective volumetric electron-phonon coupling coefficient :


To understand the spatial variation of electron-phonon coupling across a semiconductor-metal interface, we also define a local Eliashberg function as Sadasivam et al. (2015):


where denotes the phonon eigenvector, the index runs over all the atoms in the unitcell, and the index represents the vibrational degrees of freedom (, , ) for each atom. The local Eliashberg function is then used to define a local volumetric electron-phonon coupling coefficient :


where the additional factor ensures that is the local volumetric coupling coefficient around atom that occupies a volume .

The Eliashberg functions of bulk strained CoSi and the interface supercell calculated from DFPT are provided in the Supplemental Material along with a discussion on convergence with respect to k-point grid and smearing parameters. Eqs. (22) and (23) are appropriate for bulk materials in which translational periodicity is assumed in the scattering matrix elements . These equations also apply for the Si-CoSi interface supercells because these supercells represent Si-CoSi superlattices with periodicities of the order of a few nm. The Eliashberg function for the supercells physically represents the coupling between electron and phonon modes of the Si-CoSi superlattice. To ensure that the electron-phonon coupling coefficients obtained from DFT/DFPT calculations on superlattices are transferrable to transport simulations of a single Si-CoSi interface, we performed calculations on a series of Si-CoSi supercells with varying Si and CoSi slab thicknesses. Figure 6a shows the spatial variation of the electron-phonon coupling coefficient for three interface supercells (SC) of the 8B configuration with different lengths of the Si and CoSi slabs forming the interface. The local coupling coefficient on the CoSi side of the interface is averaged over one Co and two Si atoms to remove atomistic fluctuations in the local coupling coefficient. The electron-phonon coupling coefficients for bulk strained CoSi and bulk intrinsic Si are also shown in Figure 6a. The electron-phonon coupling in intrinsic bulk Si is zero since the Fermi level lies in the middle of the bandgap, and the delta functions around the Fermi surface in Eq. (22) are zero. An important observation from Figure 6a is the appearance of a non-zero coupling coefficient on the Si side of the interface. Also, the magnitude of this coupling is approximately constant in Si beyond two atomic layers from the interface for all three supercells considered in Figure 6a. The convergence of the plateau on the Si side of the interface for supercells SC2 and SC3 in Figure 6a indicates that the electronic wavefunctions of CoSi are sufficiently localized within the CoSi slab and do not tunnel across the Si slabs of the superlattice. The Si-CoSi interface forms a Schottky barrier with a p-type Schottky barrier height of 0.2 eV. Hence, the local electronic DOS at the Fermi level decays rapidly in Si away from the Si-CoSi interface (see Figure 6b). However, such a two-order-of-magnitude decay in local electronic DOS does not result in a commensurate reduction in the local electron-phonon coupling coefficient shown in Figure 6a.

This unusual result can be understood by considering the relative contributions of different types of phonon modes of the Si-CoSi interface supercell to the overall Eliashberg function. The different phonon modes of the interface supercell are classified into four types based on the spatial localization of the phonon eigenvector corresponding to the mode. The present approach is analogous to the classification of interface modes in ref. Gordiz and Henry, 2016a. The interface supercell (SC 2 in Figure 6a) shown in Figure 7a is decomposed into three regions consisting of Si, CoSi, and interfacial atoms. The criteria for classification of a phonon mode is defined as follows:


where , , and denote the norm of the phonon eigenvector within the Si, CoSi and interfacial regions respectively in Figure 7a. denotes the norm of the phonon eigenvector of the entire interface supercell and is normalized to unity. The choice of spatial extent of the interfacial region and the value 0.85 in Eq. (27) are arbitrary and used only to provide a physical understanding of the mechanism of cross-interface coupling between electrons in metal and phonons in the semiconductor.

Figure 6: a) Spatial variation of the electron-phonon coupling coefficient across a Si-CoSi interface for different supercell lengths. The coupling coefficient in bulk strained CoSi and bulk Si are also shown for comparison. b) Spatial variation of the local electron DOS at the Fermi energy across the Si-CoSi structure shown in Figure 3c.

The contribution of the different types of phonon modes to the total phonon DOS and the total Eliashberg function of the interface supercell are shown in Figures 7b,c respectively. Figure 7b indicates that phonon modes in the frequency range of rad/s are localized in the Si region of the interface. These high-frequency modes correspond to optical modes of Si and are above the maximum allowed phonon frequency of bulk CoSi. Although the optical modes of Si contribute to phonon DOS of the interface supercell, their contribution to the Eliashberg function shown in Figure 7c is negligible. This result demonstrates that modes localized in Si do not couple with metal electrons. However, the significant volumetric coupling coefficient on the Si side of the interface in Figure 6a can be attributed to delocalized modes whose vibrational energy is distributed across Si and CoSi atoms of the interface supercell. Metal electrons transfer energy to delocalized phonon modes whose vibrational patterns dictate that a portion of the energy is transferred to silicon atoms across the interface.

Hence, our results suggest that energy exchange between electrons in metal and atomic vibrations in the semiconductor is manifested primarily by the coupling between electrons and delocalized interface modes whose vibrational energy is distributed across Si and CoSi atoms. An important implication of this result is that strength of direct electron-phonon coupling is intimately tied to the strength of interfacial bonding and the phonon-phonon conductance across the interface. For an interface with weak or van der Waals bonding, the contribution of such delocalized modes to phonon DOS is expected to be much smaller, and the phonon modes will be localized on either side of the interface.

Figure 7c also suggests that the mechanism of energy transfer from metal electrons to phonons in Si is primarily mediated by acoustic delocalized phonon modes and the contribution from coupling between electrons and optical modes of Si is negligible. This result contrasts with electron-phonon coupling in bulk Si where the contributions from acoustic and optical modes are similar in magnitude (see section VI in Supplemental Material). In the interface supercell considered here, optical modes of Si are localized to the Si side of the interface where the electron DOS at Fermi level is very small (see Figure 6b). The acoustic modes in Si are delocalized with the acoustic modes of CoSi, leading to their stronger coupling with electrons in CoSi. Figure 7c also indicates that coupling between metal electrons and CoSi optical phonon modes contributes significantly to the Eliashberg function of the interface supercell. However, such coupling is localized within the metal and contributes little to energy transfer across the interface. Localized interfacial modes, i.e., modes with vibrational energy localized to a few atomic layers around the interface are observed to contribute to the Eliashberg function in a small frequency range rad/s.

Figure 7: a) Partitioning of different regions in the Si-CoSi interface supercell used for the classification of phonon modes. b,c) Contribution of Si, CoSi, interfacial, and delocalized phonon modes to the total DOS (b) and Eliashberg function (c) of the Si-CoSi interface supercell.

v.2 Effect of Electron-Phonon Coupling on Thermal Interface Conductance

In this section, results from first-principles calculations of electron-phonon coupling are incorporated into the AGF transport simulations. The details of the approach are described in Section II.4 and Supplemental Material. The primary difference between the results presented here and those in Section IV is the presence of non-zero energy fluxes in the Büttiker probes to represent the energy exchanged between electrons and phonons. Hence, the simulation results presented in this section include contributions from both anharmonic phonon scattering and electron-phonon coupling.

We consider first the case where electrons exchange energy only through the Büttiker probes in the metal, i.e., no direct coupling between metal electrons and semiconductor phonons. The Eliashberg function of bulk strained CoSi is used to calculate the energy exchange term in Eq. (16). Figure 8a shows a typical electron and lattice temperature profile obtained from such a simulation along with the heat fluxes from the various Büttiker probes in Si and CoSi (see Figure 8c). The heat fluxes in all the Büttiker probes on the Si side of the interface are zero while the heat fluxes in the Büttiker probes of CoSi decrease away from the interface. This decay in the electron-phonon energy transfer away from the interface is a consequence of the equilibrium between electrons and phonons away from the interfacial region (see temperature profile in Figure 8a).

Also, comparison of the lattice temperature profiles in Figure 5b with that in Figure 8a shows that for the same scattering rates and applied temperature difference across the Si and CoSi contacts, the lattice temperature drop in CoSi is reduced when electrons are included in the simulation. The reduced lattice temperature drop in CoSi is a consequence of electrons in metal providing a parallel heat flow path with lower resistance compared to phonons ( W/m/K, W/m/K). Hence, a significant fraction of energy in CoSi is carried by electrons that transfer energy to the lattice near the metal-semiconductor interface.

The present simulation is conceptually similar to the analytical model developed by Majumdar and Reddy Majumdar and Reddy (2004) who suggested that electron-phonon coupling within the metal effectively provides a resistance in series with the phonon-phonon resistance across the interface. Hence the interface conductance in Figure 8a is smaller than the phonon-only conductance in Figure 5b. Majumdar and Reddy’s model for the effective conductance with electron-phonon coupling is given by:


where is the effective electron-phonon coupling efficient in the metal, is the lattice thermal conductivity of the metal, and is the phonon interfacial conductance. The electron-phonon coupling coefficient in bulk CoSi is W/m/K (see Figure 6a), W/m/K, and W/m/K (see Figure 5b). Substituting these values in Eq. (28), we obtain MW/m/K which is close to the value from the simulation in Figure 8a.

Figure 8: a) Electron and lattice temperature profile across Si-CoSi interface with electron-phonon coupling inside the metal region only. b) Electron and lattice temperature profile across Si-CoSi interface with electron-phonon coupling inside the metal region and in two unit cells of Si closest to the interface. In both a) and b), The red line corresponds to a linear fit of the lattice temperature profile in Si and the green line corresponds to a linear fit of the electron temperature profile in CoSi away from the interface region. c) Heat flux distribution in the Büttiker probes across the Si-CoSi interface corresponding to the temperature profiles in a) and b). For the simulation with direct electron-phonon coupling, the first Büttiker probe in Si closest to the interface has a non-zero energy flux.

Although the temperature profiles presented so far in Figures 5 and 8 involve conditions near room temperature, similar simulations were also performed at temperatures of 100, 150, 200, and 250 K to obtain the temperature dependence of interface conductance. At each temperature, the Büttiker probe scattering rate in Si was changed to match the bulk thermal conductivity corresponding to that temperature (see Supplemental Material). Figure 9 shows a comparison of simulation predictions with experimental measurements using the time-domain thermoreflectance (TDTR) technique Ye et al. (2016).

The simulation predictions using the various models are presented to provide a quantitative understanding of the contributions from each heat transfer mechanism to the thermal interface conductance. Ballistic AGF simulations with only coherent interface scattering (black solid curve denoted by ‘A’ in Figure 9) under-predict the thermal interface conductance for all temperatures with a 33% difference at room temperature. Also, an elastic transport model does not capture the temperature dependence of the interface conductance. Experimental data suggests that the thermal interface conductance increases by 37% from 150 K to room temperature; however the AGF simulation predicts a modest 15% increase in interface conductance for the same change in temperature. The stronger dependence of the experimental data on temperature suggests the importance of inelastic scattering processes in cross-interface energy transport. The inclusion of inelastic phonon scattering (magenta curve with circles denoted by ‘B’ in Figure 9) in the AGF simulations increases the interface conductance by about 80% at room temperature, and the simulation predictions are closer to experimental data. However, if electrons in metal are also considered in the simulation with electron-phonon coupling limited to the metal region only (red curve with hexagrams denoted by ‘C’ in Figure 9), the thermal interface conductance decreases by about 30% at room temperature, and the simulation under-predicts the experimental data. We note that this simulation considers the contributions from both anharmonic phonon scattering and electron-phonon coupling within the metal.

The DFPT calculations of electron-phonon coupling presented in the previous section do not consider anharmonicity of phonon modes in the interface supercell. In a single Si-CoSi interface with semi-infinite Si and CoSi slabs on either side, the interface phonon modes will be localized around the interface. The spatial extent of these modes will depend on the anharmonic interaction strength with bulk Si and bulk CoSi modes. The local electron-phonon coupling coefficient is expected to equal the bulk values for Si and CoSi beyond the spatial extent of these interface modes. Different approximations for the extent of joint or interface phonon modes have been proposed in the literature. Huberman and Overhauser Huberman and Overhauser (1994) proposed that the joint modes extend to a distance equal to the bulk mean free path of the materials forming the interface. For Si, the average phonon mean free path is of the order of 40 nm and the use of this length predicts a large contribution to thermal transport from cross-interface electron-phonon coupling Sadasivam et al. (2015). Results from application of the analytical model developed by Huberman & Overhauser to the present Si-CoSi interface is discussed in the Supplemental Material. More recently, Lu et al. Lu et al. (2016b) argued that the extent of interfacial phonon modes should equal the distance over which the temperature profile obtained in molecular dynamics simulations is non-linear. This length is typically of the order of 1-2 nm, and this model predicts a much smaller contribution of cross-interface electron-phonon coupling to interface conductance. In the present work, we obtain an approximate estimate of this length by fitting the simulation predictions to experimental data.

Figure 9: a) Comparison of simulation predictions with experimental data (blue squares with error bars). A (black solid curve) - Phonon-only simulation with elastic interface scattering. B (magenta circles) - Phonon-only simulation with anharmonic phonon scattering in both Si and CoSi. C (red hexagrams) - Electrons and phonons considered in the simulation with electron-phonon energy transfer inside the metal region only. D (green diamonds) - Electrons and phonons considered in the simulation with electron-phonon energy transfer included in two (1.9 nm) unit cells of Si closest to the interface.

With the assumption that cross-interface electron-phonon coupling is responsible for the difference between experimental data and the simulation results represented by the red curve in Figure 9, we use the coupling coefficient on the Si side of the interface (see Figure 6a) to model energy transfer between electrons in metal and the semiconductor lattice. Curve ‘D’ in Figure 9 represents the thermal interface conductance obtained by coupling electrons in metal with two unit cells of Si closest to the interface along the transport direction. Direct coupling with two unitcells of Si, which represents a length of approximately 1.9 nm, is found to be sufficient to obtain a close match with experimental data at various temperatures. The close match to experimental data suggests that the extent of joint interface modes in Si is much smaller than the bulk mean free path of Si. The small spatial extent of joint modes is likely due to the increased anharmonicity of interfacial phonon modes as compared to the bulk phonon modes. Similar conclusions regarding increased anharmonicity of the interfacial region are discussed in ref. Gordiz and Henry, 2016b by computing the anharmonic contribution to the potential energy of interfacial atoms in Si/Ge interfaces. The temperature profile corresponding to the simulation with direct electron-phonon coupling (see Figure 8b) is similar to that obtained from the simulation with electron-phonon coupling only in the metal region (see Figure 8a). However, the non-zero energy flux in the Büttiker probe closest to the interface in Si (see Figure 8c) is indicative of direct electron-phonon energy transfer, and this effect contributes to the enhancement in thermal interface conductance.

Vi Conclusions

This work reports first-principles calculations of phonons and electron-phonon coupling at a Si-CoSi interface and compares simulation predictions of thermal interface conductance to experimental measurements using the TDTR technique. TEM imaging of the Si-CoSi interface confirms the epitaxial nature of the interface and thus enables a quantitative comparison between simulation and experiment. From a methodology standpoint, important contributions from the present work include the development of computationally efficient methods to include inelastic phonon scattering in a Green’s function transport simulation and the incorporation of results from first-principles calculations of electron-phonon coupling into the AGF framework. We also evaluate the validity of the ‘mixing rule’, a heuristic approximation to interfacial bonding at heterojunctions, using comparisons to results obtained from rigorous first-principles calculations of interfacial bonding, and find that simple averaging of interfacial force constants can result in errors of approximately 100% in thermal interface conductance at room-temperature.

Elastic scattering of phonons at an interface is the most widely used framework to understand and predict the thermal interface conductance of heterojunctions, but the need to include inelastic phonon and coupled electron-phonon processes has become apparent, largely due to the lack of agreement between models and experiments. The present work provides a rigorous evaluation of the contributions from various transport processes for a Si-CoSi interface. Importantly, the experimental results, performed across a wide temperature range, only agree well with predictions that include all transport processes: elastic and inelastic phonon scattering, electron-phonon coupling only in the metal, and electron-phonon coupling across interface. The relative contributions of the various transport mechanisms would however be specific to the metal-semiconductor interface. For example, the extent of joint phonon modes is expected to be strongly sensitive to the strength of bonding at the interface (e.g., van der Waals vs. covalent bonding). Also, the polarity of interfacial bonds could have a significant impact on the strength of direct electron-phonon coupling. An interesting possibility for future work would involve a systematic study of the effect of interfacial bonding parameters on the relative contributions from the various cross-interface thermal transport mechanisms.


SS acknowledges financial support from the Office of Naval Research (Award No: N000141211006) and Drs. Helen and Marvin Adelberg fellowship from the School of Mechanical Engineering at Purdue University.


  1. M. H. Kryder, E. C. Gage, T. W. McDaniel, W. A. Challener, R. E. Rottmayer, G. Ju, Y.-T. Hsia,  and M. F. Erden, Proc. IEEE 96, 1810 (2008).
  2. J. P. Reifenberg, D. L. Kencke,  and K. E. Goodson, IEEE Electron Device Lett. 29, 1112 (2008).
  3. S. P. Murarka, Intermetallics 3, 173 (1995).
  4. Y. A. Kosevich, Phys. Rev. B 52, 1017 (1995).
  5. G. T. Hohensee, R. Wilson,  and D. G. Cahill, Nat. Commun. 6 (2015).
  6. A. Majumdar and P. Reddy, Appl. Phys. Lett. 84, 4768 (2004).
  7. P. E. Hopkins, J. L. Kassebaum,  and P. M. Norris, J. Appl. Phys. 105, 023710 (2009).
  8. L. Guo, S. L. Hodson, T. S. Fisher,  and X. Xu, J. Heat Transfer 134, 042402 (2012).
  9. M. Huberman and A. Overhauser, Phys. Rev. B 50, 2865 (1994).
  10. A. Sergeev, Physica B: Condens. Matter 263, 217 (1999).
  11. G. Mahan, Phys. Rev. B 79, 075408 (2009).
  12. J. Ren and J.-X. Zhu, Phys. Rev. B 87, 241412 (2013).
  13. J. Lombard, F. Detcheverry,  and S. Merabia, J. Phys.: Condens. Matter 27, 015007 (2014).
  14. T. Lu, J. Zhou, T. Nakayama, R. Yang,  and B. Li, Phys. Rev. B 93, 085433 (2016a).
  15. W. Little, Can. J. Phys. 37, 334 (1959).
  16. E. T. Swartz and R. O. Pohl, Rev. Mod. Phys. 61, 605 (1989).
  17. Z. Tian, K. Esfarjani,  and G. Chen, Phys. Rev. B 86, 235304 (2012).
  18. P. E. Hopkins, J. Appl. Phys. 106, 013528 (2009).
  19. A. Sergeev, Phys. Rev. B 58, R10199 (1998).
  20. R. Stoner and H. Maris, Phys. Rev. B 48, 16373 (1993).
  21. R. J. Stevens, A. N. Smith,  and P. M. Norris, J. Heat Transfer 127, 315 (2005).
  22. R. M. Costescu, M. A. Wall,  and D. G. Cahill, Phys. Rev. B 67, 054302 (2003).
  23. R. Wilson, B. A. Apgar, W.-P. Hsieh, L. W. Martin,  and D. G. Cahill, Phys. Rev. B 91, 115414 (2015).
  24. N. Ye, J. P. Feser, S. Sadasivam, T. S. Fisher, T. Wang, C. Ni,  and A. Janotti, arXiv preprint cond-mat/1609.01776  (2016).
  25. K. Miao, S. Sadasivam, J. Charles, G. Klimeck, T. S. Fisher,  and T. Kubis, Appl. Phys. Lett. 108, 113107 (2016).
  26. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al.J. Phys.: Condens. Matter 21, 395502 (2009).
  27. W. Zhang, T. Fisher,  and N. Mingo, Numer. Heat Transfer, Part B 51, 333 (2007).
  28. S. Sadasivam, Y. Che, Z. Huang, L. Chen, S. Kumar,  and T. S. Fisher, Annu. Rev. Heat Transfer 17, 89 (2014).
  29. M. L. Sancho, J. L. Sancho, J. L. Sancho,  and J. Rubio, J. Phys. F 15, 851 (1985).
  30. F. Guinea, C. Tejedor, F. Flores,  and E. Louis, Phys. Rev. B 28, 4397 (1983).
  31. N. Mingo, Phys. Rev. B 74, 125402 (2006).
  32. R. Venugopal, M. Paulsson, S. Goasguen, S. Datta,  and M. Lundstrom, J. Appl. Phys. 93, 5613 (2003).
  33. J. Maassen, F. Zahid,  and H. Guo, Phys. Rev. B 80, 125423 (2009).
  34. A. Afzalian, J. Appl. Phys. 110, 094517 (2011).
  35. M. Luisier, Phys. Rev. B 86, 245407 (2012).
  36. D. Singh, J. Y. Murthy,  and T. S. Fisher, J. Heat Transfer 133, 122401 (2011).
  37. N. Mingo, Phys. Rev. B 68, 113308 (2003).
  38. C. Jeong, S. Datta,  and M. Lundstrom, J. Appl. Phys. 111, 093708 (2012).
  39. M. Anantram, M. S. Lundstrom,  and D. E. Nikonov, Proc. IEEE 96, 1511 (2008).
  40. C. G. Broyden, Math. Comput. 19, 577 (1965).
  41. P. B. Allen, Phys. Rev. Lett. 59, 1460 (1987).
  42. R. Stadler, D. Vogtenhuber,  and R. Podloucky, Phys. Rev. B 60, 17112 (1999).
  43. M. Wardle, J. Goss, P. Briddon,  and R. Jones, Phys. Status Solidi A 202, 883 (2005).
  44. Z. Huang, T. Fisher,  and J. Murthy, J. Appl. Phys. 109, 074305 (2011).
  45. Z.-Y. Ong and E. Pop, Phys. Rev. B 81, 155408 (2010).
  46. N. Mingo, D. A. Stewart, D. A. Broido,  and D. Srivastava, Phys. Rev. B 77, 033418 (2008).
  47. M. N. Luckyanova, J. Garg, K. Esfarjani, A. Jandl, M. T. Bulsara, A. J. Schmidt, A. J. Minnich, S. Chen, M. S. Dresselhaus, Z. Ren, E. A. Fitzgerald,  and G. Chen, Science 338, 936 (2012).
  48. X. Gu, X. Li,  and R. Yang, Phys. Rev. B 91, 205313 (2015).
  49. E. Landry and A. McGaughey, Phys. Rev. B 80, 165304 (2009).
  50. R. Bauer, A. Schmid, P. Pavone,  and D. Strauch, Phys. Rev. B 57, 11276 (1998).
  51. S. Sadasivam, U. V. Waghmare,  and T. S. Fisher, J. Appl. Phys. 117, 134502 (2015).
  52. K. Gordiz and A. Henry, J. Appl. Phys. 119, 015101 (2016a).
  53. Z. Lu, Y. Wang,  and X. Ruan, Phys. Rev. B 93, 064302 (2016b).
  54. K. Gordiz and A. Henry, Sci. Rep. 6, 23139 (2016b).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description