Lattice thermal conductivity of Bi, Sb,
and Bi-Sb alloy from first principles
The MIT Faculty has made this article openly available. Please share
how this access benefits you. Your story matters.
Citation Lee, Sangyeop, Keivan Esfarjani, Jonathan Mendoza, Mildred S.
Dresselhaus, and Gang Chen. “Lattice Thermal Conductivity of Bi,
Sb, and Bi-Sb Alloy from First Principles.” Phys. Rev. B 89, no. 8
(February 2014). © 2014 American Physical Society
As Published http://dx.doi.org/10.1103/PhysRevB.89.085206
Publisher American Physical Society
Version Final published version
Citable link http://hdl.handle.net/1721.1/88767
Terms of Use Article is made available in accordance with the publisher's
policy and may be subject to US copyright law. Please refer to the
publisher's site for terms of use.
PHYSICAL REVIEW B 89, 085206 (2014)
Lattice thermal conductivity of Bi, Sb, and Bi-Sb alloy from first principles
Sangyeop Lee,
1
Keivan Esfarjani,
2,3
Jonathan Mendoza,
1
Mildred S. Dresselhaus,
4,5
and Gang Chen
1
1
Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
2
Department of Mechanical and Aerospace Engineering, Rutgers University, New Brunswick, NJ 08901, USA
3
Institute for Advanced Materials and Devices for Nanotechnology (IAMDN), Rutgers University, Piscataway, NJ 08854, USA
4
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
5
Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
(Received 20 December 2013; revised manuscript received 12 February 2014; published 25 February 2014)
Using first principles, we calculate the lattice thermal conductivity of Bi, Sb, and Bi-Sb alloys, which are
of great importance for thermoelectric and thermomagnetic cooling applications. Our calculation reveals that
the ninth-neighbor harmonic and anharmonic force constants are significant; accordingly, they largely affect
the lattice thermal conductivity. Several features of the thermal transport in these materials are studied: (1) the
relative contributions from phonons and electrons to the total thermal conductivity as a function of temperature are
estimated by comparing the calculated lattice thermal conductivity to the measured total thermal conductivity,
(2) the anisotropy of the lattice thermal conductivity is calculated and compared to that of the electronic
contribution in Bi, and (3) the phonon mean free path distributions, which are useful for developing nanostructures
to reduce the lattice thermal conductivity, are calculated. The phonon mean free paths are found to range from 10 to
100 nm for Bi at 100 K.
DOI: 10.1103/PhysRevB.89.085206 PACS number(s): 63.20.kg, 66.70.−f
I. INTRODUCTION
Bi and Bi-Sb alloys have long been studied for their
promising low-temperature thermoelectric applications. Bi
and Sb have a rhombohedral crystal structure, which is a Peierl
distortion of the simple cubic crystal. The small structural
distortion results in Brillouin zone folding and a small overlap
between conduction and valence bands, thereby causing
semimetallic behavior and conduction by both electrons and
holes. Since the semimetallic behavior causes cancellation of
the hole and electron contributions to the power factor, bulk
Bi is not a good thermoelectric material. However, Bi has a
large thermomagnetic effect and a large thermomagnetic figure
of merit (ZT)[1]. The thermomagnetic effect is particularly
pronounced below 10 K due to the extremely long mean free
path of the electrons in Bi [2]. In addition, Bi nanowires
become semiconducting as their diameters approach several
nanometers, thereby exhibiting a large thermoelectric power
factor [3,4]. As a conventional bulk thermoelectric material,
Bi
1−x
Sb
x
has drawn more attention than Bi, since alloying with
a small amount of Sb causes Bi
1−x
Sb
x
to become a narrow gap
semiconductor, which is advantageous for high thermoelectric
efficiency. Currently, Bi
1−x
Sb
x
(x 0.12) is the best n-type
thermoelectric material below 200 K [5].
Before discussing the lattice thermal transport, we empha-
size that electrons, in addition to phonons, carry a considerable
amount of heat in Bi, Sb, and Bi-Sb alloys. Therefore,
both phonons and electrons contribute to the total thermal
conductivity, which can be expressed as
κ
tot
= κ
ph
+ κ
e
(1)
where κ
tot
and κ
ph
are the total thermal conductivity and the
lattice thermal conductivity, respectively. The term κ
e
includes
the thermal conductivity of electrons and holes, as well as
the bipolar contribution. The κ
e
of Bi, Sb, and Bi-Sb alloys
is expected to contribute substantially to κ
tot
, since these
materials are either semimetals or semiconductors with a very
narrow band gap.
Accurate methods to obtain separate κ
ph
and κ
e
are crucial
to developing better thermoelectric materials, but separating
κ
ph
and κ
e
is experimentally nontrivial. κ
ph
can be directly
measured under a high magnetic field, because such fields
largely suppress electron transport. Previous measurements
[6,7] in practical temperature ranges (100–300 K) utilized this
method, but the prior measurements are mainly limited to
transport along the binary crystallographic direction. We could
not find any reports on κ
ph
along the trigonal direction, which
is expected to have a greater ZT than for the binary direction
and thus is of more interest. Another way to separate κ
ph
and
κ
e
is to estimate κ
e
using either the Wiedemann-Franz law
or other electron transport properties, such as the electrical
conductivity and Seebeck coefficient [8]. Such an approach
provides a reasonable qualitative analysis, but validity of
the Wiedemann-Franz law and the simple electron transport
models used in the estimation of κ
e
is sometimes questionable
for quantitative purposes [9].
In this paper, we s tudy the lattice dynamics and quantify
κ
ph
for Bi, Sb, and Bi-Sb alloys from first principles and the
Boltzmann transport equation. As shown in recent papers
[10–16], this approach provides excellent agreement with
experimental data for many pair-bonded materials, such as Si,
GaAs, and Si-Ge alloys. We follow the same approach but pay
special attention to the range of interatomic interactions. This
is because Bi, unlike pair-bonded materials, has significant
interaction strength out to large number neighbors, such as the
ninth-nearest neighbor [17,18].
II. SECOND- AND THIRD-ORDER FORCE CONSTANTS
In this paper, we calculated the second- and third-order
force constants using density functional theory. The calcula-
tion of the second-order force constants of Bi and Sb is based
on the real space approach [19]. We calculated the force exerted
1098-0121/2014/89(8)/085206(9) 085206-1 ©2014 American Physical Society
LEE, ESFARJANI, MENDOZA, DRESSELHAUS, AND CHEN PHYSICAL REVIEW B 89, 085206 (2014)
FIG. 1. Crystal structure of Bi and Sb. The void and filled atoms
represent two basis atoms. R
1
,R
2
,andR
3
are primitive lattice vectors,
and α is a rhombohedral angle. The values of α are 57°30
for Bi and
57°84
for Sb, which are close to 60° for the simple cubic structure.
on each atom when we displace one or multiple atoms in a
4 ×4 ×4 supercell (128 atoms). For the supercell calculation,
we used 30 Ry for the cutoff energy of the plane wave basis
and a 4 × 4 × 4 k-point mesh for Brillouin zone sampling,
both of which were carefully checked for convergence. The
calculation was performed with the ABINIT package [20]
and Hartwigsen-Goedecker-Hutter pseudopotentials [21]. The
valence electrons in the pseudopotential are 6s
2
6p
3
and 5s
2
5p
3
for Bi and Sb, respectively. The spin-orbit interaction is
included in all calculations because of the strong spin-orbit
interaction in Bi and Sb [22]. The second-order force constants
are then fitted to the calculated displacement-force data set
while enforcing translational and rotational invariance. In
the fitting process, we considered up to the 14th neighbors
to include the previously reported long-ranged interaction
occurring at the ninth neighbor [17,18]. The ninth neighbors
are s hown by the atom labeled C in Fig. 1, where the
origin atom is described by atom A. Bi and Sb both have a
slightly distorted simple cubic crystal structure. Due to this
small crystallographic distortion, the six first neighbors in the
cubic structure become three first neighbors and three second
neighbors. In Fig. 1, atom B is the first neighbor to atom A
and the second neighbor to atom C. The almost collinear chain
consisting of
AB and BC forms the ninth-neighbor relation,
and atom C is the ninth neighbor to atom A. In the following
discussions, the fourth and ninth neighbors are frequently men-
tioned to discuss the range of the force constants. The fourth
and ninth neighbors in the rhombohedral crystal structure of
Bi correspond to the second neighbor (separated by
√
2a)
and the fourth neighbor (separated by 2a), respectively, in the
undistorted cubic structure, where a is the lattice constant of
the simple cubic structure.
The third-order force constants were calculated by taking
finite differences of the second-order force constants [23]. We
built a 3 × 3 × 3 supercell consisting of 54 atoms, and we
displaced one of the two basis atoms along the +R
1
direction in
Fig. 1 by 0.04
˚
A. The displacement value of 0.04
˚
A was chosen
after carefully checking the convergence of third-order force
constants with respect to the displacement values. The size of
the supercell was large enough to include the significant ninth-
neighbor interaction. In addition, the large size of the supercell
minimizes the effect from the periodic images of the displaced
atom due to periodic boundary conditions. For the cal-
culation, a cutoff energy of 30 Ry and a 3 × 3 × 3 k-point
mesh are used. We then calculate the second-order force
constants using density functional perturbation theory [24,25].
All of the procedures are repeated for another supercell with
the displacement along the −R
1
direction. By taking the
finite differences of the second-order force constants of the
two different supercells, the third-order force constants with
respect to the R
1
direction are calculated. Rotational invariance
with respect to the trigonal direction is then applied to calculate
the third-order force constants with respect to the R
2
and R
3
directions. Translational invariance is applied to the third-order
force constants by modifying the self-interaction terms.
We calculated phonon dispersions and mode Gr
¨
uneisen
parameters to validate the calculated second- and third-order
force constants. In Fig. 2(a), we plot the trace of the
Trace of 2nd order FC tensors (eV/Å
2
)
distance (Å)
distance (Å)
3rd order force constants (eV/Å
3
)
9th neighbors
9th neighbors
(a)
(b)
FIG. 2. (Color online) Force constants (FCs) of Bi and Sb versus
interatomic distance: (a) Trace values of second-order FC tensors and
(b) two-body third-order FCs.
085206-2
LATTICE THERMAL CONDUCTIVITY OF Bi, Sb, AND . . . PHYSICAL REVIEW B 89, 085206 (2014)
Frequency (THz)Frequency (THz)
XK Γ T W L Γ X
XK Γ T W L Γ X
W
T
Γ
L
X
K
(b)
(a)
(c)
FIG. 3. (Color online) Phonon dispersion of (a) Bi and (b) Sb.
Dots are experimental values from [37]forBiand[38]forSb.(c)
The high symmetry points in the Brillouin zone.
second-order force constant tensors versus distance. Both Bi
and Sb have the interactions of significant magnitude occurring
at the ninth neighbors, which agree well with the previous
reports [17,18]. In Fig. 3, the calculated phonon dispersions
for Bi and Sb are compared with the experimental values. Both
calculated phonon dispersions are similar to the experimental
data, confirming the accuracy of the calculated second-order
force constants.
Since the ninth neighbors in Bi and Sb have significant
second-order force constants, the third-order force constants
at the ninth neighbors should also be of interest. In Fig. 2(b),
we plot the two-body third-order force constants as a function
of distance. Each dot represents a third-order force constant.
As seen in Fig. 2(b), the third-order force constants have
substantial values at the ninth neighbors. The importance of
the ninth-neighbor interaction on crystal anharmonicity can
be checked with the mode Gr
¨
uneisen parameters. The mode
Gr
¨
uneisen parameters are calculated with the two different
sets of third-order force constants: one includes up to the
fourth neighbors and the other includes up to the 10th
neighbors. To validate the third-order force constants, the
reference mode Gr
¨
uneisen parameters are also calculated.
For the reference mode Gr
¨
uneisen parameters, we used
density functional perturbation theory to calculate the phonon
frequencies for two different crystal volumes: a crystal at
equilibrium and one with the volume increased by 1%. We
then take the finite differences of the two different phonon
frequencies and calculate the mode Gr
¨
uneisen parameters
from the definition γ =−dlnω/dlnV , where ω and V are a
phonon frequency and a crystal volume, respectively. Shown in
Fig. 4 are the calculated acoustic mode Gr
¨
uneisen parameters.
Figure 4 shows that the acoustic mode Gr
¨
uneisen parameters
are underestimated over a range of wave vectors when the
third-order force constants are considered only up to the fourth
neighbors. Even after considering up to the eighth neighbors,
the mode Gr
¨
uneisen parameters are relatively unchanged. This
is consistent with the negligible third-order force constants
at the fifth, sixth, seventh, and eighth neighbors, as shown
in Fig. 2(b). However, when extending the range up to
the 10th neighbors, the calculated acoustic mode Gr
¨
uneisen
parameters agree reasonably well with the reference Gr
¨
uneisen
XK Γ T W L Γ X
XK Γ T W L Γ X
Grüneisen parameterGrüneisen parameter
(a)
(b)
Reference
upto 4th neighbor
upto 10th neighbor
Reference
upto 4th neighbor
upto 10th neighbor
FIG. 4. (Color online) Acoustic mode Gr
¨
uneisen parameters of
(a) Bi and (b) Sb comparing inclusion up to the fourth and 10th
neighbors to the references. The reference Gr
¨
uneisen parameters are
calculated using the difference of phonon frequencies of the two
different crystal volumes.
085206-3
LEE, ESFARJANI, MENDOZA, DRESSELHAUS, AND CHEN PHYSICAL REVIEW B 89, 085206 (2014)
parameters. This confirms that the ninth-neighbor interaction
is playing a significant role in the anharmonic properties. The
optical mode Gr
¨
uneisen parameter was also determined from
third-order force constants that included up to the fourth- and
10th-neighbor interaction terms. Both cases yielded similar
values for t he optical Gr
¨
uneisen parameter.
The significant interaction at the ninth-nearest neighbors
can be explained by the resonant bonding in Bi and Sb [26]. Bi
and Sb have very weak sp hybridization, and the s band is well
below the p band [27]. Therefore, s electrons do not participate
in the chemical bonding, and we can consider only p electrons
forming the chemical bonds. For the three p electrons per atom
in Bi or Sb to meet the six coordination number requirement
in the cubiclike crystal structure, the electrons alternate their
positions among six chemical bonds, leading the chemical
bonding called resonant bonding [28]. This resonant bonding
picture implies two important features: (1) electrons are highly
delocalized and are therefore easily polarized upon external
perturbations, and (2) the chemical bonds in Bi and Sb are
almost collinear due to the cubiclike crystal structure. The
almost collinear bonding can be found in Fig. 1, as explained
earlier. These two features result in the significant interaction
at the ninth-nearest neighbors. The electron polarization by
the displacement of the origin atom is long ranged along
the collinear bonding direction due to the large electronic
polarizability and almost collinear bonding. This long-ranged
electron polarization reaches the ninth-nearest neighbors,
giving rise to the significant interatomic interaction between
the origin and the ninth-nearest neighbor atoms.
To study the effects of alloying on κ
ph
, the virtual crystal
approximation is used [29]. The atomic mass and the force
constants of the virtual crystal were linearly interpolated
between Bi and Sb, weighted by the composition ratio of the
constituents. The lattice constant of the virtual crystal is also
averaged according to the composition ratio, which is well
justified by the fact that the Bi-Sb alloy follows Vegard’s
law [30]. Three-phonon scattering is calculated using the
virtual crystal approximation, while the atomic mass disorder
is treated as an additional elastic scattering mechanism. This
approach is successful in predicting the Si-Ge alloy thermal
conductivity [16].
III. SCATTERING RATE AND LATTICE
THERMAL CONDUCTIVITY
The lattice thermal conductivity can be calculated from
the distribution function of the phonon modes. We calculate
the distribution function by solving the linearized Boltzmann
equation with the scattering rates due to the three-phonon
process and mass disorder. The scattering rate of the three-
phonon process i s given by
W
3
1,2
= 2π|V
3
(−1,−2,3)|
2
n
0
1
n
0
2
n
0
3
+ 1
δ(−ω
1
− ω
2
+ ω
3
)
for coalescence processes (2)
W
2,3
1
= 2π|V
3
(−1,2,3)|
2
n
0
1
n
0
2
+1
n
0
3
+ 1
δ(−ω
1
+ω
2
+ω
3
)
for decay processes (3)
where 1, 2, and 3 denote phonon modes in the three-phonon
process, while n
0
and ω indicate the Bose-Einstein equilibrium
distribution function and the phonon frequency, respectively.
The three-phonon scattering matrix element V
3
is given by
V
3
(1,2,3) =
8Nω
1
ω
2
ω
3
1/2
×
b
1
b
2
b
3
αβγ
R
2
R
3
αβγ
(0b
1
,R
2
b
2
,R
3
b
3
)
×e
iq
2
·R
2
e
iq
3
·R
3
e
αb
1
e
βb
2
e
γb
3
√
m
b
1
m
b
2
m
b
3
(4)
where
αβγ
(0b
1
,R
2
b
2
,R
3
b
3
) is a third-order force constant
with Cartesian coordinates αβγ and Rb representing the
lattice vector and basis atom. Here, e
αb
denotes the phonon
eigenvector component of the basis atom b along direction
α, while N is the total number of wave vectors in the first
Brillouin zone. The mass disorder scattering rate is
W
2
1
=
π
2
gω
1
ω
2
n
0
1
n
0
2
+ 1
b
|e
∗
b
· e
b
|
2
δ(ω
1
− ω
2
)(5)
with the mass variance factor g defined by g =
i
f
i
(1 − M
i
/M
avg
)
2
, where f
i
is the fraction of element i.
Putting both scattering rates above into the Boltzmann
equation, we obtain
− v
1
·∇T
∂n
0
1
∂T
=
2,3
W
3
1,2
(
1
+
2
−
3
)
+
1
2
W
2,3
1
(
1
−
2
−
3
)
+
2
W
2
1
(
1
−
2
)(6)
where is a linearized deviation of the distribution function
from equilibrium, defined as = (n
0
− n)/(∂n
0
/∂β) with
β = ω/k
B
T . We solve the linearized Boltzmann equation
above iteratively to find /∇T . The detailed procedure is
provided in other papers [12]. In contrast to the iterative
method mentioned above, a more commonly used method
to solve the Boltzmann equation is to neglect
2
and
3
in
Eq. (6), and this approximation is known as the single-mode
relaxation time (SMRT) approximation. The SMRT assumes
only one phonon mode is ever out of equilibrium. The
time for the nonequilibrium mode to relax to equilibrium
is then calculated. We used both the full iterative method
and the SMRT approximation to calculate the lattice thermal
conductivity from the Boltzmann equation, and we compare
the results from the two methods. After solving the Boltzmann
equation, the lattice thermal conductivity can be obtained by
κ
αβ
=−
1
V
BZ
ωv
α
n
0
(n
0
+ 1)
β
∇
β
T
(7)
where α and β are Cartesian directions, V is the crystal volume,
and v is the phonon group velocity.
One of the numerical uncertainties in this calculation occurs
in the energy conservation of the scattering rate calculation.
Due to the computational limitations, the Brillouin zone is
sampled with a relatively coarse mesh. To find sets of three
085206-4
