Comparison on two ways of phonon spectral property reconstruction from thermal
spectroscopy experiment

Mojtaba Forghani1111email<EMAIL_ADDRESS>and Nicolas G.
Hadjiconstantinou2

1 Department of Mechanical Engineering, Stanford University, CA
2 Department of Mechanical Engineering, Massachusetts Institute of Technology

###### Abstract

We examine the reliability of the cumulative thermal conductivity as a
function of the free path, $F(\Lambda)$, in the context of reconstruction of
phonon transport properties from thermal spectroscopy experiments. We
specifically show that a given $F(\Lambda)$ does not correspond to a unique
relaxation-times function, in the sense that more than one distribution of
relaxation-times can result in the same $F(\Lambda)$. Since different
relaxation-time distributions will, in general, lead to different thermal
responses, $F(\Lambda)$ does not uniquely predict the material thermal
response in all transport regimes. This implies that in the context of thermal
transport at the nanoscale within the Boltzmann relaxation-time approximation
framework, the ``fundamental" property that should be reconstructed from the
thermal spectroscopy experiments is the frequency-dependent relaxation-times
function (provided group velocities are known), since it explicitly appears in
the governing equation as the input material property. Extensive global
optimization studies show that a previously proposed formulation for
reconstruction based on the frequency-dependent relaxation-times function
[Physical Review B 94, 155439 (2016)] provides numerically unique solutions.

##  1 Introduction

The study of phonon dynamics in the context of nanoscale solid-state heat
transport has received considerable attention [12, 7, 9, 11, 2, 10, 13, 5, 4,
1, 14, 6, 3, 8]. One area of significant interest is the study of phonon
relaxation-time and free path distributions [15, 12, 18, 17, 16, 13, 19, 14].
This information is required for modeling heat transport at the kinetic level,
which becomes necessary due to the failure of Fourier-based analyses at such
small scales. The applications of this study include improved heat management
in nanoelectronic circuits and devices [20, 21, 1, 14, 2, 22],
microelectromechanical sensors [6] and nano-structured materials for improved
thermoelectric conversion efficiency [23, 24, 25, 3, 8, 26].

Purely theoretical approaches such as density functional theory (DFT) have
been widely used to predict phonon properties [5, 4]. However, such approaches
have yet to rise to a stature level where they can replace experimental
results, in part due to their failure to consistently reproduce experimental
observations [4]. Another popular (alternative) class of approaches relies on
extracting material constitutive information from thermal spectroscopy
experiments [12, 7, 9, 11, 10, 13, 8]. However, the analysis of thermal
spectroscopy data remains a challenging task. One typical approach consists of
extracting the cumulative distribution function (CDF) of thermal conductivity
as a function of the free path, $F(\Lambda)$ (see section 2 for more detail),
from the experimentally measured temperature relaxation profiles, by invoking
the concept of ``effective thermal conductivity" and proceeding to match the
experimentally measured response to solutions of the heat conduction equation
with the thermal conductivity (or thermal diffusivity) treated as an
adjustable, ``effective", property [19, 12]. Unfortunately, as it was
discussed extensively before [7], this procedure implicitly assumes that heat
transport is Fourier-like, which is only justified under fairly restrictive
conditions (late times and large scales) that are not always satisfied under
experimental conditions. This is highlighted in figure 1, adapted from [7]
(figure 13 of [7]). The figure shows calculations for the one-dimensional
transient thermal grating (1D-TTG) geometry [7] in Si material for two grating
sizes. In both cases, the reconstruction based on effective thermal
conductivity approach (denoted by ``$k_{\text{eff}}$-based'') predicts an
inaccurate thermal response compared to the one predicted by the methodology
proposed in [7] and further analyzed here (denoted by
``$\tau_{\omega}$-based''). The true thermal response is denoted by
``BTE''(Boltzmann transport equation)—see section 2 for more detail. It is
clear that the error in the $k_{\text{eff}}$-based reconstruction increases as
the length scale decreases.

(a) 100 nm

(b) 400 nm

Figure 1: Thermal response in 1D-TTG experiment for 100 nm (a) and 400 nm (b)
grating sizes. Comparison between the Boltzmann transport equation solution
(denoted ``BTE''), Boltzmann transport equation solution with $\tau_{\omega}$
reconstructed from the BTE solution using the optimization method proposed in
[7, 9] and further analyzed here (denoted by ``$\tau_{\omega}$ based''),
Fourier heat conduction solution with effective thermal conductivity
calculated using the well-known suppression function for a thermal grating
[29, 10] (denoted by ``$\kappa_{\text{eff}}$ based'') and exponential response
obtained from solution of Boltzmann transport equation [7] (denoted by
``Modified $\alpha_{\text{eff}}$'').

Figure 2 shows examples of reconstructing Si properties in a 2D-dots geometry
[19] from synthetic (numerical) data, adapted from [9] (figure 10 of [9]). The
two $\kappa_{\text{eff}}(L)/\kappa$ plots (figure 2a) correspond to the
reconstructed cumulative effective heat conductivity functions, as a function
of the system length scale, based on two different Al-Si interface models
(transmissivity functions) that have been used in the numerical experiment
(Monte-Carlo solution of the BTE); neither is able to recover the material
thermal conductivity (reach the value of 1). This is due to experiment not
satisfying the late time and large scale criteria needed for the
$\kappa_{\text{eff}}$-based approach to be valid (see [9] for more detail).
Consequently, the two reconstructed $F(\Lambda)$ shown in figure 2b,
corresponding to reconstructions based on these effective heat conductivities,
have not been able to recover the true CDF.

(a)

(b)

Figure 2: Reconstruction results using the effective thermal conductivity
approach for the 2D-dots geometry, based on synthetic data generated by
solution of the BTE [9]. The effective thermal conductivities do not match if
different interface models (DMM denotes diffuse mismatch model; expt. denotes
experimentally determined transmissivities [30]) are used when generating the
synthetic data (a). As a result, the corresponding $F(\Lambda)$ do not match
and are each different from the true $F(\Lambda)$ (b).

As noted above, in response to the failure of $\kappa_{\text{eff}}$-based
methods to provide accurate reconstructions at all scales, in [7, 9, 11] we
developed and extensively validated an optimization-based methodology for
directly reconstructing phonon relaxation-times from thermal spectroscopy
data. In this work, we attempt to address some more fundamental questions,
such as the uniqueness of reconstructed quantities (e.g., $F(\Lambda)$), as
well as their ability to reproduce the experimental results. In particular, we
show that reconstructions that use $F(\Lambda)$ (as a means to defining
``effective thermal conductivities") fail to guarantee a unique thermal
response. This can be understood by noting that it is the relaxation-time and
not $F(\Lambda)$ which appears as a primary material-parameter input to the
BTE. Specifically, we show that the same $F(\Lambda)$ can be arrived at from
multiple relaxation-time distributions, which, in general, will correspond to
multiple thermal responses. This observation seriously questions the
reliability of approaches which use $F(\Lambda)$ in the reconstruction.

We also study the reliability of the solutions returned by the optimization
algorithm proposed in [7] and in particular the solution uniqueness and
sensitivity to noisy measurements. Our results show that the optimized
solution exhibits properties associated with a unique minimum; furthermore,
the solution is sufficiently robust to noisy measurements.

The remainder of the paper is organized as follows. In section 2, we briefly
overview the governing equations, the linearized BTE and related conservation
laws. In section 3, we show theoretically and numerically that more than one
relaxation-times distribution, and consequently material thermal responses,
can be related to the same $F(\Lambda)$, and hence, the latter may not be used
to study thermal behaviors in all regimes. In section 4, we investigate the
uniqueness of the solution associated with our previously proposed algorithm
which instead reconstructs the relaxation-times function; we also investigate
the sensitivity of this approach to the uncertainty in the temperature
measurements from a Bayesian perspective. Finally, in section 5 we provide a
summary of our work and suggestions for future improvements.

##  2 Boltzmann transport equation

Given the small temperature differences usually associated with thermal
spectroscopy experiments, here we consider the linearized BTE

$\frac{\partial e^{d}}{\partial
t}+\textbf{v}_{\omega}\cdot\nabla_{\textbf{x}}e^{d}=-\frac{e^{d}-(de^{\text{eq}}/dT)_{T_{\text{eq}}}\Delta\widetilde{T}}{\tau_{\omega}},$
(1)

where
$e^{d}=e^{d}(t,\mathbf{x},\omega,\mathbf{\Omega})=e-e^{\text{eq}}_{T_{\text{eq}}}=\hbar\omega(f-f^{\text{eq}}_{T_{\text{eq}}})$
is the deviational energy distribution, $\omega$ is the phonon frequency,
$\mathbf{\Omega}$ is the phonon traveling direction,
$f=f(t,\mathbf{x},\omega,\mathbf{\Omega})$ is the occupation number of phonon
modes, $\textbf{v}_{\omega}=\textbf{v}(\omega)$ is the phonon group velocity,
$\tau_{\omega}=\tau(\omega)$ is the frequency-dependent relaxation-time, also
referred to here as the ``relaxation-times function", and $\hbar$ is the
reduced Planck constant. Here and in what follows, we use $\omega$ to denote
the dependence on both frequency and polarization.

The above equation is linearized about the equilibrium temperature
$T_{\text{eq}}$, to be understood here as the experimental reference
temperature. In general, $\tau_{\omega}=\tau(\omega,T)$; however, as a result
of the linearization,
$\tau_{\omega}=\tau(\omega,T_{\text{eq}})\equiv\tau(\omega)$; in other words,
the solutions (and associated reconstruction) are valid for the experiment
baseline temperature $T_{\text{eq}}$. Also,
$(de^{\text{eq}}/dT)_{T_{\text{eq}}}=\hbar\omega(df^{\text{eq}}_{T}/dT)|_{T_{\text{eq}}}$
and $f^{\text{eq}}_{T}$ is the Bose-Einstein distribution with temperature
parameter $T$, given by

$f^{\text{eq}}_{T}(\omega)=\frac{1}{\exp(\hbar\omega/k_{B}T)-1},$ (2)

where $k_{B}$ is Boltzmann's constant. Finally,
$\Delta\widetilde{T}=\Delta\widetilde{T}(t,\mathbf{x})=\widetilde{T}-T_{\text{eq}}$
is referred to as the deviational pseudo-temperature
($\widetilde{T}(t,\mathbf{x})$ is the pseudo-temperature). Note that the
deviational pseudo-temperature, which is different from the deviational
temperature defined below, is defined by the energy conservation statement
[27]

$\int_{\mathbf{\Omega}}\int_{\omega}\left[\frac{C_{\omega}}{\tau_{\omega}}\Delta\widetilde{T}-\frac{e^{d}}{\tau_{\omega}}D_{\omega}\right]d\omega
d\mathbf{\Omega}=0,$ (3)

in which $D_{\omega}=D(\omega)$ is the density of states,
$C_{\omega}=C(\omega;T_{\text{eq}})=D_{\omega}(de^{\text{eq}}/dT)_{T_{\text{eq}}}$
is the frequency-dependent volumetric heat capacity, and
$d\mathbf{\Omega}=\sin(\theta)d\theta d\phi$ represents the differential solid
angle element such that $\theta$ and $\phi$ refer to the polar and azimuthal
angles in the spherical coordinate system, respectively. The temperature
$T(t,\mathbf{x})$ is computed from

$\int_{\mathbf{\Omega}}\int_{\omega}\left[C_{\omega}\Delta
T-e^{d}D_{\omega}\right]d\omega d\mathbf{\Omega}=0,$ (4)

where $\Delta T(t,\mathbf{x})=T-T_{\text{eq}}$ is the deviational temperature.
The frequency-dependent free path is given by

$\Lambda_{\omega}=v_{\omega}\tau_{\omega},$ (5)

where $v_{\omega}=||\mathbf{v}_{\omega}||$ is the group velocity magnitude.
Cumulative distribution function (CDF) of thermal conductivity as a function
of the free path, introduced in section 1, is defined as
$F(\Lambda):=\frac{1}{3\kappa}\int_{\omega^{*}(\Lambda)}C_{\omega}v_{\omega}^{2}\tau_{\omega}d\omega$,
where $\omega^{*}(\Lambda)$ is the set of modes such that
$\omega^{*}(\Lambda)=\\{\omega|\Lambda_{\omega}\leq\Lambda\\}$; similarly, the
corresponding probability density function (of thermal conductivity) is given
by $\mathfrak{f}=\frac{dF}{d\Lambda}$.

##  3 Limitations associated with $F(\Lambda)$

In this section, we investigate, both numerically and analytically, whether
$F(\Lambda)$ can be used to describe heat transfer in all transport regimes.
This investigation is motivated by the widespread use of $F(\Lambda)$ in
procedures for reconstructing material properties from thermal spectroscopy
experiments.

The approximations and theoretical inconsistencies associated with the use of
the concept of effective thermal conductivity, in conjunction with
$F(\Lambda)$, to reconstruct thermal spectroscopy data have been discussed
before in [7, 9], as well as section 1. Here, we take a different approach
which allows us to examine the premise of such approaches on a more general
level. First, we prove through a theorem that a given $F(\Lambda)$ may arise
as a result of more than one distribution of relaxation-times. Since different
relaxation-time distributions will, in general, lead to different material
thermal behaviors, this implies that $F(\Lambda)$ cannot uniquely determine
the material thermal behavior, which in turn questions the ability of
$F(\Lambda)$ to describe thermal responses in sub-micron regimes. The above
assertion is also demonstrated through a numerical experiment using a
realistic material model for silicon.

#### Theorem:

For a given material, if there are two phonon frequencies, $\omega_{1}$ and
$\omega_{2}$, such that
$C_{\omega}({\omega_{1}})v_{\omega}({\omega_{1}})=C_{\omega}({\omega_{2}})v_{\omega}({\omega_{2}})$,
there exist more than one relaxation-time distributions that lead to a given
$F(\Lambda)$.

The proof of theorem is left to Appendix A.

###  3.1 Numerical Demonstration

In order to highlight the importance of the above in practical terms, namely,
the existence of multiple relaxation-time functions for a given $F(\Lambda)$,
and more importantly, the non-uniqueness of the thermal response for a given
$F(\Lambda)$, we have performed a calculation that directly studies this
hypothesis. In the interest of simplicity, all calculations were performed for
the 1D-TTG geometry with the Holland model for material Si [7, 9]. Readers are
referred to Appendix B for a discussion of the functional form of the product
$C_{\omega}({\omega})v_{\omega}({\omega})$ for Si, used in our calculations.
We consider the objective function

$\min_{{\bf
U}}\int_{\Lambda}\big{|}\mathfrak{f}^{*}(\Lambda)-\mathfrak{f}(\Lambda)\big{|}d\Lambda,$
(6)

where $\mathfrak{f}^{*}$ represents the ``true" free path-dependent
distribution function of thermal conductivity and ${\bf U}$ is a set of
parameters that parameterizes the relaxation-times function. Here, without
loss of generality, we have assumed ${\bf U}$ to be a set of parameters that
allows different branches of the relaxation-times function to take the form of
a piece-wise linear function with three segments as a function of the phonon
frequency; this is the same parameterization that has been used previously in
[7, 9, 11]. This objective function attempts to find free path-dependent
thermal conductivity distribution functions $\mathfrak{f}(\Lambda)$ that match
the true thermal conductivity distribution function by searching in the
parameter space ${\bf U}$. Note that ${\bf U}$ enters the objective function
through the definition of the free path distribution
$\Lambda_{\omega}=v_{\omega}\tau_{\omega}({\bf U})$.

We perform the optimization using the Nelder-Mead (NM) algorithm [28],
repeating it multiple times starting from different initial conditions for
${\bf U}$ in order to also study the effect of different initial conditions on
the optimization process. Figure 3 shows the relaxation-times and the
$F(\Lambda)$ obtained from one of these optimization processes as well as the
temperature profile predicted by the optimized ${\bf U}$ versus the true
profile at 100 nm length scale. This profile is generated using the inverse
fast Fourier transform (IFFT) [33] as the forward simulation method [7]. As it
is expected from equation (6), the reconstructed $F(\Lambda)$ matches the true
CDF. At the same time, both the relaxation-time functions, and more
importantly, the thermal responses do not match, even though all other
material properties, as well as the geometry are kept fixed. The mismatch in
thermal responses also exists at other length scales.

(a) Relaxation-time

(b) $F(\Lambda)$

(c) Temperature

Figure 3: Comparison between properties predicted by one representative
solution of optimization (6) and true data. As expected, although the two
$F(\Lambda)$ match, the relaxation-times and thermal responses are clearly
different.

The results provided in figure 3 clearly show, from a numerical perspective,
that $F(\Lambda)$ cannot be used reliably as a material property for
predicting thermal behaviors in all regimes. While the non-uniqueness of the
relaxation-times function for a given $F(\Lambda)$, as was proven
theoretically in this section and shown numerically in figure 3(a), may not be
catastrophic on its own, the fact that these different functions predict
different thermal responses, as we have seen in figure 3(c), implies that a
given $F(\Lambda)$ provides insufficient information for predicting the
thermal response. The results provided in this section prove that $F(\Lambda)$
cannot be used to predict thermal responses in all transport regimes.

##  4 An alternative to $F(\Lambda)$

The discussions provided previously in [7, 9, 11], and more importantly, in
the previous section, show that $F(\Lambda)$, regardless of the method being
used to reconstruct it and its accuracy, is not able to predict the material
thermal behavior in all regimes, and the sub-micron regime in particular,
which is of interest presently. This is in sharp contrast with the relaxation-
times function which is related to the system thermal response through a
direct, one-to-one relationship.

In previous work we proposed a multi-stage NM algorithm for reconstruction of
relaxation-times function from thermal spectroscopy experiments [7, 9, 11] and
validated it both numerically and experimentally. In this section, we briefly
study the reliability of the solutions of this previously proposed
optimization algorithm, that is, the uniqueness of the optimization solution
as well as a Bayesian-approach inspired study of the sensitivity of the
solutions to measurement noise. Due to the complexity of the BTE and the
related inverse problems, here we have focused on the numerical study of a
specific problem, the 1D-TTG experiment.

Due to the success of the NM algorithm in providing solutions to the
reconstruction problem [7, 9, 11], here, we use a globalized version of NM
algorithm, referred to as the ``globalized bounded Nelder-Mead (GBNM)" [31],
to study the uniqueness of the solutions to the reconstruction problem
formulation used in [7, 9, 11]. The GBNM algorithm is very similar to the NM
algorithm described before; however, instead of performing the optimization
once, starting from one initial condition, it repeats the complete
optimization process many times, starting from different probabilistically
correlated initial conditions. More details can be found in [31].

###  4.1 On the uniqueness of reconstructed solutions

To assess the uniqueness of the reconstructed solutions obtained via the
optimization formulation, we follow a procedure that is similar to our
previous work [7]: we generate synthetic temperature profiles for a material
by solving the BTE, and use the generated data to infer the material
relaxation-time by solving the following optimization problem

$\min_{{\bf U}}\mathcal{L}=\min_{{\bf
U}}\left[\frac{\sum_{t,\textbf{x},L}|T_{\text{m}}(t,\mathbf{x};L)-T_{\text{BTE}}(t,\mathbf{x};L,{\bf
U})|}{N}+\alpha\Bigg{|}1-\frac{1}{3\kappa}\int_{\omega}C_{\omega}\tau_{\omega}({\bf
U})v^{2}_{\omega}d\omega\Bigg{|}\right],$ (7)

where $T_{\text{m}}(t,\mathbf{x};L)$ denotes the experimentally measured
temperature (equivalent to synthetically generated temperatures),
$T_{\text{BTE}}$ is the temperature computed from solution of the BTE (the
same temperature as in equation (4)), $N$ is the total number of (independent)
measurements available, and $L$ is the different characteristic length scales.
In the present case, we have assumed the material to be Si and considered its
response in the 1D-TTG geometry for 10 nm$<L<$100 $\mu$m, solved for using the
IFFT method [33, 7]. We also assume that the two $TA$ branches of the material
are the same, using the same piece-wise linear parameterization of [7], for
both the ab initio and the Holland models [7, 9, 11].

As stated above, the difference from previous work is that, here, the
optimization is performed using the GBNM method. Figure 4 shows the value of
the objective function for the various points that the GBNM method explored on
its way to minimizing ${\cal L}$. More detail on the procedure and parameters
used here can be found in [32]. The horizontal axis in the figure corresponds
to the distance between a ``trial" relaxation-times function (whose value of
${\cal L}$ is plotted on the ordinate) and the ``true" solution
$\tau^{*}_{\omega}$; specifically, ``distance" is defined by

$\text{distance}:=\frac{{\bigintss}_{\omega}\left|\log\left(\frac{\tau_{\omega}}{\tau^{*}_{\omega}}\right)\right|D_{\omega}d\omega}{\int_{\omega}D_{\omega}d\omega}.$
(8)

Here, $\omega$, as before, contains information of both different frequencies
and branches.

(a)

(b)

Figure 4: The value of the objective function versus the distance (see
equation (8)) between the points sampled during GBNM and the true solution for
the Holland model (a) and the ab initio model (b).

The main message from figure 4 is that the bottom left of the cluster of
datapoints in both plots represents a monotonically increasing function, which
means that the closer we are to the true solution, the smaller the value of
the objective function becomes, and more importantly, this is a monotonic
behavior as we approach the true solution. This implies that the solution is
unique; non-uniqueness would be implied by points in the lower right hand
corner of the diagrams (low value of objective function, while the distance
from the solution is large). Here we note that the number of datapoints in
each plot in figure 4 is approximately 1 million. While exploring the whole
parameter space with high resolution is not feasible due to the high dimension
of the problem studied here ($\text{dim}({\bf U})=12$; two branches, each
parameterized with three lines), the general behavior of the objective
function, such as the monotonicity described above, will likely not change if
more iterations are used (and consequently an even better exploration of the
parameter space is performed).

###  4.2 Sensitivity of reconstructed relaxation-time to measurement
uncertainty

While the results provided in figure 4 show that the behavior of the objective
function $\cal{L}$ introduced in [7] resembles that of a function with a
unique minimum, the sensitivity of this solution to noise in the measurement
is not clear. Low sensitivity of an objective function to noise plays an
important role in its practical application to real-world measurements. Here
we study this sensitivity through the Bayesian theorem, that is, using the
relation

$p\Big{(}{\bf U}\big{|}\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L}\Big{)}\propto
p\Big{(}\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L}\big{|}{\bf U}\Big{)}\pi({\bf
U}),$ (9)

where $p\Big{(}{\bf U}\big{|}\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L}\Big{)}$
is the distribution of the parameters of the relaxation-times function given
the temperature measurements (the posterior distribution) — the quantity that
we are interested in— from which we can also calculate the distribution of the
relaxation-times function itself through the $\tau_{\omega}({\bf U})$
function; $p\Big{(}\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L}\big{|}{\bf
U}\Big{)}$ is the distribution of measured temperatures, given the parameters
of the relaxation-times function (the likelihood function). The likelihood
function can be calculated by first solving the BTE for a given relaxation-
times function $\tau_{\omega}({\bf U})$, and then adding an artificial noise
to its predicted temperatures that resembles that of noisy measurements.
Finally, $\pi({\bf U})$ is the prior distribution which is usually a wide
(least informative) distribution of the quantity of interest, the relaxation-
times function; for instance, a wide symmetrical distribution around the
correct relaxation-times function. More detailed definitions and information
on equation (9) are provided in Appendix C.

Once we infer the distribution of the different components of U (the posterior
distribution), we can calculate the distribution of the relaxation-times
function as a function of its frequency $\tau_{\omega}({\bf U})$. The results
for the ab initio silicon model are provided in figure 5. We have also plotted
the true relaxation-time functions in this figure. We notice that the
distribution matches the true solution ``in the mean''. We also observe that
for $LA$ modes, at frequencies around $\omega=6\times 10^{13}$ rad/s, the
distribution is sharper than other frequencies, implying that a relatively
more accurate reconstruction of the value of relaxation-time at those
frequencies is possible. On the other hand, the uncertainty at very large
frequencies ($\omega\gtrsim 6.5\times 10^{13}$ rad/s) is more than the other
ranges of the frequency. A similar trend can be observed for $TA$ modes.
Overall, by comparing the two plots in figure 5, we conclude that the
reconstruction of $TA$ branches is more accurate (the distributions are
sharper). This could be due to the assumption that the two $TA$ branches have
been approximated by the same function, which consequently led to their
greater influence on the thermal behavior.

(a) $LA$ modes

(b) $TA$ modes

Figure 5: Contour plot of the probability density of the frequency-dependent
relaxation-times distribution of $LA$ and $TA$ modes. The red line in (a)
denotes the true $\tau_{\omega}^{LA}$, while the red line in (b) denotes the
true $\tau_{\omega}^{TA_{1}}$ and the green line in (b) denotes the true
$\tau_{\omega}^{TA_{2}}$.

In order to obtain a better picture of the shape of the distributions at
different frequencies, we have also plotted the distributions at a few
specific frequencies (equivalent to cross-sections of figure 5 at different
frequencies). Figures 6 and 7 show the prior distribution,
$\pi(\tau_{\omega})$, versus the posterior distribution,
$p(\tau_{\omega}\big{|}\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L})$, for a few
frequencies of $LA$ and $TA$ modes, respectively. We observe that in most
cases, the posterior distribution has become significantly narrower compared
to the prior distribution, implying that the chosen prior distribution has
been wide enough to avoid bias in inferring the posterior distribution. The
distributions in the $TA$ case are sharper, as can also be seen in figure 5.
The distribution at very high frequencies ($\omega=7$ rad/s in figure 6 and
$\omega=3.5$ rad/s in figure 7) is wider than the other frequencies,
consistent with our previous observations from figure 5. While the results
provided in figure 4 point to the uniqueness of the objective function minimum
point, the narrow distributions observed in figures 5–7 complement those
results by implying that good accuracy of solution in the vicinity of the
optimal solution, minimally affected by the measurement noise, is attainable.

(a)

(b)

Figure 6: Distributions at a few low frequency (a) and high frequency (b) $LA$
modes.

(a)

(b)

Figure 7: Distributions at a few low frequency (a) and high frequency (b) $TA$
modes.

##  5 Summary and outlook

We have studied the type and reliability of information that can be extracted
from thermal spectroscopy experiments. We specifically showed that inverting
the experimental data using $F(\Lambda)$ (presumably via the effective thermal
conductivity construct) leads to an ill-posed problem because any given
$F(\Lambda)$ does not correspond to a unique system thermal response; as a
result, even correct determination of $F(\Lambda)$ from the experiment is
insufficient by itself to fully characterize the system behavior in other
situations. Although the condition
$C_{\omega}({\omega_{1}})v_{\omega}({\omega_{1}})=C_{\omega}({\omega_{2}})v_{\omega}({\omega_{2}})$
in section 3 is more general, in practice, this condition is typically met due
to the presence of more than one phonon branch, which $F(\Lambda)$ is unable
to distinguish between. Our previous results in [7, 9, 11, 32]—also briefly
reviewed in section 1— additionally showed that the Fourier heat conduction
equation with effective thermal conductivity (or diffusivity) can be used for
the reconstruction of $F(\Lambda)$ only if certain late time and large scale
assumptions are met (see [7, 9] and figures 1 and 2 of the present
manuscript), and is thus not applicable at arbitrarily small length scales.

Our results in section 4 are consistent with our previous work, showing that
the algorithm proposed in [7] can be reliably used for future reconstruction
purposes on real experimental data (as it has already been performed in [11])
which may include noise in their measurements, with no concerns about the non-
uniqueness of the solution or low sensitivity of the relaxation-times function
to the measured thermal responses. Note that although finding this unique
solution requires a global optimization algorithm which is significantly more
expensive, possibly by orders of magnitude, than the multi-stage NM
optimization algorithm used in [7, 9, 11], this is not required for regular
applications. In applications, the multi-stage algorithm proposed in [7]
should provide a solution that is sufficiently close to the global optimum at
a fraction of the cost. Specifically, 1 million forward simulations (the
number of datapoints shown in figure 4) were used to find the approximate
location of the global minimum using GBNM, while as it was discussed in [7],
our multi-stage algorithm requires at most about 1000 forward simulations to
obtain this solution.

While in the present work the numerical analysis of the reliability and
uniqueness of our previously proposed algorithm [7] has been focused on the
1D-TTG experiment, it is important to note that many thermal spectroscopy
experiments also include solid-solid interfaces in their setup, which can add
to the complexity of the analysis, for instance, due to lack of sufficient
information regarding the interface conductance and transmissivities. However,
in a previous work [9], we have studied the reliability of our proposed
reconstruction process on the 2D-dots geometry [19] and have observed that not
knowing the interface transmissivities does not influence the reconstruction
process significantly, and, in fact, the phonon relaxation-times function can
be reconstructed accurately without knowing any information of the
transmissivity functions. In future, we will also perform a study similar to
present work that aims at uniqueness and sensitivity of the reconstructed
relaxation-times in the presence of solid-solid interface with unknown
properties.

The processes being described in section 4 introduce new complementary tools
for analyses of thermal spectroscopy measurements and open new avenues for
future similar studies by combining them with the available theoretical tools,
leading to a unified numerical-analytical approach for categorization of
different problems based on the uniqueness of their inverse problems, in the
context of phonon transport. This will allow researchers to explore the
uniqueness of the solution for a particular geometry before performing the
experiment, and thus helping them to focus on experiments whose solution to
their inverse BTE equations is well-posed. Similarly, the Bayesian approach
can help with visualization of sensitivity of the solution around the true
solution, thus guiding the researchers towards geometries that reveal high
sensitivity in the vicinity of the solution, and consequently, are more robust
to measurement noise.

## Acknowledgment

The authors would like to thank Y. M. Marzouk and G. Chen for their comments
and suggestions. This work was supported by the Solid-State Solar-Thermal
Energy Conversion Center (S3TEC), an Energy Frontier Research Center funded by
the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under
Award# DE-SC0001299 and DE-FG02-09ER46577.

## References

  * [1] E. Pop, S. Sinha, and K. E. Goodson, Heat generation and transport in nanometer-scale transistors, proceedings of the IEEE 94, 1587-1601 (2006).
  * [2] D. G. Cahill, K. Goodson, and A. Majumdar, Thermometry and thermal transport in micro/nanoscale solid-state devices and structures, Journal of Heat Transfer 124, 223-241 (2001).
  * [3] Z. Tian, S. Lee, G. Chen, Heat transfer in thermoelectric materials and devices, Journal of Heat Transfer 135, 061605 (2013).
  * [4] K. Esfarjani, G. Chen, and H. T. Stokes, Heat transport in silicon from first-principles calculations, Physical Review B 84, 085204 (2011).
  * [5] A. Ward and D. A. Broido, Intrinsic phonon relaxation times from first-principles studies of the thermal conductivities of Si and Ge, Physical Review B 81, 085205 (2010).
  * [6] D. G. Cahill, P. V. Braun, G. Chen, D. R. Clarke, S. Fan, K. E. Goodson, P. Keblinski, W. P. King, G. D. Mahan, A. Majumdar, H. J. Maris, S. R. Phillpot, E. Pop, and L. Shi, Nanoscale thermal transport. II. 2003-2012, Applied Physics Reviews 1, 011305 (2014).
  * [7] M. Forghani, N. G. Hadjiconstantinou, and J.-P. M. Péraud, Reconstruction of the phonon relaxation times using solutions of the Boltzmann transport equation, Physical Review B 94, 155439 (2016).
  * [8] M. Zebarjadi, K. Esfarjani, M. S. Dresselhaus, Z. F. Ren, and G. Chen, Perspectives on thermoelectrics: from fundamentals to device applications, Energy and Environmental Science 5, 5147-5162 (2012).
  * [9] M. Forghani and N. G. Hadjiconstantinou, Reconstruction of phonon relaxation times from systems featuring interfaces with unknown properties, Physical Review B 97, 195440 (2018).
  * [10] J. A. Johnson, A. A. Maznev, J. Cuffe, J. K. Eliason, A. J. Minnich, T. Kehoe, C. M. Sotomayor Torres, G. Chen, and K. A. Nelson, Direct measurement of room-temperature nondiffusive thermal transport over micron distances in a silicon membrane, Physical Review Letters 110, 025901 (2013).
  * [11] M. Forghani and N. G. Hadjiconstantinou, Phonon relaxation time reconstruction from transient thermal grating experiments and comparison with density functional theory predictions, Applied Physics Letters 114, 026103 (2019).
  * [12] A. J. Minnich, J. A. Johnson, A. J. Schmidt, K. Esfarjani, M. S. Dresselhaus, K. A. Nelson, and G. Chen, Thermal conductivity spectroscopy technique to measure phonon mean free paths, Physical Review Letters 107, 095901 (2011).
  * [13] A. J. Minnich, Determining phonon mean free paths from observations of quasiballistic thermal transport, Physical Review Letters 109, 205901 (2012).
  * [14] E. Pop, Energy dissipation and transport in nanoscale devices, Nano Research 3, 147-169 (2010).
  * [15] D. P. Sellan, J. E. Turney, A. J. H. McGaughey, and C. H. Amon, Cross-plane phonon transport in thin films, Journal of Applied Physics 108, 113524 (2010).
  * [16] F. Yang and C. Dames, Heating-frequency-dependent thermal conductivity: An analytical solution from diffusive to ballistic regime and its relevance to phonon scattering measurements, Physical Review B 91, 165311 (2015).
  * [17] K. T. Regner, A. J. H. McGaughey, and J. A. Malen, Analytical interpretation of nondiffusive phonon transport in thermoreflectance thermal conductivity measurements, Physical Review B 90, 064302 (2014).
  * [18] A. A. Maznev, J. A. Johnson, and K. A. Nelson, Onset of nondiffusive phonon transport in transient thermal grating decay, Physical Review B 84, 195206 (2011).
  * [19] Y. Hu, L. Zeng, A. J. Minnich, M. S. Dresselhaus, and G. Chen, Spectral mapping of thermal conductivity through nanoscale ballistic transport, Nature Nanotechnology 10, 701-706 (2015).
  * [20] J.-K. Yu, S. Mitrovic, D. Tham, J. Varghese, and J. R. Heath, Reduction of thermal conductivity in phononic nanomesh structures, Nature Nanotechnology 5, 718-721 (2010).
  * [21] M. C. Wingert, Z. C. Y. Chen, E. Dechaumphai, J. Moon, J.-H. Kim, J. Xiang, and R. Chen, Thermal conductivity of Ge and GeSi coreshell nanowires in the phonon confinement regime, Nano Letters 11, 5507-5513 (2011).
  * [22] S.-M. Lee and D. G. Cahill, Heat transport in thin dielectric films, Journal of Applied Physics 81, 2590 (1997).
  * [23] K. Biswas, J. He, I. D. Blum, C.-I Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid, and M. G. Kanatzidis, High-performance bulk thermoelectrics with all-scale hierarchical architectures, Nature 489, 414-418 (2012).
  * [24] A. I. Hochbaum, R. Chen, R. D. Delgado, W. Liang, E. C. Garnett, M. Najarian, A. Majumdar, and P. Yang, Enhanced thermoelectric performance of rough silicon nanowires, Nature 451, 163-167 (2008).
  * [25] A. I. Boukai, Y. Bunimovich, J. Tahir-Kheli, J.-K. Yu, W. A. Goddard III, and J. R. Heath, Silicon nanowires as efficient thermoelectric materials, Nature 451, 168-171 (2008).
  * [26] D. Kraemer, B. Poudel, H.-P. Feng, J. C. Caylor, B. Yu, X. Yan, Y. Ma, X. Wang, D. Wang, A. Muto, K. McEnaney, M. Chiesa, Z. Ren, and G. Chen, High-performance flat-panel solar thermoelectric generators with high thermal concentration, Nature Materials 10, 532-538 (2011).
  * [27] J.-P. M. Péraud, C. D. Landon, and N. G. Hadjiconstantinou, Monte Carlo methods for solving the Boltzmann transport equation, Annual Review of Heat Transfer 17, 205-265 (2014).
  * [28] J. A. Nelder and R. Mead, A simplex method for function minimization, The Computer Journal 7, 308-313, (1965).
  * [29] C. Hua and A. J. Minnich, Transport regimes in quasiballistic heat conduction, Physical Review B 89, 094302 (2014).
  * [30] C. Hua, X. Chen, N. K. Ravichandran, and A. J. Minnich, Experimental metrology to obtain thermal phonon transmission coefficients at solid interfaces, Physical Review B 95, 205423 (2017).
  * [31] M. A. Luersen and R. L. Riche, Globalized Nelder-Mead method for engineering optimization, Computers and Structures 82, 2251-2260 (2004).
  * [32] M. Forghani, An inverse problem framework for reconstruction of phonon properties using solutions of the Boltzmann transport equation, PhD thesis, Massachusetts Institute of Technology, Cambridge, MA (2019).
  * [33] E. O. Brigham, The fast Fourier transform, Prentice-Hall, New Jersey (1988).
  * [34] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, and A. H. Teller, Equation of state calculations by fast computing machines, Journal of Chemical Physics 21, 1087 (1953).
  * [35] G. O. Roberts, A. Gelman, and W. R. Gilks, Weak convergence and optimal scaling of random walk Metropolis algorithms, The Annals of Applied Probability 7, 110-120 (1997).

## Appendix A: Proof of theorem

Proof: We consider two different relaxation-time functions, $\tau_{\omega}$
and $\tau^{\prime}_{\omega}$, with free path-dependent distribution function
of thermal conductivity $\mathfrak{f}(\Lambda)$ and
$\mathfrak{f}^{\prime}(\Lambda)$, respectively; $\tau_{\omega}$ and
$\tau^{\prime}_{\omega}$ are identical
($\tau_{\omega}=\tau^{\prime}_{\omega}$) for all
$\omega\neq\omega_{1},\omega_{2}$, while at $\omega_{1}$, $\omega_{2}$ they
satisfy the relations

$\tau^{\prime}_{\omega}(\omega_{1})=\frac{\tau_{\omega}(\omega_{2})v_{\omega}(\omega_{2})}{v_{\omega}(\omega_{1})}\
\text{,}\
\tau^{\prime}_{\omega}(\omega_{2})=\frac{\tau_{\omega}(\omega_{1})v_{\omega}(\omega_{1})}{v_{\omega}(\omega_{2})},$
(10)

where we assume $v_{\omega}(\omega_{2})\neq v_{\omega}(\omega_{1})$.

From the relaxation-times function $\tau_{\omega}$ one can derive the free
path $\Lambda_{\omega}$ and frequency-dependent heat conductivity
$\kappa_{\omega}$
($\kappa_{\omega}:=C_{\omega}v_{\omega}^{2}\tau_{\omega}/3$), considered here
as a pair
$(\Lambda_{\omega},\kappa_{\omega})=(v_{\omega}\tau_{\omega},C_{\omega}v_{\omega}^{2}\tau_{\omega}/3)$.
Since the two relaxation-times functions are identical at all
$\omega\neq\omega_{1},\omega_{2}$, their corresponding
$(\Lambda_{\omega},\kappa_{\omega})$ pairs are also identical at all
$\omega\neq\omega_{1},\omega_{2}$. To show that their pairs are identical for
all frequencies, and thus they have the same CDF of thermal conductivity, we
now show that

$(\Lambda_{\omega}(\omega_{1}),\kappa_{\omega}(\omega_{1}))=(\Lambda_{\omega}^{\prime}(\omega_{2}),\kappa^{\prime}_{\omega}(\omega_{2}))$
(11)

and

$(\Lambda_{\omega}(\omega_{2}),\kappa_{\omega}(\omega_{2}))=(\Lambda_{\omega}^{\prime}(\omega_{1}),\kappa^{\prime}_{\omega}(\omega_{1}))$
(12)

Here, we have used the definition
$(\Lambda_{\omega}(\omega_{i}),\kappa_{\omega}(\omega_{i})):=(v_{\omega}(\omega_{i})\tau_{\omega}(\omega_{i}),C_{\omega}(\omega_{i})v_{\omega}^{2}(\omega_{i})\tau_{\omega}(\omega_{i})/3)$,
for $i=1,2$ (and similarly for the
$(\Lambda^{\prime}_{\omega}(\omega_{i}),\kappa^{\prime}_{\omega}(\omega_{i}))$
pair).

The above can be shown by writing

$\tau^{\prime}_{\omega}(\omega_{1})v_{\omega}(\omega_{1})=\left[\frac{\tau_{\omega}(\omega_{2})v_{\omega}(\omega_{2})}{v_{\omega}(\omega_{1})}\right]v_{\omega}(\omega_{1})=\tau_{\omega}(\omega_{2})v_{\omega}(\omega_{2}),$
(13)

where we have used equation (10) to obtain the second relationship. The heat
conductivity value corresponding to this free path is

$\displaystyle\frac{1}{3}C_{\omega}(\omega_{1})v_{\omega}^{2}(\omega_{1})\tau^{\prime}_{\omega}(\omega_{1})$
$\displaystyle=\frac{1}{3}C_{\omega}(\omega_{1})v_{\omega}(\omega_{1})\left[v_{\omega}(\omega_{1})\tau^{\prime}_{\omega}(\omega_{1})\right]$
$\displaystyle=\frac{1}{3}C_{\omega}(\omega_{2})v_{\omega}(\omega_{2})\left[v_{\omega}(\omega_{2})\tau_{\omega}(\omega_{2})\right]$
$\displaystyle=\frac{1}{3}C_{\omega}(\omega_{2})v_{\omega}^{2}(\omega_{2})\tau_{\omega}(\omega_{2}),$
(14)

where we have used both the assumption $v_{\omega}(\omega_{2})\neq
v_{\omega}(\omega_{1})$ and equation (13) to obtain the second equality.
Relations (13) and (14) are equivalent to relation (12).

Using the same process as (13) and (14) for $\omega_{2}$ (replacing
$\omega_{1}$ by $\omega_{2}$), we obtain relation (11).

These results (relations (11) and (12)) imply that
$\mathfrak{f}(\Lambda)=\mathfrak{f}^{\prime}(\Lambda)$ for all free paths
while there are frequencies for which
$\tau_{\omega}(\omega)\neq\tau^{\prime}_{\omega}(\omega)$ $\blacksquare$.

## Appendix B: Silicon material properties

Figure 8 shows the plot of $C_{\omega}v_{\omega}$ for two different sets of
DFT-based silicon material properties, referred to as the Holland and the ab
initio models [7, 9, 11, 32]. Longitudinal acoustic and transverse acoustic
branches are denoted by $LA$ and $TA$, respectively. We observe that in both
cases, considering all branches, a number of modes have the same
$C_{\omega}v_{\omega}$ value for different values of $\omega$. In particular,
for the ab initio model, even considering one branch, there are many modes
that have the same $C_{\omega}v_{\omega}$ value. The existence of many such
modes implies that a number of other relaxation-time distributions exist which
could have the same $F(\Lambda)$ as the ab initio model. Note that the results
provided in section 3.1 are based on Holland model.

(a)

(b)

Figure 8: The plot of $C_{\omega}v_{\omega}$ versus frequency for the ab
initio (a) and the Holland model (b).

## Appendix C: Bayesian inference

The likelihood function that we considered in our study is defined as

$p\Big{(}\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L}\big{|}{\bf U}\Big{)}={\cal
N}\left(\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L};\\{T_{\text{BTE}}(t,{\bf
x}=0;L)\\}_{t,L},\sigma_{m}^{2}\right),$ (15)

in which $\sigma_{\text{m}}=0.01$ K. The prior distribution $\pi({\bf U})$ of
the parameter U — its distribution before the observations— here is taken to
be a Laplace distribution,

$\pi({\bf U})=\prod_{i}\frac{1}{2{\bf\Sigma}_{i}}\exp{\left(-\frac{|{\bf
U}_{i}-{\bf\mu}_{i}|}{{\bf\Sigma}_{i}}\right)},$ (16)

whose variance for each component ${\bf U}_{i}$ is $2{\bf\Sigma}_{i}^{2}$ and
$\mu_{i}$ is its mean (corresponding to the correct solution). The variance
and mean values for different components are taken to be similar to those in
[32]. Here, $\\{T_{\text{m}}(t,{\bf x}=0;L)\\}_{t,L}$ are the set of all
temperature measurements (observations). Equation (15) implies that here we
are assuming that the temperature measurements are corrupted by a noise of
Gaussian form with the standard deviation of $\sigma_{\text{m}}=0.01$ K. The
value of $\sigma_{m}$ is consistent with the noise in the common thermal
spectroscopy experiments [11]. The choice of Laplace distribution for the
prior is motivated by two observations:

  * •

Here we are more interested in the sensitivity of the reconstructed property,
the relaxation-times function, to the noise in the temperature measurement.
Therefore, instead of the least-informative uniform distribution, we are
assuming a prior with higher probability around the correct solution. This is
motivated by our observation that this optimization (global optimization as
well as optimizations in [7, 9, 11]) converges to a solution close to the
correct function. In other words, the Bayesian analysis is a complement to the
uniqueness analyses, in the sense that while the results provided in figure 4
show that the objective function has a global minimum, in the Bayesian
analysis we study how sensitive the objective function is to the measurement
noise in the vicinity of that global minimum.

  * •

Due to our choice of objective function ${\cal L}$ ($L$-1 norm) we used here
and previously in [7, 9, 11], the Laplace distribution is expected to be a
more accurate representation of the behavior of the objective function
compared to other common distributions such as Gaussian distribution (compare
the similarity between Laplace distribution and the shape of the objective
function at the minimum in the supplementary document of [11]). The result
corresponding to a Gaussian distribution counterpart of this distribution can
be found in [32]

The posterior distribution was sampled using the Metropolis-Hasting (MH)
algorithm [34]. The parameters of the algorithm (e.g., the variance of the
proposal distribution) are similar to those of [32]. The acceptance rate for
the chosen proposal distribution is 0.21 which is close to the optimal value
of 0.23 for multi-dimensional variables [35]. We have used a total of about
$N_{\text{MH}}=7\times 10^{5}$ sample points; the first $10\%$ of all samples,
approximately, were used as burn-in samples.