of 7
Spectrally resolved specular reflections of thermal phonons from atomically rough
surfaces
Supplementary Material
Navaneetha K. Ravichandran
Department of Physics, Boston College, Chestnut Hill, MA 02467, USA
Hang Zhang
Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing 100190, China
Austin J. Minnich
Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA
(Dated: September 17, 2018)
2
S1. ADDITIONAL TEM IMAGES
In the main manuscript, we showed the representative transmission electron microscope (TEM) images for the
membranes M1 and M2. Figure S1 shows the TEM images for M1 and M2 at a few other locations, demonstrating
that the surface of M1 is smoother than M2.
FIG. S1. Cross-sectional TEM images of the surface of membranes M1 (
A
) and M2 (
B
). The RMS roughness on the crystalline
silicon membrane surface is estimated to be 2.5
±
0.5
̊
A for membrane M1 and 7
±
0.5
̊
A for membrane M2. Scale bar: 20
̊
A.
S2.
AB-INITIO
PROCEDURE TO EXTRACT SPECULARITY PARAMETER
In this section, we describe the
ab-initio
modeling procedure to extract the phonon specularity parameter from
grating period- and temperature-dependent thermal conductivity measurements. As described in the main manuscript,
the measured grating period-dependent thermal conductivity of the free-standing membranes at a given temperature
in the quasiballistic regime is obtained by solving the phonon Boltzmann transport equation as,
κ
(
q
) =
ζ
S
(
q
Λ
ζ
,
Λ
ζ
/d,p
λ
)
[
1
3
C
ζ
v
ζ
Λ
ζ
]
(S1)
where
q
= 2
π/δ
is the grating wavevector corresponding to a grating period
δ
,
d
is the thickness of the membrane, Λ
ζ
is the mean free path,
p
λ
is the wavelength-dependent specularity parameter,
C
ζ
is the volumetric mode specific heat
and
v
ζ
is the group velocity of a phonon mode
ζ
(
k
,j
) denoting a particular phonon wavevector
k
, polarization
j
and wavelength
λ
= 2
π/
|
k
|
. The semi-analytical form of
S
(
q
Λ
ζ
,
Λ
ζ
/d,p
λ
) is described in ref. [47] of the main
manuscript. To interpret our experiments, the solutions of the BTE (equation S1) requires
ab-initio
phonon properties
of silicon projected onto an equivalent isotropic crystal. We use a Gaussian kernel-based regression method to obtain
the isotropic phonon properties for silicon from a complete set of phonon properties in reciprocal space as described
below (full Brillouin zone phonon properties provided by Dr. Lucas Lindsay).
3
A. Gaussian Kernel Regression for Isotropic Phonon Properties
In the equivalent isotropic dataset, all phonon properties are represented only as functions of phonon frequency. In
particular, for any property
I
n
(
k
), where
k
is the wave vector of the reciprocal space and
n
is the mode number, we
are interested in obtaining:
I
n
(
ω
) =
BZ
d
k
(2
π
)
3
V
I
n
(
k
)
δ
(
ω
ω
k
)
k
I
n
(
k
)
w
(
k
)
δ
(
ω
ω
k
)
(S2)
where
V
is the volume of the supercell and
w
(
k
) are the weights on the discrete Brillouin zone grid. To evaluate the
sum/integral in eq. S2, we use the following approximation for the
δ
-function:
δ
(
ω
ω
k
)
1
2
π
exp
(
(
ω
ω
k
)
2
2Ω
2
)
(S3)
with Ω being a smearing parameter. Note that if Ω

frequency grid spacing ∆
ω
k
, then the
δ
-function is poorly
approximated and the isotropic phonon properties will produce very different macroscopic thermal properties upon
evaluation. If Ω

ω
k
, then there will be a number of zeros in the evaluation of
I
n
(
ω
) resulting in a number of
unrealistic jumps in the isotropic phonon properties. To overcome these problems, we follow the adaptive broadening
scheme for
k
-space integration introduced in ref.
1
to choose the parameter Ω. In this technique, the parameter Ω is
adaptively chosen according to,
n
k
=
a
∂ω
k
|
k
|
=
a
|
v
g
(
n,
k
)
|
|
k
|
(S4)
where
a
is a constant on the order of 1 and
v
g
(
n,
k
) is the group velocity of the phonon mode
n
at a wave vector
k
. Here, Ω is a function of the phonon mode number
n
and the wave vector
k
. Intuitively, when the phonon group
velocity is small, a large number of points gets clustered in the dispersion. The adaptive broadening scheme balances
this effect by including a fewer number of points in the average (equation S2).
B. Smoothness constraint for specularity parameter
To derive a prior distribution function that reflects the smoothness of the specularity parameter, we define smooth-
ness of the specularity profile
p
λ
(where
λ
= 1
,
2
,...
represents the different phonon modes) as,
p
i
=
1
2
(
p
i
1
+
p
i
+1
)
(S5)
To admit some uncertainty in our prior knowledge, we add perturbations to the definition of smoothness as,
p
i
=
1
2
(
p
i
1
+
p
i
+1
) +
W
i
(S6)
where
{
W
i
}∼N
(
0
2
I
)
is a new random variable with variance
γ
. In matrix formulation, we get
Lp
=
W
(S7)
where,
L
=
1
2
1 2
1
.
.
.
.
.
.
0
1 2
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1 2
1 0
.
.
.
.
.
.
1 2
1
(S8)
4
From equation S7, the distribution of
Lp
is the same as the distribution of the random variable
W
. Since
W
is
normally distributed, the prior distribution of the specularity profile
p
λ
is given by,
π
prior
exp
(
1
γ
2
Lp
2
)
(S9)
C. Metropolis Hastings Markov Chain Monte Carlo algorithm
The next step is to sample different specularity profiles from the prior probability distribution (equation S9) and
compute the posterior probability estimate from Bayes theorem, as described in the main manuscript. In this work, we
use a Metropolis Hastings Markov Chain Monte Carlo (MH-MCMC) algorithm to sample trial specularity parameters
p
tr
λ
from the prior probability distribution (equation S9). The MH-MCMC algorithm algorithm proceeds as follows:
Algorithm 1
Metropolis Hastings Markov Chain Monte Carlo algorithm
1:
Choose a proposal density
q
(
m,p
) =
1
2
πγ
2
exp
(
1
2
γ
2
m
p
2
)
.
2:
Choose initial sample
p
0
at random.
3:
for
k
= 0
,...,N
do
4:
Draw a sample
m
from the proposal density
q
(
p
k
,m
)
5:
Compute
π
prior
(
m
)
6:
Compute the acceptance probability
α
(
p
k
,m
) =
min
[
1
,
π
(
m
)
π
(
p
k
)
]
7:
Accept and set
p
k
+1
=
m
with probability
α
(
p
k
,m
). Otherwise, reject and set
p
k
+1
=
p
k
.
In this method, the parameter
γ
represents an artificial
temperature
which controls the level of perturbation on a
profile
p
k
of the
k
th
step to obtain a sample
m
at the (
k
+ 1)
th
step. Since we only deal with ratios of
π
(
p
)’s in this
technique, there is no need to know the proportionality constant in equation S9. Moreover, since the samples
m
are
generated from a proposal distribution
q
(
m,p
) for an independent and identically distributed normal random variable,
standard sampling techniques developed for a multivariate normal distribution can be employed. The information
about the coupling between different components of
p
λ
are included in the expression for the acceptance probability
α
(
p
k
,m
).
Once we have these random trial samples (
p
tr
λ
) drawn according the prior distribution function (equation S9), we
can compute the residual function
N
(
q,T
k
expt
k
BTE
(
p
tr
λ
)
2
I
)
by solving the BTE and finally invoke the
Bayes theorem to obtain the posterior probability density (
π
posterior
(
p
tr
λ
)) for every sample trial specularity profile
p
tr
λ
. While plotting the specularity profiles
p
tr
λ
as in fig. 4(E) in the main manuscript, we use this posterior probability
density (
π
posterior
(
p
tr
λ
)) as the intensity of each of the sampled trial specularity profiles
p
tr
λ
.
D. Convergence of Bayesian inference
To start the Bayesian inference scheme, we choose an initial guess
p
0
that best represents the experimental data
points by trial and error. This process reduces the search space for possible specularity profiles that fit the measure-
ments, thereby reducing computational load. To confirm uniqueness of the final solution, we perform convergence
studies by varying the artificial temperature parameter
γ
and the number of samples, N.
Figure S2 shows the converged specularity profiles for three different values of
γ
. When
γ
is too small, as in fig. S2
(A), the perturbations in the Bayesian sampling are not sufficient to sample the entire specularity parameter solution
space. As
γ
is increased, the sampling procedure spans the entire solution space, resulting in the posterior probability
density smoothly dropping down to 0, as shown in fig. S2 (C) by the low intensity specularity profiles towards the
edges of the gray region. The specularity profiles outside this converged region of solution space cannot explain
the experimental measurements adequately. This convergence behavior demonstrates the uniqueness of the solution