of 49
arXiv:1511.03312v1 [cond-mat.mtrl-sci] 10 Nov 2015
The Role of Thermalizing and Non-thermalizing Walls in Phon
on
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
(Dated: November 16, 2015)
Abstract
Phonon boundary scattering is typically treated using the F
uchs-Sondheimer theory, which as-
sumes that phonons are thermalized to the local temperature
at the boundary. However, whether
such a thermalization process actually occurs and its effect o
n thermal transport remains unclear.
Here we examine thermal transport along thin films with both t
hermalizing and non-thermalizing
walls by solving the spectral Boltzmann transport equation
(BTE) for steady state and transient
transport. We find that in steady state, the thermal transpor
t is governed by the Fuchs-Sondheimer
theory and is insensitive to whether the boundaries are ther
malizing or not. In contrast, under
transient conditions, the thermal decay rates are significa
ntly different for thermalizing and non-
thermalizing walls. We also show that, for transient transp
ort, the thermalizing boundary condition
is unphysical due to violation of heat flux conservation at th
e boundaries. Our results provide in-
sights into the boundary scattering process of thermal phon
ons over a range of heating length scales
that are useful for interpreting thermal measurements on na
nostructures.
aminnich@caltech.edu
1
I. INTRODUCTION
Engineering the thermal conductivity of nanoscale materials has be
en a topic of con-
siderable research interest over the past two decades [1]. While ap
plications such as GaN
transistors [2, 3] and light emitting diodes (LEDs) [4] require high th
ermal conductivity sub-
strates to dissipate heat, the performance of thermoelectric an
d thermal insulation devices
can be significantly enhanced by reducing their thermal conductivit
y [5, 6]. In many of these
applications, phonon boundary scattering is the dominant resistan
ce to heat flow, making
the detailed understanding of this process essential for advancin
g applications.
Phonon boundary scattering has been studied extensively both th
eoretically and exper-
imentally. The thermal conductivity reduction due to boundary sca
ttering of phonons is
conventionally treated using the Fuchs-Sondheimer theory, which
was first derived for elec-
tron boundary scattering independently by Fuchs [7] and Reuter a
nd 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 importa
nt assumption that the
diffusely scattered part of the phonon spectrum at a partially spec
ular 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 sc
attering limit of
Casimir
'
s theory [12].
Several computational works [11, 13–16] have studied the reduc
tion in thermal con-
ductivity due to phonon boundary scattering in nanostructures b
y solving the phonon
Boltzmann transport equation (BTE). These works have consider
ed either thermalizing
or non-thermalizing boundaries but have never compared the effec
t of these two different
boundary conditions on the thermal conductivity of nanostructu
res. Several experimental
works have also studied the reduction in thermal conductivity of na
nomaterials such as
nanowires [17–19], thin films [10, 20, 21] and nanopatterned struct
ures [22] due to phonon
boundary scattering. These works have used the Fuchs-Sondhe
imer theory to interpret their
measurements. However, it is not clear if the assumptions made in th
e Fuchs-Sondheimer
2
theory are necessarily applicable for these experiments. In fact,
an analysis of the effect
of the key assumption made in the Fuchs-Sondheimer theory, that
the walls are thermal-
izing, has never been investigated due to the challenges involved in so
lving the BTE for
non-thermalizing walls.
Here, we examine the role of thermalizing and non-thermalizing walls in h
eat conduction
along thin films by solving the spectral phonon Boltzmann transport
equation (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 the
rmalized or not at the bound-
aries and that Fuchs-Sondheimer theory accurately describes th
ermal transport along the
thin film. In the case of transient transport, we find that the deca
y rates are significantly
different for thermalizing and non-thermalizing walls and that Fuchs-
Sondheimer theory
accurately predicts the thermal conductivity only when the therm
al transport is diffusive.
Moreover, under transient transport conditions, we find that ph
onons cannot undergo ther-
malization at the boundaries in general due to the violation of heat flu
x conservation. Our
results provide insights into the boundary scattering process of t
hermal phonons that are
useful for interpreting thermal measurements on nanostructu
res.
II. MODELING
A. Boltzmann Transport Equation
We begin our analysis by considering the two dimensional spectral tr
ansient Boltzmann
transport equation (BTE) under the relaxation time approximation
for an isotropic crystal,
given by,
∂g
ω
∂t
+
μv
g
∂g
ω
∂z
+
v
g
p
1
μ
2
cos
φ
∂g
ω
∂x
=
g
ω
g
o
(
T
)
τ
ω
+
Q
ω
4
π
(1)
Here,
g
ω
is the phonon energy distribution function,
ω
is the phonon frequency,
v
g
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
3
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 bou
ndary
conditions enforce that the diffusely scattered phonons are ther
malized while also allowing
some phonons to be specularly reflected. Here, we generalize thes
e boundary conditions to
allow for the possibility of both partial thermalization and partial spe
cularity as:
For
μ
(0
,
1]
,
g
+
ω
(0
, μ, φ
) =
p
ω
g
ω
(0
,
μ, φ
)
+ (1
p
ω
)
σ
ω
C
ω
T
(
z
= 0)
4
π
(1
σ
ω
)
π
Z
2
π
0
Z
0
1
g
ω
(0
, μ
, φ
)
μ
d
μ
d
φ
!
For
μ
[
1
,
0)
,
g
ω
(
d, μ, φ
) =
p
ω
g
+
ω
(
d,
μ, φ
)
+ (1
p
ω
)
σ
ω
C
ω
T
(
z
=
d
)
4
π
+
(1
σ
ω
)
π
Z
2
π
0
Z
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, μ, φ
) is the 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
ω
,
p
ω
and
σ
ω
are the phonon
specularity parameter and the thermalization parameter for the t
hin film walls respectively.
The specularity parameter represents the fraction of specularly
scattered phonons at the
boundaries and the thermalization parameter represents the fra
ction of the phonon distri-
bution that is absorbed and reemitted at the local equilibrium temper
ature of the thin film
walls. For simplicity, we ignore mode conversion for non-thermalizing b
oundary condition
in our analysis.
The unknown quantities in this problem are the phonon distribution fu
nction (
g
ω
(
t, x, z, μ, φ
))
4
and the deviational temperature distribution (∆
T
(
t, x, z
)). They are related to each other
through the energy conservation requirement,
Z
ω
m
ω
=0
Z
1
μ
=
1
Z
2
π
φ
=0

g
ω
τ
ω
1
4
π
C
ω
τ
ω
T

d
φ
d
μ
d
ω
= 0
(3)
Due to the high dimensionality of the BTE, analytical or semi-analytica
l solutions are only
available in literature for either semi-infinite domains [23–25] or domain
s with simple bound-
ary and transport conditions [26] or with several approximations [
27]. For nanostructures
with physically realistic boundaries, several numerical solutions of t
he BTE have been re-
ported [11, 16, 28]. However, computationally efficient analytical or
semi-analytical solutions
for the in-plane heat conduction along even simple unpatterned films
[10, 20] are unavail-
able. To overcome this problem, we solve the BTE analytically for stea
dy state transport
(section II B) and semi-analytically for transient transport along t
hin films in the TG exper-
iment [10, 20] (section II C).
B. Steady State Heat Conduction in Thin Films
In this section, we extend the Fuchs-Sondheimer relation for ther
mal conductivity sup-
pression due to phonon boundary scattering to the general boun
dary conditions described in
equation 2. To simulate steady state transport,
Q
ω
is set to 0 in the BTE (equation 1). Fur-
thermore, we assume that a one-dimensional temperature gradie
nt exists along the thin film
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
μ
∂g
ω
∂z
+
v
g
p
1
μ
2
cos
φ
∂g
0
ω
∂x
=
g
ω
g
0
ω
τ
ω
(4)
For steady state transport, it is convenient to solve the BTE in ter
ms of the deviation from
equilibrium distribution ( ̄
g
ω
=
g
ω
g
0
ω
(∆
T
(
x
))). In this case, the BTE transforms into,
̄
g
ω
∂z
+
̄
g
ω
μ
Λ
ω
=
cos
φ
p
1
μ
2
μ
∂g
0
ω
∂x
(5)
5
The boundary conditions (equation 2) for ̄
g
ω
now become,
For
μ
(0
,
1]
,
̄
g
+
ω
(0
, μ, φ
) =
p
ω
̄
g
ω
(0
,
μ, φ
)
(1
p
ω
) (1
σ
ω
)
π
Z
2
π
0
Z
0
1
̄
g
ω
(0
, μ
, φ
)
μ
d
μ
d
φ
For
μ
[
1
,
0)
,
̄
g
ω
(
d, μ, φ
) =
p
ω
̄
g
+
ω
(
d,
μ, φ
) +
(1
p
ω
) (1
σ
ω
)
π
Z
2
π
0
Z
1
0
̄
g
+
ω
(
d, μ
, φ
)
μ
d
μ
d
φ
(6)
The general solution of the BTE (equation 5) along with the boundar
y conditions (equa-
tion 6) is given by,
̄
g
+
ω
(
z, μ, φ
) =
Λ
ω
cos
φ
p
1
μ
2
∂g
0
ω
∂x
1
exp

z
μ
Λ
ω

(1
p
ω
)
1
p
ω
exp

d
μ
Λ
ω

+
(1
p
ω
) (1
σ
ω
)
h
A
+
ω
+
p
ω
exp

d
μ
Λ
ω

A
ω
i
1
p
2
ω
exp

2
d
μ
Λ
ω

exp

z
μ
Λ
ω

|
{z
}
I
̄
g
ω
(
z,
μ, φ
) =
Λ
ω
cos
φ
p
1
μ
2
∂g
0
ω
∂x
1
exp

(
d
z
)
μ
Λ
ω

(1
p
ω
)
1
p
ω
exp

d
μ
Λ
ω

+
(1
p
ω
) (1
σ
ω
)
h
A
ω
+
p
ω
exp

d
μ
Λ
ω

A
+
ω
i
1
p
2
ω
exp

2
d
μ
Λ
ω

exp

(
d
z
)
μ
Λ
ω

|
{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 expres-
sions for ̄
g
+
ω
(
z, μ, φ
) and ̄
g
ω
(
z,
μ, φ
) (equation 7) is shown in section I of the supplementary
material. The expression for the in-plane (
x
direction) spectral heat flux is given by,
q
x,ω
=
1
d
Z
d
z
=0
Z
1
μ
=
1
Z
2
π
φ
=0
v
x
̄
g
ω
D
(
ω
)
4
π
d
φ
d
μ
d
z
=
"
1
3
C
ω
v
g
Λ
ω
#
∂T
∂x
1
3 (1
p
ω
) Λ
ω
2
d
Z
1
0
μ
μ
3

1
exp

d
μ
Λ
ω

1
p
ω
exp

d
μ
Λ
ω

d
μ
(8)
6
since the diffuse contributions to the distribution functions ̄
g
+
ω
(
z, μ, φ
) and ̄
g
ω
(
z,
μ, φ
)
(terms I and II in equation 7) are independent of the azimuthal ang
le
φ
and integrate out
to 0. Comparing equation 8 with the expression for heat flux from th
e Fourier’s law, the
spectral effective thermal conductivity of the thin film is obtained a
s a product of the bulk
spectral thermal conductivity and the well-known Fuchs-Sondhe
imer reduction factor due to
phonon boundary scattering given by,
k
ω,
eff
(
d
) =
"
1
3
C
ω
v
g
Λ
ω
#
|
{z
}
k
ω,
bulk
1
3 (1
p
ω
) Λ
ω
2
d
Z
1
0
μ
μ
3

1
exp

d
μ
Λ
ω

1
p
ω
exp

d
μ
Λ
ω

d
μ
|
{z
}
Fuchs
Sondheimer reduction factor
F
(
Λ
ω
d
)
(9)
It is interesting to observe from equation 9 that the spectral effe
ctive thermal conductivity is
independent of the thermalization parameter
σ
ω
even though a general boundary condition
(equation 2) has been used in this derivation. Thus, the steady sta
te thermal conductivity
suppression due to boundary scattering is only influenced by the re
lative extent of specular
and diffuse scattering (parameterized by the specularity paramet
er
p
ω
) and does not depend
on the type of diffuse scattering process (parameterized by the t
hermalization parameter
σ
ω
). We explicitly demonstrate this result using numerical simulations in s
ection III A.
C. Transient Heat Conduction in Thin Films
In this section, we solve the BTE (equation 1) for transient therma
l transport along a thin
film. The initial temperature profile considered in this work is identical
to that which occurs
in the Transient Grating (TG) experiment, which has been used exte
nsively 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 d
ecay of a one-dimensional
impulsive sinusoidal temperature grating on the sample at different g
rating periods. In the
large grating period limit of heat diffusion, the temporal decay is a sing
le exponential. Since
the initial temperature distribution is an infinite one-dimensional sinu
soid in the
x
direction,
the temperature distribution remains spatially sinusoidal at all later
times. Therefore, each
7
wave vector
q
in the spatially Fourier transformed BTE directly corresponds to a u
nique
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 equation 1 in the time variable
t
. With these transformations, the BTE reduces to,
iηG
ω
+
μv
g
∂G
ω
∂z
+
iqv
g
p
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 en
ergy distribution function
g
ω
.
The outline of the solution methodology for equation 10 is as follows. T
he general solution
is given by,
For
μ
(0
,
1]
, G
+
ω
(
z, μ, φ
) =
G
+
ω
(0
, μ, φ
) exp
γ
FS
μφ
μ
Λ
ω
z
!
+
exp

γ
FS
μφ
μ
Λ
ω
z

4
πμ
Λ
ω
Z
z
0
C
ω
̄
T
+
̄
Q
ω
τ
ω

exp
γ
FS
μφ
μ
Λ
ω
z
!
d
z
For
μ
[
1
,
0)
, G
ω
(
z, μ, φ
) =
G
ω
(
d, μ, φ
) exp
γ
FS
μφ
μ
Λ
ω
(
d
z
)
!
exp

γ
FS
μφ
μ
Λ
ω
z

4
πμ
Λ
ω
Z
d
z
C
ω
̄
T
+
̄
Q
ω
τ
ω

exp
γ
FS
μφ
μ
Λ
ω
z
!
d
z
where
, γ
FS
μφ
= (1 +
iητ
ω
) +
i
Λ
ω
q
p
1
μ
2
cos
φ
(11)
Here,
G
+
ω
(0
, μ, φ
) and
G
ω
(
d, μ, φ
) are determined by solving the boundary conditions (equa-
tion 2) with the following procedure. First, the angular integrals in th
e boundary conditions
are discretized using Gauss quadrature, which results in the followin
g 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
8
from the discretization.
G
+
ω
(0
, μ
i
, φ
j
) =
p
ω
G
ω
(
d,
μ
i
, φ
j
) exp
γ
FS
ij
μ
i
Λ
ω
d
!
+
p
ω
4
πμ
i
Λ
ω
Z
d
0

C
ω
̄
T
+
̄
̃
Q
ω
τ
ω

exp
γ
FS
ij
μ
i
Λ
ω
z
!
d
z
+ (1
p
ω
)
"
σ
ω
C
ω
̄
T
(
z
= 0)
4
π
+
(1
σ
ω
)
π
X
i
j
G
ω
d,
μ
i
, φ
j

exp
γ
FS
i
j
μ
i
Λ
ω
d
!
μ
i
w
μ
i
w
φ
j
+
(1
σ
ω
)
4
π
2
Λ
ω
X
i
j
Z
d
0

C
ω
̄
T
+
̄
̃
Q
ω
τ
ω

exp
γ
FS
i
j
μ
i
Λ
ω
z
!
d
z
w
μ
i
w
φ
j
#
G
ω
(
d,
μ
i
, φ
j
) =
p
ω
G
+
ω
(0
, μ
i
, φ
j
) exp
γ
FS
ij
μ
i
Λ
ω
d
!
+
p
ω
4
πμ
i
Λ
ω
Z
d
0

C
ω
̄
T
+
̄
̃
Q
ω
τ
ω

exp
γ
FS
ij
μ
i
Λ
ω
(
d
z
)
!
d
z
+ (1
p
ω
)
"
σ
ω
C
ω
̄
T
(
z
=
d
)
4
π
+
(1
σ
ω
)
π
X
i
j
G
+
ω
0
, μ
i
, φ
j

exp
γ
FS
i
j
μ
i
Λ
ω
d
!
μ
i
w
μ
i
w
φ
j
+
(1
σ
ω
)
4
π
2
Λ
ω
X
i
j
Z
d
0

C
ω
̄
T
+
̄
̃
Q
ω
τ
ω

exp
γ
FS
i
j
μ
i
Λ
ω
(
d
z
)
!
d
z
w
μ
i
w
φ
j
#
(12)
To obtain equation 12, we have substituted the general BTE solutio
n into the boundary
conditions to eliminate
G
ω
(0
, μ, φ
) and
G
+
ω
(
d, μ, φ
). Therefore, the only unknowns in the
set of linear equations (equation 12) are
G
+
ω
(0
, μ, φ
) and
G
ω
(
d, μ, φ
). By bringing the terms
containing
G
+
ω
(0
, μ, φ
) and
G
ω
(
d, μ, φ
) to the left hand side, equation 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)
where,
̄
̃
c
+
ω
(0
, μ
i
, φ
j
) and
̄
̃
c
ω
(
d, μ
i
, φ
j
) are analytical functions of the unknown temperature
9
distribution function ∆
̄
T
obtained from the right hand side of equation 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 section II A of the supplementary material. To clos
e the
problem, the expressions for
G
+
ω
(
z, μ, φ
) and
G
ω
(
z, μ, φ
) (equation 11) and the boundary
conditions (equation 14) are substituted into the energy conserv
ation equation (equation 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
) +
Z
d
0

K
(
z
, z
) ∆
̄
T
(
z
)

d
z
(15)
where the functional form of the inhomogeneous parts
f
(
z
),
h
(
z
) and the kernel
K
(
z
, z
) are
described in section II B of the supplementary material. This integra
l equation (equation 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 section II C of the supplementary material. Finally, t
he solution ∆
̄
T
(
z
)
is substituted into equation 11 to obtain expressions for
G
ω
(
z, μ, φ
) and also the thickness-
10
averaged in-plane heat flux
j
x,ω
given by,
j
x,ω
=
1
4
πd
Z
d
0
Z
2
π
0
Z
1
1
G
ω
v
g
p
1
μ
2
cos
φ
d
μ
d
φ
d
z
=
1
4
π
X
ij
μ
i
Kn
d
ω
γ
FS
ij
X
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
X
m
=1
t
m
1 + (
1)
m
m
2
π
2
+

γ
FS
ij
μ
i
Kn
d
ω

2
v
g
q
1
μ
2
i
w
μ
i
w
φ
j
(16)
where
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 th
e thermal
transport properties of the thin film is to compare the expression f
or heat flux from the BTE
solution with that expected for heat diffusion, as was done in equatio
n 8 for the steady state
Fuchs-Sondheimer theory. However, in practice, equation 16 is no
t easily reduced into the
form of Fourier’s law. To overcome this problem, the following strate
gy is adopted. The solu-
tion of the Fourier heat equation to a one-dimensional heat conduc
tion 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 vol-
umetric 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 (
R
d
0
̄
T
(
η, q, z
) d
z
) 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 f
ails, the transport is in the
strongly quasi-ballistic regime [25] and we conclude that the Fourier la
w description of the
heat conduction with an effective thermal conductivity
k
eff
is not valid for that case.
11
The semi-analytical solution of the BTE for transient transport pr
esented in this work
is computationally very efficient, taking only a few seconds on a single c
omputer processor,
while the direct Monte Carlo simulation of the BTE takes up to a few day
s on a high-
performance computer cluster executed in parallel mode. Moreov
er, it is computationally
challenging to extract the heat flux distribution directly from the Mo
nte Carlo solution,
while in our semi-analytical solution, the evaluation of heat flux distrib
ution is a single step
process (equation 16).
III. RESULTS & DISCUSSION
We now present the results of the calculations for free-standing s
ilicon thin films. To
obtain these results, we use an isotropic dispersion and intrinsic sca
ttering rates calculated
using a Gaussian kernel-based regression [31] from the
ab-initio
phonon properties of iso-
topically pure silicon. The first principles phonon properties are calcu
lated by
J. Carette &
N. Mingo
using ShengBTE [32, 33] and Phonopy [34] from the inter-atomic for
ce constants
calculated using VASP [35–38].
A. Steady State Transport in Thin Films
1. Comparison with Monte Carlo Solution
We first examine steady state heat condition along thin films. Figures
1 (a) and (b) show
the cross-plane distribution of the in-plane heat flux and the effect
ive thermal conductivity
respectively, for steady state transport through thin films comp
uted using a Monte Carlo
technique and the analytical solution from this work. The details of t
he Monte Carlo tech-
nique used in this work is described in section III of the the supplemen
tary material. For both
fully diffuse and partially specular boundary conditions, the heat flux
distribution and the
effective thermal conductivity of the thin film show excellent agreem
ent between the Monte
Carlo solutions and the analytical solution from this work over a rang
e of temperatures and
12
film thicknesses. In particular, both solutions predict identical hea
t flux and thermal con-
ductivities for thermalizing and non-thermalizing boundary condition
s at the thin film walls
since the steady state transport is insensitive to the type of diffus
e boundary scattering of
phonons, as discussed in section II B. This observation can be gene
ralized further to state
that in steady state thermal transport experiments on thin films,
it is impossible to distin-
guish between non-thermalizing 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 pho
non mean free path in
thin films is by using the Matthiessen’s rule [39] 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 the Matthiessen’s rule has been us
ed in the past
for computational [40] and experimental [41] investigations of ph
onon boundary scatter-
ing, the mathematical rigor of such an expression for effective mea
n free path is unclear.
On the other hand, the effective mean free path of phonons in thin fi
lms can also be
determined rigorously from the Fuchs-Sondheimer factor (
F
ω
/d
)), since by definition,
F
ω
/d
) =
k
ω,
eff
/k
ω,
bulk
= Λ
ω,
eff
/
Λ
ω,
bulk
. Figure 1 (c) shows the comparison of the nor-
malized 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, Matthies
sen’s rule predicts a
shorter effective phonon mean free path compared to the predict
ions of the Fuchs-Sondheimer
factor from the rigorous solution of the BTE, which is consistent wit
h the findings of another
work based on Monte Carlo sampling [42]. This result highlights the impor
tance of using the
rigorous BTE solution to estimate the extent of diffuse phonon boun
dary scattering even in
13
(a)
(b)
100 K
300 K
500 K
x
z
Thin Film Walls
Fully Diffuse
Partially
Specular
(c)
1 mm
1 nm
100 nm
10 μm
MR
FS
FIG. 1. (a) Comparison of the cross-plane distribution of th
e in-plane steady state heat flux
between analytical and Monte Carlo solutions of the BTE at 30
0 K and film thickness of 100 nm
for different boundary conditions. The geometry of the thin fil
m and the coordinate axes used
in this work are shown in the inset. (b) Comparison of the stea
dy state thermal conductivity
between analytical and Monte Carlo solutions of the BTE at di
fferent temperatures and thin film
thicknesses. For both (a) and (b), the Monte Carlo solutions
are identical for thermalizing and non-
thermalizing boundary scattering and agree well with the an
alytical solution derived in this work
for both fully diffuse and partially specular boundary condit
ions (RMS 0.1 nm). For the partially
specular boundary condition, the specularity parameter (
p
ω
) is calculated from Ziman’s specularity
model [39] for a surface RMS roughness of 0.1 nm. (c) Effective M
FPs of phonons computed using
the Matthiessen’s rule (MR) and the Fuchs-Sondheimer (FS) t
heory for different film thicknesses
and fully diffuse boundary scattering. Matthiessen’s rule un
derpredicts the effective phonon MFPs
in thin films compared to the Fuchs-Sondheimer theory, which
is a rigorous BTE solution.
simple nanostructures.
14
B. Transient Transport in Thin Films
We now examine transient thermal conduction along thin films observ
ed in the TG ex-
periment. To perform this calculation, we solve the integral equatio
n (equation 15) semi-
analytically using the same isotropic phonon properties used in stead
y state transport calcu-
lations. The source term in the BTE (equation 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.
1. Difference between Thermalizing and Non-thermalizing Bou
ndary Scattering
Figure 2 (a) shows a comparison of the time traces calculated from t
he degenerate kernel
method and the Monte Carlo method for a grating period of 20
μ
m. The transient decays
are in good agreement between the degenerate kernel and the Mo
nte 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 deca
y than the diffuse boundary
conditions since a specularly reflecting wall does not resist the flow o
f heat in the in-plane
direction. However, the transient decay for the non-thermalizing
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 n
on-thermalizing diffuse
scattering.
This observation is also evident from figure 2 (b) which shows the the
rmal conductivi-
ties obtained by fitting the time traces to an exponential decay for
different temperatures,
different grating periods and different boundary conditions. The ob
served thermal conduc-
tivity of the thin film decreases with decreasing grating period due to
the breakdown of the
Fourier’s law of heat conduction and the onset of quasiballistic therm
al transport [25] when
the grating period is comparable to phonon MFPs. Consistent with th
e findings from the
time traces, the thermal conductivity of the thin film with specular w
alls 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 F
ourier’s law, the thermal
15
conductivity of thin film with non-thermalizing diffuse walls is higher than
that of the thin
films with thermalizing diffuse walls. This observation is in stark contras
t with the steady
state condition, where there was no difference in thermal conduct
ivity between thermalizing
and non-thermalizing boundary conditions.
T = 500 K
Grating Period : 20 μm
500 K
100 K
Lines : Degenerate Kernels
Symbols : Monte Carlo
(b)
(a)
(c)
T = 100 K
Grating Period = 1 μm
T = 500 K
Grating Period = 1000 μm
FIG. 2. (a) Comparison between time traces from the Monte Car
lo (colored noisy lines) and the
degenerate kernels solutions (black lines) of the BTE for a g
rating period of 20
μm
at 500 K.
(b) Comparison of the thermal conductivity predictions fro
m the Monte Carlo (symbols) and the
degenerate kernels solutions (black lines) of the BTE for di
fferent temperatures and grating periods.
For both (a) and (b), the Monte Carlo solutions and the BTE sol
utions from this work are in very
good agreement. (c) Plot showing the heat flux conservation a
t the film boundaries. Specular
and non-thermalizing diffuse boundary conditions conserve h
eat flux to numerical precision while
thermalizing diffuse boundary condition violates heat flux co
nservation at the film wall under
quasiballistic (T = 100 K, grating period = 1
μ
m) and diffusive (T = 500 K, grating period =
1000
μ
m) transport regimes.
2. Validity of the Thermalizing and Non-thermalizing Bounda
ry Conditions
At this point, it is important to investigate the validity of the thermaliz
ing and non-
thermalizing boundary condition for the thin film walls. The non-therm
alizing boundary
16
scattering condition can be naturally derived from the conservatio
n of heat flux at the
boundary [27]. However, the thermalizing boundary condition is not d
erived from the heat
flux conservation at the boundary. Therefore, in the absence of
any external scattering
mechanisms, phonons cannot reach the local thermal equilibrium an
d simultaneously con-
serve heat flux at the boundary in general, due to the following reas
on.
Consider a boundary at
z
= 0 separating a solid at
z >
0 from vacuum in
z <
0. The
incoming phonon distribution at
z
= 0 is
g
ω
(0
, μ, φ
), which is a general phonon distribu-
tion, not necessarily at the local thermal equilibrium. According to t
he formulation of the
thermalizing diffuse boundary condition, the outgoing phonon distrib
ution, 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 he
at
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
leads to the following
relation for ∆
T
(
z
= 0).
T
(
z
= 0) = 4
P
p
R
ω
max
ω
=0
R
0
μ
=
1
R
2
π
φ
=0
g
ω
(0
, μ, φ
)
v
g
μ
d
μ
d
φ
d
ω
P
p
R
ω
max
ω
=0
C
ω
v
g
d
ω
(18)
Additionally, energy conservation (equation 3) has to be satisfied a
t all locations including
the boundaries in the absence of any other source or sink of phono
ns. This requirement
further adds constraints on ∆
T
(
z
= 0) through the relation,
T
(
z
= 0) = 2
P
p
R
ω
max
ω
=0
R
0
μ
=
1
R
2
π
φ
=0
g
ω
(0
,μ,φ
)
τ
ω
d
μ
d
φ
d
ω
P
p
R
ω
max
ω
=0
C
ω
τ
ω
d
ω
(19)
For the assumptions made in the Fuchs-Sondheimer theory under s
teady state transport con-
ditions, the integrals of the incoming and the outgoing distribution fu
nctions (equation 7)
over the azimuthal angle
φ
are 0. Therefore, there is no heat flux towards or away from the
boundary and the constraints on ∆
T
(
z
= 0) (given by equations 18 and 19) are trivially
satisfied. However, in general, these two expressions for ∆
T
(
z
= 0) are not equal, indicating
phonons cannot thermalize at the boundaries in the absence of any
external source or sink
17
of phonons.
Figure 2 (c) shows the difference between the incoming and outgoing
total heat flux at
the thin film wall (
z
= 0) as a function of the temporal frequency
η
. The specular and
non-thermalizing diffuse boundary conditions satisfy heat flux cons
ervation to numerical
precision. However, there is a significant difference between the inc
oming and the outgoing
heat flux for the thermalizing diffuse boundary condition under quas
iballistic (T = 100
K, grating period = 1
μ
m) and diffusive (T = 500 K, grating period = 1000
μ
m) transport
regimes. Nevertheless, it is still possible for inelastic (but not therm
alizing) diffuse boundary
scattering to take place as long as the following conditions for heat fl
ux are met at the thin
film boundaries:
X
p
Z
ω
max
ω
=0
g
+
ω
(
z
= 0)
v
g
d
ω
=
1
π
X
p
Z
ω
max
ω
=0
Z
0
μ
=
1
Z
2
π
φ
=0
g
ω
(
z
= 0
, μ, φ
)
v
g
μ
d
μ
d
φ
d
ω
X
p
Z
ω
max
ω
=0
g
ω
(
z
=
d
)
v
g
d
ω
=
1
π
X
p
Z
ω
max
ω
=0
Z
1
μ
=0
Z
2
π
φ
=0
g
+
ω
(
z
=
d, μ, φ
)
v
g
μ
d
μ
d
φ
d
ω
3. Comparison with Fuchs-Sondheimer Theory at Different Grati
ng Periods
We now examine if the Fuchs-Sondheimer theory can be used to expla
in transient heat
conduction in the TG experiment along thin films. If the suppression in
thermal conduc-
tivity of thin films due to phonon boundary scattering and quasiballist
ic effects in the TG
experiment are assumed to be independent, Fuchs-Sondheimer th
eory can be employed to
describe quasiballistic transport in the TG experiment using the follow
ing expression:
k
(
q, d
) =
X
p
Z
ω
max
0
F

p
ω
,
Λ
ω
d

S
(
q
Λ
ω
)

1
3
C
ω
v
g
Λ
ω

d
ω
(20)
where
F
p
ω
,
Λ
ω
d

is the Fuchs-Sondheimer suppression function from the steady st
ate trans-
port condition and
S
(
q
Λ
ω
) is the quasiballistic suppression function [25] for a grating period
q
. Recent works [20] have used a similar expression for the thermal
conductivity suppression
of the form:
k
(
q, d
) =
X
p
Z
ω
max
0
F

p
ω
,
Λ
ω
d

S

q
Λ
ω
F

p
ω
,
Λ
ω
d

1
3
C
ω
v
g
Λ
ω

d
ω
(21)
18
Henceforth, equation 20 is referred to as FS I and equation 21 is re
ferred to as FS II.
Figure 3 (a) shows the comparison of thermal conductivity obtaine
d by fitting the BTE
solution for temperature decay, and thermal conductivities from
FS I and FS II models for
fully diffuse boundary scattering. We only consider non-thermalizing
diffuse scattering as
we have shown that thermalizing diffuse scattering is unphysical for
the problem considered
here. At very long grating periods, when the transport is primarily d
iffusive, the thermal
conductivity predictions from FS I and FS II match well with the BTE s
olution from this
work, as expected. However, at the shorter grating periods com
parable to phonon MFPs,
where the transport is in the quasiballistic regime, FS I underpredict
s the thin film thermal
conductivity while FS II overpredicts it.
This observation is also evident from the magnitude of the suppress
ion function plotted at
η
= 0 for fully diffuse boundary conditions shown in figures 3 (b) and (c)
. The suppression
function for the thin film geometry is defined as
S
(
q
Λ
ω
,
Λ
ω
/d, ητ
ω
, p
ω
) =
κ
ω,
BTE
κ
ω,
Fourier
(22)
where,
κ
ω
=
j
x,ω
/
̄
T
is the conductance per phonon mode and
j
x,ω
is the thickness-averaged
in-plane heat flux defined in equation 16. In figures 3 (b) and 3 (c), t
he magnitude of the
suppression function at
η
= 0 is plotted against phonon MFP non-dimensionalized with
respect to the grating period
q
. The suppression functions from the complete BTE solution
and the models FS I and FS II are identical at high temperatures and
long grating periods,
when the transport is primarily diffusive, governed by the Fourier’s la
w of heat conduction.
However, for low temperatures and short grating periods, FS I un
derpredicts the heat flux
and FS II overpredicts the heat flux carried by phonons with very lo
ng MFPs. Moreover,
the difference between the models FS I and FS II, and the BTE solutio
n is smaller for
thinner films indicating that enhanced boundary scattering in thinne
r films delays the onset
of quasiballistic heat conduction. These observations emphasize th
e importance of using the
complete BTE solution to accurately investigate boundary scatter
ing when grating periods
are comparable to phonon MFPs.
19
(a)
(b)
(c)
T = 100 K
λ = 1 μm
T = 500 K
λ
=
100 μm
d = 2 μm
d = 500 nm
d = 100 nm
d = 10 nm
T = 100 K
d = 10 nm
d = 100 nm
d = 10 nm
d = 100 nm
FIG. 3. Comparison of the thermal conductivity (a) and the su
ppression functions ((b) and (c))
calculated from the models FS I, FS II and by solving the BTE fo
r non-thermalizing diffuse bound-
ary conditions at different temperatures, grating periods (
λ
) and film thicknesses. In figures (b)
and (c), the symbols correspond to the degenerate kernel sol
ution, the solid lines correspond to
FS I model and the dashed solid lines correspond to FS II model
. For very thin films and long
grating periods, the models FS I and FS II are in good agreemen
t with the BTE predictions. For
thicker films and shorter grating periods, FS I underpredict
s and FS II overpredicts the thermal
conductivity at short grating periods (a) and the contribut
ion of phonons with long MFP ((b) and
(c)) compared to the complete BTE solution.
IV. CONCLUSION
We have studied the effect of thermalizing and non-thermalizing boun
dary scattering of
phonons in steady state and transient heat conduction along thin fi
lms by solving the BTE
using analytical and computationally efficient semi-analytical techniq
ues. From our analysis,
we reach the following conclusions. First, under steady state tran
sport conditions, we find
that the thermal transport is governed by the Fuchs-Sondheime
r theory and is insensitive
to whether the boundaries are thermalizing or not. In contrast, u
nder transient conditions,
the decay rates are significantly different for thermalizing and non-
thermalizing walls and
the Fuchs-Sondheimer theory is only applicable in the heat diffusion re
gime. We also show
that, for transient transport, the thermalizing wall boundary co
ndition is unphysical due to
20