Cross-plane heat transport in thin films with thicknesses comparable to the phonon mean free paths is of both fundamental and practical interest for applications such as light-emitting diodes and quantum well lasers. However, physical insight is difficult to obtain for the cross-plane geometry due to the challenge of solving the Boltzmann equation in a finite domain. Here, we present a semi-analytical series expansion method to solve the transient, frequency-dependent Boltzmann transport equation that is valid from the diffusive to ballistic transport regimes and rigorously includes the frequency-dependence of phonon properties. Further, our method is more than three orders of magnitude faster than prior numerical methods and provides a simple analytical expression for the thermal conductivity as a function of film thickness. Our result enables a straightforward physical understanding of cross-plane heat conduction in thin films.

In the past two decades, thermal transport in thin solid films with thicknesses from tens of nanometers to micrometers has become a topic of considerable importance.1–4 Such films are used in applications ranging from quantum well lasers to electronic devices.5–7 For example, boundary scattering in these films leads to reduced thermal conductivity that results in the inefficient removal of heat in GaN transistors and light emitting diodes (LEDs).8–10 To address these and other problems, it is first necessary to understand the fundamental physics of heat conduction in micro-scale solid thin films.

Heat transport in thin films with thicknesses comparable to the phonon mean free paths (MFPs) is governed by the Boltzmann transport equation (BTE), which is an integro-differential equation of time, real space, and phase space. Due to its high dimensionality, it is in general very challenging to solve. For transport along the thin film, an analytical solution can be easily derived by assuming that the characteristic length scale of the thermal gradient is much longer than the phonon MFPs. Analytical solutions were derived for electron transport by Fuchs and Sondheimer with partially specular and partially diffuse boundary scattering.11,12 Later, the Fuchs-Sondheimer solutions were extended to phonon thermal transport assuming an average phonon MFP, enabling the calculation of thermal conductivity as a function of the film thickness.13,14 Mazumder and Majumdar used a Monte-Carlo method to study the phonon transport along a silicon thin film including phonon dispersion and polarizations.15 

Heat conduction perpendicular to the thin film (cross-plane direction) is much more challenging. In other fields such as neutron transport and thermal radiation, solutions to the BTE for a slab geometry have been obtained using an invariant embedding method,16,17 an iterative method,18 and an eigenfunction expansion approach.19 For heat conduction, Majumdar numerically solved the gray phonon Boltzmann transport using a discrete-ordinate method by assuming that the two surfaces of the film were black phonon emitters.20 Later, Joshi and Majumdar developed an equation of phonon radiative transfer for both steady-state and transient cases, which showed the correct limiting behavior for both purely ballistic and diffusive transport.21 Chen and Tien applied solutions from radiative heat transfer to calculate the thermal conductivity of a thin film attached to two thermal reservoirs.13 Chen obtained approximate analytical solutions of the BTE to study ballistic phonon transport in the cross-plane direction of superlattices and addressed the inconsistent use of temperature definition at the interfaces.22 Cross-plane heat conduction using a consistent temperature definition was then re-investigated by Chen and Zeng.23,24

Despite these extensive efforts to study transport in thin films based on the BTE, solutions for the cross-plane geometry are still only available with expensive numerical calculations. For example, no analogous Fuchs-Sondheimer formula for the in-plane thermal conductivity exists for the cross-plane direction. Furthermore, most of the previous approaches assumed a single phonon MFP even though recent work has demonstrated that the transport properties of phonons in solids vary widely over the broad thermal spectrum.25,26 Incorporating frequency-dependent phonon properties with these prior numerical methods is extremely computationally expensive.

In this work, we present a semi-analytical solution of the frequency-dependent transient BTE using the method of degenerate kernels, also known as a series expansion method.27 Our approach is valid from the diffusive to ballistic transport regimes, is capable of incorporating a variety of boundary conditions, and is more than three orders of magnitude faster than prior numerical approaches. Further, we obtain a simple closed-form expression for cross-plane thermal conductivity, analogous to the Fuch-Sondheimer formula for the in-plane thermal conductivity, which enables the cross-plane thermal conductivity of a thin film to be easily calculated. Our results can be applied to efficiently solve heat conduction problems in numerous practical geometries such as superlattices and the thin films used in thermoreflectance experiments while rigorously incorporating the full phonon dispersion.

The one-dimensional (1D) frequency-dependent BTE for an isotropic crystal under the relaxation time approximation is given by

gωt+vgμgωx=gωg0(T)τω+Qω4π,
(1)

where gω=ω(fω(x,t,θ)f0(T0)) is the desired deviational energy distribution function, g0(T) is the equilibrium deviational distribution function defined below, Qω is the spectral volumetric heat generation, vg is the phonon group velocity, and τω is the phonon relaxation time. Here, x is the spatial variable, t is the time, ω is the phonon frequency, T is the temperature, and μ=cos(θ) is the directional cosine of the polar angle.

Assuming a small temperature rise, ΔT=TT0, relative to a reference temperature, T0, the equilibrium deviational distribution is proportional to ΔT

g0(T)=14πωD(ω)(fBE(T)fBE(T0))14πCωΔT.
(2)

Here, is the reduced Planck constant, D(ω) is the phonon density of states, fBE is the Bose-Einstein distribution, and Cω=ωD(ω)fBET is the mode specific heat. The volumetric heat capacity is then given by C=0ωmCωdω and the thermal conductivity k=0ωmkωdω, where kω=13CωvωΛω and Λω=τωvω is the phonon MFP.

Both gω and ΔT are unknown. Therefore, to close the problem, energy conservation is used to relate gω to ΔT, given by

0ωm[gω(x,t)τω14πCωτωΔT(x,t)]dωdΩ=0,
(3)

where Ω is the solid angle in spherical coordinates and ωm is the cut-off frequency. Note that summation over phonon branches is implied without an explicit summation sign whenever an integration over phonon frequency is performed.

To solve this equation, we first transform it into an inhomogeneous first-order differential equation by applying a Fourier transform to the time variable, giving

iηg̃ω+vgμdg̃ωdx=g̃ωτω+CωτωΔT̃4π+Q̃ω4π,
(4)

where η is the temporal frequency. This equation has the general solution

g̃ω+(x)=Pωeγωμx+0xCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(xx)dx(μ(0,1]),
(5)
g̃ω(x)=Bωeγωμ(Lx)xLCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(xx)dx(μ[1,0]),
(6)

where γω=(1+iητω)/Λω, L is the distance between the two walls, and Pω and Bω are the unknown coefficients determined by the boundary conditions. Here, g̃+(x) indicates the forward-going phonons and g̃(x) the backward-going phonons. In this work, g̃+(x) is specified at one of the two walls and g̃(x) is specified at the other.

Let us assume that the two boundaries are nonblack but diffuse with wall temperature ΔT1 and ΔT2, respectively. The boundary conditions can be written as

g̃ω+(x=0)=Pω=ϵ1Cω4πΔT1+(1ϵ1)10g̃ω(x=0,μ)dμ,
(7)
g̃ω(x=L)=Bω=ϵ2Cω4πΔT2+(1ϵ2)01g̃ω+(x=L,μ)dμ,
(8)

where ϵ1 and ϵ2 are the emissivities of the hot and cold walls, respectively. When ϵ1=ϵ2=1, the walls are black and we recover Dirichlet boundary conditions. Note that while we assume a thermal spectral distribution for the boundary condition, an arbitrary spectral profile can be specified by replacing the thermal distribution with the desired distribution. Equations (7) and (8) are specific for diffuse boundary scattering; the specular case is presented in Appendix  A.

Applying the boundary conditions to Eqs. (5) and (6), we have

g̃ω+(x)=A1ωCω4πeγωμx+eγωμx0LCωΔT̃(x)+Q̃ω(x)τω4πΛω[DωE1(γω(Lx))+B1ωE1(γωx)]dx+0xCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(xx)dx(μ[0,1]),
(9)
g̃ω(x)=A2ωCω4πeγωμ(Lx)+eγωμ(Lx)0LCωΔT̃(x)+Q̃ω(x)τω4πΛω[DωE1(γωx)+B2ωE1(γω(Lx))]dx+xLCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(xx)dx(μ[0,1]),
(10)

where

A1ω=ϵ1ΔT1+(1ϵ1)ϵ2ΔT2E2(γωL)1(1ϵ1)(1ϵ2)(E2(γωL))2,A2ω=ϵ2ΔT2+(1ϵ2)ϵ1ΔT1E2(γωL)1(1ϵ1)(1ϵ2)(E2(γωL))2,B1ω=1ϵ11(1ϵ1)(1ϵ2)(E2(γωL))2,B2ω=1ϵ21(1ϵ1)(1ϵ2)(E2(γωL))2,Dω=(1ϵ1)(1ϵ2)E2(γωL)1(1ϵ1)(1ϵ2)(E2(γωL))2,

and En(x) is the exponential integral given by28 

En(x)=01μn2exμdμ.
(11)

To close the problem, we plug Eqs. (9) and (10) into Eq. (3) and obtain an integral equation for temperature as

20ωmCωτωdωΔT̃(x̂)=0ωmCωτω[A1ωE2(γ̂ωx̂)+A2ωE2(γ̂ω(1x̂))]dω+010ωmQ̃ω(x)Gω(x̂,x̂)Knωdωdx̂+01ΔT̃(x̂)0ωmCωGω(x̂,x̂)Knωτωdωdx̂,
(12)

where x̂=x/L, Knω=Λω/L is the Knudsen number, γ̂ω=1+iητωKnω and

Gω(x̂,x̂)=E2(γ̂ωx̂)[DωE1(γ̂ω(1x̂))+B1ωE1(γ̂ωx̂)]+E2(γ̂ω(1x̂))[DωE1(γ̂ωx̂)+B1ωE1(γ̂ω(1x̂))]+E1(γ̂ω|x̂x̂|).
(13)

Equation (12) can be written in the form

ΔT(x̂)=f(x̂)+01K(x̂,x̂)ΔT(x̂)dx̂,
(14)

where the kernel function K(x̂,x̂) is given by

K(x̂,x̂)=120ωmCωτωdω0ωmCωGω(x̂,x̂)Knωτωdω
(15)

and the inhomogeneous function f(x̂) is given by

f(x̂)=120ωmCωτωdω0ωmCωτω[A1ωE2(γ̂ωx̂)+A2ωE2(γ̂ω(1x̂))]dω+120ωmCωτωdω010ωmQ̃ω(x)Gω(x̂,x̂)Knωdωdx̂.
(16)

From Eq. (14), we see that the governing equation is a Fredholm integral equation of the second kind. Previously, the gray version of Eq. (12) that assumes average phonon properties has been solved numerically using integral discretization method.28 While this approach does yield the solution, it requires the filling and inversion of a large, dense matrix, an expensive calculation even for the gray case. Considering the full phonon dispersion adds additional integrations to calculate each element of the matrix, dramatically increasing the computational cost. Additionally, care must be taken to account for a singularity point at x̂=x̂ since E1(0).

Here, we solve this equation using the method of degenerate kernels,27 which is much more efficient than the integral discretization method and automatically accounts for the singularity point at x̂=x̂. This method is based on expanding all the functions in Eq. (14) in a Fourier series, then solving for the coefficients of the temperature profile. From the temperature ΔT̃(x̂), all other quantities such as the distribution and heat flux can be obtained.

To apply this method, we first need to expand the inhomogeneous function f(x̂) and kernel K(x̂,x̂) with a Fourier series. This expansion is possible because both f(x̂) and K(x̂,x̂) are continuous and continuously differentiable on the relevant spatial domains of normalized length between [0,1] and [0,1]×[0,1], respectively.27 All the necessary functions can be expanded using a linear combination of sines and cosines; however, a substantial simplification can be obtained by solving a symmetric problem in which the spatial domain is extended to include its mirror image by extending both f(x̂) and K(x̂,x̂) to [−1,1] and [−1,1]× [−1,1]. Because of the symmetry of this domain, all the coefficients for sine functions equal zero and the Fourier series for both functions reduces to a cosine expansion. f(x̂) is then approximated as

f(N)(x̂)=12f0+m=1Nfmcos(mπx̂),
(17)

where fm=201f(x̂)cos(mπx̂)dx̂. The kernel K(x̂,x̂) can be represented by a degenerate double Fourier series, given by

K(N)(x̂,x̂)=14k00+12m=1Nkm0cos(mπx̂)+12n=1Nk0ncos(nπx̂)+m=1Nn=1Nkmncos(mπx̂)cos(nπx̂),
(18)

where

kmn=40101K(x̂,x̂)cos(mπx̂)cos(nπx̂)dx̂dx̂.
(19)

Moreover, the convergence and completeness theorem of cosine functions guarantees that K(N)(x̂,x̂) and f(N)(x̂) converge to K(x̂,x̂) and f(x̂) as N.29 

Inserting Eqs. (17) and (18) into Eq. (14), we then obtain the following integral equation

m=0Nxmcos(mπx̂)=12f0+m=0Nfmcos(mπx̂)+01n=0Nxmcos(nπx̂)[14k00+12m=1Nkm0cos(mπx̂)+12n=1Nk0ncos(nπx̂)+m=1Nn=1Nkmncos(mπx̂)cos(nπx̂)]dx̂,
(20)

where xm are the desired but unknown Fourier coefficients of ΔT̃(x̂).

Using the orthogonality of cos (nπx̂) on [0,1] gives a simpler form of Eq. (20)

m=0Nxmcos(mπx̂)=12f0+m=0Nfmcos(mπx̂)+14m=0Nkm0xm+12m=1Nn=1Nkmnxncos(mπx̂).
(21)

Grouping the terms with the same index m in cosine, a system of linear equations of xm can be obtained as

A¯¯x¯=f¯,
(22)

where x¯ is the vector of unknown coefficient xm and f¯ is the vector of fm in Eq. (17). The matrix A¯¯ contains elements A00=1k004,A0n=12k0n,An0=kn04,Ann=112knn, and Anm=12knm (mn0). Expressions of the elements in A¯¯ can be obtained analytically for the specific kernel here and are given in Appendix  B for steady-state heat conduction with diffuse walls. Since there is no row or column in A¯¯ that is all zeros or a constant multiple of another row or column, it is always guaranteed that A¯¯ is non-singular and its inverse exists.

Solving the matrix system yields xm and thus the temperature ΔT̃(x̂),g̃ω+(x), and g̃ω(x) can be obtained from ΔT̃(x̂) using Eqs. (9) and (10). Finally, the spectral heat is given by

qω(x)=11gωvωμdμ=01gω+vωμdμ01gωvωμdμ,
(23)

thereby closing the problem.

We now summarize the key steps to implement our method. The first step is to determine the appropriate boundary conditions for the problem and compute the constants in Eqs. (5) and (6). Subsequently, the kernel function K(x̂,x̂) and the inhomogenous function f(x̂) can be obtained from Eq. (3), and their Fourier coefficients can be computed using Eqs. (17) and (18). The elements in A¯¯ correspond to the Fourier coefficients of kernel function K(x̂,x̂), and f¯ is a vector of the Fourier coefficients of the inhomogeneous part of Eq. (12). We emphasize that analytic expressions for all of these elements exist and can be obtained; examples of these coefficients for steady heat conduction with non-black, diffuse boundaries are given in Appendix  B. Once A¯¯ and f¯ are obtained, Eq. (22) is solved by standard matrix methods to yield the coefficients xm. Finally, ΔT̃(x̂) is given by m=0Nxmcos(mπx̂).

The primary benefit of our method is the substantial reduction in computational cost compared to the widely used integral discretization approach. Since both K(N)(x̂,x̂) and f(N)(x̂) converge to K(x̂,x̂) and f(x̂) as 1/N2, only a few terms of expansion are required for accurate calculations. In practice, we find that only 20 terms are necessary before the calculation converges, meaning the required matrix is only 20 × 20. Compared to the traditional integral discretization method that requires a matrix on the order of 1000 × 1000 before convergence is achieved, our approach is at least 3 orders of magnitude faster. Further, as we will show in Sec. III, our semi-analytical approach enables a closed-form solution for the cross-plane thermal conductivity of a thin film that is not possible to derive from the integral discretization method.

As an example calculation, we consider steady-state heat conduction between two walls that are either black or non-black. In the former case, both wall emissivities ϵ1 and ϵ2 equal 1 while in the latter case they are set to 0.5. Assuming steady state and no heat generation inside the domain, η = 0, and Qω=0. The Fourier coefficients of K(x̂,x̂) and f¯ for these two specific cases are given in Appendix  B. We perform our calculations for crystalline silicon, using the experimental dispersion in the [100] direction and assuming the crystals are isotropic. The numerical details concerning the dispersion and relaxation times are given in Ref. 30.

We calculate the deviational temperature distribution ΔT(x̂) for different film thickness at different equilibrium temperatures as shown in Fig. 1 while keeping |ΔT1|=|ΔT2|=1 K. When the averaged Knudsen number is small such that Knavg1, the temperature profile remains linear. As thin film thickness decreases such that Knavg1 or 1, we observe a similar temperature slip as discussed in Ref. 25. These calculations take approximately one second to compute on a standard laptop computer. In contrast, the integral discretization method is at least 1000 times slower, requiring on the order of one hour to arrive at the same result.

FIG. 1.

Temperature distribution ΔT(x̂) for a planar slab with black walls (solid lines) and nonblack walls (dashed lines) when (a) KnavgO(102), (b) Knavg1, and (c) KnavgO(102). As Knavg increases, temperature slip at the boundaries grows larger.

FIG. 1.

Temperature distribution ΔT(x̂) for a planar slab with black walls (solid lines) and nonblack walls (dashed lines) when (a) KnavgO(102), (b) Knavg1, and (c) KnavgO(102). As Knavg increases, temperature slip at the boundaries grows larger.

Close modal

In addition to the finite-layer geometry we consider above, the series expansion approach can be readily applied to many other thin film geometries, such as superlattices and the transducer film used in thermoreflectance experiments, by imposing different boundary conditions. Similar large reductions in computational cost can be expected for these cases.

Our semi-analytical approach also allows us to obtain a simple closed form expression for the cross-plane thermal conductivity as a function of film thickness, analogous to the Fuchs-Sondheimer expression for in-plane thermal conductivity. Such a formula allows the cross-plane thermal conductivity to be easily evaluated because the full solution of the BTE is no longer required. To derive this formula, we assume black walls and calculate the spatially averaged spectral heat flux that is integrated over the domain, defined as

01qω(x̂)dx̂=1L0Lqω(x)dx=(ΔT1ΔT22)[13CωvωKnωCωvωKnωE4(1Knω)]+Cωvω2Knω[010x̂ΔT(x̂)E2(|x̂x̂|Knω)dx̂dx̂01x̂1ΔT(x̂)E2(|x̂x̂|Knω)dx̂dx̂].
(24)

Once xm is solved from Eq. (20), we can insert the Fourier series of ΔT(x) into Eq. (24), which leads to

01qω(x̂)dx̂=[(ΔT1ΔT22)13CωvωKnωCωvωKnωE4(1Knω)]+Cωvω2Knωm=1xm[1(1)m]01(Knωμ)2(1+e1Knωμ)1+(Knωμ)2(mπ)2dμ.
(25)

According to Fourier's law, the integrated heat flux is given by

01qωf(x̂)dx̂=13CωvωKnω(ΔT1ΔT2).
(26)

The heat suppression function is defined as the ratio of the BTE and Fourier's heat flux,31 given as

S(Knω,L)=1232E4(1Knω)+32m=1xm[1(1)m]01μ2(1+e1Knωμ)1+(Knωμ)2(mπ)2dμ.
(27)

Note that the suppression function in general not only depends on Knω but also is a function of geometry through xm. The reduced or apparent thermal conductivity at a given domain thickness L is then given by

k(L)=0ωm13CωvωΛωS(Knω,L)dω.
(28)

This formula is analogous to the Fuch-Sondheimer equation for transport along thin films and allows the evaluation of the cross-plane thermal conductivity provided the expansion coefficients xm are known. However, obtaining the expansion coefficients still requires solving the BTE as described in Sec. II B. A more useful result would be a suppression function that depends only on the Knudsen number as is available for in-plane heat conduction with the Fuchs-Sondheimer formula.11,12

To overcome this difficulty, we derive a simplified form of Eq. (27) that is valid under the conditions of most experiments. Note from Fig. 1(a) that for KnavgO(102), the temperature distribution is still linear, allowing us to simplify Eq. (27) by inserting the linear temperature distribution. Doing so leads to a simplified suppression function

Ssimplified(Knω)=1+3Knω[E5(1Knω)14].
(29)

This equation depends only on the Knudsen number and hence can be used to directly evaluate the cross-plane thermal conductivity given the phonon dispersion. This equation is valid provided that the ballistic modes are only low frequency phonons that contribute little to heat capacity, a situation that occurs often in experiment because high frequency phonons have short MFPs, on the order of tens of nanometers, at temperatures exceeding 20 K.

One important observation from Fig. 2(a) is that the exact and simplified suppression functions converge to the same curve at large Knω. Also note that as the slab thickness decreases, the Knudsen number of a phonon with a particular MFP becomes larger. Therefore, in the limit of very small distance between the boundaries, the only important portion of the suppression function is at large values of Knudsen number exceeding Knω=1 because phonons possess a finite minimum MFP. This observation suggests that for practical purposes the simplified suppression can be used even outside the range in which it is strictly valid with good accuracy. This simplification is very desirable because the simplified suppression function only depends on the Knudsen number and thus can be applied without any knowledge of other material properties.

FIG. 2.

(a) Simplified (solid line) and exact (dashed line) suppression function versus Knudsen number. The exact and simplified suppression functions converge to the same curve at large Knω. (b) Example MFP reconstructions for silicon at 300 K using numerically simulated data. Plotted are the analytical MFP distribution (solid line), the numerical apparent thermal conductivities (squares), and the reconstructed MFP distribution by the exact suppression function (circles) and by the simplified suppression function (stars). The x axis corresponds to the MFP for the distributions and to the slab thickness for the thermal conductivity data. Both the exact and simplified suppression functions yield satisfactory MFP reconstruction results.

FIG. 2.

(a) Simplified (solid line) and exact (dashed line) suppression function versus Knudsen number. The exact and simplified suppression functions converge to the same curve at large Knω. (b) Example MFP reconstructions for silicon at 300 K using numerically simulated data. Plotted are the analytical MFP distribution (solid line), the numerical apparent thermal conductivities (squares), and the reconstructed MFP distribution by the exact suppression function (circles) and by the simplified suppression function (stars). The x axis corresponds to the MFP for the distributions and to the slab thickness for the thermal conductivity data. Both the exact and simplified suppression functions yield satisfactory MFP reconstruction results.

Close modal

To investigate the accuracy of this approximation, we perform a reconstruction procedure developed by Minnich31 to recover the MFP spectrum from thermal conductivity data as a function of slab thickness using both exact and simplified suppression functions. We follow the exact procedures of the numerical method as described in Ref. 31. Briefly, we synthesize effective thermal conductivities numerically using Eq. (28). Using these effective thermal conductivities and our knowledge of the suppression function, we use convex optimization to solve for the MFP spectrum. In the exact suppression function case, each slab thickness has its own suppression function given by Eq. (27) while in the simplified case Eq. (29) is used for all slab thicknesses.

As shown in Fig. 2(b), both the simplified and exact suppression functions yield satisfactory results. Even though the smallest thickness we consider here is 50 nm, close to the ballistic regime, the simplified suppression function still gives a decent prediction over the whole MFP spectrum, with a maximum of 15% deviation from the actual MFP spectrum. For practical purposes, this deviation is comparable to uncertainties in experimental measurements and therefore the simplified suppression function can be used as an excellent approximation in the reconstruction procedure. This result demonstrates that length-dependent thermal conductivity measurements like those recently reported for SiGe nanowires32 and graphene ribbons33 can be used to reconstruct the full MFP spectrum rather than only an average MFP. We perform an investigation of our approach for this purpose in Ref. 34.

We have presented a series expansion method to solve the one-dimensional, transient frequency-dependent BTE in a finite domain and demonstrated its capability to describe cross-plane heat conduction in thin films. Our solution is valid from the diffusive to ballistic regimes with a variety of boundary conditions, rigorously includes frequency dependence, and is more than three orders of magnitude faster than prior numerical approaches. We have also developed a simple analytical expression for thermal conductivity, analogous to the Fuchs-Sondheimer equation for in-plane transport, which enables the simple calculation of the cross-plane thermal conductivity as a function of film thickness. Our work provides an efficient method to solve cross-plane heat conduction problems that occurs in numerous situations such as in superlattices and thermoreflectance experiments.

This work was sponsored in part by the National Science Foundation under Grant No. CBET CAREER 1254213, and by Boeing under the Boeing-Caltech Strategic Research and Development Relationship Agreement.

Here, we derive the governing equation for the problem of nonblack, specular boundaries with wall temperatures ΔT1 and ΔT2, respectively. The boundary conditions can be written as

g̃ω+(x=0,μ)=Pω=ϵ1Cω4πΔT1+(1ϵ1)g̃ω(x=0,μ),
(A1)
g̃ω(x=L,μ)=Bω=ϵ2Cω4πΔT2+(1ϵ2)g̃ω+(x=L,μ).
(A2)

Applying the boundary conditions to Eqs. (5) and (6), we have

g̃ω+(x)=F1ΔT1Cω4πeγωμx+(1ϵ1)F2ΔT2Cω4πeγωμ(L+x)+(1ϵ1)F20LCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(x+x)dx+0xCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(xx)dx(μ[0,1]),
(A3)
g̃ω(x)=F2ΔT2Cω4πeγωμ(Lx)+(1ϵ2)F1ΔT1Cω4πeγωμ(2Lx)+(1ϵ2)F10LCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(2Lxx)dx+xLCωΔT̃(x)+Q̃ω(x)τω4πΛωμeγωμ(xx)dx(μ[0,1]),
(A4)

where F1=ϵ1ϵ1+ϵ2ϵ1ϵ2 and F2=ϵ2ϵ1+ϵ2ϵ1ϵ2.

To close the problem, we insert Eqs. (A3) and (A4) into Eq. (3) and nondimensionalize x by L. We then derive an integral equation for temperature for the specular boundary conditions, given by

20ωmCωτωdωΔT̃(x̂)=0ωmCωτωHω(x̂)dω+010ωmQ̃ω(x)Gω(x̂,x̂)Knωdωdx̂+01ΔT̃(x̂)0ωmCωGω(x̂,x̂)Knωτωdωdx̂,
(A5)

where x̂=x/L, Knω=Λω/L is the Knudsen number, γ̂ω=1+iητωKnω and

Hω(x̂)=F1ΔT1E2(γ̂ωx̂)+F2ΔT2E2(γ̂ω(1x̂))+(1ϵ1)F2ΔT2E2(γ̂ω(1+x̂))+(1ϵ2)F1ΔT1E2(γ̂ω(2x̂))
(A6)

and

Gω(x̂,x̂)=(1ϵ1)F2E1(γ̂ω(x̂+x̂))+(1ϵ2)F1E1(γ̂ω(2x̂x̂))+E1(γ̂ω|x̂x̂|).
(A7)

In this case, the inhomogeneous function becomes

f(x̂)=120ωmCωτωdω[0ωmCωτωHω(x̂)dω+010ωmQ̃ω(x)Gω(x̂,x̂)Knωdωdx̂],
(A8)

and the kernel function becomes

K(x̂,x̂)=120ωmCωτωdω0ωmCωGω(x̂,x̂)Knωτωdω.
(A9)

With these results, the problem can be solved by following the same procedures described in Sec. II B are followed to formulate a linear system of equations. The solution of this system then yields the temperature Fourier coefficients.

In this section, we derive the expansion coefficients for a thin film with nonblack, diffuse boundaries. For steady-state heat conduction between two non-black walls as studied in Sec. III, the inhomogeneous function becomes

f(x̂)=120ωmCωτωdω0ωmCωτω[A1ωE2(x̂Knω)+A2ωE2(1x̂Knω)]dω.
(B1)

Its Fourier coefficients in Eq. (17) are then given by

f0=120ωmCωτωdω0ωmCωKnωτω(A1ω+A2ω)[12E3(1Knω)]dω,
(B2)

and

fn=10ωmCωτωdω0ωm01CωτωKnωμ[A1ω+(1)nA2ω]e1Knωμ[(1)nA1ω+A2ω]1+(Knωμ)2(nπ)2,
(B3)

providing the right-hand side of Eq. (22). Under the same assumption of diffuse, non-black walls, the kernel function becomes

K(x̂,x̂)=120ωmCωτωdω0ωmCωGω(x̂,x̂)Knωτωdω,
(B4)

where

Gω(x̂,x̂)=E2(x̂Knω)[DωE1(1x̂Knω)+B1ωE1(x̂Knω)]+E2(1x̂Knω)[DωE1(x̂Knω)+B1ωE1(1x̂Knω)]+E1(|x̂x̂|Knω).
(B5)

Its Fourier coefficients kmn are given by Eq. (18), and can be evaluated as

k00=20ωmCωτωdω0ωmCωKnωτω{2Knω1+2E3(1Knω)+(2Dω+B1ω+B2ω)[12E3(1Knω)12E2(1Knω)+E3(1Knω)E2(1Knω)]}dω,
(B6)

and

km0=20ωmCωτωdω0ωmCωτω01Knωμ[(1)m+1][e1Knωμ1]1+(Knωμ)2(mπ)2dμdω+20ωmCωτωdω0ωmCωτω(Dω+B1ω)[1E2(1Knω)]01Knωμ[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμdω+20ωmCωτωdω0ωmCωτω(Dω+B2ω)[1E2(1Knω)]01Knωμ[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμdω,
(B7)

and

k0n=20ωmCωτωdω0ωmCωτω01Knωμ[(1)n+1][e1Knωμ1]1+(Knωμ)2(nπ)2dμdω+20ωmCωτωdω0ωmCωKnωτω(Dω+B1ω)[12E3(1Knω)]01[1(1)ne1Knωμ]1+(Knωμ)2(nπ)2dμdω+20ωmCωτωdω0ωmCωKnωτω(Dω+B2ω)[12E3(1Knω)]01[(1)ne1Knωμ]1+(Knωμ)2(nπ)2dμdω,
(B8)

and for mn

kmn=20ωmCωτωdω0ωmCωτω01Knωμ{e1Knωμ[(1)m+(1)n][1+(1)m+n}[1+(Knωμ)2(mπ)2][1+(Knωμ)2(nπ)2]dμdω+20ωmCωτωdω0ωmCωDωτω01Knωμ[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[(1)ne1Knωμ]1+(Knωμ)2(nπ)2dμdω+20ωmCωτωdω0ωmCωDωτω01Knωμ[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[1(1)ne1Knωμ]1+(Knωμ)2(nπ)2dμdω+20ωmCωτωdω0ωmCωB1ωτω01Knωμ[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[1(1)ne1Knωμ]1+(Knωμ)2(nπ)2dμdω+20ωmCωτωdω0ωmCωB2ωτω01Knωμ[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[(1)ne1Knωμ]1+(Knωμ)2(nπ)2dμdω,
(B9)

and for m0

kmm=20ωmCωτωdω{0ωmCωτωtan1(mπKnω)mπKnωdω+20ωmCωτω01Knωμ[e1(1)mKnωμ1][1+(Knωμ)2(mπ)2]2dμdω}+20ωmCωτωdω0ωmCωDωτω01Knωμ[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμdω+20ωmCωτωdω0ωmCωDωτω01Knωμ[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμdω+20ωmCωτωdω0ωmCωB1ωτω01Knωμ[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[1(1)me1Knωμ]1+(Knωμ)2(mπ)2dμdω+20ωmCωτωdω0ωmCωB2ωτω01Knωμ[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμ01[(1)me1Knωμ]1+(Knωμ)2(mπ)2dμdω.
(B10)

These equations specify the matrix elements of A¯¯ in Eq. (22). With the linear system specified, the coefficients of the temperature profile xm can be easily obtained by solving a linear system.

1.
D. G.
Cahill
,
P. V.
Braun
,
G.
Chen
,
D. R.
Clarke
,
S.
Fan
,
K. E.
Goodson
,
P.
Keblinski
,
W. P.
King
,
G. D.
Mahan
,
A.
Majumdar
,
H. J.
Maris
,
S. R.
Phillpot
,
E.
Pop
, and
L.
Shi
, “
Nanoscale thermal transport. II. 2003–2012
,”
Appl. Phys. Rev.
1
(
1
),
011305
(
2014
).
2.
Y.
Wang
,
H.
Huang
, and
X.
Ruan
, “
Decomposition of coherent and incoherent phonon conduction in superlattices and random multilayers
,”
Phys. Rev. B
90
,
165406
(
2014
).
3.
X.
Wang
and
B.
Huang
, “
Computational study of in-plane phonon transport in si thin films
,”
Sci. Rep.
4
,
6399
(
2014
).
4.
A. A.
Balandin
and
D. L.
Nika
, “
Phononics in low-dimensional materials
,”
Mater. Today
15
(
6
),
266
275
(
2012
).
5.
I.
Chowdhury
,
R.
Prasher
,
K.
Lofgreen
,
G.
Chrysler
,
S.
Narasimhan
,
R.
Mahajan
,
D.
Koester
,
R.
Alley
, and
R.
Venkatasubramanian
, “
On-chip cooling by superlattice-based thin-film thermoelectrics
,”
Nat. Nanotechnol.
4
(
4
),
235
238
(
2009
).
6.
Z.
Su
,
L.
Huang
,
F.
Liu
,
J. P.
Freedman
,
L. M.
Porter
,
R. F.
Davis
, and
J. A.
Malen
, “
Layer-by-layer thermal conductivities of the Group III nitride films in blue/green light emitting diodes
,”
Appl. Phys. Lett.
100
(
20
),
201106
(
2012
).
7.
A. L.
Moore
and
L.
Shi
, “
Emerging challenges and materials for thermal management of electronics
,”
Mater. Today
17
(
4
),
163
174
(
2014
).
8.
Z.
Yan
,
G.
Liu
,
J. M.
Khan
, and
A. A.
Balandin
, “
Graphene quilts for thermal management of high power gan transistors
,”
Nat. Commun.
3
,
827
(
2012
).
9.
Y. K.
Koh
,
Y.
Cao
,
F.
Liu
,
D. G.
Cahill
, and
D.
Jena
, “
Heat-transport mechanisms in superlattices
,”
Adv. Functional Mater.
19
(
4
),
610
615
(
2009
).
10.
J.
Cho
,
Y.
Li
,
W. E.
Hoke
,
D. H.
Altman
,
M.
Asheghi
, and
K. E.
Goodson
, “
Phonon scattering in strained transition layers for gan heteroepitaxy
,”
Phys. Rev. B
89
,
115301
(
2014
).
11.
F.
Fuchs
, “
The conductivity of thin metallic films according to the electron theory of metals
,”
Proc. Cambridge Philos. Soc.
34
,
100
108
(
1938
).
12.
E. H.
Sondheimer
, “
The mean free path of electrons in metals
,”
Adv. Phys.
1
,
1
42
(
1952
).
13.
G.
Chen
and
C. L.
Tien
, “
Thermal conductivities of quantum well structures
,”
J. Thermophys. Heat Transfer
7
,
311
318
(
1993
).
14.
G.
Chen
, “
Size and interface effects on thermal conductivity of superlattices and periodic thin-film structures
,”
J. Heat Transfer
119
,
220
229
(
1997
).
15.
S.
Mazumder
and
A.
Majumdar
, “
Monte Carlo study of phonon transport in solid thin films including dispersion and polarization
,”
J. Heat Transfer
123
,
749
759
(
2001
).
16.
S.
Chandrasekhar
and
G.
Münch
, “
The Theory of the Fluctuations in Brightness of the Milky way I
,”
Astrophys. J.
112
,
380
(
1950
).
17.
R.
Bellman
and
G. M.
Wing
,
An Introduction to Invariant Imbedding
(
Wiley
,
New York
,
1975
).
18.
K. M.
Case
and
P. F.
Zweifel
,
Linear Transport Theory
(
Addison-Wesley Publishing Company, Inc
,
1967
).
19.
C. C.
Lii
and
M. N.
Ziik
, “
Transient radiation and conduction in an absorbing, emitting, scattering slab with reflective boundaries
,”
Int. J. Heat Mass Transfer
15
(
5
),
1175
1179
(
1972
).
20.
A.
Majumdar
, “
Microscale heat conduction in dielectric thin films
,”
J. Heat Transfer
115
,
7
16
(
1993
).
21.
A. A.
Joshi
and
A.
Majumdar
, “
Transient ballistic and diffusive phonon heat transport in thin films
,”
J. Appl. Phys.
74
(
1
),
31
39
(
1993
).
22.
G.
Chen
, “
Thermal conductivity and ballistic-phonon transport in the cross-plane direction of superlattices
,”
Phys. Rev. B
57
,
14958
14973
(
1998
).
23.
G.
Chen
and
T. F.
Zeng
, “
Nonequilibrium phonon and electron transport in heterostructures and superlattices
,”
Microscale Thermophys. Eng.
5
,
71
88
(
2001
).
24.
T. F.
Zeng
and
G.
Chen
, “
Phonon heat conduction in thin films: Impact of thermal boundary resistance and internal heat generation
,”
J. Heat Transfer
123
,
340
347
(
2001
).
25.
K.
Esfarjani
,
G.
Chen
, and
H. T.
Stokes
, “
Heat transport in silicon from first-principles calculations
,”
Phys. Rev. B
84
(
8
),
085204
(
2011
).
26.
A.
Ward
and
D. A.
Broido
, “
Intrinsic phonon relaxation times from first-principles studies of the thermal conductivities of si and ge
,”
Phys. Rev. B
81
,
085205
(
2010
).
27.
A. D.
Polyanin
and
A. V.
Manzhirov
,
Handbook of Integral Equations
(
Chapman & Hall/CRC
,
2008
).
28.
G.
Chen
,
Nanoscale Energy Transport and Conversion
(
Oxford University Press
,
New York
,
2005
).
29.
H. F.
Weinberger
,
A First Course in Partial Differential Equations with Complex Variables and Transform Methods
(
Dover Publications, Inc.
,
1995
).
30.
A. J.
Minnich
,
G.
Chen
,
S.
Mansoor
, and
B. S.
Yilbas
, “
Quasiballistic heat transfer studied using the frequency-dependent Boltzmann transport equation
,”
Phys. Rev. B
84
,
235207
(
2011
).
31.
A. J.
Minnich
, “
Determining phonon mean free paths from observations of quasiballistic thermal transport
,”
Phys. Rev. Lett.
109
,
205901
(
2012
).
32.
T.-K.
Hsiao
,
H.-K.
Chang
,
S.-C.
Liou
,
M.-W.
Chu
,
S.-C.
Lee
, and
C.-W.
Chang
, “
Observation of room-temperature ballistic thermal conduction persisting over 8.3 μm in sige nanowires
,”
Nat. Nanotechnol.
8
,
534
538
(
2013
).
33.
X.
Xu
,
L. F. C.
Pereira
,
Y.
Wang
,
J.
Wu
,
K.
Zhang
,
X.
Zhao
,
S.
Bae
,
C. T.
Bui
,
R.
Xie
,
J. T. T. L.
Thong
,
B. H.
Hong
,
K. P.
Loh
,
D.
Donadio
,
B.
Li
, and
B.
Ozyilmaz
, “
Length-dependent thermal conductivity in suspended single-layer graphene
,”
Nat. Commun.
5
,
3689
(
2014
).
34.
H.
Zhang
,
C.
Hua
,
D.
Ding
, and
A.
Minnich
, “
Length dependent thermal conductivity measurements yield phonon mean free path spectra in nanostructures
,”
Sci. Rep.
5
,
9121
(
2015
).