Phonon transport across a Si-Ge interface: the role of inelastic bulk scattering

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.


I. INTRODUCTION
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 comprised 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.
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 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, 13-16, 18, 19, 21, 24, 28, 29], this work builds on prior studies [23,27] to further elucidate the role of inelastic bulk scattering, a process that may be important in a broader class of heterojunctions.

II. SIMULATED STRUCTURE: SI-GE INTERFACE
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 phonon transport is induced by a ∆T = 1 K.

A. McKelvey-Shockley flux method
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 [30,31], it was subsequently adapted for phonon/thermal transport [32,33]. 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 [32][33][34][35][36]. Importantly, this method supports inelastic phonon scattering [36].
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 [36]. Temperature, in steady-state, is given by where T (x, ǫ) is the temperature resolved at each phonon energy, T 0 (x) the temperature of a fictitious local equilibrium distribution (originating from the use of the relaxation time approximation), D(ǫ) = λ(ǫ)v + x (ǫ)/2 the ǫ-resolved diffusion coefficient, τ Q (ǫ) = λ(ǫ)/(2v + x (ǫ)) the ǫ-resolved heat relaxation time, C V (ǫ) = ǫ g(ǫ) [∂n BE /∂T ] T ref the ǫ-resolved heat capacity, λ(ǫ) the phonon mean-free-path for backscattering (average x-projected distance traveled before backscattering), v + x (ǫ) the average x-projected phonon velocity, g(ǫ) the phonon density of states, n BE (ǫ, T ref ) the Bose-Einstein distribution, and T ref a reference temperature (300 K in this work). λ(ǫ) and v + x (ǫ) are defined in [37]. The expression for T 0 (x) is obtained by imposing conservation of energy (see Ref. [36]).
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 [32][33][34][35], 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 [36] T This approach was shown to compare well to the more rigorous and computationally demanding BTE, particularly when solved using a full phonon dispersion [32][33][34][35][36].

B. Ideal interface model
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 [24], 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 T i→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 channel 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 [37] for the definition of M (ǫ). The transmission coefficients calculated from Eq. (8) are used in the following expression the for ǫ-resolved thermal contact resistance (see Appendix A): where R ball (ǫ) = C V (ǫ)v + x (ǫ)/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 Eqns. (8)- (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. Note that an effective thermal contact resistance is extracted using R C = ∆T (l)/I Q (l) from the ǫ-integrated, effective temperature and heat current.

C. DFT modeling of 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 [38,39] using norm-conserving Perdew-Burke-Ernzerhof [40] exchangecorrelation 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. Second-order interatomic force constants were obtained using density functional perturbation theory (DFPT) as implemented in QE for Si with a 4×4×4 q-grid and Phonopy [41] for Ge with a 5×5×5 supercell and Γ point sampling. For anharmonic third-order force constant calculations, performed with ShengBTE [42], 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.

IV. RESULTS
A. Bulk phonon properties of Si and Ge Fig. 2a shows the calculated phonon dispersions, ǫ(q), for Si and Ge, which agree well with previous calculations [43]. 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. 2b. At low energies the M Si,Ge (ǫ) scale quadratically, as expected for linear dispersions. Fig. 2c shows the energy-dependent thermal contact resistance calculated with Eqns. (9)-(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. 2d. λ(ǫ) 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. 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 [44]. These properties serve as the input for our Si-Ge interface calculations. This length is interesting since some high-energy phonons are diffusive (λ(ǫ)≪L), while simultaneously some other low-energy phonons are ballistic (λ(ǫ)≫L). Fig. 3a shows the temperature of the Si-Ge junction versus position, x, and phonon energy, ǫ. Above 34 meV, 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 lowenergy acoustic phonons, where R C (ǫ) is high (see Fig. 2c). 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. [32][33][34][35]). The high-energy Si phonons have a temperature at the x = −L/2 boundary approching T L , indicating diffusive phonon scattering. Fig. 3b 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. 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. This value is close to other predicted values using the DMM (4.6 × 10 −9 Km 2 /W) [29] and MD (∼3.1 × 10 −9 Km 2 /W) [5], but lower than the AMM (1.7 × 10 −8 Km 2 /W) [29]. 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. Fig. 4a 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. [27]. At the interface, there are select energies, near 13 meV and 26 meV, that 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. Fig. 4b 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 (this would be 8% in bulk Si, shown as 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 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 Eqns. (1),(3), we obtain: Fig. 5a 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 strong redistribution of heat, i.e. inelastic scattering. On the Si side, heat is removed from the high-energy phonons and tranferred to midenergy 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 the heat across the interface from Si to Ge, as shown in Fig. 4a, this creates a highly nonequilibrium 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 redistribuion 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 [32,33], 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. 5a, in Fig. 5b we present the energyresolved 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.  5b shows which phonon energies will be heated/cooled, which correlates with the pattern in Fig. 5a. The low-energy acoustic modes are far from the equilibrium temperature, because they 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.  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 [32]. 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 intermediate length, arises from inelastic scattering, since assuming only elastic scattering (using the same approach) results in a linear T (x) [32].
The heat current versus L is plotted in Fig. 6b, 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. Fig. 6c 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 currrent is mostly uniform in x (away from the interface and contacts).

V. CONCLUSIONS
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 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 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.
Lastly, while our treatment is not comprehensive, for example not including 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.
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. [32] for more details on these quantities). We assume small temperature variations, Next, a condition that no net heat current flows in equilibrium is imposed on Eq. (A1). Using Eq. (A1), we find where I ± Q → I ± Q,ref in equilibrium. Using I ± Q,ref = ǫM n BE /h, leads to the following condition relating T 1→2 and T 2→1 : 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)-(A4) indicate that only T 1→2 must be specified, and that T 2→1 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 I − Q (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 I − Q (l + ) and I + Q (l − ) 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 Eqns. (A7)-(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. Lastly, Eq. (A12) is applicable under non-equilibrium phonon conditions. Next we specify the particular choice in T i→f 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 This model retrives two known limits: 1) Two identical materials. In the case of two identical materials, the IIM states that T 1→2 = T 2→1 = 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 R C → R ball 1 /2 when M 1 ≪ M 2 , a well-known result from the Landauer approach for phonon [32] and electron transport [45].
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 the phonon wavenumber (see Fig. 7a). The distribution of modes for a linear dispersion is M (ǫ) = 3πǫ 2 /(h 2 v 2 ) (factor of 3 is for three polarizations) [46]. According to the IIM, T 1→2 = v 2 1 /v 2 2 and T 2→1 = 1, which gives R C = (R ball 2 − R ball 1 )/2. The reason why the transmission from 1 to 2 is below one is exemplified in Fig. 7b, 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 horizontal arrows in Fig. 7b. The maximum possible transverse wavenumber in 1 is larger than in 2, q max 1,⊥ > q max 2,⊥ . The incoming phonons from 1 can transfer into 2 if q 1,⊥ ≤ q max 2,⊥ , but are reflected if q 1,⊥ > q max 2,⊥ since there are no available final states in 2. Oppositely, all incoming phonons from 2 can transfer into 1, since q 2,⊥ < q max 1,⊥ 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. 7b. An incoming phonon with q 1,⊥ = q max 2,⊥ , will approach the interface with a moderate angle, but after crossing into material 2 will exit traveling parallel to the interface. As mentioned above, any incoming phonon with q 1,⊥ > q max 2,⊥ 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 [47] T 1→2 = M 2 /(M 1 + M 2 ) and T 2→1 = M 1 /(M 1 + M 2 ) (using the definition M in [37]), which once inserted into Eq. (A12) gives R C = (R ball 2 + R ball 1 )/2 in the case of linear dispersions (larger than the IIM value of R C = (R ball 2 − R ball 1 )/2). We note that our approach is similar to the radiation limit [48], but does not rely on a linear dispersion assumption, correctly predicts R C = 0 for two identical materials, and works in the T → 1 limit.
One approximation of this model is that it assumes both materials have states that overlap in q ⊥ (as considered in Fig. 7b). 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. Lastly, we note that the IIM does not discriminate among phonon polarizations; it focuses solely on the availability of final states for phonons.