Phonon transport across a Si–Ge interface: The role of inelastic bulk scattering

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),³ 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¹⁸ and transmission coefficients²⁵ 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.²⁸ 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.³⁸ 

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.³⁸ Temperature, in steady-state, is given by

where T(x,ϵ) is the temperature resolved at each phonon energy, T₀(x) is the temperature of a fictitious local equilibrium distribution (originating from the use of the relaxation time approximation), $D(ϵ)=λ(ϵ)vx+(ϵ)/2$ is the ϵ-resolved diffusion coefficient, $τQ(ϵ)=λ(ϵ)/(2vx+(ϵ))$ is the ϵ-resolved heat relaxation time, $CV(ϵ)=ϵg(ϵ)[∂nBE/∂T]Tref$ is the ϵ-resolved heat capacity, λ(ϵ) is the phonon mean-free-path for backscattering (average x-projected distance traveled before backscattering), $vx+(ϵ)$ 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+(ϵ)$ are defined in Ref. 40.

The expression for T₀(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³⁸ 

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²⁵ 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³⁹) 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→f(ϵ)$ 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(ϵ)=CV(ϵ)vx+(ϵ)/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.⁴¹ 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.⁴² 

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⁸ W/m², resulting in an effective R_C = 2.2 × 10⁻⁹ km²/W with the IIM (3.5 × 10⁻⁹ km²/W with the DMM). These values are close to other reported calculations using the DMM (4.6 × 10⁻⁹ km²/W),³¹ the AMM (≈5.5 × 10⁻⁹ km²/W),³¹ and MD (≈3.1 × 10⁻⁹ km²/W).⁵ 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 $∫0∞S(x,ϵ)dϵ=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) (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₀(x); heat is added when T₀(x) − T(x, ϵ) > 0 and removed when T(x, ϵ) − T₀(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).³⁴ 

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⁻¹ 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.