Efficient Ab Initio Calculations of Electron-Defect Scattering and Defect-Limited Carrier Mobility

TL;DR

This paper introduces an efficient first-principles method to calculate electron-defect interactions, relaxation times, and defect-limited mobility at low temperatures, revealing energy and defect type dependencies overlooked by traditional models.

Contribution

The authors develop a computationally efficient ab initio approach for electron-defect scattering, enabling accurate low-temperature mobility predictions for materials with defects.

Findings

Electron-defect relaxation times vary significantly with energy and defect type.

Defect-limited mobility shows temperature dependence contrary to conventional assumptions.

The new method surpasses traditional heuristic models in accuracy for defect-related charge transport.

Abstract

Electron-defect ($e$-d) interactions govern charge carrier dynamics at low temperature, where they limit the carrier mobility and give rise to phenomena of broad relevance in condensed matter physics. Ab initio calculations of $e$-d interactions are still in their infancy, mainly because they require large supercells and computationally expensive workflows. Here we develop an efficient ab initio approach for computing elastic $e$-d interactions, their associated $e$-d relaxation times (RTs), and the low-temperature defect-limited carrier mobility. The method is applied to silicon with simple neutral defects, such as vacancies and interstitials. Contrary to conventional wisdom, the computed $e$-d RTs depend strongly on carrier energy and defect type, and the defect-limited mobility is temperature dependent. These results highlight the shortcomings of widely employed heuristic models of $e$-d interactions in materials. Our method opens new avenues for studying $e$-d scattering and low-temperature charge transport from first principles.

Equations

$${\Gamma }_{n\mathbf{k}}=\frac{2\pi }{\hslash }\frac{n_{\text{at}}{C}_{\mathrm{d}}}{N_{{\mathbf{k}}^{\prime }}}\sum _{n^{\prime }{\mathbf{k}}^{\prime }}{\left|M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}\right|}^2\delta \left({\epsilon }_{n^{\prime }{\mathbf{k}}^{\prime }}-{\epsilon }_{n\mathbf{k}}\right),$$

$$M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}=⟨n^{\prime }{\mathbf{k}}^{\prime }⟩\Delta V_{\mathrm{e}-\mathrm{d}}n\mathbf{k}.$$

$$V_{\mathrm{KS}}=V_{\mathrm{L}}(\mathbf{r})+{\hat{V}}_{\mathrm{NL}}.$$

$$V_{\mathrm{L}}(\mathbf{r})=V_{\mathrm{H}}(\mathbf{r})+V_{\mathrm{XC}}(\mathbf{r})+V_{\mathrm{pp}}(\mathbf{r}).$$

$${\hat{V}}_{\text{NL}}=\sum _{s=1}\sum _{ij}{D}_{ij}^{(s)}|{\beta }_i^{(s)}⟩⟨{\beta }_j^{(s)}|,$$

$$M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}=M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\text{L}}+M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\text{NL}},$$

$$\begin{array}{rl}{M}_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\text{L}}& =⟨n^{\prime }{\mathbf{k}}^{\prime }⟩\Delta V_{\mathrm{L}}(\mathbf{r})n\mathbf{k}\\ M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\text{NL}}& =⟨n^{\prime }{\mathbf{k}}^{\prime }⟩\Delta {\hat{V}}_{\mathrm{NL}}n\mathbf{k}.\end{array}$$

$$M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\text{L}}=\sum _{\mathbf{G}}\Delta {\tilde{V}}_{\mathrm{L}}({\mathbf{k}}^{\prime }-\mathbf{k}+\mathbf{G})⟨u_{n^{\prime }{\mathbf{k}}^{\prime }}⟩e^{i\mathbf{G}\cdot \mathbf{r}}{u_{n\mathbf{k}}}_{\mathrm{uc}},$$

$$\Delta {\tilde{V}}_{\mathrm{L}}(\mathbf{q})=\frac{1}{{\Omega }_{\mathrm{uc}}}{\int }_{{\Omega }_{\text{sup}}}d\mathbf{r}\Delta V_{\mathrm{L}}\left(\mathbf{r}\right)e^{-i\mathbf{q}\cdot \mathbf{r}},$$

$$\begin{array}{rl}{M}_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\mathrm{NL}}& =⟨n^{\prime }{\mathbf{k}}^{\prime }⟩{\hat{V}}_{\mathrm{NL}}^{(\mathrm{d})}-{\hat{V}}_{\mathrm{NL}}^{(\mathrm{p})}n\mathbf{k}\\ & =⟨n^{\prime }{\mathbf{k}}^{\prime }⟩{\hat{V}}_{\mathrm{NL}}^{(\mathrm{d})}n\mathbf{k}-⟨n^{\prime }{\mathbf{k}}^{\prime }⟩{\hat{V}}_{\mathrm{NL}}^{(\mathrm{p})}n\mathbf{k},\end{array}$$

$$\begin{array}{rl}⟨n^{\prime }{\mathbf{k}}^{\prime }⟩{\hat{V}}_{\mathrm{NL}}^{(\alpha )}n\mathbf{k}=\sum _{s_{\alpha }=1}\sum _{ij}& D_{ij}^{(s_{\alpha })}⟨n^{\prime }{\mathbf{k}}^{\prime }⟩{\beta }_i^{(s_{\alpha })}⟨{\beta }_j^{(s_{\alpha })}⟩n\mathbf{k}.\end{array}$$

$$⟨{\beta }_j^{(s)}⟩n\mathbf{k}=\frac{1}{\sqrt{{\Omega }_{\mathrm{uc}}}}\sum _{\mathbf{G}}{B}_{j\mathbf{k}}^{(s)\ast }(\mathbf{G})⟨e^{i\mathbf{G}\cdot \mathbf{r}}⟩{u_{n\mathbf{k}}}_{\mathrm{uc}},$$

$$B_{j\mathbf{k}}^{(s)}(\mathbf{G})=\frac{1}{\sqrt{{\Omega }_{\mathrm{uc}}}}{\int }_{{\Omega }_{\sup }}d\mathbf{r}\left({\beta }_j^{(s)}(\mathbf{r})e^{-i\mathbf{k}\cdot \mathbf{r}}\right)e^{-i\mathbf{G}\cdot \mathbf{r}}.$$

$$M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\mathrm{L}}=⟨n^{\prime }{\mathbf{k}}^{\prime }⟩\Delta V_{\mathrm{L}}(\mathbf{r}){n\mathbf{k}}_{\sup }$$

$${\mu }_{\alpha \beta }(T)=\frac{e}{n_c}{\int }_{-\infty }^{+\infty }dE\left[-\partial f(T,E)/\partial E\right]\times {\Sigma }_{\alpha \beta }\left(E\right),$$

$${\Sigma }_{\alpha \beta }\left(E\right)=\frac{2}{{\Omega }_{\mathrm{uc}}}\sum _{n\mathbf{k}}{\tau }_{n\mathbf{k}}{\mathbf{v}}_{n\mathbf{k}}^{\alpha }{\mathbf{v}}_{n\mathbf{k}}^{\beta }\delta \left(E-{\epsilon }_{n\mathbf{k}}\right),$$

$${\Gamma }_{n\mathbf{k}}=\frac{2\pi }{\hslash }\sum _{n^{\prime }{\mathbf{k}}^{\prime }}{\left|\text{matrixelement}{\psi }_{n^{\prime }{\mathbf{k}}^{\prime }}\Delta \hat{H}{\psi }_{n\mathbf{k}}\right|}^2\delta \left({\epsilon }_{n^{\prime }{\mathbf{k}}^{\prime }}-{\epsilon }_{n\mathbf{k}}\right),$$

$$⟨\mathbf{r}⟩{\psi }_{n\mathbf{k}}=\frac{1}{\sqrt{N_{\mathbf{k}}}}⟨\mathbf{r}⟩n\mathbf{k}=\frac{1}{\sqrt{N_{\mathbf{k}}}}{u}_{n\mathbf{k}}\left(\mathbf{r}\right)e^{i\mathbf{k}\cdot \mathbf{r}},$$

$${\int }_{{\Omega }_{\mathrm{uc}}}{d}^{3}r{u}_{n\mathbf{k}}^{\ast }\left(\mathbf{r}\right)u_{n\mathbf{k}}\left(\mathbf{r}\right)=1.$$

$$\Delta \hat{H}={\hat{H}}^{(\mathrm{d})}-{\hat{H}}^{(\mathrm{p})}=\sum _{i=1}^{N_{\mathrm{d}}}\Delta V_{\mathrm{e}-\mathrm{d}}(\mathbf{r}-{\mathbf{r}}_i),$$

$${\Gamma }_{n\mathbf{k}}=\frac{2\pi }{\hslash }\frac{n_{\mathrm{at}}{C}_{\mathrm{d}}}{N_{{\mathbf{k}}^{\prime }}}\sum _{n^{\prime }{\mathbf{k}}^{\prime }}{\left|M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}\right|}^2\delta \left({\epsilon }_{n^{\prime }{\mathbf{k}}^{\prime }}-{\epsilon }_{n\mathbf{k}}\right),$$

$$M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}=\text{matrixelement}{n}^{\prime }{\mathbf{k}}^{\prime }\Delta V_{\mathrm{e}-\mathrm{d}}n\mathbf{k}.$$

$$\begin{array}{rl}& M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\mathrm{L}}=\text{matrixelement}{n}^{\prime }{\mathbf{k}}^{\prime }\Delta V_{\mathrm{L}}\left(\mathbf{r}\right)n\mathbf{k}\\ & ={\int }_{{\Omega }_{\text{BvK}}}{d}^{3}r{u}_{n^{\prime }{\mathbf{k}}^{\prime }}^{\ast }\left(\mathbf{r}\right)e^{-i{\mathbf{k}}^{\prime }\cdot \mathbf{r}}\Delta V_{\mathrm{L}}\left(\mathbf{r}\right)u_{n\mathbf{k}}\left(\mathbf{r}\right)e^{i\mathbf{k}\cdot \mathbf{r}}.\end{array}$$

$$\Delta V_{\mathrm{L}}\left(\mathbf{r}\right)=\frac{1}{N_{\mathbf{k}}}\sum _{\mathbf{q}}\Delta {\tilde{V}}_{\mathrm{L}}\left(\mathbf{q}\right)e^{i\mathbf{q}\cdot \mathbf{r}}$$

$$\begin{array}{rl}\Delta {\tilde{V}}_{\mathrm{L}}\left(\mathbf{q}\right)=& \frac{1}{{\Omega }_{\mathrm{uc}}}{\int }_{{\Omega }_{\mathrm{BvK}}}{d}^{3}r\Delta V_{\mathrm{L}}\left(\mathbf{r}\right)e^{-i\mathbf{q}\cdot \mathbf{r}}\\ \approx & \frac{1}{{\Omega }_{\mathrm{uc}}}{\int }_{{\Omega }_{\sup }}{d}^{3}r\Delta V_{\mathrm{L}}\left(\mathbf{r}\right)e^{-i\mathbf{q}\cdot \mathbf{r}},\end{array}$$

$$M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\mathrm{L}}=\sum _{\mathbf{G}}\Delta {\tilde{V}}_{\mathrm{L}}({\mathbf{k}}^{\prime }-\mathbf{k}+\mathbf{G})\text{matrixelement}{u}_{n^{\prime }{\mathbf{k}}^{\prime }}{e}^{i\mathbf{G}\cdot \mathbf{r}}{u_{n\mathbf{k}}}_{\mathrm{uc}},$$

$$\text{matrixelement}{u}_{n^{\prime }{\mathbf{k}}^{\prime }}{e}^{i\mathbf{G}\cdot \mathbf{r}}{u_{n\mathbf{k}}}_{\mathrm{uc}}={\int }_{{\Omega }_{\mathrm{uc}}}{d}^{3}r{u}_{n^{\prime }{\mathbf{k}}^{\prime }}^{\ast }\left(\mathbf{r}\right)e^{i\mathbf{G}\cdot \mathbf{r}}{u}_{n\mathbf{k}}\left(\mathbf{r}\right),$$

$$u_{n\mathbf{k}}\left(\mathbf{r}\right)=\frac{1}{\sqrt{{\Omega }_{\mathrm{uc}}}}\sum _{\mathbf{G}}{C}_{n\mathbf{k}}\left(\mathbf{G}\right)e^{i\mathbf{G}\cdot \mathbf{r}},$$

$$\begin{array}{rl}& M_{n^{\prime }{\mathbf{k}}^{\prime },n\mathbf{k}}^{\mathrm{L}}=\sum _{\mathbf{G}}\Delta {\tilde{V}}_{\mathrm{L}}({\mathbf{k}}^{\prime }-\mathbf{k}+\mathbf{G})\\ & \times \left[\sum _{{\mathbf{G}}^{{}^{\prime \prime }}}\sum _{{\mathbf{G}}^{{}^{\prime }}}{C}_{n^{\prime }{\mathbf{k}}^{\prime }}^{\ast }({\mathbf{G}}^{{}^{\prime \prime }})C_{n\mathbf{k}}({\mathbf{G}}^{{}^{\prime }}){\delta }_{{\mathbf{G}}^{{}^{\prime \prime }},{\mathbf{G}}^{{}^{\prime }}+\mathbf{G}}\right].\end{array}$$

$$\begin{array}{rl}& ⟨{\beta }_j^{(s)}⟩n\mathbf{k}={\int }_{{\Omega }_{\mathrm{BvK}}}{d}^{3}r{\beta }_j^{\ast }\left(\mathbf{r}-{\tau }_s\right)e^{i\mathbf{k}\cdot \mathbf{r}}{u}_{n\mathbf{k}}\left(\mathbf{r}\right)\\ & ={\int }_{{\Omega }_{\mathrm{BvK}}}{d}^{3}r{\left[\frac{{\beta }_j\left(\mathbf{r}-{\tau }_s\right)}{\sqrt{{\Omega }_{\mathrm{uc}}}}{e}^{-i\mathbf{k}\cdot \mathbf{r}}\right]}^{\ast }\sqrt{{\Omega }_{\mathrm{uc}}}{u}_{n\mathbf{k}}\left(\mathbf{r}\right).\end{array}$$

Full text

I Introduction

Defects in materials can scatter or capture charge carriers. They control carrier dynamics at low temperature, where phonons are frozen out, or even at room temperature in the presence of strong disorder. Impurity defects are key to controlling the carrier mobility in semiconductors Radisavljevic and Kis (2013); Wang et al. (2012); Li et al. (2016); Herz (2017); Mao et al. (2017); Sun et al. (2015), and defects are widely employed to engineer advanced materials and devices Xu et al. (2017); Koenraad and Flatté (2011); Gunlycke and White (2011); Yazyev and Louie (2010). Weak localization Bergmann (1984); Altshuler et al. (1980); Datta (1997) and conductance fluctuations Lee and Stone (1985); Akkermans and Montambaux (2007) are but two examples of subtle quantum effects that can be induced by defects. For all these reasons, understanding how defects interact with charge carriers is important in both applied and fundamental condensed matter physics.

Calculations of electron-defect ($e$-d) interactions have so far relied almost exclusively on empirical models Chattopadhyay and Queisser (1981), in which the atomic details of the defects are suppressed and simple approximations are employed for the electronic screening and band structure. For example, the perturbation induced by charged defects is typically modeled with a simple Yukawa potential, and models of $e$-d scattering due to neutral (non-charged) defects are scarce and mainly limited to spherical-well-potential approaches Erginsoy (1950); Sclar (1956); Blakemore (1980). Yet, the perturbation induced by defects in real materials can be anisotropic and complex, and its effects on electronic screening intricate. How the band structure and defect states impact $e$-d scattering is difficult to predict.

Due to their severe approximations, empirical models are not reliable to compute $e$-d interactions, and their predictions can be qualitatively and quantitatively incorrect. First-principles approaches, such as those based on density functional theory (DFT) and related methods Martin (2004), provide quantitative accuracy in many areas of materials physics, and are therefore desirable to treat $e$-d interactions. Such ab initio $e$-d calculations can take advantage of existing tools developed for electron-phonon ($e$-ph) interactions, for which accurate calculations of relaxation times (RTs) Bernardi et al. (2014, 2015a, 2015b); Tanimura et al. (2016), matrix elements and their interpolation Agapito and Bernardi (2018), and phonon-limited carrier mobility Mustafa et al. (2016); Park et al. (2014); Zhou and Bernardi (2016); Lee et al. (2018); Zhou et al. (2018) have been recently developed. These and other workflows developed for $e$-ph calculations can be generalized to treat $e$-d interactions.

However, there are open challenges specific to $e$-d calculations that currently prevent their broad applicability. First-principles $e$-d calculations need large supercells to obtain the electron wavefunctions and defect perturbation potentials; they additionally involve computationally costly $e$-d matrix elements, and require systematic convergence of the RTs with respect to supercell size and Brillouin zone (BZ) grids. Early work on first-principles $e$-d interactions Evans et al. (2005); Restrepo et al. (2009); Lordi et al. (2010) proposed an approach, called here the all-supercell method, in which the cost of computing $e$-d matrix element is prohibitively large due to the use of electron wavefunctions from large supercells. This work aims to address these open challenges and make ab initio $e$-d calculations more affordable and broadly applicable.

Here we formulate an efficient scheme to compute and converge the $e$-d matrix elements and the associated RTs and defect-limited carrier mobility. Our approach does not require the wavefunctions of large supercells, which dramatically speeds up the calculations. We apply our method to study $e$-d interactions in silicon in two cases, neutral vacancy and interstitial defects, for which we compute and converge the $e$-d RTs as a function of energy and the carrier mobility as a function of temperature below 150 K. The results show that, contrary to conventional wisdom, the RTs depend strongly on carrier energy and defect type, and the defect-limited mobility for neutral defects depends on temperature. Our results provide new microscopic insight into $e$-d scattering, and our approach, together with its future extensions, can uncover new defect physics in materials and devices for electronics, energy and quantum technologies.

The paper is organized as follows: Section II describes the new formulas and workflow for computing $e$-d matrix elements and RTs, and compares our approach with the all-supercell method. In Section III we compute and converge $e$-d RTs for electrons and holes in silicon and the associated defect-limited mobility at low temperature. Section IV discusses technical points and outlines future research directions.

II Methodology

In the following, we work under the assumption that the defects are neutral (non-charged) and that the $e$-d scattering events are independent, uncorrelated and elastic. The $e$-d scattering rate $\Gamma_{n\bm{k}}$ (and its inverse, the RT, $\tau_{n\bm{k}}=\Gamma_{n\bm{k}}^{-1}$) for a Bloch state $\ket{n\bm{k}}$, where $n$ is the band index and $\bm{k}$ the crystal momentum, is computed using lowest-order perturbation theory (see Appendix A):

[TABLE]

where $n_{\text{at}}$ is the number of atoms in a primitive cell, $C_{\rm d}$ the (dimensionless) defect atomic concentration, $N_{\bm{k}^{\prime}}$ the number of BZ $\bm{k}^{\prime}$-points used in the sum, and $\varepsilon_{n\bm{k}}$ the unperturbed energy of the state $\ket{n\bm{k}}$. The $e$-d matrix elements $M_{n^{\prime}\bm{k}^{\prime},n\bm{k}}$ are the central quantities computed in this work; they encode the probability amplitude for scattering from the unperturbed state $\ket{n\bm{k}}$ to $\ket{n^{\prime}\bm{k}^{\prime}}$ due to the perturbation potential $\Delta V_{\rm e-d}$ from a defect:

[TABLE]

Within DFT Kohn and Sham (1965), the $e$-d perturbation can be computed as the difference between the Kohn-Sham (KS) potentials $V_{\rm KS}$ of a defect-containing supercell and a pristine supercell with no defect, namely, $\Delta V_{\rm e-d}=V_{\rm KS}^{\rm(d)}-V_{\rm KS}^{\rm(p)}$. Here and below, we use superscripts (d) and (p) to denote the defect-containing and pristine systems, respectively. When using, as we do here, norm-conserving pseudopotentials in the Kleinman-Bylander (KB) form Kleinman and Bylander (1982), the Kohn-Sham potential can be written as a sum of local and nonlocal parts Martin (2004):

[TABLE]

The local potential $V_{\rm L}(\bm{r})$ comprises the Hartree and exchange-correlation potentials plus the local part of the pseudopotentials,

[TABLE]

The nonlocal potential $\hat{V}_{\text{NL}}$, which is due to the pseudopotentials, is as a sum over all atoms in the supercell of KB projectors $\ket{\beta_{i}^{(s)}}$, each localized in the core region of atom $s$:

[TABLE]

where $i$ and $j$ are orbital angular momentum quantum numbers, and $D_{ij}^{(s)}$ are KB coefficients Martin (2004).

Accordingly, we separate the matrix elements into a local and a nonlocal part,

[TABLE]

each due to the respective defect perturbation, $\Delta V_{\rm L}(\bm{r})=V_{L}^{\rm(d)}(\bm{r})-V_{L}^{\rm(p)}(\bm{r})$ for the local and $\Delta\hat{V}_{\rm NL}=\hat{V}_{\rm NL}^{\rm(d)}-\hat{V}_{\rm NL}^{\rm(p)}$ for the nonlocal part:

[TABLE]

II.1 Electron-defect matrix element computation

We develop a new approach to reduce the computational cost of the local and nonlocal matrix elements. We Fourier transform the local perturbation potential $\Delta V_{\text{L}}(\bm{r})$ and compute the local matrix elements as (see Appendix B)

[TABLE]

where $u_{n\bm{k}}(\bm{r})$ is the periodic part of the Bloch wavefunction (normalized in a primitive cell with volume $\Omega_{\text{uc}}$), $\bm{G}$ are reciprocal lattice vectors of the primitive cell, and $\Delta\widetilde{V}_{\rm L}$ are Fourier coefficients of the local perturbation potential (computed in a supercell with volume $\Omega_{\text{sup}}$):

[TABLE]

where $\bm{q}=\bm{k}^{\prime}-\bm{k}+\bm{G}$ is the transferred momentum in the $\ket{n\bm{k}}\rightarrow\ket{n^{\prime}\bm{k}^{\prime}}$ scattering process.

Through Eq. (8), we effectively separate the local matrix element calculation into two independent parts, the plane-wave matrix elements $\braket{u_{n^{\prime}\bm{k}^{\prime}}}{\,e^{i\bm{G}\cdot\bm{r}}\,}{{u_{n\bm{k}}}}_{\rm uc}$, which are easily computed by integrating over the primitive cell, and the Fourier coefficients $\Delta\widetilde{V}_{\rm L}(\bm{q})$, which are computed in the supercell. In addition, since the local perturbation potential is smooth and decays in real space over a few angstroms, its Fourier coefficients $\Delta\widetilde{V}_{\rm L}(\bm{q})$ decay rapidly in reciprocal space, and the summation over $\bm{G}$ in Eq. (8) can be truncated to a small cutoff $-$ typically, just a few reciprocal lattice vectors. For the same reason, the Fourier coefficients $\Delta\widetilde{V}_{\rm L}(\bm{q})$ can be efficiently interpolated (in our case, with a B-spline) at arbitrary $\bm{q}$ starting from their calculation at a few thousand $\bm{q}$-points in a cubic box centered at $\bm{q}=0$. A great advantage of this formulation is that one can compute the local matrix elements using only the wavefunctions of the primitive cell (as opposed to those of the supercell), and effectively interpolate the perturbation potential to fine BZ grids.

The nonlocal matrix elements $M_{n^{\prime}\bm{k}^{\prime},n\bm{k}}^{\text{NL}}$ are computed as the difference between the nonlocal potentials of a defect-containing and a pristine supercell:

[TABLE]

where the matrix elements of $\hat{V}_{\rm NL}$ for each of the two supercells (labelled by $\alpha=$ d, p), using Eq. (5), read:

[TABLE]

Similar to the local matrix elements, the computation in Eq. (11) is split into a primitive cell and a supercell calculation, by expressing the scalar products $\braket{\beta_{j}^{(s)}}{n\bm{k}}$ as (see Appendix C)

[TABLE]

where $B_{j\bm{k}}^{(s)}\left(\bm{G}\right)$ is the Fourier coefficient of the KB projector $\beta_{j}^{(s)}$ (multiplied by the phase factor $e^{-i\bm{k}\cdot\bm{r}}$) at the primitive-cell reciprocal lattice vector $\bm{G}$:

[TABLE]

Note that the nonlocal matrix elements, which are computed as the difference in Eq. (10), are nonzero both because the atomic positions change upon relaxing the structure in the defect-containing supercell and because the number and type of atoms, in general, differ in the two supercells, as is the case when considering a vacancy or impurity.

Figure 1 shows our workflow for computing properties related to $e$-d interactions from first principles. To obtain the relevant DFT inputs, the KS equations are solved in a primitive cell and separately in pristine and defect-containing supercells. The local and nonlocal matrix elements are then computed by splitting the calculation into a primitive cell and a supercell part, an approach that dramatically reduces the computational cost. Computing local matrix elements is the most expensive step, while the nonlocal matrix elements only involve relatively inexpensive reciprocal space sums. Importantly, only the KS wavefunctions, band structure and $\bm{k}$-points of the primitive cell are used, while the supercells are employed only to obtain the perturbation potential due to the defect. Once computed, the $e$-d matrix elements are employed to calculate the $e$-d RTs and the defect-limited mobility, among other quantities of interest. This approach allows us to systematically converge the RTs and other properties related to $e$-d interactions with respect to supercell size and BZ grids.

II.2 Comparison with the all-supercell method

In the all-supercell method Restrepo et al. (2009), one uses the pristine and defect-containing supercells to provide all the necessary quantities for computing the $e$-d matrix elements, including the wavefunctions, band structure, perturbation potentials, and BZ grids. However, using supercell wavefunctions makes it challenging to compute and converge the $e$-d matrix elements and RTs, and ultimately to carry out accurate $e$-d calculations, since unconverged $e$-d RTs and transport properties can differ widely from the converged results.

All $e$-d calculations need to be converged with respect to supercell size; as we discuss below, converging the RTs for a neutral defect typically requires very large supercells with hundreds of atoms. In our approach, this convergence does not constitute a challenge since the same (primitive cell) wavefunctions are employed, regardless of supercell size. Conversely, in the all-supercell method, one uses wavefunctions from the pristine supercell, and the computational cost to compute and store the wavefunctions and obtain the matrix elements increases dramatically with supercell size, making accurate convergence tests too computationally demanding.

Let us analyze the cost of the most computationally intensive step, namely obtaining the local $e$-d matrix elements, $M^{\rm L}_{n^{\prime}\bm{k}^{\prime},n\bm{k}}$. Using a uniform BZ grid with $N_{\bm{k}}$ points, one obtains $\mathcal{O}(N_{\bm{k}}^{2})$ matrix elements, each for a distinct $\ket{n\bm{k}}\rightarrow\ket{n^{\prime}\bm{k}^{\prime}}$ $e$-d scattering process. In a typical calculation, a uniform grid with at least $N_{\bm{k}}\approx 10^{6}$ points is needed to converge the RTs in the entire BZ. In a mobility calculation, one typically selects a small energy window of $\sim$100 meV near the band edges (in a semiconductor, or near the Fermi energy in a metal), which reduces the required number of $\bm{k}$-points to $N_{\bm{k}}\approx 10^{4}$.

In the all-supercell method, the local matrix elements are computed as:

[TABLE]

where the subscript (sup) denotes that both the local perturbation potential $\Delta V_{\rm L}(\bm{r})$ and the wavefunctions are obtained from a DFT calculation on a supercell. Since the cost of the DFT calculations scales with system size as roughly $N^{3}_{\rm at}$, where $N_{\rm at}$ is the number of atoms in the supercell, computing $N_{\bm{k}}$ supercell wavefunctions, from which the local matrix elements are computed on the uniform grid, costs $N_{\bm{k}}N^{3}_{\rm at}$ in the all-supercell method.

By contrast, in our method only the primitive cell wavefunctions are used, and thus the computational cost of the matrix elements does not depend on $N_{\rm at}$ through the wavefunctions. To obtain the local matrix elements on the uniform grid with our method [see Eq. (8)], the only supercell data one needs are the Fourier coefficients $\Delta\widetilde{V}_{\rm L}(\bm{q})$ of the local perturbation potential. Obtaining these coefficients at a few thousand $\bm{q}$ points $-$ from which an interpolation table can be constructed $-$ has a cost that scales as $N^{3}_{\rm at}$, but this step is required only once for a given supercell size. Therefore, computing the local matrix elements on a uniform BZ grid with $N_{\bm{k}}$ points has a cost of order $N^{3}_{\rm at}$ in our method, versus a cost of $N_{\bm{k}}N^{3}_{\rm at}$ in the all-supercell method. For the typical mobility calculation mentioned above, this represents a speed-up by a factor of $N_{\bm{k}}\approx 10,000$ over the all-supercell method. Note that our mobility calculations are already expensive (tens of thousands of core-hours), so approaches that are thousands of times more expensive are clearly impractical. An additional benefit is that in our approach the large supercell wavefunctions are never stored or loaded into memory, so the speed up is significant even for computing a single $e$-d RT.

Finally, one would like to map the $e$-d scattering processes onto the band structure of the primitive cell, as is done for $e$-ph scattering processes. This is possible in our approach due to our use of primitive cell band structures and $\bm{k}$-point grids, but impractical in the all-supercell method, where one uses the supercell band structures and $\bm{k}$-point grids, which depend on the choice of a supercell and differ from those of the primitive cell due to nontrivial BZ folding effects. Due to its computational efficiency and convenience, we thus believe that our approach solves key technical challenges that have so far prevented efficient and accurate ab initio calculations of $e$-d interactions.

Figure 2 validates our approach, by comparing the local matrix elements computed with our method [using Eq. (8)] with those obtained with the all-supercell method using Eq. (14), for vacancy defects in silicon (see below for the computational details). It is seen that for a test case of a primitive cell and a 2$\,\times\,$2$\,\times\,$2 supercell the results obtained with the two methods are in perfect agreement. This is but one of many benchmark tests we have performed.

III Results

III.1 Computational details

We apply our approach to compute the $e$-d RTs and defect-limited mobility in silicon, separately for vacancy and (tetrahedral) interstitial defects. A defect concentration of 1 ppm (one defect in 106 atoms) is assumed in both cases. The ground states of the primitive cell and of supercells with size $N\,\times\,N\,\times\,N$ (where $N$ is the number of primitive cells along each lattice vector) are computed using DFT within the local density approximation Perdew and Zunger (1981), using a plane wave basis and norm-conserving pseudopotentials Kleinman and Bylander (1982) with the Quantum Espresso package Giannozzi et al. (2009). Briefly, for the primitive cell we use a lattice constant of $5.43$ Å, a $40$ Ry kinetic energy cutoff and a $12\,\times\,12\,\times\,12$ Monkhorst-Pack $\bm{k}$-point grid Monkhorst and Pack (1976), converging the total energy to within $10$ meV/atom; a consistent lattice constant and total energy convergence criterion is employed for the supercells. In the defect-containing supercells, the atomic forces are relaxed to within 25 meV/Å to account for the structural changes induced by the defect, and the resulting KS potentials are used to compute the $e$-d matrix elements. Due to the different reference potentials in the pristine and defect-containing supercells, we employ the core-average potential alignment method Kumagai and Oba (2014) to align the local potentials of the two supercells when computing the local perturbation potential; the reference potential is taken as the average of the local potential at the atom that is farthest from the defect site. In the $e$-d RT and mobility calculations, we select the electronic states of relevance in a small ($\sim$100 meV) energy window near the band edges, and interpolate the band structure using maximally localized Wannier functions Marzari and Vanderbilt (1997) with the Wannier90 code Yates et al. (2007); Mostofi et al. (2014). All $e$-d calculations have been implemented in our perturbo code Per following the workflow in Fig. 1.

III.2 Relaxation times and their convergence

We first analyze the $e$-d RTs for vacancy defects in silicon, and study their convergence. Figure 3(a) gives the converged RTs (and their inverse, the scattering rates) of electrons and holes due to scattering with vacancy defects with a 1 ppm concentration. As the electron and hole energy increases away from the respective band edges, the scattering rates increase (and the RTs decrease) due to the increased phase space for scattering $-$ that is, the larger number of final states that can be accessed in the scattering process. The latter is quantified by the density of states (DOS), which indeed shows a trend similar to the scattering rates. Overall, the $e$-d RTs are in the ps $-$ ns range, and thus much longer than the typical $e$-ph RTs near room temperature, which are in the fs $-$ ps range Bernardi (2016). This result is consistent with the fact that $e$-ph scattering dominates near room temperature, while $e$-d scattering becomes important at low temperatures, where phonons are frozen out.

We find that the RTs depend on carrier energy, type of carrier, and type of neutral defect, at odds with simplified empirical models employed for decades to model $e$-d scattering Erginsoy (1950). Figure 3(b) shows the RTs as a function of carrier energy in silicon for vacancy and interstitial defects, both with a 1 ppm concentration. For the vacancy defects, the electron RT is roughly 1 ns at the conduction band minimum and the hole RT is roughly 300 ps at the valence band maximum; both these RTs decrease by an order of magnitude $\sim$100 meV away from the band edges. Overall, the holes relax significantly faster than the electrons, both near the band edge and at higher energies, a result we attribute to the higher density of states near the three-fold degenerate valence band maximum Jhalani et al. (2017). A similar asymmetry in the electron and hole RTs has been predicted recently in $e$-ph scattering in GaN Jhalani et al. (2017), where it was also attributed to the valence band degeneracy.

For interstitials, in contrast, we find that the RTs are similar for electrons and holes. They are 150 ps for electrons and 100 ps for holes at the respective band edges, and for both carriers the RTs approach a value of 30 ps roughly 100 meV away from the band edges. The unexpected energy, carrier type, and defect type dependence of the RTs cannot be explained by the widely used Erginsoy formula Erginsoy (1950), which predicts an energy-independent RT for $e$-d scattering due to neutral defects. Different from the Erginsoy model, our ab initio calculations take the atomic and electronic structure into account, providing accurate results that are material and defect specific.

In ab initio $e$-d calculations, there are significant challenges with converging the RTs, which have so far not been examined in detail. This convergence is crucial since many transport properties and physical observables associated with $e$-d interactions depend sensitively on the RTs. Figure 4 shows how to systematically converge the $e$-d RTs with respect to three key factors $-$ the supercell size, structural relaxation, and BZ grid used in the sum over final states in Eq. (1). This convergence study is discussed here for the vacancy case, although we find similar results for interstitial defects.

Figure 4(a) shows the convergence of the RTs with respect to supercell size. Results are given for supercell sizes ranging from $4\,\times\,4\,\times\,4$ to $8\,\times\,8\,\times\,8$, in each case containing one vacancy at the center of the supercell. To isolate the role of supercell size, the atomic structure is not relaxed in these calculations. The RTs in the $4\,\times\,4\,\times\,4$ and $6\,\times\,6\,\times\,6$ supercells are within $20$% and $5$%, respectively, of the $8\,\times\,8\,\times\,8$ supercell results, which can be considered fully converged. To verify that the RTs do not depend on defect position, we compute the RTs for an off-center vacancy that is placed away from the center of a $6\,\times\,6\,\times\,6$ supercell. For the same supercell size, the RTs of the off-center and centered vacancy match exactly, as they should.

Structural relaxation can be extensive around a defect and is expected to play an important role in accurately computing $e$-d interactions. Since it is costly to relax the structure in large supercells, an approximate scheme that retains accurate RTs is desirable. To this end, we compute the RTs in a supercell of a fixed size (here, $6\,\times\,6\,\times\,6$) in which only the atoms up to i-th nearest neighbor shell of the vacancy defect are relaxed, and those that are farther away are kept fixed. Figure 4(b) shows the RTs for structural relaxation up to the 1st, 2nd and 3rd nearest neighbor shell. We find that the RTs are almost converged for structural relaxation within the 2nd shell, and nearly identical to those for relaxation up to the 3rd shell, which can be considered converged since the atomic forces are negligible outside the 3rd nearest neighbor shell. The conclusion is that one needs to relax only a small portion of the atoms around the defect to accurately compute the $e$-d RTs.

Most critical when computing the $e$-d RTs is converging the $\bm{k}^{\prime}$-point grid in the sum over final states in Eq. (1), which is equivalent to converging the grid of transferred momenta, $\bm{q}=\bm{k}^{\prime}-\bm{k}+\bm{G}$. There is a cross-convergence effect between this grid and the energy broadening $\eta$ employed to represent the delta function in Eq. (1), which is implemented as a normalized Gaussian with broadening, $\delta_{\eta}(x)=\frac{1}{\sqrt{2\pi}\eta}e^{-x^{2}/2\eta^{2}}$. The situation is fully analogous to converging the $e$-ph scattering rates Zhou and Bernardi (2016). Briefly, the broadening has to be small enough to not alter the final result, but the smaller the broadening the denser the $\bm{k}^{\prime}$-point BZ grid needed to converge the sum in Eq. (1). Systematic convergence is achieved by starting with a small broadening (say, $\eta\approx 10$ meV) and converging the $\bm{k}^{\prime}$-point BZ grid, and then decreasing the broadening to a smaller value and converging the BZ sum again. At convergence, the RTs do not change upon decreasing the broadening and converging the BZ sum. Note that the $\bm{k}^{\prime}$-point grid can in principle be distinct from the $\bm{k}$-point grid at which the RTs are computed, but this is feasible in practice only if one has a mechanism to effectively interpolate the matrix elements. When this is possible, using random or importance sampling $\bm{k}^{\prime}$-point grids can significantly speed up the calculations Bernardi (2016). Here, in each calculation, we use the same uniform BZ grid for $\bm{k}$- and $\bm{k}^{\prime}$-points, and refer to it below as the BZ grid (a uniform $M\times M\times M$ grid will be denoted as an $M^{3}$ grid).

Figure 4(c) shows the RTs for several values of the energy broadening $\eta$ and gives the corresponding BZ grid at convergence. The BZ grid required to converge the RTs are denser for smaller values of the broadening; for $\eta$ values of 1, 5 and 15 meV, uniform BZ grids with $400^{3}$, $200^{3}$ and $150^{3}$ points are needed, respectively. For electron energies higher than $25$ meV above the conduction band minimum, a $15$ meV broadening and a $150^{3}$ BZ grid are sufficient to converge the RTs. For electron energies within 25 meV of the band edge, a 5 meV broadening with a $200^{3}$ BZ grid gives the same RTs as a smaller 1 meV broadening with a $400^{3}$ BZ grid 111Converging the RTs within a few meV of the band edges may require even smaller broadenings and denser BZ grids, but it is not necessary for mobility calculations, even at very low temperatures, because the carrier velocity vanishes at the band edges.. The broadening and BZ grid values at convergence are similar for electrons and holes, and for vacancies and interstitials.

Importantly, the RTs computed with unconverged grids can differ widely from the converged values, especially at energies near the band edges [see Fig. 4(c)], which critically contribute to charge transport. It is therefore essential to have an efficient method for computing and converging the $e$-d RTs on fine BZ grids to accurately compute charge transport at low temperature.

III.3 Defect-limited carrier mobility

At room temperature, where $e$-ph interactions typically dominate, charge transport can be accurately predicted from first principles in several families of materials Mustafa et al. (2016); Park et al. (2014); Zhou and Bernardi (2016); Lee et al. (2018); Zhou et al. (2018). However, many devices and experiments operate at low temperature, where charge transport is governed by $e$-d interactions. It is thus critically important to develop ab initio calculations that can predict carrier dynamics at low temperature in the presence of $e$-d scattering. We compute the defect-limited carrier mobility $\mu$ at temperature $T$ within the RT approximation of the Boltzmann transport equation Bernardi (2016):

[TABLE]

where $e$ is the electron charge, $n_{c}$ the carrier concentration, $f(T,E)$ the Fermi-Dirac distribution, and $\Sigma\left(E\right)$ the transport distribution function (TDF) at energy $E$,

[TABLE]

where $\alpha$ and $\beta$ are Cartesian directions. The TDF is computed with a tetrahedron integration method Zhou and Bernardi (2016), using converged $e$-d RTs and Wannier-interpolated band velocities $\mathbf{v}_{n\bm{k}}$ Yates et al. (2007); Mostofi et al. (2014). To estimate the carrier energy range contributing significantly to the mobility, we plot the integrand of the mobility formula in Eq. (15), the function $(-\partial f/\partial E)\times\Sigma\left(E\right)$, in Fig. 4(c) for temperatures of $2$, $10$ and $160$ K. As the temperature increases, the peak of the function broadens and moves away in energy from the band edges, indicating that the energy region contributing to the mobility shifts to higher carrier energies. The most stringent conditions for computing the RTs are below $10$ K, where the contribution to the mobility peaks $5$ meV away from the band edge; in this regime, BZ grids as dense as $200^{3}$ $\bm{k}^{\prime}$-points and a broadening of 5 meV are needed to accurately compute the mobility.

We compute the mobility of electrons and holes in silicon, considering separately vacancy and interstitial defects. Figure 5 shows the computed mobility curves at temperatures below 150 K. The electron and hole mobilities for vacancy defects are higher than the corresponding mobilities for interstitial defects due to the longer $e$-d RTs for vacancies. We find that, in all cases, the defect-limited mobility is roughly constant below 10 K, and decreases at higher temperatures. Note that $e$-ph interactions are not included here, so these trends are due solely to the $e$-d interactions. Above 50 K, the temperature dependence of the mobility is well approximated by a power law, $\mu\propto T^{-\alpha}$, with coefficients $\alpha$ of order 0.5. For vacancy defects, the best-fit values of $\alpha$ are 0.53 for electrons and 0.55 for holes, and for interstitials 0.46 for electrons and 0.38 for holes.

An interesting interpretation, which is particularly apt for vacancies, is that the defects can be regarded as a substance added to the pure crystal to make an “alloy”. One thus expects that the temperature dependence of the carrier mobility for defect scattering is similar to that of alloy scattering, for which a power law with $\alpha=0.5$ is expected based on existing models Vaughan et al. (2012); Hamaguchi (2010). Our results above 50 K are consistent with this interpretation, but we additionally find that the value of the exponent $\alpha$ depends on carrier and defect type. By contrast, the Erginsoy formula Erginsoy (1950) predicts a temperature independent mobility for neutral defect scattering, which is clearly inconsistent with our results, and also with experiment.

An early experiment Norton et al. (1973) on $n$-type doped silicon obtained the mobility due to neutral impurity scattering by subtracting the lattice and ionized impurity contributions. The resulting mobility, which is limited by neutral defect scattering alone, decreases above 50 K, in agreement with our results. While this trend has been attributed by the authors to inelastic $e$-d scattering Norton et al. (1973), our results show that it can be explained by accurately computing elastic $e$-d scattering. The mobility decrease is due to the energy dependence of the RTs $-$ as the temperature increases, so does the average energy of the electrons contributing to the mobility [see Fig. 4(c)], and their average RTs decrease as we have shown, causing the mobility decrease with temperature.

The mobility computed here using $e$-d interactions with only neutral point defects can be seen as an upper bound corresponding to an ideally pure material. In silicon, the mobility in the purest crystals (with impurity concentrations of $\sim$1012 cm*-3*) is roughly 106 cm2/Vs at 10 K, a value that can be explained by ionized impurity scattering Canali et al. (1975). In these pure silicon samples, the concentration of neutral defects can also be as low as 1012 cm*-3* ($C_{\rm d}\approx 10^{-10}$ in dimensionless units) Canali et al. (1975). We extend our results to this lower neutral defect concentration, and estimate a mobility limit in silicon of $\sim$1010 cm2/Vs at 10 K for an ideal scenario in which ionized impurity scattering is absent and only neutral point defects scatter the carriers. This mobility limit is higher than the value measured in samples with ionized impurity scattering, since the latter is much stronger than neutral defect scattering due to its long-range character. More extensive quantitative comparisons between computed and measured low-temperature charge transport data will be the subject of future work. We conclude that our approach is a powerful tool to compute charge transport at low temperature and estimate mobility limit values in materials.

IV Discussion

We close by discussing technical remarks and future research directions. Several improvements are possible to our calculations of $e$-d RTs. Similar to $e$-ph calculations, interpolating the $e$-d matrix elements, for example using Wannier functions or atomic orbitals Agapito and Bernardi (2018), would be desirable, since one could compute the primitive cell wavefunctions only on coarse grids rather than on the fine grids needed to converge the RTs. In addition, since the broadening needed to converge the RTs increases with increasing carrier energy, using an adaptive broadening scheme could significantly speed up the RT calculations. One could use relatively coarse BZ grids and larger broadening values at higher carrier energies to save computational time, while using finer grids (and a smaller broadening) only at low carrier energy.

The method presented in this work can be extended in several ways. While our calculations focus on neutral point defects, our method can be generalized to treat charged defects, by adding their long-range Coulomb interaction to the local perturbation potential in reciprocal space, similar to what is done in $e$-ph calculations Zhou and Bernardi (2016). The method can also be extended, using spin-polarized DFT calculations, to treat $e$-d scattering processes involving spin and spin-orbit coupling. Example applications include spin-flip processes due to magnetic impurities and defects in topological materials.

The proposed $e$-d calculations are general since they take into account the atomic structure of the material (including important structural relaxation effects around the defect) and its electronic structure. Unlike empirical models, there is no particular extension needed to treat different types of point and extended defects or different materials, provided one can afford the large DFT calculations needed to obtain the perturbation potentials. For example, our method is suitable for extended defects such as dislocations or grain boundaries, but to study them one may need supercells with thousands of atoms.

While this work focuses on the lowest order of perturbation theory, by including higher order $e$-d interactions one could investigate a wide range of low-temperature phenomena, including weak localization or antilocalization and universal conductance fluctuations Bruus and Flensberg (2004); Mahan (2000). Our method can also be a starting point for efficient inelastic $e$-d scattering calculations Barmparis et al. (2015). Ab initio calculations of $e$-d interactions are still in their infancy, and more work is needed to develop their potential and expand their scope.

V Conclusion

We have presented an efficient approach that overcomes main technical challenges for ab initio calculations of $e$-d interactions. The method is applied to compute and systematically converge the elastic $e$-d RTs and the associated defect-limited carrier mobility below $150$ K for vacancy and interstitial defects in silicon. The RTs exhibit a pronounced dependence on energy, carrier and defect type, and the defect-limited mobility is temperature dependent. These results cannot be explained using widely used empirical models of $e$-d interactions. Our approach can provide new microscopic insight into $e$-d scattering processes. It is broadly applicable and can be generalized to treat charged defects, magnetic impurities, and extended defects. We expect that this work will lay a solid foundation for efficient ab initio calculations of $e$-d interactions.

Acknowledgements.

This work was supported by the Air Force Office of Scientific Research through the Young Investigator Program Grant FA9550-18-1-0280. J.-J. Z. was supported by the National Science Foundation under Grant No. ACI- 1642443, which provided for code development. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. I-T. L. thanks Dr. Luis Agapito, Dr. Davide Sangalli, Vatsal Jhalani, Jinsoo Park and Xiao Tong for fruitful discussions.

APPENDIX A ELASTIC AND INCOHERENT ELECTRON-DEFECT SCATTERING RATE

The scattering rate $\Gamma_{n\bm{k}}$ can be written using Fermi’s golden rule within lowest-order perturbation theory:

[TABLE]

where the perturbation $\Delta\hat{H}$ is the difference between the Kohn-Sham Hamiltonian $\hat{H}^{(\rm d)}$ of a crystal containing $N_{\rm d}$ identical defects and the Hamiltonian $\hat{H}^{(\rm p)}$ of the same crystal with no defects, namely, $\Delta\hat{H}=\hat{H}^{(\rm d)}-\hat{H}^{(\rm p)}$. The crystal is made up by $N_{\bm{k}}$ primitive cells, and we apply Born-von Karman (BvK) periodic boundary conditions; the crystal volume is $\Omega_{\rm BvK}=N_{\bm{k}}\Omega_{\rm uc}$, where $\Omega_{\rm uc}$ is the volume of the primitive cell. Above, $\ket{\psi_{n\bm{k}}}$ are unperturbed Bloch wavefunctions with energy $\varepsilon_{n\bm{k}}$, which in coordinate space read:

[TABLE]

where $\ket{n\bm{k}}$ is the Bloch wavefunction without the prefactor, and $u_{n\bm{k}}\left(\bm{r}\right)$ is the periodic part of the Bloch wavefunction, normalized in the primitive cell as

[TABLE]

Since the kinetic energy is the same in the pristine and defect-containing systems, the difference of their Hamiltonians equals the sum of the perturbations due to all defects:

[TABLE]

where $\Delta V_{\rm e-d}(\bm{r}-\bm{r}_{i})$ denotes the perturbation potential due to a defect located at $\bm{r}_{i}$, and we consider non-interacting defects of the same kind. Assuming that the scattering events are independent, we can write the scattering rate for elastic and incoherent scattering processes due to all defects as

[TABLE]

where $n_{\rm at}$ is the number of atoms in the primitive cell, $C_{\rm d}$, which is formally equal to $N_{\rm d}/\left(N_{\bm{k}}\times n_{\rm at}\right)$, is in practice an assumed value of the defect concentration, and $M_{n^{\prime}\bm{k}^{\prime},n\bm{k}}$ is defined as the $e$-d matrix element for the perturbation due to a single defect:

[TABLE]

Within our approximations, the scattering rate is proportional to the defect concentration, and can be computed at any desired value of the defect concentration, provided that the scattering events remain uncorrelated and the defects non-interacting throughout the concentration range of interest.

APPENDIX B LOCAL MATRIX ELEMENTS

The local matrix elements can be written as:

[TABLE]

We define the forward and inverse Fourier transforms of the local perturbation potential, respectively, as

[TABLE]

and

[TABLE]

where in the last line we replace the crystal volume $\Omega_{\rm BvK}$ with the supercell volume $\Omega_{\rm sup}$, using the fact that the local perturbation potential vanishes at the supercell boundary, which is typically the case for supercells larger than a few primitive cells due to the localized nature of the perturbation potential. Inserting Eq. (24) into Eq. (23) and using the translational invariance of $u_{n\bm{k}}\left(\bm{r}\right)$, we have

[TABLE]

with the plane-wave matrix elements defined as

[TABLE]

where $\bm{G}$ are reciprocal lattice vectors of the primitive cell. This formula is valid for any basis set. Here we use a plane wave basis for $u_{n\bm{k}}\left(\bm{r}\right)$, and write

[TABLE]

where $C_{n\bm{k}}\left(\bm{G}\right)$ is the Fourier coefficient of $u_{n\bm{k}}\left(\bm{r}\right)$ at the reciprocal lattice vector $\bm{G}$. The local matrix element formula in a plane wave basis becomes:

[TABLE]

This is the formula implemented in our code and used in this work.

APPENDIX C NONLOCAL MATRIX ELEMENTS

To obtain an expression for the nonlocal matrix elements, we focus on the scalar product

[TABLE]

We first fix the atomic position at the origin (by setting $\tau_{s}=0$) and then generalize the result to arbitrary atomic positions. We define generalized forward and inverse Fourier transforms of the Kleinman-Bylander (KB) projectors, respectively, as:

[TABLE]

and

[TABLE]

where we replace the crystal volume $\Omega_{\rm BvK}$ with the supercell volume $\Omega_{\rm sup}$ because the KB projector is localized around the core atomic region. If a general atomic position $\mathbf{\tau}_{s}$ is chosen, the Fourier coefficient $B_{j\bm{k}}^{(s)}\left(\bm{q}\right)$ becomes, using the properties of the Fourier transforms,

[TABLE]

The scalar product in Eq. (30), after inserting Eq. (31) into Eq. (30), becomes

[TABLE]

where we used the translational invariance of $u_{n\bm{k}}\left(\bm{r}\right)$. This formula is valid for any basis set. Here we use a plane wave basis, so the matrix elements of $\hat{V}_{\rm NL}$, for each of the pristine and defect-containing supercells (labelled by $\alpha$ = d, p), read:

[TABLE]

where we use the same notation as in Eq. (11). The nonlocal matrix elements are computed as the difference in Eq. (10). This is the formula implemented in our code and used in this work.