Calculation of Kapitza resistance with kinetic equation
CCalculation of Kapitza resistance with kinetic equation
A. P. Meilakhs ∗ and B. V. Semak Ioffe Institute, 26 Politekhnicheskaya, St. Petersburg 194021, Russian Federation (Dated: June 30, 2020)A new method is introduced for calculation of interfacial thermal resistance in the case of heattransport through the interface by phonons. A unique feature of the method is taking into accountall the consequences of a non-equilibrium character of phonon distribution functions during heattransfer.We introduce a model set of transmission and reflection amplitudes of phonons at theinterface based on the most common in the literature Diffuse Mismatch Model. For the proposedmodel we derive an exact analytical solution. The problem is also solved for a set of transmissionand reflection amplitudes characterized by a free parameter.
I. INTRODUCTION
When heat flows through a boundary between twomedia a temperature jump occurs at the interface. Aproportionality coefficient between the heat flux and thetemperature jump is called interfacial thermal resistanceor Kapitza resistance [1]. Sometimes the inverse valuecalled Kapitza conductance is referred to. Investigationof the Kapitza resistance effect attracts a lot of atten-tion in recent years for it is essential for the optimizationof heat management of many applied systems includingnanostructured materials for heat sink [2, 3], nanofluids[4, 5], thermoelectric generators [6, 7], and nanoelectron-ics [8, 9].Two most important models for the analytical descrip-tion of heat transport through an interface are AcousticMismatch Model (AMM) and Diffuse Mismatch Model(DMM) [10]. In the first case, it is assumed that trans-mission and reflection amplitudes of phonons at the inter-face can be calculated with elasticity theory [11, 12]. Inthe second case, it is assumed that interfacial scatteringis very strong and the phonon incident on the interface”forgets” its initial direction and is scattered uniformlyin all directions [13]. In both models, the distributionfunction of the phonons incident on the boundary is as-sumed to be an equilibrium distribution function with thecorresponding temperature. Meanwhile, the presence ofconstant heat flux across the interface means the pres-ence of the same flux in the media, and, consequently,the phonon distribution functions differ from the equilib-rium one.In numerous contemporary papers on Kapitza resis-tance theory computer simulations of lattice dynamicsat the interface are mostly considered [14]. Such ques-tions as an effect of anharmonic scattering [15] and lat-tice symmetry [16] on phonon transmission through aninterface, accuracy of DMM assumptions [17], a criticalangle concept [18], and influence of chemical bonds onphonon transmission coefficient [19] are investigated bycomputer modeling. These works significantly improveour understanding of the lattice dynamics at the inter-face; yet, questions about the form of the phonon distri- ∗ A. P. Meilakhs: [email protected]ffe.ru bution functions are not considered. Papers on phononkinetics at the interface [20, 21] remain quite rare.However, the questions of kinetics are fundamentallyimportant for understanding the phonon heat transportthrough the interface. As it was first shown in Ref. [12],if we formally apply the AMM formulae to calculate theKapitza resistance of some crystallographic plane in anideal homogeneous crystal, the result would be non zero.It is obvious that the correct calculation result should beexactly zero for such an imaginary boundary, on whichthe reflection of phonons does not occur. This paradoxwas resolved for the first time in Ref. [22] for the one-dimensional model of two bound harmonic strings, witharguments quite similar to those used to describe electri-cal resistance of one-dimensional lattices [23]. It turnedout that the source of the paradox is that AMM doesnot consider a non-equilibrium type of phonon distribu-tion function. It is a clear manifestation that the kineticassumption of the AMM (and also the DMM) is incor-rect: the distribution function of phonons incident onthe boundary is not an equilibrium distribution functionwith the corresponding temperature.On the other hand, in Ref. [24] the question of the cor-rect definition of the temperature of phonons on oppositesides of the interface is raised, given that the vibrationalmodes include the atoms of both crystals when the crys-tals are in contact: ”Since the metal is bonded to theinsulator, the phonon normal modes in equilibrium in-volve all of the atoms in both systems. However, duringheat flow, the phonons in the insulator are at a differ-ent temperature from those in the metal. This situationcannot be described by using the normal modes of thecombined system.” In a paper [25] the definition of thetemperature of phonons on opposite sides of the interfaceis proposed.The defenition is based on the transition from the ba-sis of combined modes to the basis of ordinary phonons,that is, those that correspond to an ordinary flat waveand only in one of the crystals. Wherein, the amplitudesof reflection and transmission of phonons at the interfaceare interpreted as the matching conditions for the distri-bution functions of phonons at the interface. With thisapproach, it is possible to create a formalism that natu-rally takes into account the nonequilibrium phonon dis-tribution functions at the interface. For one-dimensionallattice this approach gives the same result as in paper a r X i v : . [ c ond - m a t . o t h e r] J un ΔT R ΔT B ΔT L ΔT x0T
FIG. 1. The temperature of the phonons near the interfacewith the heat flux from left to right. Shown are: the inter-face plane passing through zero along the coordinate x , solidcurves show the course of the temperature of the phononsnear the interface, the dashed straight lines show linear ex-trapolations of the course of the temperature to the interface.Denoted are the temperature jump at the interface associatedwith the reflection of phonons at the interface ∆ T B and theeffective contributions to the temperature jump associatedwith the nonoptimality of heat transport near the interface∆ T L,R and the full temperature jump ∆ T . [22].Another manifestation of the non-equilibrium processof heat transport through the interface is an additionalcontribution to Kapitza resistance arising from the per-turbation of the phonon distribution function by the in-terface [26]. Since the distribution function of phonons atthe interface is perturbed by interfacial scattering, heattransport near the interface is not optimal, which leadsto an additional effective contribution to the interfacialthermal resistance (Fig. 1). This contribution was in-vestigated for various systems including the interface be-tween normal and superconducting phases in type-II su-perconductors [27, 28], interface between metal and in-sulator [29], and dielectric nanolayers [30]. Since bothcontributions to Kapitza resistance, the one associatedwith the reflection of phonons at the interface and theone associated with the nonoptimality of heat transportnear the boundary, are crucially dependent on the formof the distribution function of phonons, they should becalculated simultaneously, within the framework of a uni-fied system of equations.In the present paper, the approach presented in Ref.[25] is generalized to the three-dimensional case. In theframework of the proposed method, both contributions toKapitza resistance are calculated simultaneously: the one associated with the reflection of phonons at the interfaceand the one associated with the nonoptimality of heattransport near the boundary. II. SYSTEM OF EQUATIONS
In analytic consideration of the problem of the phonontransmission through the crystal boundary, both usingthe elasticity theory [31] and taking into account theatomic structure of the interface [32–34], the solution issought, on the one side of the boundary, as a superpo-sition of incident and reflected plane waves, and on theother side, in the form of a superposition of transmittedplane waves. The solution in this form, which includesoscillations of atoms on both sides of the interface, will becalled the combined mode. Let us introduce the notationfor the parameters of the combined modes.Consider the interface of two semi-infinite crystals: theleft one ( L ) and the right one ( R ). Let a unit amplitudewave, with a frequency ω , a component of the wave vec-tor in the plane of the interface q || and a polarization j ,incident onto the interface from the left side. We denoteas α the whole set of indexes characterising the wave: α = ( q || , ω, j ). The incident wave is partly reflected withan amplitudes A Lαβ and partly transmitted with an am-plitudes B Rαβ , where the first index stands for parametersof the incident wave and the second one for the parame-ters of the departing wave (Fig. 2a).Whatever model of the phonon transmission across theinterface we consider, the law of conservation of energyflow must be fulfilled for the combined mode. In otherwords, the energy flux to the interface and from the in-terface for each mode must be equal. Otherwise, thesolution in the form of combined mode would not be sta-tionary – the interfacial atoms would have to increase ordecrease the amplitude of oscillations over time. The en-ergy flux in the direction perpendicular to the interface, x , is proportional to the wave velocity in the x direction,the density of the substance, the square of the frequency,and the square of the wave amplitude. In this paper, wewill not consider inelastic phonon scattering, so we canreduce the frequency. Thus, we have: ρ L v Lx,α = ρ L (cid:88) β | A Lαβ | v Lx,β + ρ R (cid:88) β | B Rαβ | v Rx,β ρ R v Rx,α = ρ R (cid:88) β | A Rαβ | v Rx,β + ρ L (cid:88) β | B Lαβ | v Lx,β . (1)This property of the combined mode will be used later.It is impossible to correctly determine the crystal tem-peratures using the basis of combined modes since theenergy of each combined mode includes the vibrationalenergy of both crystals. Following the method proposedin [25], we transit from the basis of the combined modesto the basis of ordinary plane waves (Fig. 2b). Herewith,the squares of the reflection and transmission amplitudes A L,Rαβ , B
L,Rαβ become the coefficients for the equation ofmatching of the distribution functions of phonons - plane A αβ β βα n β n α n α n β L B αβ R N α L N α R L ← L → R←R→ FIG. 2. a) Shown are: unit amplitude wave, incident onto the interface from the left side, with parameters α , transmitted andreflected waves with parameters β and amplitudes A αβ and B αβ correspondingly. For simplicity drawn are only one reflectedand one transmitted wave, while there can be infinitely many such waves because of scattering. b) Shown is the transition fromthe basis of combined modes with occupation numbers N Lα , N Rα to the basis of ordinary flat waves with occupation numbers n L → α , n L ← β , n R → α , n R ← β . Indexes → , ← characterize the direction of the phonon movement, from left to right and right to left,respectively. Occupation numbers of departing phonons n L ← β , n R → α are related to occupation numbers of incident phonons n L → α , n R ← α by distribution functions matching equations (2). waves. In the discrete case, these equations are writtenas: n L ← α = (cid:88) β | A Lβα | n L → β + ρ L ρ R (cid:88) β | B Lβα | n R ← β n R → α = (cid:88) β | A Rβα | n R ← β + ρ R ρ L (cid:88) β | B Rβα | n L → β . (2)The physical sense of the matching equations (2) isvery simple: phonons flying from the interface n L ← α , n R → α consist of those that fell on the interface from the sameside and were reflected and of those fallen from the otherside and transmitted n L → α , n R ← α . Since the number ofphonons in a given mode is proportional to the energyin a given mode, and the energy is proportional to thesquare of the amplitude, it is clear that it is the squares ofthe amplitudes that should be used as coefficients. Therigorous derivation of these equations is a straightforwardgeneralization of the derivation of the matching equationsin the one-dimensional case from Ref. [25].To describe the distribution function of phonons nearthe boundary, we use the conventional stationary Boltz-mann equations in the relaxation time approximation: v L,Rx ∂n L,R ∂x = − χ L,R τ L,R . (3)Here χ L,R = n L,R − n L,R is a nonequilibrium part ofthe distribution functions, n – Bose-Einstein distribu-tion function.There is an ambiguity in the division of the distributionfunctions into the sum of the equilibrium and nonequi-librium parts. To eliminate this ambiguity, we introduce the Chapman-Enskog conditions [35]: the temperatureof a nonequilibrium system is defined as the temperatureof an equilibrium system with the same energy, whichyields (cid:88) j (cid:90) d kχ L,R (cid:126) ω = 0 . (4)In order to connect the phonon distribution functionswith the heat flux, we write (cid:88) j (cid:90) d kv x (cid:126) ωχ L,R = q. (5)Since the heat flux through the boundary is a conservedquantity, the last two equations are equivalent and canbe restricted to an equation only for phonons in the leftcrystal. It is easy to verify that these equations are for-mally equivalent since the distribution functions obey thelinking equations (2) whose coefficients satisfy the flowconservation equation (1).The matching equations for the phonon distributionfunctions (2) and the Boltzmann equations for phononson both sides of the interface (3), together with theEnskog-Chapman conditions (4) and the expression forheat flux (5), make up the complete system of equationsfor the phonon distribution functions on both sides of theinterface.It is easy to verify that the result of calculation of theKapitza resistance for an imaginary boundary – someplane in a homogeneous crystal, with the system of equa-tions (2-5) is exactly zero (unlike the AMM calculations).In this case A L,Rαβ = 0 B L,Rαβ = δ αβ . (6)For non-equilibrium distribution functions, we make asubstitution that satisfies the Boltzmann equations (3): χ L = χ R = − τ v x ∂n ∂T dTdx , dTdx = q/κ. (7)Heat flux equation − (cid:88) j (cid:90) d k (2 π ) (cid:126) ωτ v x ∂n ∂T dTdx = q. (8)is satisfied because the expression in parentheses is ther-mal conductivity in the relaxation time approximation.Chapman-Enskog conditions (4) are satisfied, since v x isan odd function of k x . Thus, we have found the non-equilibrium part of the distribution function. Substitut-ing (6) into (2) we find that n L ← α = n R ← β , n R → α = n L → β ,and because of χ L = χ R we find T L = T R = T anddistribution functions n L = n R = n ( T ) − τ v x ∂n ∂T qκ aresolutions of the system of equations (2-5) with the trans-mission and reflection amplitudes defined by conditions(6).Since T L = T R there is no temperature jump associ-ated with phonon reflection ∆ T B = 0. From expression(8) it can be seen that the temperature gradient at theinterface is equal to the temperature gradient far fromthe interface q/κ , hence ∆ T L = ∆ T R = 0. The tem-perature jump on the imaginary boundary vanishes as itshould be. III. MODEL
We use a Debye model of phonon spectrum. We alsointroduce the averaged values of the phonon velocities bythe formula [36]: 3 v L,R = 2 v L,R t + 1 v L,R l . (9)Here v L,Rt is a transversal and v L,Rl is a longitudinal speedof sound in the corresponding crystal.With this in mind, we rewrite the matching equations(2) in a more convenient continuous form and considerthe reflection and transmission amplitudes as functionsof angles, θ (cid:48) – angle of incidence, θ – angle of reflectionor transmission, angles are counted from the x axis per- pendicular to the interface: n L ← = (cid:90) dcosθ (cid:48) n L → ( θ (cid:48) ) | A Lθθ (cid:48) | ++ ρ L v L ρ R v R (cid:90) dcosθ (cid:48) n R ← ( θ (cid:48) ) | B Lθθ (cid:48) | n R → = (cid:90) dcosθ (cid:48) n R ← ( θ (cid:48) ) | A Rθθ (cid:48) | ++ ρ R v R ρ L v L (cid:90) dcosθ (cid:48) n L → ( θ (cid:48) ) | B Rθθ (cid:48) | . (10)The ratio of velocity cubes before the second term inthe right-hand sides arises as the ratio of the densities ofstates in the Debye model when transiting to the contin-uous limit.Let us now choose a set of phonon reflection and trans-mission amplitudes A L,Rθθ (cid:48) , B
L,Rθθ (cid:48) . We take as a basis thesimple for calculations and often used in the literatureDMM [13]. The choice of such a model also justifies theuse of averaged velocities (9), since in this model the de-tails of the phonon spectrum are not important, but onlythe phonon densities of states are important.In DMM, it is assumed that the amplitudes of the re-flected and transmitted waves do not depend on the angleof reflection or transmission, as well as on the angle ofincidence of the wave on the boundary: | A L,Rθθ (cid:48) | = | A L,R | | B L,Rθθ (cid:48) | = | B L,R | . (11)However, we want to consider the coherent scattering ofphonons at the interface (as, for example, in [34]) and thedescription of the lattice dynamics with combined modes.It is easy to see that assumption (11) is not compatiblewith the condition of conservation of energy flow in eachcombined mode (1). Instead of the assumption of theclassical DMM, we introduce the following assumption:the fraction of the energy flux dissipated in a certaindirection does not depend on the angle of incidence ofthe wave on the interface, which yields | A L,Rθθ (cid:48) | = | A L,R | cos θ (cid:48) | B L,Rθθ (cid:48) | = | B L,R | cos θ (cid:48) (12)Let us find the set of amplitudes A L,R , B
L,R .When substituting the equilibrium distribution func-tion in system (3), the system must become an iden-tity for any set of reflection and transmission amplitudes,since at the same crystal temperatures there is no heattransfer between them, and the phonon distribution func-tions are equilibrium.1 = (cid:90) dcosθ (cid:48) | A Lθθ (cid:48) | + ρ L v L ρ R v R (cid:90) dcosθ (cid:48) | B Lθθ (cid:48) | (cid:90) dcosθ (cid:48) | A Rθθ (cid:48) | + ρ R v R ρ L v L (cid:90) dcosθ (cid:48) | B Rθθ (cid:48) | . (13)This equation is an analog to the principle of the detailedequilibrium of the classical DMM [10].For the conservation of flow equation (1), we proceedto the continuous form in the same way as we did for thematching equations: cosθ (cid:48) = (cid:90) dcosθcosθ | A Lθθ (cid:48) | + ρ R v R ρ L v L (cid:90) dcosθcosθ | B Rθθ (cid:48) | cosθ (cid:48) = (cid:90) dcosθcosθ | A Rθθ (cid:48) | + ρ L v L ρ R v R (cid:90) dcosθcosθ | B Lθθ (cid:48) | . (14)We substitute the expression (12) into the equations(13, 14) and get2 = | A L | + ρ L v L ρ R v R | B L | | A R | + ρ R v R ρ L v L | B R | | A L | + ρ R v R ρ L v L | B R | | A R | + ρ L v L ρ R v R | B L | (15)The rank of this system of equations is 3. In order toobtain a complete system of equations for determiningamplitudes A L,R , B
L,R , we add the condition of the in-dependence of the fraction of energy transmitted in acertain direction from the side of incidence (
L, R ), sameas in DMM. Thus | A L | = ρ L v L ρ R v R | B L | | A R | = ρ R v R ρ L v L | B R | (16)The rank of the system of equations (15, 16) is 4, sowe can find | A L | = 2 v R v R + v L | A R | = 2 v L v R + v L | B L | = 2 v R v R + v L ρ R v R ρ L v L | B R | = 2 v L v R + v L ρ L v L ρ R v R (17)This is the desired model set of transmission and reflec-tion amplitudes.We will also assume that the phonon relaxation timesare independent of the frequency τ L,R = const ( ω ). Thisassumption is strictly fulfilled for strongly non-perfect[37] crystals, and is thus justified for areas near theboundary. In addition, as it will be seen later, the fi-nal answer will not depend on τ L,R .We will call the system of equations (2-5) with ampli-tudes (17) and the Debye model of the phonon spectrumand the constant relaxation times ”kinetic DMM”.
IV. SOLUTION
Let us find the solution of kinetic DMM. We first findthe distribution functions of the phonons incident on theinterface by solving the stationary Boltzmann equation(3). Due to the perturbation of the distribution functionof phonons by interfacial scattering, we must considerthe non-equilibrium part of the distribution function asa function of the coordinate χ L,R ( x ). For the left crystal ∂χ L ∂x + χ L v Lx τ L = − ∂n L ∂x . (18)This is an ordinary differential equation for x , and generalsolution of the non-homogeneous equation is the sum ofa particular solution and the complementary solution ofthe associated homogeneous equation χ L = χ Lp + χ Lc .Complementary solution is: χ Lc = χ L exp ( − x/v Lx τ L ) , (19)where χ L = χ L ( x = 0). Similarly for the right crystal χ Rc = χ R exp ( x/v Rx τ R ).For incident phonons v Lx > , v Rx <
0, the solution ofthis type increases infinitely. Thus, it turns out that thecomplementary solution for the distribution function ofincident phonons does not satisfy the boundedness condi-tion, which means that the distribution function of inci-dent phonons is determined only by a particular solution.Particular solution is χ L,Rp = − τ L,R v L,Rx ∂n L,R ∂T (cid:18) dTdx (cid:19) L,R . (20)Here (cid:0) dTdx (cid:1) L,R – unknown functions of coordinate, since,due to the perturbation of the phonon distribution func-tions by the interface, the temperature gradients nearthe boundary differ from the gradients in a homogeneousmedium for the same value of the heat flux. We introducethe parameters ˜∆ T L,R = τ L,R v L,R (cid:0) dTdx (cid:1)(cid:12)(cid:12)(cid:12)
L,Rx =0 which canbe considered as preliminary estimates for ∆ T L,R , theeffective contributions to the temperature jump relatedto the nonoptimality of heat transport near the interface.We consider linear heat transport across the border,thus ∆ T B (cid:28) T L = T R + ∆ T B = T . In this case,through the entered parameters, the phonon distributionfunctions on the boundary ( x = 0) can be written as, n L → = n ( T ) + ∂n ∂T ∆ T B + ∂n ∂T ∆ ˜ T L cos θ (cid:48) n R ← = n ( T ) + ∂n ∂T ∆ ˜ T R cos θ (cid:48) . (21)Using the matching equations (2), we express the dis-tribution functions of phonons flying from the interfacethrough the distribution functions of phonons incident onthe interface. We will use the matching equations in theform (10). Substituting (21) into (10) with amplitudes(17) we find n L ← = n R → = n ( T ) + v R v R + v L ∂n ∂T ∆ T B ++ 23 v R v R + v L ∂n ∂T ∆ ˜ T L + 23 v L v R + v L ∂n ∂T ∆ ˜ T R . (22)The distribution functions of phonons near the boundaryare thus expressed through three unknown parameters:∆ ˜ T L,R , ∆ T B .In order to finally find the distribution functions ofphonons near the boundary, we substitute (21, 22) intoequations (4, 5) and obtain a system of three equationswith three unknowns, from which we find∆ T B = 8 π (cid:126) qk B v R + v L T I T D , (23)and ∆ ˜ T L,R = 16 π (cid:126) qk B v L,R T I T D , (24)where T D is the Debye temperature in that of the crys-tals, in which the Debye temperature is lower, I T T = (cid:90) T /TT /T y e y ( e y − dz. (25)is a dimensionless integral depending on tempera-ture. For the corresponding integration limits, is ex-pressed through the Debye function F D as I T D = T D / (3 T ) F D ( T D /T ).Let the Debye temperature be lower in the left crystal.To find ∆ T L , we investigate the damping of the pertur-bation of the distribution function near the interface. Tofind the complementary part of the nonequilibrium dis-tribution function at the interface, we use the expression χ L = n L ← − n ( T ) − ∂n ∂T ∆ T B − ∂n ∂T ∆ ˜ T L cos θ (cid:48) . (26)Substitute the expression for the complementary (19, 26)and particular (20) parts of the nonequilibrium distribu-tion function in the heat flow expression (5) and find: dTdx = − qκ − (cid:18) dTdx (cid:19) pert , (27)where (cid:0) dTdx (cid:1) pert is a change in temperature gradient dueto perturbation by the interface: (cid:18) dTdx (cid:19) pert = 3 (cid:90) d kv Lx (cid:126) ωχ Lh exp ( − x/v Lx τ L ) . (28)Because of ∆ T L = − (cid:90) −∞ dx (cid:18) dTdx (cid:19) pert , (29) after substituting relation (28) into (29) and integratingover x , we find ∆ T L = 18 ∆ ˜ T L . (30)The calculation of ∆ T R is performed similarly, withthe difference that for phonons with frequencies greaterthan the maximum phonon frequency in the left crystal, k B T LD / (cid:126) coefficients in the matching equation (17) mustbe changed so that B L = 0, since such phonons do notpass through the interface. We obtain∆ T R = 18 I T LD I T RD ∆ ˜ T R + 616 I T RD T LD I T RD ∆ ˜ T R , (31)where the second term is the relaxation of the distribu-tion function of those phonons that are not included inthe heat transport through the interface. V. RESULTS AND DISCUSSION
T, K G , W / K c m Exp. Al O -Au Kinetic DMM, no HFP Kinetic DMM Classic DMM FIG. 3. The dependence of the Kapitza conductance at theboundary of sapphire and gold on temperature. The pointsare experimental data from Ref. [38]. The solid line is theresult of calculations with kinetic DMM, formula (32). Thedashed line is the result of calculations with classic DMM. Thedash-dotted line is kinetic DMM without taking into accountthe contribution of the relaxation of high-frequency phonons(HFP). Data is presented on a normal scale.
Thus, the Kapitza conductance is expressed as G = q ∆ T B + ∆ T L + ∆ T R == 72 π k B (cid:126) T I T D v R + v L ) + 2 v R (cid:16) I T RD T LD /I T RD (cid:17) . (32)It is interesting that, just like in classic DMM, in kineticDMM, the Kapitza conductance turns out to be indepen-dent of the density of crystals. Also, from the expression(31) it is clear why the averaging of the velocities shouldbe carried out according to the formula (9). G , W / K c m T, K
Exp. Al O -Al Kinetic DMM Classic DMM FIG. 4. The dependence of the Kapitza conductance at theboundary of sapphire and aluminum on temperature. Thesolid line is the result of calculations with kinetic DMM, for-mula (32). The dashed line is the result of calculations withclassic DMM. The points are experimental data from Ref.[38]. Data is presented on a logarithmic scale.
The contribution of the relaxation of high-frequencyphonons 2 v R (cid:16) I T RD T LD /I T RD (cid:17) is small for a pair of substanceswith approximately equal values of the Debye tempera-ture. However, if the Debye temperatures of two sub-stances are significantly different, this contribution canbe quite significant, given that the speed of sound is usu-ally higher in the substance where the Debye tempera-ture is higher. In the case when the results of kineticDMM coincide well with the experimental results, takinginto account the contribution from the relaxation of high-frequency phonons significantly improves the coincidence(Fig. 3).If we remove from the expression (31) the contributionof the relaxation of high-frequency phonons (HFP), weobtain the formula, which differs from the classical DMMformula, when averaging (9) is applied to it, only by anumerical factor. The result of calculations using kineticDMM without HFP is approximately 1 . t = | B R | . Hav-ing done the calculations for this model, similar to those FIG. 5. The dependence of the Kapitza conductance ontemperature. The solid line is the result of calculations withkinetic DMM, formula (32). The dashed line is the result offitting with the free parameter t – the transmission coefficient,formula (34). Data is presented on a normal scale. a) Theboundary of sapphire and titanium nitride. The points areexperimental data from Ref. [39]. The optimal value of thefree parameter is t = 0 , t DMM . b) The boundary of cobaltdisilicide and silicon. The points are experimental data fromRef. [40]. The optimal value of the free parameter is t =0 , t DMM . c) The boundary of sapphire and titan. Thepoints are experimental data from [38]. The optimal value ofthe free parameter is t = 0 , t DMM . above, we find G B = q ∆ T B = 38 π k B (cid:126) T I T D t ρ L v L ρ R v R v L − ( v R + v L ) t . (33)We can see that the dependence on the density of crystalsappeared in the new expression. Thus, it turns out thatthe independence of the Kapitza conductance (31) ondensities was obtained due to the acceptance of assump-tion (16). Unexpectedly, the magnitudes of temperaturejumps associated with the nonoptimality of heat transfernear the boundary ∆ T L,R are independent of t and aredetermined by expressions (24, 29, 30).Expression (33) equals to (32) at t = ρ L v L ρ R v R v L v R + v L .We will call this value t DMM . Let us express G B interms of t DMM G B = 38 π k B (cid:126) T I T D t ( v R + v L )( t DMM − t ) . (34)It can be seen that the resistance associated with thereflection of phonons from the boundary goes to infinity,or, equivalently, the temperature jump at the boundaryvanishes at t = 7 / t DMM . A feature of the method thatyields a zero value of Kapitza resistance at the imaginaryboundary is shown in the singular nature of the formulae(33, 34). G , W / K c m T, K
Exp. Diam.-Au Only HFP Kinetic DMM Classic DMM
FIG. 6. The dependence of the Kapitza conductance at theboundary of diamond and gold on temperature. The solidline is the result of calculations with kinetic DMM, formula(32). The dashed line is the result of calculations with classicDMM. The dash-dotted line is kinetic DMM at t = 7 / t DMM ,thus taking into account the only the relaxation contribution.The points are experimental data from Ref. [38]. Data ispresented on a normal scale.
The nonlinear dependence of Kapitza conductance G B on the transmission coefficient t (33), is a pos-sible cause of a substantial discrepancy between the-ory and experiment in some cases. For example, at t = 3 / t DMM , the conductivity decreases by almost two times G B (3 / t DMM ) = 9 / G B ( t DMM ). Figure (5)shows the dependencies of G ( T ) on the temperature at t which gives the best approximation of the experimen-tal data. It turns out that for most pairs of boundariesa good approximation is achieved by varying t within60% t DMM .The contribution of the relaxation of high-frequencyphonons to the calculated value of Kapitza resistance issmall at low temperatures and is growing with the tem-perature grows. That leads to a slower growth rate ofthe overall Kapitza conductance with temperature. Insome cases, it leads to a very good coincidence withthe experimentally measured temperature dependence ofKapitza conductance (Fig. 3, Fig.5c). However, in somecases it can be seen (Fig. 5a, Fig. 5b) that in thelow-temperature region the theoretical curves pass be-low the experimental points, but at temperatures closeto room temperature, on the contrary, turn out to behigher. We can offer two possible explanations for thisdeviation of experimental data from the temperature de-pendence predicted by theory. It can be assumed thatthe DMM approach is not correct, the scattering at theboundary is rather weak, and in this case, as shown bytheoretical models [32, 34], the phonon transmission coef-ficient should decrease with increasing phonon frequency.This would explain the slower increase in the Kapitzaconductance with increasing temperature, at which high-frequency oscillations are excited. It is possible, on thecontrary, to assume that the DMM assumption about thediffuse nature of scattering is correct, and the discrep-ancy between theory and experiment is associated withthe use of the Debye approximation. In Ref. [41] theKapitza conductance was studied in the DMM model,but with a more realistic dependence of the density ofphonon states on the frequency than that in the Debyemodel. As a result, it turns out that at low frequencies,Kapitza conductance grows faster with increasing tem-perature than classical DMM predicts, and at high tem-peratures, on the contrary, it grows more slowly, whichexplains the observed discrepancy.In the framework of kinetic DMM we can not introducesuch a concept as ”Radiation limit” [42] since, as it wasshown in a previous section, the temperature jump on theimaginary boundary vanishes. However the magnitudesof temperature jumps associated with the nonoptimal-ity of heat transfer near the boundary ∆ T L,R are inde-pendant of t so we can introduce the ”Relaxation limit”which is kinetic DMM at t = 7 / t DMM so the temper-ature jump associated with the reflection of phonons atthe interface vanishes, thus taking into account the onlythe relaxation contribution. For pairs of substances withvery different sound velocities, such as, for example, di-amond - plumbum, diamond - aurum, kinetic DMM, aswell as classic DMM, predict a significantly lower Kapitzaconductance value than that measured in the experiment(Fig. 6). In such cases, experimental data cannot be fit-ted by variation of the free parameter t , because even at”Relaxation limit” the interfacial thermal resistance as-sociated with the nonoptimality of heat transport nearthe boundary is still larger than the experimentally mea-sured. We share the common point of view [24, 43] thatthis discrepancy can be explained only with the help ofa model that takes into account the direct transfer ofheat from the metal electrons to the dielectric phonons[44–46]. VI. SUMMARY
A new method for calculating the Kapitza conductancewas proposed. Matching equations for the phonon dis-tribution functions are introduced, which, together withthe Boltzmann equations for phonons and some addi-tional equations, form a complete system of equationsdescribing the kinetics of heat transfer across the bound-ary by phonons. A method of solution to such a systemwas introduced. It is shown that the use of this method does not lead to a paradox: the solution of the corre-sponding system for an imaginary boundary, that is, fora plane in a homogeneous crystal, gives zero Kapitza re-sistance. For a model set of amplitudes based on theclassical DMM, an exact analytical solution is obtained.The dependence of Kapitza conductance on temperatureis found. The problem was also solved for a set of ampli-tudes with a free parameter characterizing the fraction ofenergy transmitted by a phonon across the boundary. Afeature of the method that yields a zero value of Kapitzaresistance at the imaginary boundary is shown in thesingular character of the formula. It turned out that theselection of the corresponding parameter value leads to agood approximation of the experimental data.The authors thank A. Y. Vul for attention to the work,M. G. Glazov for helpful comments, A. S. Khorolskaya forfruitful discussions. A.P. Meilakhs thanks RSF (grant18-72-00131) for support. [1] P. L Kapitza, J. Phys. USSR , 181 (1941).[2] S. V. Kidalov, F. M. Shakhov, Materials , 46(2013).[4] S. P. Jang and S. U. S. Choi, Appl. Phys. Lett. , 4316(2004).[5] M. A. Serebryakova, A. V. Zaikovskii et. al. J. Heat andMass Transfer , 1314 (2017).[6] E. D. Eidelman, A. P. Meilakhs, Nanosyst.: Phys.,Chem., Math., , 919 (2016).[7] E. D. Eidelman, A. P. Meilakhs, B. V. Semak and F. M.Shakhov, J. Phys. D: Appl. Phys., , 464007 (2017).[8] H. Adhikari, A. F. Marshall, C. E. D. Chidsey, and P. C.McIntyre, Nano Lett. , 318 (2006).[9] Y. Zhou, J. Anaya, J. Pomeroy, H. Sun, X. Gu, A. Xie, E.Beam, M. Becker, T. A. Grotjohn, C. Lee and M. Kuball,ACS Appl. Mater. Interfaces , 34416 (2017).[10] E. T. Swartz, and R. O. Pohl, Rev. Mod. Phys. , 605(1989).[11] I. M. Khalatnikov, JETP. , 687 (1952).[12] W. A. Little, Can. J. Phys. , 334 (1959).[13] E. T. Swartz, and R. O. Pohl, Appl. Phys. Lett. , 200(1987).[14] N. Yang, T. Luo et. al. J. Comp. and Theor. Nanosci. , 168 (2015).[15] K. Saaskilahti, J. Oksanen, J. Tulkki and S. Volz, Phys.Rev. B. , 134312 (2014).[16] K. Bi, Y. Liu, C. Zhang, J. Li, M. Chen, Y. Chen, Appl.Phys. A. , 883 (2016).[17] R. R. Kakodkar and J. P. Feser, Phys. Rev. B , 125434(2017).[18] A. Alkurdi, S. Pailhs, and S. Merabia, Appl. Phys. Lett. , 093101 (2017).[19] Z. Huanga, C. Huanga, D. Wua, Z. Raoa Comp. Mat.Sci. , 316 (2018).[20] S. Merabia, K. Termentzidis, Phys. Rev. B , 094303(2012).[21] G. Varnavides, A.S. Jermyn, P. Anikeeva, P. Narang,Phys. Rev. B , 115402, (2019).[22] A. M. van den Brink, H. Dekker, Phys. Rev. B. , P.17842 (1995).[23] R. Landauer, Philosophical Mag. , (1970). [24] G. D. Mahan, Phys. Rev. B , 075408 (2009).[25] A. P. Meilakhs, Phys. Solid State, , 148 (2015).[26] A. V. Latyshev, A. A. Yushkanov, Theor. Math. Phys., , 105 (2010).[27] A G. Aronov, A. S. Ioselevich, JETP. , 1839 (1981).[28] V. I. Kozub JETP. , 1847 (1985).[29] A. Majumdar and P. Reddy, Appl. Phys. Lett. , 4768(2004).[30] K. Alaili, J. Ordonez-Miranda, Y. Ezzahri, Int. J. Ther-mal Sci. , 40 (2018).[31] L. D. Landau and E. M. Lifshitz, Course of TheoreticalPhysics, Vol. 7: Theory of Elasticity ( Pergamon, NewYork, 1986).[32] L. Zhang, P. Keblinski, J.-S. Wang and B. Li, Phys. Rev.B. , 064303 (2011).[33] D. A. Young, H. J. Maris, Phys. Rev. B. , 3685 (1988).[34] A. P. Meilakhs, Nanosyst.: Phys., Chem., Math. , 971(2016).[35] L. D. Landau and E. M. Lifshitz, Course of TheoreticalPhysics, Vol. 10: Physical Kinetics ( Pergamon, NewYork, 1981).[36] H.-K. Lyeo, D. G. Cahill, , 144301 (2006).[37] J. M. Ziman, Electrons and Phonons. The Theory ofTransport Phenomena in Solids (Oxford UniversityPress, Oxford, 2001).[38] R. J. Stoner, H. J. Maris, Phys. Rev. B , 16373 (1993).[39] R. M. Costescu, M. A. Wall, D. G. Cahill, Phys. Rev. B. , 054302 (2003).[40] N. Ye, J. P. Feser, et. al. Phys. Rev. B. , 085430 (2017).[41] P. Reddy, K. Castelino A. Majumdar, Appl. Phys. Lett. , 211908 (2005).[42] N. S. Snyder, Cryogenics , 89 (1970).[43] M. L. Huberman and A. W. Overhauser, Phys. Rev. B , 2865 (1994).[44] A. P. Meilakhs, E. D. Eidelman, JETP Lett., , 38(2013).[45] T. Lu, J. Zhou, T. Nakayama, R. Yang, B. Li, Phys. Rev.B , 085433 (2016).[46] B. Shi, X. Tang, T. Lu, T. Nakayamay, Y. Liz, J. Zhouz,Mod. Phys. Lett. B32