Thermal phonon boundary scattering in anisotropic thin films

Boundary scattering of thermal phonons in thin solid films is typically analyzed using Fuchs-Sondheimer theory, which provides a simple equation to calculate the reduction of thermal conductivity as a function of the film thickness. However, this widely-used equation is not applicable to highly anisotropic solids like graphite because it assumes the phonon dispersion is isotropic. Here, we derive a generalization of the Fuchs-Sondheimer equation for solids with arbitrary dispersion relations and examine its predictions for graphite. We find that the isotropic equation vastly overestimates the boundary scattering that occurs in thin graphite films due to the highly anisotropic group velocity, and that graphite can maintain its high in-plane thermal conductivity even in thin films with thicknesses as small as ten nanometers.

Thermal transport in thin solid films with thicknesses from tens of nanometers to micrometers is a topic of considerable importance, 1-4 with such films being used in applications ranging from quantum well lasers to LED lighting. [5][6][7] In these thin films, the phonon mean free paths (MFPs) are comparable to the film thickness, resulting in boundary scattering that reduces the thermal conductivity and inhibits heat dissipation. This thermal management problem is presently an important challenge in many devices such as GaN transistors. 6,8,9 Heat transport along the in-plane direction of thin films is typically described using the Fuchs-Sondheimer equation, which is an analytical solution of the Boltzmann transport equation (BTE) for thin films with partially specular and partially diffuse walls. This equation was originally derived for electron transport 10,11 and later extended to phonon thermal transport assuming an average phonon MFP, enabling the calculation of thermal conductivity as a function of the film thickness. 12, 13 Mazumder and Majumdar used a Monte-Carlo method to study the phonon transport along a silicon thin film including phonon dispersion and polarizations. 14 While the Fuchs-Sondheimer equation is in wide use, it cannot be applied to anisotropic transport in its typical form because it assumes the crystal of interest is isotropic. However, there are many situations in which boundary scattering occurs in thin anisotropic films, with the most familiar example being thin graphite films. Such films have been studied experimentally, 15 and few-layer graphene films have been investigated as in heat spreaders for GaN transistors. 8 Provided that the films are sufficiently thick that phonon dispersion modifications in the cross-plane direction can be neglected, mathematically describing thermal transport in anisotropic thin films requires a Fuchs-Sondheimer equation that is valid for any crystal, regardless of its anisotropy. Surprisingly, despite the simplicity of the derivation, no such equation has been reported to the best of our knowledge.
Here, we report the generalization of Fuchs-Sondheimer theory to crystals with arbitrary anisotropies. We find that highly anisotropic solids with small group velocities along certain crystallographic directions experience minimal boundary scattering because the thermal conductivity reduction depends only on the component of the MFP normal to the boundary rather than the overall MFP. As a result, thin films of anisotropic crystals like graphite maintain their high thermal conductivity even as the film thickness becomes very small. This observation has important implications for heat spreading in electronic devices and the thermal conductivity of graphite foams. 16 We begin by considering the steady BTE in a thin film under the relaxation time approximation. We let the y direction represent the cross-plane direction and assume that a thermal gradient exists along the z direction. Then, the BTE is given by: where f k is the desired distribution function, v yk and v zk are the components of the phonon group velocity along the y and z directions, τ k is the phonon relaxation time, and k is the phonon wavevector in phase space. This equation is solved for the thin film by letting the small perturbation g k = f − f 0 (T (z, t)), yielding: This result assumes that ∂g k /∂z ≈ 0 and that the length scale of the thermal gradient is long compared to the phonon MFPs so that a temperature gradient can be defined.
Equation 2 is a one-dimensional ODE with the boundary conditions g k (k y < 0, y = d) = g k (k y > 0, y = 0) = 0 corresponding to thermalizing, diffuse scattering, meaning that the phonon distribution emerging from the wall is a thermal distribution at the local boundary temperature. The general solution of this equation is well-known and given in Ref. 17 as: where η = y/Λ yk , ξ y = d/Λ yk is the nondimensional thickness defined by the y component of the MFP, and S 0 (z) = −Λ zk ∂f 0 ∂T ∂T ∂z . The thermal conductivity can be obtained by substituting this solution into the expression for heat flux: where V is the crystal volume. After evaluating this expression, the thermal conductivity of the thin film along the z direction is identified as: where the intrinsic thermal conductivity κ zz,k = C k v 2 zk τ k . S(ξ y ), which we term the boundary scattering function, is the generalized version of the Fuchs-Sondheimer expression describing the reduction in thermal conductivity that occurs due to boundary scattering.
For reference, the traditional Fuchs-Sondheimer expression obtained by integrating over all solid angles in a spherically symmetric Brillouin zone is 11 : where E n (ξ) is the exponential integral function and ξ = d/Λ ω where Λ ω is the isotropic MFP that depends only on phonon frequency. To examine this prediction, we consider transport along the basal plane of a thin film of graphite. We use a phonon dispersion calculated from an optimized Tersoff potential provided by Lucas Lindsay. 18 The actual phonon-phonon relaxation times of graphite are not readily available, and thus we must assume a general form to proceed. In a previous work for which we modeled the graphite dispersion with an anisotropic Debye model, 19 we found that isotropic relaxation times were able to explain the thermal conductivity anisotropy of graphite. In this work for which we use the actual dispersion, we find that these same relaxation times yield κ yy = 2000 W/mK and κ zz = 230 W/mK. Given that the in-and cross-plane thermal conductivities of graphite are 2000 W/mK and 6 W/mK, 20 respectively, this calculation indicates that these relaxation times cannot explain the large thermal anisotropy of graphite.
We now must find the relaxation times that best reproduce the following experimental observations. First, the thermal conductivity along each crystal axis is well-known. 20 Second, recent computational and experimental reports indicate that phonon MFPs along the c-axis are on the order of several hundred nanometers. 21,22 We satisfy these constraints by assuming relaxation times of the form 23 : where τ p (ω) is the relaxation time for branch p, A p is a constant for each branch, ω is the phonon frequency, and ω max,p is the maximum phonon frequency for each branch.
Because identifying phonon polarizations off high-symmetry points in the Brillouin zone is challenging, branches are determined by sorting the frequencies obtained from the dynamical matrix from low to high.
We can satisfy the constraints by changing A p for specific branches even though the relaxation times remain isotropic because the different branches contribute differently to thermal conductivity along each crystal axis, We take A p = A p,0 = 8 × 10 16 s −1 for all branches except 1, 2, 3, and 5. Branches 1 and 2 primarily carry heat along the c-axis and must have smaller relaxation times to yield the correct c-axis thermal conductivity. We find that A p,1 = A p,0 /100 and A p,2 = A p,0 /20. Branch 3 must have shorter relaxation times, yielding A p,3 = A p,0 /20. We also cap the maximum relaxation times at 250 ps; otherwise the c-axis MFPs are too long compared to experiment. Branch 5, the LA branch along the ab axis, must have longer relaxation times to reproduce the high ab-axis thermal conductivity, A p,3 = 3 × A p,0 , with a maximum relaxation time of 500 ps. Finally, the minimum relaxation time of all branches is enforced to be 5 ps so that the MFPs are not unphysically short. These relaxation times yield κ ab = 1500 W/mK and κ c = 6.9 W/mK, in reasonable agreement with experiment. We emphasize that these relaxation times are not necessarily those of actual graphite, but the thermal transport properties obtained using them are sufficiently close to the experimental values that they can be used for further analysis.
The accumulative thermal conductivity versus the component of the MFP along the a, b, and c crystal axes is shown in Fig. 1, showing that the chosen relaxation times result in phonons along the ab-axis having MFPs on the order of microns to tens of microns, as might be expected due to the high basal plane thermal conductivity. On the other hand, phonons along the c-axis have MFPs up to 1 micron and typically in the hundreds of nanometers range, much shorter than those along the ab axis but still considerably longer than prior estimates. 24 The calculated MFPs are, however, of the same order as those recently observed in experiment 21 and calculated with molecular dynamics. 22 The large difference in MFP between the ab and c axes affects thermal transport in thin films because the strength of boundary scattering only depends on the component of MFP normal to the boundary. In the case of transport in a thin film along the ab-axis, the far shorter c-axis MFP leads to only minimal boundary scattering even for film thicknesses as small as 10 nm.
The ab-axis thermal conductivity versus film thickness obtained using Eq. 6 is shown in Fig. 2a. The ab-axis thermal conductivity remains close to 1000 W/mK even at thicknesses of around 10 nm, and nearly reaches its ultimate value at a thickness of around 1 micron.
These observations are quite unexpected if one considers that the MFPs along the ab-axis are tens of microns as in Fig. 1.
To place these calculations in perspective, we seek to calculate the thermal conductivity of an equivalent isotropic crystal with the same ab-axis MFPs. However, we were unable to find a simple way to convert the highly anisotropic graphite dispersion to an equivalent As a result, even extremely thin graphite films maintain their high thermal conductivity despite the long MFPs along the ab-axis.
The same calculation as above but for a thermal gradient along the c-axis is shown in Fig.   2b. In this case, the graphite is infinite along the c-axis and has a finite thickness along one basal axis, similar to the geometry used in a prior experimental study. 21 The thermal gradient exists along the c-axis of the graphite. The thermal conductivity versus film thickness is presented according to Eq. 6, which uses the ab-axis MFPs, and for comparison the same calculation with the c-axis MFPs in the boundary scattering function. The calculation shows that, counterintuitively, the c-axis thermal conductivity obtained with the correct ab-axis MFPs in the boundary scattering function is actually higher than that with the c-axis MFPs, even though the ab-axis MFPs are considerably longer than those in the c-axis. The reason is that modes that primarily contribute to c-axis thermal conductivity have short in-plane MFPs and thus are minimally affected by the boundary function in Eq. 6. In contrast, using Eq. 6 with the c-axis MFPs incorrectly reduces the contribution of the long MFP c-axis modes that contribute the most to c-axis thermal conductivity, resulting in a smaller thermal conductivity than predicted by Eq. 6.
This point is illustrated in Fig. 3, which shows the spectral c-axis thermal conductivity κ c,k = C k v 2 ck τ k S(ξ y ) versus phonon frequency calculated using the ab-and c-axis components of the MFPs in the boundary scattering function S(ξ y ). The correct boundary scattering function based on the ab-axis MFPs only minimally affects the phonons that primarily contribute to c-axis heat conduction because these phonons have short ab-axis MFPs and long c-axis MFPs. The primary effect of the correct boundary scattering function is to make the