11institutetext: CERN, Esplanade des Particules 1, Geneva, Switzerland
22institutetext: RWTH Aachen University, Templergraben 55, Aachen, Germany
33institutetext: DESY, NotkestraÃe 85, Hamburg, Germany

# Physics Validation of Novel Convolutional 2D Architectures for Speeding Up
High Energy Physics Simulations

Florian Rehm 1122<EMAIL_ADDRESS>Sofia Vallecorsa 11 Kerstin
Borras and Dirk KrÃcker 223333

###### Abstract

The precise simulation of particle transport through detectors remains a key
element for the successful interpretation of high energy physics results.
However, Monte Carlo based simulation is extremely demanding in terms of
computing resources. This challenge motivates investigations of faster,
alternative approaches for replacing the standard Monte Carlo approach.

We apply Generative Adversarial Networks (GANs), a deep learning technique, to
replace the calorimeter detector simulations and speeding up the simulation
time by orders of magnitude. We follow a previous approach which used three-
dimensional convolutional neural networks and develop new two-dimensional
convolutional networks to solve the same 3D image generation problem faster.
Additionally, we increased the number of parameters and the neural networks
representational power, obtaining a higher accuracy. We compare our best
convolutional 2D neural network architecture and evaluate it versus the
previous 3D architecture and Geant4 data. Our results demonstrate a high
physics accuracy and further consolidate the use of GANs for fast detector
simulations.

##  1 Introduction

Accurate Simulations of elementary particles in High Energy Physics (HEP)
detectors are fundamental to correctly reproduce and interpret the
experimental results. Detector simulations rely on Monte Carlo-based methods,
such as the Geant4 toolkit Geant4  which suffer from high computational costs:
currently more than half of the Worldwide Large Hadron Collider (LHC) Grid
resources are used for the generation and processing of simulated data
RoadmapHEP . Unfortunately, the operational requirements related to the future
High Luminosity phase of the LHC will exceed the expected available
computational resources drastically even if taking technological improvements
into account HL-LHC . The more excessive need for simulation capacities
motivates the research on alternative fast simulation approaches. Based on
generative models several prototypes have demonstrated great potential
de2017learning ; paganini2017calogan ; Salamani2018 ; dijet ; gan_lhcb . While
orders of magnitude speed ups are promised, the most crucial goal is given by
the high level of accuracy required in order to correctly reproduce the
particle showers.

Generative Adversarial Networks (GANs) represent one promising alternative
approach which have been employed primarily for simulating calorimeter
detectors EnergyGAN . In this paper we introduce three novel GAN architectures
which use Convolutional Neural Networks (CNNs) in two-dimensional (Conv2D)
layers instead of convolutional three-dimensional (Conv3D) layers for
representing 3D images. The baseline model to which we compare is taken from
Ref. EnergyGAN . Due to the large computational complexity of three-
dimensional convolutions with respect to the two-dimensional case, we expect
to observe a speed up and a reduction of the necessary computational resources
in terms of needed memory. Additionally, we increased the number of parameters
of the Conv2D neural networks what leads to a higher representation power and
therefore to a higher physics accuracy. In the end we select the best among
our new Conv2D architectures and make a detailed physics comparison to the
baseline Conv3D model and to Geant4 data.

The following section 2 provides a brief overview on related work. The GAN
approach for calorimeter simulations is introduced in section 3. In section 4,
we explain the novel Conv2D architectures, while section 5 describes the
evaluation of the results of our tests in terms of computational resources and
physics accuracy. The last section summarizes the conclusion.

##  2 Related Work

The first successful application of CNNs to image recognition problems dates
back to the 1990s with the LeNet architecture designed to classify hand-
written digits LeNet . In 2014, VGGNet VGG  was introduced: with 140M
parameters it proved for the first time, that the network depth remain a
critical aspect to achieve good performance. CNNs are extremely successful and
today employed in a broad range of various applications such as image
classification and object tracking: CSPDarkNet-53 YOLOv4  is a recent example
for object detection with YOLOv4; ResNet50 ResNet  and Inception-v4
Inception_v4  are optimised for 2D image classification. CNNs play today a key
role in many research areas and applications, ranging from medical imaging to
autonomous driving GANS_Applications .

On the other hand, the use of 2D convolutions to solve three-dimensional
problems is however limited. In 2Dslicemethod  a classification of a 3D shape
is achieved by building 3 two-dimensional slices along the three shape’s
canonical axes. Three independent images are created and classified by another
neural network. conv2duse  expanded on the previous method by stacking a set
of Conv2D layers along the $x$-, $y$\- and $z$\- axes, thus rebuilding a 3D
image from the initial set of 2D image slices. This approach resulted in
higher memory efficiency.

Within the field of High Energy Physics simulations, Torben_Thesis  proposed a
network consisting of Conv2D layers to generate the three dimensional volume
of a calorimeter: a three-dimensional 12x15x7 image is generated by stacking
12 two-dimensional 15x7 images together. This model yields accurate physics
results and an impressive $6\,660$x speed up with respect to standard Monte
Carlo techniques, running on a NVIDIA GTX 1080 card. We compare our results to
the architecture presented in Torben_Thesis .

##  3 Calorimeter Simulation with GANs

Generative Adversarial Networks (GANs) consist of two deep neural networks
trained by employing an adversarial approach. The first network, the
generator, generates images (fake images) starting from a latent vector. The
second network represents the discriminator, which evaluates and distinguishes
the generated images (fake images) from the true images, belonging to the
training set. The GAN training is successful, when the discriminator
classification probability approaches 50% for both, the true and the fake
images. Interpreting the calorimeter output as a three-dimensional image
encourages the use of architectures consisting primarily of convolutional
layers. The model developed in EnergyGAN , which uses Conv3D layers, is
building the baseline in this study and comprises 3D images of 25x25x25
calorimeter cells (pixels). An accurate representation of the energy pattern
across the pixels is key, since this pattern carries physics information used
to reconstruct the particles’ properties, such as their type.

Therefore, the generator has to reliably reproduce the particle shower
development along all three volume axes. In the following section are the
details of the alternative network approaches which we have developed and
analysed introduced. Additionally, we discuss and evaluate the corresponding
results by comparing them to the baseline Conv3D architecture and to Geant4.

##  4 New Convolutional 2D GAN Architectures

We developed three generator architectures and one discriminator architecture
using Conv2D layers. It should be noted that each of the three Conv2D GAN
models is built using the same topology for the discriminator network in order
to efficiently compare the representational power of the different Conv2D
generators to each other. The Conv2D discriminator is presented in figure 1.
Its inputs are either generated images (fake images) or training images (true
images). The discriminator outputs three values: the first is the typical GAN
true/fake probability goodfellow  which is used to calculate a binary cross
entropy loss crossentropy . The second loss (named AUX, for AUXiliary loss)
represents the result of a regression task on the initial particle energy
$E_{p}$, that the discriminator estimates from the images using a dense layer.
It is implemented as a Mean Absolute Percentage Error (MAPE) MAPE . The third
discriminator output comes from a lambda layer, calculating the sum over the
pixels of the input image which, therefore, corresponds to the total energy of
the input image. It is entitled ECAL and uses again the MAPE loss function.

Figure 1: New Conv2D discriminator network used for training all generator
networks.

###  4.1 The G1 Generator: The Three Branches Architecture

The G1 generator implements a three-branches architecture similar to the one
used for the discriminator. It outputs is a 3D 25x25x25 pixels image. The
latent space is initialised as a set of 200 random numbers, drawn from a
uniform distribution, and it is multiplied by the primary particle energy
$E_{p}$. In addition to the Conv2D layers, the generator network includes
transposed 2D convolutional (Conv2D_transpose) layers to increase the image
size, batch normalization (BatchNorm), rectified linear units activation
function (ReLU), linear ReLU activation functions (LeakyReLU) and dropout
layers (Dropout). The initial latent vector, reshaped to a 3D grid, is
transposed along the three axes and run through the three Conv2D blocks. After
transposing it back in the opposite direction, the three branches outputs are
concatenated and input to a final Conv2D layer in order to achieve the desired
dimensions (25x25x25). These transpositions allow the Conv2D network to learn
correlations within all three image dimensions. The G1 generator architecture
is shown in figure 2.

Figure 2: The G1 generator: Three branches architecture.

###  4.2 The G2 Generator: The 25 Branches Architecture

The second generator architecture (G2) follows the structure outlined in
section 2. It consists of 25 branches, one for each of the 25 layers along the
$z$ axis. The 25 branches create 25 two-dimensional images (with a 25x25
shape), that are stacked in order to rebuild a 3D 25x25x25 volume. Finally,
two additional Conv2D layers are applied to allow some interaction between the
single dimensions. The G2 generator is indicated in figure 3.

Figure 3: The G2 generator: 25 branches architecture.

###  4.3 The G3 Generator: The Single Branch Architecture

The third architecture we developed and tested, the G3 generator, is a much
deeper network. It alternates Conv2D layers to transpositions along the three
axes, in such a way that each axis is represented once as the convolutional
layers channel dimension. In order to improve convergence in the training, we
have introduced skip connections across the layers, as shown in figure 4. Due
to its depth, the G3 model training process requires a large number of epochs
before reaching convergence.

Figure 4: The G3 generator: Single branch architecture with three skip
connections.

Table 1 summarizes the total number of parameters for the three architectures:
the initial Conv3D network has less than half the number of parameters and a
smaller number of convolutional layers, compared to the three Conv2D
architectures. Furthermore, the Conv3D model consists of most of its
parameters concentrated in the first dense layer (630k out of 753k parameters
$\hat{=}$ 84%). However, for image processing convolutional layers typically
tend to have a higher representational power than dense layers while being
able to keep the number of parameter smaller. Moreover, taking into account
the Conv3D to Conv2D differences in terms of computational resources, we
expect to be able to decrease computation time by converting neural network
architectures from Conv3D to Conv2D layers. This we expect even though the
model in EnergyGAN  contains only 4 Conv3D layers compared to the 28, 177, and
25 layers in the G1, G2 and G3 architectures.

##  5 Evaluation

In this section we evaluate the performance of the three Conv2D generators. We
assess in terms of computational resources required by the inference and
training processes and the physics accuracy of the generated images. We
compare to the initial Conv3D network architecture and to Geant4.

###  5.1 Computational Evaluation

In order to measure the inference time, we process an initial set of 20 warm-
up batches and then use 100 inference steps including 20 batches for each of
the measurements. Tests were performed on a Nvidia Tesla T4 GPU using Python
version 3.6.8, TensorFlow version 2.2.0 and a batch size of 128. Table 1
summarises the inference results for each model.

Table 1: Summarizes the number of parameters and the number of convolutional layers (Conv) for each generator model. Additionally, the inference times on the GPU are shown. The speed up is given with respect to Geant4 and in bracket versus the baseline Conv3D network. The last row shows the average GPU utilization during inference. Model | Parameters | Conv | Inference | Speed Up | Utilization
---|---|---|---|---|---
Conv3D | 752k | $4$ | 7.00 s |  $6\,200$x (1.00x) | 78.75%
G1 3-branches |  $2\,052$k | $28$ | 4.91 s |  $8\,000$x (1.29x) | 21.75%
G2 25-branches |  $2\,131$k | $177$ | 17.48 s |  $2\,480$x (0.40x) | 21.83%
G3 1-branch |  $2\,397$k | $25$ | 6.21 s |  $7\,006$x (1.13x) | 32,90%

The values clearly demonstrate, that the G1 3 branches architecture has the
lowest inference time and therefore the highest speed up (factor of 1.29x)
although it has more than twice as much parameters as the Conv3D architecture.
It should be noted, that the G1 architecture has an enormous $8\,000$x speed
up versus the Geant4 simulation, aimed at being replaced by a more efficient
GAN model. Additionally, the Conv2D G1 architecture has a considerably lower
GPU utilization compared to the Conv3D baseline model (factor of 3.6x). With
running multiple streams on the GPU the inference time could be further
decreased by a potential factor of 4.64x ($=1.29\cdot 3.6$) compared to the
baseline Conv3D model. The G2 model exhibits a higher inference time than the
Conv3D model, potentially caused by the larger number of convolutional layers.
However, it shows the same low GPU utilization as the G1 architecture. The G3
model has a marginally higher inference time as the Conv3D model, but
considerably less GPU utilization.

Using the same hardware, we have investigated the training time: results per
epoch are presented in table 2. The Conv3D network is by far the slowest and
the G1 3-branches network the fastest, demonstrating a speed up of 6.5x. The
G2 25-branches network is the slowest Conv2D model.

Table 2: Training time per epoch for the different models on a Nvidia Tesla T4 GPU and the corresponding speed up compared to the Conv3D architecture. Additionally shown, are the MSE validation index for all GAN models compared to Geant4 simulation. The lower the MSE value, the better the physics accuracy and representation. Model | Time per Epoch [min] | Speed Up | MSE
---|---|---|---
Conv3D | 258 | 1x | 0.065
G1 3-branches | 40 | 6.5x | 0.027
G2 25-branches | 71 | 3.6x | 0.048
G3 1-branches | 47 | 5.5x | 0.071

###  5.2 Physics Evaluation

Evaluating the performance of a generative model for image generation remains
a non-trivial task which is mainly done with visually inspections of the
energy patterns across the calorimeter volume and comparison to Monte Carlo
predictions. An initial visual inspection demonstrates that the Conv2D and the
Conv3D model can produce realistic images, similar to Geant4 (see figure 5).

Figure 5: Example Geant4 (left), Conv2D generator (middle), Conv3D generator
(right) 3D images for a 400 GeV electron. The primary particle enters the
calorimeter from the left in the middle of the $x$\- and $y$-axis ($x=y=13$)
at $z=0$.

In order to quantify the agreement of the GAN models with Geant4, we define a
composite accuracy value. It is calculated by building 2-dimensional
projections of the particle shower distributions along the $x$-, $y$\- and
$z$-axis (averaged over $20\,000$ samples) for the GAN models and the Monte
Carlo. We measure the total mean squared error (MSE) between Geant4 and the
corresponding GAN model (summarised in table 2). The G1 architecture shows the
lowest $\mathrm{MSE}=0.027$ which corresponds to the highest accuracy of our
models. The G2 architecture has as well a lower $\mathrm{MSE}=0.048$ compared
to the Conv3D model with a $\mathrm{MSE}=0.065$. The G3 model performs worst
with a $\mathrm{MSE}=0.071$. We compared our Conv2D models with different
metrics and in all of them the G1 architecture performs best. Therefore, we
focus in the following exclusively on the G1 architecture and compare it to
the baseline Conv3D model and Geant4.

2D projections of the energy shapes along the $y$\- and $z$\- axis (averaged
over $20\,000$ samples) are presented in figure 6 in linear (left) and
logarithmic (right) scale. The $x$-axis projections are similar to the ones of
the $y$-axis and therefore not shown. In general the GAN models are close to
Geant4. A more detailed analysis of the logarithmic $y$-axis distribution (top
right), shows that the tails of the Conv3D distribution (in blue) are off,
similar as it was already observed in EnergyGAN , while the Conv2D model
(green) improves the description of the tails.

Figure 6: Particle shower distributions along the $y-$ (top) and $z-$ (bottom)
axes in linear (left) and logarithmic (right) scale. The Geant4 prediction,
the Conv3D model and the Conv2D G1 model outputs are shown in red, blue and
green respectively.

The left panel in figure 7 shows the single cell energy deposition for
$20\,000$ showers with input energies within the full energy range of 2-500
GeV. We observe, that the high energy hits are well described by both models,
even though below 1 MeV both GAN models fail to correctly describe the
expected distribution. For this reason, we apply a minimal energy threshold of
$10^{-5}$ GeV (grey shaded region) for all further validation steps.

Figure 7: Single cell energy deposition (left). Calibrated SSIM index
($L=10^{-4}$) (right), see text for further details.

In figure 7 on the right we show the Structural Similarity Index (SSIM) SSIM
between different image sets. The SSIM measures the perceptual difference
between similar images but does not judge if one image is better than the
other. The SSIM index is sensitive to the dynamic range of the pixels and must
be adjusted using the parameter $L$ in order to maintain sensitivity. We have
scanned several $L$ values ranging over orders of magnitude from $L=0.1$ to
$L=10^{-6}$. We found that the SSIM index stabilizes starting from $L=10^{-4}$
and below. We can see in figure 7 (right) that the SSIM values for Monte Carlo
(MC) of the Geant4 tool vs. MC, MC vs. G1 and MC vs. Conv3D are close to each
other, proving that the similarity between GAN generated images and Monte
Carlo images is of the same level as the similarity measured within Monte
Carlo images. However, the generated images vs. generated images SSIM is
clearly higher proving that both GAN models generate images that are more
similar to each other than what is predicted by Monte Carlo. This result
highlights a well known weakness of GAN models which tend to produce samples
exhibiting a lower diversity than the original data set. While we plan to
continue our investigation in order to better quantify GAN performance in
terms of support space size, sample diversity and mode dropping effects, the
current result is well within the precision level we expect of fast simulation
models.

In an effort to standardise a performance evaluation of generative models
applied to calorimeter simulations, we have chosen the following validation
plots following the format proposed in Buhmann:2020pmy . Figure 8 shows on the
left side the total number of hits and on the right side the total deposited
energy, for three input particle energies (50, 200 and 500 GeV). In our case,
the generation energy spectrum is continuous, so in order to probe how the
number of hits varies with the input particle energy we build $\pm 2$ GeV bins
around the corresponding central energy value: the GAN models are capable to
reproduce both the mean and the width of the Geant4 distributions.

Figure 8:  (left) Number of hits. (right) Total deposited energy. A $10^{-5}$
GeV single cell energy threshold is applied.

Figure 9 shows additional details about the deposited energy with respect to
the incident electron energy. The GAN electrons are generated using the same
input energy as the Geant4 events to minimize statistical effects. We
calculate the root-mean-square and the standard deviation of the 90% core of
the distribution, labeled with $\mu_{90}$ and $\sigma_{90}$ respectively. Both
distributions are correctly reproduced by the GANs, although a slight
degradation is visible at the start and the end of the energy range.

Figure 9: Deposited energy mean $\mu_{90}$ (left) and relative width
$\sigma_{90}/\mu_{90}$ (right). The lower panels show the relative deviation
to Geant4.

##  6 Conclusion and Future Work

This work compares three alternative Conv2D neural network architectures for
generating 3D calorimeter shower images. We focus in the further discussion on
our best Conv2D architecture (in terms of computational resources), the G1
three branches model, by comparing its performance to the Conv3D. It achieves
the largest speed up ($8\,000$x) compared to Geant4 and a factor of 1.29x
compared to the baseline Conv3D model. Additionally, our Conv2D model utilizes
the GPU with a factor of 3.29x less than the Conv3D model, leaving still the
possibility to achieve further acceleration by using multiple streams. As far
as detailed physics validation is concerned, the G1 Conv2D model achieves
similar or better performance (in the description of the distributions tails)
than the baseline Conv3D architecture developed in EnergyGAN .

In conclusion, we exhibit how a computationally simpler model, based on two
dimensional convolutions, can be used to correctly reproduce three dimensional
images. The demonstrated model in this research is not a finished software
product which will be employed exactly in this form. Nevertheless, a similar
approach, adopted to the specific calorimeter geometry, will most presumably
be used in the future HL-LHC phase to counteract the exploding demands for
simulated data. Beyond the field of High Energy Physics, this work joins the
efforts ongoing in the computer vision community to design accurate models
suitable for cases in which computational resources are limited.

## ACKNOWLEDGEMENTS

This work has been sponsored by the Wolfgang Gentner Programme of the German
Federal Ministry of Education and Research.

## References

  * (1) S. Agostinelli et al., _GEANT4–a simulation toolkit_ (2003), Vol. 506, pp. 250–303
  * (2) A. Albrecht et al., _A Roadmap for HEP Software and Computing R &D for the 2020s_ (2019), Vol. 3, https://doi.org/10.1007/s41781-018-0018-8
  * (3) G. Apollinari, B. Alonso et al., _High-Luminosity Large Hadron Collider (HL-LHC): Technical Design Report V. 0.1_ (2017), Vol. 4/2017
  * (4) de Oliveira et al., _Learning Particle Physics by Example: Location-Aware Generative Adversarial Networks for Physics Synthesis_ (Springer Science and Business Media LLC, 2017), Vol. 1
  * (5) M. Paganini et al., _CaloGAN: Simulating 3D high energy particle showers in multilayer electromagnetic calorimeters with generative adversarial networks_ (American Physical Society (APS), 2018), Vol. 97
  * (6) A. Ghosh (ATLAS Collaboration), Tech. Rep. ATL-SOFT-PROC-2019-007, CERN, Geneva (2019), http://cds.cern.ch/record/2680531
  * (7) D. Sipio et al., _DijetGAN: a Generative-Adversarial Network approach for the simulation of QCD dijet events at the LHC_ (Springer Science and Business Media, 2019)
  * (8) C. et al, _Generative Models for Fast Calorimeter Simulation: the LHCb case_ (EDP Sciences, 2019), p. 02034, http://dx.doi.org/10.1051/epjconf/201921402034
  * (9) G. Khattak et al., _Three Dimensional Energy Parametrized Generative Adversarial Networks for Electromagnetic Shower Simulation_ (2018)
  * (10) Y. Lecun, L. Bottou, Y. Bengio, P. Haffner, _Gradient-based learning applied to document recognition_ , in _Proceedings of the IEEE_ (1998), pp. 2278–2324
  * (11) S. Liu, W. Deng, _Very deep convolutional neural network based image classification using small training sample size_ , in _2015 3rd IAPR Asian Conference on Pattern Recognition (ACPR)_ (2015), pp. 730–734
  * (12) A. Bochkovskiy, C.Y. Wang, H.y. Liao, _Yolov4: Optimal speed and accuracy of object detection_ (2020)
  * (13) K. He, X. Zhang, S. Ren, J. Sun, _Deep residual learning for image recognition_ (2015)
  * (14) C. Szegedy, S. Ioffe, V. Vanhoucke, A. Alemi, _Inception-v4, inception-resnet and the impact of residual connections on learning_ (2016), 1602.07261
  * (15) Alqahtani et al., _Applications of Generative Adversarial Networks (GANs): An Updated Review_ (2019)
  * (16) F. Gomez-Donoso, A. Garcia-Garcia, J. RodrÃguez, S. Orts, M. Cazorla, _LonchaNet: A sliced-based CNN architecture for real-time 3D object recognition_ (2017), pp. 412–418
  * (17) K. Sarkar, B. Hampiholi, K. Varanasi, D. Stricker, _Learning 3d shapes as multi-layered height-maps using 2d convolutional networks_ (2018), 1807.08485
  * (18) T. Quast, _Qualification, performance validation and fast generative modelling of beam test calorimeter prototypes for the cms calorimeter endcap upgrade_ (2020)
  * (19) I.J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, _Generative adversarial networks_ (2014), 1406.2661
  * (20) G.E. Nasr, E. Badr, C. Joun, _Cross Entropy Error Function in Neural Networks: Forecasting Gasoline Demand_ , in _FLAIRS Conference_ (2002)
  * (21) P.M. Swamidass, ed., _MAPE (mean absolute percentage error)_ (Springer US, Boston, MA, 2000), pp. 462–462, https://doi.org/10.1007/1-4020-0612-8_580
  * (22) J. Nilsson, T. Akenine-MÃller, _Understanding ssim_ (2020)
  * (23) E. Buhmann, _Getting High: High Fidelity Simulation of High Granularity Calorimeters with High Speed_ (2020)


{N_{\mathrm{W}}}<p\leq{N_{\mathrm{D},{\bm{\mathbf{k}}}}}.\end{dcases}$ (92)

The total number of states with nonzero unperturbed eigenvalues of
${\hat{P}}_{\mathrm{smooth}}$ is $N_{\bm{\mathbf{k}}}{N_{\mathrm{W}}}-N_{\rm
P}$, and the perturbed version of these states will be selected for
Wannierization. Since the mixing between the selected states is not important
because our goal here is to determine the subspace spanned by them, we only
consider the mixing between the eigenstates of
${\hat{P}}_{\mathrm{smooth}}(\lambda=0)$ with nonzero eigenvalues and those
with zero eigenvalues.

Using first-order perturbation theory in a similar manner to Eq. (A), the
perturbed eigenstates of ${\hat{P}}_{\mathrm{smooth}}$ are

$\displaystyle\ket{\phi_{p\mathbf{k}}}-\ket{\phi^{(\mathrm{0})}_{p\mathbf{k}}}=\begin{dcases}-\lambda(1-l_{p\mathbf{k}})\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})&\text{for
}p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{F},{\bm{\mathbf{k}}}}}\\\
-\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})&\text{for
}p={N_{\mathrm{F},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{W}}}\end{dcases}$
(93)

The only difference with Eq. (A) is the $\left(1-l_{p\mathbf{k}}\right)$
factor for
$p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{F},{\bm{\mathbf{k}}}}}$.

Finally, we include the states inside ${\mathcal{W}_{\mathrm{P}}}$
($\ket{\phi^{(\mathrm{0})}_{p\mathbf{k}}}$ for
$p=1,\cdots,{N_{\mathrm{P},{\bm{\mathbf{k}}}}}$) and orthogonalize them to
$\mathcal{O}(\lambda)$ as done in Eq. (73). With Eq. (67), the resulting basis
states for the disentangled subspace are

$\displaystyle\ket{\widetilde{\phi}_{p\mathbf{k}}}=\begin{dcases}\ket{\psi_{p\mathbf{k}}^{(\mathrm{0})}}+\lambda\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{W}}}\right)\ket{\psi_{p\mathbf{k}}^{(\mathrm{1})}}+{\mathcal{O}}(\lambda^{2})&\text{for
}p=1,\cdots,{N_{\mathrm{P},{\bm{\mathbf{k}}}}},\\\
\ket{\psi_{p\mathbf{k}}^{(\mathrm{0})}}+\lambda\left(\hat{\mathbb{1}}-{\hat{P}^{(0)}_{\mathrm{D}}}+l_{p\mathbf{k}}({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}})\right)\ket{\psi_{p\mathbf{k}}^{(\mathrm{1})}}+{\mathcal{O}}(\lambda^{2})&\text{for
}p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{F},{\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{F},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{W}}}.\end{dcases}$
(94)

Note the additional
$l_{p\mathbf{k}}({\hat{P}^{(0)}_{\mathrm{D}}}-{\hat{P}^{(0)}_{\mathrm{W}}})$
factor for the
$p={N_{\mathrm{P},{\bm{\mathbf{k}}}}}+1,\cdots,{N_{\mathrm{F},{\bm{\mathbf{k}}}}}$
case, which is absent in Eq. (73).

By projecting and orthogonalizing the unperturbed WFs to the disentangled
subspace, the WFPs become

$\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}}}}+\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+\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}-{\mathcal{W}_{\mathrm{P}}}\end{subarray}}l_{m\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}}}}.$
(95)

The only difference between Eq. (C) and the original formula for WFPs [Eq.
(8)] is that Eq. (C) has the third term. The factor $l_{m\mathbf{k}}$ in this
term smoothly changes from 1 to 0 as $\varepsilon^{\rm(0)}_{m\mathbf{k}}$
moves from the boundary energy of ${\mathcal{W}_{\mathrm{P}}}$ to that of
${\mathcal{W}_{\mathrm{F}}}$. Hence, this third term eliminates the
discontinuity generated by the second term of Eq. (C) while maintaining the
desired property that the WFPs should accurately represent the wavefunction
perturbation of the states inside ${\mathcal{W}_{\mathrm{P}}}$. Therefore, Eq.
(C) is an alternative way to create WFPs and is expected to give more
localized WFPs than Eq. (8). We detail how the real-space matrix elements can
be calculated with this modified definition of WFPs in Appendix F.2.

For the systems studied in this work, we find that the WFPs constructed
without smoothing the inner window are already sufficiently localized.
Smoothing the inner window will be useful in systems where the WFP
delocalization due to the discontinuity at the boundary of the inner window is
a problem.

##  Appendix D WFP correction to the Hamiltonian matrix element

In this appendix, we prove Eq. (25).

$\displaystyle H^{(\mathrm{\rm
W1})}_{i^{\prime}{\bm{\mathbf{k^{\prime}}}},i{\bm{\mathbf{k}}}}=$
$\displaystyle\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}},{\bm{\mathbf{R^{\prime}}}},{\bm{\mathbf{q}}}}\Bigg{[}\matrixelement{w^{(\mathrm{0})}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}}}}{\hat{V}_{\bm{\mathbf{q}}}}{w^{(\mathrm{0})}_{i{\bm{\mathbf{R}}}}}+\matrixelement{w^{(\mathrm{0})}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{1})}_{{i\mathbf{R}};{\bm{\mathbf{q}}}}}+\matrixelement{w^{(\mathrm{1})}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}};-{\bm{\mathbf{q}}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{0})}_{i\mathbf{R}}}\Bigg{]}{e^{-i{{\bm{\mathbf{k^{\prime}}}}}\cdot{{\bm{\mathbf{R^{\prime}}}}}}}e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle=$
$\displaystyle\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}},{\bm{\mathbf{R^{\prime}}}},{\bm{\mathbf{q}}}}\Bigg{[}\matrixelement{w^{(\mathrm{0})}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}}}}{\hat{V}_{\bm{\mathbf{q}}}}{w^{(\mathrm{0})}_{i{\bm{\mathbf{R+R^{\prime}}}}}}+\matrixelement{w^{(\mathrm{0})}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{1})}_{i{\bm{\mathbf{R+R^{\prime}}}};{\bm{\mathbf{q}}}}}+\matrixelement{w^{(\mathrm{1})}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}};-{\bm{\mathbf{q}}}}}{\hat{H}^{(\mathrm{0})}}{w^{(\mathrm{0})}_{i{\bm{\mathbf{R+R^{\prime}}}}}}\Bigg{]}{e^{-i{{\bm{\mathbf{k^{\prime}}}}}\cdot{{\bm{\mathbf{R^{\prime}}}}}}}e^{i{\bm{\mathbf{k}}}\cdot({\bm{\mathbf{R+R^{\prime}}}})}$
$\displaystyle=$
$\displaystyle\frac{1}{N_{k}}\sum_{{\bm{\mathbf{R}}},{\bm{\mathbf{R^{\prime}}}},{\bm{\mathbf{q}}}}\left(g_{i^{\prime}i{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}+\delta
g_{i^{\prime}i{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\right)e^{i{\bm{\mathbf{q}}}\cdot{\bm{\mathbf{R^{\prime}}}}}{e^{-i{{\bm{\mathbf{k^{\prime}}}}}\cdot{{\bm{\mathbf{R^{\prime}}}}}}}e^{i{\bm{\mathbf{k}}}\cdot({\bm{\mathbf{R+R^{\prime}}}})}$
$\displaystyle=$
$\displaystyle\sum_{{\bm{\mathbf{R}}},{\bm{\mathbf{q}}}}\left(g_{i^{\prime}i{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}+\delta
g_{i^{\prime}i{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}\right)e^{i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}\delta_{{\bm{\mathbf{k}}}+{\bm{\mathbf{q}}},{\bm{\mathbf{k^{\prime}}}}}$
$\displaystyle=$
$\displaystyle\widetilde{g}^{\mathrm{(W)}}_{i^{\prime}i{\bm{\mathbf{k}}};{\bm{\mathbf{k^{\prime}-k}}}}$
(96)

##  Appendix E Interpolation of matrix elements for Kubo formula

In this section, we discuss the interpolation of the $K_{mn{\bm{\mathbf{k}}}}$
matrix element defined in Eq. (37). Wannier interpolation of
$K_{mn{\bm{\mathbf{k}}}}$ can be done analogously to that of
$S_{mn{\bm{\mathbf{k}}}}$ in Eqs. (33-36). For
${n\mathbf{k}},{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$, one can use the
Wannier interpolation of the wavefunction perturbation [Eq. (II.2.3)] to find

$\displaystyle K_{mn{\bm{\mathbf{k}}}}=$
$\displaystyle\innerproduct{\psi^{(\mathrm{\rm
H1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{\psi^{\prime\rm(H1)}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}$
$\displaystyle=$
$\displaystyle\frac{1}{N_{k}}\sum_{i,i^{\prime},{\bm{\mathbf{R}}},{\bm{\mathbf{R^{\prime}}}}}e^{i{\bm{\mathbf{k}}}\cdot({\bm{\mathbf{R}}}-{\bm{\mathbf{R^{\prime}}}})}V^{\dagger{(\mathrm{0})}}_{mi^{\prime};{\bm{\mathbf{k}}}}\innerproduct{w^{(1)}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}};{\bm{\mathbf{q}}}}}{w^{(1)}_{i^{\prime}{\bm{\mathbf{R^{\prime}}}};{\bm{\mathbf{q}}}}}{w^{\prime(1)}_{{i\mathbf{R}};{\bm{\mathbf{q}}}}}V^{(\mathrm{0})}_{in;{\bm{\mathbf{k}}}}+\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{\bm{\mathbf{k+q}}}})(\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}})}$
$\displaystyle=$
$\displaystyle\sum_{i,i^{\prime},{\bm{\mathbf{R}}}}V^{\dagger{(\mathrm{0})}}_{mi^{\prime};{\bm{\mathbf{k}}}}k_{i^{\prime}i{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}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}^{\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{\bm{\mathbf{k+q}}}})(\varepsilon^{(\mathrm{H0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{H0})}_{p{\bm{\mathbf{k+q}}}})}.$
(97)

In the third equality, we used the discrete translation properties of the WFP
[Eq. (16)] and renamed ${\bm{\mathbf{R}}}-{\bm{\mathbf{R^{\prime}}}}$ to
${\bm{\mathbf{R}}}$. We also used the real-space matrix element defined in Eq.
(38).

Using Eq. (41) and Eq. (40), one can rewrite Eq. (E) as Eq. (39). Since all
quantities appearing in Eq. (39) can be calculated using a WF- or WFP-based
interpolation, $K_{mn{\bm{\mathbf{k}}}}$ can also be interpolated.

##  Appendix F Calculation of real-space matrix elements

###  F.1 Without the smoothed inner window

In this appendix, we describe how to calculate the real-space matrix elements
using the Bloch eigenstates computed at a coarse $k$ grid. For the potential
matrix element $g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, Ref. Giustino  _et
al._ (2007) gives

$g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=\frac{1}{N_{k}}\sum_{m,n,{\bm{\mathbf{k}}}}U^{(0)\dagger}_{im;{\bm{\mathbf{k+q}}}}\matrixelement{\psi^{(0)}_{m\mathbf{k}+\mathbf{q}}}{\hat{V}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}},$
(98)

which easily follows from the definition of
$g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ [Eq. (17)]. Here, we derive
analogous expressions for the newly introduced matrices $\delta
g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ and
$s_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$.

First, for $\delta g_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, we start from its
definition [Eq. (18)] and use Eqs. (1, 8) to find

$\displaystyle\delta g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\frac{1}{N_{k}}\sum_{\begin{subarray}{c}m,n,p,{\bm{\mathbf{k}}}\\\
{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}U^{(0)\dagger}_{im;{\bm{\mathbf{k+q}}}}\frac{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}+\mathbf{q}}(P^{(\mathrm{0})}_{\mathrm{D}}-P^{(\mathrm{0})}_{\mathrm{W}})_{mp;{\bm{\mathbf{k+q}}}}\matrixelement{\psi^{(0)}_{p{\bm{\mathbf{k+q}}}}}{\hat{V}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle+$
$\displaystyle\frac{1}{N_{k}}\sum_{\begin{subarray}{c}m,n,p,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}+\mathbf{q}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}U^{(0)\dagger}_{im;{\bm{\mathbf{k+q}}}}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k}+\mathbf{q}}}{\hat{V}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{p{\bm{\mathbf{k}}}}}(P^{(\mathrm{0})}_{\mathrm{D}}-P^{(\mathrm{0})}_{\mathrm{W}})_{pn;{\bm{\mathbf{k}}}}\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}+\mathbf{q}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}.$
(99)

In Eq. (F.1), all band indices are limited to the bands inside
${\mathcal{W}_{\mathrm{D}}}$. The
${\hat{Q}^{(0)}_{\mathrm{D}}}\ket{\psi_{m\mathbf{k}}^{(\mathrm{1})}}$ term in
Eq. (8) does not appear because the WFs are orthogonal to the
${\hat{Q}^{(0)}_{\mathrm{D}}}$ subspace. Thus, $\delta
g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$ can be evaluated using only the
energies and matrix elements for the states inside
${\mathcal{W}_{\mathrm{D}}}$.

For $s_{{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, using the definitions [Eqs. (1,
8, 33)], we find

$\displaystyle s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\frac{1}{N_{k}}\sum_{m,n,{\bm{\mathbf{k}}}}U^{(0)\dagger}_{im,{\bm{\mathbf{k}}}}\matrixelement{\psi^{(\mathrm{1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{{\hat{Q}^{(0)}_{\mathrm{D}}}\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}U^{(\mathrm{0})}_{nj,{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle+$
$\displaystyle\frac{1}{N_{k}}\sideset{}{{}^{\prime}}{\sum}_{{\begin{subarray}{c}m,n,p,q,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}}U^{(0)\dagger}_{im,{\bm{\mathbf{k}}}}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k}}}{\hat{V}_{\bm{\mathbf{q}}}^{\dagger}}{\psi^{(0)}_{p\mathbf{k}+\mathbf{q}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p\mathbf{k}+\mathbf{q}}}(P^{(\mathrm{0})}_{\rm
D}-P^{(\mathrm{0})}_{{\rm
W}})_{pq;{\bm{\mathbf{k+q}}}}\matrixelement{\psi^{(0)}_{q{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}.$
(100)

The first term involves the wavefunction perturbation, which can be calculated
by solving the Sternheimer equation. The second term refers only to the bands
inside the outer window.

For $k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}$, we similarly find

$\displaystyle k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\sum_{m,n,{\bm{\mathbf{k}}}}U^{(0)\dagger}_{im;{\bm{\mathbf{k}}}}\matrixelement{\psi^{(\mathrm{1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{{\hat{Q}^{(0)}_{\mathrm{D}}}}{\psi^{\prime(1)}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}{e^{-i{{\bm{\mathbf{k}}}}\cdot{{\bm{\mathbf{R}}}}}}$
$\displaystyle+$
$\displaystyle\frac{1}{N_{k}}\sum_{\begin{subarray}{c}m,n,p,q,{\bm{\mathbf{k}}}\\\
{n\mathbf{k}},{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}\end{subarray}}U^{(0)\dagger}_{im,{\bm{\mathbf{k}}}}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k}}}{\hat{V}_{\bm{\mathbf{q}}}^{\dagger}}{\psi^{(0)}_{p\mathbf{k}+\mathbf{q}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p\mathbf{k}+\mathbf{q}}}(P^{(\mathrm{0})}_{{\rm
D}}-P^{(\mathrm{0})}_{{\rm
W}})_{pq;{\bm{\mathbf{k+q}}}}\frac{\matrixelement{\psi^{(0)}_{q{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{q{\bm{\mathbf{k+q}}}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}.$
(101)

###  F.2 With the smoothed inner window

In this subsection, we detail how the formulas for the real-space matrix
elements change if the WFPs are constructed with the smoothed inner window as
described in Appendix C. For brevity, we define

$\widetilde{l}_{m\mathbf{k}}=\begin{cases}1&\text{:
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$}\\\ l_{m\mathbf{k}}&\text{:
${m\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}-{\mathcal{W}_{\mathrm{P}}}$}\end{cases}$
(102)

and rewrite Eq. (C) as

$\displaystyle\ket{w^{(1)}_{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}}}}+\frac{1}{\sqrt{N_{k}}}\sum_{\begin{subarray}{c}m,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}\end{subarray}}\widetilde{l}_{m\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}}}}.$
(103)

Then, one finds that the real-space matrix elements are

$\displaystyle\delta g_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\frac{1}{N_{k}}\sum_{\begin{subarray}{c}m,n,p,{\bm{\mathbf{k}}}\\\
{n\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}\end{subarray}}U^{(0)\dagger}_{im;{\bm{\mathbf{k+q}}}}\frac{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}+\mathbf{q}}(P^{(\mathrm{0})}_{\mathrm{D}}-P^{(\mathrm{0})}_{\mathrm{W}})_{mp;{\bm{\mathbf{k+q}}}}\matrixelement{\psi^{(0)}_{p{\bm{\mathbf{k+q}}}}}{\hat{V}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k+q}}}}}\widetilde{l}_{n\mathbf{k}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle+$
$\displaystyle\frac{1}{N_{k}}\sum_{\begin{subarray}{c}m,n,p,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}+\mathbf{q}}\in{\mathcal{W}_{\mathrm{F}}}\end{subarray}}U^{(0)\dagger}_{im;{\bm{\mathbf{k+q}}}}\widetilde{l}_{m\mathbf{k}+\mathbf{q}}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k}+\mathbf{q}}}{\hat{V}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{p{\bm{\mathbf{k}}}}}(P^{(\mathrm{0})}_{\mathrm{D}}-P^{(\mathrm{0})}_{\mathrm{W}})_{pm;{\bm{\mathbf{k}}}}\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}+\mathbf{q}}-\varepsilon^{(\mathrm{0})}_{p{\bm{\mathbf{k}}}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}},$
(104) $\displaystyle s_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\frac{1}{N_{k}}\sum_{m,n,{\bm{\mathbf{k}}}}U^{(0)\dagger}_{im,{\bm{\mathbf{k}}}}\matrixelement{\psi^{(\mathrm{1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{{\hat{Q}^{(0)}_{\mathrm{D}}}\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}U^{(\mathrm{0})}_{nj,{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}$
$\displaystyle+$
$\displaystyle\frac{1}{N_{k}}\sideset{}{{}^{\prime}}{\sum}_{{\begin{subarray}{c}m,n,p,q,{\bm{\mathbf{k}}}\\\
{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}\end{subarray}}}U^{(0)\dagger}_{im,{\bm{\mathbf{k}}}}\widetilde{l}_{m\mathbf{k}}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k}}}{\hat{V}_{\bm{\mathbf{q}}}^{\dagger}}{\psi^{(0)}_{p\mathbf{k}+\mathbf{q}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p\mathbf{k}+\mathbf{q}}}(P^{(\mathrm{0})}_{\rm
D}-P^{(\mathrm{0})}_{{\rm
W}})_{pq;{\bm{\mathbf{k+q}}}}\matrixelement{\psi^{(0)}_{q{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}},$
(105)

and

$\displaystyle k_{ij{\bm{\mathbf{R}}};{\bm{\mathbf{q}}}}=$
$\displaystyle\sum_{m,n,{\bm{\mathbf{k}}}}U^{(0)\dagger}_{im;{\bm{\mathbf{k}}}}\matrixelement{\psi^{(\mathrm{1})}_{{m\mathbf{k}};{\bm{\mathbf{q}}}}}{{\hat{Q}^{(0)}_{\mathrm{D}}}}{\psi^{\prime(1)}_{{n\mathbf{k}};{\bm{\mathbf{q}}}}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}{e^{-i{{\bm{\mathbf{k}}}}\cdot{{\bm{\mathbf{R}}}}}}$
$\displaystyle+$
$\displaystyle\frac{1}{N_{k}}\sum_{\begin{subarray}{c}m,n,p,q,{\bm{\mathbf{k}}}\\\
{n\mathbf{k}},{m\mathbf{k}}\in{\mathcal{W}_{\mathrm{F}}}\end{subarray}}U^{(0)\dagger}_{im,{\bm{\mathbf{k}}}}\widetilde{l}_{m\mathbf{k}}\frac{\matrixelement{\psi^{(0)}_{m\mathbf{k}}}{\hat{V}_{\bm{\mathbf{q}}}^{\dagger}}{\psi^{(0)}_{p\mathbf{k}+\mathbf{q}}}}{\varepsilon^{(\mathrm{0})}_{m\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{p\mathbf{k}+\mathbf{q}}}(P^{(\mathrm{0})}_{{\rm
D}}-P^{(\mathrm{0})}_{{\rm
W}})_{pq;{\bm{\mathbf{k+q}}}}\frac{\matrixelement{\psi^{(0)}_{q{\bm{\mathbf{k+q}}}}}{\hat{V}^{\prime}_{\bm{\mathbf{q}}}}{\psi^{(0)}_{n\mathbf{k}}}}{\varepsilon^{(\mathrm{0})}_{n\mathbf{k}}-\varepsilon^{(\mathrm{0})}_{q{\bm{\mathbf{k+q}}}}}\widetilde{l}_{n\mathbf{k}}U^{(\mathrm{0})}_{nj;{\bm{\mathbf{k}}}}e^{-i{\bm{\mathbf{k}}}\cdot{\bm{\mathbf{R}}}}.$
(106)

When $l_{p\mathbf{k}}=0$, which means that the inner window is not smoothed at
all, Eq. (F.2), Eq. (F.2), and Eq. (F.2) reduce to Eq. (F.1), Eq. (F.1), and
Eq. (F.1), respectively.

##  Appendix G Improved double-grid method for phonon-induced electron self-
energy

In this section, we discuss the improved double-grid method for accelerating
the convergence of the phonon-induced electron self-energy with respect to the
coarse $k$ grid size. Here, we focus on the diagonal part of the self-energy,
but the results can also be applied to the off-diagonal one Lihm and Park
(2020). Our aim is to calculate $\Sigma_{n\mathbf{k}}$ for all
${n\mathbf{k}}\in{\mathcal{W}_{\mathrm{P}}}$.

The convergence of the self-energy with respect to the $q$-point sampling is
slow. This convergence rate is dominated by the convergence of the lower Fan
self-energy Poncé  _et al._ (2015). The convergence of the upper Fan term with
respect to the $q$-point sampling is much faster because the denominator
$\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}+\mathbf{q}}$ in Eq. (50)
is never small. Therefore, one can save considerable computational cost by
calculating the rapidly convergent upper Fan self-energy at a coarse $q$-point
grid and the slowly convergent lower Fan self-energy at a finer $q$-point
grid.

$\Sigma_{n\mathbf{k}}(\eta)=\sum_{{\bm{\mathbf{q}}},\nu}^{\mathrm{coarse}}\left[\Sigma^{\rm
upper\ Fan}_{n\mathbf{k}}+\Sigma^{\rm
DW}_{n\mathbf{k}}\right]+\sum_{{\bm{\mathbf{q}}},\nu}^{\mathrm{fine}}\left[\Sigma^{\rm
lower\ Fan}_{n\mathbf{k}}(\eta)\right]$ (107)

This method is the double-grid method Brown-Altvater  _et al._ (2020). To
distinguish this method from the improved version we propose later, we call it
the “conventional double-grid method.” Note that one should also calculate the
DW self-energy at the coarse grid when using the double-grid method. The
reason for this is that the convergence of the sum of the DW and the upper Fan
self-energies is much faster than the convergence of each term because of the
approximate cancellation between the two terms Allen and Heine (1976);
Querales-Flores  _et al._ (2019); Lihm and Park (2020).

However, we found that the conventional double-grid method is not very
successful in the case of the finite-temperature self-energy of silicon.
Figure 13 shows that the sum of the upper Fan and DW self-energy (black dots)
converges as $1/N_{q}$ if an $N_{q}\times N_{q}\times N_{q}$ $q$-point grid is
used. When a $6\times 6\times 6$ coarse grid is used, the error in the band-
gap renormalization at 400 K is 14 meV, around 13% of the converged value.
Therefore, we need to speed up the convergence.

Figure 13:  Phonon-induced electron self-energy of silicon at $T$=400 K
calculated using the conventional and improved double-grid (DG) methods for
the (a) VBM, (b) CBM, and the lowest conduction band (CB) at $\Gamma$. The
fine $q$ grid size and $\eta$ for the self-energy are fixed to $125\times
125\times 125$ and 100 meV, respectively.  Figure 14:  (a) Contribution of
each phonon mode to the sum of the DW and upper Fan self-energies of silicon
for the electron state at the valence band maximum. The size and color of the
circles indicate the magnitude and sign of the phonon-mode-resolved electron
self-energy, where red (blue) means a positive (negative) contribution. (b)
Sum of the upper Fan (UF), DW, and optionally lower Fan (LF) self-energies for
the lowest-frequency ($\nu$=1) acoustic mode multiplied by
$\absolutevalue{{\bm{\mathbf{q}}}}^{2}$ for $q$ points along the
$\mathrm{\Gamma}$-L line.

By investigating the contribution of each phonon mode to the self-energy, we
find that the contribution of the low-energy acoustic phonons dominates the
slow convergence rate as shown in Fig. 14(a). The slow $1/N_{q}$ convergence
can be understood as follows. Let us focus on the acoustic modes at small
$|{\bf q}|$ satisfying $\omega_{\mathbf{q}\nu}\approx
v\absolutevalue{\mathbf{q}}$. We set the momentum cutoff $q_{\rm cut}\ll
2\pi/a$ for the low-energy modes where $a$ is the lattice parameter. We assume
that the temperature is high enough to satisfy $vq_{\rm cut}\ll
k_{\mathrm{B}}T$ so that $2n_{\mathbf{q}\nu}+1\approx
2k_{\mathrm{B}}T/\hbar\omega_{\mathbf{q}\nu}$. To make the $q$ dependence
explicit, we write $g_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{q}}}\nu}$ in terms of
$h_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{q}}}\nu}$ using Eq. (85) and approximate
$h_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{q}}}\nu}$ by
$h_{mn{\bm{\mathbf{k}}},{\bm{\mathbf{0}}}\nu}$. Note that $h$, not $g$, should
be approximated because $g$ has an additional $|{\bf q}|$ dependence in the
$1/\sqrt{\omega_{\mathbf{q}\nu}}$ factor in Eq. (85). Then, the sum of Eq.
(50) for $\absolutevalue{\mathbf{q}}<q_{\rm cut}$ becomes

$\displaystyle\Sigma^{\rm upper\ Fan,\ acoustic}_{n{\bm{\mathbf{k}}}}\approx$
$\displaystyle\sum_{m\notin{\mathcal{W}_{\mathrm{P}}}}\frac{\absolutevalue{h_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{0}}}\nu}}^{2}}{\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}}}\int_{\absolutevalue{\mathbf{q}}<q_{\rm
cut}}d{\bm{\mathbf{q}}}\frac{k_{\mathrm{B}}T}{\hbar
v^{2}\absolutevalue{\mathbf{q}}^{2}}$ $\displaystyle=$ $\displaystyle q_{\rm
cut}\frac{4\pi k_{\mathrm{B}}T}{\hbar
v^{2}}\sum_{m\notin{\mathcal{W}_{\mathrm{P}}}}\frac{\absolutevalue{h_{mn{\bm{\mathbf{k}}};{\bm{\mathbf{0}}}\nu}}^{2}}{\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}}}.$
(108)

The acoustic phonon contribution to the upper Fan self-energy is proportional
to $q_{\rm cut}$. Using a coarse $q$ grid gives an effective cutoff $q_{\rm
cut}\sim 1/N_{q}$ because the acoustic phonon modes with $q<1/N_{q}\times
2\pi/a$ are not sampled. Therefore, the conventional double-grid method
converges as $1/N_{q}$. Note that the convergence is much faster at $T=0$
because then $2n_{\mathbf{q}\nu}+1=1$ does not give a $1/q$ factor, so
$\Sigma^{\rm upper\ Fan,\ acoustic}_{n{\bm{\mathbf{k}}}}$ becomes proportional
to $q_{\rm cut}^{2}\sim 1/N_{q}^{2}$.

To speed up the convergence, we utilize the fact that the acoustic modes
around $\mathbf{q}=0$ actually have little effect on the total self-energy
because of the translational symmetry. Since the acoustic modes at
$\mathbf{q}=0$ describe a uniform translation of the whole lattice, the band
eigenvalues should not change. However, the sum of the phonon-mode-resolved
upper Fan term and the DW term is not small and even diverges at
$\mathbf{q}=0$. This divergence is canceled by the divergence of the lower Fan
term in the opposite sign, so their sum, the total self-energy, is zero at
$\mathbf{q}=0$.

To make use of this cancellation, we rewrite the total self-energy as follows.

$\displaystyle\Sigma_{n\mathbf{k}}(\eta)=$
$\displaystyle\sum_{{\bm{\mathbf{q}}},\nu}^{\mathrm{coarse}}\left[\Sigma^{\rm
upper\ Fan}_{n\mathbf{k}}+\Sigma^{\rm DW}_{n\mathbf{k}}+\Sigma^{\rm lower\
Fan}_{n\mathbf{k}}(\eta_{\mathrm{large}})\right]$ $\displaystyle+$
$\displaystyle\sum_{{\bm{\mathbf{q}}},\nu}^{\mathrm{fine}}\left[\Sigma^{\rm
lower\ Fan}_{n\mathbf{k}}(\eta)-\Sigma^{\rm lower\
Fan}_{n\mathbf{k}}(\eta_{\mathrm{large}})\right]$ (109)

Here, $\eta$ is the small parameter that should be converged to 0, while
$\eta_{\mathrm{large}}$ is a large value around 1 eV and does not need to be
converged. Equation Eq. (G) is the “improved double-grid method” we propose in
this work. The key property of Eq. (G) is that the first term converges much
faster with respect to the $q$ grid than the first term of Eq. (107), which is
just the sum of the upper Fan and DW self-energies. The reason for this is
that the diverging contribution of the low-energy acoustic modes in the upper
Fan and DW terms is suppressed by the additional lower Fan term. Figure 14(b)
shows the self-energy contributions multiplied by the phase-space volume
factor $\absolutevalue{\mathbf{q}}^{2}$. We find that the contribution of
modes with small $\absolutevalue{\mathbf{q}}$ is indeed suppressed by adding
the lower Fan term.

The cancellation becomes exact as $\eta_{\mathrm{large}}$ is decreased.
However, we need to choose $\eta_{\mathrm{large}}$ to be large enough so that
the denominator in the lower Fan term [Eq. (III.1)] is always far from zero
and does not create additional convergence issues. For example, Fig. 13(c)
shows a slow convergence of the self-energy for the lowest conduction band at
$\mathrm{\Gamma}$ for $\eta_{\mathrm{large}}=0.5$ eV. The slow convergence is
more pronounced for states far from band extrema because there are more
resonant couplings mediated by phonons to other electron states, leading to a
rapid variation in the phonon-mode-resolved lower Fan self-energy.
Practically, we find that $\eta_{\mathrm{large}}$ in [1 eV, 2 eV] is a good
choice for silicon.

Figure 13 shows that the total self-energy converges much faster with respect
to the coarse $q$ grid when using the improved double-grid method. By using
the improved double-grid method with $\eta_{\mathrm{\large}}=1~{}\mathrm{eV}$,
the absolute (relative) error in the band-gap renormalization at 400 K is only
2.4 meV (2%) with a $6\times 6\times 6$ coarse $q$ grid, in contrast to 14 meV
(13%) for the conventional double-grid method. Also, the choice of
$\eta_{\mathrm{\large}}$ gives only a minor effect on the convergence rate,
and more importantly, the convergence is consistently faster than the
conventional double-grid method for a wide range of $\eta_{\mathrm{\large}}$
in [1 eV, 2 eV].

## References

  * Kohn (1996) W. Kohn, “Density functional and density matrix method scaling linearly with the number of atoms,” Phys. Rev. Lett. 76, 3168–3171 (1996).
  * Prodan and Kohn (2005) E. Prodan and W. Kohn, “Nearsightedness of electronic matter,” Proceedings of the National Academy of Sciences 102, 11635–11638 (2005).
  * Wannier (1937) Gregory H. Wannier, “The structure of electronic excitation levels in insulating crystals,” Phys. Rev. 52, 191–197 (1937).
  * Kohn (1959) W. Kohn, “Analytic properties of Bloch waves and Wannier functions,” Phys. Rev. 115, 809–821 (1959).
  * Marzari and Vanderbilt (1997) Nicola Marzari and David Vanderbilt, “Maximally-localized generalized Wannier functions for composite energy bands,” Physical Review B 56, 12847–12865 (1997).
  * Souza  _et al._ (2001) Ivo Souza, Nicola Marzari, and David Vanderbilt, “Maximally localized Wannier functions for entangled energy bands,” Physical Review B 65, 035109 (2001).
  * Marzari  _et al._ (2012) Nicola Marzari, Arash A. Mostofi, Jonathan R. Yates, Ivo Souza, and David Vanderbilt, “Maximally localized Wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012).
  * Giustino  _et al._ (2007) Feliciano Giustino, Marvin L. Cohen, and Steven G. Louie, “Electron-phonon interaction using Wannier functions,” Physical Review B 76, 165108 (2007).
  * Zhou and Bernardi (2016) Jin-Jian Zhou and Marco Bernardi, “ _Ab Initio_ electron mobility and polar phonon scattering in GaAs,” Physical Review B 94, 201201 (2016).
  * Poncé  _et al._ (2018) Samuel Poncé, Elena R. Margine, and Feliciano Giustino, “Towards predictive many-body calculations of phonon-limited carrier mobilities in semiconductors,” Phys. Rev. B 97, 121201 (2018).
  * Noffsinger  _et al._ (2012) Jesse Noffsinger, Emmanouil Kioupakis, Chris G. Van de Walle, Steven G. Louie, and Marvin L. Cohen, “Phonon-Assisted Optical Absorption in Silicon from First Principles,” Physical Review Letters 108, 167402 (2012).
  * Margine and Giustino (2013) E. R. Margine and F. Giustino, “Anisotropic migdal-eliashberg theory using Wannier functions,” Phys. Rev. B 87, 024505 (2013).
  * Jhalani  _et al._ (2017) Vatsal A. Jhalani, Jin-Jian Zhou, and Marco Bernardi, “Ultrafast Hot Carrier Dynamics in GaN and Its Impact on the Efficiency Droop,” Nano Letters 17, 5012–5019 (2017).
  * Wang  _et al._ (2006) Xinjie Wang, Jonathan R. Yates, Ivo Souza, and David Vanderbilt, “ _Ab Initio_ calculation of the anomalous Hall conductivity by Wannier interpolation,” Physical Review B 74, 195118 (2006).
  * Yates  _et al._ (2007) Jonathan R. Yates, Xinjie Wang, David Vanderbilt, and Ivo Souza, “Spectral and Fermi surface properties from Wannier interpolation,” Physical Review B 75, 195121 (2007).
  * Ibañez-Azpiroz  _et al._ (2018) Julen Ibañez-Azpiroz, Stepan S. Tsirkin, and Ivo Souza, “ _Ab Initio_ calculation of the shift photocurrent by Wannier interpolation,” Physical Review B 97, 245143 (2018).
  * Tsirkin  _et al._ (2018) Stepan S. Tsirkin, Pablo Aguado Puente, and Ivo Souza, “Gyrotropic effects in trigonal tellurium studied from first principles,” Physical Review B 97, 035158 (2018).
  * Qiao  _et al._ (2018) Junfeng Qiao, Jiaqi Zhou, Zhe Yuan, and Weisheng Zhao, “Calculation of intrinsic spin Hall conductivity by Wannier interpolation,” Physical Review B 98, 214402 (2018).
  * Ryoo  _et al._ (2019) Ji Hoon Ryoo, Cheol-Hwan Park, and Ivo Souza, “Computation of intrinsic spin Hall conductivities from first principles using maximally localized Wannier functions,” Physical Review B 99, 235113 (2019).
  * Sternheimer (1954) R. M. Sternheimer, “Electronic polarizabilities of ions from the Hartree-Fock wave functions,” Phys. Rev. 96, 951–968 (1954).
  * Baroni  _et al._ (2001) Stefano Baroni, Stefano de Gironcoli, Andrea Dal Corso, and Paolo Giannozzi, “Phonons and related crystal properties from density-functional perturbation theory,” Rev. Mod. Phys. 73, 515–562 (2001).
  * Ge and Lu (2015) Xiaochuan Ge and Deyu Lu, “Local representation of the electronic dielectric response function,” Physical Review B 92, 241107 (2015).
  * Giustino  _et al._ (2010a) Feliciano Giustino, Steven G. Louie, and Marvin L. Cohen, “Electron-phonon renormalization of the direct band gap of diamond,” Physical Review Letters 105, 265501 (2010a).
  * Gonze  _et al._ (2011) X. Gonze, P. Boulanger, and M. Côté, “Theoretical approaches to the temperature and zero-point motion effects on the electronic band structure,” Ann. Phys. 523, 168–178 (2011).
  * Sipe and Shkrebtii (2000) J. E. Sipe and A. I. Shkrebtii, “Second-order optical response in semiconductors,” Physical Review B 61, 5337–5352 (2000).
  * (26) Jae-Mo Lihm and Cheol-Hwan Park, unpublished .
  * Xu  _et al._ (2021) Haowei Xu, Hua Wang, Jian Zhou, and Ju Li, “Pure spin photocurrent in non-centrosymmetric crystals: Bulk spin photovoltaic effect,” Nature Communications 12, 4330 (2021).
  * Mu  _et al._ (2021) Xingchi Mu, Yiming Pan, and Jian Zhou, “Pure bulk orbital and spin photocurrent in two-dimensional ferroelectric materials,” npj Computational Materials 7, 61 (2021).
  * He and Vanderbilt (2001) Lixin He and David Vanderbilt, “Exponential decay properties of wannier functions and related quantities,” Phys. Rev. Lett. 86, 5341–5344 (2001).
  * Cornean  _et al._ (2019) Horia D Cornean, David Gontier, Antoine Levitt, and Domenico Monaco, “Localised Wannier functions in metallic systems,” in  _Annales Henri Poincaré_, Vol. 20 (Springer, 2019) pp. 1367–1391.
  * Macfarlane  _et al._ (1958) G. G. Macfarlane, T. P. McLean, J. E. Quarrington, and V. Roberts, “Fine structure in the absorption-edge spectrum of si,” Phys. Rev. 111, 1245–1254 (1958).
  * Allen and Heine (1976) P B Allen and V Heine, “Theory of the temperature dependence of electronic band structures,” Journal of Physics C: Solid State Physics 9, 2305–2312 (1976).
  * Allen and Cardona (1981) P. B. Allen and M. Cardona, “Theory of the temperature dependence of the direct gap of germanium,” Physical Review B 23, 1495–1505 (1981).
  * Allen and Cardona (1983) P. B. Allen and M. Cardona, “Temperature dependence of the direct gap of Si and Ge,” Phys. Rev. B 27, 4760–4769 (1983).
  * Giustino (2017) Feliciano Giustino, “Electron-phonon interactions from first principles,” Rev. Mod. Phys. 89, 015003 (2017).
  * Lihm and Park (2020) Jae-Mo Lihm and Cheol-Hwan Park, “Phonon-induced renormalization of electron wave functions,” Physical Review B 101, 121102 (2020).
  * Poncé  _et al._ (2014) S. Poncé, G. Antonius, Y. Gillet, P. Boulanger, J. Laflamme Janssen, A. Marini, M. Côté, and X. Gonze, “Temperature dependence of electronic eigenenergies in the adiabatic harmonic approximation,” Phys. Rev. B 90, 214304 (2014).
  * Nery  _et al._ (2018) Jean Paul Nery, Philip B. Allen, Gabriel Antonius, Lucia Reining, Anna Miglio, and Xavier Gonze, “Quasiparticles and phonon satellites in spectral functions of semiconductors and insulators: Cumulants applied to the full first-principles theory and the Fröhlich polaron,” Phys. Rev. B 97, 115145 (2018).
  * Miglio  _et al._ (2020) Anna Miglio, Véronique Brousseau-Couture, Emile Godbout, Gabriel Antonius, Yang-Hao Chan, Steven G. Louie, Michel Côté, Matteo Giantomassi, and Xavier Gonze, “Predominance of non-adiabatic effects in zero-point renormalization of the electronic band gap,” npj Computational Materials 6, 167 (2020).
  * Brown-Altvater  _et al._ (2020) Florian Brown-Altvater, Gabriel Antonius, Tonatiuh Rangel, Matteo Giantomassi, Claudia Draxl, Xavier Gonze, Steven G. Louie, and Jeffrey B. Neaton, “Band gap renormalization, carrier mobilities, and the electron-phonon self-energy in crystalline naphthalene,” Physical Review B 101, 165102 (2020).
  * Kioupakis  _et al._ (2010) Emmanouil Kioupakis, Patrick Rinke, André Schleife, Friedhelm Bechstedt, and Chris G. Van de Walle, “Free-carrier absorption in nitrides from first principles,” Physical Review B 81, 241201 (2010).
  * Patrick and Giustino (2014) Christopher E. Patrick and Feliciano Giustino, “Unified theory of electron-phonon renormalization and phonon-assisted optical absorption,” Journal of Physics: Condensed Matter 26, 365503 (2014).
  * Hybertsen and Louie (1986) Mark S. Hybertsen and Steven G. Louie, “Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies,” Physical Review B 34, 5390–5413 (1986).
  * Zacharias  _et al._ (2015) Marios Zacharias, Christopher E. Patrick, and Feliciano Giustino, “Stochastic approach to phonon-assisted optical absorption,” Physical Review Letters 115, 177401 (2015).
  * Zacharias and Giustino (2016) Marios Zacharias and Feliciano Giustino, “One-shot calculation of temperature-dependent optical spectra and phonon-induced band-gap renormalization,” Physical Review B 94, 075125 (2016).
  * Monserrat and Vanderbilt (2016) Bartomeu Monserrat and David Vanderbilt, “Temperature effects in the band structure of topological insulators,” Physical Review Letters 117, 226801 (2016).
  * Gorelov  _et al._ (2020) Vitaly Gorelov, David M. Ceperley, Markus Holzmann, and Carlo Pierleoni, “Electronic structure and optical properties of quantum crystals from first principles calculations in the Born–Oppenheimer approximation,” The Journal of Chemical Physics 153, 234117 (2020).
  * Zacharias  _et al._ (2020) Marios Zacharias, Matthias Scheffler, and Christian Carbogno, “Fully anharmonic nonperturbative theory of vibronically renormalized electronic band structures,” Phys. Rev. B 102, 045126 (2020).
  * Rangel  _et al._ (2020) Tonatiuh Rangel, Mauro Del Ben, Daniele Varsano, Gabriel Antonius, Fabien Bruneval, Felipe H. da Jornada, Michiel J. van Setten, Okan K. Orhan, David D. O’Regan, Andrew Canning, Andrea Ferretti, Andrea Marini, Gian-Marco Rignanese, Jack Deslippe, Steven G. Louie, and Jeffrey B. Neaton, “Reproducibility in $G_{0}W_{0}$ calculations for solids,” Computer Physics Communications 255, 107242 (2020).
  * Antonius  _et al._ (2014) G. Antonius, S. Poncé, P. Boulanger, M. Côté, and X. Gonze, “Many-Body Effects on the Zero-Point Renormalization of the Band Structure,” Physical Review Letters 112, 215501 (2014).
  * Bludau  _et al._ (1974) W. Bludau, A. Onton, and W. Heinke, “Temperature dependence of the band gap of silicon,” Journal of Applied Physics 45, 1846–1848 (1974).
  * Hamann and Vanderbilt (2009) D. R. Hamann and David Vanderbilt, “Maximally localized Wannier functions for $GW$ quasiparticles,” Phys. Rev. B 79, 045109 (2009).
  * Note (1) 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.
  * Huang  _et al._ (2016) H. H. Huang, Xiaofeng Fan, David J. Singh, Hong Chen, Q. Jiang, and W. T. Zheng, “Controlling phase transition for single-layer MTe2 (M = Mo and W): modulation of the potential barrier under strain,” Phys. Chem. Chem. Phys. 18, 4086–4094 (2016).
  * Baroni and Resta (1986) Stefano Baroni and Raffaele Resta, “Ab initio calculation of the macroscopic dielectric constant in silicon,” Phys. Rev. B 33, 7017–7021 (1986).
  * Guo  _et al._ (2008) G. Y. Guo, S. Murakami, T.-W. Chen, and N. Nagaosa, “Intrinsic Spin Hall Effect in Platinum: First-Principles Calculations,” Physical Review Letters 100, 096401 (2008).
  * Guo  _et al._ (2005) G. Y. Guo, Yugui Yao, and Qian Niu, “Ab initio calculation of the intrinsic spin hall effect in semiconductors,” Phys. Rev. Lett. 94, 226601 (2005).
  * Schlipf  _et al._ (2018) Martin Schlipf, Samuel Poncé, and Feliciano Giustino, “Carrier Lifetimes and Polaronic Mass Enhancement in the Hybrid Halide Perovskite CH3NH3PbI3 from Multiphonon Fröhlich Coupling,” Physical Review Letters 121, 086402 (2018).
  * D’Souza  _et al._ (2020) Ransell D’Souza, Jiang Cao, José D. Querales-Flores, Stephen Fahy, and Ivana Savić, “Electron-phonon scattering and thermoelectric transport in p -type PbTe from first principles,” Physical Review B 102, 115204 (2020).
  * Vitale  _et al._ (2020) Valerio Vitale, Giovanni Pizzi, Antimo Marrazzo, Jonathan R. Yates, Nicola Marzari, and Arash A. Mostofi, “Automated high-throughput Wannierisation,” npj Computational Materials 6, 66 (2020).
  * Bistoni  _et al._ (2019) O Bistoni, P Barone, E Cappelluti, L Benfatto, and F Mauri, “Giant effective charges and piezoelectricity in gapped graphene,” 2D Materials 6, 045015 (2019).
  * Binci  _et al._ (2021) Luca Binci, Paolo Barone, and Francesco Mauri, “First-principles theory of infrared vibrational spectroscopy of metals and semimetals: Application to graphite,” Physical Review B 103, 134304 (2021).
  * Dreyer  _et al._ (2021) Cyrus E Dreyer, Sinisa Coh, and Massimiliano Stengel, “Nonadiabatic Born effective charges in metals and the Drude weight,” arXiv:2103.04425 [cond-mat] , 16 (2021), arXiv:2103.04425 .
  * Eiguren and Ambrosch-Draxl (2008) Asier Eiguren and Claudia Ambrosch-Draxl, “Wannier interpolation scheme for phonon-induced potentials: Application to bulk ${\text{MgB}}_{2}$, W, and the $(1\times{}1)$ H-covered W(110) surface,” Phys. Rev. B 78, 045124 (2008).
  * Giustino  _et al._ (2010b) Feliciano Giustino, Marvin L. Cohen, and Steven G. Louie, “$GW$ method with the self-consistent Sternheimer equation,” Phys. Rev. B 81, 115105 (2010b).
  * Umari  _et al._ (2010) P. Umari, Geoffrey Stenuit, and Stefano Baroni, “$GW$ quasiparticle spectra from occupied states only,” Phys. Rev. B 81, 115104 (2010).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R.M. Wentzcovitch, “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials,” J. Phys.: Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, D Ceresoli, M Cococcioni, N Colonna, I Carnimeo, A Dal Corso, S de Gironcoli, P Delugas, R A DiStasio, A Ferretti, A Floris, G Fratesi, G Fugallo, R Gebauer, U Gerstmann, F Giustino, T Gorni, J Jia, M Kawamura, H-Y Ko, A Kokalj, E Küçükbenli, M Lazzeri, M Marsili, N Marzari, F Mauri, N L Nguyen, H-V Nguyen, A Otero-de-la-Roza, L Paulatto, S Poncé, D Rocca, R Sabatini, B Santra, M Schlipf, A P Seitsonen, A Smogunov, I Timrov, T Thonhauser, P Umari, N Vast, X Wu, and S Baroni, “Advanced capabilities for materials modelling with Quantum ESPRESSO,” Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Hamann (2013) D. R. Hamann, “Optimized norm-conserving Vanderbilt pseudopotentials,” Physical Review B 88, 085117 (2013).
  * van Setten  _et al._ (2018) M.J. van Setten, M. Giantomassi, E. Bousquet, M.J. Verstraete, D.R. Hamann, X. Gonze, and G.-M. Rignanese, “The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table,” Computer Physics Communications 226, 39–54 (2018).
  * Perdew and Wang (1992) John P. Perdew and Yue Wang, “Accurate and simple analytic representation of the electron-gas correlation energy,” Physical Review B 45, 13244–13249 (1992).
  * Reeber and Wang (1996) Robert R. Reeber and Kai Wang, “Thermal expansion and lattice parameters of group IV semiconductors,” Materials Chemistry and Physics 46, 259 – 264 (1996).
  * Lautenschlager  _et al._ (1985) P. Lautenschlager, P. B. Allen, and M. Cardona, “Temperature dependence of band gaps in si and ge,” Phys. Rev. B 31, 2163–2171 (1985).
  * Perdew  _et al._ (1996) John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
  * Pizzi  _et al._ (2020) Giovanni Pizzi, Valerio Vitale, Ryotaro Arita, Stefan Blügel, Frank Freimuth, Guillaume Géranton, Marco Gibertini, Dominik Gresch, Charles Johnson, Takashi Koretsune, Julen Ibañez-Azpiroz, Hyungjun Lee, Jae-Mo Lihm, Daniel Marchand, Antimo Marrazzo, Yuriy Mokrousov, Jamal I Mustafa, Yoshiro Nohara, Yusuke Nomura, Lorenzo Paulatto, Samuel Poncé, Thomas Ponweiser, Junfeng Qiao, Florian Thöle, Stepan S Tsirkin, Małgorzata Wierzbowska, Nicola Marzari, David Vanderbilt, Ivo Souza, Arash A Mostofi, and Jonathan R Yates, “Wannier90 as a community code: New features and applications,” Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Deslippe  _et al._ (2012) Jack Deslippe, Georgy Samsonidze, David A. Strubbe, Manish Jain, Marvin L. Cohen, and Steven G. Louie, “BerkeleyGW: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures,” Computer Physics Communications 183, 1269–1289 (2012).
  * Deslippe  _et al._ (2013) Jack Deslippe, Georgy Samsonidze, Manish Jain, Marvin L. Cohen, and Steven G. Louie, “Coulomb-hole summations and energies for $GW$ calculations with limited number of empty orbitals: A modified static remainder approach,” Phys. Rev. B 87, 165124 (2013).
  * Poncé  _et al._ (2016) S. Poncé, E.R. Margine, C. Verdi, and F. Giustino, “EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions,” Computer Physics Communications 209, 116–133 (2016).
  * Rohlfing and Louie (2000) Michael Rohlfing and Steven G. Louie, “Electron-hole excitations and optical spectra from first principles,” Phys. Rev. B 62, 4927–4944 (2000).
  * Green and Keevers (1995) Martin A Green and Mark J Keevers, “Optical properties of intrinsic silicon at 300 K,” Progress in Photovoltaics: Research and Applications 3, 189–192 (1995).
  * Nastos and Sipe (2006) F. Nastos and J. E. Sipe, “Optical rectification and shift currents in GaAs and GaP response: Below and above the band gap,” Physical Review B 74, 035201 (2006).
  * Tsirkin (2021) Stepan S. Tsirkin, “High performance Wannier interpolation of Berry curvature and related quantities with WannierBerri code,” npj Computational Materials 7, 33 (2021).
  * Lee  _et al._ (2015) Chia-Hui Lee, Eduardo Cruz Silva, Lazaro Calderin, Minh An T. Nguyen, Matthew J. Hollander, Brian Bersch, Thomas E. Mallouk, and Joshua A. Robinson, “Tungsten ditelluride: a layered semimetal,” Scientific Reports 5, 10013 (2015).
  * Momma and Izumi (2011) Koichi Momma and Fujio Izumi, “VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data,” Journal of Applied Crystallography 44, 1272–1276 (2011).
  * Poncé  _et al._ (2015) S. Poncé, Y. Gillet, J. Laflamme Janssen, A. Marini, M. Verstraete, and X. Gonze, “Temperature dependence of the electronic structure of semiconductors and insulators,” The Journal of Chemical Physics 143, 102813 (2015).
  * Querales-Flores  _et al._ (2019) José D. Querales-Flores, Jiang Cao, Stephen Fahy, and Ivana Savić, “Temperature effects on the electronic band structure of PbTe from first principles,” Physical Review Materials 3, 055405 (2019).
  * Poncé  _et al._ (2021) Samuel Poncé, Francesco Macheda, Elena Roxana Margine, Nicola Marzari, Nicola Bonini, and Feliciano Giustino, “First-principles predictions of hall and drift mobilities in semiconductors,” Phys. Rev. Research 3, 043022 (2021).
  * Adolph  _et al._ (1996) B. Adolph, V. I. Gavrilenko, K. Tenelsen, F. Bechstedt, and R. Del Sole, “Nonlocality and many-body effects in the optical properties of semiconductors,” Phys. Rev. B 53, 9797–9808 (1996).
  * Lopez  _et al._ (2012) M. G. Lopez, David Vanderbilt, T. Thonhauser, and Ivo Souza, “Wannier-based calculation of the orbital magnetization in crystals,” Physical Review B 85, 014435 (2012).