PHYSICAL REVIEW B 90, 064302 (2014)
Analytical interpretation of nondiffusive phonon transport in thermoreflectance
thermal conductivity measurements
K. T. Regner,
1
A. J. H. McGaughey,
1,2
and J. A. Malen
1,2,*
1
Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
2
Department of Materials Science and Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
(Received 19 June 2014; revised manuscript received 1 August 2014; published 20 August 2014)
We derive an analytical solution to the Boltzmann transport equation (BTE) to relate nondiffusive thermal
conductivity measurements by thermoreflectance techniques to the bulk thermal conductivity accumulation
function, which quantifies cumulative contributions to thermal conductivity from different mean free path energy
carriers (here, phonons). Our solution incorporates two experimentally defined length scales: thermal penetration
depth and heating laser spot radius. We identify two thermal resistances based on the predicted spatial temperature
and heat flux profiles. The first resistance is associated with the interaction between energy carriers and the surface
of the solution domain. The second resistance accounts for transport of energy carriers through the solution
domain and is affected by the experimentally defined length scales. Comparison of the BTE result with that from
conventional heat diffusion theory enables a mapping of mean-free-path-specific contributions to the measured
thermal conductivity based on the experimental length scales. In general, the measured thermal conductivity will
be influenced by the smaller of the two length scales and the surface properties of the system. The result is used
to compare nondiffusive thermal conductivity measurements of silicon with first-principles-based calculations of
its thermal conductivity accumulation function.
DOI: 10.1103/PhysRevB.90.064302 PACS number(s): 65.40.−b, 44.05.+e, 63.20.−e
I. INTRODUCTION
Nondiffusive thermal transport occurs when length or time
scales of a system are on the order of the mean free paths
(MFPs) or lifetimes of the energy carriers. As a result, a local
equilibrium temperature cannot be defined and the thermal
transport properties of the system can no longer be taken
as the bulk values. When system boundaries are decreased
below energy carrier MFPs, nondiffusive transport can be
described with a reduced, effective thermal conductivity [1–5].
Heat dissipation in light emitting diodes and transistors is
adversely impacted by reductions in thermal conductivity,
while thermoelectric energy conversion devices benefit.
Determination of the relationship between system dimen-
sions and effective thermal conductivity has been a research
focus for over 20 years and requires two fundamental pieces
of information: (i) the intrinsic (i.e., bulk) MFP-dependent
contributions of energy carriers to thermal conductivity
[6–8] k
and (ii) the relationship between system dimensions
and the modified MFPs of the energy carriers [9,10]. In
semiconducting materials, the former can be described by the
thermal conductivity accumulation function for phonons [11],
k
accum
, which identifies cumulative contributions to thermal
conductivity from phonons having a MFP less than or equal to
the length scale
∗
. Under the isotropic assumption,
k
accum
(
∗
) =
∗
0
k
d =
∗
0
1
3
C
()v()d. (1)
Here, is MFP, C
is volumetric heat capacity per unit MFP,
and v is group velocity. Thermal conductivity accumulation
functions have been determined theoretically for bulk and
nanostructured materials using analytical scattering relation-
ships [ 10], molecular dynamics simulations with empirical
*
jonmalen@andrew.cmu.edu
potentials [7], and by first-principles calculations [8,12,13],
but require experimental validation.
Recent attempts have been made to experimentally measure
k
accum
by inducing nondiffusive thermal transport through
varying an experimentally controllable length scale L
c
in a
range comparable to phonon MFPs. Techniques include tran-
sient thermal grating (TTG), where L
c
is the period of a pulsed
optical grating that induces a spatially periodic temperature
profile [14] and thermoreflectance techniques including time
domain thermoreflectance (TDTR) and broadband frequency
domain thermoreflectance (BB-FDTR), where the experimen-
tal length scales are the spot size of a heating laser and the
thermal penetration depth of a temporally sinusoidal laser heat
flux [6,15–18]. An effective thermal conductivity as a function
of L
c
is found by interpreting nondiffusive measurements with
a s olution to the heat diffusion equation.
Initially, the interpretation to obtain k
accum
was that energy
carriers with >L
c
do not contribute to the experimentally
measured thermal conductivity k
exp
and energy carriers with
L
c
fully contribute, as they would in a purely diffusive
regime [6,15,16,18]. Mathematically, this assumption takes
the form
k
exp
(L
c
) =
L
c
0
k
d. (2)
This mapping between L
c
and MFP contributions to the
effective thermal conductivity leads to accumulation functions
that are consistent with first-principles predictions in silicon
and gallium arsenide [15,16,18] but lacks rigorous justifica-
tion. More generally,
k
exp
(L
c
) =
∞
0
S(,L
c
)k
d, (3)
where S(, L
c
) is known as the suppression function. In the
simple interpretation in Eq. (2), S(, L
c
) is a step function
1098-0121/2014/90(6)/064302(10) 064302-1 ©2014 American Physical Society
K. T. REGNER, A. J. H. MCGAUGHEY, AND J. A. MALEN PHYSICAL REVIEW B 90, 064302 (2014)
from 1 to 0 at = L
c
. But discrepancies between BB-FDTR
[16] and TDTR [6] results using Eq. (2) demand a deeper
understanding of the suppression function.
Comparison of analytical [19] and numerical solutions
[20–22] of the Boltzmann transport equation (BTE) to the heat
diffusion equation for TTG leads to the functional dependence
of the suppression function on L
c
and MFP and reconciles
nondiffusive TTG measurements and k
accum
. Although the
form of the suppression function has been identified for
TTG, a new analysis is required for BB-FDTR and TDTR
since the experimental setups are physically different, i.e.,
L
c
is different. Ding et al. predicted suppression due to
spot size in TDTR using a Monte Carlo–based numerical
solution to the BTE [23], but neither suppression due to
thermal penetration depth nor analytical analyses for these
experiments have been demonstrated in the literature. Three
important questions remain unresolved: (1) What is the f orm of
thermal-penetration-depth-based suppression? (2) What is its
interplay with spot-size-based suppression? (3) Under what
circumstances can BB-FDTR and TDTR measurements be
interpreted with the conventional heat diffusion equation?
In this work, we derive an analytical suppression function
for thermoreflectance techniques by solving the BTE. In ther-
moreflectance techniques, there are two experimental length
scales: (1) the thermal penetration depth L
p
=
√
2k/C,
which characterizes the exponential decay length of the
temperature amplitude into a solid with thermal conductivity k
and volumetric heat capacity C due to sinusoidal laser heating
with angular frequency at the surface, and (2) the e
−2
radius of the Gaussian laser spot, r
o
. The presence of r
o
in
thermoreflectance experiments necessitates a comparison of
length scales rather than the time scales 1/ and phonon
lifetimes. In Secs. II and III, we account for both experimental
length scales in our expression for the suppression function.
The results are used in Sec. IV to interpret nondiffusive
measurements of phonon transport in silicon by BB-FDTR
and TDTR, although our solution does not account for the
multiple time scales in TDTR that arise from using a pulsed
laser.
II. SUPPRESSION FUNCTION IN A PLANAR GEOMETRY
As shown in Fig. 1(a), we first consider a planar medium
with a temporally oscillating surface temperature with angular
frequency and amplitude T
s
= 1K, such that T (x = 0,t) =
T
s
e
it
. Because we are solving for deviations from the
mean temperature, for convenience we define the temperature
T (x →∞,t) = T
∞
= 0 K. The one-dimensional (1D) nature
of this problem will yield an analytical solution that provides
insight into the functional dependence of the suppression
function on thermal penetration depth.
We begin with the gray, 1D BTE for phonons in Cartesian
coordinates under the relaxation time approximation in an
isotropic medium [24,25],
1
v
∂n
∂t
+ μ
∂n
∂x
=
n
e
− n
τv
, (4)
where the nonequilibrium distribution function n(x,t,μ)isthe
phonon energy density per unit phonon frequency per unit
solid angle and equals ωD(ω)g(x,t,μ)/4π . Here, is the
n(x, t, μ) =
n
+
(x, t), 0 < μ ≤ 1
n
-
(x, t), -1 ≤ μ ≤ 0
εn
e
(T
s
) = εn
-
(x = 0, t)ρn
-
(x = 0, t)
T(x = 0, t) = T
s
e
iΩt
x
θ
(a)
C
ω
T
s
4π
θ
T(r = r
o
, t) = T
s
e
iΩt
n(r, t, μ) =
n
+
(r, t), 0 < μ ≤ 1
n
-
(r, t), -1 ≤ μ ≤ 0
εn
e
(T
s
) = ε
n
-
(r = r
o
, t)
ρn
-
(r = r
o
, t)
r
o
(b)
r
C
ω
T
s
4π
FIG. 1. (Color online) Schematic diagrams for ( a) the 1D planar
system (Sec. II) and (b) the spherically symmetrical system (Sec. III),
both with oscillating surface temperatures. Here, μ is the directional
cosine, μ = cosθ . The parameters ε and ρ are the phonon emissivity
and reflectivity and are discussed further in Sec. IV.
reduced Planck constant, ω is the phonon frequency, D(ω)
is the phonon density of states, g(x,t,μ) is the occupation
function, n
e
(x,t) is the equilibrium distribution function and is
specified for phonons when g is the Bose-Einstein distribution
g
BE
, τ is the gray lifetime /v, v is the frequency-independent
phonon group velocity (i.e., sound velocity), and μ is the
directional cosine ( μ = cos θ) that accounts for the velocity
of phonons traveling at an angle θ from the x direction
[see Fig. 1(a)]. For small temperature variation, n
e
(x,t) ≈
ωD(ω)
dg
BE
dT
|
x,t
T (x,t)/4π = C
ω
T (x,t)/4π, where C
ω
is the
volumetric heat capacity per unit frequency and T (x,t)is
the departure from T
∞
= 0[19,20,27]. Thus, we solve for the
deviations from the equilibrium distribution function, which
are r elated to deviations of temperature from T
∞
.
Since the oscillating surface temperature determines the
temporal behavior of the solution, we separate variables such
that n(x,t,μ) =
¯
n(x,μ)e
it
, where
¯
n is the component of n
that is only a function of x and μ. Substituting into Eq. (4)
yields
μ
d
¯
n
dx
+
i
v
+
1
τv
¯
n =
¯
n
e
τv
. (5)
The difficulty in solving Eq. (5) arises from the fact that
we must account for phonons traveling over all directions
064302-2
ANALYTICAL INTERPRETATION OF NONDIFFUSIVE . . . PHYSICAL REVIEW B 90, 064302 (2014)
μ. For TTG, Collins et al. demonstrated a Volterra integral
solution to a BTE of similar form [19], but the dependence on
μ in our formulation leads to a divergent integral. Henceforth,
we follow a two-flux procedure similar to that of the Milne-
Eddington approximation for radiative heat transfer [28]. This
method involves taking the zeroth and first moments of Eq. (5),
i.e., Eq. (5) is integrated over all directions after multiplication
with μ
0
= 1 (zeroth moment) and μ
1
= μ (first moment). The
distribution moments are defined as
¯
n
l
(x) = 2π
1
−1
¯
n(x,μ)μ
l
dμ, l = 0,1,.... (6)
Furthermore, the distribution function is assumed to be
isotropic over the upper and lower hemispheres such that
¯
n
+
(x) ≡
¯
n(x,0 <μ 1) and
¯
n
−
(x) ≡
¯
n(x, − 1 μ 0)
[see Fig. 1(a)] [28]. From Eq. (6), the zeroth and first moments
are
¯
n
0
= 2π(
¯
n
+
+
¯
n
−
) = 3
¯
n
2
and
¯
n
1
= π(
¯
n
+
−
¯
n
−
), which
can be physically related to temperature and heat flux [28].
Applying the two-flux method to Eq. (5) yields a coupled set
of equations:
d
¯
n
0
dx
+ 3
i
v
+
1
τv
¯
n
1
= 0, (7a)
d
¯
n
1
dx
+
i
v
¯
n
0
= 0. (7b)
In formulating Eqs. (7a) and (7b), we use conservation
of energy for a gray medium to determine the equilibrium
distribution
¯
n
e
in terms of
¯
n
0
as (Ref. [29])
¯
n
e
=
1
2
1
−1
¯
ndμ =
¯
n
0
4π
. (8)
This coupled set of ordinary, linear, homogeneous differ-
ential equations is an eigenvalue problem and has a solution
of the form [
¯
n
0
¯
n
1
] = c
1
v
1
e
−λx
+ c
2
v
2
e
λx
, where c
1
and c
2
are
constants to be determined by the boundary conditions, ±λ
are the eigenvalues, and v
1
and v
2
are the eigenvectors. Since
the spatial domain is semi-infinite, c
2
= 0 because
¯
n
0
and
¯
n
1
cannot increase unbounded. The boundary condition at x = 0
is depicted schematically in Fig. 1(a) and is [30,31]
¯
n
+
(x = 0) = ε
C
ω
T
s
4π
+ ρ
¯
n
−
(x = 0), (9)
where ε and ρ are the phonon emissivity and reflectivity,
both of which will be discussed in further detail in Sec. IV.
Physically, Eq. (9) states that the total energy carried by
phonons traveling in the positive x direction at the surface is
equal to the sum of the energy carried by phonons emitted due
to the induced surface temperature T
s
and the energy carried by
phonons traveling in the negative x direction that are reflected
from the surface.
By solving the system of equations and integrating over
all phonon frequencies, the spatial temperature and heat flux
profiles are found to be
¯
T
BTE
(x) =
∫
∞
0
¯
n
0
(x)dω
∫
∞
0
C
ω
dω
=
εT
s
(1 + ρ)
4iβ
3η
+ (1 − ρ)
exp
−
η
L
p
x
, (10a)
¯
q
BTE
(x) =
∞
0
v
¯
n
1
(x)dω
=
εT
s
Cv
2(1 + ρ) + (1 − ρ)
3η
2iβ
exp
−
η
L
p
x
, (10b)
where η =
√
2i − 2τ,ß= /L
p
, and k
bulk
=
1
3
Cv
2
τ . Since
we use the gray approximation,
¯
n
0
and
¯
n
1
are independent of ω
and the integral over ω only changes C
ω
to the total volumetric
heat capacity, i.e.,
∞
0
C
ω
dω = C. To generate figures in
this section and Sec. III, we use properties of bulk silicon
(C = 1.65×10
6
Jm
−3
K
−1
,k
bulk
= 145 W m
−1
K
−1
, and
v = 8430 ms
−1
)[32,33] and determine L
p
using k
bulk
.
The magnitudes of the spatial temperature profiles from the
diffusion solution [
¯
T
diff
(x) = T
s
exp(−
√
2ix
L
p
)] and BTE solution
for ε = 1−ρ = 1 and /L
p
= 1 are shown in Fig. 2(a).
The spatial temperature profile from the diffusion solution
is a continuous exponential decay where the diffusive thermal
resistance can be defined as R
diff,x
= (T
s
− T
∞
)/
¯
q
diff
(
x = 0
)
.
The real part of the exponential in Eqs. (10a) and (10b)
represents the BTE prediction of penetration depth L
p-BTE
,
(b)
10
−2
10
−1
10
0
10
1
10
−11
10
−10
10
−9
10
−8
10
−7
10
−4
10
−2
10
0
10
2
Λ/L
p
τΩ
Thermal Resistance
(m
2
-K/W)
|R
diff,x
|
|R
BTE,x
|
|R
i,x
|
|R
ε
|
(a)
ΔT
ε
ΔT
i
|T
diff
(x)|
|T
BTE
(x)|
T
∞
= 0
Temperature (K)
T
s
= 1
x/L
p
012345
Λ /
L
p
= 1
FIG. 2. (Color online) 1D planar geometry with temporally os-
cillating surface temperature and ε = 1−ρ = 1. (a) Magnitude of the
spatial temperature profiles from the diffusion and BTE solutions for
/L
p
= 1. The BTE solution has two distinct regions that correspond
to two distinct thermal resistances. (b) Magnitude of the thermal
resistances R
diff,x
and R
BTE,x
= R
ε
+ R
i,x
plotted as a function of
/L
p
and τ.
064302-3
K. T. REGNER, A. J. H. MCGAUGHEY, AND J. A. MALEN PHYSICAL REVIEW B 90, 064302 (2014)
which can be written as
L
p-BTE
=
L
p
(1 + τ
2
2
)
1/2
− τ
. (11)
When τ 1,L
p-BTE
= L
p
and
¯
T
BTE
(x) collapses to
¯
T
diff
(x),
but when τ 1,L
p-BTE
→∞, which indicates purely
ballistic transport. Thus, as /L
p
increases, the temperature
decay rate predicted by the BTE decreases.
The spatial temperature profile from the BTE solution
indicates two distinct regions: a surface temperature jump of
T
ε
and a spatial temperature decay spanning T
i
. When ε
= 1−ρ, the total thermal resistance from the BTE solution
R
BTE,x
is comprised of two parts,
R
ε
=
T
ε
¯
q
BTE
(
x = 0
)
=
4 − 2ε
εCv
, (12a)
R
i,x
=
T
i
¯
q
BTE
(
x = 0
)
=
L
2
p
2i
+
2
3
k
bulk
, (12b)
such that
R
BTE,x
= R
ε
+ R
i,x
. (12c)
The thermal resistances in Eqs. (12a), (12b), and (12c)
are complex. Complex thermal resistances are analogous to
impedance in alternating current circuit analysis. In the plots
throughout this paper, we plot the magnitude of such complex
thermal resistances.
The magnitude of the terms R
ε
, R
i,x
, and R
BTE,x
are plotted
in Fig. 2(b) as a function of /L
p
and τwith ε = 1−ρ
= 1 and are compared to the magnitude of R
diff,x
.Theterm
R
ε
is a resistance that arises from the interaction between
the surface and ballistic phonons originating within one MFP
of the surface and is associated with the surface temperature
jump in BTE [31,34–36] and radiative transfer [37] problems.
The term R
ε
is independent of any experimentally controllable
length scale but is always present. The term R
i,x
is intrinsic to
the material and accounts for transport of phonons associated
with two length scales: L
p
and . It should be noted that R
i,x
says nothing about the surface properties (i.e., R
i,x
is not a
function of ε). Thus, when /L
p
1,R
BTE,x
= R
ε
+ R
i,x
=
R
ε
+ L
p
/(
√
2ik
bulk
) ≈ L
p
/(
√
2ik
bulk
) and the BTE thermal
resistance converges to the diffusive thermal resistance be-
cause R
i,x
dominates R
ε
. However, as the phonon MFP
approaches L
p
, the second term in R
i,x
and the R
ε
term become
non-negligible and the BTE thermal resistance becomes larger
than the diffusive thermal resistance. In the ballistic limit
(/L
p
1),R
BTE,x
= R
ε
+ R
i,x
= R
ε
+ /(
√
3k
bulk
) and
becomes independent of . It should be noted that the total
thermal resistance is independent of whether a temporally
oscillating surface temperature or heat flux is imposed, the
latter of which is more consistent with the experiments.
As in the analysis of the experimental measurements, we
can now determine the effective thermal conductivity k
eff
that
equates the complex diffusive thermal resistance (R
diff,x
=
1/
√
iCk
eff
) to the complex thermal resistance determined
by the BTE,R
BTE,x
[21,31]. Since, by definition, T
s
is identical
in both systems, this procedure is equivalent to equating
surface heat fluxes from the diffusion and BTE solutions.
Furthermore, similar functional forms of t he BTE and diffusion
solutions suggest that interpreting nondiffusive transport with
an effective, suppressed k is reasonable. We define the
suppression function for this planar geometry S
x
(,L
p
,ε,ρ)
as the fractional contribution to thermal conductivity made by
a phonon with a MFP of in a thermoreflectance experiment
with , ε, and ρ and is
S
x
(,L
p
,ε,ρ) =
k
eff
k
bulk
=
9ε
2
2iβ
2
2(1 +ρ) +(1 −ρ)
3η
2iβ
2
. (13)
It should be noted that S
x
(,L
p
,ε,ρ) is complex. Thus
the phase angle of the suppression function influences the
observed phase angle in thermoreflectance experiments, ulti-
mately influencing the value of thermal conductivity obtained.
In plots of the suppression function throughout the paper, we
plot its magnitude.
In Figs. 3(a) and 3(b), we plot the magnitude of the thermal
resistance of the system from the BTE and diffusion solutions
(a)
(b)
p
10
−2
10
−1
10
0
10
1
10
−11
10
−10
10
−9
10
−8
10
−7
10
−4
10
−2
10
0
10
2
Λ/L
τΩ
Diffusion Solution
ε = 0.1
ε = 0.5
ε = 1
|R
diff,x
| or |R
BTE,x
| (m
2
-K/W)
10
−2
10
−1
10
0
10
1
0
0.2
0.4
0.6
0.8
1.0
10
−4
10
−2
10
0
10
2
Λ/L
p
τΩ
ε = 1 Eq. (2)
P
1
, || Plates
ε = 0.5
ε = 0.1
|S
x
|
FIG. 3. (Color online) 1D planar geometry with temporally os-
cillating surface temperature and ε = 1−ρ = 1, 0.5, and 0.1. (a)
Magnitude of the thermal resistance from the diffusion and BTE
solutions vs /L
p
and τ. The BTE predicts a higher thermal
resistance than the diffusion solution, which can be accounted for by
reducing the effective thermal conductivity in the diffusion solution.
(b) Magnitude of the suppression function plotted as a function of
/L
p
and τ. These curves are compared to the P
1
solution to the
BTE for parallel, black, isothermal plates and to the step function
suppression function [Eq. (2)] [6,16,18].
064302-4
ANALYTICAL INTERPRETATION OF NONDIFFUSIVE . . . PHYSICAL REVIEW B 90, 064302 (2014)
and the magnitude of S
x
(,L
p
,ε,ρ) as a function of /L
p
and
τfor ε = 1−ρ = 1, 0.5, and 0.1. The suppression function
[Fig. 3(b)] accounts f or the increase in thermal resistance
compared to the diffusion solution [Fig. 3(a)], and reduces the
effective thermal conductivity of the material. The suppression
function is different than that previously assumed [i.e., a step
function; see Eq. (2)] [6,16,18] in that phonons with /L
p
< 1
contribute less and phonons with /L
p
> 1 contribute more
near /L
p
= 1.
The effect of changing ε is highlighted in Figs. 3(a) and
3(b). In our BTE solution, the resistance associated with
the surface temperature jump R
ε
= (4 − 2ε)/εCv (for ρ =
1−ε) is independent of any experimentally controllable length
scale, i.e., L
p
. Consequently, this resistance is always present
and of the same magnitude but only becomes non-negligible
when R
i,x
is sufficiently small, which happens when the
penetration depth is on the order of or smaller than the MFP.
Decreasing ε increases R
ε
, increasing the surface temperature
jump, and hastening the onset of suppression. This fact can be
qualitatively understood with an analogy to radiative transfer,
i.e., the energy transfer rate from an isothermal gray surface
will be less than that from an isothermal black surface at a
given surface temperature. Reducing the phonon emissivity
reduces the number of phonons emitted from the surface and
hence reduces the energy transfer away from the surface,
increasing the thermal resistance and reducing the effective
thermal conductivity of the material in the nondiffusive regime.
Furthermore, it is reasonable that emissivity is related to the
interface resistance between the transducer and substrate in
a thermoreflectance experiment considering that emissivity
affects the size of the surface temperature jump [38]. The
effect of changing ε and ρ will be revisited in Sec. IV.
To verify the behavior of our suppression function, we
compare it to a solution to the gray BTE for two infinite,
parallel, black (ε = 1), isothermal plates. This scenario is
similar to our problem except that we consider an oscillating
surface temperature that defines our length scale L
p
.The
solution to this problem is obtained using the P
1
approximation
and is plotted against the ratio of and plate separation
distance in Fig. 3(b). A similar trend instills confidence in
our solution and suggests that although L
p
is not a physical
boundary, it similarly suppresses contributions of phonons to
thermal transport.
III. SUPPRESSION FUNCTION IN
A SPHERICAL GEOMETRY
In BB-FDTR and TDTR experiments, there are two relevant
length scales: the thermal penetration depth and the spot size of
the heating laser. Thus, in order to obtain an accurate suppres-
sion function for relating thermoreflectance measurements to
k
accum
, both length scales should be incorporated. The most
accurate solution would involve solving the spectral BTE in
cylindrical coordinates, under conditions of radially symmetric
Gaussian surface heating. While other studies have reached
numerical solutions to similar problems [23,39], it is our goal
to reach an analytical solution for a simpler problem.
As depicted in Fig. 1(b), we consider a sphere with
radius r
o
embedded in an infinite medium with tempera-
ture T
(
r →∞,t
)
= T
∞
= 0 K and a temporally oscillating
surface temperature at the sphere-medium interface. Solving
the BTE within the medium will provide a solution that is
dependent on L
p
due to the periodic nature of the surface
temperature as well as the effect of spot size, which can be
captured by varying the radius of the embedded sphere. We
note that Chen solved a similar problem for a sphere with
steady-state heating [31]. While this geometry is not an exact
representation of a thermoreflectance experiment, the spherical
symmetry (1D in the radial direction) of the problem allows
us to derive an analytical solution for the suppression function
that is dependent on L
p
and r
o
.
We begin with the 1D, gray BTE under the relaxation time
approximation in spherical coordinates in the radial direction
r (Ref. [24]),
1
v
∂n
∂t
+ μ
∂n
∂r
+
1 − μ
2
r
∂n
∂μ
=
n
e
− n
τv
. (14)
The μ-dependence in Eq. (14) can be eliminated using the
method of spherical harmonics (P
N
approximation), which is a
generalization of the Milne-Eddington approximation and has
been thoroughly studied in spherically symmetrical geometries
in radiative transfer [28,40–42]. The method involves reducing
the governing equation into a set of N simpler partial
differential equations by taking advantage of the orthogonality
of spherical harmonics. Applying the P
N
approximation to
Eq. (14) yields
1
v
∂n
l
∂t
+
l + 1
2l + 1
∂n
l+1
∂r
+
l
2l + 1
∂n
l−1
∂r
+
(l + 1)(l + 2)
r(2l + 1)
n
l+1
−
l(l − 1)
r(2l + 1)
n
l−1
+
n
l
τv
=
n
0
τv
δ
0l
, (15)
where l = 0,1,2,...,N and δ
0l
is the Kronecker delta. In the
limit where N → , the exact solution is obtained. Here
we use the P
1
approximation, which is accurate for scattering
media at large optical thicknesses with decreasing accuracy
as the optical thickness is decreased [28]. For our problem,
large optical thicknesses correspond to L
p
.UsingtheP
1
approximation and separating variables in a similar fashion as
Eq. (5), Eq. (14) reduces to
d
¯
n
0
dr
+ 3
i
v
+
1
τv
¯
n
1
= 0, (16a)
d
¯
n
1
dr
+
i
v
¯
n
0
+
2
r
¯
n
1
= 0. (16b)
By employing an analogous boundary condition as used for
the planar solution, i.e.,
¯
n
+
(
r = r
o
)
= ε
C
ω
T
s
4π
+ ρ
¯
n
−
(
r = r
o
)
,
we obtain closed-form solutions for the spatial temperature
and heat flux profiles for r r
o
,
¯
T
BTE
(r) =
r
o
εT
s
r
(1 + ρ)
4i(βη+)
3η
2
+ (1 − ρ)
×exp
−
η
L
p
(r − r
o
)
, (17a)
¯
q
BTE
(r) =
βη
r
r
o
+
βη +
r
2
o
εT
s
Cv
r
2
2(1 + ρ) + (1 − ρ)
3η
2
2i(βη+)
×exp
−
η
L
p
(r − r
o
)
, (17b)
064302-5
![FIG. 5. (Color online) Spherical particle embedded in an infinite medium with oscillating temperature at the surface of the sphere (r = ro) with ε = 1−ρ = 1. (a) Magnitude of the thermal resistance from the diffusive and BTE solutions vs /ro for different values of /Lp. (b) Magnitude of the suppression function plotted as a function of /ro for different values of /Lp. (c) Comparison of Eq. (18) when /Lp = 0 for a particle with radius ro and a particle with radius 3ro, the exact solution for a sphere with steady-state heating from Ref. [31], and the suppression function found numerically by solving the spectral BTE for a Gaussian-shaped laser spot from Ref. [23]. Using a particle radius of 3ro in our analytical result compares well with numerical results from Ref. [23].](/figures/fig-5-color-online-spherical-particle-embedded-in-an-1hj7ltu7.webp)
![FIG. 6. (Color online) Comparison of thermal conductivity measurements and predicted kexp for silicon. (a) kaccum from first-principles calculations (solid lines) is transformed into predicted kexp vs 3ro (dashed lines) by Eq. (3) using the suppression function from Eq. (18). A value of ε = 0.88 results in the best fit to TDTR measurements at T = 300 K and T = 80 K from Ref. [23] with a heating frequency of 106 Hz. Predicted kexp vs 3ro for a heating frequency of 107 Hz with ε = 0.88 is shown for comparison. In TDTR, the transducer/substrate interface is Al/Si. (b) kaccum from first-principles calculations is transformed into kexp vs Lp by Eq. (3) using the suppression function from Eq. (18). A value of ε = 0.6 results in the best fit to BB-FDTR measurements (circles) at T = 300 K and T = 80 K from Ref. [16]. In BB-FDTR, the transducer/substrate interface is Cr/Si. Predicted kexp vs Lp with ε = 0.88 is compared to TDTR measurements from Ref. [6] (diamonds).](/figures/fig-6-color-online-comparison-of-thermal-conductivity-1ad29827.webp)
![FIG. 3. (Color online) 1D planar geometry with temporally oscillating surface temperature and ε = 1−ρ = 1, 0.5, and 0.1. (a) Magnitude of the thermal resistance from the diffusion and BTE solutions vs /Lp and τ . The BTE predicts a higher thermal resistance than the diffusion solution, which can be accounted for by reducing the effective thermal conductivity in the diffusion solution. (b) Magnitude of the suppression function plotted as a function of /Lp and τ . These curves are compared to the P1 solution to the BTE for parallel, black, isothermal plates and to the step function suppression function [Eq. (2)] [6,16,18].](/figures/fig-3-color-online-1d-planar-geometry-with-temporally-2zolvk3n.webp)