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 transport10,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

vykfky+vzkfkz=fkf0(T(z,t))τk,
(1)

where fk is the desired distribution function, vyk and vzk 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 gk=ff0(T(z,t)), yielding

Λykgky+gk=Λzkf0TTz.
(2)

This result assumes that gk/z0 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

g+=S0(z)(eη1),
(3)
g=S0(z)(eξyη1),
(4)

where η=y/Λyk,ξy=d/Λyk is the nondimensional thickness defined by the y component of the MFP, and S0(z)=Λzkf0TTz. The thermal conductivity can be obtained by substituting this solution into the expression for heat flux

Jq=0d1Vkωgkvzkdy,
(5)

where V is the crystal volume. After evaluating this expression, the thermal conductivity of the thin film along the z direction is identified as

κzz(d)=ky>0κzz,k(11eξyξy)=ky>0κzz,iS(ξy),
(6)

where the intrinsic thermal conductivity κzz,k=Ckvzk2τ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 is11 

S(ξ)=138ξ[ 14E3(ξ)+4E5(ξ) ],
(7)

where En(ξ) is the exponential integral function and ξ=d/Λω where Λω 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 κ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 form23 

τp(ω)1=Ap(ωωmax,p)2,

where τp(ω) is the relaxation time for branch p, Ap 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 Ap 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×1016 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×Ap,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 μ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 experiment21 and calculated with molecular dynamics.22 

FIG. 1.

Normalized, accumulated thermal conductivity versus component of MFP along the transport direction for each crystal axis using the actual graphite phonon dispersion and semi-empirical relaxation times. The c-axis MFPs are on the order of hundreds of nanometers, in agreement with experiment and considerably shorter than those along the ab-axis.21 

FIG. 1.

Normalized, accumulated thermal conductivity versus component of MFP along the transport direction for each crystal axis using the actual graphite phonon dispersion and semi-empirical relaxation times. The c-axis MFPs are on the order of hundreds of nanometers, in agreement with experiment and considerably shorter than those along the ab-axis.21 

Close modal

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.

FIG. 2.

(a) Actual ab-axis thermal conductivity using Eq. (6) (solid line) versus film thickness and comparison using the same equation with the ab-axis rather than c-axis component of the MFPs in the boundary scattering function (dashed line). The actual thermal conductivity remains high even for thicknesses as small as 10 nm due to the small MFP component along the thickness direction of the film. (b) Corresponding plot for c-axis thermal conductivity. Counterintuitively, the correct use of the ab-axis MFPs in Eq. (6) leads to a higher thermal conductivity than using the c-axis MFPs, even though the ab-axis MFPs are larger than those along the c-axis.

FIG. 2.

(a) Actual ab-axis thermal conductivity using Eq. (6) (solid line) versus film thickness and comparison using the same equation with the ab-axis rather than c-axis component of the MFPs in the boundary scattering function (dashed line). The actual thermal conductivity remains high even for thicknesses as small as 10 nm due to the small MFP component along the thickness direction of the film. (b) Corresponding plot for c-axis thermal conductivity. Counterintuitively, the correct use of the ab-axis MFPs in Eq. (6) leads to a higher thermal conductivity than using the c-axis MFPs, even though the ab-axis MFPs are larger than those along the c-axis.

Close modal

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 κc,k=Ckvck2τkS(ξ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 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.

FIG. 3.

Spectral c-axis thermal conductivity versus phonon frequency using Eq. (6) (crosses) and with the c-axis MFPs in the boundary scattering function (circles). Even though the ab-axis MFPs are longer than those in the c-axis, the correct boundary scattering function only minimally affects phonons with long c-axis MFPs, resulting in a higher thermal conductivity than if the shorter c-axis MFPs are used in the boundary scattering function.

FIG. 3.

Spectral c-axis thermal conductivity versus phonon frequency using Eq. (6) (crosses) and with the c-axis MFPs in the boundary scattering function (circles). Even though the ab-axis MFPs are longer than those in the c-axis, the correct boundary scattering function only minimally affects phonons with long c-axis MFPs, resulting in a higher thermal conductivity than if the shorter c-axis MFPs are used in the boundary scattering function.

Close modal

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.

1.
D. G.
Cahill
,
P. V.
Braun
,
G.
Chen
,
D. R.
Clarke
,
S.
Fan
,
K. E.
Goodson
,
P.
Keblinski
,
W. P.
King
,
G. D.
Mahan
,
A.
Majumdar
,
H. J.
Maris
,
S. R.
Phillpot
,
E.
Pop
, and
L.
Shi
, “
Nanoscale thermal transport. ii. 2003–2012
,”
Appl. Phys. Rev.
1
(
1
),
011305
(
2014
).
2.
Y.
Wang
,
H.
Huang
, and
X.
Ruan
, “
Decomposition of coherent and incoherent phonon conduction in superlattices and random multilayers
,”
Phys. Rev. B
90
,
165406
(
2014
).
3.
X.
Wang
and
B.
Huang
, “
Computational study of in-plane phonon transport in si thin films
,”
Sci. Rep.
4
,
6399
(
2014
).
4.
A. A.
Balandin
and
D. L.
Nika
, “
Phononics in low-dimensional materials
,”
Mater. Today
15
(
6
),
266
275
(
2012
).
5.
I.
Chowdhury
,
R.
Prasher
,
K.
Lofgreen
,
G.
Chrysler
,
S.
Narasimhan
,
R.
Mahajan
,
D.
Koester
,
R.
Alley
, and
R.
Venkatasubramanian
, “
On-chip cooling by superlattice-based thin-film thermoelectrics
,”
Nat. Nanotechnol.
4
(
4
),
235
238
(
2009
).
6.
Z.
Su
,
L.
Huang
,
F.
Liu
,
J. P.
Freedman
,
L. M.
Porter
,
R. F.
Davis
, and
J. A.
Malen
, “
Layer-by-layer thermal conductivities of the group III nitride films in blue/green light emitting diodes
,”
Appl. Phys. Lett.
100
(
20
),
201106
(
2012
).
7.
A. L.
Moore
and
L.
Shi
, “
Emerging challenges and materials for thermal management of electronics
,”
Mater. Today
17
(
4
),
163
174
(
2014
).
8.
Z.
Yan
,
G.
Liu
,
J. M.
Khan
, and
A. A.
Balandin
, “
Graphene quilts for thermal management of high power GaN transistors
,”
Nat. Commun.
3
,
827
(
2012
).
9.
J.
Cho
,
Y.
Li
,
W. E.
Hoke
,
D. H.
Altman
,
M.
Asheghi
, and
K. E.
Goodson
, “
Phonon scattering in strained transition layers for GaN heteroepitaxy
,”
Phys. Rev. B
89
,
115301
(
2014
).
10.
F.
Fuchs
, “
The conductivity of thin metallic films according to the electron theory of metals
,”
Proc. Cambridge Philos. Soc.
34
,
100
108
(
1938
).
11.
E. H.
Sondheimer
, “
The mean free path of electrons in metals
,”
Adv. Phys.
1
,
1
42
(
1952
).
12.
G.
Chen
and
C. L.
Tien
, “
Thermal conductivities of quantum well structures
,”
J. Thermophys. Heat Transfer
7
,
311
318
(
1993
).
13.
G.
Chen
, “
Size and interface effects on thermal conductivity of superlattices and periodic thin-film structures
,”
J. Heat Transfer
119
,
220
229
(
1997
).
14.
S.
Mazumder
and
A.
Majumdar
, “
Monte Carlo study of phonon transport in solid thin films including dispersion and polarization
,”
J. Heat Transfer
123
,
749
759
(
2001
).
15.
W.
Jang
,
Z.
Chen
,
W.
Bao
,
C. N.
Lau
, and
C.
Dames
, “
Thickness-dependent thermal conductivity of encased graphene and ultrathin graphite
,”
Nano Lett.
10
(
10
),
3909
3913
(
2010
).
16.
M. T.
Pettes
,
H.
Ji
,
R. S.
Ruoff
, and
L.
Shi
, “
Thermal transport in three-dimensional foam architectures of few-layer graphene and ultrathin graphite
,”
Nano Lett.
12
(
6
),
2959
2964
(
2012
).
17.
G.
Chen
,
Nanoscale Energy Transport and Conversion
(
Oxford University Press
,
New York
,
2005
).
18.
L.
Lindsay
and
D. A.
Broido
, “
Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and grapheme
,”
Phys. Rev. B
81
(
20
),
205441
(
2010
).
19.
A. J.
Minnich
, “
Phonon heat conduction in layered anisotropic crystals
,”
Phys. Rev. B
91
(
8
),
085206
(
2015
).
20.
Y. S.
Touloukian
,
Thermophysical Properties of Matter: Thermal Conductivity: Nonmetallic Solids
(
IFI/Plenum
,
1970
), ISBN: 9780306670206.
21.
Q.
Fu
,
J.
Yang
,
Y.
Chen
,
D.
Li
, and
D.
Xu
, “
Experimental evidence of very long intrinsic phonon mean free path along the c-axis of graphite
,”
Appl. Phys. Lett.
106
(
3
),
031905
(
2015
).
22.
Z.
Wei
,
J.
Yang
,
W.
Chen
,
K.
Bi
,
D.
Li
, and
Y.
Chen
, “
Phonon mean free path of graphite along the c-axis
,”
Appl. Phys. Lett.
104
(
8
),
081903
(
2014
).
23.
M. T.
Pettes
,
M. M.
Sadeghi
,
H.
Ji
,
I.
Jo
,
W.
Wu
,
R. S.
Ruoff
, and
L.
Shi
, “
Scattering of phonons by high-concentration isotopic impurities in ultrathin graphite
,”
Phys. Rev. B
91
(
3
),
035429
(
2015
).
24.
T.
Tanaka
and
H.
Suzuki
, “
The thermal diffusivity of pyrolytic graphite at high temperatures
,”
Carbon
10
(
3
),
253
257
(
1972
).