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 10 nm.

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 lighting with light-emitting diodes lighting.^{5–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.

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 and assuming that the film is sufficiently thick that the phonon dispersion and relaxation times of the thin film are accurately described by those of the bulk solid. 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 $fk$ is the desired distribution function, $vyk$ and $vzk$ are the components of the phonon group velocity along the *y* and *z* directions, $\tau 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 $gk=f\u2212f0(T(z,t))$, yielding

This result assumes that $\u2202gk/\u2202z\u22480$ 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 ordinary differential equation with the boundary conditions $gk(ky<0,y=d)=gk(ky>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 $\eta =y/\Lambda yk,\u2009\xi y=d/\Lambda yk$ is the nondimensional thickness defined by the *y* component of the MFP, and $S0(z)=\u2212\Lambda zk\u2202f0\u2202T\u2202T\u2202z$. 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 $\kappa zz,k=Ckvzk2\tau k$. $S(\xi 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 $En(\xi )$ is the exponential integral function and $\xi =d/\Lambda \omega $ where $\Lambda \omega $ is the isotropic MFP that depends only on phonon frequency.

Comparing Eqs. (6) and (7), we see that the traditional Fuchs-Sondheimer expression averages out the solid angle dependence in phase space so that the final expression depends only on the overall MFP. The generalized expression does not make any assumptions about the symmetry of the crystal, and so the expression depends on the components of the MFPs along the crystal axis normal to the boundary. If this MFP component is different than the value along the in-plane axis, the actual in-plane thermal conductivity obtained with Eq. (6) may be very different than the prediction of the traditional expression.

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 Lindsay *et al.*^{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 the 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 $\kappa yy=2000$ W/mK and $\kappa 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 $\tau p(\omega )$ is the relaxation time for branch *p*, *A _{p}* is a constant for each branch, $\omega $ is the phonon frequency, and $\omega 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 $Ap=Ap,0=8\xd71016$ 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 $Ap,1=Ap,0/100$ and $Ap,2=Ap,0/20$. Branch 3 must have shorter relaxation times, yielding $Ap,3=Ap,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 longitudinal acoustic branch along the ab axis, must have longer relaxation times to reproduce the high ab-axis thermal conductivity, $Ap,3=3\xd7Ap,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 $\kappa ab=1500$ W/mK and $\kappa 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 *μ*m 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. In the following discussion, we consider only film thicknesses greater than 10 nm so that possible modifications of the phonon dispersion in ultrathin films can be neglected.

The ab-axis thermal conductivity versus film thickness obtained using Eq. (6) is shown in Fig. 2(a). 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 *μ*m. 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 isotropic one. As a means of comparison, we instead calculate the thermal conductivity reduction using the longer ab-axis MFPs in the boundary scattering function rather than c-axis MFPs in the same figure. The result shows the trend expected of an isotropic solid, with the thermal conductivity approaching zero for very thin films and slowly approaching the bulk value as the film thickness increases. This prediction differs greatly from the correct calculation because of the form of Eq. (6): the thermal conductivity reduction depends only on the component of the MFP along the thickness direction, regardless of the ab-axis MFPs. 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. 2(b). In this case, the graphite is of infinite extent 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 $\kappa c,k=Ckvck2\tau kS(\xi y)$ versus phonon frequency calculated using the ab- and c-axis components of the MFPs in the boundary scattering function $S(\xi y)$. The correct boundary scattering function based on the ab-axis MFPs only minimally affects the phonons that primarily contribute to the 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 MFPs of phonons with small c-axis MFPs even smaller, thereby resulting in only a small reduction to the thermal conductivity. On the other hand, the incorrect boundary scattering function significantly reduces the MFPs of phonons with long c-axis MFPs, resulting in an anomalously large reduction in the thermal conductivity. These calculations show the importance of considering the dispersion anisotropies of a crystal when considering boundary scattering in thin films.

In summary, we have derived a generalized Fuchs-Sondheimer theory for crystals with arbitrary anisotropies and applied it to thin graphite films. We find that the large velocity anisotropy in graphite causes boundary scattering to have only a minimal effect on the ab-axis thermal conductivity in films as thin as 10 nm. This observation has important implications for heat-spreading in electronic devices using thin graphite films and heat conduction in graphite foams.

The author thanks Lucas Lindsay for providing the graphite dispersion and for useful discussions. This work was sponsored in part by the National Science Foundation under Grant No. CBET CAREER 1254213 and by Boeing under the Boeing-Caltech Strategic Research & Development Relationship Agreement.