PHYSICAL REVIEW B
93
, 035314 (2016)
Role of thermalizing and nonthermalizing walls in phonon heat conduction along thin films
Navaneetha K. Ravichandran and Austin J. Minnich
*
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, California 91125, USA
(Received 5 November 2015; revised manuscript received 31 December 2015; published 28 January 2016)
Phonon boundary scattering is typically treated using the Fuchs-Sondheimer theory, which assumes that
phonons are thermalized to the local temperature at the boundary. However, whether such a thermalization
process actually occurs and its effect on thermal transport remains unclear. Here we examine thermal transport
along thin films with both thermalizing and nonthermalizing walls by solving the spectral Boltzmann transport
equation for steady state and transient transport. We find that in steady state, the thermal transport is governed by
the Fuchs-Sondheimer theory and is insensitive to whether the boundaries are thermalizing or not. In contrast,
under transient conditions, the thermal decay rates are significantly different for thermalizing and nonthermalizing
walls. We also show that, for transient transport, the thermalizing boundary condition is unphysical due to violation
of heat flux conservation at the boundaries. Our results provide insights into the boundary scattering process of
thermal phonons over a range of heating length scales that are useful for interpreting thermal measurements on
nanostructures.
DOI:
10.1103/PhysRevB.93.035314
I. INTRODUCTION
Engineering the thermal conductivity of nanoscale mate-
rials has been a topic of considerable research interest over
the past two decades [
1
]. While applications such as GaN
transistors [
2
,
3
] and light emitting diodes (LEDs) [
4
] require
high thermal conductivity substrates to dissipate heat, the
performance of thermoelectric and thermal insulation devices
can be significantly enhanced by reducing their thermal
conductivity [
5
,
6
]. In many of these applications, phonon
boundary scattering is the dominant resistance to heat flow,
making the detailed understanding of this process essential for
advancing applications.
Phonon boundary scattering has been studied extensively
both theoretically and experimentally. The thermal conduc-
tivity reduction due to boundary scattering of phonons is
conventionally treated using the Fuchs-Sondheimer theory,
which was first derived for electron boundary scattering
independently by Fuchs [
7
] and Reuter and Sondheimer [
8
]
and was later extended to phonon boundary scattering in
several works [
9
–
11
]. Fuchs-Sondheimer theory is widely used
to interpret experiments but makes an important assumption
that the diffusely scattered part of the phonon spectrum at a
partially specular wall is at a local thermal equilibrium with the
wall—the thermalizing boundary condition. The thermalizing
boundary condition is also a key assumption in the diffuse
boundary scattering limit of Casimir’s theory [
12
].
Several computational works [
11
,
13
–
16
] have studied the
reduction in thermal conductivity due to phonon boundary
scattering in nanostructures by solving the phonon Boltzmann
transport equation (BTE). These works have considered either
thermalizing or nonthermalizing boundaries but have never
compared the effect of these two different boundary condi-
tions on the thermal conductivity of nanostructures. Several
experimental works have also studied the reduction in thermal
conductivity of nanomaterials such as nanowires [
17
–
19
], thin
films [
10
,
20
,
21
], and nanopatterned structures [
22
] due to
phonon boundary scattering. These works have used the Fuchs-
*
aminnich@caltech.edu
Sondheimer theory to interpret their measurements. However,
it is not clear whether the assumptions made in the Fuchs-
Sondheimer theory are necessarily applicable for these exper-
iments. In fact, an analysis of the effect of the key assumption
made in the Fuchs-Sondheimer theory, that the walls are
thermalizing, has never been investigated due to the challenges
involved in solving the BTE for nonthermalizing walls.
Here, we examine the role of thermalizing and nonthermal-
izing walls in heat conduction along thin films by solving the
spectral phonon BTE for a suspended thin film under steady
state and transient transport conditions. We find that steady
state transport is insensitive to whether phonons are thermal-
ized or not at the boundaries and that the Fuchs-Sondheimer
theory accurately describes thermal transport along the thin
film. In the case of transient transport, we find that the
decay rates of the initial temperature distribution [defined by
γ
=
4
π
2
k
eff
/
(
Cλ
2
), where
C
is the volumetric specific heat
of the material and
k
eff
is the effective thermal conductivity
of the thin film at a heating length scale
λ/
2] are significantly
different for thermalizing and nonthermalizing walls, and that
the Fuchs-Sondheimer theory accurately predicts the thermal
conductivity only when the thermal transport is diffusive.
Moreover, under transient transport conditions, we find that
phonons cannot undergo thermalization at the boundaries in
general due to the violation of heat flux conservation. Our
results provide insights into the boundary scattering process
of thermal phonons that are useful for interpreting thermal
measurements on nanostructures.
II. MODELING
A. Boltzmann transport equation
We begin our analysis by considering the two-dimensional
spectral transient BTE under the relaxation time approxima-
tion for an isotropic crystal, given by
∂g
ω
∂t
+
μv
ω
∂g
ω
∂z
+
v
ω
1
−
μ
2
cos
φ
∂g
ω
∂x
=−
g
ω
−
g
o
(
T
)
τ
ω
+
Q
ω
4
π
.
(1)
2469-9950/2016/93(3)/035314(10)
035314-1
©2016 American Physical Society
NAVANEETHA K. RAVICHANDRAN AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
93
, 035314 (2016)
Here,
g
ω
is the phonon energy distribution function,
ω
is the
phonon frequency,
v
ω
is the phonon group velocity,
τ
ω
is the
phonon relaxation time,
x
and
z
are the spatial coordinates,
t
is
the time variable,
g
0
(
T
) is the equilibrium phonon distribution
function at a deviational temperature
T
=
T
0
+
T
from
an equilibrium temperature
T
0
,μ
is the direction cosine,
φ
is the azimuthal angle, and
Q
ω
is the rate of volumetric
heat generation for each phonon mode. As the in-plane (
x
)
direction is infinite in extent, we require boundary conditions
only for the cross-plane (
z
) direction. In the traditional
Fuchs-Sondheimer problem, the boundary conditions enforce
that the diffusely scattered phonons are thermalized while
also allowing some phonons to be specularly reflected. Here,
we generalize these boundary conditions to allow for the
possibility of both partial thermalization and partial specularity
as
for
μ
∈
(
0
,
1
]
,
g
+
ω
(
0
,μ,φ
)
=
p
ω
g
−
ω
(
0
,
−
μ,φ
)
+
(
1
−
p
ω
)
σ
ω
C
ω
T
(
z
=
0
)
4
π
−
(
1
−
σ
ω
)
π
2
π
0
0
−
1
g
−
ω
(0
,μ
,φ
)
μ
dμ
dφ
,
for
μ
∈
[
−
1
,
0
)
,
g
−
ω
(
d,μ,φ
)
=
p
ω
g
+
ω
(
d,
−
μ,φ
)
+
(
1
−
p
ω
)
σ
ω
C
ω
T
(
z
=
d
)
4
π
+
(
1
−
σ
ω
)
π
2
π
0
1
0
g
+
ω
d,μ
,φ
μ
dμ
dφ
,
(2)
where
d
is the thickness in the cross-plane direction,
g
+
ω
(0
,μ,φ
) is the phonon distribution leaving the cross-
plane wall at
z
=
0
,g
−
ω
(0
,μ,φ
) is the phonon distribution
approaching the cross-plane wall at
z
=
0
,g
+
ω
(
d,μ,φ
)isthe
phonon distribution approaching the cross-plane wall at
z
=
d,g
−
ω
(
d,μ,φ
) is the phonon distribution leaving the cross-
plane wall at
z
=
d,C
ω
is the specific heat of a phonon mode
with frequency
ω
, and
p
ω
and
σ
ω
are the phonon specularity
parameter and the thermalization parameter for the thin film
walls, respectively. The specularity parameter represents the
fraction of specularly scattered phonons at the thin film walls
and the thermalization parameter represents the fraction of the
phonon distribution that is absorbed and reemitted at the local
equilibrium temperature of the thin film walls. For simplicity,
we ignore mode conversion for nonthermalizing boundary
condition in our analysis. The simulation domain and the
boundary conditions [Eq. (
2
)] are pictorially represented in
Fig.
1
.
The unknown quantities in this problem are the phonon
distribution function [
g
ω
(
t,x,z,μ,φ
)] and the deviational
temperature distribution [
T
(
t,x,z
)]. They are related to each
other through the energy conservation requirement,
ω
m
ω
=
0
1
μ
=−
1
2
π
φ
=
0
g
ω
τ
ω
−
1
4
π
C
ω
τ
ω
T
dφdμdω
=
0
.
(3)
Due to the high dimensionality of the BTE, analytical or
semianalytical solutions are only available in the literature
for either semi-infinite domains [
23
–
25
] or domains with
simple boundary and transport conditions [
26
] or with sev-
eral approximations [
27
]. For nanostructures with physically
realistic boundaries, several numerical solutions of the BTE
have been reported [
11
,
16
,
28
]. However, computationally
efficient analytical or semianalytical solutions for the in-plane
heat conduction along even simple unpatterned films [
10
,
20
]
are unavailable. To overcome this problem, we solve the
BTE analytically for steady state transport (Sec.
II B
) and
semianalytically for transient transport along thin films in the
transient grating (TG) experiment [
10
,
20
] (Sec.
II C
).
B. Steady state heat conduction in thin films
In this section, we extend the Fuchs-Sondheimer relation
for thermal conductivity suppression due to phonon boundary
scattering to the general boundary conditions described in
Eq. (
2
). To simulate steady state transport,
Q
ω
issetto0inthe
BTE [Eq. (
1
)]. Furthermore, we assume that a one-dimensional
Nonthermalizing Wall (g
diff
≠ g
0
)
Thermalizing Wall (g
diff
= g
0
)
g
in
g
in
g
spec
g
spec
g
diff
g
diff
z
x
d
FIG. 1. Spatial distribution of the temperature profile and picto-
rial representation of the boundary conditions [Eq. (
2
)] used in this
work. The thin film is assumed to be infinite in extent along the
in-plane (
x
) direction and has a finite thickness (
d
) in the cross-plane
(
z
) direction. For steady state transport calculations, the temperature
gradient exists only in the in-plane (
x
) direction. For transient
transport, the initial temperature distribution is an in-plane sinusoidal
distribution with a period
λ
, but can develop a cross-plane temperature
gradient at later times. In the case of the nonthermalizing wall, the
diffusely reflected component [
g
diff
] of the distribution function is
isotropic but away from local thermal equilibrium at the boundary
while for the thermalizing wall, [
g
diff
] is equal to the local equilibrium
distribution function [
g
0
]. For both thermalizing and nonthermalizing
boundary conditions, the specular reflection component [
g
spec
]hasits
direction reversed compared to the incoming distribution [
g
in
]and
is also away from the local thermal equilibrium. In our work, we
consider either thermalizing or nonthermalizing boundary conditions
for both walls of the thin films at a time.
035314-2
ROLE OF THERMALIZING AND NONTHERMALIZING . . .
PHYSICAL REVIEW B
93
, 035314 (2016)
temperature gradient exists along the thin film (Fig.
1
)
and
∂g
ω
∂x
≈
∂g
0
ω
∂x
. These assumptions are consistent with the
conditions under which typical steady state thermal transport
measurements are conducted on nanostructures [
19
,
29
,
30
].
Under these assumptions, the BTE is simplified as
v
ω
μ
∂g
ω
∂z
+
v
ω
1
−
μ
2
cos
φ
∂g
0
ω
∂x
=−
g
ω
−
g
0
ω
τ
ω
.
(4)
For steady state transport, it is convenient to solve the
BTE in terms of the deviation from equilibrium distribution
(
̄
g
ω
=
g
ω
−
g
0
ω
(
T
(
x
))). In this case, the BTE transforms
into
∂
̄
g
ω
∂z
+
̄
g
ω
μ
ω
=−
cos
φ
1
−
μ
2
μ
∂g
0
ω
∂x
.
(5)
The boundary conditions [Eq. (
2
)] for
̄
g
ω
now become
for
μ
∈
(
0
,
1
]
,
̄
g
+
ω
(
0
,μ,φ
)
=
p
ω
̄
g
−
ω
(
0
,
−
μ,φ
)
−
(
1
−
p
ω
)(
1
−
σ
ω
)
π
2
π
0
0
−
1
̄
g
−
ω
(0
,μ
,φ
)
μ
dμ
dφ
,
for
μ
∈
[
−
1
,
0
)
,
̄
g
−
ω
(
d,μ,φ
)
=
p
ω
̄
g
+
ω
(
d,
−
μ,φ
)
+
(
1
−
p
ω
)(
1
−
σ
ω
)
π
2
π
0
1
0
̄
g
+
ω
(
d,μ
,φ
)
μ
dμ
dφ
.
(6)
The general solution of the BTE [Eq. (
5
)] along with the boundary conditions [Eq. (
6
)] is given by
̄
g
+
ω
(
z,μ,φ
)
=−
ω
cos
φ
1
−
μ
2
∂g
0
ω
∂x
1
−
exp
−
z
μ
ω
(
1
−
p
ω
)
1
−
p
ω
exp
−
d
μ
ω
+
(
1
−
p
ω
)(
1
−
σ
ω
)
A
+
ω
+
p
ω
exp
−
d
μ
ω
A
−
ω
1
−
p
2
ω
exp
−
2
d
μ
ω
exp
−
z
μ
ω
I
,
̄
g
−
ω
(
z,
−
μ,φ
)
=−
ω
cos
φ
1
−
μ
2
∂g
0
ω
∂x
1
−
exp
−
(
d
−
z
)
μ
ω
(
1
−
p
ω
)
1
−
p
ω
exp
−
d
μ
ω
+
(
1
−
p
ω
)(
1
−
σ
ω
)
A
−
ω
+
p
ω
exp
−
d
μ
ω
A
+
ω
1
−
p
2
ω
exp
−
2
d
μ
ω
exp
−
(
d
−
z
)
μ
ω
II
(7)
for
μ
∈
(0
,
1]. Here, the terms
A
+
ω
and
A
−
ω
only depend on phonon frequency. In particular, they are independent of the angular
coordinates
μ
and
φ
. The derivation of the final expressions for
̄
g
+
ω
(
z,μ,φ
) and
̄
g
−
ω
(
z,
−
μ,φ
)[Eq.(
7
)] is shown in Sec. I of the
Supplemental Material [
31
]. The expression for the in-plane (
x
direction) spectral heat flux is given by
q
x,ω
=
1
d
d
z
=
0
1
μ
=−
1
2
π
φ
=
0
v
x
̄
g
ω
D
(
ω
)
4
π
dφdμdz
=−
1
3
C
ω
v
ω
ω
∂T
∂x
1
−
3
(
1
−
p
ω
)
ω
2
d
1
0
(
μ
−
μ
3
)
1
−
exp
−
d
μ
ω
1
−
p
ω
exp
−
d
μ
ω
dμ
(8)
since the diffuse contributions to the distribution functions
̄
g
+
ω
(
z,μ,φ
) and
̄
g
−
ω
(
z,
−
μ,φ
) [terms I and II in Eq. (
7
)] are independent
of the azimuthal angle
φ
and integrate out to 0. Comparing Eq. (
8
) with the expression for heat flux from the Fourier’s law, the
spectral effective thermal conductivity of the thin film is obtained as a product of the bulk spectral thermal conductivity and the
well-known Fuchs-Sondheimer reduction factor due to phonon boundary scattering given by
k
ω,
eff
(
d
)
=
1
3
C
ω
v
ω
ω
k
ω,
bulk
1
−
3
(
1
−
p
ω
)
ω
2
d
1
0
(
μ
−
μ
3
)
1
−
exp
−
d
μ
ω
1
−
p
ω
exp
−
d
μ
ω
dμ
Fuchs-Sondheimer reduction factor
F
(
ω
d
)
.
(9)
It is interesting to observe from Eq. (
9
) that the spectral effective thermal conductivity is independent of the thermalization
parameter
σ
ω
even though a general boundary condition [Eq. (
2
)] has been used in this derivation. Thus, the steady state thermal
conductivity suppression due to boundary scattering is only influenced by the relative extent of specular and diffuse scattering
(parametrized by the specularity parameter
p
ω
) and does not depend on the type of diffuse scattering process (parametrized by
the thermalization parameter
σ
ω
). We explicitly demonstrate this result using numerical simulations in Sec.
III A
.
035314-3
NAVANEETHA K. RAVICHANDRAN AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
93
, 035314 (2016)
C. Transient heat conduction in thin films
In this section, we solve the BTE [Eq. (
1
)] for transient thermal transport along a thin film. The initial temperature profile
considered in this work is identical to that of the TG experiment, which has been used extensively to study heat conduction in
suspended thin films [
10
,
20
]. In the TG experiment, the thermal transport properties of the sample are obtained by observing the
transient decay of a one-dimensional impulsive sinusoidal temperature grating on the sample at different grating periods. In the
large grating period limit of heat diffusion, the temporal decay is a single exponential. Since the initial temperature distribution
is an infinite one-dimensional sinusoid in the
x
direction, the temperature distribution remains spatially sinusoidal at all later
times. Therefore, each wave vector
q
in the spatially Fourier transformed BTE directly corresponds to a unique grating period
λ
=
2
π/q
. Unlike in the steady state case, here we solve for the absolute phonon distribution
g
ω
rather than the deviation
̄
g
ω
=
g
ω
−
g
0
ω
. Furthermore, the BTE is solved in the frequency domain (
η
) by Fourier transforming Eq. (
1
) in the time variable
t
. With these transformations, the BTE reduces to
iηG
ω
+
μv
ω
∂G
ω
∂z
+
iqv
ω
1
−
μ
2
cos
φG
ω
=−
G
ω
τ
ω
+
1
4
π
C
ω
τ
ω
̄
T
+
̄
Q
ω
4
π
,
(10)
where the substitution
G
0
(
T
)
=
1
4
π
C
ω
̄
T
has been made and
G
ω
represents the spatial (in-plane axis) and temporal Fourier
transform of absolute phonon energy distribution function
g
ω
.
The outline of the solution methodology for Eq. (
10
) is as follows. The general solution is given by
for
μ
∈
(0
,
1]
,G
+
ω
(
z,μ,φ
)
=
G
+
ω
(
0
,μ,φ
)
exp
−
γ
FS
μφ
μ
ω
z
+
exp
−
γ
FS
μφ
μ
ω
z
4
πμ
ω
z
0
(
C
ω
̄
T
+
̄
Q
ω
τ
ω
)exp
γ
FS
μφ
μ
ω
z
dz
,
for
μ
∈
[
−
1
,
0
)
,G
−
ω
(
z,μ,φ
)
=
G
−
ω
(
d,μ,φ
)
exp
γ
FS
μφ
μ
ω
(
d
−
z
)
−
exp
−
γ
FS
μφ
μ
ω
z
4
πμ
ω
d
z
(
C
ω
̄
T
+
̄
Q
ω
τ
ω
)exp
γ
FS
μφ
μ
ω
z
dz
,
where
γ
FS
μφ
=
(
1
+
iητ
ω
)
+
i
ω
q
1
−
μ
2
cos
φ.
(11)
Here,
G
+
ω
(0
,μ,φ
) and
G
−
ω
(
d,μ,φ
) are determined by solving the boundary conditions [Eq. (
2
)] with the following procedure.
First, the angular integrals in the boundary conditions are discretized using Gauss quadrature with weights
w
μ
i
and
w
φ
j
for
the variables
μ
and
φ
respectively, which results in the following set of linear equations in the variables
G
+
ω
(0
,μ
i
,φ
j
) and
G
−
ω
(
d,
−
μ
i
,φ
j
) for every
{
μ
i
,φ
j
}∈
(0
,
1]
×
[0
,
2
π
] doublet from the discretization:
G
+
ω
(0
,μ
i
,φ
j
)
=
p
ω
G
−
ω
(
d,
−
μ
i
,φ
j
)exp
−
γ
FS
ij
μ
i
ω
d
+
p
ω
4
πμ
i
ω
d
0
(
C
ω
̄
T
+
̄
̃
Q
ω
τ
ω
)exp
−
γ
FS
ij
μ
i
ω
z
dz
+
(
1
−
p
ω
)
σ
ω
C
ω
̄
T
(
z
=
0
)
4
π
+
(
1
−
σ
ω
)
π
i
j
G
−
ω
(
d,
−
μ
i
,φ
j
)exp
−
γ
FS
i
j
μ
i
ω
d
μ
i
w
μ
i
w
φ
j
+
(
1
−
σ
ω
)
4
π
2
ω
i
j
d
0
(
C
ω
̄
T
+
̄
̃
Q
ω
τ
ω
)exp
−
γ
FS
i
j
μ
i
ω
z
dz
w
μ
i
w
φ
j
,
G
−
ω
(
d,
−
μ
i
,φ
j
)
=
p
ω
G
+
ω
(0
,μ
i
,φ
j
)exp
−
γ
FS
ij
μ
i
ω
d
+
p
ω
4
πμ
i
ω
d
0
(
C
ω
̄
T
+
̄
̃
Q
ω
τ
ω
)exp
−
γ
FS
ij
μ
i
ω
(
d
−
z
)
dz
+
(
1
−
p
ω
)
σ
ω
C
ω
̄
T
(
z
=
d
)
4
π
+
(
1
−
σ
ω
)
π
i
j
G
+
ω
(0
,μ
i
,φ
j
)exp
−
γ
FS
i
j
μ
i
ω
d
μ
i
w
μ
i
w
φ
j
+
(
1
−
σ
ω
)
4
π
2
ω
i
j
d
0
(
C
ω
̄
T
+
̄
̃
Q
ω
τ
ω
)exp
−
γ
FS
i
j
μ
i
ω
(
d
−
z
)
dz
w
μ
i
w
φ
j
.
(12)
To obtain Eq. (
12
), we have substituted the general BTE solution into the boundary conditions to eliminate
G
−
ω
(0
,μ,φ
) and
G
+
ω
(
d,μ,φ
). Therefore, the only unknowns in the set of linear equations [Eq. (
12
)] are
G
+
ω
(0
,μ,φ
) and
G
−
ω
(
d,μ,φ
). By
bringing the terms containing
G
+
ω
(0
,μ,φ
) and
G
−
ω
(
d,μ,φ
) to the left-hand side, Eq. (
12
) can be written in a concise matrix
form:
U
+
kk
U
−
kk
D
+
kk
D
−
kk
G
+
ω
(0
,μ
i
,φ
j
)
G
−
ω
(
d,
−
μ
i
,φ
j
)
=
̄
̃
c
+
ω
(0
,μ
i
,φ
j
)
̄
̃
c
−
ω
(
d,μ
i
,φ
j
)
,
(13)
035314-4
ROLE OF THERMALIZING AND NONTHERMALIZING . . .
PHYSICAL REVIEW B
93
, 035314 (2016)
where,
̄
̃
c
+
ω
(0
,μ
i
,φ
j
) and
̄
̃
c
−
ω
(
d,μ
i
,φ
j
) are analytical functions of the unknown temperature distribution function
̄
T
obtained
from the right-hand side of Eq. (
12
). The solution to this set of linear equations can be represented as
G
+
ω
(0
,μ
i
,φ
j
)
G
−
ω
(
d,
−
μ
i
,φ
j
)
=
T
+
kk
T
−
kk
B
+
kk
B
−
kk
̄
̃
c
+
ω
(0
,μ
i
,φ
j
)
̄
̃
c
−
ω
(
d,μ
i
,φ
j
)
,
(14)
where
k
is the index which represents the doublet
{
μ
i
,φ
j
}
. The details of the simplification of the boundary conditions and the
evaluation of
T
+
kk
,T
−
kk
,B
+
kk
,B
−
kk
,
̄
̃
c
+
ω
(0
,μ
i
,φ
j
), and
̄
̃
c
−
ω
(
d,μ
i
,φ
j
) are described in Sec. II A of the Supplemental Material [
31
].
To close the problem, the expressions for
G
+
ω
(
z,μ,φ
) and
G
−
ω
(
z,μ,φ
)[Eq.(
11
)] and the boundary conditions [Eq. (
14
)] are
substituted into the energy conservation equation [Eq. (
3
)] and an integral equation in the variable
z
for
̄
T
(
z
) at each
η
and
q
is obtained, which has the form
̄
T
(
z
)
=
h
(
z
)
+
f
(
z
)
+
d
0
[
K
(
z
,z
)
̄
T
(
z
)]
dz
,
(15)
where the functional form of the inhomogeneous parts
f
(
z
)
,h
(
z
) and the kernel
K
(
z
,z
) are described in Sec. II B of the
Supplemental Material [
31
]. This integral equation [Eq. (
15
)] is then solved using the method of degenerate kernels for each
η
and
q
to obtain the frequency domain solution
̄
T
(
z
) for every
η
and
q
. The details of the degenerate kernel calculations
are described in Sec. II C of the Supplemental Material [
31
]. Finally, the solution
̄
T
(
z
) is substituted into Eq. (
11
) to obtain
expressions for
G
ω
(
z,μ,φ
) and also the thickness-averaged in-plane heat flux
j
x,ω
given by
j
x,ω
=
1
4
πd
d
0
2
π
0
1
−
1
G
ω
v
ω
1
−
μ
2
cos
φdμdφdz
=
1
4
π
ij
⎛
⎝
μ
i
Kn
d
ω
γ
FS
ij
i
j
[(
T
+
kk
+
B
+
kk
)
̄
̃
c
+
ω
(0
,μ
i
,φ
j
)
+
(
T
−
kk
+
B
−
kk
)
̄
̃
c
−
ω
(
d,μ
i
,φ
j
)]
1
−
exp
−
γ
FS
ij
μ
i
Kn
d
ω
+
2
4
πγ
FS
ij
C
ω
t
0
2
+
̄
Q
ω
τ
ω
1
−
μ
i
Kn
d
ω
γ
FS
ij
1
−
exp
−
γ
FS
ij
μ
i
Kn
d
ω
−
C
ω
1
−
exp
−
γ
FS
ij
μ
i
Kn
d
ω
4
πμ
i
Kn
d
ω
N
m
=
1
t
m
1
+
(
−
1
)
m
m
2
π
2
+
γ
FS
ij
μ
i
Kn
d
ω
2
⎞
⎟
⎠
v
ω
1
−
μ
2
i
w
μ
i
w
φ
j
,
(16)
where the
t
i
’s are the Fourier coefficients for the expansion
of
̄
T
in the cross-plane (
z
) direction and Kn
d
ω
=
ω
/d
is
the Knudsen number. The conventional approach to describe
the thermal transport properties of the thin film is to compare
the expression for heat flux from the BTE solution with that
expected for heat diffusion, as was done in Eq. (
8
)forthe
steady state Fuchs-Sondheimer theory. However, in practice,
Eq. (
16
) is not easily reduced into the form of Fourier’s law. To
overcome this problem, we adopt the following strategy. The
solution of the Fourier heat equation to a one-dimensional
heat conduction with an instantaneous spatially sinusoidal
heat source is a simple exponential decay
T
(
t,x
=
0)
=
T
0
exp (
−
γt
), where the decay rate (
γ
) is related to the
effective thermal conductivity (
k
eff
) and the volumetric heat
capacity of the solid (
C
)as
γ
=
k
eff
q
2
/C
. Therefore, to
obtain the effective thermal conductivity from our calculations,
we perform an inverse Fourier transform of the temperature
distribution averaged in the
z
direction [
!
d
0
̄
T
(
η,q,z
)
dz
]
with respect to the variable
η
, fit the resulting solution to an
exponentially decaying function
T
0
exp (
−
γt
), and extract
the thermal conductivity from the fit. If the fitting fails, the
transport is in the strongly quasiballistic regime [
25
] and
we conclude that the Fourier law description of the heat
conduction with an effective thermal conductivity
k
eff
is not
valid for that case.
The semianalytical solution of the BTE for transient
transport presented in this work is computationally very
efficient, taking only a few seconds on a single computer
processor, while the direct Monte Carlo simulation of the BTE
takes up to a few days on a high-performance computer cluster
executed in parallel mode. Moreover, it is computationally
challenging to extract the heat flux distribution directly from
the Monte Carlo solution, while in our semianalytical solution,
the evaluation of heat flux distribution is a single step process
[Eq. (
16
)].
III. RESULTS AND DISCUSSION
We now present the results of the calculations for free-
standing silicon thin films. To obtain these results, we use an
isotropic dispersion and intrinsic scattering rates calculated
using a Gaussian kernel-based regression [
32
]fromthe
ab
initio
phonon properties of isotopically pure silicon. The
first-principles phonon properties are calculated by Carrete and
Mingo using ShengBTE [
33
,
34
] and Phonopy [
35
] from the
interatomic force constants calculated using VASP [
36
–
39
].
035314-5
NAVANEETHA K. RAVICHANDRAN AND AUSTIN J. MINNICH
PHYSICAL REVIEW B
93
, 035314 (2016)
FIG. 2. (a) Comparison of the cross-plane distribution of the in-plane steady state heat flux between analytical and Monte Carlo solutions
of the BTE at 300 K and film thickness of 100 nm for different boundary conditions. The geometry of the thin film and the coordinate axes
used in this work are shown in the inset. (b) Comparison of the steady state thermal conductivity between analytical and Monte Carlo solutions
of the BTE at different temperatures and thin film thicknesses. For both (a) and (b), the Monte Carlo solutions are identical for thermalizing
and nonthermalizing boundary scattering and agree well with the analytical solution derived in this work for both fully diffuse and partially
specular boundary conditions (RMS 0.1 nm). For the partially specular boundary condition, the specularity parameter (
p
ω
) is calculated from
Ziman’s specularity model [
40
] for a surface rms roughness of 0.1 nm. (c) Effective MFPs of phonons computed using Matthiessen’s rule (MR)
and the Fuchs-Sondheimer (FS) theory for different film thicknesses and fully diffuse boundary scattering. Matthiessen’s rule underpredicts
the effective phonon MFPs in thin films compared to the Fuchs-Sondheimer theory, which is a rigorous BTE solution.
A. Steady state transport in thin films
1. Comparison with Monte Carlo solution
We first examine steady state heat condition along thin
films. Figures
2(a)
and
2(b)
show the cross-plane distribution
of the in-plane heat flux and the effective thermal conductivity,
respectively, for steady state transport through thin films
computed using a Monte Carlo technique and the analytical
solution from this work. The details of the Monte Carlo
technique used in this work are described in Sec. III of the
Supplemental Material [
31
]. For both fully diffuse and partially
specular boundary conditions, the heat flux distribution and the
effective thermal conductivity of the thin film show excellent
agreement between the Monte Carlo solutions and the ana-
lytical solution from this work over a range of temperatures
and film thicknesses. In particular, both solutions predict
identical heat flux and thermal conductivities for thermalizing
and nonthermalizing boundary conditions at the thin film
walls since the steady state transport is insensitive to the
type of diffuse boundary scattering of phonons, as discussed
in Sec.
II B
. This observation can be generalized further to
state that in steady state thermal transport experiments on thin
films, it is impossible to distinguish between nonthermalizing
and
any
type of inelastic diffuse scattering of phonons at
boundaries.
2. Effective phonon mean-free path
We also examine the effective mean-free path (MFP) of
phonons within the thin film for various film thicknesses. An
approach to estimate the effective phonon mean-free path in
thin films is by using Matthiessen’s rule [
40
] given by
1
ω,
eff
=
1
ω,
bulk
+
1
−
p
ω
1
+
p
ω
1
d
,
(17)
where
d
is the thickness of the thin film and
ω,
bulk
is the
intrinsic phonon mean-free path in the bulk material. Although
Matthiessen’s rule has been used in the past for computational
[
41
] and experimental [
42
] investigations of phonon boundary
scattering, the mathematical rigor of such an expression for
the effective mean-free path is unclear. On the other hand,
the effective mean-free path of phonons in thin films can also
be determined rigorously from the Fuchs-Sondheimer factor
[
F
(
ω
/d
)], since by definition,
F
(
ω
/d
)
=
k
ω,
eff
/k
ω,
bulk
=
ω,
eff
/
ω,
bulk
. Figure
2(c)
shows the comparison of the
normalized effective phonon mean-free paths obtained from
the Fuchs-Sondheimer factor and Matthiessen’s rule for
different film thicknesses. Matthiessen’s rule underpredicts
phonon MFPs comparable to the thickness of the film. Even for
phonons with intrinsic mean-free path an order of magnitude
smaller than the film thickness, Matthiessen’s rule predicts
a shorter effective phonon mean-free path compared to the
predictions of the Fuchs-Sondheimer factor from the rigorous
solution of the BTE, which is consistent with the findings of
another work based on Monte Carlo sampling [
43
]. This result
highlights the importance of using the rigorous BTE solution
to estimate the extent of diffuse phonon boundary scattering
even in simple nanostructures.
B. Transient transport in thin films
We now examine transient thermal conduction along thin
films observed in the TG experiment. To perform this calcula-
tion, we solve the integral equation [Eq. (
15
)] semianalytically
using the same isotropic phonon properties used in steady
state transport calculations. The source term in the BTE
[Eq. (
10
)] is assumed to follow a thermal distribution given
by
Q
ω
=
C
ω
T
0
, where
C
ω
is the volumetric specific heat of
the phonon mode.
035314-6
ROLE OF THERMALIZING AND NONTHERMALIZING . . .
PHYSICAL REVIEW B
93
, 035314 (2016)
FIG. 3. (a) Comparison between time traces from the Monte Carlo (colored noisy lines) and the degenerate kernels solutions (black lines)
of the BTE for a grating period of 20
μ
m at 500 K. (b) Comparison of the thermal conductivity predictions from the Monte Carlo (symbols)
and the degenerate kernels solutions (black lines) of the BTE for different temperatures and grating periods. For both (a) and (b), the Monte
Carlo solutions and the BTE solutions from this work are in very good agreement. (c) Plot showing the difference between the incoming and
outgoing heat flux normalized by the incoming heat flux at the thin film boundaries as a function of the temporal frequency normalized by
the maximum temporal frequency at which the simulations were performed (
η
max
). Specular and nonthermalizing diffuse boundary conditions
conserve heat flux to numerical precision while the thermalizing diffuse boundary condition violates heat flux conservation at the film wall
under quasiballistic (
T
=
100 K, grating period
=
1
μ
m) and diffusive (
T
=
500 K, grating period
=
1000
μ
m) transport regimes.
1. Difference between thermalizing and nonthermalizing
boundary scattering
Figure
3(a)
shows a comparison of the time traces calculated
from the degenerate kernel method and the Monte Carlo
method for a grating period of 20
μ
m. To obtain the thickness
averaged time traces using the degenerate kernels method, we
used 6 Gauss quadrature points for the
μ
variable, 10 Gauss
quadrature points for the
φ
variable, and 1 term in the Fourier
cosine expansion to achieve a convergence threshold of 10
−
6
on the relative change in the solution for an increase in the
number of discretization points. The transient decays are in
good agreement between the degenerate kernel and the Monte
Carlo solutions over a wide range of temperatures and different
boundary conditions. As expected, the solution for the specular
boundary condition results in a faster transient decay than
the diffuse boundary conditions since a specularly reflecting
wall does not resist the flow of heat in the in-plane direction.
However, the transient decay for the nonthermalizing diffuse
boundary condition is faster that the thermalizing diffuse
boundary condition, indicating that the thermalizing boundary
condition offers higher resistance to heat flow than the
nonthermalizing diffuse scattering.
This observation is also evident from Fig.
3(b)
which
shows the thermal conductivities obtained by fitting the time
traces to an exponential decay for different temperatures,
different grating periods, and different boundary conditions.
The observed thermal conductivity of the thin film decreases
with decreasing grating period due to the breakdown of
Fourier’s law of heat conduction and the onset of quasiballistic
thermal transport [
25
] when the grating period is comparable
to phonon MFPs. Consistent with the findings from the time
traces, the thermal conductivity of the thin film with specular
walls is higher than that of the thin film with diffuse walls.
Moreover, even for very long grating periods compared to
phonon MFPs, where the thermal transport is diffusive and
obeys Fourier’s law, the thermal conductivity of thin film with
nonthermalizing diffuse walls is higher than that of the thin
films with thermalizing diffuse walls. This observation is in
stark contrast with the steady state condition, where there was
no difference in thermal conductivity between thermalizing
and nonthermalizing boundary conditions.
2. Validity of the thermalizing and nonthermalizing
boundary conditions
At this point, it is important to investigate the validity of
the thermalizing and nonthermalizing boundary condition for
the thin film walls. The nonthermalizing boundary scattering
condition can be naturally derived from the conservation of
heat flux at the boundary [
27
]. However, the thermalizing
boundary condition is not derived from the heat flux con-
servation at the boundary. Therefore, in the absence of any
external scattering mechanisms, phonons cannot reach the
local thermal equilibrium and simultaneously conserve heat
flux at the boundary in general, because of the following
reason.
Consider a boundary at
z
=
0 separating a solid at
z>
0
from vacuum in
z<
0. The incoming phonon distribution at
z
=
0is
g
−
ω
(0
,μ,φ
), which is a general phonon distribution,
not necessarily at the local thermal equilibrium. According
to the formulation of the thermalizing diffuse boundary
condition, the outgoing phonon distribution, in the case of
fully diffuse boundary scattering, is given by
g
0
(
T
(
z
=
0)
)
≈
C
ω
4
π
T
(
z
=
0), where
C
ω
is the heat capacity of the phonon
mode and
T
(
z
=
0) is the local equilibrium temperature at
the boundary
z
=
0. Since the boundary separates a solid from
vacuum, all of the heat flux incident on the boundary has to
be reflected back into the solid. This constraint on the incident
and reflected heat flux at the thermalizing diffuse boundary
035314-7