Thermal boundary resistance in semiconductors by non-equilibrium thermodynamics

We critically address the problem of predicting the thermal boundary resistance at the interface between two semiconductors by atomistic simulations. After reviewing the available models, lattice dynamics calculations and molecular dynamics simulation protocols, we reformulate this problem in the language of non-equilibrium thermodynamics, providing an elegant, robust and valuable theoretical framework for the direct calculation of the thermal boundary resistance through molecular dynamics simulations. The foundation of the method, as well as its subtleties and the details of its actual implementation are presented. Finally, the Si/Ge interface showcase is discussed as the prototypical example of semiconductor heterojunction whose thermal properties are paramount in many front-edge nanotechnologies. Graphical Abstract


Introduction
The prediction of the thermal transport properties of a bulk semiconductor material (where the main heat carriers are phonons or, more generally, vibrational modes) usually proceeds through the calculation of its thermal resistivity R qq (or, equivalently, its inverse quantity κ qq , namely the thermal conductivity) [1,2]. While the most fundamental theory to accomplish this task is provided by the Boltzmann transport equation [3] (which is now numerically solvable either in the usual single-mode relaxation-time approximation [3,4] or even exactly [5][6][7]), molecular dynamics (MD) simulations [8] represent in fact the most popular tool for computing the thermal transport coefficients, due to their versatility, ease of implementation and comparatively small computational effort. The resistivity is typically evaluated by assuming the Fourier law ∂T/∂z = −R qq J q and by explicitly computing the heat flux J q and the temperature gradient ∂T/∂z in a steady-state condition of thermal transport (for further convenience we remark that for a sample with total thickness L z and cross section S the bulk thermal resistance is defined as R qq = R qq L z /S). This approach is largely exploited by non equilibrium molecular dynamics (NEMD) methods where the flux J q is imposed and the resulting temperature gradient ∂T/∂z is calculated [9,10], or viceversa [11,12]. Alternatively, the thermal diffusivitȳ κ qq = κ qq /ρC v (where ρ and C v are the system mass density and specific heat, respectively) is evaluated during the system approach to equilibrium (AEMD) in a microcanonical evolution from an initial configuration set by a non-uniform temperature profile [13]. Finally, equilibrium statistical mechanics offers a direct way to compute thermal conductivity through the general Green-Kubo theory of transport coefficients (equilibrium MD, EMD) [12,14]. Overall, this panoply of methods allows for computing, in a large variety of (in principle) equivalent ways, the thermal transport coefficients as function of temperature, atomic structure, defect-induced disorder, or chemical composition, thus providing a full characterization of materials thermal properties. It is also possible to include as well quantum features [15][16][17] if the system is simulated well below its Debye temperature. However, concerns have been raised about the actual reliability of the proposed methods [18]. The case of heat transport across an interface is intriguingly different and much more subtle [19][20][21][22].
Here a temperature gradient is applied across a heterojunction (HT) between two semiconductors which behave differently as for their heat transport properties. Because of this and of interface-specific properties, across the interface region a sudden temperature drop T occurs, giving rise to a thermal boundary resistance (TBR) R TBR qq which is phenomenologically described as R TBR qq = T/J q , assuming a steady-state condition. This phenomenon is nowadays widely referred to as Kapitza resistance [20][21][22], although this expression originally referred to the thermal resistance occurring at the interface when heat flows from a solid into a liquid. While casting the problem in this way could seem straightforward, actually there are a number of subtleties indeed requiring a deeper understanding, which basically defines the main goal of the present short review (incidentally, we warn against a widely common abuse of notation: while using the same term 'resistance' for both the bulk and interface case, the two quantities have in fact different units: K/W and Km 2 /W, respectively).
This issue is not just academic but rather a really key feature in nanoscience. In fact, interfaces are everywhere at the scale where most of the present-day semiconductor technology deploys: their role can hardly be underestimated since the fabrication methods of the electronics industry are continuously refined to produce ever smaller devices in the nm size range. In this context, therefore, atomically well-defined interfaces represent the most important structural feature affecting the transfer of energy through thermal exchange. So, the heat flux across interfaces plays a vital role in many front-end applications, including electronics (cooling of nanodevices and optimal control of their thermal budget) [23,24], information technology (phononics) [25,26] and energy harvesting/production (superlattices for thermoelectric conversion) [27][28][29][30][31][32]. In conclusion, exploring heat-related phenomena occurring at semiconductor interfaces will facilitate the precise control of their properties and, ultimately, their best-tailored processing for applications in the above technologies.
While the above scenario stands for the need of a detailed understanding of interface thermal properties (possibly including thermal rectification phenomena, as well [33][34][35][36]), the theoretical and conceptual paraphernalia usually underlying direct calculations of the TBR is somewhat oversimplified, if not even rudimentary. The most commonly used theory frames, namely the acoustic (AMM) [21,37] and diffuse (DMM) [21,38] mismatch models, in fact neglect the actual atomic-scale structure of the interface and also estimate, under the Debye approximation, the phonon dispersion branches as linear. They further assume that the phonon interface scattering is purely elastic, either if it is guessed to be diffusive (DMM) or specular (AMM). None of the assumptions is completely fulfilled by real HT and, therefore, these models fail in predicting quantitatively the TBR in most cases. Occasionally, improved versions of AMM and DMM have been conceived [39][40][41], including some of the features neglected in the original models, but no major step forward in understanding the underlying physics was reached without including an all-atom treatment of the interface.
A rather different approach is based on lattice dynamics calculations of TBR, where either the spectral density of phonon transmitted across the interface is computed [42] or the phonon interface scattering is calculated directly (see, e.g. Ref. [43] and references therein). This is typically done under the assumption of elastic scattering (i.e. anharmonicities are disregarded) and neglecting the actual junction width. This amounts to approximate the TBR as a 'junction thermal resistance' and, therefore, both the left and right material segments emit phonons to the junction, as well as they absorb phonons from it. A corresponding net heat flux can be calculated by following a Laundauer-like approach, either by assuming that the phonon distributions are equilibrium (i.e. Bose-Einstein) ones or by assuming non-equilibrium (but bulk-like) expressions for them. Although more fundamental and superior than AMM and DMM (for instance, it directly accounts for quantum effects in the phonon population), even this approach does not provide a satisfactory quantitative prediction of TBR, as found by confronting its predictions for a symmetrically strained Si/Ge interface to the results of a MD simulation with no guess about the phonon scattering events [43]. In addition, in case of non-abrupt interfaces the assumption of specular scattering is questionable. More recently, an anharmonic nonequilibrium Green function approach was developed; this approach overcomes most of the above limitations, but it has been applied to single-molecule junction rather than to solid-solid interfaces [44].
A third and last approach is entirely based on MD simulations, making no direct use of the phonon language, i.e. not requiring any explicit calculation of phonon frequencies, populations, scattering rates or lifetimes. Basically, the TBR of the HT of interest is described as a series of thermal resistances, corresponding to two materials leads A and B embedding an interface layer, whose morphology is fully defined by the sample preparation stage [45,46,48]. The resulting TBR is therefore written as R TBR qq is the overall thermal conductivity of the HT of total thickness L z , while κ A,B qq (l A,B ) are the thermal conductivities of the leads at their actual thickness l A,B . The factor N = 1 or 1/2 reflects, respectively, the fact that the HT is nonperiodic or, rather that periodic conditions are imposed along the z direction. While this approach makes no a priori assumptions about the behaviour of phonons, it is affected by finite-size effects similar to those ones found in any EMD, NEMD or AEMD calculation of transport coefficients [9,[12][13][14]. As a matter of fact, a set of three different calculations is required, namely the l-dependent conductivities of the two materials forming the leads, as well as the overall thermal conductivity of the simulated sample. This could likely result into a heavy computational effort. In any case, the estimated value of TBR actually depends on L z and, therefore, different calculations must be repeated for increasing L z so as to properly extrapolate the TBR value for two semi-infinite leads as, outlined, for instance, in Ref. [46]. Furthermore, in some instances the application of the NEMD protocol is problematic because it does not provide a sizeable T at the interface from which to compute TBR. In these cases, AEMD has proven to be a valuable alternative [47].
In this short review, we highlight an alternative formulation of the TBR problem, based on the non-equilibrium thermodynamics of transport phenomena as for its theoretical formulation, and non-equilibrium MD as for its actual implementation. To our way of thinking, this approach combines at best the merits of a general theory to a robust numerical tool which, albeit operating at the true atomic scale, does not require any simplifying assumption about the interface morphology nor about the physics of the scattering of the thermal energy carriers at such boundary. The method also presents a practical advantage, reducing the computational effort. In Section 2, we briefly outline the basics of non-equilibrium thermodynamics by specifically addressing the problem of TBR. In the next Section 3, we discuss the implementation of this method within NEMD. Finally, in Section 4 we discuss the showcase of a Si/Ge interface, namely the prototypical example of semiconductor HT of paramount important in nanotechnology.

Non-equilibrium thermodynamics theory for heat transport
In order to properly define the conceptual framework for a predictive thermodynamical theory of heat transport across an interface, it is useful to consider at first a homogeneous system (not containing any interface), subject to a temperature gradient. We will further assume that no mass or charge transport phenomena occur, as well as chemical reactions, without any loss of generality since we are focusing on pure heat transport in solid semiconductor materials. If such a homogeneous system is in a steady-state condition of thermal transport and we assume a linear response regime, then local equilibrium holds anywhere and, therefore, all equations of thermodynamics can be cast in local form and all relevant quantities can be given in units of volume (i.e. as a density).
By selecting a volume region in the system, the corresponding change of entropy density ∂s/∂t is given by the sum of the net flow of entropy J s in and out that volume element and an entropy production term σ provided by any possible source inside the same volume (entropy continuity equation) The rate of generation σ is usually expressed in the Onsager form σ = i J i X i [49], i.e. as a sum of products between the ith flux J i and the corresponding generalized force X i . Here the index i spans the different transport mechanisms occurring in the system. By using the Gibbs equation du = Tds (where u is the internal energy density and a constant-volume situation is depicted) and the energy balance equation where T is the temperature. By comparing Equations (1 and 2), we can immediately identify the entropy density production term for the case here considered (no mass or charge transport), which turns out to be given by the product between the heat flux J q and its generalized force X q = ∂ 1/T /∂z. Such a force, in principle, is a linear combination of all fluxes J i occurring in the system. However, in the present case we simply get X q = i r qi J i = r qq J q since, as assumed, no transport phenomena other than the heat current are present. The r qi terms are called Onsager resistivity coefficients: in this context they describe all transport mechanisms and, in Figure 1. Schematic illustration of the appearance of TBR at an interface across which a temperature gradient is established. This cartoon conceptualizes a steady-state condition and, therefore, the temperature profile T (z) far away the interface region is linear (with different slopes at left and right, mimicking an heterojunction between two materials with unlike thermal properties). particular, r qq is ascribed to pure heat conduction (in other words, we remark that the non-diagonal resistivity terms r qi with i = q would describe the possible coupling of heat to charge and mass transport and, therefore, they are null due to the present assumptions). By inserting in the flux-force equation the explicit form of the generalized force provided by Equation (2), we obtain the key result ∂ ∂z where it is shown by very general arguments that the actual thermodynamic driving force for thermal transport is the gradient of an inverse temperature [50]. Although somewhat surprising, this result is consistent with the Fourier equation ∂T/∂z = −R qq J q,z (normally used in predicting thermal transport features in homogeneous materials) by simply developing the z-derivative of the inverse temperature and obtaining ∂T/∂z = −T 2 r qq J q . This unveils the link between the Osanger r qq and the ordinary R qq thermal resistivity, namely T 2 r qq = R qq .

Thermal boundary resistance
The application of the above theory to the case of an interface is not trivial, since the system is no longer homogeneous (see Figure 1). When an interface is present we need, at first, to unambiguously define its location and thickness and, then, properly define any relevant interface thermodynamical quantities. Let us consider the prototypical situation of two semi-infinite material leads meeting at a nominal interface. By selecting any suitable property P(z) having two different values in the bulk-like regions far away from the interface, it is possible to draw its variation along the direction normal to it, as shown in Figure 2, top. This will clearly identify both the left and right interface boundaries (and, therefore, the actual thickness of the interface region) once that it is assumed to define such an interface as that region where P(z) differs from the pure material values (respectively: P left and P right ) by some arbitrary amount. This procedure is named the Gibbs construction [49]. For instance, in the case an interface between two lattice mismatched semiconductors P(z) could be the interplanar lattice constant along z. Alternatively, P(z) could represent the actual content of a dopant, or a given chemical species, or any other structural defect. For a crystalline/amorphous interface P(z) could, finally, represent the local average atomic coordination. As schematically shown in Figure 2, while the P(z) is normally well behaved in the Gibbs interface region, there is in principle no reason for such portion being symmetrically extended into the two facing leads nor to be centred at their nominal interface. The Gibbs construction, therefore, provides a robust definition of the interface which is treated as an autonomous thermodynamical system [49,51].
If we now apply a temperature gradient along z and plot the temperature profile T(z) once the steady-state regime has been reached, we typically obtain what is pictured in Figure 2, bottom. This figure directly defines three relevant temperatures in our problem, namely the left T left and right T right temperatures just at the boundaries of the interface region, as well as the interface temperature T s , defined as the average value in the same volume. This is a quantity easy to compute within a MD simulation using the kinetic energy of the atoms placed in the interface region.
Once that the interface is identified, Equation (1) is recast in the form where the superscript 'int' indicates an interface quantity, while J s,left and J s,left are the entropy density flux through the left and right boundary of the interface, respectively. By developing the same algebra as in the previous section and exploiting the fact that we are addressing a steady-state condition of thermal transport, we eventually obtain the entropy density production term σ int for the interface This result is quite interesting since it shows that two different flux-force equations are indeed necessary to correctly describe the interface problem, namely each one defining its own Onsager resistivity coefficient. Therefore, the total TBR R TBR qq is in fact a series of two Onsager resistances This approach is very clean and robust, since: (i) it is based on very general and elegant thermodynamical arguments; (ii) it does not imply any guess or assumption or approximation about the atomic-scale mechanisms ruling over the thermal energy exchange at the interface; (iii) it does not rely on such a concept as phonon, which only stems from crystals, and therefore can deal with systems lacking of order; (iv) it allows, through the Gibbs construction, to unambiguously define where the interface indeed occurs and how large it is. It is worth clarifying the relationship between the standard definition of Kapitza resistance R TBR qq,Kapitza = (T left − T right )/J q given in Section 1 and the corresponding Onsager definition R TBR qq provided by the above development. As a matter of fact, Equation (7) can be recast in the form where it is clearly shown that the two definitions only differ by the term T 2 s /T left T right . It is now convenient to distinguish between the opposite cases of thick and thin interface. When the interface has a non-vanishing width, like in the case of a rough boundary or when interdiffusion of some chemical species indeed occurs, the ratio T 2 s /T left T right may be significantly different from unity since, as explained above, T s is evaluated through the average kinetic energy of the entire boundary region. Also, if the interface is large enough it could possibly host an additional heat source or a sink, which contributes to the actual value of T s (definitely no longer related to T left or T right in this configuration). The case of a sharp interface is more subtle. For an infinitesimally thin boundary, we can assume with no loss of generality that According to Equation (9) we can state that (11) and therefore (12) indicates that even in the case of a sharp interface, the factor T 2 s /T left T right is not unity when T is an appreciable fraction of the interface temperature, which in turn falls between T left and T right .

Addressing TBR through computer experiments
The numerical implementation of the procedure described in the previous section basically requires a threefold task: (i) the Gibbs construction for the interface; (ii) the set up of a steady state of thermal conduction; and (iii) the evaluation of heat flux in this condition.
The first task is really straightforward and simply requires the calculation of the property P(z). Likely, such a selected property is a structural one and this simply implies that some care must be devoted in preparing the computational sample in a fully relaxed configuration prior to any further calculation.
As for the setting up of a steady-state thermal transport condition, it can be generated by coupling the left and right terminal ends of the system to two heat reservoirs set at different temperatures [12,14,45]. By MD simulation, the system is so aged for a long enough time to reach the steady state which is assessed by a constant-in-time temperature profile (see Figure 2, bottom). We warn against the fact that some details of the profile established across the simulated sample slightly depend on the kind of heat bath used: this, in principle, could somewhat affect the estimation of the temperature drop at the interface (and, therefore, the estimation of T left , T right , and T s ). In general, Langevin thermostatting is recommended since it provides more consistent results with experiments for a large set of simulation parameters [52].
The calculation of the heat current vector is not at all a trivial matter. J q is needed in Equation (7) and it is defined as J q ≡ d/dt( i r i E i ), where r i is the position of the ith particle and E i its energy. By using empirical potentials, it is possible to elaborate non equivalent heat current formulas for the same many-body force field, because of the ambiguity in defining the onsite energy E i . This problem has been recently solved [53] by working out a general pairwise force expression valid for any potential, which avoids the partition of the potential itself into arbitrary single-particle contributions. The same difficulty in uniquely decomposing an ab initio total energy functional (as typically provided by density functional theory) into individual contributions from each atom is usually reported; however, such a misconception has been eventually clarified and a computable expression of the heat current is now available for ab initio MD calculations [54] as well. A different solution to this problem [45,55] consists in calculating instead the work W hot and W cold spent by the hot and cold thermostat, respectively, and evaluating the corresponding heat fluxes as J hot,cold q = (1/S)(∂W hot,clod /∂t), where S is the cross section of the simulated sample. The steady-state condition is now proclaimed when J hot q = J cold q to within the accepted numerical error. In this way, there is no need to make use of any atomic-scale formulation for the heat current.
Whatever solution is adopted to calculate the key ingredients of Equation (7), we remark that as many as O(10 6 ) MD steps could be needed for reaching the steady-state condition [45,55]. Therefore, the direct calculation of R TBR qq through non-equilibrium thermodynamics, while simple in principle, is made non trivial by such a heavy computational demand. Nevertheless, the method outlined in Section 2.2 offers the advantage of requiring the calculation of the relevant quantities just for the HT system, without need to calculate the corresponding properties for the two leads. This translates into a dramatic reduction in the overall computational workload as compared, for instance, to the method based on the treatment of the TBR problem as a series of thermal resistances.

TBR at the Si/Ge interface
Thermal transport across Si/Ge HT can be tailored by engineering their superperiodicity [56], as well as the stoichiometry of the barrier layers [57,58]. While high-frequency phonons are efficiently scattered by Si-Ge alloying, midand low-frequency ones are affected by a suitable distribution of Si/Ge interfaces. Overall, this state of affairs makes SiGe superlattices systems with tunable thermal conductivity, a feature useful for the design of thermoelectric generators [30][31][32] or nanocooling in Si-based devices [23,24]. The thermal resistance at the Si/Ge interface represents in this framework the key feature, which is here addressed in order to show the potential of the theory outlined in the previous section in a case of great practical interest.
In Figure 3, we report the Gibbs construction for a planar Si/Ge abrupt interface. Inspired by experimental work [32], we have modelled the growth of a (001)-oriented Si/Ge HT on crystalline Ge; therefore, the in-plane lattice constant was set at a Ge 0 = 5.6567 Å, namely the bulk-like value predicted for Ge by the adopted Tersoff force field [59]. The cross section of the sample was 5a 0 × 5a 0 . The Si slab is modelled as pseudomorphic (p-Si), meaning that its interplanar spacing is given by a where α is the ratio between the lattice constant in bulk Ge and bulk Si, while C's are the elastic moduli of bulk Si. After a careful conjugated gradient structural relaxation, the resulting sample length was 108.48 nm. For this configuration, the first nextneighbour distance d fn was calculated along the z direction as an average taken over a passing window as wide as a Ge 0 , corresponding for the present case to the property P(z) discussed in Section 2.2. Far away from the nominal (or, equivalently, chemical) interface, d fn recovers the Ge (left) and p-Si (right) bulklike values, as expected, while the central segment where such a distance deviated by more than two standard deviations from the reference values was selected as the interface region. This procedure proves that the interface, while chemically abrupt, has in fact a finite thickness of about 17 Å; interestingly enough, it is also observed that the interface region is not symmetrically spread in the two facing leads but, rather, it is almost entirely hosted by the Si one. This is a system where the chemical interface does not necessarily overlap the thermodynamical one, as anticipated in Section 2.2. All data needed for computing the TBR are contained in Figure 4 where the corresponding temperature profile is shown, as calculated during a NEMD simulation lasted for 5 ns (the first 2 ns are used to set up the steady-state condition). In particular, we get T left = 296.8 K, T right = 271.7 K, and T s = 284.3 K; the calculated stationary heat current is J q = 8.32 GW/m 2 for a nominal temperature offset of 300 K between the hot (Ge side) and cold (Si side) thermostats. Overall, through Equation (7), they provide R TBR qq = 3.02 m 2 K/GW. This result is nicely consistent with previous calculations [43,46] based on the more conventional NEMD prediction performed with the same interatomic potential, although in somewhat different structural conditions (in Ref. [43] symmetrically strained interfaces were considered, while in Ref. [46] no interface relaxation was allowed, while clamping atomic planes at the position predicted by continuum elasticity for a pseudomorphic configuration).
The application of the Gibbs construction and the corresponding calculation of the TBR is very robust, since it does not depend on the actual structure of the investigated systems, nor on the thermal bias conditions. In order to prove this, we investigated a Si/Ge interface occurring in a nanowire with total length of 40 nm and diameter of 5 nm, where the hot and cold thermostats are coupled to the Si and Ge end, respectively. The nominal temperature offset is now larger, i.e. 400 K. The NW system, at variance with the previous bulk-like interface, can laterally accomodate the lattice mismatch between the component materials and, therefore, a situation inherently different than the pseudmorphic one is here experienced. Figure 5 summarises both the Gibbs construction (see inset) and the resulting temperature profile in the steady state. The resulting TBR is now R TBR qq = 2.55 m 2 K/GW. The difference with respect to the bulk-like case is easily accounted for by considering that the interface temperature T s is 284.3 K and 323.4 K for the bulk-like and nanowire system, respectively. The key point here is that, thanks to the Gibbs construction above, the interface corresponds to an autonomous thermodynamical system; therefore, its thermal resistance can be treated as a system variable, only depending upon the interface temperature T s (calculated, once again, through the Gibbs construcution). In other words, as discussed in [51], T s in fact represents a tunable interface property suitable to fully engineer the thermal resistance at a Si/Ge boundary. Interface alloying, here not considered, provides an additional tool [45,46,58].

Conclusions
In this short review, we have readdressed the problem of predicting the TBR at a semiconductor/semiconductor interface by computer simulations. We have shown that the most general paradigm to formulate the problem is provided by non-equilibrium thermodynamics which also suggests the details of an effective protocol for its actual calculation. The Gibbs construction in fact allows for the determination of the interface position and thickness for any given HT, while the calculation of the temperature values at its boundaries and the average value within it allocates all the quantities needed to calculate TBR in a steady state of thermal conduction (the heat flux is needed, as well, similarly to other methods). The present formulation draws a direct link between the Onsager and the phenomenological resistivity and allows to attribute to the interface region its own temperature, making in fact such an interface an autonomous thermodynamical system according to nonequilibrum thermodynamics. While here not addressed, this issue could be important when studying heat transport across an interface in a non-stationary regime and/or for assessing the dependence of boundary resistance upon the local temperature at the interface.

Note
1. In bulk systems with transport of heat, mass and electric charge the three processes are, in general, not adequately described by the simple flux equations usually known as Fourier, Fick and Ohm law, respectively. Rather, a system of flux-(generalized)force is derived within the Onsager theory of transport phenomena, where all transport mechanisms are in fact coupled. By labelling the heat, mass and charge transport by q, m and c, respectively, a set of three transport coefficients L qq , L qm and L qc are accordingly obtained, from which restitivities or conductivities can be derived straightforwardly. In this short review, we will focus on a very special case, namely pure heat transport, as indeed adequate for solid non-metallic materials. We have nevertheless indicated all the relevant heat transport coefficients by the double qq label in order to better frame the present discussion into the more general Onsager theory.