Length Dependent Thermal Conductivity Measurements Yield Phonon Mean Free Path Spectra in Nanostructures

Thermal conductivity measurements over variable lengths on nanostructures such as nanowires provide important information about the mean free paths (MFPs) of the phonons responsible for heat conduction. However, nearly all of these measurements have been interpreted using an average MFP even though phonons in many crystals possess a broad MFP spectrum. Here, we present a reconstruction method to obtain MFP spectra of nanostructures from variable-length thermal conductivity measurements. Using this method, we investigate recently reported length-dependent thermal conductivity measurements on SiGe alloy nanowires and suspended graphene ribbons. We find that the recent measurements on graphene imply that 70% of the heat in graphene is carried by phonons with MFPs longer than 1 micron.

T hermal transport in nanostructures has been a topic of intense interest in recent years 1-3 . When the characteristic dimensions of nanostructures such as the diameter of a nanowire approach phonon mean free paths (MFPs), the thermal conductivity can be substantially smaller than the bulk value due to scattering from sample boundaries. Significant thermal conductivity reductions have been observed in a number of nanoscale systems, including nanowires 4-6 , nanotubes 7 , thin Si membranes 8 , and micron size beams at cryogenic temperatures 9 . This concept has been widely adopted in thermoelectrics applications [10][11][12][13] .
Understanding and engineering the thermal conductivity reduction in nanostructures requires knowledge of phonon scattering mechanisms in the form of the phonon MFPs. The MFP accumulation function, which we term the MFP spectrum in this work, has been demonstrated to be a particularly useful quantity to describe the values of the MFPs relevant for heat conduction 14 . In several works, information about MFPs was obtained by measuring the thermal conductivity over variable lengths of nanostructures such as nanotubes 7 , graphene ribbons 15 and SiGe nanowires 16 . If phonons have MFPs exceeding the distance between the heat source and sink their contribution to thermal conductivity is reduced compared to that in the bulk material, and thus the deviations of the measured thermal conductivity from the bulk value provide information on the phonon MFPs. However, prior studies extracted only an average MFP despite the fact that recent works have demonstrated that in many solids phonon MFPs vary over orders of magnitude, making the approximation of an average MFP for all phonons quite poor 17,18 .
In principle, information about the full MFP spectrum should be contained in these variable-length thermal conductivity measurements, just as the MFP spectrum can be obtained from thermal conductivity measurements performed over variable thermal length scales in MFP spectroscopy 19 . In particular, the method proposed by Minnich based on convex optimization 20 should be applicable to the present situation provided that the suppression function that describes the effect of the finite length on the thermal conductivity can be identified. Li et al obtained the phonon MFP spectrum of graphite along the c-axis from thickness dependent thermal conductivities obtained with molecular dynamics simulations 21 , but their suppression function was not rigorously obtained from the Boltzmann Transport Equation (BTE).
In this report, we present a reconstruction approach to obtain MFP spectra from variable-length thermal conductivity measurements. We use a recently reported analytical solution of the BTE, along with efficient numerical simulations, to identify a suppression function that describes the discrepancy between the actual heat flux and that predicted by Fourier's law for a finite length domain. The MFP spectrum is then obtained using the convex optimization method described in Ref. 20. We apply this approach to SiGe nanowires and graphene ribbons. The measurements on graphene ribbons imply that MFPs are exceedingly long, with 70% of the heat being carried by phonons with MFPs longer than 1 micron.

Theory
Our goal is to relate experimentally measured thermal conductivities to the MFP spectrum, or the accumulated thermal conductivity as a function of MFP 14 . Following the approach of Ref. 20, we therefore seek an equation of the form: where Kn v 5 L v /L is the Knudsen Number, L v is the MFP, L is the sample length along the direction of the temperature gradient, k denotes thermal conductivity as a function of length L, f(L v ) and F(L v ) are differential and accumulative MFP spectra related by ÞdL, and S is the heat flux suppression function that equals the ratio of actual heat flux to the Fourier's law prediction.
The kernel K is defined as K(Kn v ) 5 2dS/dKn v . The inputs to this equation are a finite number of measured thermal conductivities k as a function of lengths L. To close the problem, we must identify the suppression function by solving the BTE, given by 22 : where e v is the desired distribution function, v is the angular frequency, e 0 v is the distribution function at the equilibrium state, v is the group velocity of phonons, and t v is the relaxation time of phonons at certain frequency.
We obtain this function using two distinct approaches: a semianalytic method and a numerical Monte Carlo (MC) method. First, we use a recently reported semi-analytical solution for steady heat conduction through a crystal of thickness L with two blackbody boundaries 23 . In this solution, the BTE is linearized and solved using a series expansion method. The full details are given in Ref. 24. The final result for the suppression function and kernel are: where E n (x) is the exponential integral function, given by: This equation was derived by neglecting temperature slip between the black walls. Physically, this assumption is similar to the weakly quasiballistic regime described in Ref. 23 and implies that the ballistic phonons are low frequency modes with a small heat capacity. The assumption has been shown to be quite accurate for experimentally accessible length scales 24 .
We plot this result in Fig. 1a. For extremely short MFPs compared to the sample length L, the suppression function equals unity, indicating these phonons are diffusive and their heat flux contribution equals the Fourier's law prediction. As the Knudsen number increases, the suppression function decreases and eventually approaches zero, indicating that phonons contribute a smaller amount to the heat flux than predicted by Fourier's law. Physically, this suppression occurs because phonons cannot travel a full MFP before being absorbed by the blackbody boundary.  We additionally solve the BTE numerically to validate the calculations above as well as to consider more complex situations such as when boundary scattering occurs. For this calculation, we use a linearized deviational Monte Carlo method to solve the adjoint BTE as described by Péraud et al 26 . This technique solves BTE by simulating advection and scattering of particles that represent phonons traveling inside the simulation domain. Substantial reductions in computational cost are achieved through a number of simple changes to the original MC algorithm. First, the deviational algorithm simulates only the deviation from a known equilibrium Bose-Einstein distribution, thereby incorporating deterministic information and reducing the variance. Further, for small temperature differences, the collision term in the BTE can be linearized, allowing particles to be simulated completely independently and without spatial and temporal discretization 27 . Next, we use a variable local equilibrium temperature method that closely matches the steady-state temperature profile.
Finally, we solve the adjoint BTE rather than the traditional BTE 28,29 . In the original algorithms of Refs. 26, 27, the probability for a certain phonon mode to be sampled is proportional to the density of states. Therefore, low frequency phonons are rarely sampled even though they contribute substantially to thermal conductivity, leading to large stochastic noise. The adjoint method overcomes this limitation by drawing particles with equal probability among all phonon modes and correcting the bias introduced by this sampling when thermal properties are calculated. With these advances in numerical approach, we are able to solve the BTE in a 100 micron long domain in minutes on a desktop computer. Further, this numerical approach can incorporate boundary scattering mechanisms for arbitrary geometries exactly, unlike the analytical treatment.
To validate the code, we calculate the MFP spectrum for an infinite planar slab with two blackbody boundary conditions, the same problem solved by the semi-analytical method. For this calculation, we use an isotropically averaged dispersion for Si. The original dispersion and relaxation times were calculated by density functional theory (DFT) by Jesús Carrete and N. Mingo with ShengBTE 30,31 and Phonopy 32 , from interatomic force constants obtained with VASP [33][34][35][36] . We reduce computational cost by taking advantage of the cubic symmetry of Si and computing an isotropic equivalent dispersion as described in Ref. 37. Using this dispersion, we calculate the MFP spectrum for variable lengths, as in Fig. 1b. We observe that decreasing the length of the domain results in the suppression of long MFP phonons to thermal conductivity compared to the bulk spectrum. The ratio of the differential MFP spectrum for a finite length to that for an infinite length yields the suppression function and is plotted for two lengths in Fig. 1a, demonstrating that the function obtained from our numerical approach exactly agrees with the analytical result.
With these tools, we now demonstrate the principal result of this work by using the suppression function to reconstruct the MFP spectrum from variable-length thermal conductivity measurements on a Si slab. We synthesized thermal conductivities as a function of length as shown in Fig. 1c. With these length dependent thermal conductivities and the suppression function, we used the same convex optimization method introduced in Ref. 20 to reconstruct the MFP spectrum. As in Figure 1d, the reconstructed accumulative MFP distribution is in excellent agreement with the actual accumulative MFP distribution, which is obtained from the DFT calculation, demonstrating that our approach can accurately reconstruct the MFP spectrum from length-dependent thermal conductivities of Si slabs.
We now numerically demonstrate that our approach can be applied to more general problems than Si slabs with a single scattering mechanism. We consider a SiGe nanowire in which point defects and the nanowire boundaries scatter phonons in addition to the intrinsic phonon-phonon scattering mechanism. The phononphonon relaxation times are taken to be the same as those of pure Si 30-36 , while the mass defect scattering rate is given by ÞAv 438 , where A is a constant of 3.01 3 10 241 s 3 for Si 12x Ge x , which is obtained from Ref. 38. This model predicts a thermal conductivity of 14 W/mK for bulk Si 0.9 Ge 0.1 , which is consistent with other models by DFT calculation and experimental result 38,39 . These scattering rates are combined using Matthiessen's rule. We incorporate boundary scattering by explicitly simulating phonon trajectories inside a nanowire with a square cross-section of size 100 nm by 100 nm. We use Ziman's specularity parameter, p 5 exp(216p 2 s 2 /l 2 ), to determine the probability of specular or diffuse scattering, where s is surface roughness and l is the phonon wavelength 40 .
With this framework, we use our MC simulations to calculate the length-dependent thermal conductivities for a Si 0.9 Ge 0.1 nanowire with surface roughness of s 5 0.1 nm, a 100 nm by 100 nm square cross-section, over lengths from L 5 5 nm to 16 mm as would be obtained in an experiment. Then, using only our knowledge of these thermal conductivities and the suppression function, we perform the reconstruction procedure to obtain the MFP spectrum of the nanowire. This result is shown in Fig. 2a. The reconstructed spectrum is in good agreement with the actual one without requiring any knowledge of the scattering mechanisms in the nanowire. Thus, the success of the reconstruction of MFP spectrum of the nanowire demonstrates the self-consistency of our approach.

Discussion
We now use our approach to examine two recent reports of lengthdependent thermal conductivities in nanostructures. First, we consider SiGe alloy nanowires as investigated by Hsiao et al 16 . These nanowires were reported to have ballistic heat conduction persisting over approximately 8 microns. To investigate this experimental report, we calculate the length-dependent thermal conductivities of Si 0.9 Ge 0.1 nanowires, which has the approximately the same crosssectional area, 100 nm by 100 nm cross section, as that of the nanowires in Ref. 16, whose diameters range from 50 nm to 180 nm. Simulated data are shown in Fig. 2b for nanowires with extremely rough boundaries, s R ', and smooth boundaries, with s 5 0.1 nm. We observe a discrepancy between the trends of experimental data and our simulations. The experimental data suggest that the thermal conductivities of these SiGe nanowires are mainly due to phonons with a narrow MFP spectrum around ,8.3 mm, but our simulations indicate that even for very smooth nanowires (s 5 0.1 nm), phonons within this range only contribute ,15% of the total thermal conductivity (Fig. 2a). Most of the heat is carried by MFPs less than 1 micron, with some heat being carried by longer MFPs in the partially specular case. Additionally, most actual nanowires have surface roughness higher than 0.1 nm 41,42 , and in this case long MFP phonons contribute only a small amount to heat conduction. Due to the relatively large diameter of the nanowire, changes to the phonon dispersion for thermal phonons due to phonon confinement are unlikely 43 . The experimental measurements thus do not agree with our self-consistent calculations. Due to the lack of appropriate experimental reports, we are unable to apply our approach to other nanowire data sets. Further experimental investigation is necessary to address this discrepancy.
Next, we consider recent measurements on graphene 15 . In this work, Xu et al 15 performed thermal conductivity measurements over variable lengths on suspended single-layer graphene ribbons to infer an average MFP of 240 nm at room temperature. Using our approach, we can use these same measurements to obtain the MFP spectrum of graphene. Due to the difficulty of fabrication for very long suspended graphene devices and vulnerability of these devices during measurement, the authors of Ref. 15 didn't obtain the saturated thermal conductivity from an extremely long, or ''bulk'', graphene sample. Therefore, we evaluated the saturation value as ,2,000 W/mK with the same extrapolating method in Ref. 21, which is in the range of previous reported experimental results [44][45][46] .
There are two subtleties that require discussion before applying our reconstruction approach to graphene. First, our derivation is based on the relaxation time approximation (RTA) of the BTE. It is well known that the RTA with the computed relaxation times from DFT underpredicts the thermal conductivity of graphene due to the importance of normal processes 47 . However, while our derivation is based on the RTA, we do not make any assumption of the values of the relaxation times, but rather only that an effective relaxation time for each phonon mode can be identified. Our approach can be applied to graphene provided that we regard the MFP variable as an effective MFP that represents the average propagation length for a particular phonon frequency as determined by both normal and Umklapp processes.
Second, we have derived our suppression function for an isotropic material with a three-dimensional phase space. While graphene can be reasonably modeled as isotropic for the in-plane directions, the phase space is two-dimensional. This dimensionality change requires a modification of the form of the BTE for 2D materials. Repeating the derivation in Ref. 24 for a 2D phase space yields the suppression function S and kernel K as:  Figure 3a plots the 2D and 3D suppression functions and kernels, demonstrating that the two are quite similar. Using the 2D kernel, we apply our method to obtain the MFP spectrum of a graphene ribbon as in Fig. 3b. This result shows that MFPs in graphene span a large range from ,100 nm to ,10 mm. In addition, we observe that a large portion of heat in graphene is carried by long MFP phonons: ,70% of thermal conductivity are from phonons with MFPs greater than 1 mm. Further, as reported in Ref. 15, the widths of these graphene ribbons are only between 2 to 4 microns. From Fig. 3b, phonons with MFPs longer than 4 mm still carry 13% of heat, which means that some of these long MFP phonons must be specularly reflected at the edges of the graphene ribbons. This observation further confirms the report in Ref. 15 of weak width-dependent thermal conductivities of suspended graphene ribbons when widths are larger than 1.5 mm. Our MFP reconstruction approach has thus provided valuable insights into the intrinsic and edge scattering mechanisms in graphene ribbons that are difficult to obtain from knowledge of only the average MFP.

Summary
We have presented a reconstruction method that allows MFP spectra of nanostructures to be obtained from length-dependent thermal conductivity measurements. Our approach requires no prior knowledge of the scattering mechanisms in the nanostructure. By applying our approach to recent measurements on graphene ribbons, we find that more than half of the heat in graphene is contributed by phonons with MFPs exceeding 1 micron.