Calculation of Kapitza resistance with kinetic equation

TL;DR

This paper presents a novel analytical method for calculating interfacial thermal resistance due to phonons, accounting for non-equilibrium phonon distributions and validated against experimental data.

Contribution

It introduces an exact analytical solution for Kapitza resistance considering non-equilibrium phonon effects using a model based on the Diffusive Mismatch Model.

Findings

Results agree well with experimental data

The method accounts for non-equilibrium phonon distributions

Analytical solutions are derived for different transmission models

Abstract

A new method is introduced for calculation of interfacial thermal resistance in the case of heat transport through the interface by phonons. A unique feature of the method is taking into account all the consequences of a non-equilibrium character of phonon distribution functions during the heat transport. We introduce a model set of transmission and reflection amplitudes of phonons at the interface based on the most common in the literature Diffusive Mismatch Model. For the proposed model we derive an exact analytical solution. The problem is also solved for a set of transmission and reflection amplitudes characterized by a free parameter. We found that the calculation results are in a good agreement with the experimental data.

Figures

b). Herewith, the squares of the reflection and transmission amplitudes $A_{\alpha\beta}^{L,R},B_{\alpha\beta}^{L,R}$ become the coefficients for the equation of matching of the distribution functions of phonons - plane waves. In the discrete case, these equations are written as:

c). However, in some cases it can be seen (Fig. 5a, Fig. 5b) that in the low-temperature region the theoretical curves pass below the experimental points, but at the temperatures close to room temperature the theoretical curves turn out to be, on the contrary, higher. We can offer two possible explanations of this deviation of experimental data from the temperature dependence predicted by the theory. It can be assumed that the DMM approach is not correct, the scattering at the boundary is rather weak, and in this case, as shown by theoretical models [[33](#bib.bib33), [35](#bib.bib35)], the ph

Equations

$${\rho }^L{v}_{x,\alpha }^L={\rho }^L\sum _{\beta }|A_{\alpha \beta }^L{|}^2{v}_{x,\beta }^L+{\rho }^R\sum _{\beta }|B_{\alpha \beta }^R{|}^2{v}_{x,\beta }^R$$

$${\rho }^R{v}_{x,\alpha }^R={\rho }^R\sum _{\beta }|A_{\alpha \beta }^R{|}^2{v}_{x,\beta }^R+{\rho }^L\sum _{\beta }|B_{\alpha \beta }^L{|}^2{v}_{x,\beta }^L.$$

$$n_{\alpha }^{L\leftarrow }=\sum _{\beta }|A_{\beta \alpha }^L{|}^2{n}_{\beta }^{L\to }+\frac{{\rho }^L}{{\rho }^R}\sum _{\beta }|B_{\beta \alpha }^L{|}^2{n}_{\beta }^{R\leftarrow }$$

$$n_{\alpha }^{R\to }=\sum _{\beta }|A_{\beta \alpha }^R{|}^2{n}_{\beta }^{R\leftarrow }+\frac{{\rho }^R}{{\rho }^L}\sum _{\beta }|B_{\beta \alpha }^R{|}^2{n}_{\beta }^{L\to }.$$

$$v_x^{L,R}\frac{\partial n^{L,R}}{\partial x}=-\frac{{\chi }^{L,R}}{{\tau }^{L,R}}.$$

$$n_0^{L,R}=\frac{1}{\exp (\hslash \omega /k_B{T}^{L,R})+1}.$$

$$\sum _j\int \frac{d^{3}k}{(2\pi {)}^3}{\chi }^{L,R}\hslash \omega =0.$$

$$\sum _j\int \frac{d^{3}k}{(2\pi {)}^3}{v}_x\hslash \omega {\chi }^{L,R}=q.$$

$$A_{\alpha \beta }^{L,R}$$

$$B_{\alpha \beta }^{L,R}$$

$${\chi }^{L,R}=-{\tau }^{L,R}{v}_x^{L,R}\frac{\partial n_0}{\partial T}\frac{dT}{dx},\frac{dT}{dx}=q/\kappa .$$

$$-\left(\sum _j\int \frac{d^{3}k}{(2\pi {)}^3}\hslash \omega \tau v_x^2\frac{\partial n_0}{\partial T}\right)\frac{dT}{dx}=q.$$

$$\frac{3}{v^{L,R^2}}=\frac{2}{v_t^{L,R^2}}+\frac{1}{v_l^{L,R^2}}.$$

$$n^{L\leftarrow }={\int }_0^{1}dcos{\theta }^{\prime }{n}^{L\to }({\theta }^{{}^{\prime }})|A_{\theta {\theta }^{{}^{\prime }}}^L{|}^2+$$

$$+\frac{{\rho }^L{v}^{L^3}}{{\rho }^R{v}^{R^3}}{\int }_0^{1}dcos{\theta }^{\prime }{n}^{R\leftarrow }({\theta }^{{}^{\prime }})|B_{\theta {\theta }^{{}^{\prime }}}^L{|}^2$$

$$n^{R\to }={\int }_0^{1}dcos{\theta }^{\prime }{n}^{R\leftarrow }({\theta }^{{}^{\prime }})|A_{\theta {\theta }^{{}^{\prime }}}^R{|}^2+$$

$$+\frac{{\rho }^R{v}^{R^3}}{{\rho }^L{v}^{L^3}}{\int }_0^{1}dcos{\theta }^{\prime }{n}^{L\to }({\theta }^{{}^{\prime }})|B_{\theta {\theta }^{{}^{\prime }}}^R{|}^2$$

$$|A_{\theta {\theta }^{\prime }}^{L,R}{|}^2=|A^{L,R}{|}^2$$

$$|B_{\theta {\theta }^{\prime }}^{L,R}{|}^2=|B^{L,R}{|}^2.$$

$$|A_{\theta {\theta }^{\prime }}^{L,R}{|}^2=|A^{L,R}{|}^2\cos {\theta }^{\prime }$$

$$|B_{\theta {\theta }^{\prime }}^{L,R}{|}^2=|B^{L,R}{|}^2\cos {\theta }^{\prime }.$$

$$1={\int }_0^{1}dcos{\theta }^{\prime }|A_{\theta {\theta }^{\prime }}^L{|}^2+\frac{{\rho }^L{v}^{L^3}}{{\rho }^R{v}^{R^3}}{\int }_0^{1}dcos{\theta }^{\prime }|B_{\theta {\theta }^{\prime }}^L{|}^2$$

$$1={\int }_0^{1}dcos{\theta }^{\prime }|A_{\theta {\theta }^{\prime }}^R{|}^2+\frac{{\rho }^R{v}^{R^3}}{{\rho }^L{v}^{L^3}}{\int }_0^{1}dcos{\theta }^{\prime }|B_{\theta {\theta }^{\prime }}^R{|}^2.$$

$$cos{\theta }^{{}^{\prime }}={\int }_0^{1}dcos\theta cos\theta |A_{\theta {\theta }^{\prime }}^L{|}^2+\frac{{\rho }^R{v}^R}{{\rho }^L{v}^L}{\int }_0^{1}dcos\theta cos\theta |B_{\theta {\theta }^{\prime }}^R{|}^2$$

$$cos{\theta }^{{}^{\prime }}={\int }_0^{1}dcos\theta cos\theta |A_{\theta {\theta }^{\prime }}^R{|}^2+\frac{{\rho }^L{v}^L}{{\rho }^R{v}^R}{\int }_0^{1}dcos\theta cos\theta |B_{\theta {\theta }^{\prime }}^L{|}^2.$$

$$2=|A^L{|}^2+\frac{{\rho }^L{v}^{L^3}}{{\rho }^R{v}^{R^3}}|B^L{|}^2$$

$$2=|A^R{|}^2+\frac{{\rho }^R{v}^{R^3}}{{\rho }^L{v}^{L^3}}|B^R{|}^2$$

$$2=|A^L{|}^2+\frac{{\rho }^R{v}^R}{{\rho }^L{v}^L}|B^R{|}^2$$

$$2=|A^R{|}^2+\frac{{\rho }^L{v}^L}{{\rho }^R{v}^R}|B^L{|}^2.$$

$$|A^L{|}^2=\frac{{\rho }^L{v}^L}{{\rho }^R{v}^R}|B^L{|}^2$$

$$|A^R{|}^2=\frac{{\rho }^R{v}^R}{{\rho }^L{v}^L}|B^R{|}^2.$$

$$|A^L{|}^2$$

$$|A^R{|}^2$$

$$|B^L{|}^2$$

$$|B^R{|}^2$$

$$\frac{\partial {\chi }^L}{\partial x}+\frac{{\chi }^L}{v_x^L{\tau }^L}=-\frac{\partial n_0^L}{\partial x}.$$

$${\chi }_c^L={\chi }_0^L\exp (-x/v_x^L{\tau }^L),$$

$${\chi }_p^{L,R}=-{\tau }^{L,R}{v}_x^{L,R}\frac{\partial n_0^{L,R}}{\partial T}{\left(\frac{dT}{dx}\right)}^{L,R}.$$

$$\tilde{\Delta T^{L,R}}=\tau^{L,R}v^{L,R}\left(\frac{dT}{dx}\right)\Bigr{|}_{x=0}^{L,R}.$$

Full text

I Introduction

When heat flows through a boundary between two media a temperature jump occurs at the interface. A proportionality coefficient between the heat flux and the temperature jump is called interfacial thermal resistance or Kapitza resistance [1]. The inverse value called Kapitza conductance is also widely used. Investigation of the Kapitza resistance effect attracts a lot of attention in recent years for it is essential for the optimization of heat management of many applied systems including nanostructured materials for heat sink [2, 3], nanofluids [4, 5], thermoelectric generators [6, 7], and nanoelectronics [8, 9].

Two most important models for the analytical description of heat transport through an interface are Acoustic Mismatch Model (AMM) and Diffuse Mismatch Model (DMM) [10]. In the first case, it is assumed that transmission and reflection amplitudes of phonons at the interface can be calculated with elasticity theory [11, 12]. In the second case, it is assumed that interfacial scattering is very strong and the phonon incident on the interface ”forgets” its initial direction and is scattered uniformly in all directions [13]. In both models, the distribution function of the phonons incident on the boundary is assumed to be an equilibrium distribution function with the corresponding temperature. Meanwhile, the presence of constant heat flux across the interface means the presence of the same flux in the media, and, consequently, the phonon distribution functions differ from the equilibrium one.

Computer simulation of lattice dynamics at the interface is the most common approach to tackle the problems of Kapitza resistance theory described in the contemporary research manuscripts [14]. Such questions as an effect of anharmonic scattering [15] and lattice symmetry [16] on phonon transmission through an interface, accuracy of DMM assumptions [17], a critical angle concept [18], and influence of chemical bonds on phonon transmission coefficient [19] are investigated by computer modelling. These works significantly improve our understanding of the lattice dynamics at the interface; yet, questions about the form of the phonon distribution functions had not been previously considered. Papers on phonon kinetics at the interface [20, 21] remain quite rare.

However, the questions of kinetics are fundamentally important for understanding the phonon heat transport through the interface. As it was first shown in Ref. [12], formal application of the AMM formulae to calculation of the Kapitza resistance of some crystallographic plane in an ideal homogeneous crystal, produces a non-zero result. It is obvious that the correct calculation result should be exactly zero for such an imaginary boundary, on which the reflection of phonons does not occur. This paradox was resolved for the first time in Ref. [22] where the one-dimensional model of two bound harmonic strings was used. In that research the authors applied the arguments quite similar to those used to describe electrical resistance of one-dimensional lattices [23]. It turned out that the source of the paradox is that AMM does not consider a non-equilibrium type of phonon distribution function. It is a clear manifestation that the kinetic assumption of the AMM (and also the DMM) is incorrect: the distribution function of phonons incident on the boundary is not an equilibrium distribution function with the corresponding temperature.

On the other hand, in Ref. [24] the question of the correct definition of the temperature of phonons on opposite sides of the interface is raised, given that the vibrational modes include the atoms of both crystals when the crystals are in contact: ”Since the metal is bonded to the insulator, the phonon normal modes in equilibrium involve all of the atoms in both systems. However, during the heat flow, the phonons in the insulator are at a different temperature from those in the metal. This situation cannot be described by using the normal modes of the combined system.”  The definition of the temperature of phonons on opposite sides of the interface was proposed in [25].

The defenition is based on the transition from the basis of combined modes to the basis of ordinary phonons, that is, those that correspond to an ordinary flat wave and only in one of the crystals. Wherein, the amplitudes of reflection and transmission of phonons at the interface are interpreted as the matching conditions for the distribution functions of phonons at the interface. With this approach, it is possible to create a formalism that naturally takes into account the nonequilibrium phonon distribution functions at the interface. For one-dimensional lattice this approach gives the same result as described previously [22].

Another manifestation of the non-equilibrium process of heat transport through the interface is an additional contribution to Kapitza resistance arising from the perturbation of the phonon distribution function by the interface [26]. Since the distribution function of phonons at the interface is perturbed by interfacial scattering, the distribution function of phonons near the interface is not equal to distribution function describing heat flow in a homogenous media. As it is known, the calculated thermal conductivity is maximal when it is calculated with the exact distribution function describing heat flow in the homogenous media [27]. Thus we can assume that thermal conductivity near the interface is smaller than in the homogenous media. And with the same heat flux, the temperature gradient at the interface is higher which leads to an additional effective contribution to the interfacial thermal resistance (Fig. 1). This contribution was investigated for various systems including the interface between normal and superconducting phases in type-II superconductors [28, 29], interface between metal and insulator [30], and dielectric nanolayers [31]. Since both contributions to Kapitza resistance, the one associated with the reflection of phonons at the interface and the one associated with smaller thermal conductivity near the interface, are crucially dependent on the form of the distribution function of phonons, they should be calculated simultaneously, within the framework of a unified system of equations.

In the present paper, the approach presented in Ref. [25] is generalized to the three-dimensional case. We introduce the full set of equations describing phonon kinetics at the interface. For the first time, we resolve the paradox mentioned in Ref. [12] for a three-dimensional system. We use the DMM-based set of transmission and reflection amplitudes of phonons at the interface, because of its popularity and simplicity. For this simple set of amplitudes, we find the exact solution. For the first time, we simultaneously calculate both contributions to Kapitza resistance, the one associated with the reflection of phonons at the interface and the one associated with smaller thermal conductivity near the interface.

II System of equations

In analytic consideration of the problem of the phonon transmission through the crystal boundary, both using the elasticity theory [32] and taking into account the atomic structure of the interface [33, 34, 35], the solution is sought, on the one side of the boundary, as a superposition of incident and reflected plane waves, and on the other side, in the form of a superposition of transmitted plane waves. The solution in this form, which includes oscillations of atoms on both sides of the interface, will be called the combined mode. Let us introduce the notation for the parameters of the combined modes.

Consider the interface of two semi-infinite crystals: the left one ($L$) and the right one ($R$). Let a unit amplitude wave, with a frequency $\omega$, a component of the wave vector in the plane of the interface $q_{||}$ and a polarization $j$, incident onto the interface from the left side. We denote as $\alpha$ the whole set of indexes characterising the wave: $\alpha=(q_{||},\omega,j)$. The incident wave is partly reflected with an amplitudes $A_{\alpha\beta}^{L}$ and partly transmitted with an amplitudes $B_{\alpha\beta}^{R}$, where the first index stands for parameters of the incident wave and the second one for the parameters of the departing wave (Fig. 2a).

Whatever model of the phonon transmission across the interface we consider, the law of conservation of energy flow must be fulfilled for the combined mode. In other words, the energy flux to the interface and from the interface for each mode must be equal. Otherwise, the solution in the form of combined mode would not be stationary – the interfacial atoms would have to increase or decrease the amplitude of oscillations over time. The energy 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, we will not consider inelastic phonon scattering, so we can reduce the frequency. Thus, we have:

[TABLE]

This property of the combined mode will be used later.

It is impossible to correctly determine the crystal temperatures using the basis of combined modes since the energy of each combined mode includes the vibrational energy of both crystals. Following the method proposed in [25], we transit from the basis of the combined modes to the basis of ordinary plane waves (Fig. 2b). Herewith, the squares of the reflection and transmission amplitudes $A_{\alpha\beta}^{L,R},B_{\alpha\beta}^{L,R}$ become the coefficients for the equation of matching of the distribution functions of phonons - plane waves. In the discrete case, these equations are written as:

[TABLE]

The physical sense of the matching equations (2) is very simple: phonons flying from the interface $n_{\alpha}^{L\leftarrow},n_{\alpha}^{R\rightarrow}$ consist of those that fell on the interface from the same side and were reflected and of those fallen from the other side and transmitted $n_{\alpha}^{L\rightarrow},n_{\alpha}^{R\leftarrow}$. Since the number of phonons in a given mode is proportional to the energy in a given mode, and the energy is proportional to the square of the amplitude, it is clear that it is the squares of the amplitudes that should be used as coefficients. The rigorous derivation of these equations is a straightforward generalization of the derivation of the matching equations in the one-dimensional case from Ref. [25].

To describe the distribution function of phonons near the boundary, we use the conventional stationary Boltzmann equations in the relaxation time approximation:

[TABLE]

Here $\chi^{L,R}=n^{L,R}-n_{0}^{L,R}$ is a nonequilibrium part of the distribution functions, $n_{0}$ – is equilibrium part, that is Bose-Einstein distribution function.

Here we can introduce the definition of temperatures with phonon occupation numbers:

[TABLE]

In our method, we find $n_{0}^{L},n_{0}^{R}$, and then $T^{L},T^{R}$ is derived with equation (4). The naive method, that leads to the paradox described in the Introduction, defines $T^{L},T^{R}$ by the same formula but with occupation numbers of combined modes $N^{L},N^{R}$. The elimination of the paradox with this changed definition of temperatures for the one-dimensional case is described in detail in Ref. [25].

There is an ambiguity in the division of the distribution functions into the sum of the equilibrium and nonequilibrium parts. To eliminate this ambiguity, we introduce the Chapman-Enskog conditions [27]: the temperature of a nonequilibrium system is defined as the temperature of an equilibrium system with the same energy, which yields

[TABLE]

In order to connect the phonon distribution functions with the heat flux, we write

[TABLE]

Since the heat flux through the boundary is a conserved quantity, the last two equations are equivalent and can be restricted to an equation only for phonons in the left crystal. It is easy to verify that these equations are formally equivalent since the distribution functions obey the matching equations (2) where coefficients satisfy the flow conservation equation (1).

The matching equations for the phonon distribution functions (2) and the Boltzmann equations for phonons on both sides of the interface (3), together with the Enskog-Chapman conditions (5) and the expression for heat flux (6), make up the complete system of equations for the phonon distribution functions on both sides of the interface.

It is easy to verify that the result of calculation of the Kapitza resistance for an imaginary boundary – some plane in a homogeneous crystal, with the system of equations (2, 3, 5, 6) is exactly zero (unlike the AMM calculations). In this case

[TABLE]

For non-equilibrium distribution functions, we make a substitution that satisfies the Boltzmann equations (3):

[TABLE]

And $\chi^{L}=\chi^{R}$ because in a homogeneous media $\tau^{L}=\tau^{R}=\tau$ and $v_{x}^{L}=v_{x}^{R}=v_{x}$.

Heat flux equation (6) is satisfied:

[TABLE]

because the expression in parentheses is exactly the thermal conductivity in the relaxation time approximation $\kappa$. Chapman-Enskog conditions (5) are satisfied, since $v_{x}$ is an odd function of $k_{x}$ hence the integral on $k_{x}$ equals zero. Thus, we have found the non-equilibrium part of the distribution function.

Substituting (7) into (2) we find that $n_{\alpha}^{L\leftarrow}=n_{\beta}^{R\leftarrow},n_{\alpha}^{R\rightarrow}=n_{\beta}^{L\rightarrow}$, and because $\chi^{L}=\chi^{R}$ we find $n_{0}^{L}=n_{0}^{R}$. By definition (4) it meens $T^{L}=T^{R}=T$.

So distribution functions $n^{L}=n^{R}=n_{0}(T)-\tau v_{x}\frac{\partial n_{0}}{\partial T}\frac{q}{\kappa}$ are solutions of the system of equations (2, 3, 5, 6) with the transmission and reflection amplitudes defined by conditions (7).

Since $T^{L}=T^{R}$ there is no temperature jump associated with phonon reflection $\Delta T^{B}=0$. From expression (9) it can be seen that the temperature gradient at the interface is equal to the temperature gradient far from the interface $q/\kappa$, hence $\Delta T^{L}=\Delta T^{R}=0$. The temperature jump on the imaginary boundary vanishes as it should be.

III Model

Now that we derived a set of equations describing phonon kinetics at the interface in a general form, we want to formulate a specific model, i.e. a set of parameters which we want to substitue into equations (2, 3, 5, 6).

We use a Debye model of phonon spectrum. We also introduce the averaged values of the phonon velocities by the formula [36]:

[TABLE]

Here $v^{L,R}_{t}$ is a transversal and $v^{L,R}_{l}$ is a longitudinal speed of sound in the corresponding crystal.

With this in mind, we rewrite the matching equations (2) in a more convenient continuous form and consider the reflection and transmission amplitudes as functions of angles, $\theta^{\prime}$ – angle of incidence, $\theta$ – angle of reflection or transmission, angles are counted from the $x$ axis perpendicular to the interface:

[TABLE]

The ratio of velocity cubes before the second term in the right-hand sides arises as the ratio of the densities of states in the Debye model when transiting to the continuous limit.

We will now choose a set of phonon reflection and transmission amplitudes $A^{L,R}_{\theta\theta^{\prime}},B^{L,R}_{\theta\theta^{\prime}}$. We take as a basis the simple for calculations and often used in the literature DMM [13]. The choice of such a model also justifies the use of averaged velocities (10), since in this model the details of the phonon spectrum are not important, but only the phonon densities of states are important.

In DMM, it is assumed that the amplitudes of the reflected and transmitted waves do not depend on the angle of reflection or transmission, as well as on the angle of incidence of the wave on the boundary:

[TABLE]

However, we want to consider the coherent scattering of phonons at the interface (as, for example, in [35]) and the description of the lattice dynamics with combined modes. It is easy to see that assumption (12) is not compatible with the condition of conservation of energy flow in each combined mode (1). Instead of the assumption of the classical DMM, we introduce the following assumption: the fraction of the energy flux dissipated in a certain direction does not depend on the angle of incidence of the wave on the interface. Energy flux of the incident wave in the direction perpendicular to the interface is proportional to $v\cos\theta^{\prime}$. So this assumption yields

[TABLE]

Let us find the set of amplitudes $A^{L,R},B^{L,R}$.

When substituting the equilibrium distribution function into system (2), the system must become an identity for any set of reflection and transmission amplitudes, since at the same crystal temperatures there is no heat transfer between them, and the phonon distribution functions are equilibrium distribution functions.

[TABLE]

This equation is analogous to the principle of the detailed equilibrium of the classical DMM [10].

For the conservation of flow equation (1), we proceed to the continuous form in the same way as we did for the matching equations:

[TABLE]

We substitute the expression (13) into the equations (14, 15) and get

[TABLE]

The rank of this system of equations is 3. In order to obtain a complete system of equations for determining amplitudes $A^{L,R},B^{L,R}$, we add the condition of the independence of the fraction of energy transmitted in a certain direction from the side of incidence $(L,R)$, same as in DMM. Thus

[TABLE]

The rank of the system of equations (16, 17) is 4, so we can find

[TABLE]

This is the desired model set of transmission and reflection amplitudes.

We will also assume that the phonon relaxation times are independent of the frequency $\tau^{L,R}=const(\omega)$. This assumption is strictly fulfilled for strongly non-perfect [37] crystals, and is thus justified for areas near the boundary. In addition, as it will be seen later, the final answer will not depend on $\tau^{L,R}$.

We will call the system of equations (2, 3, 5, 6) with amplitudes (18) and the Debye model of the phonon spectrum and the constant relaxation times ”kinetic DMM”.

IV Solution

Let us find the solution of kinetic DMM. We first find the distribution functions of the phonons incident on the interface by solving the stationary Boltzmann equation (3). Due to the perturbation of the distribution function of phonons by interfacial scattering, we must consider the non-equilibrium part of the distribution function as a function of the coordinate $\chi^{L,R}(x)$. For the left crystal

[TABLE]

This is an ordinary differential equation for $x$, and general solution of the non-homogeneous equation is the sum of a particular solution and the complementary solution of the associated homogeneous equation, $\chi^{L}=\chi^{L}_{p}+\chi^{L}_{c}$. Complementary solution is:

[TABLE]

where $\chi^{L}_{0}=\chi^{L}(x=0)$. Similarly for the right crystal $\chi^{R}_{c}=\chi^{R}_{0}\exp{(x/v^{R}_{x}\tau^{R})}$.

For incident phonons $v^{L}_{x}>0,v^{R}_{x}<0$, the solution of this type increases infinitely. Thus, it turns out that the complementary solution for the distribution function of incident phonons does not satisfy the boundedness condition, which means that the distribution function of incident phonons is determined only by a particular solution.

Particular solution is

[TABLE]

Here $\left(\frac{dT}{dx}\right)^{L,R}$ – unknown functions of coordinate, since, due to the perturbation of the phonon distribution functions by the interface, the temperature gradients near the interface differ from the gradients in a homogeneous media for the same value of the heat flux. So we have three unknown parameters, the temperature jump at the interface $\Delta T^{B}$, and temperature gradients on both sides of the interface. For these parameters to have the same dimension, instead of gradients we introduce:

[TABLE]

Because the considered problem is linear, we can assume the effective contributions to the temperature jump related to smaller thermal conductivity near the interface are proportional to temperature gradiens, so $\Delta T^{L,R}\sim\tilde{\Delta T^{L,R}}$.

In the case of linear heat transport across the interface, temperature jump at the interface is much smaller than temperatures of media on both sides of the interface. We expand destribution function of phonons on one side of the interface into the Tailor series:

[TABLE]

We substitute the expression (22) into (21), and express the distribution functions of incedent phonons at the interface $(x=0)$ with parameters $\tilde{\Delta T^{L,R}},\Delta T^{B}$:

[TABLE]

Here we also expanded $v_{x}=v\cos\theta^{\prime}$.

Using the matching equations (2), we express the distribution functions of receding phonons with the distribution functions of incident phonons. We will use the matching equations in the form (11). Substituting (24) into (11) with amplitudes (18) we find

[TABLE]

The distribution functions of phonons near the boundary are thus expressed through three unknown parameters: $\Delta\tilde{T}^{L,R},\Delta T^{B}$.

In order to finally find the distribution functions of phonons near the boundary, we substitute (24, 25) into equations (5, 6) and obtain a system of three equations with three unknowns. We use $\chi^{L}=n^{L}-n_{0}(T^{L})$. Here is the result of substitution into equation (6) for phonons in the left crystal:

[TABLE]

Hear $k_{||}$ – stands for the components of wave vector parallel to the interface.The multiplier 3 appears as a result of summation on three branches of the phonon spectrum.

We perform integration over angular components of integral, change the variable of integration from the wave vector to the frequency, and simplify the expression:

[TABLE]

Here $\omega_{D}$ is Debye frequency.

We make a substitution $z=\hbar\omega/k_{B}T$ for the integral. Now

[TABLE]

$T_{D}$ is Debye temperature in that of the crystals, in which the Debye temperature is lower. Let the Debye temperature be lower in the left crystal.

One can see that the right part of the expression (28) is the Debye function. The Debuy temperature dependence naturally arises in the calculation of Kapitza conductance for the Debye model of phonon spectrum. However, later we will need this integral with a different lower limit. So we introduce the notation

[TABLE]

We perform the same substitution for equations (5) and obtain

[TABLE]

Equations (27, 30) are system of three equations with three unknowns $\Delta T^{B},\Delta\tilde{T}^{L,R}$. The solution is:

[TABLE]

and

[TABLE]

Note that $\Delta\tilde{T}^{R}$ has a minus sign.

We have found the phonon distribution function at the interface (24, 25) expressed with three parameters (31, 32). Now we want to investigate the relaxation of the distribution function of phonons disturbed by the interface to undisturbed function. That is a distribution function of phonons for constant heat flow in a homogenous media. So we will find the effective contributions to the temperature jump associated with smaller thermal conductivity near the interface $\Delta T^{L,R}$.

To find the complementary part of the nonequilibrium distribution function at the interface, we use the expression

[TABLE]

Substitute the expression for the complementary (20, 33) and particular (21) parts of the nonequilibrium distribution function in the heat flow expression (6) and find:

[TABLE]

where $\left(\frac{dT}{dx}\right)_{pert}$ is a change in temperature gradient due to perturbation by the interface:

[TABLE]

Because of

[TABLE]

after substituting relation (35) into (36) and integrating, we find

[TABLE]

The calculation of $\Delta T^{R}$ is performed similarly, with the difference that for phonons with frequencies greater than the Debye frequency in the left crystal, coefficients in the matching equation (18) should be changed so that $B^{L}=0$, since such phonons do not pass through the interface. We obtain

[TABLE]

where the second term is the relaxation of the distribution function of those phonons that are not included in the heat transport through the interface. Since $\Delta\tilde{T}^{R}$ is a negative quantity, $\Delta T^{R}$ is positive.

V Results and discussion

We have found all the contributions to temperature jump at the interface (31, 32, 37, 38). Now, the Kapitza conductance is expressed as

[TABLE]

It is interesting that, just like in classic DMM, in kinetic DMM, the Kapitza conductance turns out to be independent of the density of crystals. Also, from the expression (39) it is clear why the averaging of the velocities should be carried out according to the formula (10).

The contribution of the relaxation of high-frequency phonons $2v^{R^{2}}\left(I_{T^{L}_{D}}^{T^{R}_{D}}/I_{0}^{T^{R}_{D}}\right)$ is small for a pair of substances with approximately equal values of the Debye temperature. However, if the Debye temperatures of two substances are significantly different, this contribution can be quite significant, given that the speed of sound is usually higher in the substance where the Debye temperature is higher. Considering the contribution from the relaxation of high-frequency phonons significantly improves the coincidence in the case when the results of kinetic DMM coincide well with the experimental results (Fig. 3).

By removal the relaxation of high-frequency phonons factor (HFP) from the expression (39), we obtain the formula, which differs from the classical DMM formula, when averaging (10) is applied to, only by a numerical factor. The result of calculations using kinetic DMM without HFP is approximately $1.8$ times greater. In many cases, classical DMM tends to overestimate the value of the Kapitza conductance. Kinetic DMM has an even worse agreement with experiment in such cases (Fig. 4).

For a better understanding of the proposed method and model, we also consider a model with a free parameter. In order to introduce a free parameter, we discard a rather arbitrary condition of independence of the fraction of the energy transmitted in a certain direction from the side of phonon incidence (17) and introduce the transmission coefficient $t=|B^{R}|^{2}$. Having performed the calculations for this model, similar to those above, we find

[TABLE]

We can see that the dependence on the density of crystals appeared in the new expression. Thus, it turns out that the independence of the Kapitza conductance (39) on densities was obtained due to the acceptance of assumption (17). Unexpectedly, the magnitudes of temperature jumps associated with smaller thermal conductivity near the interface $\Delta T^{L,R}$ are independent of $t$ and are determined by expressions (32, 37, 38).

Expression (40) equals to (39) at $t=\frac{\rho^{L}v^{L}}{\rho^{R}v^{R}}\frac{v^{L^{2}}}{v^{R^{2}}+v^{L^{2}}}$. We will call this value $t_{DMM}$. Let us express $G^{B}$ in terms of $t_{DMM}$

[TABLE]

It can be seen that the resistance associated with the reflection of phonons from the boundary goes to infinity, or, equivalently, the temperature jump at the boundary vanishes at $t=7/4t_{DMM}$. The calculated dependency on the parameter $t$ is quite similar to the one calculated in Ref. [25], which is $t/(1-t)$. Because in this paper we consider the more complex problem, the value 1 is replaced with a quantity $7/4t_{DMM}$ that characterizes the chosen pair of materials.

One can easily see that introducing the parameter $t$ to classic DMM in the same manner leads to linear dependence of Kapitza conductance on $t$. The non-linear dependence is a feature of the presented method that yields a zero value of Kapitza resistance at the imaginary boundary. This feature is shown in the singular nature of the formulae (40, 41).

The nonlinear dependence of Kapitza conductance $G^{B}$ on the transmission coefficient $t$ (41), is a possible cause of the substantial discrepancy between the theory and an experiment in some cases. For example, at $t=3/4t_{DMM}$, the conductivity decreases by almost two times $G^{B}(3/4t_{DMM})=9/16G^{B}(t_{DMM})$. Figure (5) shows the dependencies of $G(T)$ on the temperature at $t$ which gives the best approximation of the experimental data. It turns out that for most pairs of boundaries a good approximation is achieved by varying $t$ within $60\%\,t_{DMM}$.

The contribution of the relaxation of high-frequency phonons to the calculated value of Kapitza resistance is small at the low temperatures and is growing when the temperature grows. That leads to a slower growth rate of the overall Kapitza conductance with temperature. In some cases, it leads to a very good coincidence with the experimentally measured temperature dependence of Kapitza conductance (Fig. 3, Fig. 5c). However, in some cases it can be seen (Fig. 5a, Fig. 5b) that in the low-temperature region the theoretical curves pass below the experimental points, but at the temperatures close to room temperature the theoretical curves turn out to be, on the contrary, higher. We can offer two possible explanations of this deviation of experimental data from the temperature dependence predicted by the theory. It can be assumed that the DMM approach is not correct, the scattering at the boundary is rather weak, and in this case, as shown by theoretical models [33, 35], the phonon transmission coefficient should decrease with increasing phonon frequency. This would explain the slower increase in the Kapitza conductance with increasing temperature, at which high-frequency oscillations are excited. It is possible, on the contrary, to assume that the DMM assumption about the diffuse nature of scattering is correct, and the discrepancy between theory and experiment is associated with the use of the Debye approximation. In Ref. [41] the Kapitza conductance was studied in the DMM model, but with a more realistic dependence of the density of phonon states on the frequency than that in the Debye model. As a result, it turns out that at low frequencies, Kapitza conductance grows faster with increasing temperature than classical DMM predicts, and at high temperatures, on the contrary, it grows more slowly, which explains the observed discrepancy.

In the framework of kinetic DMM we can not introduce such a concept as ”Radiation limit” [42] since, as it was shown in a previous section, the temperature jump on the imaginary boundary vanishes. However the magnitudes of temperature jumps associated with smaller thermal conductivity near the interface $\Delta T^{L,R}$ are independant of $t$ so we can introduce the ”Relaxation limit” which is kinetic DMM at $t=7/4t_{DMM}$ so the temperature jump associated with the reflection of phonons at the interface vanishes, thus taking into account the only the relaxation contribution. The kinetic DMM, as well as classic DMM, predicts a significantly lower Kapitza conductance value than that measured in the experiment for the pairs of substances with very different sound velocities, such as, for example, diamond - lead, and diamond - gold (Fig. 6). In such cases, experimental data cannot be fitted by variation of the free parameter $t$, because even at ”Relaxation limit” the interfacial thermal resistance associated with smaller thermal conductivity near the interface is still larger than the experimentally measured. We share the common point of view [24, 43] that this discrepancy can be explained only with the help of a model that considers the direct transfer of heat from the metal electrons to the dielectric phonons [44, 45, 46], or a model that takes into account both phonon heat transport and electron-phonon coupling [47, 48].

VI Conclusion

A new method for calculating the Kapitza conductance was proposed. Matching equations for the phonon distribution functions are introduced, which, together with the Boltzmann equations for phonons and some additional equations, form a complete system of equations (2, 3, 5, 6) describing the kinetics of heat transfer across the boundary by phonons. A method of solution to such a system was introduced. It is shown that the use of this method does not lead to a paradox: the solution of the corresponding system for an imaginary boundary, that is, for a plane in a homogeneous crystal, gives zero Kapitza resistance. For a model set of amplitudes based on the classical DMM, an exact analytical solution is obtained. Moreover, the algorithm for solving the equations (2, 3, 5, 6) which can be applied to any set of transmission and reflection amplitudes is found. The dependence of Kapitza conductance on temperature is found. The problem was also solved for a set of amplitudes with a free parameter characterizing the fraction of energy transmitted by a phonon across the boundary. A feature of the method that yields a zero value of Kapitza resistance at the imaginary boundary is shown in the singular character of the formula (40). It turned out that the selection of the corresponding parameter value leads to a good approximation of the experimental data.

Acknowledgements.

The authors thank A. Y. Vul for attention to the work, M. G. Glazov for helpful comments, A. S. Khorolskaya for fruitful discussions. A.P. Meilakhs thanks RSF (grant 18-72-00131) for support.