Thermal Transport Across Metal SilicideSilicon Interfaces: FirstPrinciples Calculations and Green’s Function Transport Simulations
Abstract
Heat transfer across metalsemiconductor interfaces involves multiple fundamental transport mechanisms such as elastic and inelastic phonon scattering, and electronphonon 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 firstprinciples 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 CoSiSi 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 crossinterface 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. Firstprinciples calculations of electronphonon coupling reveal that crossinterface 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 CoSiSi interfaces using the timedomain 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, electronphonon coupling in the metal, and electronphonon 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 metalsemiconductor 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 highpower electronic devices. Understanding of thermal transport through metalsemiconductor 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, metalsemiconductor interfaces also provide a material system in which various physical mechanisms of heat transport such as elastic interfacial phonon scattering, inelastic phonon scattering, and electronphonon coupling coexist. 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 threephonon 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 Pbdiamond Kosevich (1995); Hohensee et al. (2015).
Electrons are the primary energy carriers in metals while phonons are dominant in intrinsic semiconductors. Hence, electronphonon coupling can be another important energy transfer mechanism that affects thermal interface conductance in metalsemiconductor interfaces. Electronphonon coupling within the metal provides an additional resistance to heat transfer Majumdar and Reddy (2004). However, electronphonon coupling across an interface, i.e., coupling between metal electrons and semiconductor phonons, provides a parallel heat flow path in addition to phononphonon heat transfer across the interface. Time domain thermoreflectance (TDTR) experiments in the literature Hopkins et al. (2009); Guo et al. (2012) suggest that direct electronphonon coupling can contribute significantly to heat transport across metalsemiconductor 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 metalsemiconductor interface are summarized in Figure 1.
Simplified empirical models are commonly used to interpret experimental thermal conductance data for metalsemiconductor 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 electronphonon 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 firstprinciples determination of the strength of coupling between electrons and joint interface phonon modes at a metalsemiconductor 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 metalsemiconductor 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 firstprinciples 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 latticematched 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 metalsemiconductor interfaces, the methods developed and reported here are also expected to be useful for a broad class of problems that use the nonequilibrium Green’s function (NEGF) method for atomistic transport simulations. From a methodology standpoint, we report a framework that combines firstprinciples 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 phononphonon 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 electronphonon 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 metalsemiconductor junctions. Rigorous firstprinciples calculations indicate that elastic phonon transport underpredicts the experimental data over a wide temperature range. Analysis of the crossinterface 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. Firstprinciples calculations of electronphonon 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 crossinterface energy transfer between metal electrons and the semiconductor lattice. We also obtain an effective length scale of electronphonon 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.
Ii Methods
ii.1 Firstprinciples Calculations
All firstprinciples 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). RappeRabeKaxirasJoannopoulos (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 PerdewBurkeErnzerhof (PBE) functional form. Three sets of firstprinciples calculations are performed for the results reported in this paper; these involve calculations on bulk Si (6 atom nonprimitive 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 kpoint grids used for DFT calculations on all three systems. Structural relaxation is carried out to reduce the HellmannFeynman 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 inplane direction to match its lattice with Si, and hence the inplane stresses are not relaxed.
Phonons are analyzed using density functional perturbation theory (DFPT) where the dynamical matrices are obtained on a qpoint 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 qpoint grid for the interface supercell. The realspace interatomic 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 kpoint grid  
Phonon qpoint grid 
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 semiinfinite ‘leads’. The device Green’s function is given by:
(1) 
where is the force constant matrix corresponding to the device region, and , are the contact selfenergies. The contact selfenergies are obtained from the surface Green’s functions , as follows:
(2) 
where , represent the force constant matrices for interaction between atoms in the device region and the semiinfinite contacts. and are the surface Green’s functions of the contacts that are obtained using the SanchoRubio method Sancho et al. (1985); Guinea et al. (1983). The phonon transmission function across the device is obtained from the Caroli formula:
(3) 
where denotes the imaginary part of the contact selfenergies 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 selfenergy, and the transmission function on the transverse phonon wavevector (for structures with periodicity in the transverse or inplane direction) and frequency is implicitly assumed. After obtaining the transmission function, the thermal interface conductance can be obtained using the Landauer formula:
(4) 
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 threedimensional 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 selfconsistent 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 selfenergy to the device Green’s function (in addition to the selfenergies due to the real contacts):
(5) 
In the present formulation, the Büttiker probe selfenergy is assumed to be a diagonal matrix whose diagonal elements are of the form:
(6) 
where denotes the Büttiker probe selfenergy at atom and vibrational direction (, , ). Similar to the matrices , , we also define that represents the imaginary part of the Büttiker probe selfenergy. 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 .
(7) 
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:
(8) 
where and denotes the spectral function. In the above equation, the matrices and are blockdiagonal. Hence, only the blockdiagonals 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 blockdiagonal 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 leftconnected is given by Anantram et al. (2008):
(9) 
where , . Our terminology follows that of ref. Anantram et al., 2008 where denotes the leftconnected retarded Green’s function, and denotes the lesser selfenergy. The derivative of with respect to the Büttiker probe temperature is given by:
(10) 
The equation for is given by:
(11) 
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).
(12)  
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:

, ,

, for
The algorithm proceeds as follows.

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

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

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

Compute the Jacobian matrix whose element is given by the following equation:
(13) where is the Kronecker Delta function.

Update the temperature of Büttiker probes using the Newton equation:
(14) 
If , go back to Step 1 with the new guess for Büttiker probe temperatures as .
An alternative to the NewtonRaphson 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 memorytime 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.
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 crossinterface electron tranport through the CoSiSi interface and consider diffusive transport of electrons in CoSi. Electrons are included in the AGF simulation within the framework of a twotemperature model that is commonly used to interpret ultrafast 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 electronelectron and phononphonon collisions that enable local equilibrium separately within electron and phonon subsystems. The electron and phonon subsystems exchange energy through electronphonon coupling that is expressed in terms of the Eliashberg function.The steady state Fourier diffusion equation for electrons in the metal is given by:
(15) 
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:
(16) 
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:
(17) 
where is the length that the Büttiker probe occupies along the transport direction. In a phonononly 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 electronphonon 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 BoseEinstein distribution at temperature and the local phonon density of states:
(18) 
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:
(19) 
Eqs. (15), (17), (18) constitute a set of coupled nonlinear 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 NewtonRaphson 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 SiCoSi interface from ballistic phonon transport calculations (i.e., Büttiker probe scattering turned off). Different interfacial atomic configurations are possible for the SiCoSi 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 SiCoSi stacking is continuous while the ‘B’ orientation occurs when the CoSi crystal is rotated by 180 about the [111] direction. Previous firstprinciples 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 firstprinciples 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 interatomic 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 crossinterface 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 SiCoSi interface formed by semiinfinite 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 crossinterface 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:
(20) 
where , denote atom indices while , denote the directions. The spatial range of interatomic interactions is artificially truncated in DFPT by the finite qpoint grid used in the calculations. Although the neglected longrange 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 postprocessing 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.
Crossinterface 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 SiGe 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 crossinterface 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 crossinterface 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 crossinterface 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 crossinterface IFCs are obtained from DFPT on the interface supercell shown in Figure 3a. In the averaging approach, the crossinterface CoSi and SiSi IFCs are obtained by averaging the interactions in bulk Si and bulk CoSi. The averaging approximation is found to overestimate 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 crossinterface 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 contactdevice interfaces in addition to the conductance of the SiCoSi interface in the middle of the device region. To obtain the conductance of the SiCoSi 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):
(21) 
where is the interface conductance computed from Eq. (4), and is the thermal interface conductance of a single SiCoSi 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 crossinterface IFCs overestimates 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 overestimation 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.
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 SiCoSi 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 crossinterface 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 highenergy 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).
V Effect of ElectronPhonon Coupling on Thermal Interface Conductance
v.1 FirstPrinciples Calculations
The results from firstprinciples calculations of electronphonon coupling, both in bulk strained CoSi and SiCoSi interface supercells, are reported in this section. The phonon linewidth due to electronphonon scattering is given by Allen (1987); Bauer et al. (1998):
(22) 
where is the electronphonon 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 electronphonon 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 electronphonon coupling:
(23) 
The spectral Eliashberg function can be used to obtain an effective volumetric electronphonon coupling coefficient :
(24) 
To understand the spatial variation of electronphonon coupling across a semiconductormetal interface, we also define a local Eliashberg function as Sadasivam et al. (2015):
(25)  
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 electronphonon coupling coefficient :
(26) 
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 kpoint 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 SiCoSi interface supercells because these supercells represent SiCoSi 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 SiCoSi superlattice. To ensure that the electronphonon coupling coefficients obtained from DFT/DFPT calculations on superlattices are transferrable to transport simulations of a single SiCoSi interface, we performed calculations on a series of SiCoSi supercells with varying Si and CoSi slab thicknesses. Figure 6a shows the spatial variation of the electronphonon 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 electronphonon coupling coefficients for bulk strained CoSi and bulk intrinsic Si are also shown in Figure 6a. The electronphonon 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 nonzero 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 SiCoSi interface forms a Schottky barrier with a ptype Schottky barrier height of 0.2 eV. Hence, the local electronic DOS at the Fermi level decays rapidly in Si away from the SiCoSi interface (see Figure 6b). However, such a twoorderofmagnitude decay in local electronic DOS does not result in a commensurate reduction in the local electronphonon 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 SiCoSi 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:
(27) 
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 crossinterface coupling between electrons in metal and phonons in the semiconductor.
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 highfrequency 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 electronphonon coupling is intimately tied to the strength of interfacial bonding and the phononphonon 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 electronphonon 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.
v.2 Effect of ElectronPhonon Coupling on Thermal Interface Conductance
In this section, results from firstprinciples calculations of electronphonon 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 nonzero 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 electronphonon 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 electronphonon 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 metalsemiconductor interface.
The present simulation is conceptually similar to the analytical model developed by Majumdar and Reddy Majumdar and Reddy (2004) who suggested that electronphonon coupling within the metal effectively provides a resistance in series with the phononphonon resistance across the interface. Hence the interface conductance in Figure 8a is smaller than the phonononly conductance in Figure 5b. Majumdar and Reddy’s model for the effective conductance with electronphonon coupling is given by:
(28) 
where is the effective electronphonon coupling efficient in the metal, is the lattice thermal conductivity of the metal, and is the phonon interfacial conductance. The electronphonon 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.
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 timedomain 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) underpredict 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 crossinterface 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 electronphonon 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 underpredicts the experimental data. We note that this simulation considers the contributions from both anharmonic phonon scattering and electronphonon coupling within the metal.
The DFPT calculations of electronphonon coupling presented in the previous section do not consider anharmonicity of phonon modes in the interface supercell. In a single SiCoSi interface with semiinfinite 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 electronphonon 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 crossinterface electronphonon coupling Sadasivam et al. (2015). Results from application of the analytical model developed by Huberman & Overhauser to the present SiCoSi 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 nonlinear. This length is typically of the order of 12 nm, and this model predicts a much smaller contribution of crossinterface electronphonon coupling to interface conductance. In the present work, we obtain an approximate estimate of this length by fitting the simulation predictions to experimental data.
With the assumption that crossinterface electronphonon 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 electronphonon coupling (see Figure 8b) is similar to that obtained from the simulation with electronphonon coupling only in the metal region (see Figure 8a). However, the nonzero energy flux in the Büttiker probe closest to the interface in Si (see Figure 8c) is indicative of direct electronphonon energy transfer, and this effect contributes to the enhancement in thermal interface conductance.
Vi Conclusions
This work reports firstprinciples calculations of phonons and electronphonon coupling at a SiCoSi interface and compares simulation predictions of thermal interface conductance to experimental measurements using the TDTR technique. TEM imaging of the SiCoSi 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 firstprinciples calculations of electronphonon 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 firstprinciples 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 roomtemperature.
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 electronphonon 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 SiCoSi 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, electronphonon coupling only in the metal, and electronphonon coupling across interface. The relative contributions of the various transport mechanisms would however be specific to the metalsemiconductor 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 electronphonon 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 crossinterface thermal transport mechanisms.
Acknowledgements
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.
References
 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).
 J. P. Reifenberg, D. L. Kencke, and K. E. Goodson, IEEE Electron Device Lett. 29, 1112 (2008).
 S. P. Murarka, Intermetallics 3, 173 (1995).
 Y. A. Kosevich, Phys. Rev. B 52, 1017 (1995).
 G. T. Hohensee, R. Wilson, and D. G. Cahill, Nat. Commun. 6 (2015).
 A. Majumdar and P. Reddy, Appl. Phys. Lett. 84, 4768 (2004).
 P. E. Hopkins, J. L. Kassebaum, and P. M. Norris, J. Appl. Phys. 105, 023710 (2009).
 L. Guo, S. L. Hodson, T. S. Fisher, and X. Xu, J. Heat Transfer 134, 042402 (2012).
 M. Huberman and A. Overhauser, Phys. Rev. B 50, 2865 (1994).
 A. Sergeev, Physica B: Condens. Matter 263, 217 (1999).
 G. Mahan, Phys. Rev. B 79, 075408 (2009).
 J. Ren and J.X. Zhu, Phys. Rev. B 87, 241412 (2013).
 J. Lombard, F. Detcheverry, and S. Merabia, J. Phys.: Condens. Matter 27, 015007 (2014).
 T. Lu, J. Zhou, T. Nakayama, R. Yang, and B. Li, Phys. Rev. B 93, 085433 (2016a).
 W. Little, Can. J. Phys. 37, 334 (1959).
 E. T. Swartz and R. O. Pohl, Rev. Mod. Phys. 61, 605 (1989).
 Z. Tian, K. Esfarjani, and G. Chen, Phys. Rev. B 86, 235304 (2012).
 P. E. Hopkins, J. Appl. Phys. 106, 013528 (2009).
 A. Sergeev, Phys. Rev. B 58, R10199 (1998).
 R. Stoner and H. Maris, Phys. Rev. B 48, 16373 (1993).
 R. J. Stevens, A. N. Smith, and P. M. Norris, J. Heat Transfer 127, 315 (2005).
 R. M. Costescu, M. A. Wall, and D. G. Cahill, Phys. Rev. B 67, 054302 (2003).
 R. Wilson, B. A. Apgar, W.P. Hsieh, L. W. Martin, and D. G. Cahill, Phys. Rev. B 91, 115414 (2015).
 N. Ye, J. P. Feser, S. Sadasivam, T. S. Fisher, T. Wang, C. Ni, and A. Janotti, arXiv preprint condmat/1609.01776 (2016).
 K. Miao, S. Sadasivam, J. Charles, G. Klimeck, T. S. Fisher, and T. Kubis, Appl. Phys. Lett. 108, 113107 (2016).
 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).
 W. Zhang, T. Fisher, and N. Mingo, Numer. Heat Transfer, Part B 51, 333 (2007).
 S. Sadasivam, Y. Che, Z. Huang, L. Chen, S. Kumar, and T. S. Fisher, Annu. Rev. Heat Transfer 17, 89 (2014).
 M. L. Sancho, J. L. Sancho, J. L. Sancho, and J. Rubio, J. Phys. F 15, 851 (1985).
 F. Guinea, C. Tejedor, F. Flores, and E. Louis, Phys. Rev. B 28, 4397 (1983).
 N. Mingo, Phys. Rev. B 74, 125402 (2006).
 R. Venugopal, M. Paulsson, S. Goasguen, S. Datta, and M. Lundstrom, J. Appl. Phys. 93, 5613 (2003).
 J. Maassen, F. Zahid, and H. Guo, Phys. Rev. B 80, 125423 (2009).
 A. Afzalian, J. Appl. Phys. 110, 094517 (2011).
 M. Luisier, Phys. Rev. B 86, 245407 (2012).
 D. Singh, J. Y. Murthy, and T. S. Fisher, J. Heat Transfer 133, 122401 (2011).
 N. Mingo, Phys. Rev. B 68, 113308 (2003).
 C. Jeong, S. Datta, and M. Lundstrom, J. Appl. Phys. 111, 093708 (2012).
 M. Anantram, M. S. Lundstrom, and D. E. Nikonov, Proc. IEEE 96, 1511 (2008).
 C. G. Broyden, Math. Comput. 19, 577 (1965).
 P. B. Allen, Phys. Rev. Lett. 59, 1460 (1987).
 R. Stadler, D. Vogtenhuber, and R. Podloucky, Phys. Rev. B 60, 17112 (1999).
 M. Wardle, J. Goss, P. Briddon, and R. Jones, Phys. Status Solidi A 202, 883 (2005).
 Z. Huang, T. Fisher, and J. Murthy, J. Appl. Phys. 109, 074305 (2011).
 Z.Y. Ong and E. Pop, Phys. Rev. B 81, 155408 (2010).
 N. Mingo, D. A. Stewart, D. A. Broido, and D. Srivastava, Phys. Rev. B 77, 033418 (2008).
 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).
 X. Gu, X. Li, and R. Yang, Phys. Rev. B 91, 205313 (2015).
 E. Landry and A. McGaughey, Phys. Rev. B 80, 165304 (2009).
 R. Bauer, A. Schmid, P. Pavone, and D. Strauch, Phys. Rev. B 57, 11276 (1998).
 S. Sadasivam, U. V. Waghmare, and T. S. Fisher, J. Appl. Phys. 117, 134502 (2015).
 K. Gordiz and A. Henry, J. Appl. Phys. 119, 015101 (2016a).
 Z. Lu, Y. Wang, and X. Ruan, Phys. Rev. B 93, 064302 (2016b).
 K. Gordiz and A. Henry, Sci. Rep. 6, 23139 (2016b).