of 9
Long-range quadrupole electron-phonon interaction from first principles
Jinsoo Park,
1
Jin-Jian Zhou,
1
Vatsal A. Jhalani,
1
Cyrus E. Dreyer,
2, 3
and Marco Bernardi
1,
1
Department of Applied Physics and Materials Science,
California Institute of Technology, Pasadena, California 91125, USA.
2
Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11794-3800
3
Center for Computational Quantum Physics, Flatiron Institute, 162 Fifth Avenue, New York, New York 10010
Lattice vibrations in materials induce perturbations on the electron dynamics in the form of
long-range (dipole and quadrupole) and short-range (octopole and higher) potentials. The dipole
Fr ̈ohlich term can be included in current first-principles electron-phonon (
e
-ph) calculations and
is present only in polar materials. The quadrupole
e
-ph interaction is present in both polar and
nonpolar materials, but currently it cannot be computed from first principles. Here we show an
approach to compute the quadrupole
e
-ph interaction and include it in
ab initio
calculations of
e
-ph matrix elements. The accuracy of the approach is demonstrated by comparing with direct
density functional perturbation theory calculations. We apply our method to silicon as a case of a
nonpolar semiconductor and tetragonal PbTiO
3
as a case of a polar piezoelectric material. In both
materials we find that the quadrupole term strongly impacts the
e
-ph matrix elements. Analysis of
e
-ph interactions for different phonon modes reveals that the quadrupole term mainly affects optical
modes in silicon and acoustic modes in PbTiO
3
, although the quadrupole term is needed for all
modes to achieve quantitative accuracy. The effect of the quadrupole
e
-ph interaction on electron
scattering processes and transport is shown to be important. Our approach enables accurate studies
of
e
-ph interactions in broad classes of nonpolar, polar and piezoelectric materials.
I. INTRODUCTION
Electron-phonon (
e
-ph) interactions are key to under-
standing electrical transport, nonequilibrium dynamics,
and superconductivity [1]. First-principles calculations
can provide microscopic insight into
e
-ph scattering pro-
cesses and are rapidly emerging as a quantitative tool
for investigating charge transport and ultrafast carrier
dynamics in materials [2–13]. The typical workflow com-
bines density functional theory (DFT) [14] calculations
of the ground state and band structure with density func-
tional perturbation theory (DFPT) [15] for phonon dis-
persions and
e
-ph perturbation potentials. As DFPT can
compute the electronic response to periodic lattice per-
turbations (phonons) with arbitrary wave-vector
q
, the
DFPT framework can capture both short- and long-range
e
-ph interactions.
However, a key challenge is that DFPT is too compu-
tationally demanding to be carried out on the fine Bril-
louin zone grids needed to compute electron scattering
rates and transport properties. The established approach
in first-principles
e
-ph studies [9] is to carry out DFPT
calculations on coarse Brillouin zone grids with of order
10
×
10
×
10
q
-points, followed by interpolation of the
e
-ph
matrix elements with a localized basis set such as Wan-
nier functions or atomic orbitals [16]. As the perturba-
tion potential can be non-analytic near
q
= 0 [16] or even
exhibit a divergence for certain phonon modes, interpo-
lation is particularly challenging and less reliable in the
region between
q
= 0 and its nearest-neighbor
q
-points
in the coarse DFPT grid. This small-
q
region is critical
as it is dominated by long-range
e
-ph interactions, whose
treatment can affect the quality of the interpolation even
at larger values of
q
in the Brillouin zone.
A multipole expansion of the
e
-ph perturbation po-
tential shows that in the long-wavelenght limit (phonon
wave-vector
q
0) the long-range dipole Fr ̈ohlich term
diverges as 1
/q
, the quadrupole term approaches a con-
stant value and the short-range octopole and higher
terms vanish [17, 18]. These trends in momentum space
are due to the spatial decay of the
e
-ph interactions, with
1
/r
2
trend for the dipole, 1
/r
3
for the quadrupole, and
1
/r
4
or faster for the short-range part. The open question
is how one can carry out the
e
-ph matrix element interpo-
lation in the region near
q
= 0 using analytical expressions
for the long-range dipole and quadrupole terms. These
expressions have been obtained by Vogl [18], but need
to be rewritten in the
ab initio
formalism and computed
with first-principles quantities such as the atomic dynam-
ical dipoles [15] and quadrupoles [19–21] induced by lat-
tice vibrations, which can be computed with DFPT. For
each atom
κ
, one can obtain Born charge (
Z
κ
) and dy-
namical quadrupole (
Q
κ
) tensors, which, once contracted
with the phonon eigenvector, give the atomic contribu-
tions to the dipole and quadrupole
e
-ph interactions.
The dipole Fr ̈ohlich term has been derived following
this strategy [22, 23] and employed in electron scattering
rate and transport calculations [4, 7]. The quadrupole
term has not yet been derived or implemented in first-
principles calculations and its important effect on the
e
-
ph matrix elements has been overlooked. To understand
the role of long-range dipole and quadrupole
e
-ph inter-
actions, it is useful to consider separately the interactions
for different phonon modes in the long-wavelength limit,
discerning the effect of longitudinal and transverse, and
acoustic and optical modes. Analytical models of
e
-ph
interactions rely on such an intuition for the role of dif-
ferent phonon modes in various materials [24].
In ionic and polar covalent crystals (here and below,
denoted as polar materials), the dipole Fr ̈ohlich term is
dominant as
q
0 due to its 1
/q
trend. This
e
-ph inter-
arXiv:2003.13782v1 [cond-mat.mtrl-sci] 30 Mar 2020
2
action is due to longitudinal optical (LO) phonons and it
dominates small-
q
scattering. For other phonon modes in
polar materials, the dipole term vanishes, and the domi-
nant long-range
e
-ph interaction is the quadrupole term,
which is particularly important for acoustic phonons
in piezoelectric (polar noncentrosymmetric) materials.
In nonpolar semiconductors such as silicon and germa-
nium, the dipole Fr ̈ohlich interaction vanishes and the
quadrupole term is a key contribution for all modes in
the long-wavelength limit. The quadrupole
e
-ph interac-
tion is thus expected to play an important role in many
classes of materials, making a compelling case for its in-
clusion in the first-principles framework.
Here we show
ab initio
calculations of the long-range
quadrupole
e
-ph interaction and an approach to include
it in the
e
-ph matrix elements. The accuracy of our
method is confirmed by comparing the
e
-ph matrix el-
ements with direct DFPT calculations. We find that the
quadrupole contribution is significant for most phonon
modes in both nonpolar and polar materials. In silicon,
a nonpolar semiconductor, the quadrupole term has a
large effect on the
e
-ph coupling for optical modes, but
is negligible for acoustic modes in the long-wavelength
limit. In tetragonal PbTiO
3
, a polar piezoelectric ma-
terial, the quadrupole corrections are substantial for all
phonon modes and particularly important for acoustic
modes, which contribute to the piezoelectric
e
-ph inter-
action. Including only the long-range Fr ̈ohlich interaction
and neglecting the quadrupole term leads to large errors
in PbTiO
3
, while adding the quadrupole term leads to
e
-
ph matrix elements that accurately reproduce the DFPT
benchmark results for all phonon modes in the entire Bril-
louin zone. We investigate the impact of the quadrupole
e
-ph interaction on the electron scattering rates and mo-
bility in silicon and PbTiO
3
, finding mobility corrections
of order 10% in silicon and 20% in PbTiO
3
at 100 K
(and smaller corrections at 300 K) when the quadrupole
term is included. The correction on the scattering rate
at low electron energy in PbTiO
3
is substantial. Taken
together, our results highlight the need to include the
quadrupole term in all materials to correctly capture the
long-range
e
-ph interactions. In turn, this development
enables more precise calculations of electron dynamics
and scattering processes from first principles.
FIG. 1. Schematic of the dipole and quadrupole charge con-
figurations giving rise to long-range
e
-ph interactions.
II. THEORY
The electron distribution changes in response to a dis-
placement of an atom from its equilibrium position. The
cell-integrated charge response to a displacement of atom
κ
due to a phonon with wave-vector
q
0 can be written
as a multipole expansion [21]:
C
q
κ,α
=
iZ
κ,αβ
q
β
1
2
Q
κ,αβγ
q
β
q
γ
+
...,
(1)
where summation over the Cartesian indices
β
and
γ
is
implied. This polarization response defines the Born ef-
fective charge
Z
κ
, a rank-2 tensor associated with the
dipole term, and the dynamical quadrupole
Q
κ
, the rank-
3 tensor in the quadrupole term; both tensors can be
computed in the DFPT framework [15, 19]. Each of the
dipole and quadrupole responses generates macroscopic
electric fields and corresponding long-range
e
-ph interac-
tions in semiconductors and insulators [18, 25], while in
metals they are effectively screened out.
In a field-theoretic treatment of the
e
-ph interactions,
one computes the dipole and quadrupole perturbation
potentials ∆
V
ν
q
due to a phonon with mode index
ν
and
wave-vector
q
, and the corresponding
e
-ph matrix ele-
ments [2]
g
mnν
(
k
,
q
) =
(
~
2
ω
ν
q
)
1
2
m
k
+
q
|
V
ν
q
|
n
k
,
(2)
which quantify the probability amplitude of an electron
in a Bloch state
|
n
k
with band index
n
and crystal mo-
mentum
k
to scatter into a final state
|
m
k
+
q
by emit-
ting or absorbing a phonon with energy
~
ω
ν
q
.
A. Dipole and quadrupole
e
-ph interactions
To derive the dipole and quadrupole perturbation po-
tentials, we consider a Born-von Karman (BvK) crys-
tal [26] with
N
unit cells and volume
N
Ω. The potential
due to a dipole configuration with dipole moment
p
cen-
tered at position
τ
in the crystal (see Fig. 1) can be
written as [22, 23]
V
dip
(
r
;
τ
) =
i
e
N
ε
0
q
G
6
=
q
p
·
(
q
+
G
)
e
i
(
q
+
G
)
·
(
r
τ
)
(
q
+
G
)
·

·
(
q
+
G
)
,
(3)
where

is the dielectric tensor of the material, the
phonon wave-vector
q
belongs to a regular Brillouin zone
grid with
N
points, and
G
are reciprocal lattice vectors.
This result is derived by adding together the potentials
generated in the crystal by two point charges of opposite
sign with distance
u
0, resulting in a dipole
p
[4, 23].
The potential in Eq. (3) is readily extended to the
case of an atomic dynamical dipole
p
κ,
R
from atom
κ
in
the unit cell at Bravais lattice vector
R
, due to the dis-
placement induced by a phonon with mode index
ν
and
3
wave-vector
q
. The resulting atomic dynamical dipole
is
p
κ,
R
= (
e
Z
κ
)
̃
e
(
κ
)
ν
q
e
i
q
·
R
, where the phonon eigenvec-
tor projected on atom
κ
is defined as
̃
e
(
κ
)
ν
q
=
e
(
κ
)
ν
q
/
M
κ
,
with
e
ν
q
the eigenvector of the dynamical matrix at
q
and
M
κ
the mass of atom
κ
. Summing over the con-
tributions from all atoms
κ
at lattice vectors
R
with
positions
τ
κ
R
=
τ
κ
+
R
in the BvK supercell, the to-
tal
e
-ph dipole interaction due to the phonon mode is
V
dip
ν
q
(
r
) =
κ
R
V
dip
(
r
;
τ
κ
R
).
Using the identity
1
N
R
e
i
q
·
R
=
δ
q
,
0
, we obtain:
V
dip
ν
q
(
r
) =
i
e
2
ε
0
κ
M
1
/
2
κ
G
6
=
q
(
Z
κ
e
(
κ
)
ν
q
)
·
(
q
+
G
)
e
i
(
q
+
G
)
·
(
r
τ
)
(
q
+
G
)
·

·
(
q
+
G
)
(4)
The
ab initio
Fr ̈ohlich
e
-ph coupling is obtained by eval-
uating the matrix elements with this potential:
g
dip
mnν
(
k
,
q
) =
i
e
2
ε
0
κ
(
~
2
ω
ν
q
M
κ
)
1
/
2
G
6
=
q
(
Z
κ
e
(
κ
)
ν
q
)
·
(
q
+
G
)
(
q
+
G
)
·

·
(
q
+
G
)
m
k
+
q
|
e
i
(
q
+
G
)
·
(
r
τ
)
|
n
k
.
(5)
The potential due to the dynamical quadrupole re-
sponse can be derived with a similar strategy.
We
first consider the potential generated by a quadrupole
charge configuration consisting of two equal and oppo-
sitely oriented dipoles
p
and
p
, centered at positions
τ
±
u
2
respectively (see Fig. 1). The configuration, with
quadrupole moment [27]
M
αβ
=
p
α
u
β
, gives a potential:
V
quad
(
r
;
τ
) = lim
u
0
[
V
dip
(
r
;
τ
+
u
2
)
V
dip
(
r
;
τ
u
2
)]
=
e
N
ε
0
q
G
6
=
q
(
q
+
G
)
·
M
·
(
q
+
G
)
(
q
+
G
)
·

·
(
q
+
G
)
e
i
(
q
+
G
)
·
(
r
τ
)
,
(6)
where to obtain the second line we used ∆
V
dip
(
r
;
τ
) in
Eq. (3) and expanded the first line to first order in
u
.
Similar to the dipole case, the potential from atomic
quadrupoles (
M
κ,
R
)
αβ
=
1
2
(
e
Q
κ
)
αβγ
e
(
κ
)
ν
q
e
i
q
·
R
due to
the displacement induced by a phonon is obtained as
V
quad
ν
q
(
r
) =
κ
R
V
quad
(
r
;
τ
κ
R
).
Following steps
analogous to the dipole case, we find:
V
quad
ν
q
(
r
) =
e
2
ε
0
κ
M
1
/
2
κ
G
6
=
q
1
2
(
q
+
G
)
·
(
Q
κ
e
(
κ
)
ν
q
)
·
(
q
+
G
)
(
q
+
G
)
·

·
(
q
+
G
)
×
e
i
(
q
+
G
)
·
(
r
τ
κ
)
.
(7)
The corresponding
e
-ph matrix elements due to the
quadrupole perturbation potential are:
g
quad
mnν
(
k
,
q
) =
e
2
ε
0
κ
(
~
2
ω
ν
q
M
κ
)
1
2
G
6
=
q
1
2
(
q
+
G
)
α
(
Q
κ,αβγ
e
(
κ
)
ν
q
)(
q
+
G
)
β
(
q
+
G
)
α

αβ
(
q
+
G
)
β
m
k
+
q
|
e
i
(
q
+
G
)
·
(
r
τ
κ
)
|
n
k
.
(8)
Note that in the
q
0 limit the Fr ̈ohlich
e
-ph matrix
elements are of order 1
/q
and the quadrupole matrix el-
ements of order
q
0
, thus approaching a constant value;
both quantities are non-analytic as
q
0. Octopole and
higher electronic responses in Eq. (1) lead to potentials
that vanish as
q
0 and can be grouped together into a
short-range
e
-ph interaction, commonly referred to as the
“deformation potential” in analytic
e
-ph theories [18].
B. Interpolation scheme for
e
-ph interactions
The total
e
-ph matrix elements
g
(here we omit
the band and mode indices) can be formed by adding
together the short-range part
g
S
and the dipole and
quadrupole interactions, which can be combined into a
long-range part
g
L
. Therefore,
g
=
g
S
+
g
L
=
g
S
+
g
dip
+
g
quad
.
(9)
We start from a set of
e
-ph matrix elements
g
(
k
,
q
)
computed with DFPT on a regular coarse grid of
k
-
and
q
-points [15].
The short-range part is obtained
by subtracting the long-range terms on the coarse grid,
g
S
(
k
,
q
) =
g
(
k
,
q
)
g
dip
(
k
,
q
)
g
quad
(
k
,
q
). The short-
range
e
-ph matrix elements decay rapidly in real space,
and thus are ideal for interpolation using a localized basis
set such as Wannier functions [28] or atomic orbitals [16].
After interpolating the short-ranged part [9] on fine
k
-
and
q
-point grids, we add back the long-range dipole and
quadrupole matrix elements, computed using Eqs. (5)
and (8) directly at the fine-grid
k
and
q
-points.
4
As DFPT accurately captures the long-range dipole
and quadrupole
e
-ph interactions [15], the matrix ele-
ments obtained from DFPT can be used as a benchmark
for the interpolated results. For this comparison, follow-
ing Ref. 22 we compute the gauge-invariant
e
-ph coupling
strength,
D
ν
tot
(
q
), which is proportional to the absolute
value of the
e
-ph matrix elements:
D
ν
tot
(
q
) =
2
ω
ν
q
M
uc
~
2
mn
|
g
mnν
(
k
= Γ
,
q
)
|
2
N
b
,
(10)
where
M
uc
is the mass of the unit cell and the band
indices
n
and
m
run over the
N
b
bands selected for the
analysis.
C. Computational details
We investigate the effect of the quadrupole
e
-ph inter-
action in silicon, a nonpolar semiconductor, and tetrago-
nal PbTiO
3
, a polar piezoelectric material. Calculations
on GaN are shown in the companion work [29]. The
ground state and band structure are obtained using DFT
in the local density approximation with a plane-wave ba-
sis using the
Quantum ESPRESSO
code [30]. Kinetic
energy cutoffs of 40 Ry for silicon and 76 Ry for PbTiO
3
are employed, together with scalar-relativistic norm-
conserving pseudopotentials from Pseudo Dojo [31]. The
calculations employ lattice constants of 10.102 bohr for
silicon and 7.275 bohr (with aspect ratio
c/a
= 1.046)
for PbTiO
3
. We use the dynamical quadrupole tensors
computed in Ref. 21. The phonon dispersions and
e
-
ph perturbation potentials on coarse
q
-point grids are
computed with DFPT [15]. We employ the
perturbo
code [9] to compute the
e
-ph matrix elements on coarse
Brillouin zone grids with 10
×
10
×
10
k
- and
q
-points for
silicon and 8
×
8
×
8
k
- and
q
-points for PbTiO
3
. The
Wannier functions are computed with the
Wannier90
code [28] and employed in
perturbo
[9] to interpolate
the short-range
e
-ph matrix elements.
We compute the scattering rates and electron mobility
using the
perturbo
code [9]. Briefly, the band- and
k
-
dependent
e
-ph scattering rate Γ
n
k
is obtained as
Γ
n
k
=
2
π
~
q
|
g
mnν
(
k
,
q
)
|
2
[(
N
ν
q
+ 1
f
m
k
+
q
)
δ
(
ε
n
k
ε
m
k
+
q
~
ω
ν
q
)
+ (
N
ν
q
+
f
m
k
+
q
)
δ
(
ε
n
k
ε
m
k
+
q
+
~
ω
ν
q
)]
,
(11)
where
ε
n
k
and
~
ω
ν
q
are the electron and phonon en-
ergies, respectively, and
f
n
k
and
N
ν
q
the corresponding
temperature-dependent occupations. The scattering rate
can be further divided into the long-range part [4], Γ
L
n
k
,
by replacing
|
g
|
2
in Eq. (11) with
g
L
2
. The carrier mo-
bility is computed using
μ
=
σ/
(
n
c
e
), where
σ
is the
electrical conductivity and
n
c
is the carrier concentration.
0
.
5
L
Γ
0
.
5
K
1
2
3
4
D

t
o
t
(
q
)
(
e
V
/
Å
)
S
i
l
i
c
o
n
L
A
L
A
L
O
T
O
D
F
P
T
W
a
n
n
.
o
n
l
y
W
a
n
n
.
+
q
u
a
d
.
FIG. 2. Mode-resolved
e
-ph coupling strength [see Eq. (10)]
in silicon, computed using the lowest valence band. The elec-
tron momentum
k
is fixed at the Γ point and the phonon
wave-vector
q
is varied along high-symmetry lines in the Bril-
louin zone. Benchmark results from DFPT (black circles) are
compared with Wannier interpolation with the quadrupole
e
-
ph interaction included (orange line) or neglected (blue line).
The coarse-grid
q
-points are indicated with vertical lines.
The electrical conductivity
σ
is computed within the re-
laxation time approximation of the Boltzmann transport
equation [9, 32]:
σ
αβ
=
e
2
+
−∞
dE
(
∂f/∂E
αβ
(
E,T
)
,
(12)
where Σ
αβ
(
E,T
) is the transport distribution function
at energy
E
,
Σ
αβ
(
E,T
) =
s
N
k
n
k
τ
n
k
(
T
)
v
α
n
k
v
β
n
k
δ
(
E
ε
n
k
)
,
(13)
which is computed in
perturbo
using the tetrahedron
integration method [33]. Above,
s
is the spin degeneracy,
N
k
is the number of
k
-points,
v
n
k
is the band velocity,
and
τ
n
k
= (Γ
n
k
)
1
is the relaxation time. The mobility
is computed with non-degenerate electron concentrations
of 10
15
cm
3
for silicon and 10
17
cm
3
for PbTiO
3
. To
fully converge the scattering rates and mobility, we use
e
-ph matrix elements evaluated on fine Brillouin zone
grids with 200
×
200
×
200
k
-points and 8
×
10
6
random
q
-points.
III. RESULTS
A. Quadrupole effect on the
e
-ph matrix elements
The long-range quadrupole
e
-ph interaction is present
in a wide range of semiconductors and insulators, where
the atomic dynamical quadrupoles are in general non-
zero. We illustrate this point by studying silicon, a sim-
ple nonpolar semiconductor in which the Born charges
5
and thus the Fr ̈ohlich interaction
vanish and the pres-
ence of long-range interactions is not immediately obvi-
ous. Figure 2 shows the
e
-ph coupling strength,
D
ν
tot
(
q
)
in Eq. (10), computed directly using DFPT as a bench-
mark and compared with Wannier interpolation with and
without inclusion of the quadrupole term. The DFPT
benchmark
e
-ph matrix elements for optical modes ap-
proach a constant value as
q
0, as we show for the LO
mode in the Γ
L
direction and the transverse optical
(TO) mode along Γ
K
. This trend is distinctive of the
quadrupole
e
-ph interaction, which is of order
q
0
in the
long-wavelength limit.
If the quadrupole term is neglected and all
e
-ph in-
teractions are treated as short-ranged, the
e
-ph matrix
elements for optical modes in silicon incorrectly vanish
as
q
0. The interpolated values for optical modes are
underestimated between the Γ point, where the error is
greatest, and its nearest-neighbor
q
-points in the coarse
grid, where the error vanishes. Outside this
q
-point re-
gion close to Γ, the interpolated matrix elements without
the quadrupole interaction still deviate from the DFPT
result, although the error is smaller than near Γ. When
the quadrupole term is included, the long-range
e
-ph in-
teractions for the optical modes are captured correctly,
as can be seen for the Wannier plus quadrupole curves
in Fig. 2. The root-mean-square deviation of
D
ν
tot
(
q
)
from DFPT, for the optical branches shown in Fig. 2, is
0.78 eV
/
̊
A when the quadrupole term is neglected ver-
sus 0.03 eV
/
̊
A when the quadrupole term is included in
the interpolation. This result highlights the importance
of the quadrupole term to correctly capture long-range
e
-ph interactions in nonpolar semiconductors.
Observe also how for acoustic modes in silicon the
quadrupole term has a nearly negligible effect, as we show
for the longitudinal acoustic (LA) mode in Fig. 2. As
contracting the dynamical quadrupoles
Q
κ
with a rigid
shift of the lattice leads to a vanishing quadrupole contri-
bution [18], one can obtain the quadrupole acoustic sum
rule
α
Q
κ,αβγ
= 0 for nonpolar materials [18]. This
sum rule, which is satisfied by the dynamical quadrupole
values we employ for silicon [21], leads to a negligible
quadrupole correction for acoustic modes in the long-
wavelength limit. Though we focus on silicon in this
work, on the basis of our results we expect sizable
quadrupole contributions for optical modes, and negli-
gible for acoustic modes, in all nonpolar semiconductors.
The quadrupole
e
-ph interaction is particularly critical
in piezoelectric materials, as discussed here for tetrago-
nal PbTiO
3
, a prototypical piezoelectric insulator. Piezo-
electric materials are polar noncentrosymmetric systems
with non-zero Born charges. As a result, the dipole
Fr ̈ohlich interaction is dominant for LO modes near
q
0 due to its 1
/q
divergence. The quadrupole contri-
bution is expected to be important for TO and acoustic
modes (the quadrupole acoustic sum rule does not hold
for polar noncentrosymmetric crystals).
Figure 3 shows the
e
-ph coupling strength,
D
ν
tot
(
q
) in
Eq. (10), for the DFPT benchmark in tetragonal PbTiO
3
,
1
2
3
4
D
F
P
T
D
i
p
o
l
e
o
n
l
y
D
i
p
.
+
q
u
a
d
.
D

t
o
t
(
q
)
(
e
V
/
Å
)
L
A
T
A
L
A
T
A
L
O
P
b
T
i
O
3
0
.
2
M
Γ
Γ
0
.
2
X
Γ
0
.
2
Z
L
A
L
O
FIG. 3. Mode-resolved
e
-ph coupling strength [see Eq. (10)]
in tetragonal PbTiO
3
, computed using the lowest conduction
band. The initial electron momentum is fixed at the Γ point
and the phonon wave-vector
q
is varied along high-symmetry
lines in the Brillouin zone. Benchmark results from DFPT
(black circles) are compared with Wannier interpolation plus
the Fr ̈ohlich interaction (blue line) and Wannier interpolation
plus the Fr ̈ohlich and quadrupole interactions (orange line).
and compares it with interpolated results that include
only the Fr ̈ohlich dipole interaction or both the Fr ̈ohlich
and the quadrupole interactions. The short-range inter-
actions are included through Wannier interpolation in
both cases. When only the Fr ̈ohlich dipole interaction
is included, the
e
-ph matrix elements deviate dramati-
cally from the DFPT results. The values are either over-
estimated or underestimated depending on the phonon
mode considered, with deviations from DFPT that de-
pend strongly on the direction in which
q
approaches
Γ due to the non-analytic character of the long-range
e
-ph interactions. When the quadrupole
e
-ph interac-
tion is taken into account, the interpolated
e
-ph coupling
strength matches the DFPT result very accurately for
all
phonon modes. For LO modes, the quadrupole correction
is moderate due to the dominant Fr ̈ohlich term near
q
= 0.
For other optical and acoustic modes with a finite
e
-ph
coupling at
q
= 0, the quadrupole term removes the large
error in the dipole-only results (up to an order of mag-
nitude) and gives
e
-ph matrix elements in nearly exact
agreement with DFPT. For the branches shown in Fig. 3,
the root-mean-square deviation of
D
ν
tot
(
q
) from DFPT is
0.46 eV
/
̊
A for dipole-only results versus 0.03 eV
/
̊
A for
our dipole plus quadrupole interpolation scheme. It is
clear that the quadrupole term is essential in piezoelec-
tric materials for all phonon modes.
Contrary to silicon and nonpolar materials, the
quadrupole term has a large effect for acoustic modes
in piezoelectric materials, where it is one of the two
contributions to the so-called piezoelectric
e
-ph interac-
tion [24]. Expanding the phonon eigenvectors at
q
0
as
e
ν
q
e
(0)
ν
q
+
i
q
·
e
(1)
ν
q
, one finds two contributions of
order
q
0
[18]. One is from the Born charges,
Z
κ
e
(1)
ν
q
, and
6
FIG. 4. Room temperature scattering rate versus electron
energy (referenced to the conduction band minimum) in (a)
silicon and (b) PbTiO
3
. For silicon, we plot the quadrupole
contribution multiplied by 100 (orange) and the total scat-
tering rate (black), which includes the short-range and the
quadrupole contributions. For PbTiO
3
, we show the long-
range scattering rate computed using only the Fr ̈ohlich inter-
action (blue) or both the Fr ̈ohlich and quadrupole interactions
(orange).
is a dipole-like interaction generated by atoms with a net
charge experiencing different displacements due to strain
from an acoustic mode. The other is from the dynamical
quadrupoles,
Q
κ
e
(0)
ν
q
, and is associated with a clamped-
ion electronic polarization [34]. The
ab initio
Fr ̈ohlich
interaction includes only the former term, namely the
strain component of the piezoelectric
e
-ph interaction,
and thus the dipole-only scheme leads to large errors for
acoustic phonons in PbTiO
3
(see Fig. 3) as it neglects
the important electronic quadrupole contribution. Un-
til now, the
ab initio
Fr ̈ohlich term has been mistakenly
thought to fully capture piezoelectric
e
-ph interactions.
Our results demonstrate that both dipole and quadrupole
terms are essential for accurate acoustic mode
e
-ph in-
teractions in piezoelectric materials [35]. The relative
magnitude of the strain and quadrupole contributions is
material dependent
the two terms can nearly cancel
each other out, as we have shown elsewhere for GaN [29],
or their ratio can be mode and phonon wave-vector de-
pendent, as we find in PbTiO
3
.
B. Quadrupole contribution to the scattering rate
Because the quadrupole interaction has a significant
effect on the
e
-ph matrix elements, we expect that it also
plays a role in calculations of the
e
-ph scattering rate
and mobility. Figure 4(a) shows both the quadrupole
contribution and the total
e
-ph scattering rate in silicon
at 300 K for electron energies near the conduction band
minimum. We find that the quadrupole contribution to
the scattering rate is about 1% of the total scattering
rate at temperatures between 100
400 K. At electron
energies below the optical phonon emission threshold in
silicon (
~
ω
O
65 meV relative to the conduction band
minimum), absorption and emission of acoustic phonons
dominate the scattering processes, and thus we find a
small correction due to the quadrupole interaction, which
minimally affects acoustic modes in silicon. Since the
quadrupole acoustic sum rule holds only in the long wave-
length limit, the quadrupole interaction can still con-
tribute to finite-
q
acoustic scattering, as is shown by the
fact that the quadrupole scattering rate at energy be-
low
~
ω
O
is proportional to the total scattering rate. The
quadrupole contribution increases sharply above the op-
tical emission threshold because the quadrupole term is
greater for optical modes in silicon. For the same reason,
the relative contribution of the quadrupole term increases
slightly with temperature in the 100
400 K range, vary-
ing from 1% of the total scattering rate at 100 K to 1.5%
at 400 K.
The effect of the quadrupole interaction on the scat-
tering rates is greater in PbTiO
3
. Our analysis focuses
on the
e
-ph scattering rate due to the long-range
e
-ph
interactions, although similar conclusions hold for the
total scattering rate. Figure 4(b) shows the long-range
e
-ph scattering rate in PbTiO
3
at 300 K as a function
of electron energy, comparing results that include only
the dipole Fr ̈ohlich interaction with results from our ap-
proach including both the dipole and quadrupole terms.
The scattering rate from the long-range
e
-ph interactions
is lower at all energies when the quadrupole term is taken
into account. The difference is greatest near the band
edge, where the scattering rate due to the dipole inter-
action alone is 0.075 fs
1
versus a 50% smaller value of
0.050 fs
1
for dipole plus quadrupole.
These trends can be understood on the basis of the
e
-
ph matrix element analysis in Fig. 3. The errors found
when neglecting the quadrupole term in
D
ν
tot
(
q
), which
is proportional to the absolute value of the matrix ele-
ments [see Eq. (10)], are amplified in calculations of the
scattering rate, which is proportional to the square of the
matrix elements. The largest errors we find for
D
ν
tot
(
q
)
are in the
q
0 limit, especially for the acoustic modes.
For example, for the LA mode in the Γ
M
and Γ
X
directions, the value of
D
ν
tot
(
q
) from the dipole-only cal-
culation is 0.17 eV
/
̊
A compared to a twice-greater value
of 0.40 eV
/
̊
A when the quadrupole term is included. This
leads to a four-fold increase of the LA mode scattering
rate upon including the quadrupole interaction. Oppo-
site to the silicon case, in PbTiO
3
the relative magnitude
of the quadrupole correction is greater at lower temper-
atures because the quadrupole interaction is stronger for
acoustic modes. Near the band edge, we find quadrupole
corrections to the long-range scattering rate ranging from
97% at 100 K to 38% at 400 K. Given that low-energy
electronic states near the band edge govern transport
properties, including the quadrupole term is critical to
accurately computing electronic transport.
7
1
0
0
2
0
0
3
0
0
4
0
0
0
5
0
0
0
1
0
0
0
0
1
5
0
0
0
1
0
0
2
0
0
3
0
0
4
0
0
0
5
0
1
0
0
1
5
0
2
0
0
2
5
0
M
o
b
i
l
i
t
y
(
c
m
2
/
V
s
)
T
e
m
p
e
r
a
t
u
r
e
(
K
)
W
a
n
n
i
e
r
o
n
l
y
W
a
n
n
i
e
r
+
q
u
a
d
.
S
i
l
i
c
o
n
T
e
m
p
e
r
a
t
u
r
e
(
K
)
D
i
p
o
l
e
o
n
l
y
D
i
p
o
l
e
+
q
u
a
d
.
P
b
T
i
O
3
(
a
)
(
b
)
FIG. 5. Computed temperature-dependent electron mobility
in (a) silicon and (b) tetragonal PbTiO
3
. The plot compares
the mobility obtained when the quadrupole
e
-ph interaction
is included (orange squares) or neglected (blue circles). The
PbTiO
3
results are for transport in the basal
xy
plane.
C. Quadrupole contribution to the mobility
The effect of the quadrupole
e
-ph interaction on the
mobility is noteworthy. Figure 5(a) shows the temper-
ature dependent electron mobility in silicon computed
with and without the quadrupole term. Including the
quadrupole interaction reduces the computed mobility by
approximately 5
10% due to the increased
e
-ph coupling
strength and scattering rates. For example, the com-
puted mobility at 300 K is 1390 cm
2
/
Vs when including
the quadrupole interaction versus a value of 1473 cm
2
/
Vs
with the conventional interpolation approach in which all
e
-ph interactions in silicon are treated as short-ranged.
This discrepancy is due to the underestimation of the
e
-ph coupling strength for optical modes in the conven-
tional approach, especially at small values of
q
as shown
in Fig. 2. As a result, intravalley scattering due to optical
modes is underestimated without the quadrupole term,
leading to an artificially high mobility. Although we focus
on silicon, we expect that these trends apply in general
to nonpolar semiconductors because small-
q
optical
e
-ph
coupling will consistently be underestimated without the
quadrupole term. The long-range quadrupole
e
-ph in-
teraction is thus surprisingly manifest in the transport
properties of nonpolar materials.
We find an opposite trend in PbTiO
3
, in which includ-
ing the quadrupole interaction increases the mobility by
10
25% between 100
400 K, as seen in Fig. 5(b). The
quadrupole term gives a larger correction at lower tem-
peratures, reaching values up to
25% at 100 K. This
result is due to the dominant acoustic mode contribution
at low temperatures together with the large quadrupole
correction for acoustic modes in piezoelectric materials.
At higher temperatures, where optical mode scattering
is dominant and acoustic scattering less important, the
quadrupole contribution is smaller, only about 10% at
400 K. Due to differences in the quadrupole interaction
for different phonon modes and to varying mode con-
tributions to the mobility as a function of temperature,
including the quadrupole term corrects the temperature
dependence of the mobility [29] and is essential in piezo-
electric materials.
IV. DISCUSSION
We briefly discuss a technical aspect of the
e
-ph matrix
element interpolation. The treatment of long-wavelength
perturbations with wave-vector
q
0 in DFPT is crit-
ical in semiconductors and insulators [16, 36].
The
lattice-periodic part of the phonon perturbation poten-
tial, ∆
v
q
(
r
), is the sum of a Coulomb and an exchange-
correlation contribution,
v
q
(
r
) = ∆
v
q
,C
(
r
) + ∆
v
q
,XC
(
r
)
.
(14)
The Coulomb contribution ∆
v
q
,C
(
r
) combines the vari-
ation of the Hartree and electron-nuclei interactions. Its
integral over the unit cell [36],
∆(
q
) =
1
d
r
v
q
,C
(
r
)
,
(15)
is well-behaved for insulators (and semiconductors) at fi-
nite
q
values, but is ill-defined at
q
= 0. First-principles
codes such as
Quantum ESPRESSO
[30] subtract ∆(
q
)
from the perturbation potential at
q
= 0, thus making it
discontinuous at
q
= 0. Therefore, due to both the dis-
continuity at
q
= 0 and the non-analytic behavior near
q
= 0, the
e
-ph matrix elements are challenging to inter-
polate in the long-wavelength limit.
In our scheme, we identify the quadrupole interaction
as the key long-range term in nonpolar materials, and
remove the non-analytic behavior near
q
= 0 on the
coarse grid by subtracting the quadrupole term. This
strategy improves the interpolation near
q
= 0 in non-
polar materials, at once capturing the correct physics
and smoothing the coarse-grid matrix element to be in-
terpolated. Due to the non-analytic behavior, denser
DFPT grids cannot fully remove the interpolation error
if the quadrupole term is not subtracted on the coarse
grid [37]. For polar materials such as PbTiO
3
, the non-
analytic behavior is due to both the dipole (Fr ̈ohlich) and
quadrupole long-range
e
-ph interactions. By subtracting
both terms in our scheme in polar materials, the coarse-
grid matrix elements to be interpolated are made smooth
and the interpolation approach more reliable. The non-
analytic behavior of the Coulomb potential is correctly
reconstructed by adding back the dipole (in polar ma-
terials) and quadrupole (in all insulators) contributions
after interpolation.