Intrinsic lattice thermal conductivity of semiconductors
from first principles
D. A. Broido
a兲
Department of Physics, Boston College, Chestnut Hill, Massachusetts 02467, USA
M. Malorny and G. Birner
Department of Physics, University of Regensburg, D-93040 Regensburg, Germany
Natalio Mingo
CEA-Grenoble, 17 Rue des Martyrs, Grenoble 38054, France
and Department of Electrical Engineering, University of California, Santa Cruz, California 95064, USA
D. A. Stewart
Cornell Nanoscale Facility, Cornell University, Ithaca, New York, 14853, USA
共Received 30 October 2007; accepted 19 November 2007; published online 7 December 2007兲
We present an ab initio theoretical approach to accurately describe phonon thermal transport in
semiconductors and insulators free of adjustable parameters. This technique combines a Boltzmann
formalism with density functional calculations of harmonic and anharmonic interatomic force
constants. Without any fitting parameters, we obtain excellent agreement 共⬍5% difference at room
temperature兲 between the calculated and measured intrinsic lattice thermal conductivities of silicon
and germanium. As such, this method may provide predictive theoretical guidance to experimental
thermal transport studies of bulk and nanomaterials as well as facilitating the design of new
materials. © 2007 American Institute of Physics. 关DOI: 10.1063/1.2822891兴
Thermal conductivity is a fundamental transport param-
eter that is commonly used to characterize a broad range of
materials and systems. In semiconductors and insulators,
heat is carried by the vibrating lattice. A predictive theoreti-
cal approach to calculate the lattice thermal conductivity in
these materials is of tremendous importance for modern sci-
ence and technology. It would facilitate understanding of
heat dissipation in microelectronics and nanoelectronics
1
as
well as assist in material design for efficient thermoelectric
refrigeration and power generation.
2
Above a few tens of Kelvin the behavior of the lattice
thermal conductivity of semiconductors is usually dominated
by phonon-phonon scattering, which arises because of the
anharmonicity of the interatomic potential. Unlike phonon
scattering by defects, impurities or boundaries, phonon-
phonon scattering is an intrinsic resistive process. A micro-
scopic description of the intrinsic lattice thermal conductivity
共i兲
of semiconductors and insulators was first formulated
theoretically by Peierls
3
in 1929 through what has become
known as the phonon Boltzmann equation 共PBE兲.
While the general framework to describe thermal con-
ductivity is well known, the development of a predictive the-
oretical approach to calculate
共i兲
has been hindered by the
significant complexity inherent in describing 共i兲 interatomic
forces between atoms, and 共ii兲 the inelastic phonon-phonon
scattering processes. This is true even for common materials
such as silicon and germanium.
1
A tremendous simplification
is achieved in the calculation of
共i兲
in bulk
semiconductors
4,5
and nanomaterials
6
by using the relaxation
time approximation to solve the PBE. However, the relax-
ation time approximation is derived under the assumption of
elastic scattering, but the anharmonic phonon-phonon scat-
tering is an inelastic scattering process.
4
Furthermore, be-
cause the relaxation time approximation introduces adjust-
able parameters that are fit to existing experimental data, it
has limited predictive power. In recent years, molecular dy-
namics approaches to calculate the lattice thermal conductiv-
ity of materials
7,8
have been introduced. These approaches
have the advantage that they treat anharmonicity to all or-
ders. However, they typically rely on empirical interatomic
potentials 共EIPs兲 that are fit to experimental properties of
materials 共i.e., crystal structure, elastic constants, point de-
fects, etc.兲. Also, the atomic motion is treated classically so
molecular dynamics approaches are most appropriate for
temperatures typically much higher than 300 K.
In this letter, we present a predictive theoretical ap-
proach to calculate
共i兲
that invokes no adjustable param-
eters and is valid over a wide temperature range around room
temperature. This approach implements an exact iterative so-
lution of the PBE for phonon transport,
9,10
which explicitly
incorporates the quantum mechanical phonon-phonon scat-
tering processes and solves self-consistently for the phonon
distribution function. The only inputs required for the exact
solution of the PBE are the harmonic and anharmonic inter-
atomic force constants 共IFCs兲, and these are determined from
first principles using density functional perturbation
theory.
11,12
A key initial test of the predictive capability of
this approach is how the calculated
共i兲
’s compare to the
corresponding measured lattice thermal conductivity in com-
monly studied materials such as silicon and germanium. We
find that our calculated
共i兲
’s show excellent agreement with
the experimentally determined
共i兲
for isotopically enriched
Si
13
and Ge.
14
We begin by considering a perfect bulk crystal free of
defects or impurities and take silicon or germanium atoms to
reside on a diamond lattice. The lowest order scattering pro-
cesses are between three phonons.
4,15
These processes
are constrained to satisfy conservation of energy,
j
共q兲±
j
⬘
共q
⬘
兲=
j
⬙
共q
⬙
兲, and quasimomentum,
a兲
Author to whom correspondence should be addressed. Electronic mail:
broido@bc.edu.
APPLIED PHYSICS LETTERS 91, 231922 共2007兲
0003-6951/2007/91共23兲/231922/3/$23.00 © 2007 American Institute of Physics91, 231922-1
Downloaded 10 Dec 2007 to 128.253.198.179. Redistribution subject to AIP license or copyright; see http://apl.aip.org/apl/copyright.jsp
q±q
⬘
=q
⬙
+K, where q, j, and
j
共q兲 are the phonon momen-
tum, branch index, and frequency, and K is a reciprocal lat-
tice vector that is zero for normal processes and nonzero for
umklapp processes. A small temperature gradient ⵜT is taken
to perturb the phonon distribution function n
=n
0
+n
1
,
where is a short hand for 共q , j兲, n
0
⬅n
0
共
兲 is the equi-
librium 共Bose兲 phonon distribution function, and the non-
equilibrium part n
1
produces the thermal current. The PBE
is
4,9,10
v
· ⵜT
n
0
T
=
兺
⬘
⬙
冋
W
⬘
⬙
+
共⌿
⬙
− ⌿
⬘
− ⌿
兲
+
1
2
W
⬘
⬙
−
共⌿
⬙
+ ⌿
⬘
− ⌿
兲
册
, 共1兲
where the left hand side of Eq. 共1兲 represents phonon diffu-
sion induced by the thermal gradient and the right hand side
represents the collision term for three phonon interactions. In
Eq. 共1兲, v
is the phonon velocity in mode, , W
⬘
⬙
±
are the
three-phonon scattering rates, ⌿
=n
1
/ 关n
0
共n
0
+1兲兴, and the
sums on the right-hand side of Eq. 共1兲 are over the phase
space of energy and momentum conserving normal and um-
klapp three-phonon processes. The presence of ⌿
⬘
and ⌿
⬙
conveys the inelastic nature of the phonon-phonon scatter-
ing. Replacement of collision integral by −n
1
/
gives the
commonly used relaxation time approximation.
4,5
For
phonon-phonon umklapp scattering, the relaxation time is
typically taken
4–6
to depend on frequency and temperature as
⬃
−2
T
−3
although this form was derived assuming low
frequency and low temperature
16
where umklapp scattering
is typically weak.
The three-phonon scattering rates W
⬘
⬙
±
determined
from Fermi’s golden rule are
9,10
W
⬘
⬙
±
=
ប
4N
共n
0
+1兲共n
0
⬘
+1/2±1/2兲n
0
⬙
⬘
⬙
兩V
±
共,
⬘
,
⬙
兲兩
2
⫻
␦
共
±
⬘
−
⬙
兲. 共2兲
Here, N is the number of unit cells and the delta function
ensures energy conservation. The phonon frequencies 兵
其
and eigenmodes are determined by diagonalizing the dy-
namical matrix, which depends on the harmonic IFCs,
⌽
␣
共0
;ᐉ
⬘
⬘
兲, that connect pairs of atoms 0
and ᐉ
⬘
⬘
. The
notation ᐉ
specifies the
th atom in the ᐉth unit cell, and
␣
and

are Cartesian components. The three-phonon scatter-
ing matrix elements, V
±
共 ,
⬘
,
⬙
兲=V共j ,−q , j
⬘
; ⫿q
⬘
; j
⬙
,q
⬙
兲
measure the strength of the scattering events and are given
by
V共j,q , j
⬘
;q
⬘
; j
⬙
,q
⬙
兲
=
兺
兺
ᐉ
⬘
⬘
兺
ᐉ
⬙
⬙
兺
␣␥
⌽
␣␥
共0
,ᐉ
⬘
⬘
,ᐉ
⬙
⬙
兲e
iq
⬘
·R
ᐉ
⬘
e
iq
⬙
·R
ᐉ
⬙
⫻
e
␣
j
共q兲e

⬘
j
⬘
共q
⬘
兲e
␥
⬙
j
⬙
共q
⬙
兲
冑
M
M
⬘
M
⬙
共3兲
where R
ᐉ
is a lattice vector and M
is the mass of the
th
atom. The ⌽
␣␥
共0
,ᐉ
⬘
⬘
,ᐉ
⬙
⬙
兲 are third order anharmonic
IFCs for the indicated triplets of atoms, and the e’s are pho-
non eigenvectors.
Equation 共1兲 is solved using an iterative approach
9,10
to
obtain the nonequilibrium distribution functions, ⌿
. These
are used to calculate the phonon thermal conductivity tensor
␣
共i兲
, which relates an applied temperature gradient in the

direction to the resulting heat current per unit area in the
␣
direction through Fourier’s law, J
␣
=−兺

␣
T/
x

. Here,
␣
共i兲
can be expressed as
4
␣
共i兲
=
兺
C
v
␣
v


, 共4兲
where C
=k
B
共ប
/ k
B
T兲
2
n
0
共n
0
+1兲 is the contribution per
mode 共j, q兲 to the specific heat, and the scattering times
␣
are directly proportional to the ⌿
.
10
The only inputs required to solve the PBE are the har-
monic and the anharmonic IFCs, ⌽
␣
共0
;ᐉ
⬘
⬘
兲 and
⌽
␣␥
共0
,ᐉ
⬘
⬘
,ᐉ
⬙
⬙
兲. In a previous study,
10
we calculated
the lattice thermal conductivity of bulk silicon combining a
full iterative solution of the PBE with harmonic and anhar-
monic IFCs obtained from three commonly used EIPs.
17
However, the resulting values of
共i兲
were several times
higher than the measured values.
10
This large disagreement
can be traced to inaccurate representation of both harmonic
and anharmonic IFCs obtained from the EIPs.
Since accurate harmonic and anharmonic IFCs are cru-
cial elements required to calculate
共i兲
, in this work we em-
ploy a first principles approach using density functional per-
turbation theory. Density functional perturbation theory
approaches have demonstrated good predictive ability for
phonon dispersions
11
and phonon lifetimes,
18
and we use
them here to calculate the harmonic IFCs, 关⌽
␣
共0
;ᐉ
⬘
⬘
兲兴.
The third order anharmonic IFCs, 关⌽
␣␥
共0
,ᐉ
⬘
⬘
,ᐉ
⬙
⬙
兲兴,
are determined using the 2n + 1 theorem.
19
This theorem
states that if the derivatives of the wave function up to order
n are known, then it is possible to calculate the energy de-
rivatives for the system up to order 2n+ 1. We first use the
Fourier transform expression
⌽
␣␥
⬘
⬙
共q,q
⬘
,q
⬙
兲
=
␦
q±q
⬘
,q
⬙
+K
兺
ᐉ
⬘
,ᐉ
⬙
⌽
␣␥
共0
,ᐉ
⬘
⬘
,ᐉ
⬙
⬙
兲e
±iq
⬘
·R
ᐉ
⬘
e
iq
⬙
·R
ᐉ
⬙
共5兲
and the ⌽
␣␥
⬘
⬙
共q , q
⬘
,q
⬙
兲 are calculated as described in Ref.
12. Previous calculations of the anharmonic IFCs focused on
decay of longitudinal optic and transverse optic phonons at
the ⌫ point where q = 0 and q
⬘
=−q
⬙
.
18
For calculation of the
thermal conductivity, this is insufficient. In our solution of
the PBE, hundreds of thousands of energy and momentum
conserving three-phonon scattering events throughout the
Brillouin zone are required. To our knowledge, Ref. 12 is the
first to calculate the ⌽
␣␥
⬘
⬙
共q , q
⬘
,q
⬙
兲 for general 共q,q
⬘
,q
⬙
兲,
and we apply this approach here for Si and Ge. The
⌽
␣␥
⬘
⬙
共q , q
⬘
,q
⬙
兲 are evaluated directly on a 4 ⫻ 4⫻4 mesh
for pairs of vectors 兵q, q
⬘
其 with q
⬙
fixed by the translational
symmetry condition q
⬙
=q± q
⬘
−K. Crystal symmetry re-
duces the total number of pairs that need to be evaluated to
42. The linear equations, Eq. 共5兲, are then solved for the
⌽
␣␥
共0
,ᐉ
⬘
⬘
,ᐉ
⬙
⬙
兲, with triplet interactions computed out
to seventh nearest neighbors. For all density functional cal-
culations, an 8 ⫻8⫻ 8 Monkhorst-Pack mesh was used in the
Brillouin Zone.
20
An energy cutoff of 24 Ry was used for the
plane wave expansion. The pseudopotentials were generated
based on the approach of von Barth and Car.
21
231922-2 Broido et al. Appl. Phys. Lett. 91, 231922 共2007兲
Downloaded 10 Dec 2007 to 128.253.198.179. Redistribution subject to AIP license or copyright; see http://apl.aip.org/apl/copyright.jsp
We calculate the thermal conductivity,
共i兲
=
zz
共i兲
,ofSi
and Ge for heat current parallel to the direction of the tem-
perature gradient, which is taken to lie along the 关001兴共z兲
direction of the diamond lattice. The cubic symmetry is veri-
fied numerically by assuring that
zz
共i兲
=
xx
共i兲
=
yy
共i兲
and that
xy
共i兲
=
yz
共i兲
=
xz
共i兲
=0. Contributions to the intrinsic lattice thermal
conductivity for different phonon branches are obtained by
summing the integrand in Eq. 共4兲 over all wave vectors. We
find that the contribution is dominated by the acoustic
branches. Specifically, at room temperature, the acoustic
branches provide 95% of the contribution to the thermal con-
ductivity for Si and 92% for Ge. This, however, does not
mean that the optic phonons are unimportant. In fact, to the
contrary, they provide a dominant scattering channel for the
heat-carrying acoustic phonons which, if removed would
precipitate a dramatic increase in the thermal
conductivity.
9,10
The calculated intrinsic lattice thermal conductivities for
silicon and germanium between 100 and 300 K are com-
pared with measured values
13,14
in Fig. 1. The high isotopic
purity 共99.9%
28
Si and 99.9%
70
Ge兲 and quality of these
samples ensure that the dominant mechanism that limits the
thermal conductivity is phonon-phonon scattering. The
agreement between theory and experiment seen in Fig. 1 is
exceptionally good. For example, at room temperature the
calculated value for Si 共Ge兲 is only 2% 共5%兲 different than
the corresponding measured value. This is particularly no-
table given the poor agreement obtained by us using the EIP
models for the IFCs.
10
Previous molecular dynamics calcu-
lations of
共i兲
using different EIPs for Si obtained room tem-
perature values of about 2.4 共Ref. 7兲 and 1.4 W/ cm K,
8
about 60% higher and 10% lower, respectively, than the mea-
sured value.
13
This wide variation could be due to the sensi-
tivity of the result to the EIP used in each calculation and
that molecular dynamics approaches are expected to be more
accurate at high temperatures.
The computational rigor of the present approach cur-
rently precludes consideration of bulk materials with very
large unit cells.
22
However, in some nanosystems this is not
the case. For example, single-walled carbon nanotubes
whose unit cells may contain several tens of atoms could
nevertheless be treated using the present theory because the
phonon wave vectors are only one dimensional so that the
phase space for phonon-phonon scattering is dramatically re-
ducing compared to a bulk system with similar sized unit
cell.
To our knowledge, the theoretical results presented here
are the only ones to date to provide excellent agreement with
measured lattice thermal conductivities of any material over
a substantial temperature range using a first-principles-based
adjustable-parameter-free approach. The demonstrated accu-
racy of our theory applied to silicon and germanium suggests
its predictive potential for calculating the
共i兲
of many other
bulk and nanostructured materials of scientific and techno-
logical interest. Calculations of
共i兲
for other bulk materials
are currently underway to further verify the predictive capa-
bility of the theory.
Part of this work was done on the Intel Computing Clus-
ter at the Cornell Nanoscale Facility, part of the National
Nanotechnology Infrastructure Network supported by the
National Science Foundation. Acknowledgment is also made
to the Donors of the American Chemical Society Petroleum
Research Fund and to the National Science Foundation for
support of this research. We thank A. Inyushkin for kindly
providing us with his measured lattice thermal conductivity
data for silicon and germanium and A. Mayer, G. Deinzer,
and D. Strauch for useful communications on this topic.
1
D. G. Cahill, W. K. Ford, K. E. Goodson, G. D. Mahan, A. Majumdar, H.
J. Maris, R. Merlin, and S. R. Phillpot, J. Appl. Phys. 93,793共2003兲.
2
G. D. Mahan, B. Sales, and J. Sharp, Phys. Today 50,42共1997兲;F.Di
Salvo, Science 285,703共1999兲.
3
R. Peierls, Ann. Phys. 3, 1055 共1929兲.
4
R. Peierls, Quantum Theory of Solids 共Clarendon, Oxford, 1955兲,p.40;J.
M. Ziman, Electrons and Phonons 共Oxford University Press, London,
1960兲, p. 298.
5
M. Asen-Palmer, K. Bartkowski, E. Gmelin, M. Cardona, A. P. Zhernov,
A. V. Inyushkin, A. Taldenkov, V. I. Ozhogin, K. M. Itoh, and E. E. Haller,
Phys. Rev. B 56, 9431 共1997兲.
6
See, for example, J. X. Cao, X. H. Yan, Y. Xiao, and J. W. Ding, Phys.
Rev. B 69, 073407 共2004兲; A. Khitun and K. L. Wang, Appl. Phys. Lett.
79, 851 共2001兲.
7
S. G. Volz and G. Chen, Phys. Rev. B 61, 2651 共2000兲; P. K. Schelling, S.
R. Phillpot, and P. Keblinski, ibid. 65, 144306 共2002兲.
8
L. Sun and J. Y. Murthy, Appl. Phys. Lett. 89, 171919 共2006兲.
9
M. Omini and A. Sparavigna, Phys. Rev. B 53, 9064 共1996兲;M.Omini
and A. Sparavigna, Nuovo Cimento D 19, 1537 共1997兲.
10
D. A. Broido, A. Ward, and N. Mingo, Phys. Rev. B 72, 014308 共2005兲.
11
S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod.
Phys. 73,515共2001兲.
12
G. Deinzer, G. Birner, and D. Strauch, Phys. Rev. B 67, 144304 共2003兲.
13
A. V. Inyushkin, A. N. Taldenkov, A. M. Cibin, A. V. Gusev, and H.-J.
Pohl, Phys. Status Solidi C 1, 2995 共2004兲; R. K. Kremer, K. Graf, M.
Cardona, G. G. Devyatykh, A. V. Gusev, A. M. Gibin, A. V. Inysushkin,
A. N. Taldenkov, and H.-J. Pohl, Solid State Commun. 131, 499 共2004兲.
14
V. I. Ozhogin, A. V. Inyushkin, A. N. Taldenkov, A. V. Tikhomirov, and
G. E. Popov, JETP Lett. 63,490共1996兲.
15
A simple model comparison of three and four-phonon scattering rates
shows the latter to be considerably weaker than the former even at high
temperatures. See D. J. Ecsedy and P. G. Klemens, Phys. Rev. B 15,5957
共1977兲.
16
P. G. Klemens, Proc. R. Soc. London, Ser. A A208, 108 共1951兲.
17
F. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 共1985兲; J. Tersoff,
ibid. 38, 9902 共1988兲; J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov,
and S. Yip, ibid. 58, 2539 共1998兲.
18
A. Debernardi, S. Baroni, and E. Molinari, Phys. Rev. Lett. 75, 1819
共1995兲; G. Lang, K. Karch, M. Schmitt, P. Pavone, A. P. Mayer, R. K.
Wehner, and D. Strauch, Phys. Rev. B 59, 6182 共1999兲.
19
X. Gonze and J.-P. Vigneron, Phys. Rev. B 39, 13120 共1989兲.
20
H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 共1976兲.
21
U. von Barth and R. Car, 共unpublished兲; for a brief description of this
method, see A. Dal Corso, S. Baroni, R. Resta and S. de Gironcoli, Phys.
Rev. B 47, 3588 共1993兲.
22
Calculation of the phase space of energy and momentum-conserving three-
phonon processes provides both time and memory bottlenecks. The size of
the phase space search increases dramatically with the number phonon
branches 3m, where m is the number of atoms in the unit cell.
FIG. 1. 共Color online兲 Lattice thermal conductivity
共i兲
as a function of
temperature. The red line and solid squares are the calculated and measured
thermal conductivities of silicon, respectively, while the blue line and the
solid circles are the corresponding quantities for germanium.
231922-3 Broido et al. Appl. Phys. Lett. 91, 231922 共2007兲
Downloaded 10 Dec 2007 to 128.253.198.179. Redistribution subject to AIP license or copyright; see http://apl.aip.org/apl/copyright.jsp
