Analytical interpretation of nondiffusive phonon transport in thermoreflectance thermal conductivity measurements

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.


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][2][3][4][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 dimensions 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][7][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, 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 relationships [10], molecular dynamics simulations with empirical * jonmalen@andrew.cmu.edupotentials [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 transient 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 experimental length scales are the spot size of a heating laser and the thermal penetration depth of a temporally sinusoidal laser heat flux [6,[15][16][17][18].An effective thermal conductivity as a function of L c is found by interpreting nondiffusive measurements with a solution 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 ( 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 justification.More generally, where S( , L c ) is known as the suppression function.In the simple interpretation in Eq. ( 2), S( , L c ) is a step function 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][21][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 form 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 thermoreflectance 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 = 1 K, such that T (x = 0,t) = T s e i t .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], where the nonequilibrium distribution function n(x,t,μ) is the phonon energy density per unit phonon frequency per unit solid angle and equals ωD(ω)g(x,t,μ)/4π .Here, is the 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 related 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 i t , where n is the component of n that is only a function of x and μ.Substituting into Eq.( 4) yields The difficulty in solving Eq. ( 5) arises from the fact that we must account for phonons traveling over all directions μ.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 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 n0 = 2π ( n+ + n− ) = 3 n2 and n1 = π ( 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: In formulating Eqs.(7a) and (7b), we use conservation of energy for a gray medium to determine the equilibrium distribution ne in terms of n0 as (Ref.[29]) This coupled set of ordinary, linear, homogeneous differential equations is an eigenvalue problem and has a solution of the form [ n0 n1 ] = 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 n0 and n1 cannot increase unbounded.The boundary condition at x = 0 is depicted schematically in Fig. 1(a) and is [30,31] n+ 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 where η = √ 2i − 2τ , ß = /L p , and k bulk = 1 3 Cv 2 τ .Since we use the gray approximation, n0 and n1 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 [32,33] and determine L p using k bulk .
The magnitudes of the spatial temperature profiles from the diffusion solution ] 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 Temperature (K) which can be written as When τ 1,L p-BTE = L p and TBTE (x) collapses to Tdiff (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, 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 .The term 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][35][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 bulk ) and the BTE thermal resistance converges to the diffusive thermal resistance because 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 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/ √ i Ck 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 the 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 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, ultimately influencing the value of thermal conductivity obtained.In plots of the suppression function throughout the paper, we plot its magnitude.
In Figs.  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 for 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 suppression 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 temperature 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]), 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][41][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 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 .Using the P 1 approximation and separating variables in a similar fashion as Eq. ( 5), Eq. ( 14) reduces to 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 , where = /r o .The suppression function is found by determining k eff of the infinite medium that equates the complex thermal resistance from the diffusion solution [43] to the complex thermal resistance defined by the BTE, which is equivalent to equating surface heat fluxes [21,31] and is where In Fig. 4(a), the magnitude of the complex thermal resistance is plotted as a function of /L p and τ at different values of /r o with ε = 1−ρ = 1 for both the diffusion and BTE solutions.The thermal resistance from the diffusion equation (solid lines) highlights the interplay between r o and L p .When L p r o , the solution converges to the planar solution [Fig.3(a)], and when L p r o , the diffusive thermal resistance becomes independent of /L p .
Similar to the planar solution, the total thermal resistance from the BTE is the sum of a surface component R ε = (4 − 2ε)/εCv (for ε = 1−ρ), which is the same as for the planar solution, and an intrinsic component R i,r , As in the planar solution, R i,r includes no effect from the surface properties (R i,r is not a function of ε when ε = 1−ρ).For a given value of ε, R i,r converges to the diffusion solution when /L p 1 and asymptotes to /( √ 3k bulk ) when /L p 1.But since the diffusive resistance decreases with increasing /r o when /L p 1, R ε becomes nonnegligible, and even dominates, when r o is commensurate or smaller than the MFP.Because the total thermal resistance is the sum of R ε and R i,r , the BTE and diffusion solutions do not converge when /L p 1 at larger values of /r o .When /r o = 0, the BTE solution converges to the planar solution from Eq. ( 10), as shown in Fig. 3(a).
The magnitude of the suppression function S r ( ,L p ,r o ,ε,ρ) is plotted in Fig. 4(b) for ε = 1−ρ = 1.In the limit when r o → ∞, the solution converges to the planar solution given in Eq. ( 13) and shown in Fig. 3(b).Changes in the suppression function with /r o over all /L p illustrate the interactions between the two length scales.In general, the smaller of L p or r o dominates suppression.For example, when /L p 1, suppression is solely due to decreasing particle radius and is consistent with the TDTR experimental measurements by Minnich et al. of k vs r o that were independent of heating frequency [15] and Chen in the case of steady-state heating [31].According to the BTE solution, if either L p or r o are much smaller than the phonon MFP, that phonon will not contribute to k exp .Under these circumstances, BB-FDTR and TDTR are inadequate for measuring the bulk thermal conductivity of a material.
The magnitude of the thermal resistance from the BTE and diffusion solutions and the magnitude of S r ( ,L p ,r o ,ε,ρ) are plotted in Figs.5(a) and 5(b) as a function of /r o for different values of /L p with ε = 1−ρ = 1.As heating frequency increases ( /L p increases), additional suppression occurs from L p , even at very low /r o .In Fig. 5(c), we compare our analytical solution for S r in the low limit ( /L p = 0) with Chen's exact solution from Ref. [31] for a sphere with steady-state heating and the suppression function from Ref. [23] found numerically by solving the spectral BTE for a Gaussian-shaped laser spot.Due to our use of the P 1 approximation, we find that Eq. ( 18) and the exact solution for a sphere with steady-state heating from Ref. [31] differ by a factor of 2 on the horizontal axis.We assert that this factor is not significant considering that the range of MFP spans four orders of magnitude in typical crystalline semiconductors [7,8,12].We also find that using a value of 3r o in Eq. ( 18) yields a suppression function that compares well with the suppression function from Ref. [23].We expect that there should be a correction factor to the spot size in our suppression function as a result of the geometry we have   18) when /L p = 0 for a particle with radius r o and a particle with radius 3r o , 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 3r o in our analytical result compares well with numerical results from Ref. [23].
chosen to generate an analytical solution, i.e., we approximate our spot as a finite sphere in an infinite medium while the actual experimental geometry is a Gaussian spot incident on a semi-infinite medium.

IV. RELATING EXPERIMENTS AND k accum USING THE SUPPRESSION FUNCTION
The suppression function can be used to relate experimental measurements to k accum by mapping length scales to phonon MFPs.For example, k accum can be obtained using Eq. ( 3) with thermoreflectance thermal conductivity measurements and the suppression function from Eq. ( 18) as inputs to the solution of an inverse problem, which was done by Minnich for TTG using convex optimization [21].Alternatively, as we will do here, the experimental measurement can be predicted given k as an input, which can be obtained from models (e.g., Callaway, Born-von Karman-Slack, first-principles, etc.) [8,10,18].This approach is less mathematically complex and allows for a direct comparison to the measurements.
We compare experimental measurements on silicon made by TDTR and BB-FDTR with predicted k exp in Figs.6(a) and 6(b).The solid lines are the predicted accumulation functions from first-principles calculations for silicon plotted as a function of MFP at temperatures of 80 and 300 K [44].Using Eq. ( 3) with the suppression function from Eq. ( 18), we transform this data into a predicted k exp as would be measured by BB-FDTR or TDTR, shown as the dashed lines in Figs.6(a) and 6(b).To make this transformation, we use a spot size of 3r o , which is found by comparing Eq. ( 18) to the suppression function for a Gaussian-shaped spot from Ref. [23] [see Fig. 5(c)] and a temperature-independent value of ε = 1−ρ in Eq. ( 18) that yields the best fit between experimental data and Eq. ( 3), i.e., ε is used as a fitting parameter.
In Fig. 6(a), we transform the k accum vs data from the first-principles calculations into predicted k exp vs 3r o using Eqs.(3) and (18).We find that a value of ε = 1−ρ = 0.88 fits the TDTR measurements from Ref. [23] at temperatures of 80 and 300 K well.Here, we used a heating frequency of 10 6 Hz to determine /L p .We note that the interface between the transducer and substrate for the TDTR data presented is Al/Si.It is reasonable that the value of ε obtained by fitting is related to the properties of this interface.We also plot predicted k exp vs 3r o for a heating frequency of 10 7 Hz with ε = 1−ρ = 0.88 to show how increased TDTR heating frequency is expected to further suppress k exp .
In Fig. 6(b) we transform k accum vs data from the firstprinciples calculations into predicted k exp vs L p using Eqs.( 3) and (18).Here, L p is determined using predicted k exp instead of k bulk to be consistent with our previous presentation of the experimental measurements [16].We find that a value of ε = 1−ρ = 0.6 best describes the BB-FDTR measurements from Ref. [16] at temperatures of 80 and 300 K.In the BB-FDTR results presented, the interface between the transducer and substrate is Cr/Si rather than Al/Si, and it is reasonable that there is a difference in the fitted value of ε for BB-FDTR compared to TDTR.
For silicon at a temperature of 300 K, the predicted k exp vs L p for TDTR in Fig. 6(b) shows L p -dependence over the measurement range although the experimental measurements show no L p -dependence.The TDTR spot size used is the average of the range given in Ref. [6] (3r o = 32.25 μm).For BB-FDTR (3r o = 10.2 μm), the prediction compares well to experimental results at smaller L p .The experimental measurements should plateau at larger L p due to the effect  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 10 6 Hz.Predicted k exp vs 3r o for a heating frequency of 10 7 Hz with ε = 0.88 is shown for comparison.In TDTR, the transducer/substrate interface is Al/Si.(b) k accum from first-principles calculations is transformed into k exp vs L p 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 k exp vs L p with ε = 0.88 is compared to TDTR measurements from Ref. [6] (diamonds). of spot size, but this effect is not observed.More suppression is observed in BB-FDTR relative to TDTR for the available range of TDTR data because a smaller spot size was used and the surface emissivity is lower.At T = 80 K, Eqs.(3) and (18) compare well with BB-FDTR experimental results over all L p .At this temperature, phonons have longer MFPs and are significantly suppressed by the finite spot size, i.e., even for very large L p , k exp will only attain approximately 30% of k bulk due to the spot size restriction.In Figs.6(a) and 6(b) we compare only against the overarching modulation frequency and have neglected the multiple time scales in TDTR that arise from using a pulsed laser.
To generate predicted k exp vs L p in Figs.6(a) and 6(b), we used a suppression function derived from the gray BTE [Eq.(18)] and applied it to the full phonon spectrum, where the MFP is frequency-dependent.A similar approach was used by Collins et al. for TTG, in which a gray suppression function was applied to the full phonon spectrum to obtain predictions of thermal diffusivity as a function of grating period [19].The results were compared to predictions of thermal diffusivity as a function of grating period calculated from the spectral BTE using phonon properties from first-principles calculations.The authors found favorable comparison in that the predicted effective thermal diffusivity varied by less than 7% over grating periods from 10 −1 to 10 6 nm compared to the full spectral models of Si and PbTe at a temperature of 300 K.
The parameters ε and ρ in Eq. ( 18) arise from the analogy with radiative transfer and describe the ability of the surface to emit and reflect phonons.In our comparisons with experimental results we use ε as a fitting parameter but propose that ε is related to the properties of the transducer/substrate interface in BB-FDTR and TDTR experiments.One interpretation is that the phonon emissivity is equal to the transmission coefficient of phonons from the transducer into the substrate [30].
Phonon transmission coefficients are used in the Landauer formulation to make predictions of interface thermal resistance.Following Ref. [33], the total interface resistance R total = 2R T + R L , where R T and R L are the contributions from transverse and longitudinal acoustic phonons, can be derived in a similar manner as thermal conductivity.Beginning with Eqs.(2.10) and (2.11) in Ref. [46] and using a truncated Debye dispersion and Debye density of states, where k B is the Boltzmann constant, v T and v L are the transverse and longitudinal speeds of sound, θ T and θ L are the temperatures associated with the transverse and longitudinal Brillouin zone edge frequencies, y = ω/k B T , and α = ε is the transmission coefficient.Using Eq. ( 21) with values of v T , v L , θ T , and θ L from Ref. [33] and our best-fit values for α = ε, we find that R total = 3.85 m 2 K GW −1 for a Cr/Si interface and R total = 2.63 m 2 K GW −1 for an Al/Si interface at T = 300 K.These values compare well with measured values reported in Refs.[16,15] at T = 300 K (R total = 4.76 m 2 K GW −1 for Cr/Si interface and R total = 2.78 m 2 K GW −1 for Al/Si interface).Furthermore, because ε influences the onset of suppression, we hypothesize that the interface properties contribute to the discrepancy between room temperature BB-FDTR and TDTR heating frequency-dependence for silicon [see Fig. 6(b)], though a spectral phonon model that includes the transducer may be required to reconcile this unresolved question.
It is important to note that previous nondiffusive measurements have been solely attributed to reduced thermal conductivity, i.e., the interface resistance between the substrate and the transducer is assumed constant in diffuse interpretations of the experiments [6,15,16,18].In our formulation, we are comparing a thermal resistance from the BTE that includes a surface temperature drop in the BTE domain to a diffusion solution that does not account for an interface thermal resistance (no surface temperature drop).As a result, we are including the effect of R ε [Eq.(12a)] in our definition of k eff .To generate a suppression function that does not include the surface temperature drop, one can equate the appropriate complex diffusive thermal resistance to R i,x or R i,r , resulting in S i = 1/(1 + iτ ).This result is equivalent for both the planar and spherical geometries and is independent of the particle radius.In Ref. [16], a suppression function was determined from a numerical solution to the 1D, gray BTE for phonons traveling in the positive and negative x directions (μ = 1 or −1).The result is related to S i ; the difference being a factor of π/2 on the x axis, which stems from considering −1 μ 1 when determining S i [47].
To include heating frequency-dependent interface resistance between the transducer and the substrate, a BTE formulation that explicitly includes an interface could be considered and compared to a diffusion solution including an interface.How the transducer affects nondiffusive transport, which is important in interpreting the experiments, has not been explicitly addressed, though it may contribute to the discrepancy between heating frequency-dependent measurements of silicon by BB-FDTR and TDTR.

V. CONCLUSIONS
An analytical suppression function for a system geometrically similar to a thermoreflectance experiment was obtained by solving the BTE for a gray medium.The result accounts for the two dominant length scales in thermoreflectance experiments: thermal penetration depth and heating laser spot radius.We used the suppression function to predict k exp vs L p and k exp vs 3r o to make a direct comparison to experimental measurements by both BB-FDTR and TDTR.Our results corroborate the use of BB-FDTR and TDTR as tools for identifying k accum by generating nondiffusive transport and provide insight and understanding of the measurements.Furthermore, our results suggest that if either L p or r o are much smaller than the phonon MFPs that dominate k, BB-FDTR and TDTR are inadequate for measuring the bulk thermal conductivity of a material.The phonon surface properties ε and ρ affect suppression and may explain discrepancies between TDTR and BB-FDTR measurements of similar samples with different transducers.It is clear that powerful new insight is offered by nondiffusive thermal transport measurements paired with the experiment-specific suppression function to map data into real energy carrier properties.

FIG. 1 .
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.
part of the exponential in Eqs.(10a) and (10b) represents the BTE prediction of penetration depth L p-BTE ,

FIG. 2 .
FIG. 2. (Color online) 1D planar geometry with temporally oscillating 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 τ .
3(a) and 3(b), we plot the magnitude of the thermal resistance of the system from the BTE and diffusion solutions

FIG. 3 .
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 /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].

FIG. 4 .
FIG. 4. (Color online) Spherical particle embedded in an infinite medium with oscillating temperature at the surface of the sphere (r = r o ) with ε = 1−ρ = 1.(a) Magnitude of the thermal resistance from the diffusion and BTE solutions vs /L p and τ .(b) Magnitude of the suppression function plotted as a function of /L p and τ for different values of /r o .For /r o = 0, the results collapse to the 1D planar case shown in Figs.3(a) and 3(b).

FIG. 5 .
FIG. 5. (Color online) Spherical particle embedded in an infinite medium with oscillating temperature at the surface of the sphere (r = r o ) with ε = 1−ρ = 1.(a) Magnitude of the thermal resistance from the diffusive and BTE solutions vs /r o for different values of /L p .(b) Magnitude of the suppression function plotted as a function of /r o for different values of /L p .(c) Comparison of Eq. (18) when /L p = 0 for a particle with radius r o and a particle with radius 3r o , 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 3r o in our analytical result compares well with numerical results from Ref.[23].

FIG. 6 .
FIG.6.(Color online) Comparison of thermal conductivity measurements and predicted k exp for silicon.(a) k accum from first-principles calculations (solid lines) is transformed into predicted k exp vs 3r o (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 10 6 Hz.Predicted k exp vs 3r o for a heating frequency of 10 7 Hz with ε = 0.88 is shown for comparison.In TDTR, the transducer/substrate interface is Al/Si.(b) k accum from first-principles calculations is transformed into k exp vs L p 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 k exp vs L p with ε = 0.88 is compared to TDTR measurements from Ref.[6] (diamonds).