PHYSICAL REVIEW B
95
, 205423 (2017)
Experimental metrology to obtain thermal phonon transmission coefficients at solid interfaces
Chengyun Hua,
1
Xiangwen Chen,
2
Navaneetha K. Ravichandran,
3
and Austin J. Minnich
2
,
*
1
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
2
California Institute of Technology, Pasadena, California 91125, USA
3
Boston College, Chestnut Hill, Massachusetts 02467, USA
(Received 23 December 2016; revised manuscript received 17 March 2017; published 17 May 2017)
Interfaces play an essential role in phonon-mediated heat conduction in solids, impacting applications ranging
from thermoelectric waste heat recovery to heat dissipation in electronics. From the microscopic perspective,
interfacial phonon transport is described by transmission coefficients that link vibrational modes in the materials
composing the interface. However, direct experimental determination of these coefficients is challenging because
most experiments provide a mode-averaged interface conductance that obscures the microscopic detail. Here, we
report a metrology to extract thermal phonon transmission coefficients at solid interfaces using
ab initio
phonon
transport modeling and a thermal characterization technique, time-domain thermoreflectance. In combination
with transmission electron microscopy characterization of the interface, our approach allows us to link the atomic
structure of an interface to the spectral content of the heat crossing it. Our work provides a useful perspective on
the microscopic processes governing interfacial heat conduction.
DOI:
10.1103/PhysRevB.95.205423
I. INTRODUCTION
Interfaces play an essential role in phonon-mediated heat
conduction in solids [
1
,
2
]. Material discontinuities lead to ther-
mal phonon reflections that are manifested on a macroscopic
scale as a thermal boundary resistance (TBR), also called
Kapitza resistance
R
k
, that relates the temperature drop at the
interface to the heat flux flowing across it. TBR exists at the
interface between any dissimilar materials due to differences
in phonon states on each side of the interface [
3
]. Typical
interfaces often possess defects or roughness which can lead
to additional phonon reflections and hence higher TBR.
TBR plays an increasingly important role in applications,
particularly as device sizes decrease below the intrinsic mean
free paths (MFPs) of thermal phonons [
2
]. At sufficiently
small length scales, TBR can dominate the total thermal
resistance. For instance, the effective thermal conductivity
of a superlattice can be orders of magnitude smaller than
that of the constituent materials due to TBR [
4
–
7
]. This
physical effect has been used to realize thermoelectrics with
high efficiency [
8
,
9
] and dense solids with exceptionally
low thermal conductivity [
10
]. On the other hand, TBR
can lead to significant thermal management problems [
11
–
13
] in applications such as LEDs [
14
,
15
] and high power
electronics [
13
,
16
].
Numerous works over several decades have investigated the
microscopic origin of TBR at solid-solid interfaces, starting
with studies performed at low temperatures (
∼
1 K), in which
heat is carried predominantly by phonons with frequencies
less than 1 THz [
17
,
18
]. At these low temperatures and for
pristine, ordered interfaces, transmission coefficients can be
obtained from continuum elastic theory in an analogy with
Snell’s law for light; this model is known as the acoustic
mismatch model (AMM) [
19
,
20
]. The AMM was shown to
explain the experimentally measured values of TBR at various
*
Corresponding author: aminnich@caltech.edu
solid-solid interfaces [
18
]. At higher temperatures (above 1 K),
heat transport across the interfaces was found to be sensitive
to surface roughness. For the limit of completely diffuse
scattering in which transmitted and reflected phonons cannot
be distinguished, Swartz constructed the diffuse mismatch
model (DMM) [
1
]. Despite the success of these models at
explaining TBR at low temperatures, they generally fail at
temperatures larger than 40 K and are unable to account for
the atomistic structure of the interface.
Recent works have focused on remedying these deficien-
cies. Optical methods enable the routine measurement of TBR
over a wide range of temperatures for various metal-dielectric
interfaces [
21
–
25
] as well as at interfaces with variable
bonding strength [
26
,
27
]. Other works have examined the
temperature dependence of the thermal conductivity [
28
]in
nanocrystalline samples. Computational atomistic methods
such as molecular dynamics [
29
–
37
] and atomistic Green’s
functions [
38
–
42
] have been extensively applied to obtain the
transmission coefficients at interfaces with realistic atomic
structure. These calculations generally predict the coefficients
to decrease with increasing phonon frequency due to reflec-
tions of short wavelength phonons by atomistic roughness, a
trend that is supported by basic wave physics and indirectly
by experiment [
28
,
43
]. However, a direct determination of the
spectral transmission coefficients at an actual interface has not
yet been reported.
Here, we report a metrology to extract the thermal phonon
transmission coefficients at a solid interface. Our approach,
based on combining experimental observations with
ab initio
phonon transport modeling, exploits quasiballistic transport
near the interface to significantly narrow the possible trans-
mission coefficient profiles at a solid interface compared to the
bounds obtained from traditional approaches. Applying our ap-
proach in conjunction with transmission electron microscopy
(TEM), we are able to directly link atomic structure to the
spectral content of heat crossing the interface. Our approach is
a useful tool to elucidate the microscopic transport properties
of thermal phonons at solid interfaces.
2469-9950/2017/95(20)/205423(21)
205423-1
©2017 American Physical Society
HUA, CHEN, RAVICHANDRAN, AND MINNICH
PHYSICAL REVIEW B
95
, 205423 (2017)
δ
t
i
e
Q
0
(a)
(b)
FIG. 1. (a) Schematic of the principle underlying the measurement of transmission coefficients. If the characteristic length scale of the
thermal transport is much longer than the phonon MFPs, information about the interfacial distribution is lost due to strong scattering. If some
MFPs are comparable to the thermal length scale, the nonequilibrium distribution at the interface propagates into the substrate where it can be
detected. (b) 2D schematic of the experimental configuration subject to a modulated heating source: a double layer structure of a transducer
film on a substrate (sample).
Q
0
is the amplitude of the heating source,
η
is angular modulation frequency,
δ
is the optical penetration depth of
the heating source, and
x
is the cross-plane transport direction.
x
1
and
x
2
correspond to the coordinate systems used in transducer and substrate
accordingly.
II. OVERVIEW OF APPROACH
Our approach is based on interpreting data from the TDTR
experiment with
ab initio
phonon transport model. Briefly,
TDTR is a widely used optical pump-probe technique that is
used to characterize thermal properties. In this experiment,
a sample consists of a metal transducer film on a substrate.
A pulsed laser beam from an ultrafast oscillator is split into a
pump and a probe beam. The pump pulse train is modulated at a
frequency from 1 to 15 MHz to enable lock-in detection, and is
then used to impulsively heat the metal film coated on the sam-
ple. The transient temperature decay
Z
(
t
) at the surface is de-
tected as a change in optical reflectance by the probe beam [
44
].
In the traditional TDTR approach, this transient signal
is related to the desired thermal properties by a macro-
scopic transfer function based on a multilayer heat diffusion
model [
2
,
45
]. This function maps thermal properties such as
substrate thermal conductivity and metal-substrate interface
conductance to the TDTR signal, and thus these properties
are obtained by varying these parameters until the simulated
results match the measured data sets. Put another way, one
must solve an inverse problem that links the data sets to the
unknown parameters; this calculation is often performed using
a nonlinear least squares algorithm.
This approach is widely used and has provided important in-
sights into a wide range of metal-semiconductor interfaces. A
drawback, however, is that the microscopic information about
the interface is averaged into a single number, the interface
conductance, obscuring the microscopic detail. Hopkins
et al.
used TDTR measurements on a variety of metal films with
varying phonon cutoff frequencies to extract spectral informa-
tion about phonon transmission [
23
]. However, determining
transmission coefficients is still challenging due to variations
in phonon density of states in each metal.
In this work, we aim to directly extract the transmission
coefficients from TDTR data by replacing the macroscopic
transfer function based on Fourier’s law with a microscopic
transfer function based on
ab initio
phonon transport modeling.
Just as in the traditional approach, we seek to identify the
parameters that best fit the TDTR data sets. Although this
fitting process is much more complex than the traditional
approach, in principle it is the same widely used procedure.
There are several conditions that must be satisfied to suc-
cessfully extract transmission coefficients using this approach.
First, it is essential that part of the nonequilibrium phonon
distribution emerging from the interface propagate into the
substrate ballistically. As illustrated in Fig.
1(a)
, when MFPs
are much shorter than the characteristic length scale of the
thermal gradient, information about the phonon distribution at
the interface is lost due to scattering. In this case, the TDTR
signal is largely insensitive to the transmission coefficient
profile so long as the overall interface conductance remains
fixed.
On the other hand, if some phonons have sufficiently long
MFPs, the nonequilibrium phonon distribution penetrates into
the substrate. In this case, as we will show in subsequent
sections, the TDTR signal depends not only on the magnitude
of the interface conductance but on the spectral profile of
the transmission coefficients. It is this sensitivity that we will
exploit to retrieve the transmission coefficients from the TDTR
data sets.
This discussion implies that not every substrate will be
suitable for our approach as sufficiently long phonons MFPs
compared to the induced thermal gradient are required.
Fortunately, many experimental reports have demonstrated
clear evidence of this quasiballistic heat transport regime
in different material systems [
46
–
51
]. In particular, MFPs
in Si are generally accepted to exceed one micron at room
temperature [
52
]. Considering that thermal penetration depths
in Si are on the same order in TDTR, Si is a suitable substrate
for this work.
Second, we must determine the microscopic transfer func-
tion that maps the transmission coefficients directly to the
TDTR signal without any artificial fitting parameters. This step
is challenging due to the difficulty of efficiently solving the
BTE for the TDTR experiment. A number of simplified models
[
51
,
53
–
58
] have been proposed, but these models make various
approximations that limit their predictive capability. In this
work, we overcome this challenge using two recent advances
we reported for rigorously solving the spectral BTE under the
relaxation time approximation (RTA) that yield a factor of 10
4
speedup compared to existing methods and allows the first
ab
initio
phonon transport modeling of TDTR free of artificial
205423-2
EXPERIMENTAL METROLOGY TO OBTAIN THERMAL . . .
PHYSICAL REVIEW B
95
, 205423 (2017)
parameters or simplifications of the phonon dispersion. The
solution of the BTE in the substrate uses a Green’s function
method for a semi-infinite domain [
59
], and the thin film
is treated with a cosine expansion approach for a finite
domain [
60
].
Before deriving the necessary functions, we briefly summa-
rize the necessary steps in our approach. First, we obtain TDTR
data for a given metal/substrate sample. We then calculate the
intrinsic phonon dispersion and lifetimes for each material that
are then inserted into the
ab initio
transfer function derived
in the next section. Finally, the transmission coefficients are
varied until the optimized transmission profiles are identified
that best match the experimental data. Thus the procedure
is identical to that used in traditional TDTR measurements
excepting the complications of deriving and inverting the
microscopic transfer function.
III. DERIVATION OF TRANSFER FUNCTION
We now describe the derivation of the microscopic transfer
function
Z
(
T
12
(
ω
)
,t
) that maps transmission coefficients to
the TDTR amplitude and phase data as a function of time. The
result of the derivation is a function for which the only inputs
are the phonon dispersions and lifetimes for each material
composing the interface, and the only unknown parameters
are the spectral transmission coefficients
T
12
(
ω
). The output
of the function is the TDTR amplitude and phase signal versus
delay time. For this work, the phonon dispersion and lifetimes
for Si were calculated from first principles with no adjustable
parameters by Lucas Lindsay [
61
]. We have also used first
principles inputs from Mingo
et al.
[
62
] with no appreciable
difference in results or conclusions.
Let us first briefly review the signal formation in TDTR.
Since the thermal response given by the BTE is a linear time-
invariant system, the output transient signal
Z
(
t
)ofTDTRcan
be represented in terms of frequency response solution through
the following equation [
45
]:
Z
(
T
12
(
ω
)
,t
)
=
∞
∑
n
=−∞
H
(
T
12
(
ω
)
,η
0
+
nη
s
)
e
inη
s
t
,
(1)
where
η
0
is the reference angular frequency of the periodic
heating,
η
s
is the angular sampling frequency set by the
repetition rate of the laser pulses, and
H
(
T
12
(
ω
)
,η
) is surface
temperature response subject to a periodic heating at frequency
η
given the transmission coefficients
T
12
(
ω
). Therefore iden-
tifying
Z
(
T
12
(
ω
)
,t
) requires computing
H
(
T
12
(
ω
)
,η
), or the
surface temperature frequency response to a periodic heating
in a double-layer structure of a metal film on a substrate as
shown in Fig.
1(b)
.
Thermal transport in an isotropic crystal, assuming only
cross-plane heat conduction, is described by the one-
dimensional (1D) spectral Boltzmann transport equation
(BTE) under the relaxation time approximation (RTA) [
59
,
63
],
∂g
ω
∂t
+
μv
ω
∂g
ω
∂x
=−
g
ω
+
f
0
(
T
0
)
−
f
0
(
T
)
τ
ω
+
Q
ω
4
π
,
(2)
f
0
(
T
)
=
1
4
π
̄
hωD
(
ω
)
f
BE
(
T
)
≈
f
0
(
T
0
)
+
1
4
π
C
ω
T ,
(3)
where
g
ω
=
̄
hωD
(
ω
)[
f
ω
(
x,t,μ
)
−
f
0
(
T
0
)] is the deviational
distribution function,
f
0
is the equilibrium Bose-Einstein (BE)
distribution function,
μ
=
cos(
θ
) is the directional cosine,
v
ω
is the phonon group velocity,
τ
ω
is the phonon relaxation
time, and
Q
ω
(
x,t
) is the spectral volumetric heat generation.
Assuming a small temperature rise,
T
=
T
−
T
0
, relative
to a reference temperature,
T
0
, the equilibrium distribution
is proportional to
T
,asshowninEq.(
3
). Here, ̄
h
is the
reduced Planck constant,
ω
is the phonon frequency,
D
(
ω
)
is the phonon density of states,
f
BE
is the Bose-Einstein
distribution, and
C
ω
=
̄
hωD
(
ω
)
∂f
BE
∂T
is the mode specific heat.
The volumetric heat capacity is then given by
C
=
∫
ω
m
0
C
ω
dω
and the Fourier thermal conductivity
k
=
∫
ω
m
0
k
ω
dω
, where
k
ω
=
1
3
C
ω
v
ω
ω
and
ω
=
τ
ω
v
ω
is the phonon MFP. To close
the problem, energy conservation is used to relate
g
ω
to
T
,
given by
∫∫
ω
m
0
[
g
ω
(
x,t
)
τ
ω
−
1
4
π
C
ω
τ
ω
T
(
x,t
)
]
dωd
=
0
,
(4)
where
is the solid angle in spherical coordinates and
ω
m
is the cutoff frequency. Note that summation over phonon
branches is implied without an explicit summation sign
whenever an integration over phonon frequency or MFP is
performed.
We now divide our discussion into three parts: transducer
film, substrate, and interface. The BTE in the transducer film
can be reformulated as a Fredholm integral equation of the
second kind [
60
]; the solution in the substrate can be treated as
the Green’s function to the BTE [
59
]. The solutions in the two
layers depend on each other through the interface conditions
that enforce conservation of heat flux.
A. Transducer film
The metal thin film serves as an optical transducer that
absorbs the incident optical energy while also enabling the ob-
servation of temperature decay through the thermoreflectance
coefficient. In our work, we neglect electrons and consider
that heat is only carried by phonons in Al. Justification for this
approximation is given in Appendix
B
.
Since the system is modulated at a given frequency
η
,we
can assume that both
g
1
ω
and
T
1
are of the form
e
iηt
to define
g
1
ω
=
̃
g
1
ω
(
x
1
,μ
)
e
iηt
and
T
1
=
̃
T
1
(
x
1
)
e
iηt
. The volumetric
heat generation in thin film is given by
Q
ω
=
Q
0
ω
e
iηt
e
−
x
1
/δ
,
where the amplitude of heating source
Q
0
=
∫
ω
m
0
Q
0
ω
dω
.We
also assume that phonons are specularly reflected at
x
1
=
0,
i.e.
̃
g
1
ω
(
x
1
=
0
,μ
)
=
̃
g
1
ω
(
x
1
=
0
,
−
μ
). Substituting the definition
of
̃
g
1
ω
and
̃
T
1
(
x
1
) and the specular boundary condition at
x
1
=
0 into Eq. (
2
) leads to a first-order ODE for
̃
g
1
ω
(
x
1
,μ
). Its
solution is given by
̃
g
+
1
ω
(
x
1
,μ
)
=
B
ω
e
−
γ
1
ω
μ
(
d
+
x
1
)
+
∫
d
0
C
1
ω
̃
T
(
x
1
)
+
Q
0
ω
e
−
x/δ
τ
1
ω
4
π
1
ω
μ
×
e
γ
1
ω
μ
(
x
1
+
x
1
)
dx
1
+
∫
x
1
0
C
1
ω
̃
T
(
x
1
)
+
Q
0
ω
e
−
x/δ
τ
1
ω
4
π
1
ω
μ
×
e
γ
1
ω
μ
(
x
1
−
x
1
)
dx
1
(
μ
∈
(0
,
1])
,
(5)
205423-3
HUA, CHEN, RAVICHANDRAN, AND MINNICH
PHYSICAL REVIEW B
95
, 205423 (2017)
̃
g
−
1
ω
(
x
1
,μ
)
=
B
ω
e
γ
1
ω
μ
(
d
−
x
1
)
−
∫
d
x
C
1
ω
̃
T
(
x
1
)
+
Q
0
ω
e
−
x/δ
τ
1
ω
4
π
1
ω
μ
×
e
γ
1
ω
μ
(
x
1
−
x
1
)
dx
1
(
μ
∈
[
−
1
,
0])
,
(6)
where
γ
1
ω
=
(1
+
iητ
1
ω
)
/
1
ω
,
d
is the film thickness, and
B
ω
are the unknown coefficients determined by the interface
condition at
x
1
=
d
. Here,
̃
g
+
1
ω
(
x
1
,μ
) indicates the forward-
going phonons and
̃
g
−
1
ω
(
x
1
,μ
) the backward-going phonons.
To close the problem, we plug Eqs. (
5
) and (
6
) into Eq. (
4
)
and obtain an integral equation for temperature as
̃
T
(
̂
x
1
)
−
∫
1
0
̃
T
(
̂
x
1
)
K
(
̂
x
1
,
̂
x
1
)
d
̂
x
1
=
∫
ω
m
0
B
ω
F
1
ω
(
̂
x
1
)
dω
+
F
2
(
̂
x
1
)
,
(7)
where
̂
x
1
=
x
1
/d
. The kernel function
K
(
̂
x
1
,
̂
x
)isgivenby
K
(
̂
x
1
,
̂
x
1
)
=
1
2
∫
ω
m
0
C
1
ω
τ
1
ω
dω
∫
ω
m
0
C
1
ω
τ
1
ω
Kn
1
ω
{
E
1
[
̂
γ
1
ω
(
̂
x
1
+
̂
x
1
)]
+
E
1
[
̂
γ
1
ω
|
̂
x
1
−
̂
x
1
|
]
}
dω
(8)
and the two inhomogeneous functions are given by
F
1
ω
(
̂
x
1
)
=
2
π
∫
ω
m
0
C
1
ω
τ
1
ω
dω
1
τ
1
ω
{
E
2
[
̂
γ
1
ω
(1
+
̂
x
1
)]
+
E
2
[
̂
γ
1
ω
(1
−
̂
x
1
)]
}
,
(9)
F
2
(
̂
x
1
)
=
2
π
∫
ω
m
0
C
1
ω
τ
1
ω
dω
∫
1
0
∫
ω
m
0
Q
0
ω
e
−
ρ
̂
x
1
Kn
1
ω
{
E
1
[
̂
γ
1
ω
(
̂
x
1
+
̂
x
1
)]
+
E
1
[
̂
γ
1
ω
|
̂
x
1
−
̂
x
1
|
]
}
dωd
̂
x
,
(10)
where Kn
1
ω
=
1
ω
/d
is the Knudsen number,
̂
γ
1
ω
=
1
+
iητ
1
ω
Kn
1
ω
,
and
E
n
(
x
) is the exponential integral given by [
64
]
E
n
(
x
)
=
∫
1
0
μ
n
−
2
e
−
x
μ
dμ.
(11)
Recently, we have developed a spectral method to efficiently
solve Eq. (
7
)inRef.[
60
]. Briefly, the functions in Eq. (
7
) can
be expanded as a finite cosine series, such as
̃
T
1(
N
)
(
̂
x
1
)
≈
N
∑
n
=
0
t
n
cos(
nπ
̂
x
1
)
(12)
and
K
(
N
)
(
̂
x,
̂
x
)
=
1
4
k
00
+
1
2
N
∑
m
=
1
k
m
0
cos(
mπ
̂
x
)
+
1
2
N
∑
n
=
1
k
0
n
cos(
nπ
̂
x
)
,
+
N
∑
m
=
1
N
∑
n
=
1
k
mn
cos(
mπ
̂
x
) cos(
nπ
̂
x
)
,
(13)
where
N
is the truncated basis number, and
t
n
’s and
k
nm
’s are
the Fourier coefficient. Similarly,
F
1
ω
(
̂
x
1
) and
F
2
(
̂
x
1
)arealso
expanded in term of cosines. Following the steps in the above
reference, we can express the temperature as
̃
T
1
(
̂
x
1
)
=
[
A
−
1
(
f
1
B
+
f
2
)]
T
φ
(
x
)
,
(14)
where
φ
(
x
)
=
[1 cos(
π
̂
x
) cos(2
π
̂
x
)
···
cos(
Nπ
̂
x
)]
T
, and
the matrix
A
contains elements
A
00
=
1
−
k
00
4
,
A
0
n
=−
1
2
k
0
n
,
A
n
0
=−
k
n
0
4
,
A
nn
=
1
−
1
2
k
nn
and
A
nm
=−
1
2
k
nm
(
m
=
n
=
0).
B
is a
N
ω
column vector of the unknown coefficients
B
ω
,
where
N
ω
is the number of discretization in phonon frequency.
f
1
is a
N
×
N
ω
matrix, consisting of the Fourier coefficients
of
F
1
ω
i
(
̂
x
1
) evaluated at each phonon frequency
ω
i
, and
f
2
is
a N column vector, consisting of the Fourier coefficients of
F
2
(
̂
x
1
). Then,
̃
g
+
1
ω
(
x
1
,μ
) and
̃
g
−
1
ω
(
x
1
,μ
) can be expressed in
terms of the unknown coefficients
B
ω
by plugging Eq. (
14
)
into Eqs. (
5
) and (
6
).
B. Substrate
The substrate can be treated as a semi-infinite region subject
to a surface heat flux. Therefore the BTE for the substrate
becomes
iη
̃
g
2
ω
+
μv
2
ω
∂
̃
g
2
ω
∂x
2
=−
̃
g
2
ω
τ
2
ω
+
C
2
ω
4
πτ
2
ω
̃
T
(
x
2
)
+
1
2
P
ω
v
2
ω
|
μ
|
δ
(
x
2
)
,
(15)
where the unknown coefficients
P
ω
’s are determined through
the interface conditions.
We then apply the Green’s function method given in
Ref. [
59
]. The unknown distribution function in spatial
frequency domain is then written as
̃
g
2
ω
(
η,ξ
2
)
=
C
2
ω
4
π
̃
T
2
(
η,ξ
2
)
+
1
2
P
ω
2
ω
|
μ
|
/C
2
ω
1
+
iητ
2
ω
+
iμξ
2
2
ω
,
(16)
and the temperature profile
̃
T
2
(
η,ξ
2
)
=
∫
ω
m
0
P
ω
v
2
ω
1
+
iητ
2
ω
(
2
ω
ξ
2
)
2
log
[
1
+
(
2
ω
ξ
2
1
+
iητ
2
ω
)
2
]
dω
∫
ω
m
0
C
2
ω
2
πτ
2
ω
[
1
−
1
2
ω
ξ
2
tan
−
1
(
ω
ξ
2
1
+
iητ
2
ω
)]
dω
,
(17)
where
ξ
2
is the Fourier variable of
x
2
. Again, to express
̃
g
2
ω
only in terms of unknown coefficients
P
ω
, we simply plug
Eq. (
17
) into Eq. (
16
).
C. Interface condition
The unknown coefficients in the solutions of transducer film
and substrate are obtained by applying appropriate interface
conditions. Here, we use the elastic transmission interface
condition with mode conversion, closely following the work by
Minnich
et al.
[
54
]. Briefly, for a given mode
i
, the heat fluxes
outgoing from the interface,
q
i
−
1
ω
and
q
i
+
2
ω
, must be equal to the
reflected and transmitted heat fluxes incident to the interface,
q
i
+
1
ω
and
q
i
−
2
ω
. By assuming elastic and diffuse scattering, the
transmission and reflection process for each phonon frequency
is treated independently and the heat flux equality condition
must be satisfied for each frequency and polarization.
205423-4
EXPERIMENTAL METROLOGY TO OBTAIN THERMAL . . .
PHYSICAL REVIEW B
95
, 205423 (2017)
The interface conditions are
∫
1
0
g
i
+
2
ω
v
i
2
ω
μdμ
=
∑
j
T
ji
12
(
ω
)
∫
1
0
g
j
+
1
ω
v
j
1
ω
μdμ
+
∑
j
R
ji
21
(
ω
)
∫
1
0
g
j
−
2
ω
v
j
2
ω
μdμ,
(18)
∫
1
0
g
i
−
1
ω
v
i
1
ω
μdμ
=
∑
j
T
ji
21
(
ω
)
∫
1
0
g
j
−
2
ω
v
j
2
ω
μdμ
+
∑
j
R
ji
12
(
ω
)
∫
1
0
g
j
+
1
ω
v
j
1
ω
μdμ,
(19)
where
T
ji
12
(
ω
) is the transmission coefficient of mode
j
at
frequency
ω
from side 1 to side 2 as mode
i
,
R
ji
21
(
ω
)isthe
reflection coefficient of mode
j
at frequency
ω
from side 2
back into side 2 as mode
i
, and so on.
The next question is how
T
ij
12
(
ω
) is related to the other
reflection and transmission coefficients. The reflection coef-
ficients are related to the transmission coefficients by energy
conservation given by
∑
j
T
ij
12
(
ω
)
+
R
ij
12
(
ω
)
=
1
(20)
and
∑
j
T
ij
21
(
ω
)
+
R
ij
21
(
ω
)
=
1
.
(21)
T
ji
21
(
ω
) is related to
T
ij
12
(
ω
) through the principle of detailed
balance, which requires that no net heat flux can transmit
across the interface when both materials are at an equilibrium
temperature
T
. Applying this condition to every phonon
mode on each side of the interface for each polarization and
frequency gives
T
ij
12
(
ω
)
C
i
1
ω
v
i
1
ω
=
T
ji
21
(
ω
)
C
j
2
ω
v
j
2
ω
.
(22)
Therefore we need to specify
T
ij
12
(
ω
),
R
ij
12
(
ω
), and
R
ij
21
(
ω
).
Let us first consider a special case where no mode
conversion is allowed (
T
ij
12
(
ω
),
T
ij
21
(
ω
),
R
ij
12
(
ω
),
R
ij
21
(
ω
)
=
0
for
i
=
j
). Then, the interface conditions become
∫
1
0
g
i
+
2
ω
v
2
ω
μdμ
=
T
ii
12
(
ω
)
∫
1
0
g
i
+
1
ω
v
i
1
ω
μdμ
+
R
ii
21
(
ω
)
∫
1
0
g
i
−
2
ω
v
i
2
ω
μdμ,
(23)
∫
1
0
g
i
−
1
ω
v
1
ω
μdμ
=
T
ii
21
(
ω
)
∫
1
0
g
i
−
2
ω
v
i
2
ω
μdμ
+
R
ii
12
(
ω
)
∫
1
0
g
i
+
1
ω
v
i
1
ω
μdμ,
(24)
and the detail balance becomes
T
ii
12
(
ω
)
C
i
1
ω
v
i
1
ω
=
T
ii
21
(
ω
)
C
i
2
ω
v
i
2
ω
.
(25)
Therefore, once
T
ii
12
(
ω
) is specified, all the other transmission
and reflection coefficients are determined. For now, we only
consider this special case and neglect the mode conversion in
our BTE simulations. Later, we show that the mode specific
transmission coefficients cannot be resolved by the TDTR
measurements and the measurable quantity is
∑
j
T
ij
12
(
ω
)
instead of individual transmission coefficients. For simplicity,
we will use
T
12
(
ω,p
) rather than the summation.
D. Justification of approximations
In this section, we justify the approximations made in the
derivation of the transfer function. First, we have assumed
only cross-plane transport since the pump diameter in the
experiments is 60
μ
m. We verify this assumption by extending
our model to include in-plane transport and comparing the full
3D TDTR signal to the 1D calculations. We observe no appre-
ciable deviations. Further details are given in Appendix
A
.
Second, we have neglected the role of electrons, either
in carrying heat in the metal or in carrying heat across the
interface by direct coupling to substrate phonons. For the first
point, in Appendix
B
, we solve the BTE explicitly including
electron-phonon coupling and show that it has negligible effect
on the TDTR signal and hence the transmission coefficients.
For the second point, our approach cannot rule out this
mechanism. However, the available evidence in the literature
[
21
,
65
,
66
] shows that the metal electron-substrate phonon
coupling effect is too small to be observed for several material
systems. Therefore little evidence exists to support a large
contribution from this mechanism. As a result, we assume that
heat is carried across the interface solely by phonons.
The third assumption is the neglect of mode conversion,
or the change of phonon polarization as phonons transmit
or reflect at the interface. In Appendix
C
, we explicitly
calculate the transfer function including this mechanism in
a Al/Si system, showing that it has very little impact on
the results. Therefore fitting the data with or without mode
conversion yields identical results in Al/Si systems. Note that
this assumption must be checked for other material systems.
E. Overview of calculation
We now describe how these pieces fit together to provide
the transfer function
Z
(
T
12
(
ω
)
,t
). At any given frequency
η
, the analytical expression of the unknown distributions
from both sides,
g
±
1
ω
and
g
±
2
ω
can be obtained according to
Secs.
III A
and
III B
, which then are evaluated at
x
1
=
d
and at
x
2
=
0, respectively. Given the values of
T
12
(
ω
), the values of
R
12
(
ω
),
T
21
(
ω
) and
R
21
(
ω
) can be inferred, and the unknown
coefficients
P
ω
and
B
ω
’s are obtained by plugging Eqs. (
5
),
(
6
), and (
16
) into Eqs. (
23
) and (
24
) and solving the linear
system.
Once
B
ω
is known, the temperature profile in the transducer
film is solved using Eq. (
14
) as well as the surface temperature
H
(
η
). At a given modulation frequency
ω
0
,
η
is chosen to be
ω
0
+
nω
L
, where integer
n
is typically ranging from
−
50 to 50,
and
ω
L
is the laser repetition frequency. Using Eq. (
1
) yields
Z
(
T
12
(
ω
)
,t
), which is directly compared to experimental data
Z
exp
(
t
). The values of all the constants used in the calculations
are tabulated in Table
I
.
IV. SENSITIVITY OF TDTR SIGNAL
TO TRANSMISSION COEFFICIENTS
With the transfer function obtained, it is useful to re-
examine the conditions that are required to successfully extract
205423-5
HUA, CHEN, RAVICHANDRAN, AND MINNICH
PHYSICAL REVIEW B
95
, 205423 (2017)
TABLE I. All the constants appearing in the BTE models and the
fitting process are given in the following table.
Bulk thermal properties
Al heat capacity (J m
−
3
K
−
1
)2
.
41
×
10
6
Al lattice thermal conductivity (W m
−
1
K
−
1
)
123
Al total thermal conductivity (W m
−
1
K
−
1
)
230
Si heat capacity (J m
−
3
K
−
1
)1
.
63
×
10
6
Si thermal conductivity (W m
−
1
K
−
1
)
155
SiGe heat capacity (J m
−
3
K
−
1
)1
.
63
×
10
6
SiGe thermal conductivity (W m
−
1
K
−
1
)51
Electronic thermal properties in Al
Heat capacity (J m
−
3
K
−
1
)4
.
11
×
10
4
Thermal conductivity (W m
−
1
K
−
1
)
203
Electron-phonon coupling coefficient
g
(W m
−
3
K
−
1
)2
.
1
×
10
17
Constants in Tamura formula
Volume per Si atom
V
0
(nm
3
)
0.02
Measure of the mass disorder
m
0
0.0568
Transducer film thickness
Al/Si with a clean interface (nm)
69
Al/SiGe with a clean interface (nm)
72
Al/Si with a native oxidized interface (nm)
70
Al/Si with a thermally grown oxidized interface (nm)
70
Other constants
Optical penetration depth
δ
(nm)
10
Laser repetition frequency (MHz)
76
transmission coefficients from TDTR data. First, the number
of unknowns should be on the order of the number of data
points. As an estimate, consider that we take TDTR data at
five modulation frequencies. As in Eq. (
1
), in the frequency
domain, each TDTR signal consists of surface temperature
responses at different frequencies, the lowest of which is
the modulation frequency. Typically, among those surface
temperature responses, the responses at the lowest four to five
frequencies contain most of the thermal information [
67
]. With
around five TDTR data sets containing amplitude and phases at
five different modulation frequencies, we obtain 40–50 unique
data points.
Now consider the number of unknowns. If we take the
materials to be isotropic, the transmission coefficients depend
only on phonon frequency. In our typical discretization, we
find that there around one hundred transmission coefficients
that must be fit. However, these coefficients are not all
independent due to a smoothness constraint—physically, the
transmission profile cannot fluctuate arbitrarily, an intuition
supported by atomistic simulations [
68
,
69
]. Qualitatively,
this requirement implies that each coefficient depends on
the adjacent coefficients, decreasing the effective number of
unknowns by a value on the order of two to three. Therefore
the number of unknowns is comparable to the number of data
points.
Ample data points are a necessary but not sufficient
condition to enable extraction of transmission coefficients. The
last requirement is that the TDTR signal should be sensitive
to the shape of the transmission coefficient profile. In the case
of heat diffusion, this requirement is not satisfied: so long
as the interface conductance is unchanged, the TDTR signal
will not change because scattering in the substrate obscures
the interfacial phonon distribution. To demonstrate this point,
we simulated TDTR signals using two different transmission
coefficient profiles as shown in Fig.
2(a)
with the scattering
rates of Si modified such that no MFP exceeds 50 nm.
The chosen transmission profiles possess the same interface
conductance. The amplitude and phase of the simulated TDTR
signals for these two transmission coefficient profiles are
shown in Fig.
2(b)
. The figure shows that they are nearly
identical even though the two transmission coefficient profiles
are completely different, demonstrating that in the diffusion
regime the TDTR signal only depends on the magnitude of
interface conductance.
On the other hand, in the quasiballistic regime the TDTR
signal is sensitive to the shape of the transmission coefficient
profile, enabling the coefficients to be obtained with tight
constraints. We demonstrate this sensitivity in Fig.
2(c)
.For
this calculation, the same transmission coefficient profiles
Time (ns)
2468
Phase (deg)
-55
-50
-45
-40
-35
-30
-25
Increasing
Constant
Mod. freq: 1.43 MHz
Substrate MFPs <= 50 nm
Time (ns)
2468
Phase (deg)
-55
-50
-45
-40
-35
-30
-25
-20
Constant
Increasing
Exp.
Mod. freq: 1.43 MHz
Al/Si with a clean interface at 300 K
Phonon frequency (THz)
0246810
Transmission coefficient
0.0
0.2
0.4
0.6
0.8
1.0
Constant
Increasing
LA
Si to Al
(a)
(b)
(c)
FIG. 2. (a) Two artificial transmission coefficient profiles from Si to Al for the longitudinal branch. The two profiles have different trend: a
constant value (dashed line) and an increasing trend (dotted-dash line), but both give the same interface conductance. (b) The amplitude of the
simulated TDTR signals using a constant transmission coefficient (dashed line) and a increasing transmission coefficient profile (dotted-dash
line) with silicon’s dispersion and a cutoff scattering rates such that there no MFP exceeding 50 nm. (c) The amplitude of the measured TDTR
signal (solid line) for Al/Si at 300 K and the simulated TDTR signals using the transmission coefficient profiles shown in (a). (Insets) The
phase of the corresponding TDTR signals. The two signals are almost identical under diffusion transport, while, in the quasiballistic regime,
the two simulated TDTR signals are no longer identical and do not match with the measured signal either.
205423-6
EXPERIMENTAL METROLOGY TO OBTAIN THERMAL . . .
PHYSICAL REVIEW B
95
, 205423 (2017)
are used but the bulk scattering rates of Si from Lindsay
et al.
are used without modification. The longest MFP in
silicon is on the order of micrometers [
52
,
62
], and the
characteristic length scale of a TDTR experiment with a
modulation frequency around 1 MHz is also on the order of a
micrometer. Therefore the transport in the TDTR experiment
for Al/Si is quasiballistic. In this case, the TDTR signals using
two different transmission coefficient profiles are no longer
identical as shown in Fig.
2(c)
, demonstrating that in the
quasiballistic regime the TDTR signal is sensitive to both the
magnitude of interface conductance and the spectral profile.
Figure
2(c)
plots a measured TDTR signal. From the plot,
one can immediately tell the transmission coefficient profiles
used in the calculations are not correct as the experimental
and calculation signals do not match. The correct transmission
coefficient profile will reproduce the measured TDTR signal.
Our procedure to identify this profile is described in the next
section.
V. SOLUTION OF INVERSE PROBLEM
The final step is to solve the inverse problem that identifies
the transmission coefficients that best explain the observed
data. From the BTE model, we obtain a surface temperature
decay curve as a function of time just like the one measured in
the experiments. For a given sample, the actual transmission
coefficient profile as a function of phonon frequency will
minimize the difference between the simulation curves and
experimental TDTR traces at all modulation frequencies.
We solve the inverse problem using a particle swarm
optimization (PSO) method to search for the optimal profile.
The goal of the PSO method is to minimize the objective
function defined as
f
=
α
|
g
ab initio
(
T
12
(
ω
))
−
g
measured
|
+
(1
−
α
)
∫
(
d
2
T
12
dω
2
)
2
dω.
(26)
The first part of the equation evaluates the norm of the
difference between the experimentally measured and BTE-
simulated TDTR signals given a transmission profile profile
T
12
(
ω
). The second part of the equation evaluates the second
derivative of the transmission coefficient profile, serving as
the smoothness penalty function. Note that the smoothness of
the profiles is the only constraint we impose in the objective
function. The smoothing parameter
α
determines the relative
importance of the second part to the first part. If
α
=
1, then
no smoothness constraint is imposed. Here, we use
α
=
∫
(
d
2
T
0
12
dω
2
)
2
dω
∣
∣
g
ab initio
(
T
0
12
(
ω
)
)
−
g
measured
∣
∣
,
(27)
where
T
0
12
(
ω
) is the initial profile. The formula is chosen such
that the first and second parts of the equation have the same
order of magnitude.
To search for the optimal profile that minimizes the
objective function, the PSO algorithm randomly initializes a
collection of transmission coefficient profiles and evolves them
in steps throughout the phase space which contains all possible
transmission coefficient profiles. At each step and for each
profile, the algorithm evaluates the objective function defined
as above. After this evaluation, the algorithm decides how each
profile should evolve according to the current best profile. The
profile evolves, then the algorithm reevaluates. The algorithm
stops when the objective function reaches the desired value.
The transmission coefficient profile that achieves the minimum
value of the objective function is the optimal profile that
explains the data.
However, since the inverse problem is ill-posed, a unique
solution does not exist. We generate a probability density
plot for the transmission coefficients using Gibbs sampling
to explore adjacent regions of the optimal transmission
coefficient profile. We first randomly generate around 1000
profiles by perturbing the optimal profile with a smooth
function defined using the following formula:
δ
=
A
[
r
1
cos(2
πω/ω
max
r
2
+
2
πr
3
)
+
r
4
sin(2
πω/ω
max
r
5
+
2
πr
6
)]
,
(28)
where the amplitude of the perturbation
A
is 0.1, and
r
1
,
r
2
,
r
3
,
r
4
,
r
5
, and
r
6
are random numbers between 0 to 1. We
evaluate the objective function at all the perturbed profiles and
recorded the values. Then, we start the Gibbs sampling process.
At each iteration, we randomly draw a profile,
a
,fromthe
stored population and compare the value of its corresponding
objective function
f
n
to the one from the previous step,
f
n
−
1
,
evaluated at profile
b
.If
f
n
is less than
f
n
−
1
, we accept
a
and
kept
f
n
. If not, a random number
r
is drawn and compared to
u
=
p/
(1
+
p
), where
p
=
exp
(
f
n
−
f
n
−
1
T
0
)
.
(29)
If
r
is smaller than
u
, then we accept
a
and kept
f
n
.If
not, we reject
a
and update
f
n
to be
f
n
−
1
. The system
temperature
T
0
is chosen such that the stationary distribution
is gradually changing. Here,
T
0
is set to be the mean value
of the objective functions of all the perturbed samples. We
keep track of how many times each profile was chosen at
each iteration and generated a histogram of the occurrence
frequency of each profile. We stop the sampling process when
the histogram becomes stationary. This occurrence frequency
is also called the likelihood of the transmission coefficient
profiles. The higher the value of a profile’s likelihood is, the
better the fit with the experimentally measured TDTR signals
at different modulation frequencies. Thus, by combining the
PSO method with the Gibbs sampling algorithm, we are able
to determine the most likely transmission coefficients at the
interface between Si and Al.
VI. RESULTS
A. Phonon transmission coefficients
We demonstrate our transmission coefficient measurements
on an Al film on Si substrate with the native oxide removed by
hydrofluoric acid prior to Al deposition, yielding a clean inter-
face. The TEM image in Fig.
3(a)
shows the interface thickness
is less than 0.5 nm. The amplitude and phase of signals
from the lock-in amplifier at different modulation frequencies
are given in Fig.
3
. For reference, solving the usual inverse
205423-7
HUA, CHEN, RAVICHANDRAN, AND MINNICH
PHYSICAL REVIEW B
95
, 205423 (2017)
FIG. 3. Experimental TDTR data (symbols) on this sample at
T
=
300 K for modulation frequencies
f
=
2
.
68, 5.51, and 9.79 MHz along
with the (a) amplitude and (b) phase fit to the data from the BTE simulations (shaded regions), demonstrating excellent agreement between
simulation and experiment. The shaded stripes denoted BTE simulations correspond to the likelihood of the measured transmission coefficients
possessing a certain value as plotted. Inset: TEM image showing the clean interface of an Al/Si sample with the native oxide removed.
The interface thickness is less than 0.5 nm. Transmission coefficients of longitudinal phonons T
Si
→
Al
(
ω
) (blue shaded region) vs (c) phonon
frequency and (d) phonon wavelength, along with the DMM transmission coefficient profile (green dashed line) that gives the same interface
conductance as the measured T
Si
→
Al
(
ω
). The intensity of the shaded region corresponds to the likelihood that the transmission coefficient
possesses a given value.
problem with the macroscopic transfer function on this data
set yields
G
≈
280 MW m
−
2
K
−
1
and
k
≈
140 W m
−
1
K
−
1
,
in good agreement with prior works and literature values for
the thermal conductivity of Si [
43
,
48
]. Although the good
agreement is often taken as evidence that the macroscopic
transfer function is valid for Si, this conclusion is incompatible
with several independent
ab initio
calculations that clearly
show that heat is carried by phonons with MFPs exceeding the
thermal penetration depth of TDTR [
62
,
70
]. This prediction
has recently been experimentally confirmed by Cuffe
et al.
using thermal measurements on variable thickness silicon
membranes [
52
]. This fact implies that quasiballistic transport
should be readily observable in a typical TDTR experiment on
Si, despite the seemingly correct thermal properties measured.
This apparent contradiction is resolved by observing that the
signal measured in TDTR strongly depends on the spectral
profile of the transmission coefficients in the quasiballistic
regime as shown in Fig.
2
.
We represent the transmission coefficient as a probability
density plot, with the color intensity indicating the likelihood
that a single transmission coefficient curve passing through
a particular point at a given phonon frequency is able to
simultaneously explain all of the data in Fig.
3
, without
any other adjustable parameters. The result is shown in
Fig.
3(c)
. The figure shows that the transmission coefficient
from Si to Al for longitudinal phonons, T
Si
→
Al
(
ω
), starts at
unity, its maximum possible value, and decreases steadily
to near zero for high phonon frequencies (
∼
10 THz). The
transmission coefficient profiles for the other polarizations
have similar shapes, and so throughout the paper we plot
only the longitudinal transmission coefficients for simplicity.
The transmission coefficients from Al to Si, T
Al
→
Si
(
ω
)are
calculated by satisfying the principle of detailed balance; the
relationship between T
Si
→
Al
(
ω
) and T
Al
→
Si
(
ω
) reflects the
differences in density of states and group velocity between
the two materials. The transmission coefficients for each side
of the interface and for the other polarizations are given in
Appendix
D
.
B. Comparison with conventional models
Our measured transmission coefficient profile thus indicates
that longitudinal phonons with frequencies less than 6 THz are
transmitted to the maximum extent allowed by the principle of
detailed balance, while longitudinal phonons with frequencies
larger than 8 THz are nearly completely reflected at the
interface. We now examine this result in context with the
common models for transmission coefficients. The AMM
205423-8