Quasi-harmonic theory for phonon thermal boundary conductance at high temperatures

Heterogeneous material interfaces are ubiquitous in novel material composites that drive current and future technologies from nano- to macroscales. Exemplary applications include silicon contacts on the nanometer scale used in logic,^1,2 phase change material interfaces used in non-volatile memory and photonic circuits,^3–5 wide- and ultrawide-bandgap material junctions used in amplifier and power devices,^6,7 and polycrystalline and multiphase composites used as environmental and thermal barrier coatings in engines or leading edges during high Mach number flight,^8–13 to name a few. As characteristic length scales of materials decrease, resulting in increases in densities of phases and heterogeneous interfaces, the temperature drops across these interfaces become the sources of the dominant thermal resistances in the composite material. This thermal boundary conductance then becomes the critical thermophysical property needed to understand and mitigate heat flow to ensure optimum functionality of functional materials and technologies from nano to macroscales.^14–16 

The phonon thermal boundary conductance across heterogeneous interfaces has been of interest to a range of communities dating back to Kapitza’s original work in 1941 on the temperature drop across solid/liquid helium interfaces at cryogenic temperatures.¹⁷ Several reviews over the past 50 years have discussed both theoretical and experimental advances that have furthered our understanding of this phonon interfacial process, including the importance and necessity in understanding heat flow across heterogeneous interfaces at non-cryogenic temperatures, i.e., room temperature and above.^{14–16,18–21} However, the theoretical understanding of phonon thermal boundary conductance is relatively limited at elevated temperatures, where materials are in the classical limit in terms of the thermal statistics of the phonons’ populations as compared to that at cryogenic and near-room temperatures.

This under-developed understanding of thermal boundary conductance in the classical limit is partially driven by the lack of experimental measurements at these temperatures. Classical molecular dynamics (MD) simulations of thermal boundary conductance across interfaces have provided invaluable insight into this phononic process in the classical limit; however, as we discuss in detail in Sec. II, the temperature dependence of thermal boundary conductance in several of these MD studies, which have been ascribed to anharmonic processes, are not described by standard theories that predict phonon thermal boundary conductance. Some of the most widely applied theoretical models for phonon thermal boundary conductance are rooted in harmonic thermodynamic theory, which immediately presents a discrepancy between models and molecular dynamics simulations.

Here, we derive a theoretical model for phonon thermal boundary conductance across solid interfaces in the classical limit using quasi-harmonic thermodynamics,²² an approach that accounts for phonon energy density changes via thermal expansion. This theory provides insight into the temperature dependence of the phonon thermal boundary conductance at high temperatures by accounting for anharmonic effects on the thermal boundary conductance via an increased internal energy of the material through an additional quasi-harmonic term that includes the material’s Grüneisen parameter. This anharmonicity-driven phonon energy density change cannot be predicted by traditionally employed harmonic thermodynamic-based theories, highlighting the importance of evoking quasi-harmonic thermodynamics formalisms to understand phonon transport across interfaces.

In the Dulong–Petit limit, at temperatures around and above the Debye temperature (⁠$ΘD$⁠), the volumetric heat capacity of a harmonic solid at constant volume is given by $CCL,H=3nkB$⁠, where $n$ is the atomic density and $kB$ is Boltzmann’s constant. Here, we use the subscript “H” to denote harmonic formalisms. In this high temperature classical limit (denoted by “CL”), however, volumetric thermal expansion driven by crystal anharmonicity leads to a temperature dependence that is not captured by this purely harmonic Dulong–Petit expression. This additional heat capacity, $CQH$⁠, driven by the difference between the heat capacity at constant pressure and that at constant volume, is given by $CCL,QH=9Bα2T$ in the classical limit, where $B$ is the bulk modulus, $α$ is the coefficient of thermal expansion, and $T$ is the absolute temperature. Thus, in the quasi-harmonic approximation (denoted by the subscript “QH”), the volumetric heat capacity of a solid in the high temperature classical limit is

Note, our use of the subscript “CL” is to indicate that we assume classical limit assumptions for the population distributions. Eq. (1) is readily derived from quasi-harmonic thermodynamic relationships²² and strictly valid in the classical limit. The thermal expansion coefficient can be further defined in terms of the material’s Grüneisen parameter, $γ$⁠, via the quasi-harmonic thermodynamics relationship $γ=3Bα/CH=mv2α/kB$ so that Eq. (1) becomes

where $m$ is the mass of the atom so that $mn=ρ$⁠, with $ρ$ being the mass density of the solid. These quasi-harmonic expressions for heat capacity are readily used to predict the volumetric heat capacity of solids but have not been invoked in theoretical semi-classical treatments of thermal boundary conductance, $h$⁠.

The goal of this work is to develop an analytical model that can provide insight into why the thermal boundary conductance between two solids increases with increasing temperature in the classical limit (e.g., when the phonon spectrum can be described by classical statistical mechanics and the heat capacity derived under harmonic assumptions is $CH=3nkB$⁠). Note, in this work, we approximate the classical limit as the temperature regime in which the temperature is greater than the Debye temperature; as shown in Fig. 1, when $T>ΘD$⁠, there is less than 8% error in approximating the temperature derivative Bose–Einstein distribution for the phonon population with that of the classical distribution, $kB/(ℏω)$⁠, with this difference decreasing with increasing temperature.

The increasing trend in thermal boundary conductance with increasing temperature has been well documented from a number of classical molecular dynamics (MD) studies,^23–35 where the nature of these simulations invoke classical dynamics to determine heat flow across interfaces. In this classical limit, the harmonic heat capacity, $CH$⁠, is a constant value as previously mentioned. Traditional and commonly employed models that are used to describe the thermal boundary conductance across solid interfaces (e.g., various models anchored in acoustic or diffuse mismatch theory)^19,36–49 are rooted in the assumption that the thermal boundary conductance is proportional to $CH$⁠. Thus, regardless of the inputs and assumptions used in calculating the thermal boundary conductance using any variations of these mismatch models, this approach fails to capture the classical limit temperature trend observed in classical MD and, therefore, does not properly describe the underlying phonon physics that are observed in these aforementioned classical MD results.^23–35 Recent computational works have used anharmonic nonequilibrium Green’s function (NEGF) formalisms to show that an increase in thermal boundary conductance driven by anharmonicity can also occur at temperatures below the Debye temperature,^50,51 which again is a trend not captured in harmonic NEGF simulations.^50–55 

An explanation for an unexpected increase in thermal boundary conductance with an increase in temperature that cannot be explained by these aforementioned harmonic-based mismatch theories is commonly attributed to multiple phonon processes (i.e., inelastic scattering) that is driven by anharmonicity at or near the interface. Experimental measurements suggesting these processes were first reviewed by Swartz and Pohl¹⁹ and subsequently been the focus of several works that have evoked this inelastic scattering postulate to explain experimental data.^56–61 However, it is important to note that none of these experimental measurements of thermal boundary conductance were taken at temperatures above both materials’ Debye temperatures. Thus, the mechanisms that are contributing to the disagreement between the experimental measurements and the harmonic-based mismatch theories may not necessarily translate to the phononic processes driving the increase in thermal boundary conductance with temperature that are observed in the MD and anharmonic NEGF simulations.

For example, in many of our prior works, we have implemented energy balances among multiple phonon states in both materials comprising the interface to account for multiple phonon inelastic scattering that can explain the aforementioned experimental results.^62–64 However, these models are rooted in harmonic theory. Thus, when both materials are at temperatures in the classical limit, these harmonic theories predict a temperature-independent thermal boundary conductance, in disagreement with observed molecular dynamic trends. Indeed, even when assuming maximum phonon transmission and inelastic radiation limits (i.e., allowing all inelastic processes to contribute to thermal boundary conductance),^59,65–67 harmonic-based models always predict a temperature-independent thermal boundary conductance, highlighting a void between traditional, and commonly employed mismatch models and a wide body of molecular dynamics works.

Motived by this, in the remainder of this work, we derive an analytical expression for phonon transmission and thermal boundary conductance in the classical limit, $ζCL$ and $hCL$⁠, respectively, using quasi-harmonic thermodynamics that accounts for the temperature dependence in $CCL$ in the classical limit via $CQH$⁠. Our theory suggests that the temperature dependence in $hCL$ is due in part to thermal expansion in the materials that change the phonon energy density as the temperature is increased. Our model accounts for these phononic changes through an additional quasi-harmonic term in $ζCL$ and $hCL$ with the materials’ Grüneisen parameters, $γ$⁠. In addition to our model predicting an increase in $hCL$ with temperature beyond the harmonic predictions (⁠$hCL,H$⁠), thus providing an explanation and model to predict anharmonic effects on thermal boundary conductance, this approach is also consistent with prior MD and anharmonic NEGF simulations that suggest that anharmonic effects on thermal boundary conductance is driven by phononic processes and anharmonicity in materials near the interface and not at the interface.^30,50,55,68

Indeed, over a decade ago, we first theorized this⁶⁹ by rederiving the diffuse mismatch model and accounting for intrinsic three-phonon scattering events near the interface, successfully predicting the classical MD-simulated thermal boundary conductance across a Stillinger–Weber silicon/germanium interface²³ and its increase in thermal boundary conductance vs temperature in the classical limit. This model, however, is still rooted in harmonic thermodynamics (i.e., the predicted heat capacities of the materials are constant in the classical limit, which is inconsistent with MD simulations)⁷⁰ but accounted for inelastic scattering through a three-phonon scattering rate that was included in the calculation of $hCL$⁠. Further, the model requires knowledge of three-phonon scattering rates in each material comprising the interface, which is not as widely accessible as bulk thermodynamic properties that can be inputs to the model derived in this work.

It is important to note that to date, there have been no experimental measurements of $hCL$ across solid interfaces in the classical limit, as previously mentioned. However, given the classical nature of molecular dynamics (MD) simulations, we benchmark out theory against thermal boundary conductance determined across a Lennard-Jones heavy Ar/Ar interface that we simulate with MD in this work. We show these comparisons of our model to MD simulations in Sec. IV after we derive our quasi-harmonic thermodynamic-based model for $hCL$ in Sec. III. Additionally, we compare our model to recent anharmonic NEGF simulations of thermal boundary conductance across Si/Ge interfaces⁵⁰ by extending our model to lower temperatures below the Debye temperature, which demonstrates the utility and potential of our quasi-harmonic thermal boundary conductance formalism for capturing anharmonic phonon effects on thermal boundary conductance even at lower temperatures.

Along these lines, while the derivation of our model in the classical limit thus allows for simplified analytical expressions to calculate $h$ including quasi-harmonic effects in the classical limit, the formalism derived herein can also be implemented at lower temperatures by accounting for quantum effects in the phonon distribution function. We thus present a more generalized form of our model to account for lower temperature phonon distributions and more generic phonon dispersions in the  Appendix.

As a final note before deriving our model, it is important to differentiate the anharmonicity that is captured in our quasi-harmonic model from that which is accounted for in three-phonon scattering processes in materials. Our model derived in this work accounts for the anharmonicity-driven energy density changes to the phononic spectra in the materials comprising the interface. This manifests itself as an increased temperature dependence in the phononic heat flux, which also impacts the phonon transmission coefficient. In our phonon transmission coefficient calculations, we assume all phonon modes in each material can exchange energy upon diffusive interface scattering via diffuse mismatch theory, an analytical formalism originally proposed by Dames and Chen⁶⁷ (note, a similar processes involving specular interface scattering can also be captured via acoustic mismatch theory as presented by Prasher and Phelan).³⁹ Thus, our quasi-harmonic model accounts for both energy density changes and inelastic phonon conversion driven by thermal expansion on thermal boundary conductance across interfaces via an energy balance. While decoupling these two processes is not our focus here, it is important to note that the inelastic scattering processes that we account for in the phonon energy conversion across an interface are different than the intrinsic three-phonon processes (e.g., Umklapp processes) that drive phonon–phonon scattering in a homogeneous material. Inelastic effects on the thermal conductance between two materials are driven by the phonon energy densities in both materials adjacent to the interface, where intrinsic three-phonon scattering processes are occurring only in one material. Or, stated differently in the context of our model, inelastic effects on phonon conversion across interfaces will involve the Grüneisen parameters of both materials at the interface, where intrinsic three-phonon effects will only account for anharmonicity and the Grüneisen parameter in a single material. Our current formalism is derived to capture the former process, where our prior formalism only accounted for the later.⁶⁹ Depending on the control volume assumed to analyze the thermal boundary conductance across interfaces (e.g., how far away from the interface one considers when defining the interfacial region temperature drop), both processes could be playing a role in the calculated or measured thermal boundary conductance.

We root our theoretical derivation and calculations of $hCL$ in the quasi-harmonic approximation and semi-classical formalisms, and for analytical simplicity, we assume an isotropic Brillouin zone. This allows for useful closed form expressions of $hCL$ in the classical limit permitting useful relationships between phonon velocities, Grüneisen parameters and thermal boundary conductances across material interfaces. For general applicability, we show the derivation and equations without these isotropy assumptions in the  Appendix, but all of the calculations in this work are based on the isotropic Brillouin zone assumption.

This isotropic Brillouin zone approach is common in prior derivations of $CH$ and $hH$ in the harmonic limit. For example, isotropy permits us to assume a spherical Brillouin zone so that the harmonic heat capacity is given by

where $N$ is the Bose–Einstein distribution, which in the classical limit, becomes $N=kBT/(ℏω)$⁠, so that

When assuming three modes and a maximum phonon wavevector of $kmax=2π/a′$⁠, where $a′$ is an effective lattice constant that accounts for the correct volume of the Brillouin zone, Eq. (4) reduces to $CCL,H=4πkB/a′3$⁠. The correct choice of effective lattice constant, $a′$⁠, is necessary to ensure the correct number of modes depending on the assumption of the Brillouin zone symmetry. For a spherical Brillouin zone, $a′=a0(4π/(3b))(1/3)$⁠, $a0$ is the equilibrium lattice constant and $b$ is the number of atoms in the unit cell, so that $n=b/a03$⁠. This leads to $CCL,H=3nkB$⁠.

We show calculations of Eqs. (3) and (4) in Fig. 2 for silicon assuming six modes (three acoustic and three optical) and a maximum wavevector of $kmax=2π/a0$⁠, where $a0=5.43×10−10$ m for silicon.⁷² For input into Eq. (3), we determine $ωk$ by fitting a fourth order polynomial⁴⁹ to the acoustic and optical modes along the $Γ→X$ direction of the Si phonon dispersion calculated by Weber.⁷³ We make no adjustment to the lattice parameter to account for inaccuracies in our Brillouin zone volume based on our isotropic assumption due to the fact that we are using realistic dispersion relations along a specific Brillouin zone direction; this results in an error in $CCL,H$ calculation of $(CCL,H/(3nkB))−1=π/3−1∼5%$ for a diamond cubic crystal structure. We also include measured reference data in Fig. 2 for the volumetric heat capacity of Si.⁷¹ Consistent with our discussion in reference to Fig. 1, the discrepancy between $CH$ (which agrees well with the experimental data at temperatures up to the Debye temperature, cf., Fig. 2 inset) and $CCL,H$ is $∼8%$⁠.

In this vein, the phonon heat flux in a material that is incident and perpendicular to some plane in $z$ (such as an interface) is $q=∑j∑kx,j=−∞∞∑ky,j=−∞∞∑kz,j=0∞Uk,jvkz,j$⁠, where $Uj$ is the internal energy of phonon polarization $j$⁠, $vkz,j$ is the phonon group velocity in the direction perpendicular to the interface, and $dk=dkxdkydkz$⁠. In the classical limit under the harmonic Dulong–Petit approximation and assumptions of Brillouin zone isotropy, $qH$ is given by

where $vk=vkz/|cos⁡(θ)|$ is the phonon group velocity.⁷⁴ Note, assuming an isotropic Debye solid, Eq. (5) becomes $qCL,H=3nkBTv/4$⁠. Recognizing the familiar Landauer relationship⁷⁵ that $h=∫kζk∂qk/∂Tdk$ results in the common expression for the thermal boundary conductance based on the transmission of phonon energy in the classical limit,

We discuss the transmission coefficient, $ζ$ in more detail below, but first we turn our attention to the fact that this expression is derived under the typically assumed harmonic assumptions so that, assuming a Debye solid and following standard derivations in the Dulong–Petit limit, $hH=3nkBvζCL,H/4=CHvζH/4$⁠.⁷⁶ Clearly, the assumption of a temperature-independent heat capacity in the classical limit manifests itself into predictions of $hCL$ that are also temperature independent, which is inconsistent with the aforementioned several works that report on a temperature dependent thermal boundary conductance using classical molecular dynamics simulations.^23–35 

The first consideration in accounting for anharmonicity in Eq. (6) is to properly account for changes in the heat flux, $q$⁠. Thus, we must first consider both the harmonic contribution to phonon internal energy, $UH$⁠, the temperature derivative of which results in $CH$⁠, and the anharmonic portion of the internal energy, $UQH$⁠. The quasi-harmonic thermodynamic approximation for $UQH$ accounts for changes in mode energy driven from thermal expansion effects through the Grüneisen parameter, so that

where $γk,j$ is the wavevector dependent Grüneisen parameter of phonon polarization $j$⁠. In the limit of a Debye system with three identical modes and a spectrally independent Grüneisen parameter, Eq. (7) simplifies to $UCL=3nkBT+9nkB2γ2T2/(2mv2)$ so that $CCL=∂U/∂T$ is equivalent to Eq. (2).

Note, the quantity in the denominator of Eq. (7) is related to the bulk modulus, $B$⁠, only defined on a modal basis with high frequency dispersive phonon modes. Strictly speaking, the bulk modulus is defined as a continuum mechanical property that we define from the continuum relationship $B=ρv2$⁠, so we caution that our definition of this spectrally dependent $B$ only converges with the measured elastic property when assuming Debye modes. For example, if we assume different longitudinal and transverse modes in this Debye continuum limit, the denominator of Eq. (7) becomes $B=Mn(vL2+2vT2)/3$ and predicts $B=99.5$ GPa for the bulk modulus of silicon, assuming $M=4.66×10−26$ kg, longitudinal and transverse phonon velocities as constant values of $vL=8190$ m s$−1$ and $vT=5510$ m s$−1$⁠, respectively, and a maximum wavevector of $kmax=(2π/a0)33b/(4π)$⁠, where $b=8$ atoms and $a0$ is the lattice constant for Si.⁷³ This predicted value is in close agreement with the accepted measured bulk modulus for silicon (⁠$B=100$ GPa).⁷⁷ 

In Fig. 2, we show calculations of $C=CH+CQH$⁠, which is the temperature derivative of Eq. (7) assuming $N=(exp⁡[ℏω/(kBT)]−1)−1$⁠. Calculations of $C$ when including $CQH$ capture the correct heat capacity temperature trend for silicon in the classical limit due to the inclusion of quasi-harmonic thermodynamics in its derivation. Note, there is a $<5%$ overestimate in $CCL$ compared to the calculated $C$ without using classical statistics at $T=ΘD$⁠, with this difference decreasing with increasing temperature into the classical limit.

The spectrally dependent expression of the internal energy defined in Eq. (7) allows us to now include the role of phonon energy density changes due to thermal expansion on the phonon heat flux via

From Eq. (8), we can now define a thermal boundary conductance in the classical limit that accounts for anharmonicity through thermal expansion via the quasi-harmonic treatment given by

where now we define $ζCL$⁠, a transmission coefficient in the classical limit that accounts for both the harmonic flux and the flux accounting for anharmonicity via the quasi-harmonic treatment. We note that we assumed a spectrally independent phonon transmission coefficient in Eq. (9), which facilitates calculations of $ζCL$ under the assumptions of detailed balance among harmonic and anharmonic fluxes, described below. While the role of crystal anharmonicity on the spectrally dependent phonon transmission should be investigated (specifically deconvolving the harmonic and anharmonic contributions for each mode), this is beyond the scope of this first-time quasi-harmonic theoretical derivation of phonon thermal boundary conductance.

Analytically determining the phonon transmission coefficient across interfaces has been a topic of many studies,^{39–43,45–49,62–64,78,79} and requires evoking assumptions of how modes interact with other modes across interfaces. At the elevated temperatures when materials are in the classical limit, the dominance of relatively shorter wavelength, higher frequency phonons, and highly dispersive, low group velocity optical modes supports the assumption of diffusive phonon scattering as the process driving thermal boundary conductance across material interfaces.⁴⁹ For this work, we thus derive (⁠$ζCL$⁠) based on the widely applied diffusive mismatch model (DMM).¹⁹ In this approach, the definition of diffusive scattering coupled with the evocation of detailed balance, which is rooted in the assumption of black boundaries and subsequent complete phonon thermalization, leads to a transmission coefficient taking the form of $ζ=q2/(q1+q2)$⁠. For example, assuming harmonic thermodynamics, the elastic phonon transmission coefficient becomes

where the subscripts “1” and “2” refer to materials 1 and 2. Following this, the thermal boundary conductance predicted via the DMM assuming only elastic scattering is given by

For the purpose of this work, we assume this elastic transmission coefficient is frequency independent up to the maximum phonon frequency of the material with the lower maximum phonon frequency,⁸⁰ an assumption that we have rigorously vetted previously against other spectral assumptions in calculating $ζCL,H$⁠.⁸¹ Thus, we only consider wavevectors of phonons with energies that are elastically accessible when calculating $qCL,H$ and $ζCL,H$ in this work.

Evoking the same assumptions while including both the harmonic and quasi-harmonic fluxes yields the phonon transmission coefficient that accounts for anharmonicity via volume expansion, given as

which includes contributions from both harmonic and quasi-harmonic fluxes in the classical limit, as defined in Eq. (8), consistent with the application of detailed balance used to derive the interfacial phonon transmission coefficient. From this, the thermal boundary conductance in the classical limit that now accounts for anharmonicity via volume expansion through quasi-harmonic thermodynamics is given by Eq. (9) with Eq. (12). It is important to emphasize that when calculating $qH$ in $ζCL$ and $hCL$⁠, only wavevectors of phonons with energies that are elastically accessible should be considered. However, for the $qQH$ contributions, all phonon modes should be considered and integrated over when calculating this term in $ζCL$ and $hCL$⁠. Thus, as discussed in Sec. II, our quasi-harmonic model accounts for both energy density changes and inelastic phonon conversion driven by thermal expansion on thermal boundary conductance across interfaces.

To validate our quasi-harmonic DMM, we first turn to Guo et al.’s⁵⁰ anharmonic NEGF simulations reporting on the thermal boundary conductance across Si/Ge interfaces as a function of temperature, as plotted in Fig. 3 (square and circle symbols). The advantage of comparing our model to these simulations is that the NEGF formalism is rooted in Landauer theory so that the fundamental assumptions used for calculating thermal boundary conductance based on phonon transmission are similar to those evoked in our DMM approach. Similarly, the force constants used to simulate Si and Ge in Guo et al.’s work were derived via ab initio approaches⁸² so that we can use experimentally measured parameters on Si and Ge for input to our model, as opposed to modified parameters that would be more representative of empirical potentials. The disadvantage, however, is that the anharmonic NEGF simulations⁵⁰ were conducted in a regime that was not in the classical limit (i.e., the simulations were conducted at temperatures below the Debye temperature), and, thus, we modify Eqs. (8) and (9) to account for a Bose–Einstein distribution at temperatures below the classical limit (see  Appendix). Further, the force constants and lattice spacing in their simulated Si and Ge systems were assumed identical (albeit derived via ab initio approaches), so the validity of using experimentally derived input to our model still could introduce discrepancies.

We plot model calculations for the thermal boundary conductance assuming only harmonic thermodynamics (“Harmonic model”) and accounting for anharmonicity through our quasi-harmonic thermodynamic derivation (“Quasi-harmonic model”). For input into our models, we determine $vk,j$ by fitting a fourth order polynomial⁴⁹ to the acoustic and optical modes along the $Γ→X$ direction of the Si and Ge phonon dispersion calculated by Weber.⁷³ Atomic mass and number density are taken from references.^72,77 The last unknown in our model is the Grüneisen parameter for Si and Ge, which dictates the relative contribution of the anharmonic heat flux in both the phonon transmission and thermal boundary conductance calculations in the quasi-harmonic model shown in Fig. 3. The Grüneisen parameter of materials is a mode and frequency dependent quantity, and, as we address below, rigorous input into our model should account for this spectral variation in $γk,j$⁠.^83–93 However, for this comparison to the anharmonic NEGF results from Guo et al.,⁵⁰ we assume a gray Grüneisen parameter in the calculations of the quasi-harmonic model (i.e., $γ$ is independent of mode and frequency). We assume $γ=0.6$⁠, which we determine from our fit of our quasi-harmonic heat capacity model to the experimental data, shown in Fig. 2. This gray Grüneisen parameter may be an acceptable assumption for Si and Ge since prior works have shown that $γk,j$ is relatively constant as a function of wavevector for the longitudinal modes in these materials.^83,84,86,87 However, the transverse modes’ Grüneisen parameters can exhibit more spectral dependencies than those of the longitudinal modes in Si and Ge^83,84,86,87 and so thus this introduces a source of assumption-driven uncertainty in our model.

Our calculations for thermal boundary conductance across the Si/Ge interface using our quasi-harmonic model qualitatively capture the increase in thermal boundary conductance with an increase in temperature approaching the Debye temperature of Si found by Guo et al.’s⁵⁰ in their anharmonic NEGF simulations. This temperature trend is not captured by thermal boundary conductance predictions using harmonic thermodynamic-based models, also shown in Fig. 3 and further exemplified in Guo et al.’s⁵⁰ harmonic NEGF simulations. This comparison demonstrates that our quasi-harmonic model can, at least qualitatively, capture of the role of anharmonicity on thermal boundary conductance across interfaces. Further, given our model is based on changes in material properties adjacent to the interface, and cannot capture interfacial phonon processes at the interface, our model supports the conclusions from Guo et al.⁵⁰ NEGF simulations and other prior MD works^30,55,68 that anharmonicity in thermal boundary conductance can be captured by phononic processes near the interface (i.e., phonon energy density changes in the material adjacent to the interface) as opposed to processes localized at the interface.

We note here that while our model captures the qualitative trends with the Si/Ge NEGF simulations from Guo et al.’s work,⁵⁰ there have been other MD or NEGF simulations of thermal boundary conductance across Si/Ge interfaces in the literature that report different values than those from reported by Guo et al.⁵⁰ For example, prior anharmonic NEGF simulations from Dai and Tian⁵¹ predict a qualitatively similar trend to those of Guo et al.⁵⁰ and to our model predictions in this work, however, the magnitude of the anharmonic contribution in Dai and Tian’s⁵¹ Si/Ge thermal boundary conductance is $∼3$ times larger. The differing results between Guo et al.’s⁵⁰ and Dai and Tian’s⁵¹ reports could be due to procedural differences in their simulation domain, NEGF implementation, or ab initio inputs, but this is speculation as an analysis of the reasons for these differences between two works in the literature from other groups is not in the scope of this current manuscript. Additionally, we currently do not understand why our model would agree better with Guo et al.’s⁵⁰ calculations as opposed to Dai and Tian’s⁵¹ unless the use of our literature values for the Grüneisen parameter for Si and Ge are better captured in Guo et al.’s⁵⁰ calculations. We address this source of uncertainty in more detail below.

As another example, Landry and McGaughey²³ used MD simulations to determine the thermal boundary conductance across interfaces between Stillinger–Weber Si and Ge, and found values that are not only higher than those reported by Guo et al.⁵⁰ but also showed a much steeper increase in $h$ as a function of temperature. In our previous work,⁶⁹ we developed a model that captured these MD predicted trends from Landry and McGaughey²³ by accounting for intrinsic three-phonon scattering events near the interface. Our model for these data accounted for phonon scattering at some finite distance away from the interface, which led to increases in the thermal boundary conductance. For this process to be taking place, the material control volume in which thermal boundary conductance is considered must be large enough to include phonon scattering events several unit cells away from the interface (depending on the phonon wavelength). Our current formalism in this work is not derived to capture these effects, as we only consider phonon transmission across the interface and not phonon scattering away from the interface. However, it is interesting to consider future work where both of these processes are accounted for in a thermal boundary conductance model. In fact, recent work by Thomas and Srivastava⁹⁴ accounted for both interface scattering and three-phonon anharmonic scattering in the thermal transport in Si/Ge superlattices using a semi-ab initio approach.⁹⁵ Their calculated thermal boundary conductances across the interfaces in the Si/Ge superlattices exhibit trends more indicative of purely harmonic interactions, and they observed no increase in $h$ with temperature as temperature increased approaching and exceeding the Ge and Si Debye temperatures. We posit two potential possibilities for the differences between these results by Thomas and Srivastava⁹⁴ and those from the other works discussed here. First, atomic potential and procedural differences between this semi-ab initio approach and the other aforementioned NEGF and MD studies will certainly play a role. Second, the work by Thomas and Srivastava⁹⁴ is simulating a Si/Ge superlattice, and the interface spacing in the structure could be small enough that interface scattering is dominating over intrinsic anharmonic scattering in the individual Si or Ge layers. Thus, if the aforementioned delocalized phonon scattering events in our previous model⁶⁹ were to play a role in thermal conductance across the Si/Ge interfaces, then due to the high density of interfaces, the scattering processes at the Si/Ge interface could be dominating over delocalized scattering processes near the interface, which would lead to additional differences between the superlattice modeled by Thomas and Srivastava⁹⁴ and the individual isolated interfaces that are modeled by NEGF and MD and discussed in the model in this work.

A glaring shortcoming in our model validation approach above is the assumption of a gray Grüneisen parameter for Si and Ge. Our model developed in Eq. (9) is derived assuming spectrally mode dependent $γk,j$⁠. Thus, without knowledge of the spectral and modal dependence of $γk,j$⁠, we cannot fully validate our theory. Additionally, we cannot rule out discrepancies between our model and the NEGF simulations based on the fact that Guo et al.⁵⁰ assumed identical force constants and lattice spacings in their simulated Si and Ge systems. This different model input from what is assumed in the NEGF simulations would impact the calculated phononic properties, including our assumption of $γ$⁠.^90–93 

Thus, to overcome this and show the general applicability of our quasi-harmonic theoretical model to predict $hCL$⁠, we compare our theoretical predictions to $hCL$ calculated from MD simulations for an interface comprised of two mass-mismatched Lennard-Jones (LJ) solids. All interatomic interactions in our MD simulations are described by the 6–12 LJ potential, $U(r)=4ε[(σ/r)12−(σ/r)6]$⁠, where $r$ is the interatomic separation, and $σ$ and $ε$ are the LJ length and energy parameters, respectively. The potential is parameterized for Ar,⁹⁶ where $ε=10.3$ meV and $σ=3.4$ Å. The phononic spectra mismatch resulting in an appreciable temperature drop (⁠$ΔT$⁠) at the interface between the LJ solids is introduced by setting the mass of one of the solids to be four times higher than that of LJ Ar. The size of the computational domain is set to $10a0×10a0×60a0$⁠, which ensures that the cross-sectional area and the length of the simulation domain have negligible influence on the calculated $hCL$⁠. We enforce periodic boundary conditions on the $x$- and $y$-directions, whereas fixed boundaries with four monolayers of atoms at each end are prescribed in the $z$-direction. The LJ-based structures are allowed to relax under the Nose–Hoover thermostat followed by an isothermal–isobaric ensemble for a total of 2 ns at 0 bar pressure. Following equilibration, we apply a heat flux across the computational domain in the $z$-direction and determine the $ΔT$ from the steady-state temperature profiles of the LJ solids calculated from the linear regression analysis. All MD simulations are performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package.⁹⁷ In order to calculate the phonon dispersions and the spectral Grüneisen parameters for our Ar and mass-heavy Ar solids, we utilize lattice dynamics (LD) calculations carried out with the General Utility Lattice Program (GULP).⁹⁸ We then use $ωk,j$ and $γk,j$ calculated along the $Γ→X$ direction from these LD simulations as direct input into our model for $hCL$⁠, thus providing a more direct assessment of our model. Given the classical nature of the MD simulations, we use Eq. (9) to calculate $hCL$ and $qCL$ in Eq. (8) as input to calculate $ζCL$ via Eq. (12).

Figure 4 shows the MD-simulated thermal boundary conductance for mass-heavy Ar/Ar interfaces compared to $hCL$ predicted from our model. Without any adjustable parameters, our model for $hCL$ captures the linear trend in thermal boundary conductance and agrees well with the absolute values from the MD simulations; both values and trends in the MD results cannot be captured with harmonic thermodynamic-based theory, exemplified by our calculations of $hCL,H$ [Eq. (11)]. While we certainly expect error in our predictions due to the assumption of a spherical/isotropic Brillouin zone, which is of the order of $<10%$ as addressed in Sec. III, we also expect additional disagreement between our model and the MD results based on our application and assumptions of the DMM transmission (which have been addressed in detail in prior works),^48,81 and our assumption of a temperature-independent Grüneisen parameter.

To understand potential inaccuracies due to our assumption of an isotropic Brillouin zone and the use of the dispersion along the $Γ→X$ direction, we repeat the calculations for $hCL$ and $hCL,H$ for the mass-heavy Ar/Ar interfaces using the dispersion along the $Γ→K$ direction to represent the phonon density of states in the spherical Brillouin zone. The non-degenerate transverse branches and $6%$ larger Brillouin zone dimension along this direction of high symmetry make this comparison a pertinent test of uncertainty due to our assumed phonon dispersion to represent the density of states in the isotropic Brillouin zone. We plot the calculations of $hCL$ normalized by $hCL,H$ using both $Γ→X$ and $Γ→K$ in Fig. 5 (lines labeled as “$X$” and “$K$⁠,” respectively) and compare these calculations to the data from Fig. 4 normalized to the lowest temperature data point. While the use of different phonon dispersions certainly changes the predicted magnitude of the quasi-harmonic contribution to the predicted thermal boundary conductance, the general temperature trends remain the same, and the agreement between the model and the MD simulations remains acceptable. Similarly, we test our assumption of using a temperature-independent Grüneisen parameter by recalculating our model calculations for $hCL$ normalized by $hCL,H$ using the phonon dispersion along the $Γ→X$ direction, but assuming $20%$ higher value for $γk,j$⁠, consistent with the increase in $γ$ that has been reported in solid Ar by Tilford and Swenson.⁹⁹ The results for this calculation are shown by the line labeled as “$γ$” in Fig. 5, and, similar to the result when we change the phonon dispersion, the general temperature trends remain the same, and the agreement between the model and the MD simulations remains acceptable. This further supports that our quasi-harmonic formalism derived in this work can properly account for the phonon dynamics driving the temperature dependent thermal boundary conductance in the classical limit.

Finally, the relatively small disagreement between $hCL,H$ and the MD results approaching 0 K in Fig. 4 could be due to assumptions related to the DMM formalism, and any shift in $hCL,H$ to better agree with the low temperature limit of the MD results would result in even better agreement between our $hCL$ calculations and the temperature dependent MD results. For example, recent work by Thomas and Srivastava⁹⁴ have studied thermal boundary conductance across interfaces in superlattices using a partially diffusive and partially specular scattering assumption. They found that specular scattering at interfaces as captured with an acoustic mismatch model (AMM)^36,38 formalism can lead to a predicted thermal boundary conductance across these interfaces that decreases as compared to the DMM predictions. In this light, the diffusive assumption imposed in our DMM approach to calculating $h$ could lead to an over prediction in our model as compared to the MD simulations in Fig. 4. Most likely, a similar partially diffusive–partially specular scattering assumption at interfaces coupled with non-local intrinsic phonon scattering and accurate phonon potentials and Brillouin zone dimensions, as detailed in the formalism outline for phonon transport in superlattices by Thomas and Srivastava,⁹⁴ would widen the general applicability of our quasi-harmonic thermal boundary conductance model, and should be the focus of future works.

This work presents a theoretical model for phonon thermal boundary conductance across solid interfaces in the high temperature classical limit using quasi-harmonic thermodynamics, an approach that naturally accounts for phonon energy density changes in a solid driven by thermal expansion. Our theory suggests that an increase in thermal boundary conductance with an increase in temperature above the Debye temperature is due to the volumetric energy density changes driven by anharmonic effects in the materials adjacent to the interface. Our model accounts for this through an additional quasi-harmonic term that includes the material’s Grüneisen parameter, $γ$⁠. Materials with larger values of $γ$ and thus larger thermal expansion coefficients will have increased contributions from anharmonic processes in heat flow across interfaces and thus increased thermal boundary conductance in the classical limit. In addition to our model predicting an increase in $hCL$ with temperature, thus providing an explanation and model to predict the aforementioned NEGF and MD studies, our results also support the results from prior MD and anharmonic NEGF studies that suggest that inelastic scattering effects on thermal boundary conductance are driven by phononic processes and anharmonicity in materials near the interface and not at the interface.^30,50,55,68

The utility in this model also comes in its simplicity. The role of anharmonicity in $hCL$ can be assessed through knowledge of bulk material properties of the materials comprising the interface. As a concluding example, assuming a Debye solid and a gray Grüneisen parameter, we can rewrite Eq. (9) as

Given that $γ=mv2α/kB$⁠, Eq. (13) provides a tractable and semi-quantitative method to predict the role of anharmonicity on the thermal boundary conductance across solid interfaces at high temperatures with only knowledge of phonon velocity, density, and thermal expansion coefficient. Thus, this model can help in screening materials for high interface density composites to increase thermal conductance and mitigate temperature in a range of applications.

We appreciate support from the National Science Foundation (NSF) (Grant No. 1921973) and the Office of Naval Research (Grant Nos. N00014-21-1-2622 and N00014-18-1-2429). P.E.H. acknowledges support from the Alexander von Humboldt Foundation while completing this manuscript.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

A more general form of the thermal boundary conductance that accounts for energy density changes via thermal expansion and applicable to any phonon dispersion and distribution is derived here. We begin with the phonon heat flux incident and perpendicular to the interface, given by

Generalizing the expression for phonon internal energy from Eqs. (7) and (A1) becomes

From this, the thermal boundary conductance becomes

where, following Eq. (12), the transmission coefficient $ζ$ is given by

where $q$ is defined in Eq. (A2).