Understanding phonon transport across heterojunctions is important to achieve a wide range of thermal transport properties. Using the McKelvey-Shockley flux method with first-principles modeling, we theoretically investigate the phonon transport properties of a Si–Ge interface with a focus on the role of inelastic bulk phonon processes. We observe significant inelastic scattering near the interface that redistributes the heat among the phonons as a result of non-equilibrium effects driven by the junction. These effects are most pronounced when the length of the junction is comparable to the average phonon mean-free-path. What controls these inelastic processes is elucidated.

Developing material heterojunctions with a wide range of thermal transport properties is the focus of intense research, which can directly impact applications related to electronics, heat management, and thermoelectrics.^{1,2} Although this topic has a long history, there still remain open fundamental questions regarding what microscopic processes control heat transport and dictate the thermal contact resistance, *R*_{C}.

Focusing on the case of interfaces composed of semiconductors and/or insulators, in which phonons are the sole heat carriers, there have been many theoretical studies using a variety of methods to understand phonon transport and predict *R*_{C}. Simple models include the acoustic mismatch model (AMM) and diffuse mismatch model (DMM),^{3} where the AMM assumes wave coherence and specular scattering off a smooth interface and the DMM assumes fully diffuse scattering off a rough interface. Thus, together both provide a range of *R*_{C} values that can provide information about the quality of the interface when compared to experiment, although measurements do not always fall within this range.

Interfacial phonon transport has also been investigated using molecular dynamics (MD),^{4–21} atomistic Green’s function (AGF),^{22–26} and the Boltzmann transport equation (BTE).^{27–29} Each approach has its advantages and has been used to identify various phonon transport processes across heterojunctions. For example, MD captures all-order anharmonic phonon scattering and was used to demonstrate the role of inelastic multi-phonon processes across the interface.^{13,15,19} AGF and MD naturally treat the interfacial phonon properties, which lead to non-bulk interface modes that have been proposed to mediate phonon transport.^{16,17} Detailed mode-resolved phonon transport properties^{18} and transmission coefficients^{25} have been obtained using AGF and MD. The BTE showed that inelastic bulk scattering can result in non-trivial temperature and heat current distributions near an interface.^{28} This may be a general feature of heterojunctions, but what exactly controls such inelastic processes and how this affects phonon transport have yet to be fully understood.

When there is a mismatch in the phonon energies, for example, with Si and Ge as depicted in Fig. 1, inelastic phonon scattering is likely an important process. For example, since the maximum Si phonon energy (*ϵ*_{Si,max}) is nearly double that of Ge (*ϵ*_{Ge,max}), the high-energy Si phonons cannot conduct heat without inelastic scattering.

In this study, we investigate a Si–Ge interface focusing on the role and origin of inelastic bulk scattering in Si and Ge, using the McKelvey-Shockley flux method combined with density functional theory (DFT) and an ideal interface model. Our approach has the advantage of being computationally efficient and physically transparent, which helps provide new insights. We find that the high-energy Si phonons do carry significant heat (more than in bulk) due to strong inelastic processes activated by non-equilibrium effects driven by the junction. While Si–Ge interfaces have been extensively investigated,^{4–8,14–17,19,20,22,25,29–31} this work builds on prior studies^{24,28} to further elucidate the role of inelastic bulk scattering, a process that may be important in a broader class of heterojunctions.

We consider a single Si–Ge interface of length *L*, with *L*/2 of both Si and Ge, as shown in Fig. 1. The interface is ideal, meaning atomically smooth with no strain or roughness present. The left side of Si and the right side of Ge are joined by two ideal contacts that are in thermodynamic equilibrium at different temperatures, *T*_{L} and *T*_{R}. The contacts are assumed to be perfectly absorbing; any phonon in Si or Ge reaching a contact will not be reflected. In this study, we take *T*_{L} = 301 K and *T*_{R} = 300 K such that the phonon transport is induced by a Δ*T* = 1 K.

The treatment of phonon transport inside the Si and Ge is carried out using the McKelvey-Shockley flux (McK-S) method. Originally developed for electron transport,^{32,33} it was subsequently adapted for phonon/thermal transport.^{34,35} The McK-S approach is a simple version of the BTE that captures ballistic and non-equilibrium effects and can treat transport from the nano-scale to the macro-scale.^{34–38} Importantly, this method supports inelastic phonon scattering.^{38}

The McK-S describes the spatial and temporal evolution of the forward- and reverse-moving phonon fluxes (along the transport direction, taken here as *x*). Instead of dealing with directed fluxes, it is possible to rewrite the governing equations directly in terms of temperature and heat current.^{38} Temperature, in steady-state, is given by

where *T*(*x*,*ϵ*) is the temperature resolved at each phonon energy, *T*_{0}(*x*) is the temperature of a fictitious local equilibrium distribution (originating from the use of the relaxation time approximation), $D(\u03f5)=\lambda (\u03f5)vx+(\u03f5)/2$ is the *ϵ*-resolved diffusion coefficient, $\tau Q(\u03f5)=\lambda (\u03f5)/(2vx+(\u03f5))$ is the *ϵ*-resolved heat relaxation time, $CV(\u03f5)=\u03f5\u2009g(\u03f5)\u2009[\u2202nBE/\u2202T]Tref$ is the *ϵ*-resolved heat capacity, *λ*(*ϵ*) is the phonon mean-free-path for backscattering (average *x*-projected distance traveled before backscattering), $vx+(\u03f5)$ is the average *x*-projected phonon velocity, *g*(*ϵ*) is the phonon density of states, *n*_{BE}(*ϵ*, *T*_{ref}) is the Bose-Einstein distribution, and *T*_{ref} is the reference temperature (300 K in this work). *λ*(*ϵ*) and $vx+(\u03f5)$ are defined in Ref. 40.

The expression for *T*_{0}(*x*) is obtained by imposing conservation of energy (see Ref. 38). In the form given by Eq. (2), inelastic scattering is permitted (phonons can change energy during scattering). In some previous studies with the McK-S method,^{34–37} only elastic scattering was permitted (phonons could not change energy during scattering); a discussion comparing elastic and inelastic scattering can be found in Ref. 38. The treatment of scattering in this framework is at the level of the relaxation time approximation, which requires a scattering time to be specified; in this case, 3-phonon anharmonic scattering is calculated from first-principles.

Heat current, in steady-state, is expressed as

where *κ*(*ϵ*) = *C*_{V}(*ϵ*)*D*(*ϵ*) is the *ϵ*-resolved bulk thermal conductivity.

The boundary conditions (BCs) needed to solve for *T*(*x*, *ϵ*) at the left/right contacts are

These mixed BCs are necessary to capture ballistic and non-equilibrium phonon effects^{34–37} and retrieve the traditional BCs [i.e., *T*(−*L*/2, *ϵ*) = *T*_{L} and *T*(*L*/2, *ϵ*) = *T*_{R}] in the diffusive limit. Separate BCs are adopted for the Si–Ge junction (discussed below).

The effective, energy-integrated temperature and heat current are obtained using^{38}

This approach was shown to compare well to the more rigorous and computationally demanding BTE, particularly when solved using a full phonon dispersion.^{34–38}

To treat the Si–Ge interface, we adopt a simple interface model intended to capture the minimal thermal contact resistance between two dissimilar materials, assuming phonons transfer elastically across the junction. The motivation of this model, hereinafter referred to as the ideal interface model (IIM), is not to achieve better agreement with experiment but rather to (approximately) include the fundamental interfacial scattering that should exist in an ideal, smooth junction. A derivation, provided in Appendix A, is based on minimal assumptions: conservation of energy and detailed balance. AGF can provide a more accurate assessment of the fundamental transmission coefficient^{25} but is more computationally expensive, particularly if both materials are not well lattice matched. Non-idealities (e.g., interfacial strain or defects) and other processes not considered here (e.g., inelastic interface scattering^{39}) will alter the transport properties. We refer interested readers to Appendix A for further details of the IIM and how it relates to the AMM and DMM.

With the IIM, the transmission coefficients are simply

where $Ti\u2192f(\u03f5)$ is the transmission probability of a phonon with energy *ϵ* traveling from material *i* to *f*, and *M*_{i/f}(*ϵ*) is the distribution of modes (also referred to as the number of conducting channels per cross-sectional area) of the phonon states in material *i*/*f*. Phonons traveling from a material with few channels to a material with more channels (*M*_{f} > *M*_{i}) have a transmission of one, since there are enough final phonon states to accommodate all incoming phonons. However, the transmission is below one for the phonons impinging on a material with fewer number of channels (*M*_{f} < *M*_{i}), since there are not sufficient final states. See Ref. 40 for the definition of *M*(*ϵ*).

The transmission coefficients calculated from Eq. (8) are used in the following expression for the *ϵ*-resolved thermal contact resistance (see Appendix A):

where $Rball(\u03f5)=CV(\u03f5)vx+(\u03f5)/2$ is the ballistic thermal resistance. Equation (9) is relatively general, assuming only elastic transport across the interface and detailed balance, and handles non-equilibrium phonon distributions. The IIM retrieves two known limits: *R*_{C} = 0 in the case of two identical materials and *R*_{C}(*ϵ*) = *R*^{ball}/2 when a material is joined with an ideal Landauer contact.

Using the IIM, given by Eqs. (8) and (9), we obtain the following boundary conditions for the Si–Ge interface, assumed here to be at *x* = *l*:

which states that the contact resistance will result in a discrete temperature drop across the interface Δ*T*(*l*, *ϵ*) = *T*(*l*^{−}, *ϵ*) − *T*(*l*^{+}, *ϵ*) when a heat current flows and that the heat current is continuous at the interface, respectively. The condition of continuous heat current across the interface, Eq. (11), implies that there is no heat generation at the interface. Note that the effective thermal contact resistance is not obtained by an integration of *R*_{C}(*ϵ*) (or its inverse) over energy but rather is extracted using *R*_{C} = Δ*T*(*l*)/*I*_{Q}(*l*) from the *ϵ*-integrated, effective temperature and heat current.

DFT was utilized to compute the phonon dispersions and 3-phonon anharmonic scattering rates for Si and Ge, which serve as input for our Si–Ge interface calculations (details of the DFT modeling are found in Appendix B). Figure 2(a) shows the calculated phonon dispersions, *ϵ*(*q*), for Si and Ge, which agree well with previous calculations.^{41} The maximum phonon energies in Si and Ge are *ϵ*_{Si,max} = 62 meV and *ϵ*_{Ge,max} = 34 meV, respectively. From the *ϵ*(*q*), we compute the number of conducting channels per cross section, *M*_{Si,Ge}(*ϵ*), presented in Fig. 2(b). At low energies, the *M*_{Si,Ge}(*ϵ*) scale quadratically, as expected for linear dispersions. Figure 2(c) shows the energy-dependent thermal contact resistance calculated using Eqs. (9) and (8). *R*_{C}(*ϵ*) is largest for the low-energy acoustic modes and decreases sharply until ∼10 meV above which the resistance is relatively similar. The mean-free-path for backscattering, *λ*(*ϵ*), is extracted and presented in Fig. 2(d). *λ*(*ϵ*) is longest for the acoustic phonons, where 3-phonon scattering is less likely, reaching ∼30 *μ*m before decaying rapidly to values below 100 nm for *ϵ* > 10 meV. A comparison of this DFT computed *λ*(*ϵ*) for Si to the previously reported *λ*(*ϵ*) in Ref. 34, based on phenomenological expressions for scattering, shows that both agree well and yield similar magnitude and energy dependence. The bulk thermal conductivities are calculated to be 160 W/m-K for Si and 51 W/m-K for Ge at 300 K, which compare well to previous simulations.^{42}

The case of *L* = 200 nm is interesting since some high-energy phonons are diffusive (*λ*(*ϵ*) ≪ *L*), while simultaneously some other low-energy phonons are ballistic (*λ*(*ϵ*) ≫ *L*). Figure 3(a) shows the temperature of the Si–Ge junction versus position, *x*, and phonon energy, *ϵ*. Above 34 meV, the temperature on the Ge side is not defined since there are no phonon states. Most of the temperature drop occurs on the Ge side compared to the Si side, since the former has a lower thermal conductivity. At the interface, there is a discrete temperature drop that varies with energy; it is largest for the low-energy acoustic phonons, where *R*_{C}(*ϵ*) is high [see Fig. 2(c)]. The low-energy acoustic phonons have a flat temperature profile, accompanied with a large temperature drop at the left/right contacts, which is a signature of ballistic phonon transport (as discussed in Refs. 34–37). The high-energy Si phonons have a temperature at the *x* = −*L*/2 boundary approaching *T*_{L}, indicating diffusive phonon scattering. Figure 3(b) shows the temperature drop in various sections of the structure versus energy. At low *ϵ*, Δ*T* is mostly localized at the interface/contacts (ballistic transport), while at higher *ϵ*, much of the Δ*T* occurs across the Si and mostly the Ge (diffusive transport).

Using Eq. (6), we calculate the effective, energy-integrated temperature profile, plotted as a solid black line in Fig. 3(a). We find a temperature drop of 0.24 K, which is relatively large compared to *T*_{L} − *T*_{R} = 1 K. The total heat current is *x*-independent and equal to *I*_{Q} = 1.1 × 10^{8} W/m^{2}, resulting in an effective *R*_{C} = 2.2 × 10^{−9} km^{2}/W with the IIM (3.5 × 10^{−9} km^{2}/W with the DMM). These values are close to other reported calculations using the DMM (4.6 × 10^{−9} km^{2}/W),^{31} the AMM (≈5.5 × 10^{−9} km^{2}/W),^{31} and MD (≈3.1 × 10^{−9} km^{2}/W).^{5} The IIM yields the lowest *R*_{C}. Surprisingly, the reported AMM value is larger than that of the DMM. A traditional solution of the 1D heat equation gives a linear temperature profile; however, we find a non-linear *T*(*x*), particularly across the Ge. As discussed later, this is a result of strong inelastic scattering near the interface.

Figure 4(a) shows the heat current versus *x* and *ϵ*. *I*_{Q}(*x*, *ϵ*) flows at different energies in Si and Ge, which is only possible via inelastic scattering; elastic scattering would give an *x*-independent *I*_{Q}(*ϵ*). Similar features were observed in Ref. 28. At the interface, there are select energies, near 13 meV and 26 meV, which carry much of the heat from Si to Ge. These energies are found to have roughly equal number of conducting channels in both materials, *M*_{Si}(*ϵ*) ≈ *M*_{Ge}(*ϵ*), which results in small *R*_{C}(*ϵ*).

Focusing on the high-energy Si phonons, a non-zero *I*_{Q}(*x*, *ϵ*) is observed near 50 meV. Figure 4(b) shows the fraction of the total heat current carried above the maximum Ge energy, which reaches over 20% at the left boundary and decreases to zero at the interface (the fraction of heat current flowing above *ϵ*_{Ge,max} in bulk Si is 8%, shown as a dashed line). Thus, a significant amount of heat is carried by the high-energy Si phonons, even though they are completely reflected at the interface. The fact that *I*_{Q}(*x*, *ϵ*) varies with *x* indicates that inelastic processes are rearranging the heat among the different phonons (i.e., the occupation of the phonon states), which we will illustrate next.

With no external sink/source removing or adding heat to the structure, we obey the continuity equation for heat: *C*_{V} *∂T*(*x*)/*∂t* = −*∂I*_{Q}(*x*)/*∂x*. We can also write a similar equation that is phonon energy resolved

in which we include a “sink-source” term, *S*(*x*, *ϵ*), describing how heat is being removed or added at each phonon energy, *ϵ*. It can be shown that $\u222b0\u221eS(x,\u03f5)\u2009d\u03f5=0$, stating that heat is only redistributed among the phonons and that no net heat is removed or added. In steady-state, the left-hand side of Eq. (12) is zero, giving *S*(*x*, *ϵ*) = *∂I*_{Q}(*x*, *ϵ*)/*∂x*. Using Eqs. (1) and (3), we obtain

Figure 5(a) presents *S*(*x*, *ϵ*) calculated using Eq. (13). The positive (red) and negative (blue) values correspond to heat being added and removed, respectively. In the vicinity of the interface, there is a strong redistribution of heat, i.e., inelastic scattering. On the Si side, heat is removed from the high-energy phonons and transferred to the mid-energy phonons (∼30 meV) that can travel into the Ge. Roughly speaking, as the high-energy phonons approach, they “sense” the interface and pump heat down to lower-energy phonons that can cross the junction. Since certain phonons transport most of the heat across the interface from Si to Ge, as shown in Fig. 4(a), this creates a highly non-equilibrium distribution in the Ge side of the junction. As an attempt to drive the system towards equilibrium, in the Ge, we observe strong inelastic scattering filling depleted states and emptying others. The redistribution of heat near an interface occurs over a length scale comparable to *λ*(*ϵ*). We also see similar features near the left/right contacts. As previously discussed,^{34,35} the phonons injected from the contacts are in equilibrium with their originating thermal reservoir, but the phonons impinging on the contacts have a distribution controlled by the scattering that occurred inside the structure. The resulting phonon population near the contacts may be far from equilibrium, thus triggering inelastic processes.

To better understand what controls the heating/cooling pattern in Fig. 5(a), we present in Fig. 5(b) the energy-resolved temperature, *T*(0, *ϵ*) (solid lines), and local equilibrium temperature, *T*_{0}(0) (dashed lines), on both sides of the interface (*x* = 0^{−} for Si and *x* = 0^{+} for Ge). According to Eq. (13), heat redistribution occurs whenever *T*(*x*, *ϵ*) ≠ *T*_{0}(*x*); heat is added when *T*_{0}(*x*) − *T*(*x*, *ϵ*) > 0 and removed when *T*(*x*, *ϵ*) − *T*_{0}(*x*) > 0. With this insight, Fig. 5(b) shows which phonon energies will be heated/cooled, which correlates with the pattern in Fig. 5(a). The low-energy acoustic modes are far from the equilibrium temperature because these are near ballistic and do not efficiently scatter; however, their value of *C*_{V}(*ϵ*)/*τ*_{Q}(*ϵ*) is small. Note that under near-equilibrium conditions *T*(*x*, *ϵ*) is independent of *ϵ*; thus, an energy dependence is indicative of a non-equilibrium phonon distribution. In summary, interfaces drive the phonon population out of equilibrium, which in turn activate inelastic scattering.

Figure 6(a) shows the temperature profile across the Si–Ge junction with *L* ranging from 20 nm to 200 *μ*m. With *L* = 20 nm, most of the temperature drop occurs at the left/right contacts and the interface, indicating near-ballistic transport. With *L* = 200 *μ*m, we retrieve the traditional heat equation solution, a linear *T*(*x*) with vanishing temperature drops at the interface and contacts. Interestingly, a distinct non-linearity is observed with *L* = 200 nm, since this is the length scale over which inelastic scattering is strongest. With smaller *L*, scattering is sparse leading to a linear and flat *T*(*x*) when approaching the ballistic limit, as discussed in Ref. 34. With larger *L*, most of the scattering occurs close to the interface/contacts and the phonons inside each material are near equilibrium resulting in a linear *T*(*x*), as predicted by the classical heat equation. The non-linearity in *T*(*x*), observed at an intermediate length, arises from inelastic scattering, since assuming only elastic scattering (using the same approach) results in a linear *T*(*x*).^{34}

The heat current versus *L* is plotted in Fig. 6(b), comparing the IIM to the DMM. For *L* > 20 *μ*m, the heat current follows a *L*^{−1} trend, as expected in the diffusive limit. Below *L* = 1 *μ*m, the heat current rolls off, falling below the traditional scaling, as phonon transport becomes quasi-ballistic and approaches the ballistic limit. Figure 6(c) presents the energy-resolved heat current for various *L*. Note how each plot is visually different until reaching ∼20 *μ*m, above which *I*_{Q}(*x*, *ϵ*) for a given *x* is the same as bulk. With *L* = 20 nm, the heat current flows at roughly the same energies in both Si and Ge, since inelastic scattering occurs over a longer length scale (as discussed above). With *L* = 200 nm, we observe the transition of heat current carried by certain energies in Si to others in Ge. For longer *L*, the energy-dependence of the heat current is mostly uniform in *x* (away from the interface and contacts).

In summary, we theoretically analyzed phonon transport across a Si–Ge interface using the McKelvey-Shockley flux method combined with DFT materials modeling, which captures ballistic, non-equilibrium, and inelastic scattering effects. This study focused mainly on the role of inelastic bulk phonon scattering that occurs near a heterojunction. We started by investigating a 200 nm long Si–Ge junction, in which we simultaneously observed diffusive transport of high-energy phonons and ballistic transport of low-energy phonons. Our analysis showed that strong inelastic scattering exists near the interface and the contacts. This originates from the non-equilibrium phonon distribution that is induced by the presence of the interface; heat is conducted across the junction at select phonon energies while other phonons are completely reflected (e.g., the high-energy Si phonons), which drive the system away from equilibrium.

As a result of these inelastic processes, there is a significant redistribution of heat among the phonons near an interface. This enables the high-energy Si phonons to carry significant heat, which are “cooled” within roughly one mean-free-path of the junction by pumping heat down to the mid-energy Si phonons that can transfer into the Ge. What activates and controls these processes was elucidated. Moreover, a non-linear temperature profile is observed in the regions where inelastic scattering is significant. We also demonstrated how the phonon transport properties vary with the length of the junction, *L*, from 20 nm to 200 *μ*m. We observe the ballistic to diffusive transition and show how this impacts the temperature and heat current distributions.

For our treatment of the Si–Ge interface, we adopted a simple interface model intended to (approximately) capture the minimum thermal resistance of an ideal smooth junction. The physics of this model, and its connection to others including the AMM and DMM, were discussed, and an expression for thermal contact resistance was provided.

Finally, while our treatment is not comprehensive, for example, excluding interfacial phonon modes or inelastic phonon transport across the interface, it does highlight the importance of inelastic bulk scattering and non-equilibrium physics, which are processes that likely play a role in many other heterointerfaces beyond Si–Ge.

This work was partially supported by DARPA MATRIX (Award No. HR0011-15-2-0037) and NSERC (Discovery Grant No. RGPIN-2016-04881). Computational resources were provided by Compute Canada.

### APPENDIX A: THERMAL CONTACT RESISTANCE AND THE IDEAL INTERFACE MODEL

Consider a heterojunction at *x* = *l*, with the material to the left (*x* ≤ *l*^{−}) labeled as 1 and the material to the right (*x* ≥ *l*^{+}) labeled as 2. Assuming phonons travel elastically across the interface, a scattering matrix relating the incoming heat fluxes to the outgoing heat fluxes can be constructed

where $IQ\xb1(x,\u03f5)$ are the directed heat fluxes (calculated with the McK-S method in this study) and $Ti\u2192f(\u03f5)$ is the transmission probability of a phonon traveling from the initial material *i* to the final material *f*. Note that in this section the explicit *ϵ* dependence of quantities will often be dropped for clarity (but is implicitly assumed). Using this scattering matrix, we first derive an expression for thermal contact resistance, *R*_{C}, and then discuss the particular choice in transmission coefficients that make up the ideal interface model (IIM).

We begin by relating the directed heat fluxes to the net heat current and temperature, and deriving several useful properties. The directed heat fluxes are written as $IQ\xb1(x)=IQ,ref\xb1+\delta IQ\xb1(x)$, where $IQ,ref+=IQ,ref\u2212$ are the heat fluxes of an equilibrium phonon distribution at a constant reference temperature *T*_{ref} and $\delta IQ\xb1(x)$ are the corrections to the background fluxes. The heat current is $IQ(x)=IQ+(x)\u2212IQ\u2212(x)$, and the temperature is *T*(*x*) = (*T*^{+}(*x*) + *T*^{−}(*x*))/2, where $T\xb1(x)=[IQ\xb1(x)\u2212IQ,ref\xb1]/Kball+Tref$ is the temperature of the forward/reverse phonon streams, $Kball(\u03f5)=CVvx+/2=\u03f5M(\u2202nBE/\u2202T)|Tref/h$ is the ballistic thermal conductance, *C*_{V}(*ϵ*) is the heat capacity, $vx+(\u03f5)$ is the average *x*-projected velocity (*x* being the transport direction), *M*(*ϵ*) is the distribution of modes, and *n*_{BE}(*ϵ*, *T*) is the Bose-Einstein distribution (see Ref. 34 for more details on these quantities). We assume small temperature variations such that |*T*(*x*) − *T*_{ref}|≪ *T*_{ref}.

Next, a condition that no net heat current flows in equilibrium is imposed on Eq. (A1). Using Eq. (A1), we find

where $IQ\xb1\u2192IQ,ref\xb1$ in equilibrium. Using $IQ,ref\xb1=\u03f5MnBE/h$ leads to the following condition relating $T1\u21922$ and $T2\u21921$:

or equivalently,

Note that this same condition is obtained if we had started with *I*_{Q}(*l*^{−}) = 0 instead of *I*_{Q}(*l*^{+}) = 0. Equations (A3) and (A4) indicate that only $T1\u21922$ must be specified and that $T2\u21921$ is automatically determined by detailed balance (or vice versa).

By adding both equations defined by (A1), one finds *I*_{Q}(*l*^{−}) = *I*_{Q}(*l*^{+}), stating that the heat current is continuous across the interface. By subtracting $IQ\u2212(l+)$ from the first equation in (A1), the following equivalent expressions for heat current are found:

where Eq. (A4) was used.

We can now derive an expression for thermal contact resistance, *R*_{C}. By adding $IQ\u2212(l+)$ and $IQ+(l\u2212)$ to both sides of the first and second equations in (A1), respectively, and manipulating, we arrive at expressions for temperature on both sides of the junction

The discrete temperature drop across the interface is Δ*T* = *T*(*l*^{−}) − *T*(*l*^{+}). Subtracting Eq. (A9) from Eq. (A10), and using Eqs. (A7) and (A8), we find

Defining the thermal contact resistance as *R*_{C} = Δ*T*/*I*_{Q}(*l*^{±}), we obtain our final expression

where *R*^{ball}(*ϵ*) = 1/*K*^{ball}(*ϵ*) is the ballistic thermal resistance. Note that Eq. (A12) is rather general; its derivation is based on the assumption of elastic transfer of phonons across the interface and the principle of detailed balance. No assumption on the phonon dispersion was made, and this expression can be easily evaluated with a first-principles full phonon dispersion. Finally, Eq. (A12) is applicable under non-equilibrium phonon conditions.

Next we specify the particular choice in $Ti\u2192f$ that will constitute the IIM. The purpose of this model is to (approximately) capture the minimal thermal contact resistance between two dissimilar materials under the assumption of an ideal, smooth interface. For this purpose, we ask ourselves what is the largest possible interfacial transmission coefficient? A transmission of one is allowed, but only for phonons flowing along one direction. Assuming that *M*_{1} < *M*_{2}, then $T1\u21922=1$ is permitted, but according to Eq. (A3), $T2\u21921=M1/M2<1$ by detailed balance. $T2\u21921=1$ is not allowed, since this would result in $T1\u21922=M2/M1>1$. If we now assume that *M*_{1} > *M*_{2}, then $T2\u21921=1$ and $T1\u21922=M2/M1<1$. When *M*_{1} = *M*_{2}, both transmissions are one, $T1\u21922=T2\u21921=1$.

This model retrieves two known limits: (1) *Two identical materials*. In the case of two identical materials, the IIM states that $T1\u21922=T2\u21921=1$, and Eq. (A12) gives *R*_{C} = 0, as one would physically expect. (2) *Material interfaced with a Landauer contact*. An ideal thermal reservoir, sometimes referred to as a Landauer contact, is assumed to have a very large number of channels compared to the material. Labeling the material as 1 and the Landauer contact as 2, a careful evaluation of Eq. (A12) yields $RC\u2192R1ball/2$ when *M*_{1} ≪ *M*_{2}, a well-known result from the Landauer approach for phonon^{34} and electron transport.^{43}

As an example, let us consider an interface formed by two materials with linear phonon dispersions, *E*(*q*) = *ℏvq*, where *v* is the group velocity (here we assume *v*_{1} < *v*_{2}) and *q* is the phonon wavenumber [see Fig. 7(a)]. The distribution of modes for a linear dispersion is *M*(*ϵ*) = 3*πϵ*^{2}/(*h*^{2}*v*^{2}) (a factor of 3 is for three polarizations).^{44} According to the IIM, $T1\u21922=v12/v22$ and $T2\u21921=1$, which gives $RC=(R2ball\u2212R1ball)/2$. The reason why the transmission from 1 to 2 is below one is exemplified in Fig. 7(b), showing the iso-energy surface of the phonon dispersion in *q*-space for a specific energy *ϵ*_{0}. With an ideal interface, transverse momentum/wavenumber (*q*_{⊥}) is conserved; such transitions are depicted by the horizontal arrows in Fig. 7(b). The maximum possible transverse wavenumber in 1 is larger than in 2, $q1,\u22a5max>q2,\u22a5max$. The incoming phonons from 1 can transfer into 2 if $q1,\u22a5\u2264q2,\u22a5max$ but are reflected if $q1,\u22a5>q2,\u22a5max$ since there are no available final states in 2. Contrarily, all incoming phonons from 2 can transfer into 1, since $q2,\u22a5<q1,\u22a5max$ always holds.

The phonon momentum vector is altered when traveling from one material to another. In the linear dispersion example, phonon propagation occurs normal to any point on the iso-energy surface of the *ϵ*(*q*) shown in Fig. 7(b). An incoming phonon with $q1,\u22a5=q2,\u22a5max$ will approach the interface with a moderate angle, but after crossing into material 2 it will exit traveling parallel to the interface. As mentioned above, any incoming phonon with $q1,\u22a5>q2,\u22a5max$ will be backscattered. Hence, we can define a critical angle above which phonons are reflected and below which phonons are refracted. Thus, this model is similar to the AMM, with the main difference being that the IIM gives a transmission of one up to the critical angle (which once averaged over all possible angles gives an average transmission less than one). The IIM should thus produce lower *R*_{C} than the AMM, as well as the DMM. For example, the DMM transmission coefficients are^{45} $T1\u21922=M2/(M1+M2)$ and $T2\u21921=M1/(M1+M2)$ (using the definition *M* in Ref. 40), which once inserted into Eq. (A12) gives $RC=(R2ball+R1ball)/2$ in the case of linear dispersions [larger than the IIM value of $RC=(R2ball\u2212R1ball)/2$]. We note that our approach is similar to the radiation limit^{46} but does not rely on a linear dispersion assumption, correctly predicts *R*_{C} = 0 for two identical materials, and works in the $T\u21921$ limit.

One approximation of this model is that it assumes that both materials have states which overlap in *q*_{⊥} [as considered in Fig. 7(b)]. At a given *ϵ*, the phonon states in both materials may be in different regions of the Brillouin zone, which could result in a transmission of zero if they do not overlap in *q*_{⊥}. However, the IIM would predict a finite transmission since it only considers the number of conducting channels on both sides. Finally, we note that the IIM does not discriminate among phonon polarizations; it focuses solely on the availability of final states for phonons.

### APPENDIX B: DFT MODELING DETAILS (BULK PHONONS)

The phonon dispersions and scattering rates for Si and Ge were computed using DFT. Structural relaxation was performed with the Quantum Espresso (QE) code^{47,48} using norm-conserving Perdew-Burke-Ernzerhof^{49} exchange-correlation pseudopotentials. The converged plane-wave energy cutoffs are 60 Ry and 70 Ry for Si and Ge, respectively. With a Brillouin zone sampling of 11 × 11 × 11, the converged lattice parameters are 5.574 Å and 5.772 Å, respectively, which are 2.6% and 2.0% larger than experiment. The second-order interatomic force constants were calculated as follows: With Si, density functional perturbation theory was employed, as implemented in QE, using a 4 × 4 × 4 *q*-grid. With Ge, the finite-displacement method was adopted, using Phonopy,^{50} with force calculations carried out on a 5 × 5 × 5 supercell and with a *k*-grid sampling the Γ-point. For anharmonic third-order force constant calculations, performed with ShengBTE,^{51} a 4 × 4 × 4 supercell was employed in both materials. Interactions of up to fourth (Si) and seventh (Ge) nearest neighbors were considered. The phonon dispersions and 3-phonon anharmonic scattering rates, used to calculate *M*_{Si,Ge}(*ϵ*), *λ*_{Si,Ge}(*ϵ*), and *g*_{Si,Ge}(*ϵ*), were computed on a *q*-grid of 60 × 60 × 60.

## REFERENCES

*M*(

*ϵ*), mean-free-path for backscattering,

*λ*(

*ϵ*), and average

*x*-projected velocity, $vx+(\u03f5)$, are defined as:

^{35,36}