11footnotetext: See Supplemental Material, which includes Refs. Poncé  _et
al._ (2021); Adolph  _et al._ (1996); Lopez  _et al._ (2012), at [URL will be
inserted by publisher] for the details of the calculation of the velocity and
related matrix elements and the convergence study.

# Wannier Function Perturbation Theory:
Localized Representation and Interpolation of Wavefunction Perturbation

Jae-Mo Lihm<EMAIL_ADDRESS>Cheol-Hwan Park<EMAIL_ADDRESS>Center
for Correlated Electron Systems, Institute for Basic Science, Seoul 08826,
Korea  Department of Physics and Astronomy, Seoul National University, Seoul
08826, Korea  Center for Theoretical Physics, Seoul National University, Seoul
08826, Korea

###### Abstract

Thanks to the nearsightedness principle, the low-energy electronic structure
of solids can be represented by localized states such as the Wannier
functions. Wannier functions are actively being applied to a wide range of
phenomena in condensed matter systems. However, the Wannier-function-based
representation is limited to a small number of bands and thus cannot describe
the change of wavefunctions due to various kinds of perturbations, which
require sums over an infinite number of bands. Here, we introduce the concept
of the Wannier function perturbation, which provides a localized
representation of wavefunction perturbations. Wannier function perturbation
theory allows efficient calculation of numerous quantities involving
wavefunction perturbation, among which we provide three applications. First,
we calculate the temperature-dependent indirect optical absorption spectra of
silicon near the absorption edge nonadiabatically, i.e., differentiating
phonon-absorption and phonon-emission processes, and without arbitrary
temperature-dependent shifts in energy. Second, we establish a theory to
calculate the shift spin conductivity without any band-truncation error.
Unlike the shift charge conductivity, an exact calculation of the shift spin
conductivity is not possible within the conventional Wannier function methods
because it cannot be obtained from geometric quantities for low-energy bands.
We apply the theory to monolayer WTe2. Third, we calculate the spin Hall
conductivity of the same material again without any band-truncation error.
Wannier function perturbation theory is a versatile method that can be readily
applied to calculate a wide range of quantities related to various kinds of
perturbations.

##  I Introduction

Most physical and chemical systems are “nearsighted” Kohn (1996); Prodan and
Kohn (2005): Their properties can be studied by inspecting one local region at
a time and the long-range effects can be easily approximated. A localized
representation of a system is highly desirable because it enables one to
exploit the nearsightedness to its full extent. Wannier functions (WFs)
provide a localized representation of wavefunctions Wannier (1937); Kohn
(1959). The invention of the maximally localized Wannier function (MLWF)
method Marzari and Vanderbilt (1997); Souza  _et al._ (2001) enabled the
calculation of the localized WFs for electronic structures computed using
density functional theory (DFT). Using the MLWFs, the low-energy eigenvalues
and eigenstates can be interpolated very efficiently.

Applications of MLWFs are stretched over all fields of condensed matter
physics and materials science Marzari  _et al._ (2012). MLWFs are used to
efficiently compute electron-phonon coupling (EPC) at a dense grid of electron
and phonon wavevectors Giustino  _et al._ (2007). Wannier interpolation of EPC
is actively being applied to study transport Zhou and Bernardi (2016); Poncé
_et al._ (2018), phonon-assisted optical absorption Noffsinger  _et al._
(2012), superconductivity Margine and Giustino (2013), and ultrafast carrier
dynamics Jhalani  _et al._ (2017). Various linear and nonlinear responses of
solids to dc and ac electric fields are also routinely calculated using
Wannier interpolation Wang  _et al._ (2006); Yates  _et al._ (2007); Ibañez-
Azpiroz  _et al._ (2018); Tsirkin  _et al._ (2018); Qiao  _et al._ (2018);
Ryoo  _et al._ (2019).

However, one of the most severe limitations of WF interpolation is the limited
number of bands being interpolated. Consequently, WFs cannot represent the
wavefunction perturbation, which is the sum of contributions from an infinite
number of bands. Obtaining MLWFs for the high-energy, plane-wave-like bands is
a highly nontrivial problem. Even assuming that the MLWFs are obtained,
calculations in the MLWF basis will be much more costly, negating the
efficiency of Wannier interpolation. In DFT calculations, the wavefunction
perturbation can be calculated without referring to the high-energy states by
solving a linear equation known as the Sternheimer equation Sternheimer
(1954). This method forms the basis of the calculation of the linear response
properties of solids using density functional perturbation theory (DFPT)
Baroni  _et al._ (2001). However, the Sternheimer equation is not directly
applicable to WFs since the WFs do not accurately represent the wavefunction
perturbation. Therefore, WF methods cannot deal with properties that require
the calculation of wavefunction perturbations. Ge et al. Ge and Lu (2015)
developed a local representation of wavefunction perturbation for the
dielectric response of insulators. However, this formalism cannot be applied
to WFs for entangled bands Souza  _et al._ (2001), while such generalization
is necessary for many applications.

A representative example where the interpolation of wavefunction perturbations
is needed is the real part of the phonon-induced electron self-energy. This
self-energy causes the temperature dependence and zero-point renormalization
of electron band structures Giustino  _et al._ (2010a); Gonze  _et al._
(2011). The state-of-the-art method for the perturbative calculation of the
real part of the self-energy is to solve the Sternheimer equation using a
complete basis set such as the plane waves Gonze  _et al._ (2011). However,
such a calculation is too costly to be done on a large number of $k$ points.
The high computational cost prohibits the consideration of phonon-induced
renormalization at a dense $k$-point mesh in the study of optical and
transport properties of materials. Also, band disentanglement Souza  _et al._
(2001) is necessary even for the valence band of insulators because the
wavefunctions and energy eigenvalues of the low-lying conduction bands are
required.

Another example of a quantity involving the wavefunction perturbation is the
shift current Sipe and Shkrebtii (2000). For the shift charge current, the
wavefunction perturbation can alternatively be expressed as a derivative with
respect to the crystal momentum ${\bm{\mathbf{k}}}$, allowing an exact WF-
based calculation Ibañez-Azpiroz  _et al._ (2018). However, such an
alternative expression cannot be made for the shift spin current in spin-orbit
coupled systems Lihm and Park . Previous calculations of shift spin current in
spin-orbit coupled systems truncated the infinite sum over bands Xu  _et al._
(2021); Mu  _et al._ (2021) or used the tight-binding approximation Lihm and
Park . Similarly, conventional Wannier interpolation schemes for the spin Hall
conductivity Qiao  _et al._ (2018); Ryoo  _et al._ (2019) resorts to the band
truncation. To our knowledge, the effect of band truncation on the shift spin
or spin Hall conductivity has not been investigated.

In this work, we develop Wannier function perturbation theory (WFPT), which
enables a localized representation of wavefunction perturbation, alleviating
the conceptual and practical limitations of the WF method. We define the
Wannier function perturbation (WFP) as the change of the WF due to a
perturbation. The WFPs are spatially localized and thus can be used to
efficiently interpolate the perturbation of eigenstate wavefunctions and the
related matrix elements. We provide three applications of WFPT. The first
application is the temperature-dependent indirect optical absorption of
silicon. We calculate the phonon renormalization of the electron band
structure using WFPT. Using the renormalized band structures, we perform a
predictive calculation of the temperature dependence of the phonon-assisted
optical absorption spectra near the absorption edge. Our calculation takes
into account nonadiabatic effects and does not require arbitrary shifts in
energy depending on the temperature. The second and third applications are the
shift spin and spin Hall conductivity of monolayer WTe2, respectively. In
these applications, we show that WFPT is necessary to use an ultradense
$k$-point grid and eliminate the band-truncation error at the same time.

##  II Formalism of Wannier function perturbation theory

Before discussing the formalism of WFPT in detail, let us briefly summarize
the key ideas. First, we Wannierize both the unperturbed and perturbed
systems. We call the corresponding WFs the unperturbed WFs and the perturbed
WFs, respectively. To study the perturbative properties of the WFs, the
perturbed WFs must converge smoothly to the unperturbed WFs in the limit of a
small perturbation. To achieve this property, we use the unperturbed WFs as
the initial guesses and construct projection-only WFs Souza  _et al._ (2001)
of the perturbed system. The resulting first-order derivatives of the WFs with
respect to the perturbation strength are the WFPs. Next, we show that the
perturbation of the eigenstate wavefunction can be efficiently interpolated
using WFPs. Finally, we study the interpolation of matrix elements that
involve wavefunction perturbations.

###  II.1 Wannier function perturbation

In this subsection, we first introduce the setting for our theory. Then, we
derive the expression for the WFs of the perturbed system and the WFPs.
Finally, we argue that the WFPs are spatially localized.

####  II.1.1 Setting

The WFs of a perfect crystal can be written as

$\ket{w^{(\mathrm{0})}_{i\mathbf{R}}}=\frac{1}{\sqrt{N_{k}}}\sum_{n,{\bm{\mathbf{k}}}}\ket{\psi^{(\mathrm{0})}_{n\mathbf{k}}}U^{(\mathrm{0})}_{ni;{\bm{\mathbf{k}}}}{e^{-i{{\bm{\mathbf{k}}}}\cdot{{\bm{\mathbf{R}}}}}}.$
(1)

Here, $\ket{\psi^{(\mathrm{0})}_{n\mathbf{k}}}$ is the unperturbed electron
eigenstate for band $n$ and crystal momentum ${\bm{\mathbf{k}}}$,
$i=1,\cdots,{N_{\mathrm{W}}}$ the WF index, ${\bm{\mathbf{R}}}$ the unit cell
position, and $N_{k}$ the size of the Born-von Karman supercell. The semi-
unitary gauge transformation matrix $U^{(\mathrm{0})}_{ni;{\bm{\mathbf{k}}}}$
is chosen to make $\ket{w^{(\mathrm{0})}_{i\mathbf{R}}}$ spatially localized.
The unperturbed eigenstates satisfy

$\innerproduct{\psi^{(\mathrm{0})}_{m\mathbf{k}}}{\psi^{(\mathrm{0})}_{n\mathbf{k^{\prime}}}}=\delta_{mn}\delta_{{\bm{\mathbf{k}}}{\bm{\mathbf{k^{\prime}}}}}.$
(2)

We apply a perturbation $\lambda\hat{V}$, where $\lambda$ is a small
perturbation parameter. The perturbation $\hat{V}$ may or may not break the
periodicity of the lattice but should respect the Born-von Karman periodic
boundary condition. The only additional condition we impose on the
perturbation $\hat{V}$ is that
$\matrixelement{{\bm{\mathbf{r}}}}{\hat{V}}{{\bm{\mathbf{r^{\prime}}}}}$
should be small if
$\absolutevalue{{\bm{\mathbf{r}}}-{\bm{\mathbf{r^{\prime}}}}}$ is large. This
condition is necessary for the perturbed system to be nearsighted. Any local
potentials, phonon perturbation potentials that contain contributions from
nonlocal pseudopotentials, velocity operator
$\hat{v}^{a}=[\hat{r}^{a},\hat{H}^{(0)}]/i\hbar$, where $\hat{H}^{(0)}$ is the
unperturbed Hamiltonian, and spin-velocity operator
$\hat{j}^{s;a}=\frac{1}{2}\anticommutator{\hat{S}_{s}}{\hat{v}^{a}}$ all
fulfill this condition. Here, $\hat{S}_{s}$ is the spin operator for spin
polarization along $s$.

Eigenstates of the perturbed system can be expanded as

$\ket{\psi_{n\mathbf{k}}}=\ket{\psi^{(\mathrm{0})}_{n\mathbf{k}}}+\lambda\ket{\psi^{(\mathrm{1})}_{n\mathbf{k}}}+{\mathcal{O}}(\lambda^{2}),$
(3)

where the wavefunction perturbation is

$\ket{\psi^{(\mathrm{1})}_{n\mathbf{k}}}=\sideset{}{{}^{\prime}}{\sum}_{{{n^{\prime}\mathbf{k^{\prime}}}}}\ket{\psi^{(0)}_{n^{\prime}\mathbf{k^{\prime}}}}\frac{\matrixelement{\psi^{(0)}_{n^{\prime}\mathbf{k^{\prime}}}}{\hat{V}}{\psi^{(0)}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}}.$
(4)

Hereafter, the primed sum denotes that the terms making the denominator zero
are excluded. For degenerate states, the true perturbed eigenstates in the
$\lambda\rightarrow 0$ limit will, in general, converge to the unperturbed
eigenstates that diagonalize $\hat{V}$ in the degenerate subspace, not to
$\ket{\psi^{(0)}_{n\mathbf{k}}}$. In this case, we linearly combine the
perturbed eigenstates so that they are no longer eigenstates of the perturbed
Hamiltonian but satisfy Eqs. (3, 4). This operation is allowed because our
goal is to construct WFs of the perturbed system, and there is gauge freedom
in constructing the WFs that allows such a linear combination. We note that
this treatment is needed only to make the derivation more concise, while, in
practice, one does not have any issue with degeneracy because of the careful
use of energy windows that will be described in the following.

####  II.1.2 WFs of the perturbed system and the WFPs

We aim to accurately describe the wavefunction perturbation of the low-energy
states whose eigenvalues are inside the “perturbation window”
${\mathcal{W}_{\mathrm{P}}}$. We denote the inner (frozen) and outer
(disentanglement) windows for the unperturbed system by
${\mathcal{W}_{\mathrm{F}}}$ and ${\mathcal{W}_{\mathrm{D}}}$, respectively
Souza  _et al._ (2001). The following inclusion relation should hold:

${\mathcal{W}_{\mathrm{P}}}\subseteq{\mathcal{W}_{\mathrm{F}}}\subseteq{\mathcal{W}_{\mathrm{D}}}$
(5)

This relation is illustrated in Fig. 1. In principle,
${\mathcal{W}_{\mathrm{P}}}$ can be chosen to be identical to
${\mathcal{W}_{\mathrm{F}}}$. Using ${\mathcal{W}_{\mathrm{P}}}$ slightly
narrower than ${\mathcal{W}_{\mathrm{F}}}$ is useful when the perturbation
mixes the states right above and below the boundary of the inner window for
the unperturbed system. For later use, we define the projection operators onto
the subspaces spanned by the unperturbed eigenstates inside
${\mathcal{W}_{\mathrm{P}}}$ and ${\mathcal{W}_{\mathrm{D}}}$ as
${\hat{P}^{(0)}_{\mathrm{P}}}$ and ${\hat{P}^{(0)}_{\mathrm{D}}}$,
respectively, and we also define
${\hat{Q}^{(0)}_{\mathrm{D}}}=\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{D}}}$.
Additionally, $N_{X}$ and $N_{X,{\bm{\mathbf{k}}}}$ denote the number of
unperturbed states inside ${\mathcal{W}}_{X}$, and those with crystal momentum
${\bm{\mathbf{k}}}$, respectively, for $X=\mathrm{P},\mathrm{F},\mathrm{D}$.
We only consider small $\lambda$ such that no eigenvalues cross the boundary
of the windows as $\lambda$ increases from 0. Another important subspace is
the unperturbed Wannier subspace, of which the projection operator is

${\hat{P}^{(0)}_{\mathrm{W}}}=\sum_{{\bm{\mathbf{R}}}}\sum_{i=1}^{{N_{\mathrm{W}}}}\ket{w_{i\mathbf{R}}^{(\mathrm{0})}}\bra{w_{i\mathbf{R}}^{(\mathrm{0})}}.$
(6)

By construction,
${\hat{P}^{(0)}_{\mathrm{P}}}\subseteq{\hat{P}^{(0)}_{\mathrm{W}}}\subseteq{\hat{P}^{(0)}_{\mathrm{D}}}$
holds, as illustrated in Fig. 1. We define “subspaces 1, 2, 3, and 4” as shown
in Fig. 1.

Figure 1:  Schematic illustration of the subspaces and energy windows. The
vertical lines indicate the separation of the full Hilbert space into
subspaces. Solid vertical lines separating different regions indicate the
separation by energy eigenvalues. The dashed vertical line separates the
subspace spanned by the unperturbed WFs from its complement and does not
indicate an energetic separation. The “high-energy” states may also include
deep semicore states which are excluded in the Wannierization process.

In general, the perturbed WFs need not make any reference to the unperturbed
WFs. However, to study the perturbative property of the WFs, the WFs must
change smoothly around $\lambda=0$. Our key idea to make the WFs smooth
functions of $\lambda$ is to use the unperturbed WFs as the initial guesses
for the projection-only WFs of the perturbed system. For Wannierization of the
perturbed system, the inner and outer windows are set to ${\mathcal{W}}_{\rm
P}$ and ${\mathcal{W}}_{\rm D}$, respectively. Using the perturbation window
as the inner window is natural because only the states inside the perturbation
window need to be accurately represented. The outer window is set to make the
conditions for Wannierization of the perturbed system as similar as possible
to that of the unperturbed system. This condition guarantees that the
perturbed WFs are close to the unperturbed WFs; hence, it makes the WFPs well
localized. See Appendix B for the numerical verification that this choice is
indeed beneficial.

Now, we explicitly construct the perturbed WFs. First, we select the subspace
to be Wannierized following the initialization recipe described in Sec. III.G
of Ref. Souza  _et al._ (2001). We take the perturbed eigenstates inside the
inner window ${\mathcal{W}_{\mathrm{P}}}$ directly. For the states not
included in ${\mathcal{W}_{\mathrm{P}}}$, we project
${\hat{P}^{(0)}_{\mathrm{W}}}$, the projection operator for the unperturbed
Wannier subspace, onto the subspace
${\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}$ and diagonalize it.
Since we are Wannierizing a supercell, we need
$N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}$ WFs. Hence, we select
$N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}-N_{\rm P}$ states inside
${\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}$ with the largest
eigenvalues Souza  _et al._ (2001) together with $N_{\rm P}$ states inside
${\mathcal{W}_{\mathrm{P}}}$ and use them as the basis states for the
perturbed WFs.

After we obtain the disentangled subspace, we compute the perturbed WFs.
Following the scheme of Ref. Marzari and Vanderbilt (1997), we project the
initial guess functions, the unperturbed WFs, onto the disentangled subspace
and perform Löwdin orthonormalization. The resulting WFs of the perturbed
system are (see Appendix A for the detailed derivation)

$\displaystyle\ket{w_{i\mathbf{R}}(\lambda)}=\ket{w^{(0)}_{i\mathbf{R}}}+\lambda\ket{w^{(1)}_{i\mathbf{R}}}+{\mathcal{O}}(\lambda^{2}),$
(7)

where

$\displaystyle\ket{w^{(1)}_{i\mathbf{R}}}=$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{m,{\bm{\mathbf{k}}}}{\hat{Q}^{(0)}_{\mathrm{D}}}\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi,{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(8) $\displaystyle+$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}\left({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}.$

We call $\ket{w^{(1)}_{i\mathbf{R}}}$, the first-order derivative of a WF at
$\lambda=0$, the “Wannier function perturbation.” The concept of WFP and its
explicit form [Eq. (8)] is the first main theoretical result of this paper.

The WFPs can be computed without explicitly referring to the states outside
the outer window. The first term of Eq. (8) can be evaluated by solving the
Sternheimer equation to obtain the wavefunction perturbation for states
${m\mathbf{k}}$ inside ${\mathcal{W}_{\mathrm{D}}}$. The second term of Eq.
(8) can be rewritten using

$\displaystyle\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}$
$\displaystyle\left({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}=\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}\sum_{m^{\prime}}\sideset{}{{}^{\prime}}{\sum}_{{n^{\prime},{\bm{\mathbf{k^{\prime}}}}}}\ket{\psi^{(0)}_{m^{\prime}\mathbf{k^{\prime}}}}$
$\displaystyle\times\left(P^{(\mathrm{0})}_{\rm D}-P^{(\mathrm{0})}_{\rm
W}\right)_{m^{\prime}n^{\prime};{\bm{\mathbf{k^{\prime}}}}}\frac{\matrixelement{\psi^{(0)}_{n^{\prime}\mathbf{k^{\prime}}}}{\hat{V}}{\psi^{(0)}_{m\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}}.$
(9)

Evaluating this expression requires only the eigenenergies and matrix elements
of $\hat{V}$ and ${\hat{P}^{(0)}_{\mathrm{W}}}$ for states inside
${\mathcal{W}_{\mathrm{D}}}$.

Figure 2(a) schematically illustrates the definition of the WFP as the change
of the WF due to a perturbation. Figures 2(b) and 2(c) show the actual WF of
silicon and the corresponding WFP for an optical phonon perturbation,
respectively. They visually demonstrate that the WFP is the change of the WF
due to the perturbation and that the WFP is spatially localized.

The formula for the WFPs [Eq. (8)] can be intuitively understood from the
partition of subspaces shown in Fig. 1. First, consider subspaces 1 and 2,
whose union is the unperturbed Wannier subspace. The wavefunction perturbation
in this subspace is not included in the WFPs because it can be reconstructed
from the perturbation theory of the WF-based tight-binding Hamiltonian.
Second, consider subspace 3, which includes the states inside the outer window
but not selected for Wannierization. To accurately describe the wavefunction
perturbation, these states should be included in the WFPs. However, the
contribution of these states to the perturbation of eigenstates in subspace 2
might diverge because the separation of subspaces 2 and 3 is done by
Wannierization not by the energy eigenvalues. In Eq. (8), this problem is
absent as the contribution of subspace 3 is included only if
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$. From Eq. (II.1.2), one can find
that this condition ensures that the denominator
$\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}$
is not zero because states ${m\mathbf{k}}$ and
${n^{\prime}\mathbf{k^{\prime}}}$ are inside subspaces 1 and 3, respectively.
The last subspace, subspace 4, is energetically separated from other states by
the outer window ${\mathcal{W}_{\mathrm{D}}}$. Its contribution to the
wavefunction perturbation is included for all states as shown in the first
term of Eq. (8). The separation of energy ranges by the windows is the key
idea in WFPT which makes the WFPs spatially localized.

####  II.1.3 Spatial localization of the WFPs

The most important property of the WFPs is their spatial localization. We
provide two arguments for the localization of the WFPs. First, the
localization of WFPs follows from the localization of the perturbed WFs. There
exists much empirical evidence that the WFs are well localized in real space
if a proper initial guess is used He and Vanderbilt (2001); Yates  _et al._
(2007); Marzari  _et al._ (2012). The initial guesses we used for the WFs of
the perturbed system are the MLWFs of the unperturbed system, which is an
ideal guess in the $\lambda\rightarrow 0$ limit. Thus, it is very likely that
the perturbed WFs are also localized, at least for small $\lambda$. Since the
WFPs are the derivative of the WFs at $\lambda=0$, the localization of the WFs
for a finite range of $\lambda$ including $\lambda=0$ implies the localization
of the WFPs.

The second argument is that the summands of Eq. (8) are almost smooth
functions of ${\bm{\mathbf{k}}}$. The summand being smooth in reciprocal space
implies localization in the real space. We emphasize that the condition
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$ in the second term of Eq. (8)
plays a crucial role. If this condition is absent, the denominator
$\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}$
in Eq. (II.1.2) can be zero, leading to a divergence in the summand and
rendering the WFPs delocalized. In our construction, the denominator never
becomes zero because $\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}$ is forced to
be inside ${\mathcal{W}_{\mathrm{P}}}$ (subspace 1), while
$\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}$ is forced to be
in subspace 3.

Although these two arguments are not rigorous in a mathematical sense, they
provide a solid justification that the WFPs are spatially localized. We note
that the localization of disentangled MLWFs is also an empirical fact without
a mathematical proof Marzari  _et al._ (2012); Cornean  _et al._ (2019).
Later, we numerically demonstrate that the WFPs and the matrix elements
involving WFPs are localized.

One source of singularity in Eq. (8) is the discontinuity that occurs when
$\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}$ crosses
${\mathcal{W}_{\mathrm{P}}}$. In this case, the summand may change
discontinuously because of the condition
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$, although it does not diverge.
This discontinuity can be removed by smoothing the inner window, as described
in Appendix C. However, for the systems studied in this work, we find that the
WFPs are already sufficiently localized without this additional step. Hence,
we use the WFPs as defined in Eq. (8).

We also note that if two bands lie close to the boundary of the outer window
from inside and outside, respectively, at some $k$ point, the first term of
Eq. (8) can give a large contribution rapidly changing with
${\bm{\mathbf{k}}}$. In such a rare case, one can locally, i.e., at that
particular $k$ point, extend the outer window to include the band outside the
original outer window to make the summand change slowly. The local extension
of the outer window can be done just for the construction of WFPs without
affecting the unperturbed WFs. In the applications studied in this paper, we
did not perform this modification as the WFPs were already sufficiently
localized.

Figure 2:  (a) Schematic illustration of the definition of a WFP. (b)
Isosurface of a WF of silicon. Isovalues are $+0.21$ and $-0.21$ times the
maximum absolute value for the yellow and cyan isosurfaces, respectively. (c)
Isosurface of a WFP for an optical phonon perturbation whose vibration mode is
indicated by the arrows. Isovalues are $+0.13$ and $-0.13$ times the maximum
absolute value for the yellow and cyan isosurfaces, respectively.

###  II.2 Interpolation of wavefunction perturbation

Next, we study how to interpolate the wavefunction perturbation of Bloch
eigenstates using the WFPs. In this subsection, we first construct the
perturbed tight-binding model using the matrix elements of the perturbed
Hamiltonian in the perturbed WF basis. Then, we use first-order perturbation
theory to obtain the eigenvector of this tight-binding model. Finally, we
calculate the wavefunction perturbation.

####  II.2.1 Perturbed tight-binding Hamiltonian

Let us calculate the Hamiltonian matrix elements of the perturbed system in
the perturbed WF basis. The Hamiltonian matrix element is

$\displaystyle H_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}=$
$\displaystyle\matrixelement{w_{i^{\prime}\mathbf{R^{\prime}}}}{(\hat{H}^{(\mathrm{0})}+\lambda\hat{V})}{w_{i\mathbf{R}}}$
$\displaystyle=$ $\displaystyle
H^{(\mathrm{0})}_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}+\lambda
H^{(\mathrm{1})}_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}+{\mathcal{O}}(\lambda^{2}),$
(10)

where

$H^{(\mathrm{0})}_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}=\matrixelement{w^{(\mathrm{0})}_{i^{\prime}\mathbf{R^{\prime}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{0})}_{{i\mathbf{R}}}}$
(11)

is the unperturbed Hamiltonian matrix element and

$\displaystyle
H^{(\mathrm{1})}_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}=$
$\displaystyle\matrixelement{w^{(\mathrm{0})}_{i^{\prime}\mathbf{R^{\prime}}}}{\hat{V}}{w^{(\mathrm{0})}_{{i\mathbf{R}}}}$
(12) $\displaystyle+$
$\displaystyle\matrixelement{w^{(\mathrm{1})}_{i^{\prime}\mathbf{R^{\prime}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{0})}_{{i\mathbf{R}}}}+\matrixelement{w^{(\mathrm{0})}_{i^{\prime}\mathbf{R^{\prime}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{1})}_{{i\mathbf{R}}}}$

is the perturbative correction. Note that
$H^{(\mathrm{1})}_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}$ contains a
direct contribution from $\hat{V}$ as well as indirect contributions from the
change of the WFs.

We write the perturbation $\hat{V}$ as a sum of monochromatic perturbations

$\hat{V}=\sum_{\bm{\mathbf{q}}}\hat{V}_{\bm{\mathbf{q}}},$ (13)

where each monochromatic perturbation satisfies

$\matrixelement{{\bm{\mathbf{r+R}}}}{\hat{V}_{\bm{\mathbf{q}}}}{{\bm{\mathbf{r^{\prime}+R}}}}={e^{i{{\bm{\mathbf{q}}}}\cdot{{\bm{\mathbf{R}}}}}}\matrixelement{{\bm{\mathbf{r}}}}{\hat{V}_{\bm{\mathbf{q}}}}{{\bm{\mathbf{r^{\prime}}}}}$
(14)

for every lattice vector ${\bm{\mathbf{R}}}$. Correspondingly, the WFPs can be
written as

$\ket{w^{(\mathrm{1})}_{i\mathbf{R}}}=\sum_{\bm{\mathbf{q}}}\ket{w^{(\mathrm{1})}_{{i\mathbf{R}};{\bm{\mathbf{q}}}}},$
(15)

where

$\innerproduct{{\bm{\mathbf{r+R^{\prime}}}}}{w^{(\mathrm{1})}_{i{\bm{\mathbf{R+R^{\prime}}}};{\bm{\mathbf{q}}}}}={e^{i{{\bm{\mathbf{q}}}}\cdot{{\bm{\mathbf{R^{\prime}}}}}}}\innerproduct{{\bm{\mathbf{r}}}}{w^{(\mathrm{1})}_{i{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}}.$
(16)

We also define real-space matrix elements

$g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\equiv\matrixelement{w^{(\mathrm{0})}_{i{\bm{\mathbf{0}}}}}{\hat{V}_{\bm{\mathbf{q}}}}{w^{(\mathrm{0})}_{j{\bm{\mathbf{R}}}}},$
(17) $\delta
g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\equiv\matrixelement{w^{(\mathrm{0})}_{i{\bm{\mathbf{0}}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{1})}_{{j\mathbf{R}};{\bm{\mathbf{q}}}}}+\matrixelement{w^{(\mathrm{1})}_{i{\bm{\mathbf{0}}};{\bm{\mathbf{-q}}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{0})}_{j\mathbf{R}}},$
(18)

and

$\widetilde{g}_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\equiv
g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}+\delta
g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}.$ (19)

Here, $g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ is the usual WF-based matrix
element of $\hat{V}_{\bm{\mathbf{q}}}$, and $\delta
g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ is a correction originating from the
change of the WFs. This correction term is not needed for calculating the
scattering matrix elements, but it is necessary when calculating the perturbed
wavefunctions. We define the Fourier transformation of
$g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ as

$g^{\mathrm{(W)}}_{ij{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}=\sum_{{\bm{\mathbf{R}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}},$
(20)

and analogously, we define
$\widetilde{g}^{\mathrm{(W)}}_{{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}$.

####  II.2.2 Perturbed tight-binding eigenvector

The Hamiltonian of the perturbed system in the Bloch representation is

$H^{(\mathrm{\rm
W})}_{i^{\prime}{\bm{\mathbf{k^{\prime}}}},i{\bm{\mathbf{k}}}}\equiv\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}},{\bm{\mathbf{R^{\prime}}}}}H_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}{e^{-i{{\bm{\mathbf{k^{\prime}}}}}\cdot{{\bm{\mathbf{R^{\prime}}}}}}}{e^{i{{\bm{\mathbf{k}}}}\cdot{{\bm{\mathbf{R}}}}}}.$
(21)

Since the unperturbed system has lattice translation symmetry, the zeroth-
order part of $H^{(\mathrm{\rm W})}$ is block diagonal:

$H^{(\mathrm{\rm
W0})}_{i^{\prime}{\bm{\mathbf{k^{\prime}}}},i{\bm{\mathbf{k}}}}=\delta_{{\bm{\mathbf{k}}},{\bm{\mathbf{k^{\prime}}}}}H^{(\mathrm{\rm
W0})}_{i^{\prime}i;{\bm{\mathbf{k}}}}=\delta_{{\bm{\mathbf{k}}},{\bm{\mathbf{k^{\prime}}}}}\sum_{\bm{\mathbf{R}}}H^{(\mathrm{0})}_{i^{\prime}{\bm{\mathbf{0}}},i{\bm{\mathbf{R}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}.$
(22)

The unperturbed Hamiltonian $H^{(\mathrm{\rm W0})}_{{\bm{\mathbf{k}}}}$ can be
diagonalized as

$(V^{(0)\dagger}_{{\bm{\mathbf{k}}}}H^{(\mathrm{\rm
W0})}_{\bm{\mathbf{k}}}V^{(\mathrm{0})}_{\bm{\mathbf{k}}})_{mn}=\varepsilon^{(\mathrm{\rm
H0})}_{m\mathbf{k}}\delta_{m,n},$ (23)

where $V^{(\mathrm{0})}_{\bm{\mathbf{k}}}$ is an
${N_{\mathrm{W}}}\times{N_{\mathrm{W}}}$ unitary matrix. In the Born-von
Karman supercell, the unperturbed eigenvector is

$V^{(\mathrm{0})}_{{i\mathbf{R}},{n\mathbf{k}}}=\frac{1}{\sqrt{N_{k}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}.$
(24)

Next, we consider the perturbative correction. Using Eqs. (14, 16), one can
show that the Fourier transformation of
$H^{(\mathrm{1})}_{{i^{\prime}\mathbf{R^{\prime}}},{i\mathbf{R}}}$ is

$\displaystyle H^{(\mathrm{\rm
W1})}_{i^{\prime}{\bm{\mathbf{k^{\prime}}}},i{\bm{\mathbf{k}}}}=$
$\displaystyle\widetilde{g}^{\mathrm{(W)}}_{i^{\prime}i{\bm{\mathbf{k}}};{\bm{\mathbf{k^{\prime}-k}}}}$
(25)

(see Appendix D for a proof). The matrix element of $H^{(\mathrm{\rm W1})}$ in
the unperturbed eigenstate basis is

$H^{(\mathrm{H1})}_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}=\left(V^{{(\mathrm{0})}\dagger}_{{\bm{\mathbf{k^{\prime}}}}}\widetilde{g}^{\mathrm{(W)}}_{{\bm{\mathbf{k}}};{\bm{\mathbf{k^{\prime}-k}}}}V^{(\mathrm{0})}_{{\bm{\mathbf{k}}}}\right)_{n^{\prime}n}\equiv\widetilde{g}^{(\mathrm{H})}_{n^{\prime}n{\bm{\mathbf{k}}};{\bm{\mathbf{k^{\prime}-k}}}}.$
(26)

We also define $g^{(\mathrm{H})}$ analogously.

We now apply first-order perturbation theory to the tight-binding Hamiltonian.
The perturbed eigenvector in the WF basis is

$\displaystyle
V_{{i\mathbf{R}},{n\mathbf{k}}}=\frac{1}{\sqrt{N_{k}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}$
(27)
$\displaystyle+\lambda\frac{1}{\sqrt{N_{k}}}\sum_{\bm{\mathbf{k^{\prime}}}}\sideset{}{{}^{\prime}}{\sum}_{{n^{\prime}=1}}^{{N_{\mathrm{W}}}}e^{i{\bm{\mathbf{k^{\prime}}}}\cdot{\bm{\mathbf{R}}}}V^{(\mathrm{0})}_{in^{\prime};{\bm{\mathbf{k^{\prime}}}}}\frac{\widetilde{g}^{(\mathrm{H})}_{n^{\prime}n{\bm{\mathbf{k}}};{\bm{\mathbf{k^{\prime}-k}}}}}{\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{n^{\prime}\mathbf{k^{\prime}}}}+{\mathcal{O}}(\lambda^{2}).$

####  II.2.3 Calculation of wavefunction perturbation

We can now calculate the wavefunction perturbation as we know the basis
functions (the perturbed WF $\ket{w_{i\mathbf{R}}}$) as well as the
coefficients (the perturbed eigenvector $V_{{i\mathbf{R}},{n\mathbf{k}}}$) for
the perturbed wavefunction. The perturbed wavefunction reads

$\displaystyle\ket{\psi^{(\mathrm{\rm H})}_{n\mathbf{k}}}=$
$\displaystyle\sum_{i,{\bm{\mathbf{R}}}}\ket{w_{i\mathbf{R}}}V_{{i\mathbf{R}},{n\mathbf{k}}}$
$\displaystyle=$ $\displaystyle\ket{\psi^{(\mathrm{\rm
H0})}_{n\mathbf{k}}}+\lambda\sum_{\bm{\mathbf{q}}}\ket{\psi^{(\mathrm{\rm
H1})}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}+{\mathcal{O}}(\lambda^{2}),$ (28)

where

$\ket{\psi^{\mathrm{(H0)}}_{n\mathbf{k}}}=\frac{1}{\sqrt{N_{k}}}\sum_{i,{\bm{\mathbf{R}}}}\ket{w^{\mathrm{(0)}}_{i\mathbf{R}}}V^{(0)}_{in,{\bm{\mathbf{k}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(29)

and

$\displaystyle\ket{\psi^{(\mathrm{\rm
H1})}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}=$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{i,{\bm{\mathbf{R}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}\ket{w^{(\mathrm{1})}_{{i\mathbf{R}};{\bm{\mathbf{q}}}}}V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}$
$\displaystyle+$
$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{{m=1}}^{{N_{\mathrm{W}}}}\ket{\psi^{(\mathrm{H0})}_{m{\bm{\mathbf{k+q}}}}}\frac{\widetilde{g}^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}}{\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{m{\bm{\mathbf{k+q}}}}}.$
(30)

Equation (II.2.3) contains contributions from both the WFP and the
perturbation of the tight-binding eigenvector
$V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}$.

Equation (II.2.3) is the interpolation equation for the first-order
wavefunction perturbation. This equation is the second main theoretical result
of this paper. The contribution of the states outside
${\hat{P}^{(0)}_{\mathrm{W}}}$, the unperturbed Wannier subspace, to the
wavefunction perturbation is taken into account by the WFPs and is the first
term of Eq. (II.2.3). The contribution of the states inside
${\hat{P}^{(0)}_{\mathrm{W}}}$ is the second term, which is derived by the
perturbation theory in the tight-binding model. Since the WFPs are spatially
localized, the wavefunction perturbation converges quickly with the number of
lattice vectors used in the sum over ${\bm{\mathbf{R}}}$.

###  II.3 Interpolation of matrix elements

Finally, we demonstrate how the matrix elements that involve wavefunction
perturbations can be interpolated. We consider interpolation of two different
types of quantities, $S_{mn{\bm{\mathbf{k}}}}$ [Eq. (II.3.1)] and
$K_{mn{\bm{\mathbf{k}}}}$ [Eq. (37)].

Table 1: Key quantities of WFPT. Description | Real space, Wannier basis | Reciprocal space, eigenstate basis
---|---|---
Wavefunction perturbation |  $\ket{w^{(\mathrm{1})}_{i{\bm{\mathbf{R}}}}}$, Eq. (8) |  $\ket{\psi^{(\mathrm{H1})}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}$, Eq. (II.2.3)
| Perturbation to the
---
Hamiltonian matrix elements
| WF contribution: $g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, Eq. (17)
---
WFP contribution: $\delta g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, Eq. (18)
Total: $\widetilde{g}_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, Eq. (19)
$H^{(\mathrm{H1})}_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}$, Eq. (26)
Sternheimer matrix elements |  $s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, Eq. (33) |  $S_{mn{\bm{\mathbf{k}}}}$, Eqs. (II.3.1, II.3.1)
Kubo matrix elements |  $k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, Eq. (38) |  $K_{mn{\bm{\mathbf{k}}}}$, Eqs. (37, 39)

####  II.3.1 Interpolation of $S_{mn\mathbf{k}}$

First, we consider quantities of the form

$\displaystyle S_{mn{\bm{\mathbf{k}}}}=$
$\displaystyle\matrixelement{\psi^{(\mathrm{1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(\mathrm{0})}_{n\mathbf{k}}}$
$\displaystyle=$
$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{{p}}\frac{\matrixelement{\psi^{(\mathrm{0})}_{m\mathbf{k}}}{\hat{V}_{\bm{\mathbf{q}}}^{\dagger}}{\psi^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}}\matrixelement{\psi^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(\mathrm{0})}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}},$
(31)

where $\hat{V}^{\prime}_{\bm{\mathbf{q}}}$ is a monochromatic perturbation
that may or may not equal $\hat{V}_{\bm{\mathbf{q}}}$. If
$\hat{V}_{\bm{\mathbf{q}}}$ is a phonon potential, interpolation of
$S_{{\bm{\mathbf{k}}}}$ enables an efficient calculation of the phonon-induced
renormalization of the electron band structure. For
$\hat{V}_{\bm{\mathbf{0}}}$ being the spin-velocity operator and
$\hat{V}_{{\bm{\mathbf{q}}}\neq{\bm{\mathbf{0}}}}=0$, $S_{\bm{\mathbf{k}}}$ is
relevant to the shift spin current.

For ${m\mathbf{k}},{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$, one can use
the Wannier interpolation of the unperturbed wavefunction and the wavefunction
perturbation [Eqs. (29, II.2.3)] to find

$\displaystyle S_{mn{\bm{\mathbf{k}}}}=$
$\displaystyle\matrixelement{\psi^{(\mathrm{H1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(\mathrm{H0})}_{{n\mathbf{k}}}}$
$\displaystyle=$
$\displaystyle\frac{1}{N_{k}}\sum_{i,i^{\prime},{\bm{\mathbf{R}}},{\bm{\mathbf{R^{\prime}}}}}V^{(0)\dagger}_{mi^{\prime},{\bm{\mathbf{k}}}}\matrixelement{w^{(\mathrm{1})}_{{i^{\prime}\mathbf{R^{\prime}}};{\bm{\mathbf{q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{w^{(\mathrm{0})}_{i\mathbf{R}}}V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}e^{i{\bm{\mathbf{k}}}\cdot({\bm{\mathbf{R-R^{\prime}}}})}+\sideset{}{{}^{\prime}}{\sum}_{{p=1}}^{{N_{\mathrm{W}}}}\frac{\left(\widetilde{g}^{(\mathrm{H})}_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\right)^{*}\matrixelement{\psi^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(\mathrm{H0})}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}}}$
$\displaystyle=$
$\displaystyle\sum_{i,i^{\prime},{\bm{\mathbf{R}}}}V^{(0)\dagger}_{mi^{\prime},{\bm{\mathbf{k}}}}\matrixelement{w^{(\mathrm{1})}_{i^{\prime}{\bm{\mathbf{0}}};{\bm{\mathbf{q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{w^{(\mathrm{0})}_{{i\mathbf{R}}}}V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}+\sideset{}{{}^{\prime}}{\sum}_{{p=1}}^{{N_{\mathrm{W}}}}\frac{\left(\widetilde{g}^{(\mathrm{H})}_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\right)^{*}{g^{\prime}}^{(\mathrm{H})}_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}}{\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}}}$
(32)

In the third equality, we use the discrete translation properties of
$\hat{V}^{\prime}_{\bm{\mathbf{q}}}$ and the WFPs [Eq. (14, 16)].

We define the real-space matrix element

$s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\equiv\matrixelement{w^{(\mathrm{1})}_{i{\bm{\mathbf{0}}};{\bm{\mathbf{q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{w^{(\mathrm{0})}_{{j\mathbf{R}}}}.$
(33)

Since the WFs and the WFPs are localized,
$s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ rapidly decays as a function of
$\absolutevalue{{\bm{\mathbf{R}}}}$. Defining

$s^{(\mathrm{W})}_{ij{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\equiv\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}}}s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(34)

and

$s^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\equiv\left(V^{(0)\dagger}_{{\bm{\mathbf{k}}}}s^{(\mathrm{W})}_{{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}V^{(0)}_{{\bm{\mathbf{k}}}}\right)_{mn},$
(35)

one can rewrite Eq. (II.3.1) as

$\displaystyle
S_{mn{\bm{\mathbf{k}}}}=s^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}+\sideset{}{{}^{\prime}}{\sum}_{{p=1}}^{{N_{\mathrm{W}}}}\frac{\left(\widetilde{g}^{(\mathrm{H})}_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\right)^{*}{g^{\prime}}^{(\mathrm{H})}_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}}{\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}}}.$
(36)

Since all quantities appearing in Eq. (36) can be calculated using a WF- or
WFP-based interpolation, $S_{mn{\bm{\mathbf{k}}}}$ can also be interpolated.

####  II.3.2 Interpolation of $K_{mn\mathbf{k}}$

Next, we show that the matrix element

$\displaystyle K_{mn{\bm{\mathbf{k}}}}=$
$\displaystyle\innerproduct{\psi^{(\mathrm{1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{\psi^{\prime(1)}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}$
(37) $\displaystyle=$
$\displaystyle\sideset{}{{}^{\prime}}{\sum}_{{p}}\frac{\matrixelement{\psi^{(\mathrm{0})}_{m\mathbf{k}}}{\hat{V}_{\bm{\mathbf{q}}}^{\dagger}}{\psi^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}}\matrixelement{\psi^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(\mathrm{0})}_{n\mathbf{k}}}}{(\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}})(\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}})}.$

can similarly be interpolated using WFPT. This matrix element is relevant to
the evaluation of the Kubo formula for noninteracting quasiparticles. When
${\bm{\mathbf{q}}}=0$ and $\hat{V}_{\bm{\mathbf{0}}}$ and
$\hat{V}^{\prime}_{\bm{\mathbf{0}}}$ are the velocity and spin-velocity
operators, respectively, $K_{nn{\bm{\mathbf{k}}}}$ is directly related to the
spin Berry curvature and the spin Hall conductivity.

Defining the real-space matrix element

$k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\equiv\innerproduct{w^{(1)}_{i{\bm{\mathbf{0}}};{\bm{\mathbf{q}}}}}{w^{(1)}_{i{\bm{\mathbf{0}}};{\bm{\mathbf{q}}}}}{w^{\prime(1)}_{{j\mathbf{R}};{\bm{\mathbf{q}}}}},$
(38)

we show in Appendix E that $K_{mn{\bm{\mathbf{k}}}}$ can be interpolated as

$\displaystyle
K_{mn{\bm{\mathbf{k}}}}=k^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}+\sideset{}{{}^{\prime}}{\sum}_{{p=1}}^{{N_{\mathrm{W}}}}\frac{\left(\widetilde{g}^{\rm(H)}_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\right)^{*}\widetilde{g}^{\prime\rm(H)}_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}}{(\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p\mathbf{k}+\mathbf{q}})(\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p\mathbf{k}+\mathbf{q}})}.$
(39)

Here, we define

$k^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\equiv\left(V^{(0)\dagger}_{{\bm{\mathbf{k}}}}k^{(\mathrm{W})}_{{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}V^{(0)}_{{\bm{\mathbf{k}}}}\right)_{mn}$
(40)

and

$k^{(\mathrm{W})}_{ij{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}\equiv\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}}}k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}\,.$
(41)

In the actual evaluation of response quantities, positive infinitesimal $\eta$
that ensures the causality of the response is included in the denominator of
the low-energy-band contribution, as in Eq. (III.1) and Eq. (61). Numerically,
the value of $\eta$ is decreased until the calculated quantity converges.

In Appendix F, we detail how to calculate the real-space matrix elements
$g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$,
$\widetilde{g}_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$,
$s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, and
$k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ from the Bloch eigenstates at the
coarse $k$ grid. Table 1 summarizes the key quantities of WFPT.

###  II.4 Proof-of-principles example: Berry curvature

Figure 3:  Band-resolved Berry curvature of monolayer WTe2. VBM stands for the
highest-energy valence band, or valence band maximum. The plane-wave
calculation is done by solving the Sternheimer equation, Eq. (II.4), to
directly evaluate Eq. (42).

Now, let us present a proof-of-principles example of WFPT interpolation by
interpolating the band-resolved Berry curvature. Since the Berry curvature
strongly oscillates in momentum space due to the contribution from nearby
bands, we exclude the intra-valence-band contribution as follows:

$\widetilde{\Omega}_{m\mathbf{k}}^{ab}=-2\imaginary\matrixelement{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{a}}}{{\hat{Q}}_{{\bm{\mathbf{k}}},{\rm
valence}}}{\frac{\partial u_{m\mathbf{k}}}{\partial k^{b}}}\,.$ (42)

Here,
$\ket{u_{m\mathbf{k}}}=e^{-i{\bm{\mathbf{k}}}\cdot\hat{{\bm{\mathbf{r}}}}}\ket{\psi_{m\mathbf{k}}}$
and ${\hat{Q}}_{{\bm{\mathbf{k}}},{\rm
valence}}=\hat{\mathbb{1}}-{\hat{P}}_{{\bm{\mathbf{k}}},{\rm valence}}$, with
${\hat{P}}_{{\bm{\mathbf{k}}},{\rm valence}}=\sum_{n\in{\rm
valence}}\ket{u_{n\mathbf{k}}}\bra{u_{n\mathbf{k}}}$ the projection operator
to the valence bands at ${\bm{\mathbf{k}}}$. By “valence bands,” we mean only
the occupied bands inside the outer window for Wannierization,
${\mathcal{W}_{\mathrm{D}}}$: The low-lying bands that are outside
${\mathcal{W}_{\mathrm{D}}}$ are not included in
${\hat{P}}_{{\bm{\mathbf{k}}},{\rm valence}}$.

The Berry curvature can be calculated using WFPT interpolation. Since the
velocity operator is the derivative of the Hamiltonian with respect to
${\bm{\mathbf{k}}}$, the derivative of the wavefunction is equivalent to the
wavefunction perturbation due to a velocity operator perturbation. Hence, by
choosing the velocity operators as the perturbations for the WFPs,
$\hat{V}_{\bm{\mathbf{q}}}=\hat{v}^{a}\delta_{{\bm{\mathbf{q}}},{\bm{\mathbf{0}}}}$
and
$\hat{V}^{\prime}_{\bm{\mathbf{q}}}=\hat{v}^{b}\delta_{{\bm{\mathbf{q}}},{\bm{\mathbf{0}}}}$,
we have

$\ket{\psi^{(1)}_{m\mathbf{k}}}=e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{r}}}}\ket{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{a}}}\text{ and
}\ket{\psi^{{}^{\prime}(1)}_{m\mathbf{k}}}=e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{r}}}}\ket{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{b}}}.$ (43)

Using Eq. (37) and Eq. (39), we find

$\displaystyle\widetilde{\Omega}_{m\mathbf{k}}^{ab}=$
$\displaystyle-2\imaginary\left(\innerproduct{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{a}}}{\frac{\partial u_{m\mathbf{k}}}{\partial
k^{b}}}-\sum_{n\in{\rm valence}}\innerproduct{\frac{\partial
u_{m\mathbf{k}}}{\partial
k^{a}}}{u_{n\mathbf{k}}}\innerproduct{u_{n\mathbf{k}}}{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{b}}}\right)$ $\displaystyle=$
$\displaystyle-2\imaginary\left(\innerproduct{\psi^{(1)}_{m\mathbf{k}}}{\psi^{{}^{\prime}(1)}_{m\mathbf{k}}}-\sum_{n\in{\rm
valence}}\innerproduct{\psi^{(1)}_{m\mathbf{k}}}{\psi^{(0)}_{n\mathbf{k}}}\innerproduct{\psi^{(0)}_{n\mathbf{k}}}{\psi^{{}^{\prime}(1)}_{m\mathbf{k}}}\right)$
$\displaystyle=$
$\displaystyle-2\imaginary\left(K_{mm{\bm{\mathbf{k}}}}-\sideset{}{{}^{\prime}}{\sum}_{{{\begin{subarray}{c}n=1\\\
n\in{\rm
valence}\end{subarray}}}}^{{N_{\mathrm{W}}}}\frac{\left(g^{\rm(H)}_{nm{\bm{\mathbf{k}}}}\right)^{*}g^{\prime{\rm(H)}}_{nm{\bm{\mathbf{k}}}}}{(\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}})^{2}}\right)$
$\displaystyle=$
$\displaystyle-2\imaginary\left(k^{\mathrm{(H)}}_{mm{\bm{\mathbf{k}}}}+\sideset{}{{}^{\prime}}{\sum}_{{{\begin{subarray}{c}n=1\\\
n\notin{\rm
valence}\end{subarray}}}}^{{N_{\mathrm{W}}}}\frac{\left(\widetilde{g}^{{\rm(H)}}_{nm{\bm{\mathbf{k}}}}\right)^{*}\widetilde{g}^{\prime{\rm(H)}}_{nm{\bm{\mathbf{k}}}}}{(\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}})^{2}}\right).$
(44)

This equation is the WFPT interpolation formula for
$\widetilde{\Omega}_{m\mathbf{k}}^{ab}$. In the last equality, we used the
fact that in the limit of exact interpolation,
$g^{\rm(H)}_{nm{\bm{\mathbf{k}}}}=\widetilde{g}^{\rm(H)}_{nm{\bm{\mathbf{k}}}}$
holds for ${m\mathbf{k}},{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$, which
can easily be shown to follow from the
$P^{(\mathrm{0})}_{\mathrm{D}}-P^{(\mathrm{0})}_{\mathrm{W}}$ terms in Eq.
(F.1). Because of the interpolation error, $\delta
g^{(\mathrm{H})}_{nm{\bm{\mathbf{k}}}}=0$ holds exactly only if
${\bm{\mathbf{k}}}$ belongs to the coarse $k$-point grid for Wannier
interpolation. But, the deviation at a generic ${\bm{\mathbf{k}}}$ is small
and can be diminished by using a denser coarse $k$-point grid for
Wannierization, as with all other Wannier-interpolated quantities.

We also calculate the Berry curvature directly in the plane-wave basis by
solving the following Sternheimer equation to calculate
$\hat{Q}_{{\bm{\mathbf{k}}},\rm{valence}}\ket{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{a}}}$:

$\displaystyle\left(\hat{H}({\bm{\mathbf{k}}})-\varepsilon_{m\mathbf{k}}\right)\hat{Q}_{{\bm{\mathbf{k}}},\rm{valence}}\ket{\frac{\partial
u_{m\mathbf{k}}}{\partial k^{a}}}$ $\displaystyle=$
$\displaystyle-\left(\sum_{\begin{subarray}{c}n=1\\\ n\notin{\rm
valence}\end{subarray}}^{\infty}\ket{u_{n{\bm{\mathbf{k}}}}}\bra{u_{n{\bm{\mathbf{k}}}}}\right)\frac{\partial\hat{H}({\bm{\mathbf{k}}})}{\partial
k^{a}}\ket{u_{m\mathbf{k}}},$ (45)

where
$\hat{H}({\bm{\mathbf{k}}})=e^{-i{\bm{\mathbf{k}}}\cdot\hat{{\bm{\mathbf{r}}}}}\hat{H}e^{i{\bm{\mathbf{k}}}\cdot\hat{{\bm{\mathbf{r}}}}}$.
The conventional Wannier interpolation method is also exact for the
calculation of Berry curvature Wang  _et al._ (2006). There,
$\widetilde{\Omega}_{m\mathbf{k}}^{ab}$ is calculated by first calculating the
full Berry curvature and then subtracting the intra-valence-band contribution.

Figure 3 shows the Berry curvature of four valence bands of monolayer WTe2
calculated using WFPT interpolation, direct plane-wave calculations, and exact
conventional Wannier interpolation. The WFPT result agrees perfectly with the
plane-wave and Wannier interpolation results, demonstrating the validity of
the WFPT interpolation method.

##  III Results

###  III.1 Application 1: Temperature-dependent electron band structure and
indirect absorption

Figure 4:  (a, b) Band structure of silicon around the (a) VBM and (b)
conduction band minimum (CBM) computed without and with the EPC at $T$=1000 K.
The red shades around the red curves indicate the phonon-induced broadening,
with the width the same as the imaginary part of the self-energy. (c, d)
Indirect absorption spectra of silicon computed (c) without and (d) with the
EPC-induced band renormalization, and the experimental data from Ref.
Macfarlane  _et al._ (1958). The curves correspond to temperatures 20, 77, 90,
112, 170, 195, 249, 291, 363, and 415 K, from right to left. Figure 5:
Spatial decay of (a) Hamiltonian, (b) EPC, and (c) WFP matrix elements. Each
data point corresponds to the maximum absolute value over all WF pairs, phonon
momentum, and atomic displacements for each lattice vector
${\bm{\mathbf{R}}}$. The $y$ axes are on a logarithmic scale.  Table 2: Time
required to calculate the phonon-induced electron self-energy for 379 $k$
points using direct plane-wave calculations and WFPT interpolation.

| Time (cpu$\cdot$h)
---|---
| Plane wave | WFPT
SCF, DFPT (coarse $q$) | 269 | 269
NSCF + Sternheimer (coarse $k$, coarse $q$) | - | 424
Wannierization | - | 171
WFPT interpolation (fine $k$, coarse $q$) | - | 23
NSCF + Sternheimer (fine $k$, coarse $q$) | 1,241∗ | -
Lower Fan (fine $k$, fine $q$) | 4,450,591∗ | 2,846
Total | 4,452,101∗ | 3,733

∗ Estimated from a calculation with smaller numbers of $k$ and $q$ points.

As the first application of WFPT, we calculate the temperature-dependent
electronic structure and indirect optical absorption spectra of silicon. The
electron self-energy in the dynamical Allen-Heine-Cardona formalism Allen and
Heine (1976); Allen and Cardona (1981, 1983) is the sum of the Fan and the
Debye-Waller (DW) terms Giustino (2017):

$\Sigma_{nn^{\prime}{\bm{\mathbf{k}}}}(\omega)=\frac{1}{N_{q}}\sum_{{\bm{\mathbf{q}}},\nu}\left[\Sigma^{\rm
Fan}_{nn^{\prime}{\bm{\mathbf{k}}};{\mathbf{q}\nu}}(\omega)+\Sigma^{\rm
DW}_{nn^{\prime}{\bm{\mathbf{k}}};{\mathbf{q}\nu}}\right],$ (46)
$\displaystyle\Sigma^{\rm
Fan}_{nn^{\prime}{\bm{\mathbf{k}}};{\mathbf{q}\nu}}(\omega)$
$\displaystyle=\sum_{m}g_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}^{*}\,g_{mn^{\prime}{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}$
$\displaystyle\times\sum_{\pm}\frac{n_{\mathbf{q}\nu}+[1\pm(2f_{m\mathbf{k}+\mathbf{q}}-1)]/2}{\omega-\varepsilon_{m\mathbf{k}+\mathbf{q}}\pm\omega_{\mathbf{q}\nu}+i\eta},$
(47) $\displaystyle\Sigma^{\rm
DW}_{nn^{\prime}{\bm{\mathbf{k}}};{\mathbf{q}\nu}}=$
$\displaystyle\frac{n_{\mathbf{q}\nu}+\frac{1}{2}}{2\omega_{{\mathbf{q}\nu}}}\sum_{\begin{subarray}{c}m,\kappa,a\\\
\kappa^{\prime},a^{\prime}\end{subarray}}\mathcal{D}_{nn^{\prime}}^{\kappa
a\kappa^{\prime}a^{\prime}}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})U_{\kappa
a,\nu}^{*}({\bm{\mathbf{q}}})U_{\kappa^{\prime}a^{\prime},\nu}({\bm{\mathbf{q}}}).$
(48)

Here, $U_{\kappa a,\nu}({\bm{\mathbf{q}}})$ is the displacement pattern for
the phonon eigenmode with branch index $\nu$ and wavevector
${\bm{\mathbf{q}}}$ in units of inverse square root of mass,
$n_{\mathbf{q}\nu}$ the Bose-Einstein occupation factor,
$f_{m\mathbf{k}+\mathbf{q}}$ the Fermi-Dirac occupation factor, and
$\mathcal{D}_{nn^{\prime}}^{\kappa
a\kappa^{\prime}a^{\prime}}({\bm{\mathbf{k}}},{\bm{\mathbf{q}}})$ the DW
matrix element [see Eq. (5) of Ref. Lihm and Park (2020) for the definition].
In practice, the DW matrix element is calculated within the rigid-ion
approximation Allen and Heine (1976); Poncé  _et al._ (2014).

Under the on-the-mass-shell approximation Nery  _et al._ (2018); Miglio  _et
al._ (2020), the diagonal part of the self-energy of the state ${n\mathbf{k}}$
is
$\Sigma_{n\mathbf{k}}=\Sigma_{nn{\bm{\mathbf{k}}}}(\varepsilon_{n\mathbf{k}})$.
In this work, we only consider the diagonal part of the self-energy because
the off-diagonal part has little effect on silicon, which has a band gap much
larger than the band energy renormalization Lihm and Park (2020).

In practical calculations, the Fan term is separated into the lower Fan and
the upper Fan terms, which are the contributions from the low-energy and high-
energy bands, respectively. We write the set of low-energy states as
${\mathcal{W}}_{\mathrm{Lower}}$. The upper bound of the energy window for
${\mathcal{W}}_{\mathrm{Lower}}$ is chosen to be a few eV above
$\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}$. We set
${\mathcal{W}}_{\mathrm{Lower}}$ to be inside the perturbation window, which
is necessary because the lower Fan term [Eq. (III.1)] requires an explicit
Wannier interpolation of the EPC matrix elements for bands inside
${\mathcal{W}}_{\mathrm{Lower}}$.

In the upper Fan term, one neglects the phonon frequency
$\omega_{{\bm{\mathbf{q}}}\nu}$ and $\eta$ in the denominator, which is a good
approximation because
$\absolutevalue{\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}+\mathbf{q}}}$
is always much larger than $\omega_{{\bm{\mathbf{q}}}\nu}$ and $\eta$ as one
is only interested in the self-energy of low-energy states whose energy is
much lower than the upper energy bound for ${\mathcal{W}}_{\mathrm{Lower}}$.
The lower and upper Fan terms read

$\displaystyle\Sigma^{\rm lower\
Fan}_{n{\bm{\mathbf{k}}};{\mathbf{q}\nu}}(\eta)$
$\displaystyle=\sum_{\begin{subarray}{c}m\\\
m{\bm{\mathbf{k+q}}}\in{\mathcal{W}}_{\mathrm{Lower}}\end{subarray}}\absolutevalue{g_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}}^{2}$
$\displaystyle\times\sum_{\pm}\frac{n_{\mathbf{q}\nu}+[1\pm(2f_{m\mathbf{k}+\mathbf{q}}-1)]/2}{\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}+\mathbf{q}}\pm\omega_{\mathbf{q}\nu}+i\eta}\,,$
(49)

and

$\displaystyle\Sigma^{\rm upper\
Fan}_{n{\bm{\mathbf{k}}};{\mathbf{q}\nu}}=(2n_{\mathbf{q}\nu}+1)\sum_{\begin{subarray}{c}m\\\
m{\bm{\mathbf{k+q}}}\notin{\mathcal{W}}_{\mathrm{Lower}}\end{subarray}}\frac{\absolutevalue{g_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}}^{2}}{\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}+\mathbf{q}}}\,,$
(50)

respectively.

The DW self-energy under the rigid-ion approximation Allen and Heine (1976);
Poncé  _et al._ (2014) can be easily Wannier interpolated without any infinite
sum by exploiting the recently developed momentum-operator-based
representation of the DW matrix elements Lihm and Park (2020). The lower Fan
term [Eq. (III.1)] can also be calculated using the conventional Wannier
interpolation method.

Calculation of the upper Fan term [Eq. (50)] is where the WF interpolation
fails and WFPT becomes necessary. The reason is that the upper Fan term
involves a sum over an infinite number of bands. Typically, a few hundred
bands are needed to reach convergence Giustino  _et al._ (2010a); Gonze  _et
al._ (2011). Wannier interpolation of such a large number of bands is
impractical, if not impossible.

Using WFPT, the upper Fan term can be easily interpolated without explicitly
performing the sum over bands. By choosing $\hat{V}_{\bm{\mathbf{q}}}$ and
$\hat{V}^{\prime}_{\bm{\mathbf{q}}}$ as the perturbation potential for the
same phonon mode $\nu$ in the definition of $S_{mn{\bm{\mathbf{k}}}}$ [Eq.
(II.3.1)] and using the WFPT interpolation formula [Eq. (36)], we find

$\displaystyle\sum_{\begin{subarray}{c}m\\\
{m\mathbf{k}+\mathbf{q}}\notin{\mathcal{W}}_{\mathrm{Lower}}\end{subarray}}\frac{\absolutevalue{g_{mn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}}^{2}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{m\mathbf{k}+\mathbf{q}}}$
$\displaystyle=$ $\displaystyle
S_{nn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}-\sum_{\begin{subarray}{c}m=1\\\
{m\mathbf{k}+\mathbf{q}}\in{\mathcal{W}}_{\mathrm{Lower}}\end{subarray}}^{{N_{\mathrm{W}}}}\frac{\left(g^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}\right)^{*}g^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}}{\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}+\mathbf{q}}}$
$\displaystyle=$ $\displaystyle
s^{(\mathrm{H})}_{nn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}+\sum_{\begin{subarray}{c}m=1\\\
{m\mathbf{k}+\mathbf{q}}\notin{\mathcal{W}}_{\mathrm{Lower}}\end{subarray}}^{{N_{\mathrm{W}}}}\frac{\left(\widetilde{g}^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}\right)^{*}g^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};\nu{\bm{\mathbf{q}}}}}{\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}+\mathbf{q}}}.$
(51)

Equation (III.1) is the WFPT interpolation formula for the upper Fan self-
energy. In the second equality, we use the fact that in the limit of exact
interpolation $\delta g^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}}}=0$ and
$g^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}}}=\widetilde{g}^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}}}$
if ${m\mathbf{k}},{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$, as explained in
Sec. II.4.

We calculate the upper Fan and DW terms on a coarse $6\times 6\times 6$
$q$-point grid, as it converges much faster than the lower Fan term. We
develop and use an improved double-grid method to speed up the convergence
with respect to the $q$ grid used to calculate the upper Fan and the DW self-
energies. The improved double-grid method utilizes the fact that the low-
energy acoustic phonons should give a small contribution to the total self-
energy due to the translational invariance. With this improved double-grid
method, one can achieve convergence at a considerably coarser $q$ grid than
the conventional double-grid method Brown-Altvater  _et al._ (2020); Lihm and
Park (2020) (see Appendix G for details).

Figures 4(a) and 4(b) show the renormalized band structure of silicon at
$T$=1000 K computed using WFPT. For comparison, we also show the renormalized
band energies calculated directly in the plane-wave basis. WFPT gives an
accurate result at a much smaller computational cost.

Since WFPT enables a cheap calculation of the eigenvalue renormalization, we
can now calculate quantities that require renormalized eigenvalues at a large
number of $k$ points. As a prototypical example, we study the temperature
dependence of the indirect absorption spectrum. Without the phonon-induced
band-structure renormalization, the perturbative expression for the phonon-
assisted absorption coefficient is Kioupakis  _et al._ (2010); Noffsinger  _et
al._ (2012)

$\displaystyle\alpha(\Omega)$ $\displaystyle=\frac{\pi
q^{2}}{c\epsilon_{0}V\Omega
n(\Omega)}\sum_{\nu,m,n,{\bm{\mathbf{k}}},{\bm{\mathbf{q}}},\pm}\left(n_{\mathbf{q}\nu}+\frac{1}{2}\pm\frac{1}{2}\right)$
$\displaystyle\times(f_{n\mathbf{k}}-f_{m\mathbf{k}+\mathbf{q}})\absolutevalue{{\bm{\mathbf{\lambda}}}\cdot({\bm{\mathbf{S}}}_{1}+{\bm{\mathbf{S}}}_{2})}^{2}$
$\displaystyle\times\delta(\varepsilon_{m\mathbf{k}+\mathbf{q}}-\varepsilon_{n\mathbf{k}}-\hbar\Omega\pm\hbar\omega_{\mathbf{q}\nu}),$
(52)

where

${\bm{\mathbf{S}}}_{1}=\sum_{p=1}^{\infty}\frac{{\bm{\mathbf{v}}}_{np{\bm{\mathbf{k}}}}g_{pm{\bm{\mathbf{k}}};-{\mathbf{q}\nu}}}{\varepsilon_{p\mathbf{k}}-\varepsilon_{m\mathbf{k}+\mathbf{q}}\mp\hbar\omega_{\mathbf{q}\nu}+i\Gamma_{p\mathbf{k}}}$
(53)

and

${\bm{\mathbf{S}}}_{2}=\sum_{p=1}^{\infty}\frac{g_{np{\bm{\mathbf{k}}};-{\mathbf{q}\nu}}{\bm{\mathbf{v}}}_{pm{\bm{\mathbf{k+q}}}}}{\varepsilon_{p{\bm{\mathbf{k+q}}}}-\varepsilon_{n\mathbf{k}}\pm\hbar\omega_{\mathbf{q}\nu}+i\Gamma_{p{\bm{\mathbf{k+q}}}}}.$
(54)

Here, $\epsilon_{0}$ is the vacuum permittivity, $\Omega$ the frequency of the
incident photon, $c$ the speed of light, ${\bm{\mathbf{\lambda}}}$ the light
polarization, and $n(\Omega)$ the frequency-dependent refractive index of the
material. Note that we use the energy conservation imposed by the delta
function to rewrite the denominator of Eq. (53) in a form different from Eq.
(2) of Ref. Kioupakis  _et al._ (2010).

When the phonon-induced band-structure renormalization is included, the
absorption coefficient becomes Patrick and Giustino (2014)

$\displaystyle\alpha(\Omega)=2\frac{4\pi^{2}q^{2}}{\Omega
c\,n(\Omega)}\frac{1}{V}\sum_{\nu,m,n,{\bm{\mathbf{k}}},{\bm{\mathbf{q}}},\pm}\left(n_{\mathbf{q}\nu}+\frac{1}{2}\pm\frac{1}{2}\right)$
(55)
$\displaystyle\times(f_{m\mathbf{k}}-f_{n\mathbf{k}+\mathbf{q}})\absolutevalue{{\bm{\mathbf{\lambda}}}\cdot({\bm{\mathbf{S}}}_{1}+{\bm{\mathbf{S}}}_{2})}^{2}$
$\displaystyle\times\delta(\varepsilon_{m\mathbf{k}+\mathbf{q}}+\real\Sigma_{m\mathbf{k}+\mathbf{q}}-\varepsilon_{n\mathbf{k}}-\real\Sigma_{n\mathbf{k}}-\hbar\Omega\pm\hbar\omega_{\mathbf{q}\nu}).$

The phonon-induced real part of the self-energy appears only in the delta
function while ${\bm{\mathbf{S}}}_{1}$ and ${\bm{\mathbf{S}}}_{2}$ are
unaffected Patrick and Giustino (2014); Giustino (2017).

Often, the absorption spectra are insensitive to the imaginary part of the
self-energy of the intermediate state (see Fig. LABEL:fig:s_conv_indabs(d)).
Following Ref. Noffsinger  _et al._ (2012), we approximate
$\Gamma_{p\mathbf{k}}$ as a constant:

$\Gamma_{p\mathbf{k}}=\Gamma_{p{\bm{\mathbf{k+q}}}}=\Gamma.$ (56)

Also, we broaden the delta function by a Gaussian:

$\delta(\varepsilon)\approx\frac{1}{\gamma\sqrt{\pi}}e^{-\varepsilon^{2}/\gamma^{2}},$
(57)

where $\gamma$ is the broadening parameter.

We correct the band structure and the velocity matrix elements using the
$G_{0}W_{0}$ method while calculating the EPC and the phonon-induced
renormalization at the DFT level. The $G_{0}W_{0}$-corrected band gap is 1.28
eV, which is consistent with a previous calculation Hybertsen and Louie
(1986).

Figures 4(c) and 4(d) show the indirect absorption spectra near the onset
calculated without [Fig. 4(c)] and with [Fig. 4(d)] the phonon-induced band
renormalization. The spectra calculated using the WFPT-interpolated band
renormalization successfully reproduce the experimental spectra for the entire
range of temperatures, from 20 K to 415 K. Previous perturbative calculation
of the indirect absorption spectra used an arbitrary temperature-dependent
shift of the calculated spectra Noffsinger  _et al._ (2012) because the phonon
contribution to the electron self-energy was not considered. Our calculation
is free from such temperature-dependent parameters and thus is an accurate and
predictive calculation of the phonon-assisted indirect absorption.

An alternative method to study EPC is the adiabatic supercell method, where
supercells with distorted lattice structures are used Zacharias  _et al._
(2015); Zacharias and Giustino (2016); Monserrat and Vanderbilt (2016);
Gorelov  _et al._ (2020); Zacharias  _et al._ (2020). The adiabatic supercell
methods are limited in that the phonon emission and absorption processes are
not separated. Consequently, such methods cannot reproduce the characteristic
two-slope behavior near the onset of indirect absorption due to the different
onsets of phonon emission and absorption Noffsinger  _et al._ (2012). Also,
the absorption onset is hard to resolve due to the large computational cost of
supercell calculations and the need for a fine $q$-point sampling. For
infrared-active materials, adiabatic supercell methods also miss the
nonadiabatic effect of the Fröhlich interaction, which is the predominant
mechanism for band renormalization in many materials Miglio  _et al._ (2020).

We note that we used two empirical parameters regarding quasiparticle energies
to compute the spectra in Figs. 4(c) and 4(d). First, we shifted all
theoretical spectra by $-54$ meV to match the experimental ones, which
reflects the typical accuracy of $G_{0}W_{0}$ calculations Rangel  _et al._
(2020). The spectra computed without EPC were additionally shifted by $-93$
meV to match the experimental ones at low temperatures. Second, we scaled the
EPC-induced self-energy by a factor of 1.36, correcting the underestimation
due to the neglect of electron-electron interactions in calculating EPC
Antonius  _et al._ (2014). The scaling factor was obtained by dividing the
experimental slope of the temperature-dependent indirect band-gap
renormalization at the high-temperature region Macfarlane  _et al._ (1958);
Bludau  _et al._ (1974) by the calculated one. We emphasize that these
parameters are temperature independent. The need for these parameters is the
limitation of semilocal DFT and the $G_{0}W_{0}$ method. More sophisticated
methods to treat electron-electron interactions may allow a completely
parameter-free calculation but this is beyond the scope of this study.

The real-space matrix elements calculated using WFs and WFPs are well
localized, as shown in Fig. 5. This fact explains the success of the WF- and
WFP-based interpolation.

We summarized the computational cost for calculating the phonon-induced
electron self-energy in Table 2. We find a speedup of several orders of
magnitude. The major speedup of the interpolation scheme comes from the
calculation of the lower Fan term, where both the $k$ and $q$ points are
finely sampled using electron and phonon WFs. Although this speedup comes from
the conventional WF interpolation, WFPT interpolation enables one to
consistently calculate the self-energy in the WF basis. If WFPT is not used,
the upper and lower Fan parts of the self-energy should have been calculated
in the plane-wave and the WF bases, respectively. Assigning band indices to
those two calculations when the bands are nearly degenerate is a nontrivial
issue. The WFPT interpolation of the upper Fan self-energy for the fine $k$
points (23 cpu$\cdot$h) also gives a considerable speedup of above 50 times
compared to solving the Sternheimer equation in the plane-wave basis (1,241
cpu$\cdot$h). In our calculation, the number of fine $k$ points is quite small
(379), as we sampled only the states inside a small energy window and the
crystal is highly symmetric. The speedup of the WFPT interpolation will be
much greater if a wider energy window is considered, or if the valence and
conduction bands have a larger effective mass. Even for silicon, the speedup
will reach the above-mentioned factor of over 50 as the number of points in
the fine $k$ grid is increased.

We note that the WFP interpolation of EPC-induced electron self-energy is
conceptually completely different from the Wannier interpolation of the $GW$
quasiparticle band structure Hamann and Vanderbilt (2009). In the latter, the
$GW$ quasiparticle energy is first calculated for all bands to be Wannierized
and then Wannier interpolated. However, in the WFP interpolation, the matrix
elements are Wannier interpolated, and then the self-energy is calculated. We
briefly discuss why one should not directly interpolate the phonon-induced
electron self-energy. First, one may try to calculate the self-energy on the
coarse $k$ grid before Wannierization. However, doing so is very inefficient
because one needs to sample $q$ points on a very fine grid to converge the
self-energy. Second, one may try to interpolate only the sum of the upper Fan
and DW terms, as it converges at a relatively coarse $q$ grid. However, to do
so, one must set ${\mathcal{W}}_{\mathrm{Lower}}$ to be outside the outer
window ${\mathcal{W}_{\mathrm{D}}}$ because the upper Fan term should be
calculated for all states inside ${\mathcal{W}_{\mathrm{D}}}$ for
Wannierization. Then, the lower Fan terms cannot be calculated using Wannier
interpolation because the contribution from the states that are excluded from
Wannierization, which are outside ${\mathcal{W}_{\mathrm{F}}}$ but inside
${\mathcal{W}}_{\mathrm{Lower}}$, is not captured. In our method, we get
around all these problems by using WFPT to interpolate the upper Fan self-
energy and by using electron and phonon WFs to efficiently calculate the EPC
on a fine $k$ and $q$ grid.

###  III.2 Application 2: Shift spin current

Figure 6:  (a) Band structure of monolayer WTe2. (b, c) Matrix element of the
(b) generalized spin-velocity derivative and (c) generalized derivative of the
velocity operator. The conduction band $m$ is indicated by the red dashed
curve in panel (a). The valence bands $n$ for panels (b) and (c) are indicated
by the green dotted and blue dash-dotted curves in panel (a), respectively.
Figure 7:  (a) Shift spin and (b) shift charge current of monolayer WTe2
calculated using WFPT, WF with band truncation, and the difference between the
two. In panel (b), the shift current spectrum calculated using an exact WF-
based interpolation Ibañez-Azpiroz  _et al._ (2018) is also shown.  Figure 8:
(a) Isosurface of the spin-up part of a WF. Isovalues are $+0.087$ and
$-0.087$ times the maximum absolute value for the yellow and cyan isosurfaces,
respectively. (b) Isosurface of the spin-down part of a WFP for the spin-
velocity perturbation $j^{x;x}$. Isovalues are $+0.25$ and $-0.25$ times the
maximum absolute value for the yellow and cyan isosurfaces, respectively.
(c)-(e) Spatial decay of the (c) Hamiltonian, (d) spin-velocity, and (e) WFP
matrix elements. Each data point corresponds to the maximum absolute value
over all WF pairs and the directions $a$, $b$, and $s$ for each
${\bm{\mathbf{R}}}$. The $y$ axes are on a logarithmic scale.  Table 3: Time
required to calculate the shift spin current of the WTe2 monolayer using
direct plane-wave calculations and WFPT interpolation.

| Time (cpu$\cdot$h)
---|---
| Plane wave | WFPT
SCF | 1.9 | 1.9
NSCF + Sternheimer (coarse $k$) | - | 47.2
Wannierization | - | 0.2
WFPT interpolation (fine $k$) | - | 3.0
NSCF + Sternheimer (fine $k$) | 52454.9∗ | -
Total | 52456.8∗ | 52.4

∗ Estimated from a calculation with a smaller number of $k$ points.

The second application of WFPT is the shift spin current, which is a second-
order dc spin current response to an ac electric field. The shift spin and
shift charge conductivity can be written as Lihm and Park

$\displaystyle\sigma^{S_{s},a;bc}_{\mathrm{shift}}(\Omega)$
$\displaystyle=\frac{i\pi
q^{3}}{2\hbar^{2}V\Omega^{2}}\sum_{{\bm{\mathbf{k}}},m,n}(d^{s,b;a}_{mn{\bm{\mathbf{k}}}}v^{c}_{nm{\bm{\mathbf{k}}}}-d^{s,c;a}_{nm{\bm{\mathbf{k}}}}v^{b}_{mn{\bm{\mathbf{k}}}})$
$\displaystyle\times(f_{m\mathbf{k}}-f_{n\mathbf{k}})\delta(\Omega+(\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}})/\hbar)$
(58)

where $\Omega$ is the frequency of the incident photon, $s$ and $a$ the spin
polarization axis and flow direction of the current, and $b$ and $c$ the
directions of the external electric field. In an actual calculation, we
broaden the delta function by a Gaussian as in Eq. (57). The matrix element
$d^{s,b;a}_{mn{\bm{\mathbf{k}}}}$ is called the generalized spin-velocity
derivative Lihm and Park  and is defined as

$d^{s;b,a}_{mn{\bm{\mathbf{k}}}}=j^{s;ab}_{mn{\bm{\mathbf{k}}}}+\sideset{}{{}^{\prime}}{\sum}_{{p}}\frac{j^{s;a}_{mp{\bm{\mathbf{k}}}}v^{b}_{pn{\bm{\mathbf{k}}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}+\sideset{}{{}^{\prime}}{\sum}_{{p}}\frac{v^{b}_{mp{\bm{\mathbf{k}}}}j^{s;a}_{pn}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}},$
(59)

where
$\hat{j}^{s;ab}=-\anticommutator{\hat{S}_{s}}{[\hat{r}^{a},[\hat{r}^{b},\hat{H}^{(\mathrm{0})}]]}/2\hbar^{2}$.
We define $\hat{S}_{0}=\hat{\mathbb{1}}$ so that $s=0$ represents the charge
current. In contrast to its charge counterpart, the generalized spin-velocity
derivative is not geometric: It cannot be written as a derivative with respect
to ${\bm{\mathbf{k}}}$. Thus, the shift spin current cannot be calculated
using the conventional WF interpolation method.

By choosing
$\hat{V}_{{\bm{\mathbf{q}}}}=\hat{j}^{s;a}\delta_{{\bm{\mathbf{q}}},{\bm{\mathbf{0}}}}$
and
$\hat{V}^{\prime}_{{\bm{\mathbf{q}}}}=\hat{v}^{b}\delta_{{\bm{\mathbf{q}}},{\bm{\mathbf{0}}}}$
in the definition of $S_{mn{\bm{\mathbf{k}}}}$ [Eq. (II.3.1)], we find

$d^{s;b,a}_{mn{\bm{\mathbf{k}}}}=j^{s;ab}_{mn{\bm{\mathbf{k}}}}+S_{mn{\bm{\mathbf{k}}}}+S_{nm{\bm{\mathbf{k}}}}^{*}.$
(60)

The matrix elements of $\hat{j}^{s;ab}$ can easily be Wannier interpolated via
a direct calculation of the matrix elements in the plane-wave basis for coarse
$k$ points Note (1). WFPT interpolation of $S_{mn{\bm{\mathbf{k}}}}$ enables
the interpolation of $d^{s;b,a}_{mn{\bm{\mathbf{k}}}}$. With the matrix
elements interpolated to a dense grid of $k$ points, one can calculate the
shift spin and shift charge current spectra using Eq. (III.2).

We calculate the shift spin and shift charge current of monolayer WTe2 in the
2H phase Huang  _et al._ (2016). In Fig. 6, we show the interpolated matrix
elements for the shift spin and shift charge currents. For comparison, we also
evaluate Eq. (59) exactly by solving the Sternheimer equation in the plane-
wave basis and approximately by truncating Eq. (59) to include only the
Wannierized bands. We find that this band truncation gives incorrect matrix
elements. In contrast, the WFPT interpolation accurately reproduces the exact
plane-wave basis result. We also find that the matrix element calculated with
an exact WF-based interpolation of the generalized derivative of the interband
dipole matrix element Ibañez-Azpiroz  _et al._ (2018) agrees with the result
of the WFPT interpolation for the charge case. Such an exact WF-based
interpolation is not possible for the spin case.

By solving the Sternheimer equation, one can compute the shift spin current
without any band-truncation error. WFPT enables such a band-truncation-error-
free calculation at a much lower computational cost. To our knowledge, the
shift spin current has not been studied using the Sternheimer equation-based
method, which is devoid of the band-truncation error. Here, we use WFPT to
calculate the shift spin conductivity without the truncation error and study
the effect of band truncation.

Figure 7 shows the calculated shift spin and shift charge conductivity. We
again find that the truncation of Eq. (59) gives considerable error in the
shift current spectra. Also, for the shift charge current, the result of the
exact WF-based interpolation Ibañez-Azpiroz  _et al._ (2018) agrees with the
result of our WFPT interpolation.

The WF and WFP for WTe2 are shown in Figs. 8(a) and 8(b). The WF shown is a
$d_{z^{2}}$-orbital-like state of the W atom, with its spin almost fully
polarized (99%) along the $+z$ direction. The WFP shown is for the spin-
velocity perturbation $\hat{j}^{x;x}$. The result can be qualitatively
understood by considering the spin and orbital parts separately. The spin
operator $\hat{S}_{x}$ rotates the spin polarization by 180∘ around the $x$
axis so that the spin of the WFP is almost fully polarized (97%) along the
$-z$ direction. The velocity operator acts on the orbital part. Wavefunction
perturbation due to the velocity operator is equivalent to the application of
the position operator to the wavefunction Baroni and Resta (1986). Hence, the
orbital part of the WFP will approximately be the orbital part of WF times
$x$, which is indeed the case in Fig. 8(b). Figures 8(c)-8(e) shows that the
real-space matrix elements are spatially localized.

We summarize the computational cost for calculating the shift current in Table
3. The WFPT interpolation method speeds up the shift current calculation by 3
orders of magnitude. We note that the calculation can be further optimized by
using crystal symmetry to compute the Sternheimer matrix elements only in the
irreducible wedge of the Brillouin zone, for both the coarse and the fine $k$
grids. We leave such an optimization for future work.

###  III.3 Application 3: Spin Hall conductivity

Finally, we apply WFPT to the calculation of spin Hall conductivity. The key
quantity for calculating the spin Hall conductivity is the spin Berry
curvature Guo  _et al._ (2008)

$\Omega^{s,ab}_{{m\mathbf{k}}}=2\sum_{n\neq
m}\frac{\imaginary\left[j^{s,a}_{mn{\bm{\mathbf{k}}}}v^{b}_{nm{\bm{\mathbf{k}}}}\right]}{(\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}})^{2}+\eta^{2}}.$
(61)

Here, $\eta$ is a positive infinitesimal, and $a$, $b$, and $s$ are the
directions of the external field, current flow, and the spin polarization,
respectively. The intrinsic dc spin Hall conductivity is given by the Kubo
formula as Guo  _et al._ (2005, 2008)

$\sigma^{s}_{ab}=\frac{q}{\hbar
V}\sum_{{\bm{\mathbf{k}}},m}f_{m\mathbf{k}}\Omega^{s,ab}_{{m\mathbf{k}}}.$
(62)

Here, $V$ is the volume of the unit cell, $q$ the electron charge, and
$f_{m\mathbf{k}}$ the Fermi-Dirac occupation factor. Contrary to the Berry
curvature, the spin Berry curvature is not a geometric quantity. Thus, in
principle, an infinite number of bands should be included to evaluate Eq.
(61). For example, Ref. Guo  _et al._ (2005) included states up to 50 eV above
the valence band maximum to converge the spin Hall conductivity of Si, Ge,
GaAs, and AlAs. Moreover, a very fine $k$-point grid is required to achieve
convergence, which calls for efficient interpolation methods such as Wannier
interpolation. However, conventional Wannier interpolation methods for
computing spin Hall conductivity Qiao  _et al._ (2018); Ryoo  _et al._ (2019)
truncate the sum over band $n$ in Eq. (61) to include only the Wannier-
interpolated bands.

By solving the Sternheimer equation, one can accurately calculate the spin
Berry curvature and the spin Hall conductivity without explicitly including
the high-energy bands. To our knowledge, the spin Hall conductivity has not
been studied using the Sternheimer-equation-based method. However, the
computational cost of the Sternheimer approach is much higher than Wannier
interpolation for the same $k$-point grid. WFPT enables an efficient Wannier-
interpolation calculation of the spin Hall conductivity without any band-
truncation error, taking the advantages of both the Sternheimer method and the
WF method.

The spin Berry curvature can be calculated by WFPT interpolation similarly to
the case of the spin-velocity derivative. We set
$\hat{V}_{\bm{\mathbf{q}}}=\hat{j}^{s,a}\delta_{{\bm{\mathbf{q}}},{\bm{\mathbf{0}}}}$
and
$\hat{V}^{\prime}_{\bm{\mathbf{q}}}=\hat{v}^{b}\delta_{{\bm{\mathbf{q}}},{\bm{\mathbf{0}}}}$
and use Eq. (39) to find

$\displaystyle\frac{1}{2}\Omega^{s,ab}_{{m\mathbf{k}}}\approx$ $\displaystyle
K_{mm{\bm{\mathbf{k}}}}-\beta_{{m\mathbf{k}}}(0)+\beta_{{m\mathbf{k}}}(\eta)=k^{(\mathrm{H})}_{mm{\bm{\mathbf{k}}}}+\beta_{{m\mathbf{k}}}(\eta)$
(63)

where

$\beta_{{m\mathbf{k}}}(\eta)=\sideset{}{{}^{\prime}}{\sum}_{{n=1}}^{{N_{\mathrm{W}}}}\frac{\imaginary\left[\widetilde{j}^{s,a\rm(H)}_{mn{\bm{\mathbf{k}}}}\widetilde{v}^{b\rm(H)}_{nm{\bm{\mathbf{k}}}}\right]}{\left(\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}\right)^{2}+\eta^{2}}.$
(64)

Here, we approximate that for ${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$ and
${n\mathbf{k}}\notin{\mathcal{W}_{\mathrm{F}}}$,
$\absolutevalue{\varepsilon_{m\mathbf{k}}-\varepsilon_{n\mathbf{k}}}\gg\eta$,
so the $\eta$ in the denominator can be ignored. In practice, we avoid the
sums over pairs of fully occupied states ${m\mathbf{k}}$ and ${n\mathbf{k}}$
in Eq. (61) and Eq. (64) since such contributions get canceled in Eq. (62)
Wang  _et al._ (2006).

We calculate the spin Hall conductivity of monolayer WTe2 in the 2H phase. In
Fig. 9, we compare spin Hall conductivity at $T=50$ K computed using the WFPT
interpolation method, conventional Wannier interpolation with band truncation
Ryoo  _et al._ (2019), and the slow but exact Sternheimer method in the plane-
wave basis. The WFPT result agrees perfectly with the exact results, while the
WF interpolation with band truncation shows some discrepancy, albeit small.
Interestingly, the band-truncation error is present even for the hole-doped
case, where all the relevant bands are separated from the boundary of the
inner window by more than 3 eV. Table 4 shows the speedup of WFPT
interpolation compared to the plane-wave calculation for a $200\times
200\times 1$ $k$-point grid. WFPT gives a speedup of over 2 orders of
magnitude, without lowering the accuracy by band truncation.

Figure 9:  Spin Hall conductivity of monolayer WTe2 computed using WFPT, using
WF with band truncation, and by solving the Sternheimer equation in the plane-
wave basis. The dashed vertical lines indicate the valence band maximum and
the conduction band minimum energies. The temperature is set to 50 K. The sum
over $k$ points in Eq. (62) is performed over an $N_{k}\times N_{k}\times 1$
grid, and $\eta$ is the broadening parameter in Eq. (61). The plane-wave
calculation is not performed for $N_{k}=200$ due to the high computational
cost.  Table 4: Time required to calculate the spin Hall conductivity of the
WTe2 monolayer using a $200\times 200\times 1$ $k$-point grid by direct plane-
wave calculations and WFPT interpolation.

| Time (cpu$\cdot$h)
---|---
| Plane wave | WFPT
SCF | 1.9 | 1.9
NSCF + Sternheimer (coarse $k$) | - | 47.2
Wannierization | - | 0.2
WFPT interpolation (fine $k$) | - | 0.1
NSCF + Sternheimer (fine $k$) | 13113.7∗ | -
Total | 13115.6∗ | 49.4

∗ Estimated from a calculation with a smaller number of $k$ points.

##  IV Discussion

If one is only interested in the valence bands of an insulator, one may choose
all the windows, ${\mathcal{W}_{\mathrm{P}}}$, ${\mathcal{W}_{\mathrm{F}}}$,
and ${\mathcal{W}_{\mathrm{D}}}$, to include only the occupied bands. If WFPT
is applied to the electronic dielectric response of such insulators, the WFP
derived in this work simplifies to the result of Ref. Ge and Lu (2015), i.e.,
Eq. (8) without the second term on the right-hand side. We note, however, that
even the Wannierization of only the valence bands may not be possible in the
case of narrow-gap semiconductors or topological insulators where the orbital
characters of the valence bands change rapidly with k. In the case of such
semiconductors, where the low-energy conduction bands need to be interpolated
and thus disentanglement is required, or that of metals, the simplified
results do not apply. Our formulation of WFPT is applicable to both isolated
and entangled bands and provides a clear physical understanding of why the
WFPs are spatially localized. The generalization of WFPT from isolated bands
to entangled bands requires the deliberate use of energy windows and subspaces
to separate the energy ranges. In addition, the changes in the Wannier basis
functions and the resulting changes in the Hamiltonian matrix [Eq. (12)] must
be addressed to correctly interpolate matrix elements involving the sum over
an infinite number of bands as done in Eqs. (II.3.1,39). For WFs constructed
without disentanglement, the Hamiltonian matrix elements do not change because
the WFPs are spanned by eigenstates outside the outer window so that the
Hamiltonian matrix elements between the WFPs and the unperturbed WFs [Eq.
(18)] are zero. Only as these highly nontrivial problems are solved, the
methods for accurately calculating the proper sum over an infinite number of
bands can be developed (as summarized in Table 1), and WFPT enables several
new kinds of calculations that were beyond the computational reach.

WFPT could find wide applications. For EPC, WFPT opens up a way to
systematically study the effect of phonon-induced band renormalization to
quantities that require integration over fine $k$ points. Phonon-assisted
optical properties are one example. We emphasize that the predictive
calculation of the indirect absorption onset (Fig. 4(d)) is currently possible
only with the WFPT interpolation. Another example is electronic transport. In
transport calculations, the EPC-induced band renormalization is ignored or
included only partially as a correction to the effective mass Schlipf  _et
al._ (2018); D’Souza  _et al._ (2020). WFPT interpolation will allow the
complete consideration of the renormalization by calculating the band
renormalization at every $k$ point in the fine grid for transport
calculations. Also, a self-consistent calculation of the phonon-induced
electron self-energy Brown-Altvater  _et al._ (2020) without neglecting the
$k$ dependence of the self-energy may become feasible.

Zero-point renormalization and temperature dependence of band structures are
also interesting by themselves. Calculation of the zero-point renormalization
for a large set of materials began only very recently Miglio  _et al._ (2020).
Also, most calculations are limited to the correction of the valence band
maximum and conduction band minimum, partly due to the large computational
cost of calculations for a large number of $k$ points. WFPT will considerably
reduce the barrier of the calculation of the zero-point and temperature-
dependent renormalization for a large number of materials or for many $k$
points. Combined with the automatic generation of Wannier functions Vitale
_et al._ (2020), an efficient, automatized high-throughput calculation will
also be possible.

Similar developments can also be made for the shift spin current. The
interpolation scheme based on WFPs is necessary to calculate converged shift
spin current spectra at a reasonable computational cost. Considering the non-
geometric property of the spin-transport-related matrix elements, WFPT will be
useful for studying a wide range of spin transport phenomena.

Another promising application of WFPT is the calculation of the nonadiabatic
Born effective charge, a concept recently introduced in Refs. Bistoni  _et
al._ (2019); Binci  _et al._ (2021); Dreyer  _et al._ (2021). Calculation of
the nonadiabatic Born effective charge of metals and doped semiconductors
requires a fine $k$-point sampling of wavefunction perturbations for phonon
perturbations Dreyer  _et al._ (2021). WFPT interpolation allows such a
calculation for all kinds of materials, not just for the special case where
the low-energy conduction bands are isolated from the high-energy bands.

Finally, we discuss two ways the WFPT formalism itself can be expanded. First,
in this work, we did not consider the localized representation of the
perturbations. The concept of WFPs may be generalized to localized
perturbations, for example by using phonon Wannier functions Giustino  _et
al._ (2007); Eiguren and Ambrosch-Draxl (2008). Such a generalization may
enable the study of phonon-induced band renormalization in complex
heterostructures by stitching WF- and WFP-based real-space matrix elements.
Second, considering time-periodic perturbations will allow the localized
representation of time-dependent wavefunction perturbations. This
generalization will be useful for the study of dynamic properties such as the
dielectric function and the higher-harmonic generation. Also, these two
generalizations may enable the application of WFPT to $GW$ calculations,
especially in the calculation of frequency- and momentum-dependent screened
Coulomb interaction using the Sternheimer equation Giustino  _et al._ (2010b)
or the Lanczos method Umari  _et al._ (2010).

##  V Conclusion

In conclusion, we introduced the concept of WFPT, which is the extension of
the nearsightedness principle to wavefunction perturbations. We defined the
WFP and showed that it gives an accurate localized representation of
wavefunction perturbations. Using WFPs, various matrix elements involving the
wavefunction perturbation can be efficiently interpolated. We applied WFPT to
the calculation of the temperature-dependent electron band structures, the
shift spin conductivity, and the spin Hall conductivity. The significant
speedup gained thanks to WFPT enabled the study of phenomena such as the
temperature-dependent onset of the indirect optical absorption and the shift
spin currents, which were beyond the reach of previous methods. As density-
functional perturbation theory has significantly expanded the applications of
density functional theory, WFPT will also widen the scope of Wannier function
methods for the theoretical and computational study of condensed matter.

##  VI Methods

We used the Quantum ESPRESSO package to perform plane-wave pseudopotential DFT
calculations Giannozzi  _et al._ (2009, 2017). We used optimized norm-
conserving Vanderbilt (ONCV) pseudopotentials Hamann (2013) taken from the
PseudoDojo library (v0.4) van Setten  _et al._ (2018). For silicon, we used
the Perdew-Wang functional Perdew and Wang (1992) and a scalar-relativistic
pseudopotential. We used a kinetic energy cutoff of 70 Ry and a
14$\times$14$\times$14 $k$-point grid for DFT and DFPT calculations. DFPT
calculation was done on a $6\times 6\times 6$ $q$-point grid. The experimental
lattice parameter at 300 K ($a$=5.43 Å) was used Reeber and Wang (1996).
Thermal expansion was not taken into account as its effect on the indirect
band gap is negligible Lautenschlager  _et al._ (1985); Noffsinger  _et al._
(2012). For the WTe2 monolayer, we used the Perdew-Burke-Ernzerhof functional
Perdew  _et al._ (1996) and included spin-orbit coupling using a fully
relativistic pseudopotential. We used a kinetic energy cutoff of 70 Ry and a
$12\times 12\times 1$ $k$-point grid. The lattice structure was optimized
until the internal stresses and forces were below 0.01 kbar and $2\times
10^{5}$ Ry/Bohr, respectively. The optimized lattice parameters of a
rectangular in-plane supercell were $a$=3.56 Å and $b$=6.17 Å, consistent with
a previous calculation Huang  _et al._ (2016). The supercell length along the
nonperiodic direction was set to 15 Å.

To construct the MLWFs, we used the Wannier90 package Pizzi  _et al._ (2020).
For silicon, we constructed 10 MLWFs from 18 bands to accurately interpolate
the bands near the conduction band minimum. We used $sp^{3}$ orbitals centered
at the Si atom and $s$ orbitals centered at the two interstitial sites as the
initial guesses. The upper bounds of the perturbation window
(${\mathcal{W}_{\mathrm{P}}}$), the inner window of the unperturbed system
(${\mathcal{W}_{\mathrm{F}}}$), and the outer window
(${\mathcal{W}_{\mathrm{D}}}$) were 5.6 eV, 6.1 eV, and 18.5 eV above the
conduction band minimum, respectively. A $12\times 12\times 12$ $k$-point grid
was used. For WTe2, we constructed 28 MLWFs from 32 bands, using W $s$ and
$d$, and Te $s$ and $p$ orbitals as initial guesses. We excluded the semicore
states in the Wannierization. The upper bound of the perturbation window and
the inner window of the unperturbed system were 2.8 eV and 3.3 eV above the
valence band maximum, respectively. The outer window was not set so that all
32 bands are used for Wannierization. A $12\times 12\times 1$ $k$-point grid
was used for the shift current and spin Hall applications. In the Berry
curvature calculation, we used a $24\times 24\times 1$ grid to reduce the
finite-difference error in the position matrix elements Marzari and Vanderbilt
(1997). For plotting the isosurfaces of the WFs and WFPs for WTe2, we used a
coarser $6\times 6\times 6$ $k$-point grid due to the high computational cost.
For all Wannier interpolation, we used the minimal-distance replica selection
based on the Wigner-Seitz supercell that improves the quality of interpolation
Yates  _et al._ (2007).

The lower Fan term in the phonon-induced electron self-energy of silicon was
calculated using a $125\times 125\times 125$ $q$-point grid. The denominator
was smoothed with $\eta=$10 meV [Eq. (III.1)]. We used the on-the-mass-shell
approximation (Rayleigh-Schrödinger perturbation) which is known to give a
better estimate of the zero-point renormalization than the Dyson-Migdal
approach Nery  _et al._ (2018); Brown-Altvater  _et al._ (2020). The upper Fan
and the DW self-energies were calculated on a coarse $6\times 6\times 6$
$q$-point grid without interpolation over $q$ points. To speed up the
convergence of the sum of the upper Fan and the DW self-energies, we used the
improved double-grid method as detailed in Appendix G. See Sec.
LABEL:sec:s_conv_selfen of the Supplementary Material Note (1) for the
convergence study of the lower Fan term.

The $G_{0}W_{0}$ calculation of the quasiparticle band structure of silicon
was done using BerkeleyGW Deslippe  _et al._ (2012), within the plasmon-pole
approximation Hybertsen and Louie (1986). We converged the direct and indirect
band gaps to within 5 meV. We used a $12\times 12\times 12$ $q$-point grid, a
cutoff of 30 Ry and 402 bands for the dielectric matrix, and 308 bands with
the static remainder correction Deslippe  _et al._ (2013) for the self-energy.
The quasiparticle energies were calculated on a $12\times 12\times 12$
$k$-point grid and then Wannier interpolated Hamann and Vanderbilt (2009)
using the EPW code Poncé  _et al._ (2016). The velocity matrix elements for
the indirect optical absorption were renormalized following Ref. Rohlfing and
Louie (2000).

The indirect absorption spectra of silicon were calculated using the EPW code
Poncé  _et al._ (2016). We used a 16-band WF model to include a sufficient
number of intermediate bands. We checked that the indirect absorption spectra
change only slightly when more bands are included. We used $60\times 60\times
60$ $k$\- and $q$-point grids. We sampled the initial and final states only in
the $\pm$1 eV window around the center of the band gap. Using the lattice
symmetry, we only sampled the $k$ points in the irreducible wedge of the
Brillouin zone. In total, 379 $k$ points were sampled. The imaginary part of
the self-energy of intermediate electronic states was set to 50 meV. Delta
functions for energy conservation were broadened using a Gaussian function of
width 15 meV. We used the experimental frequency-dependent refractive index at
300 K Green and Keevers (1995). See Sec. LABEL:sec:s_conv_indabs of the
Supplementary Material Note (1) for the convergence study.

The shift conductivity of monolayer WTe2 was calculated using a $400\times
400\times 1$ $k$-point grid. Delta functions for energy conservation were
broadened using a Gaussian function of width 10 meV. Following the scheme of
Refs. Nastos and Sipe (2006); Ibañez-Azpiroz  _et al._ (2018), we regularized
the denominators including an intermediate-state energy with $\eta=$10 meV,
which is in the range where the spectrum calculated using the given $k$-point
grid remains stable Nastos and Sipe (2006); Ibañez-Azpiroz  _et al._ (2018).
See Sec. LABEL:sec:s_conv_shift of the Supplementary Material Note (1) for the
convergence study.

The spin Hall conductivity of monolayer WTe2 was calculated using a $200\times
200\times 1$ $k$-point grid. Note that $\eta=$10 meV was used in the
calculation of the spin Berry curvature [Eq. (61)]. We used the Smoother
module of the WannierBerri package Tsirkin (2021) to compute the spin Hall
conductivity at a finite temperature. The response was scaled using an
effective single-layer thickness of 7.1 $\mathrm{\AA}$, which is taken from
the bulk lattice parameter Lee  _et al._ (2015). See Sec. LABEL:sec:s_conv_shc
of the Supplementary Material Note (1) for the convergence study.

The WFPs and the corresponding matrix elements in the coarse $k$ grid were
calculated using an in-house modified version of the ph.x program of Quantum
ESPRESSO Giannozzi  _et al._ (2017). We used in-house modified versions of the
EPW Poncé  _et al._ (2016) and Wannier90 Pizzi  _et al._ (2020) code for the
EPC, and the shift spin and spin Hall current applications, respectively. In
the EPC calculation, we used the crystal symmetry to efficiently unfold the
electron-phonon, DW, and Sternheimer matrix elements calculated for the $q$
points in the irreducible wedge to the full coarse $q$-point grid. We also
implemented the optimization Tsirkin (2021) of the minimal-distance replica
selection method that significantly reduces the computational cost in the in-
house versions of EPW and Wannier90. Isosurface plots were created using VESTA
Momma and Izumi (2011).

###### Acknowledgements.

This work was supported by the Creative-Pioneering Research Program through
Seoul National University, the NRF of Korea No-2020R1A2C1014760, and the
Institute for Basic Science (No. IBSR009-D1). Computational resources have
been provided by KISTI (KSC-2020-INO-0078).

##  Appendix A Derivation of the Wannier function perturbation

In this appendix, we derive Eq. (8), the formula for the WFP. First, we find
the disentanglement subspace following the initialization scheme of Ref. Souza
_et al._ (2001). Concretely, we project ${\hat{P}^{(0)}_{\mathrm{W}}}$, the
projection operator for the unperturbed Wannier subspace, onto the subspace
${\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}$ and diagonalize it. We
select $N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}-N_{\rm P}$ states with the largest
eigenvalues Souza  _et al._ (2001) and use them as the basis states for the
perturbed WFs, together with $N_{\mathrm{P}}$ states inside
${\mathcal{W}_{\mathrm{P}}}$. From Eqs. (1, 3, 4), we find that for states
${n\mathbf{k}},{n^{\prime}\mathbf{k^{\prime}}}\in{\mathcal{W}_{\mathrm{D}}}$,
the matrix element of ${\hat{P}^{(0)}_{\mathrm{W}}}$ is

$\displaystyle P_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}=$
$\displaystyle\matrixelement{\psi_{n^{\prime}\mathbf{k^{\prime}}}}{{\hat{P}^{(0)}_{\mathrm{W}}}}{\psi_{n\mathbf{k}}}=\sum_{{\bm{\mathbf{R}}}}\sum_{i=1}^{{N_{\mathrm{W}}}}\innerproduct{\psi_{n^{\prime}\mathbf{k^{\prime}}}}{w_{i\mathbf{R}}^{(\mathrm{0})}}\innerproduct{w_{i\mathbf{R}}^{(\mathrm{0})}}{\psi_{n\mathbf{k}}}$
$\displaystyle=$
$\displaystyle\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}}}\sum_{i=1}^{{N_{\mathrm{W}}}}U_{n^{\prime}i;{\bm{\mathbf{k^{\prime}}}}}^{(\mathrm{0})}U_{in;{\bm{\mathbf{k}}}}^{(0)\dagger}e^{i({\bm{\mathbf{k-k^{\prime}}}})\cdot{\bm{\mathbf{R}}}}$
$\displaystyle+$
$\displaystyle\lambda\frac{1}{N_{k}}\sum_{\bm{\mathbf{R}}}\sideset{}{{}^{\prime}}{\sum}_{{m,{\bm{\mathbf{k^{\prime\prime}}}}}}\sum_{i=1}^{{N_{\mathrm{W}}}}\left[U_{n^{\prime}i;{\bm{\mathbf{k^{\prime}}}}}^{(\mathrm{0})}U_{im;{\bm{\mathbf{k^{\prime\prime}}}}}^{(0)\dagger}e^{i({\bm{\mathbf{k^{\prime\prime}-k^{\prime}}}})\cdot{\bm{\mathbf{R}}}}\frac{\matrixelement{\psi_{m{\bm{\mathbf{k^{\prime\prime}}}}}^{(\mathrm{0})}}{\hat{V}}{\psi^{(\mathrm{0})}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{m{\bm{\mathbf{k^{\prime\prime}}}}}}+({n\mathbf{k}}\leftrightarrow{n^{\prime}\mathbf{k^{\prime}}})^{*}\right]+{\mathcal{O}}(\lambda^{2})$
$\displaystyle=$ $\displaystyle P_{{\rm
W},n^{\prime}n;{\bm{\mathbf{k}}}}^{(\mathrm{0})}\delta_{{\bm{\mathbf{k}}},{\bm{\mathbf{k^{\prime}}}}}+\lambda\sideset{}{{}^{\prime}}{\sum}_{{m}}\left(P_{{\rm
W},n^{\prime}m;{\bm{\mathbf{k^{\prime}}}}}^{(\mathrm{0})}\frac{\matrixelement{\psi_{m\mathbf{k^{\prime}}}^{(\mathrm{0})}}{\hat{V}}{\psi_{n\mathbf{k}}^{(\mathrm{0})}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{m\mathbf{k^{\prime}}}}-\frac{\matrixelement{\psi_{n^{\prime}\mathbf{k^{\prime}}}^{(\mathrm{0})}}{\hat{V}}{\psi_{m\mathbf{k}}^{(\mathrm{0})}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}}P_{{\rm
W},mn;{\bm{\mathbf{k}}}}^{(\mathrm{0})}\right)+{\mathcal{O}}(\lambda^{2}).$
(65)

Here, we defined

$P_{{\rm
W},n^{\prime}n;{\bm{\mathbf{k}}}}^{(\mathrm{0})}=\sum_{i=1}^{{N_{\mathrm{W}}}}U_{n^{\prime}i;{\bm{\mathbf{k}}}}^{(\mathrm{0})}U_{in;{\bm{\mathbf{k}}}}^{(0)\dagger}.$
(66)

We can obtain the eigenstates of
$P_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}$ using first-order
perturbation theory. The zeroth-order matrix is $P_{{\rm
W},n^{\prime}n;{\bm{\mathbf{k}}}}^{(\mathrm{0})}\delta_{{\bm{\mathbf{k}}},{\bm{\mathbf{k^{\prime}}}}}$.
We write its eigenvectors as
$\widetilde{U}^{(\mathrm{0})}_{np;{\bm{\mathbf{k}}}}$, where the eigenvalues
for $p=1,\cdots,{N_{\mathrm{W}}}$ are 1 (in subspace 1 and 2) and those for
$p={N_{\mathrm{W}}}+1,\cdots,N_{\mathrm{D},{\bm{\mathbf{k}}}}$ are 0 (in
subspace 3). Since ${\mathcal{W}_{\mathrm{F}}}$ is the inner window for
unperturbed WFs, one can choose the eigenvectors in
${\mathcal{W}_{\mathrm{F}}}$ to satisfy

$\widetilde{U}^{(\mathrm{0})}_{np;{\bm{\mathbf{k}}}}=\delta_{n,p}\quad{\rm
for}\quad p=1,\cdots,{N_{\mathrm{F},{\bm{\mathbf{k}}}}}\,.$ (67)

Since the basis states are the perturbed wavefunctions, the zeroth-order
eigenstates of $P_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}$ are

$\ket{\phi^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}=\sum_{n=1}^{N_{\mathrm{D},{\bm{\mathbf{k}}}}}\ket{\psi_{{n\mathbf{k}}}}\widetilde{U}^{(\mathrm{0})}_{np;{\bm{\mathbf{k}}}}.$
(68)

The $N_{\rm P}$ energy eigenstates inside ${\mathcal{W}_{\mathrm{P}}}$, the
inner window for the perturbed system, are chosen without further
disentanglement:

$\left\\{\ket{\phi^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}:p=1,\cdots,{N_{\mathrm{P},{\bm{\mathbf{k}}}}}\right\\}\,.$
(69)

Next, we choose the remaining $N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}-N_{\rm P}$
states for the disentangled subspace as follows. We apply first-order
perturbation theory to $\ket{\phi^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}$ for
$p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{W}}}$, which are
the zeroth-order eigenstates of
$P_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}$ inside
${\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}$ with eigenvalue 1.
This choice maximizes the perturbed eigenvalues because the unperturbed
eigenvalues of the states that are not selected are 0. The mixing between the
selected states is not important because our goal here is to determine the
subspace spanned by them. The WFs that are obtained by projection depend only
on the selected subspace, not on its basis Souza  _et al._ (2001). Hence, we
use nondegenerate perturbation theory and only consider the mixing of the
eigenstates of $P_{{\rm W};{\bm{\mathbf{k}}}}^{(\mathrm{0})}$ with eigenvalue
1 with those with eigenvalue 0. Writing the first-order matrix elements in Eq.
(A) as $P^{(\mathrm{1})}_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}$, we
find that $\ket{\phi_{p{\bm{\mathbf{k}}}}^{(\mathrm{0})}}$ is perturbed as

$\displaystyle\ket{\phi_{p{\bm{\mathbf{k}}}}}-\ket{\phi^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}$
$\displaystyle=\lambda\sum_{{\bm{\mathbf{k^{\prime}}}}}\sum_{q={N_{\mathrm{W}}}+1}^{N_{\mathrm{D},{\bm{\mathbf{k^{\prime}}}}}}\sum_{n,n^{\prime}}\ket{\phi^{(\mathrm{0})}_{q{\bm{\mathbf{k^{\prime}}}}}}\innerproduct{\phi^{(\mathrm{0})}_{q{\bm{\mathbf{k^{\prime}}}}}}{\psi_{n^{\prime}\mathbf{k^{\prime}}}}P^{(\mathrm{1})}_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}\innerproduct{\psi_{n\mathbf{k}}}{\phi^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}+{\mathcal{O}}(\lambda^{2})$
$\displaystyle=\lambda\sideset{}{{}^{\prime}}{\sum}_{{n,m,n^{\prime},{\bm{\mathbf{k^{\prime}}}}}}\sum_{q={N_{\mathrm{W}}}+1}^{N_{\mathrm{D},{\bm{\mathbf{k^{\prime}}}}}}\ket{\phi^{(\mathrm{0})}_{q{\bm{\mathbf{k^{\prime}}}}}}\widetilde{U}^{(0)\dagger}_{qn^{\prime};{\bm{\mathbf{k^{\prime}}}}}\left(P_{{\rm
W},n^{\prime}m;{\bm{\mathbf{k^{\prime}}}}}^{(\mathrm{0})}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k^{\prime}}}}{\hat{V}}{\psi^{(0)}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{m\mathbf{k^{\prime}}}}-\frac{\matrixelement{\psi^{(0)}_{n^{\prime}\mathbf{k^{\prime}}}}{\hat{V}}{\psi^{(0)}_{m\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{n^{\prime}\mathbf{k^{\prime}}}}P_{{\rm
W},mn;{\bm{\mathbf{k}}}}^{(\mathrm{0})}\right)\widetilde{U}^{(\mathrm{0})}_{np;{\bm{\mathbf{k}}}}+{\mathcal{O}}(\lambda^{2})$
$\displaystyle=-\lambda\sum_{m}({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}})\ket{\psi_{m\mathbf{k}}^{(\mathrm{1})}}\widetilde{U}^{(\mathrm{0})}_{mp;{\bm{\mathbf{k}}}}+{\mathcal{O}}(\lambda^{2}).$
(70)

In the third equality, we used the fact that for
$p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,...,{N_{\mathrm{W}}}$ and
$q={N_{\mathrm{W}}}+1,...,{N_{\mathrm{D},{\bm{\mathbf{k}}}}}$,
$\widetilde{U}^{(\mathrm{0})}_{np;{\bm{\mathbf{k}}}}$ and
$\widetilde{U}^{(0)\dagger}_{qn^{\prime};{\bm{\mathbf{k^{\prime}}}}}$ are the
right and left eigenvectors of
$P_{\mathrm{W};{\bm{\mathbf{k}}}}^{(\mathrm{0})}$ with eigenvalues 1 and 0,
respectively.

The unperturbed state $\ket{\phi_{p{\bm{\mathbf{k}}}}^{(\mathrm{0})}}$ defined
in Eq. (68) also contains an ${\mathcal{O}}(\lambda^{1})$ contribution because
the basis states are $\ket{\psi_{n\mathbf{k}}}$’s, the perturbed eigenstates.
Using Eq. (4), the perturbative expansion of the disentangled state
$\ket{\phi_{p{\bm{\mathbf{k}}}}}$ becomes

$\ket{\phi_{p{\bm{\mathbf{k}}}}}=\sum_{n}\ket{\psi^{(0)}_{n\mathbf{k}}}\widetilde{U}^{(\mathrm{0})}_{np;{\bm{\mathbf{k}}}}+\lambda\sum_{m}\left(\hat{\mathbb{1}}-({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}})\right)\ket{\psi_{m\mathbf{k}}^{(\mathrm{1})}}\widetilde{U}^{(\mathrm{0})}_{mp;{\bm{\mathbf{k}}}}+{\mathcal{O}}(\lambda^{2})$
(71)

for $p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,...,{N_{\mathrm{W}}}$.

Combining results in Eq. (67) and Eq. (71), the final disentangled subspace
for the perturbed system is spanned by

$\bigcup_{{\bm{\mathbf{k}}}}\left(\left\\{\ket{\phi^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}:p=1,\cdots,{N_{\mathrm{P},{\bm{\mathbf{k}}}}}\right\\}\cup\left\\{\ket{\phi_{p{\bm{\mathbf{k}}}}}:p=N_{{\rm
P},{\bm{\mathbf{k}}}}+1,\cdots,{N_{\mathrm{W}}}\right\\}\right).$ (72)

The dimension of the disentangled subspace is
$N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}$, as desired.

Next, to make computation simpler, we linearly combine the basis states and
write them in a simpler form. For $\lambda=0$, the disentangled subspace [Eq.
(72)] is identical to the unperturbed Wannier subspace. Hence, we can linearly
combine the basis states to make their first-order terms orthogonal to the
unperturbed Wannier subspace. Mathematically, this orthogonalization is
equivalent to applying
$\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)$ to the first-
order terms:

$\displaystyle\ket{\widetilde{\phi}_{p\mathbf{k}}}=\begin{dcases}\sum_{m}\left[\ket{\psi_{m\mathbf{k}}^{(\mathrm{0})}}+\lambda\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi_{m\mathbf{k}}^{(\mathrm{1})}}\right]\widetilde{U}^{(\mathrm{0})}_{mp;{\bm{\mathbf{k}}}}+{\mathcal{O}}(\lambda^{2})&\text{for
}p=1,\cdots,{N_{\mathrm{P},{\bm{\mathbf{k}}}}},\\\
\sum_{m}\left[\ket{\psi_{m\mathbf{k}}^{(\mathrm{0})}}+\lambda\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{D}}}\right)\ket{\psi_{m\mathbf{k}}^{(\mathrm{1})}}\right]\widetilde{U}^{(\mathrm{0})}_{mp;{\bm{\mathbf{k}}}}+{\mathcal{O}}(\lambda^{2})&\text{for
}p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{W}}}.\end{dcases}$
(73)

This application of
$\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)$ guarantees that
the basis states for the disentangled subspace,
$\ket{\widetilde{\phi}_{p\mathbf{k}}}$’s, are orthonormal to each other with
errors of order ${\mathcal{O}}(\lambda^{2})$. Although the final form of the
WFs of the perturbed system is not affected by this orthonormalization, it
makes the following analysis simpler.

Using Eq. (67) and the fact that $1\leq
m\leq{N_{\mathrm{P},{\bm{\mathbf{k}}}}}$ is equivalent to
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$, we find

$\displaystyle\ket{\widetilde{\phi}_{p\mathbf{k}}}=$
$\displaystyle\sum_{m}\Big{[}\ket{\psi_{m\mathbf{k}}^{(\mathrm{0})}}+\lambda{\hat{Q}^{(0)}_{\mathrm{D}}}\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}\Big{]}\widetilde{U}^{(\mathrm{0})}_{mp;{\bm{\mathbf{k}}}}+\lambda\sum_{\begin{subarray}{c}m\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}\Big{[}\left({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}\Big{]}\widetilde{U}^{(\mathrm{0})}_{mp;{\bm{\mathbf{k}}}}+{\mathcal{O}}(\lambda^{2})\,.$
(74)

Finally, we compute the perturbed WFs. We follow the scheme of Ref. Marzari
and Vanderbilt (1997), where the initial guess functions are projected onto
the disentangled subspace and then orthonormalized by Löwdin
orthogonalization. This procedure does not change the subspace but fixes the
unitary transformation within the subspace to give localized WFs. The overlap
matrix $A$ between the initial guesses and the disentangled basis states is

$\displaystyle
A_{{n\mathbf{k}},{i\mathbf{R}}}\equiv\innerproduct{\widetilde{\phi}_{n{\bm{\mathbf{k}}}}}{w_{i\mathbf{R}}^{(\mathrm{0})}}=\frac{1}{\sqrt{N_{k}}}(\widetilde{U}^{(0)\dagger}_{{\bm{\mathbf{k}}}}U^{(\mathrm{0})}_{{\bm{\mathbf{k}}}})_{ni}{e^{-i{{\bm{\mathbf{k}}}}\cdot{{\bm{\mathbf{R}}}}}}+{\mathcal{O}}(\lambda^{2}).$
(75)

The first-order term of $A_{{n\mathbf{k}},{i\mathbf{R}}}$ vanishes because we
linearly combined the disentangled basis functions so that their first-order
term is orthogonal to the unperturbed Wannier subspace [Eq. (73)]. We perform
Löwdin orthogonalization to find the gauge matrix and multiply it by
$\ket{\widetilde{\phi}_{n\mathbf{k}}}$ to obtain the WFs for the perturbed
system:

$\displaystyle\ket{w_{i\mathbf{R}}(\lambda)}=\sum_{{n\mathbf{k}}}\ket{\widetilde{\phi}_{n\mathbf{k}}}\left[A(A^{\dagger}A)^{-1/2}\right]_{{n\mathbf{k}},{i\mathbf{R}}}.$
(76)

From Eq. (75), one can show that the gauge matrix satisfies

$\displaystyle\left[A(A^{\dagger}A)^{-1/2}\right]_{{n\mathbf{k}},{i\mathbf{R}}}=A_{{n\mathbf{k}},{i\mathbf{R}}}+{\mathcal{O}}(\lambda^{2}).$
(77)

Using Eqs. (74-77), we obtain Eq. (7) and Eq. (8).

##  Appendix B Necessity of all terms in the WFP

In this section, we discuss four possible modifications to the WFP formula
[Eq. (8)] that might seem reasonable but turn out to be undesirable. The two
most important properties of the WFPs are that they accurately represent the
wavefunction perturbations and that they are spatially localized. We
demonstrate that the modifications make the WFPs lose either one of these two
properties.

Let us consider the following cases.
_Case 0_ : The true WFPs as defined in Eq. (8).
_Case 1_ : Keep only the first term in Eq. (8). This case is equivalent to
using no inner window for the Wannierization of the perturbed system, or, in
other words, ${\mathcal{W}_{\mathrm{P}}}=\emptyset$.

$\ket{w^{(1),\
\mathrm{Case1}}_{i\mathbf{R}}}=\frac{1}{\sqrt{N_{k}}}\sum_{m,{\bm{\mathbf{k}}}}{\hat{Q}^{(0)}_{\mathrm{D}}}\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi,{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(78)

_Case 2_ : Keep only the second term in Eq. (8). In other words, ignore the
contribution from states in subspace 4, which are outside the outer window
${\mathcal{W}_{\mathrm{D}}}$.

$\ket{w^{(1),\
\mathrm{Case2}}_{i\mathbf{R}}}=\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}\left({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(79)

_Case 3_ : Do not impose the condition
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$ in the second term of Eq. (8).

$\displaystyle\ket{w^{(1),\ \mathrm{Case3}}_{i\mathbf{R}}}=$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{m,{\bm{\mathbf{k}}}}{\hat{Q}^{(0)}_{\mathrm{D}}}\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi,{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}+\frac{1}{\sqrt{N_{k}}}\sum_{m,{\bm{\mathbf{k}}}}\left({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{m,{\bm{\mathbf{k}}}}\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(80)

_Case 4_ : Impose the condition ${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$
not only in the second term, but also in the first term of Eq. (8). This case
is equivalent to not using an outer window for the Wannierization of the
perturbed system to include all bands, or in other words,
${\mathcal{W}_{\mathrm{D}}}=(-\infty,\infty)$.

$\displaystyle\ket{w^{(1),\ \mathrm{Case4}}_{i\mathbf{R}}}=$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}{\hat{Q}^{(0)}_{\mathrm{D}}}\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi,{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}+\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}\left({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi^{(\mathrm{1})}_{m\mathbf{k}}}U^{(\mathrm{0})}_{mi;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
(81)

To study whether the modified WFPs can accurately and efficiently interpolate
the wavefunction perturbation, we calculate the following matrix element which
we call the “smoothed Sternheimer matrix element”:

$\displaystyle X_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\sum_{\begin{subarray}{c}p=1\\\
p{\bm{\mathbf{k+q}}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}^{\infty}\frac{(h_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu})^{*}h_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}}{\varepsilon^{\rm(0)}_{m\mathbf{k}}-\varepsilon^{\rm(0)}_{p{\bm{\mathbf{k+q}}}}}\times
f(\varepsilon^{\rm(0)}_{m\mathbf{k}}-\varepsilon^{\rm(0)}_{p{\bm{\mathbf{k+q}}}})$
$\displaystyle+$ $\displaystyle\sum_{\begin{subarray}{c}p=1\\\
p{\bm{\mathbf{k+q}}}\notin{\mathcal{W}_{\mathrm{P}}}\end{subarray}}^{\infty}\frac{(h_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu})^{*}h_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}}{\varepsilon^{\rm(0)}_{m\mathbf{k}}-\varepsilon^{\rm(0)}_{p{\bm{\mathbf{k+q}}}}},$
(82)

where

$f(x)=1-e^{-(x/\sigma)^{2}}.$ (83)

Here, we defined

$h_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{q}}}\nu}=\sum_{\kappa,a}\matrixelement{u_{m{\bm{\mathbf{k+q}}}}}{\partial_{{\bm{\mathbf{q}}}\kappa
a}\hat{H}_{\mathrm{KS}}}{u_{n\mathbf{k}}}U_{\kappa a,\nu}({\bm{\mathbf{q}}}),$
(84)

where $\partial_{{\bm{\mathbf{q}}}\kappa a}\hat{H}_{\mathrm{KS}}$ is the
derivative of the Kohn-Sham Hamiltonian with respect to a monochromatic
displacement of atom $\kappa$ along Cartesian direction $a$ with wavevector
${\bm{\mathbf{q}}}$ and $U_{\kappa a,\nu}({\bm{\mathbf{q}}})$ the displacement
pattern for phonon eigenmode $\nu$ in units of the inverse square root of
mass. The $h$ matrix element is related to the electron-phonon matrix element
$g$, which is in units of energy, as

$g_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{q}}}\nu}=\frac{1}{\sqrt{2\omega_{{\mathbf{q}\nu}}}}h_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{q}}}\nu}.$
(85)

The form of the matrix element $X_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}$ in
Eq. (B) is similar to the static Fan self-energy except that the smoothing
factor
$f(\varepsilon^{\rm(0)}_{m\mathbf{k}}-\varepsilon^{\rm(0)}_{p{\bm{\mathbf{k+q}}}})$
suppresses the contribution from states $p{\bm{\mathbf{k+q}}}$ that are
energetically close to ${m\mathbf{k}}$. Since our goal is to study whether the
contribution of the high-energy bands is accurately captured, we apply this
factor to suppress the contribution from the low-energy bands. In the
following, we set $\sigma$=0.5 eV.

Figure 10:  Band structure of silicon. Matrix elements shown in Fig. 11 are
for the states on the band denoted with the red dashed curve.  Figure 11:
Smoothed Sternheimer matrix elements for silicon obtained from an exact
calculation using the plane-wave basis and from interpolations using different
modifications of the WFP formula. The vertical dashed lines indicate $k$
points on the coarse $k$ grid where the interpolation should be exact.  Figure
12:  Spatial decay of the WFP matrix elements for silicon. Each data point
corresponds to the maximum absolute value over all WF pairs and atomic
displacements for each lattice vector ${\bm{\mathbf{R}}}$. The $y$ axes are on
a logarithmic scale. Only the ${\bm{\mathbf{q}}}=0$ case is considered.

We calculate the smoothed Sternheimer matrix elements using the WFPT
interpolation and the exact plane-wave evaluation of Eq. (B). For the WFPT
calculation, we use Eq. (36) to find

$\displaystyle X_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\sum_{\begin{subarray}{c}p=1\\\
p{\bm{\mathbf{k+q}}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}^{{N_{\mathrm{W}}}}\frac{(\widetilde{h}^{(\mathrm{H})}_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu})^{*}h^{(\mathrm{H})}_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}}{\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}}}\times
f(\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}})$
$\displaystyle+$ $\displaystyle\sum_{\begin{subarray}{c}p=1\\\
p{\bm{\mathbf{k+q}}}\notin{\mathcal{W}_{\mathrm{P}}}\end{subarray}}^{{N_{\mathrm{W}}}}\frac{(\widetilde{h}^{(\mathrm{H})}_{pm{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu})^{*}h^{(\mathrm{H})}_{pn{\bm{\mathbf{k}}};{\bm{\mathbf{q}}}\nu}}{\varepsilon^{(\mathrm{H0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}}}+s^{(\mathrm{H})}_{mn{\bm{\mathbf{k}}};{\mathbf{q}\nu}}.$
(86)

Note that $\widetilde{h}$ is equal to $h$ for states inside
${\mathcal{W}_{\mathrm{P}}}$, just as we discussed for $\widetilde{g}$ and $g$
in Sec. II.4. For the plane-wave calculation, we solve the following
Sternheimer equation to evaluate the infinite sum in the second term of Eq.
(B), in which we use the fact that
$p{\bm{\mathbf{k+q}}}\notin{\mathcal{W}_{\mathrm{P}}}$ if $p\geq 19$.

$\displaystyle\left(\hat{H}_{\mathrm{KS}}-\varepsilon_{m\mathbf{k}}\right)\ket{\delta_{{\bm{\mathbf{q}}}\kappa
a}\psi_{m\mathbf{k}}}$ $\displaystyle=$
$\displaystyle-\left(\sum_{p=19}^{\infty}\ket{\psi_{p{\bm{\mathbf{k+q}}}}}\bra{\psi_{p{\bm{\mathbf{k+q}}}}}\right)\left(\partial_{{\bm{\mathbf{q}}}\kappa
a}\hat{H}_{\mathrm{KS}}\right)\ket{\psi_{m\mathbf{k}}}$ (87)

The contribution of the bands with $p\leq 18$ to Eq. (B) is obtained by
directly calculating the eigenvalues and the matrix elements.

We calculate the smoothed Sternheimer matrix element of silicon for
${\bm{\mathbf{q}}}=0$ and $m=n$ along the band path shown in Fig. 10. Figure
11 shows the comparison for the five cases and for all six phonon modes. For
case 0, the Wannier interpolation is accurate for all $k$ points and phonon
modes. For cases 1 and 2, the Wannier interpolation does not reproduce the
exact result, even at the coarse $k$ points denoted by the vertical dashed
lines. This result indicates that the two terms in Eq. (8) are both required
to make the WFPs accurately represent the perturbation of the wavefunctions.
For cases 3 and 4, the Wannier interpolation is accurate at the coarse $k$
points, indicating that the WFPs at least represent the wavefunction
perturbations. This property can be easily shown from Eq. (B) and Eq. (B) as
the summands are exact wavefunction perturbations for
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$. However, the interpolation
displays wiggles [Figs. 11(v-ab)], which imply that the WFPs and the matrix
elements are not localized enough.

The spatial decay of the WFP matrix elements for cases 0, 3, and 4 are shown
in Fig. 12. (We did not plot the decay for cases 1 and 2 because the WFPs for
cases 1 and 2 fail to represent the wavefunction perturbation so they are of
no use even if they are spatially localized.) We find that the matrix elements
in cases 3 and 4 are much more delocalized than those in case 0, consistent
with the interpolation quality in Fig. 11.

This result can be understood as follows. For case 3, because the condition
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$ is removed, the denominator in
Eq. (II.1.2) can be close to zero making the summand a rapidly varying and
diverging function of ${\bm{\mathbf{k}}}$. Thus, the corresponding WFPs are
highly delocalized. For case 4, the outer window for Wannierization is set for
the unperturbed system but not for the perturbed system. Hence, the difference
between the unperturbed and perturbed WFs will be greater than in case 0 where
the outer window is consistently used. Therefore, the WFPs constructed with
case 4 are more delocalized than case 0.

In passing, we note that the good quality of the interpolation for the
acoustic modes for case 3 [Figs. 11(s-u)] is an exceptional case resulting
from the acoustic sum rule. Because of the translational invariance, the
matrix elements of the acoustic phonon modes at ${\bm{\mathbf{q}}}=0$ satisfy

$h_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{0}}}\nu}\sim\matrixelement{\psi_{m\mathbf{k}}}{i\commutator{\hat{H}}{\hat{p}}}{\psi_{n\mathbf{k}}}=i(\varepsilon^{\rm(0)}_{m\mathbf{k}}-\varepsilon^{\rm(0)}_{n\mathbf{k}})p_{mn{\bm{\mathbf{k}}}}$
(88)

where $p_{mn{\bm{\mathbf{k}}}}$ is the matrix element of the momentum operator
Lihm and Park (2020). Thus, when the denominator in Eq. (II.1.2) becomes close
to zero, the numerator also becomes small.

To summarize, all four modifications to the definition of the WFPs considered
in this section either fail to produce WFPs that accurately represent the
wavefunction perturbations or generate delocalized WFPs.

##  Appendix C Smoothed inner window for WFPT

In this section, we discuss how to remove the discontinuity due to the
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$ condition in Eq. (8) by smoothing
the inner window in the disentanglement step for constructing the projection-
only WFs. The issue we aim to solve is that the inner window is applied in an
“all-or-nothing” manner. So, the solution is to assign weights to the
eigenstates that are 1 for states inside the inner window and continuously
decay to 0 outside the inner window.

Concretely, we define a coefficient $l_{n\mathbf{k}}$ for each state
${n\mathbf{k}}$ inside
${\mathcal{W}_{\mathrm{F}}}-{\mathcal{W}_{\mathrm{P}}}$. The coefficient
should decay from 1 to 0 as the unperturbed eigenvalue moves away from the
boundary energy for ${\mathcal{W}_{\mathrm{P}}}$. In Appendix A, in the
disentanglement step, we projected ${\hat{P}^{(0)}_{\mathrm{W}}}$ to the
perturbed eigenstates inside
${\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}$ and diagonalized the
resulting matrix $P_{{n^{\prime}\mathbf{k^{\prime}}},{n\mathbf{k}}}$ [Eq.
(A)]. Here, instead of projecting to
${\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}$, we use the following
operator:

$\displaystyle{\hat{Q}}_{\mathrm{smooth}}=$
$\displaystyle\sum_{{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{F}}}}\ket{\psi_{n\mathbf{k}}}\bra{\psi_{n\mathbf{k}}}$
$\displaystyle+$
$\displaystyle\sum_{{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}-{\mathcal{W}_{\mathrm{P}}}}\frac{1}{1-l_{n\mathbf{k}}}\ket{\psi_{n\mathbf{k}}}\bra{\psi_{n\mathbf{k}}}.$
(89)

For each state ${n\mathbf{k}}$ inside
${\mathcal{W}_{\mathrm{F}}}-{\mathcal{W}_{\mathrm{P}}}$, the weight of the
projector $\ket{\psi_{n\mathbf{k}}}\bra{\psi_{n\mathbf{k}}}$ in
${\hat{Q}}_{\mathrm{smooth}}$ increases as
$\varepsilon^{\rm(0)}_{n\mathbf{k}}$ approaches the boundary energy of
${\mathcal{W}_{\mathrm{P}}}$ and $l_{n\mathbf{k}}$ approaches 1. Then, we
diagonalize

${\hat{P}}_{\mathrm{smooth}}={\hat{Q}}_{\mathrm{smooth}}{\hat{P}^{(0)}_{\mathrm{W}}}{\hat{Q}}_{\mathrm{smooth}}$
(90)

and take the $N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}-N_{\rm P}$ states with the
largest eigenvalues. If all $l_{n\mathbf{k}}=0$,

${\hat{Q}}_{\mathrm{smooth}}=\sum_{{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{D}}}-{\mathcal{W}_{\mathrm{P}}}}\ket{\psi_{n\mathbf{k}}}\bra{\psi_{n\mathbf{k}}}$
(91)

holds and the original scheme of Ref. Souza  _et al._ (2001) is obtained.

Now, we calculate the eigenstates of ${\hat{P}}_{\mathrm{smooth}}$ using
perturbation theory. Writing the zeroth-order eigenstates as
$\ket{\phi^{(\mathrm{0})}_{p\mathbf{k}}}$ for
$p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{D},{\bm{\mathbf{k}}}}}$,
and using Eq. (67), we have

${\hat{P}}_{\mathrm{smooth}}(\lambda=0)\ket{\phi^{(\mathrm{0})}_{p\mathbf{k}}}=\begin{dcases}\frac{\ket{\phi^{(\mathrm{0})}_{p\mathbf{k}}}}{(1-l_{p\mathbf{k}})^{2}}&:\
{N_{\mathrm{P},{\bm{\mathbf{k}}}}}<p\leq{N_{\mathrm{F},{\bm{\mathbf{k}}}}}\\\
\ket{\phi^{(\mathrm{0})}_{p\mathbf{k}}}&:\
{N_{\mathrm{F},{\bm{\mathbf{k}}}}}<p\leq{N_{\mathrm{W}}}\\\ 0&:\