PHYSICAL REVIEW B
102
, 125203 (2020)
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, USA
3
Center for Computational Quantum Physics, Flatiron Institute, 162 Fifth Avenue, New York 10010, USA
(Received 30 March 2020; accepted 29 May 2020; published 21 September 2020)
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öhlich 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.
DOI:
10.1103/PhysRevB.102.125203
I. INTRODUCTION
Electron-phonon (e-ph) interactions are key to under-
standing electrical transport, nonequilibrium dynamics, and
superconductivity [
1
]. First-principles calculations can pro-
vide microscopic insight into e-ph scattering processes and
are rapidly emerging as a quantitative tool for investigating
charge transport and ultrafast carrier dynamics in materials
[
2
–
13
]. The typical workflow combines density functional
theory (DFT) [
14
] calculations of the ground state and band
structure with density functional perturbation theory (DFPT)
[
15
] for phonon dispersions and e-ph perturbation potentials.
As DFPT can compute the electronic response to periodic
lattice perturbations (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 computation-
ally demanding to be carried out on the fine Brillouin 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 Bril-
louin zone grids with order of 10
×
10
×
10
q
points, followed
by interpolation of the e-ph matrix elements with a localized
basis set such as Wannier functions or atomic orbitals [
16
].
As the perturbation potential can be nonanalytic near
q
=
0
[
16
] or even exhibit a divergence for certain phonon modes,
interpolation is particularly challenging and less reliable in
the region between
q
=
0 and its nearest-neighbor
q
points
*
Corresponding author: bmarco@caltech.edu
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 potential
shows that in the long-wavelength limit (phonon wave-vector
q
→
0) the long-range dipole Fröhlich term diverges as 1
/
q
,
the quadrupole term approaches a constant 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 interpolation 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
formal-
ism and computed with first-principles quantities such as
the atomic dynamical dipoles [
15
] and quadrupoles [
19
–
21
]
induced by lattice vibrations, which can be computed with
DFPT. For each atom
κ
, one can obtain Born charge (
Z
κ
) and
dynamical quadrupole (
Q
κ
) tensors, which, once contracted
with the phonon eigenvector, give the atomic contributions to
the dipole and quadrupole e-ph interactions.
The dipole Fröhlich 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 interactions, it is useful to consider
2469-9950/2020/102(12)/125203(8)
125203-1
©2020 American Physical Society
JINSOO PARK
et al.
PHYSICAL REVIEW B
102
, 125203 (2020)
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 different phonon modes in various materials [
24
].
In ionic and polar covalent crystals (here and below, de-
noted as polar materials), the dipole Fröhlich term is dominant
as
q
→
0 due to its 1
/
q
trend. This e-ph interaction 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 dominant long-range e-ph inter-
action is the quadrupole term, which is particularly important
for acoustic phonons in piezoelectric (polar noncentrosym-
metric) materials. In nonpolar semiconductors such as sili-
con and germanium, the dipole Fröhlich interaction vanishes
and the quadrupole term is the dominant long-range e-ph
interaction for all modes. The quadrupole e-ph interaction is
thus expected to play an important role in many classes of
materials, making a compelling case for its inclusion 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 elements 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 piezo-
electric material, the quadrupole corrections are substantial
for all phonon modes and are particularly important for
acoustic modes, which contribute to the piezoelectric e-ph
interaction. Including only the long-range Fröhlich interaction
and neglecting the quadrupole term leads to large errors in
PbTiO
3
, while adding the quadrupole term leads to e-ph ma-
trix elements that accurately reproduce the DFPT benchmark
results for all phonon modes in the entire Brillouin zone.
We investigate the impact of the quadrupole e-ph interaction
on the electron scattering rates and mobility 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 cap-
ture the long-range e-ph interactions. In turn, this development
enables more precise calculations of electron dynamics and
scattering processes from first principles.
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)
FIG. 1. Schematic of the dipole and quadrupole charge configu-
rations giving rise to long-range e-ph interactions.
where summation over the Cartesian indices
β
and
γ
is
implied. This polarization response defines the Born effec-
tive 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 corre-
sponding long-range e-ph interactions 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 elements [
2
]
g
mn
ν
(
k
,
q
)
=
(
̄
h
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 momentum
k
to scatter into a final state
|
m
k
+
q
by emitting or absorbing
a phonon with energy ̄
h
ω
ν
q
.
A. Dipole and quadrupole e-ph interactions
To derive the dipole and quadrupole perturbation poten-
tials, we consider a Born–von Karman (BvK) crystal [
26
] with
N
unit cells and volume
N
. The potential due to a dipole
configuration with dipole moment
p
centered at position
τ
in
the crystal (see Fig.
1
) can be written as [
22
,
23
]
V
dip
(
r
;
τ
)
=
i
e
N
ε
0
∑
q
∑
G
=−
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 displacement induced by a
phonon with mode index
ν
and wave-vector
q
. The resulting
atomic dynamical dipole is
p
κ,
R
=
(
e
Z
κ
)
̃
e
(
κ
)
ν
q
e
i
q
·
R
, where the
phonon eigenvector projected on atom
κ
is defined as
̃
e
(
κ
)
ν
q
=
e
(
κ
)
ν
q
/
√
M
κ
, with
e
ν
q
the eigenvector of the dynamical matrix at
q
and
M
κ
themassofatom
κ
. Summing over the contributions
125203-2
LONG-RANGE QUADRUPOLE ELECTRON-PHONON ...
PHYSICAL REVIEW B
102
, 125203 (2020)
from all atoms
κ
at lattice vectors
R
with positions
τ
κ
R
=
τ
κ
+
R
in the BvK supercell, the total 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
=−
q
×
(
Z
κ
e
(
κ
)
ν
q
)
·
(
q
+
G
)
e
i
(
q
+
G
)
·
(
r
−
τ
)
(
q
+
G
)
·
·
(
q
+
G
)
.
(4)
The
ab initio
Fröhlich e-ph coupling is obtained by evaluating
the matrix elements with this potential:
g
dip
mn
ν
(
k
,
q
)
=
i
e
2
ε
0
∑
κ
(
̄
h
2
ω
ν
q
M
κ
)
1
/
2
∑
G
=−
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 response
can be derived with a similar strategy. We first consider
the potential generated by a quadrupole charge configuration
consisting of two equal and oppositely 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
=−
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
;
τ
)inEq.(
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 dis-
placement 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
=−
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
∑
κ
(
̄
h
2
ω
ν
q
M
κ
)
1
/
2
∑
G
=−
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öhlich e-ph matrix elements
are of order 1
/
q
and the quadrupole matrix elements of
order
q
0
, thus approaching a constant value; both quantities
are nonanalytic 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 interac-
tion, 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
) com-
puted 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.
As DFPT accurately captures the long-range dipole and
quadrupole e-ph interactions [
15
], the matrix elements ob-
tained from DFPT can be used as a benchmark for the inter-
polated results. For this comparison, following 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
̄
h
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 interaction
in silicon, a nonpolar semiconductor, and tetragonal PbTiO
3
,
a polar piezoelectric material. Calculations on GaN are shown
in the companion work [
29
]. The ground state and band struc-
ture are obtained using DFT in the local density approxima-
tion with a plane-wave basis using the Q
UANTUM
E
SPRESSO
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
]. We have verified that spin-orbit coupling has
a negligible effect. The calculations employ relaxed lattice
constants of 10.102 bohrs for silicon and 7.275 bohrs (with
aspect ratio
c
/
a
=
1
.
046) for PbTiO
3
. We use the dynamical
quadrupole tensors computed in Refs. [
21
,
32
]. The phonon
dispersions and e-ph perturbation potentials on coarse
q
-
point grids are computed with DFPT [
15
]. We employ the
125203-3
JINSOO PARK
et al.
PHYSICAL REVIEW B
102
, 125203 (2020)
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 W
ANNIER
90
code [
28
] and employed in
PERTURBO
[
9
] to interpolate the
short-range e-ph matrix elements.
We compute the scattering rates and electron mobility us-
ing the
PERTURBO
code [
9
]. Briefly, the band- and
k
-dependent
e-ph scattering rate
n
k
is obtained as
n
k
=
2
π
̄
h
∑
m
ν
q
|
g
mn
ν
(
k
,
q
)
|
2
×
[(
N
ν
q
+
1
−
f
m
k
+
q
)
δ
(
ε
n
k
−
ε
m
k
+
q
−
̄
h
ω
ν
q
)
+
(
N
ν
q
+
f
m
k
+
q
)
δ
(
ε
n
k
−
ε
m
k
+
q
+
̄
h
ω
ν
q
)]
,
(11)
where
ε
n
k
and ̄
h
ω
ν
q
are the electron and phonon energies, re-
spectively, and
f
n
k
and
N
ν
q
are 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 mobility is computed using
μ
=
σ/
(
n
c
e
), where
σ
is the electrical conductivity and
n
c
is the carrier concentration. The electrical conductivity
σ
is
computed within the relaxation time approximation of the
Boltzmann transport equation [
9
,
33
]:
σ
αβ
=
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 integra-
tion method [
34
]. 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 [
35
].
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 nonzero. We
illustrate this point by studying silicon, a simple nonpolar
semiconductor in which the Born charges—and thus the
Fröhlich interaction—vanish and the presence of long-range
interactions is not immediately obvious. Figure
2
shows the
e-ph coupling strength
D
ν
tot
(
q
)inEq.(
10
), computed directly
using DFPT as a benchmark and compared with Wannier
interpolation with and without inclusion of the quadrupole
term. The DFPT benchmark e-ph matrix elements for optical
modes approach 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
FIG. 2. Mode-resolved e-ph coupling strength [see Eq. (
10
)]
in silicon, computed using the lowest valence band. The electron
momentum
k
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 inter-
polation with the quadrupole e-ph interaction included (orange line)
or neglected (blue line). The coarse-grid
q
points are indicated with
vertical lines.
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 interactions
are treated as short ranged, the e-ph matrix elements for opti-
cal modes in silicon incorrectly vanish as
q
→
0. The interpo-
lated 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 region 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 interactions
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
,is0.78eV
/
Å when the quadrupole
term is neglected versus 0.03 eV
/
Å when the quadrupole term
is included in the interpolation. This result highlights the im-
portance of the quadrupole term to correctly describe 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 contribution [
18
], one can
obtain the quadrupole acoustic sum rule
∑
α
Q
κ,αβγ
=
0for
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 negligible for
acoustic modes, in all nonpolar semiconductors.
125203-4