Introduction

The interactions between electrons and defects control charge and spin transport at low temperature, and give rise to a range of quantum transport phenomena1,2,3,4,5. Understanding such electron–defect (e–d) interactions from first principles can provide microscopic insight into carrier dynamics in materials in the presence of point and extended defects, and accelerate materials discovery. Currently, e–d interactions can be computed ab initio mainly through the all-electron Korringa–Kohn–Rostoker Green’s function method6, which is based on multiple-scattering theory. Several properties have been computed with this approach, including the magnetoresistance due to impurity scattering7, the residual resistivity of metals and alloys8, and spin relaxation9,10,11,12 and the spin Hall effect13,14,15. In contrast, ab initio e–d calculations using pseudopotentials or projector augmented waves16,17,18 have progressed more slowly, mainly due to the high computational cost of obtaining the e–d interaction matrix elements needed for perturbative calculations19.

We recently developed an ab initio method18 to compute efficiently the e–d interactions and the associated matrix elements. Our approach uses only the wave functions of the primitive cell, thus significantly reducing computational cost compared to e–d calculations that use supercell wave functions16,17. As our method uses a plane-wave (PW) basis set and pseudopotentials, it is compatible with widely used density functional theory (DFT) codes. A different method developed by Kaasbjerg et al.20 uses an atomic orbital (AO) basis to compute the e–d matrix elements; its advantage is that one can compute the e–d matrix elements using only a small set of AOs, although one is limited by the quality and completeness of the AO basis set21.

To benefit from both the completeness and accuracy of the PW basis and the versatility of a small localized basis set, approaches combining PWs and AOs22 or Wannier functions (WFs)23,24 have been developed for electron–phonon (e–ph) interactions. They have enabled efficient interpolation of the e–ph matrix elements and have been instrumental to advancing carrier dynamics calculations25,26,27. To date, such an interpolation scheme does not exist for e–d interactions to our knowledge; thus, performing demanding Brillouin zone (BZ) integrals needed to compute e–d relaxation times (RTs) and defect-limited charge transport remains an open problem. Interpolating the interaction matrix elements to uniform, random, or importance-sampling fine BZ grids is key to systematically converging the RTs and transport properties25,28, and it has been an important development in first-principles calculations of e–ph interactions and phonon-limited charge transport.

In this work, we develop a method for interpolating the e–d interaction matrix elements using WFs. Through a generalized double-Fourier transform, our approach can efficiently transform the matrix elements from a Bloch representation on a coarse BZ grid to a localized WF representation and ultimately to a Bloch representation on an arbitrary fine BZ grid. We show the rapid spatial decay of the e–d interactions in the WF basis, which is crucial to the accuracy and efficiency of the method. Using our approach, we investigate e–d interactions due to charge-neutral vacancies in silicon and copper. In both cases, we can accurately interpolate the e–d matrix elements and converge the e–d scattering rates and defect-limited carrier mobility or resistivity. In copper, we map the e–d RTs directly on the Fermi surface and show their dependence on electronic state. We demonstrate computations of e–d matrix elements on random and uniform BZ grids as dense as 600 × 600 × 600 points, whose computational cost would be prohibitive for direct computation. We expect that our efficient e–d computation and interpolation approaches will enable a wide range of studies of e–d interactions in materials ranging from metals to semiconductors and insulators, and for applications such as electronics, nanodevices, spintronics, and quantum technologies.

Results

Theory

The perturbation potential ΔVe−d introduced by a point defect in a crystal couples different Bloch eigenstates of the unpertured (defect-free) crystal. The matrix elements associated with this e–d interaction are defined as

$${M}_{\mathrm {mn}}({\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}})=\left\langle m{\boldsymbol{k}}^{\prime} \ |\Delta {V}_{{{e}}-{\rm{d}}}| \ n{\boldsymbol{k}}\right\rangle ,$$
(1)

where \(\left|n{\boldsymbol{k}}\right\rangle\) is the Bloch state with band index n and crystal momentum k. To handle these e–d interactions, one needs to store and manipulate a matrix Mmn of size \({N}_{\mathrm {b}}^{2}\) (Nb is the number of bands) for each pair of crystal momenta \({\boldsymbol{k}}^{\prime}\) and k in the BZ. Within DFT, the perturbation potential ΔVe−d can be computed as the difference between the Kohn–Sham potential of a defect-containing supercell and that of a pristine supercell with no defect18.

We compute the e–d matrix elements in Eq. (1) with the method we developed in ref. 18, which uses only the Bloch wave functions of the primitive cell and does not require computing or manipulating the wave functions of the supercell, thus significantly reducing computational cost. The Bloch states can be expressed in terms of maximally localized WFs using

$$\left|n{\boldsymbol{k}}\right\rangle ={\sum \limits_{j{\boldsymbol{R}}}}{e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{R}}}{U}_{jn,{\boldsymbol{k}}}^{\dagger }\left|j{\boldsymbol{R}}\right\rangle ,$$
(2)

where \(\left|j{\boldsymbol{R}}\right\rangle\) is the WF with index j centered at the Bravais lattice vector R. The unitary matrices U in Eq. (2) maximize the spatial localization of the WFs29

$$\left|j{\boldsymbol{R}}\right\rangle =\frac{1}{{N}_{{\boldsymbol{k}}}}{\sum \limits_{n{\boldsymbol{k}}}}{e}^{-i{\boldsymbol{k}}\cdot {\boldsymbol{R}}}{U}_{nj,{\boldsymbol{k}}}\left|n{\boldsymbol{k}}\right\rangle ,$$
(3)

where Nk is the number of k-points in the BZ. The e–d matrix element \({M}_{ij}({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})\) between two WFs centered at the lattice vectors \({\boldsymbol{R}}^{\prime}\) and R in the Wannier representation is defined as

$${M}_{ij}({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})=\left\langle i{\boldsymbol{R}}^{\prime} \ |\Delta {V}_{{{e}}-{\rm{d}}}| \ j{\boldsymbol{R}}\right\rangle .$$
(4)

If the center of the perturbation potential lies in the unit cell at the origin, the absolute value of the e–d matrix elements, \(| {M}_{ij}({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})|\), decays rapidly (within a few lattice constants) for increasing values of the lattice vectors \({\boldsymbol{R}}^{\prime}\) and R due to the short-ranged nature of the perturbation potential from the defect, which is assumed to be charge-neutral in this work. As a result, only a small number of lattice vectors R, which we arrange in a Wigner–Seitz (WS) supercell centered at the origin, is needed to compute the e–d matrix elements in the Wannier representation.

Using Eqs. (1)–(4), the e–d matrix elements in the Wannier representation can be written as a generalized double-Fourier transform of the matrix elements in the Bloch representation, which are first computed on a coarse BZ grid with points kc:

$$M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})={\left(\frac{1}{{N}_{{{\boldsymbol{k}}}_{{\rm{c}}}}}\right)}^{2}{\sum \limits_{{{\boldsymbol{k}}}_{{\rm{c}}}^{\prime}{{\boldsymbol{k}}}_{{\rm{c}}}}}{e}^{i({{\boldsymbol{k}}}_{{\rm{c}}}^{\prime}\cdot {\boldsymbol{R}}^{\prime} -{{\boldsymbol{k}}}_{{\rm{c}}}\cdot {\boldsymbol{R}})}{U}_{{{\boldsymbol{k}}}_{{\rm{c}}}^{\prime}}^{\dagger }M({{\boldsymbol{k}}}_{{\rm{c}}}^{\prime},{{\boldsymbol{k}}}_{{\rm{c}}}){U}_{{{\boldsymbol{k}}}_{{\rm{c}}}}.$$
(5)

Here and below, we omit all band indices for clarity. Through the inverse transform, we can interpolate the e–d matrix elements to any desired pair of fine BZ grid points \({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime}\) and kf, using

$$M({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime},{{\boldsymbol{k}}}_{{\rm{f}}})={\sum \limits_{{\boldsymbol{R}}^{\prime} {\boldsymbol{R}}}}{e}^{-i({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime}\cdot {\boldsymbol{R}}^{\prime} -{{\boldsymbol{k}}}_{{\rm{f}}}\cdot {\boldsymbol{R}})}{U}_{{{\boldsymbol{k}}}_{{\rm{f}}}^{\prime}}M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}}){U}_{{{\boldsymbol{k}}}_{{\rm{f}}}}^{\dagger }.$$
(6)

Although \({U}_{{{\boldsymbol{k}}}_{{\rm{c}}}}\) in Eq. (5) is the coarse-grid unitary matrix used to construct the WFs [see Eq. (3)], the unitary matrix on the fine grid, \({U}_{{{\boldsymbol{k}}}_{{\rm{f}}}}\) in Eq. (6), is obtained by diagonalizing the fine-grid Hamiltonian,

$$H({{\boldsymbol{k}}}_{{\rm{f}}})={\sum \limits_{{\boldsymbol{R}}}}{e}^{i{{\boldsymbol{k}}}_{{\rm{f}}}\cdot {\boldsymbol{R}}}H({\boldsymbol{R}}),$$
(7)

where H(R) is the electronic Hamiltonian in the Wannier basis. These equations are analogous to those used for interpolating the e–ph matrix elements22,23,24, except that here the lattice vectors \({\boldsymbol{R}}^{\prime}\) and R are both associated with electronic states. The lattice vectors \({\boldsymbol{R}}^{\prime}\) and R in the WS supercell are determined—through the periodic boundary conditions—by the \({{\boldsymbol{k}}}_{{\rm{c}}}^{\prime}\) and kc coarse grids, respectively. In practice, we choose a uniform coarse BZ grid and the size of the WS supercell is equal to the size of this coarse grid. It is noteworthy that our e–d matrix elements in the Wannier representation require WFs only for the primitive cell, whereas a method developed in ref. 30 requires WFs for the defect-containing supercells.

Similar to the e–ph case22, the rapid spatial decay of the e–d matrix elements in the WF basis is crucial to reducing the computational cost, as it puts an upper bound to the number of lattice sites \({\boldsymbol{R}}^{\prime}\) and R at which \(M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})\) needs to be computed. In particular, computing \(M({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime},{{\boldsymbol{k}}}_{{\rm{f}}})\) at small \({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime}\) and kf vectors would in principle require summing the Fourier transform in Eq. (6) up to correspondingly large lattice vectors of length \(| {\boldsymbol{R}}^{\prime} | =2\pi /| {{\boldsymbol{k}}}_{{\rm{f}}}^{\prime}|\) and ∣R∣ = 2π∕∣kf∣; in practice, this is not needed due to the rapid spatial decay. The choice of a WS supercell and its relation to the DFT supercell are discussed in detail in the Supplementary Information.

Interpolation workflow and validation

The workflow for interpolating the e–d matrix elements to a fine grid with points kf consists of several steps (see Fig. 1) as follows: (i) compute the e–d matrix elements in the Bloch representation on a coarse BZ grid with points kc using Eq. (1); (ii) obtain the e–d matrix elements in the Wannier representation using Eq. (5); (iii) interpolate the Hamiltonian using Eq. (7) and diagonalize it to obtain the fine-grid unitary matrices \({U}_{{{\boldsymbol{k}}}_{{\rm{f}}}}\); (iv) interpolate the e–d matrix elements to any desired pair of fine-grid points \({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime}\) and kf using the matrix elements in the Wannier representation and the fine-grid unitary matrices [see Eq. (6)].

Fig. 1: Workflow for interpolating the e–d matrix elements using WFs.
figure 1

The steps are numbered as in the text.

Full size image

We validate our WF-based interpolation method using vacancy defects in silicon as an example (see Methods). The relaxed symmetry of the vacancy in our calculation is Td, whereas previous work and experiment find a D2d symmetry31. The reason for this inconsistency is that we do not randomly displace the atoms before relaxation, which is needed to break the symmetry and obtain the lower-energy D2d vacancy structure. Although our calculations use a vacancy with Td rather than D2d symmetry, the results we present are not affected by this choice. Figure 2a compares the e–d matrix elements calculated directly using Eq. (1) with the same matrix elements obtained by interpolation starting from two different coarse BZ grids with respectively 43 and 103 points kc (here and below, we denote an N × N × N uniform grid as N3). The interpolated results can qualitatively reproduce the direct computation for both coarse grids, but the results from the 103 coarse grid achieve a superior quantitative accuracy as the interpolated matrix elements agree with the directly computed ones within 1% over the entire BZ. The accuracy can be systematically improved by increasing the size of the coarse kc-grid (see Table 1). This trend implies that the e–d perturbation potential decays to a negligible value over >4 but <10 lattice constants.

Fig. 2: Matrix elements and their spatial decay.
figure 2

a Absolute value of the e–d matrix elements, computed along high-symmetry BZ lines. The initial state is set to the lowest valence band at Γ, whereas the final states span all four valence bands and possess crystal momenta chosen along the high-symmetry lines shown in figure. Panels b and c show the spatial decay of the e–d matrix elements in the Wannier basis, \(| | M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})| |\), which are plotted in b as a function of \(| {\boldsymbol{R}}^{\prime} |\) for R = 0 and in c as a function of ∣R∣ for \({\boldsymbol{R}}^{\prime} ={\bf{0}}\). The highest value of \(| | M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})| |\) is normalized to 1 in both cases and the plots use a logarithmic scale.

Full size image
Table 1 The mean and maximum difference between the interpolated and directly computed e–d matrix elements, given for several coarse kc-grids.
Full size table

The spatial decay of the matrix elements in the WF basis is essential for the accuracy of our approach. To analyze the spatial behavior of the matrix elements in the WF basis, we define for each pair of lattice vectors \({\boldsymbol{R}}^{\prime}\) and R the maximum absolute value of the e–d matrix elements as \(| | M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})| | ={\max }_{ij}| {M}_{ij}({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})|\). Figure 2b, c show the spatial behavior of \(| | M({\boldsymbol{R}}^{\prime} ,{\boldsymbol{R}})| |\) for a vacancy defect in silicon as a function of \(| {\boldsymbol{R}}^{\prime} |\) while keeping R = 0 and as a function of ∣R∣ while keeping \({\boldsymbol{R}}^{\prime} ={\bf{0}}\), respectively. Note that \(| M({\boldsymbol{R}}^{\prime} ,0)|\) and ∣M(0, R)∣ in Fig. 2b, c are identical, because the defect perturbation potential is centered at the origin of the WS supercell. This result does not hold in general for an arbitrary position of the defect in the WS supercell. We find that the matrix elements in the WF basis decay exponentially over a few unit cells, thus confirming that the WFs are a suitable basis set for interpolating the e–d matrix elements.

Relaxation times and defect-limited transport

The e–d RTs and the defect-limited carrier mobility are key to characterizing carrier dynamics at low temperatures and also near room temperature in highly doped or disordered materials. We compute the e–d RTs, τnk, associated with elastic carrier-defect scattering using the lowest-order Born approximation3,18:

$${\tau }_{n{\boldsymbol{k}}}^{-1}=\frac{2\pi }{\hslash }\frac{{n}_{\text{at}}{C}_{{\rm{d}}}}{{N}_{{\boldsymbol{k}}^{\prime} }}{\sum \limits_{m{\boldsymbol{k}}^{\prime}} }{\left|{M}_{mn}({\boldsymbol{k}}^{\prime} ,{\boldsymbol{k}})\right|}^{2}\delta ({\varepsilon }_{m{\boldsymbol{k}}^{\prime} }-{\varepsilon }_{n{\boldsymbol{k}}}),$$
(8)

where ℏ is the reduced Planck constant, nat the number of atoms in a primitive cell, Cd the (dimensionless) defect concentration (defined as the number of defects divided by the number of atoms), \({N}_{{\boldsymbol{k}}^{\prime} }\) the number of k-points used in the summation, and εnk the unperturbed energy of the Bloch state \(\left|n{\boldsymbol{k}}\right\rangle\) in the primitive cell. The delta function is implemented as a normalized Gaussian with a small broadening η, \({\delta }_{\eta }(x)={e}^{-{x}^{2}/2{\eta }^{2}}/\sqrt{2\pi }\eta\). The e–d RTs are proportional to the defect concentration because our approach assumes that the scattering events are independent and uncorrelated18.

For defect-limited carrier transport, we first compute the conductivity tensor σ(T) at temperature T using25

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

where e is the electron charge and E the electron energy, f(T, E) the Fermi-Dirac distribution, and \(\Sigma \left(E\right)\) the transport distribution function (TDF) at energy E, defined as

$${\Sigma }_{\alpha \beta }\left(E\right)=\frac{2}{{\Omega }_{{\rm{uc}}}}{\sum \limits_{n{\boldsymbol{k}}}}{\tau }_{n{\boldsymbol{k}}}{{\bf{v}}}_{n{\boldsymbol{k}}}^{\alpha }{{\bf{v}}}_{n{\boldsymbol{k}}}^{\beta }\delta \left(E-{\varepsilon }_{n{\boldsymbol{k}}}\right),$$
(10)

where α and β are Cartesian directions, and Ωuc the volume of the primitive cell. The TDF is computed with a tetrahedron integration method25, using our calculated e–d RTs and Wannier-interpolated band velocities vnk32,33. The mobility is obtained as μ = σ∕nce, where nc is the carrier concentration, whereas the resistivity is obtained by inverting the conductivity tensor.

We first study the e–d RTs and hole carrier mobility in silicon with vacancy defects (see Methods). We use interpolated e–d matrix elements and focus on the accuracy and convergence of our interpolation method. The results given here assume a defect concentration of one vacancy in 106 silicon atoms (1 p.p.m. concentration), but results for different defect concentrations can be obtained by rescaling these reference RTs to a different defect concentration [using Eq. (8)]; this approach is valid only within the concentration range in which the Born approximation holds. Figure 3a compares the RTs obtained from directly computed and interpolated e–d matrix elements. The two sets of RTs are in close agreement with each other for both electrons and holes, confirming the accuracy of our interpolation method. Figure 3b shows the convergence of the defect-limited hole mobility with respect to the size of the fine BZ grids, for BZ grids ranging from 403 to 6003 points; the convergence is studied at two temperatures (10 K and 100 K) and for two types of grids, random and uniform. At 10 K, the mobilities are fully converged for fine grids with 2003 points, for both random and uniform grids. We observe a similar trend at 100 K. The converged values of the mobility are consistent with our previous calculations using directly computed (rather than interpolated) matrix elements18.

Fig. 3: Carrier relaxation times and the hole mobility for vacancy defects in silicon.
figure 3

a Electron–defect RTs obtained from directly computed and interpolated e–d matrix elements. Here, εc is the conduction band minimum and εv the valence band maximum. b Convergence of the hole mobility with respect to the size of the fine BZ grid used for interpolation, shown for both uniform and random grids. The computed hole mobilities at 10 K and 100 K from ref. 18 are also shown. A reference vacancy concentration of 1 p.p.m. is used in all calculations.

Full size image

The interpolation method allows us to use extremely dense BZ grids with up to 6003 points due to its superior computational efficiency. Let us briefly analyze the overall speed-up of the interpolation method for a carrier mobility calculation. To converge the mobility at low temperature, one needs to consider only fine BZ grid points in a small energy window, roughly within 100 meV of the band edges in semiconductors18,25 or of the Fermi energy in metals; these are the only states contributing to the conductivity in Eq. (9). In this small energy window, the number of k-points is a small fraction α of the total number of points \({N}_{{{\boldsymbol{k}}}_{{\rm{f}}}}\) in the entire fine BZ grid. In the direct computation, one computes the e–d matrix elements \(M({{\boldsymbol{k}}}_{{\rm{f}}}^{\prime},{{\boldsymbol{k}}}_{{\rm{f}}})\) between all crystal momentum pairs and thus a number of matrix elements of order \({(\alpha {N}_{{{\boldsymbol{k}}}_{{\rm{f}}}})}^{2}\). In the interpolation approach, the most time-consuming step is directly computing the \({N}_{{{\boldsymbol{k}}}_{{\rm{c}}}}^{2}e\)-d matrix elements on the coarse grid [step (i) in Fig. 1], whereas interpolating the matrix elements per se is orders of magnitude less computationally expensive. In our machine, the average central processing unit (CPU) time to directly compute one matrix element is ~0.2 s for silicon (with an electronic kinetic energy cutoff of 40 Ry), whereas the same calculation done with our interpolation method requires only ~80 μs. The speed-up of the interpolation scheme for computing matrix elements is thus around three orders of magnitude. As a result, the overall speed-up of the interpolation approach over the direct computation is \(\sim\! {(\alpha {N}_{{{\boldsymbol{k}}}_{{\rm{f}}}})}^{2}/{N}_{{{\boldsymbol{k}}}_{{\rm{c}}}}^{2}\). The typical value of \({N}_{{{\boldsymbol{k}}}_{{\rm{f}}}}\) is around 103 – 104\({N}_{{{\boldsymbol{k}}}_{{\rm{c}}}}\). For our silicon calculations, the value of α is of order 10−2, so the interpolation approach speeds up the mobility calculation by at least two to four orders of magnitude.

The method for directly computing the e–d matrix elements we developed in ref. 18 and the interpolation method shown here are general, and can be applied to metals, semiconductors, and insulators. As an example, we show a calculation on a metal, copper, containing vacancy defects (see Methods). In metals, the fine grids required to compute the e–d RTs near the Fermi energy and the resistivity are a major challenge for direct e–d calculations without interpolation. Figure 4a shows the e–d RTs computed at k-points on the Fermi surface of copper, using interpolated e–d matrix elements obtained from a moderate size (8 × 8 × 8) coarse kc-grid. The e–d RTs for a reference vacancy concentration of 1 p.p.m. are between 0.3 and 1.4 ps. Interestingly, these e–d RTs are orders of magnitude shorter than the electron RTs in silicon for the same vacancy concentration, which suggests that scattering due to vacancy defects in copper is significantly stronger than in silicon (see the Supplementary Information). The e–d RTs in copper are strongly state-dependent—we find values of order 0.3 ps on the majority of the Fermi surface and values as large as 1.4 ps near the regions of the Fermi surface close to the X points of the BZ. Similar state-dependent RTs due to impurity scattering in copper have been predicted using the all-electron Korringa–Kohn–Rostoker Green’s function method9,10. These results show that our first-principles approach can access microscopic details of the e–d scattering processes.

Fig. 4: Relaxation times and the defect-limited resistivity for vacancy defects in copper.
figure 4

a RTs mapped on the Fermi surface, obtained using interpolated e–d matrix elements for a reference vacancy concentration of 1 p.p.m. b Defect-limited resistivity as a function of temperature for an assumed reference vacancy concentration of 1 p.p.m., compared with experimental data from ref. 34.

Full size image

Figure 4b shows the calculated defect-limited resistivity for vacancy defects in copper, at low temperatures between 2 and 50 K (see Methods). The calculated resistivity is independent of temperature, in agreement with experimental results34,35, even though the conductivity formula we use [Eq. (9)] depends on temperature via the Fermi-Dirac distribution. The low-temperature defect-limited resistivity for a reference 1 p.p.m. vacancy concentration is ~10−9 Ω⋅m. However, the equilibrium vacancy concentration (see Methods) of copper at 50 K is negligible (of order 10−119), and the corresponding resistivity is of order 10−122 Ω⋅m. This value is negligible in comparison with the measured resistivity of ~10−12 Ω⋅m in a highly pure copper Cu(7N) sample at low temperature, where 7N means 99.99999% purity (see Fig. 4b). We conclude that the resistivity of real copper samples at low temperature is not limited by intrinsic vacancy defects, but rather is controlled by impurities. This is a well-known result36,37,38,39,40 and our calculations are consistent with it.

Our method can predict a lower bound of the residual resistivity due to intrinsic defects in an ideally pure material. Alternatively, if the main type of defect or impurity is known from experiment, our approach can estimate the defect concentration present in the sample. The so-called residual resistivity ratio (RRR) between the low temperature and room temperature resistivities is used as a figure of merit for sample quality, and a large collection of data exists41 for RRR in metals. As at room temperature the resistivity is usually phonon-limited, combined with e–ph calculations25,42, our approach allows one to compute RRR for a wide range of materials and defect types. Taken together, these capabilities expand the tool box of first-principles methods for investigating carrier dynamics in complex materials.

Discussion

Our e–d interpolation method can be extended to charged defects, for which the long-range Coulomb interactions can be added in reciprocal space similar to what is done for e–ph interactions25,43. By including spin–orbit coupling, our approach can also be extended to study spin-flip processes and e–d interactions in magnetic and topological materials. Future work will also attempt to include multiple e–d scattering events and higher-order e–d interactions, e.g., using the T-matrix formalism20,44. Other applications include investigating carrier scattering due to extended defects such as dislocations and grain boundaries, a topic of prime relevance to materials science. In summary, the applications of the method discussed in this work are broad and so are its possible future extensions.

In conclusion, we developed a WF-based interpolation approach to efficiently compute e–d interactions and the associated matrix elements on fine BZ grids. We have shown that the interpolation method is accurate and that it can effectively compute demanding BZ integrals requiring up to 108–109 k-points. The ability to efficiently interpolate e–d matrix elements starting from moderate BZ coarse grids is a stepping stone toward perturbative calculations of defect-limited charge and spin transport, and to investigate quantum transport regimes governed by e–d interactions.

Methods

DFT calculations

The ground state of a primitive cell and of supercells with size N × N × N (where N is the number of primitive cells along each lattice vector) are computed using DFT within the local density approximation. We use a PW basis set and norm-conserving pseudopotentials45 with the QUANTUM ESPRESSO code46. The total energy is converged to within 10 meV/atom in all structures. In the defect-containing supercells, the atomic forces are relaxed to within 25 meV/Å to account for structural changes induced by the defect. For silicon, we use an experimental lattice constant of 5.43 Å and a PW kinetic energy cutoff of 40 Ry. For copper, we use an experimental lattice constant of 3.61 Å and a PW kinetic energy cutoff of 90 Ry. We employ a 12 × 12 × 12 k-point grid47 for the primitive cells of both materials to converge the charge density and total energy, and interpolate their band structures using maximally localized WFs29 with the Wannier90 code32,33. Coarse kc-grids between 4 × 4 × 4 and 10 × 10 × 10, as given in the text, are used in the non-self-consistent DFT calculations and in the wannierization procedure. The dense kf-grid used in the matrix element interpolation or relaxation time (RT) calculations is unrelated to the WF generation. In silicon, we wannierize the four highest valence and four lowest conduction bands together, using sp3 orbitals centered on the silicon atoms as the initial guess, to compute the e–d matrix elements used for electron and hole RTs. A different wannierization is employed to compute the e–d matrix elements along high-symmetry BZ lines in Fig. 2; in this case, we wannierize only the four highest valence bands, using four s orbitals centered at the fractional coordinates (−1/8, 3/8, −1/8), (−1/8, −1/8, −1/8), (3/8, −1/8, −1/8), and (−1/8, −1/8, 3/8). In copper, we exclude the four lowest (core) 3sp-bands and wannierize the next seven bands using five d orbitals centered at the copper atom and two s orbitals centered at the fractional coordinates (1/4, 1/4, 1/4) and (−1/4, −1/4, −1/4) as the initial guess.

Electron–defect matrix elements, RTs, and resistivity calculations

The methods to directly compute the e–d matrix elements and from them obtain the RTs, mobility, and conductivity (or resistivity) are described in detail in ref. 18. Briefly, we compute the coarse-grid e–d matrix elements using the wave functions of the primitive cell and obtain the perturbation potential due to a vacancy defect using a 6 × 6 × 6 supercell with a 2 × 2 × 2 k-point grid. The atomic positions around the vacancy are relaxed up to the third nearest-neighbor shell in both silicon and copper. The potential alignment for the supercell containing the relaxed vacancy is chosen as the core-averaged potential of the farthest atom from the vacancy in the same (but unrelaxed) supercell. In Fig. 3a, the directly computed matrix elements are calculated using the wave functions obtained from non-self-consistent calculations on a dense grid with 3003 kf-points, whereas the interpolated matrix elements are calculated on the same dense grid starting from a coarse grid with 103 kc-points. In the e–d RT and defect-limited mobility calculations in silicon, we use only electronic states in a small (~100 meV) energy window near the band edges, as these are the only states contributing to the mobility25; similarly, in copper we use only states within 100 meV of the Fermi energy. In silicon, we use a broadening value η = 5 meV to compute the delta function in Eq. (8) and a uniform BZ grid with 3003 points for the RTs; for the mobility, we use a 1 meV broadening and e–d matrix elements interpolated from a 103 coarse BZ grid. In copper, we compute the RTs and resistivity on a fine BZ grid with 2403 points, using a 1 meV broadening and e–d matrix elements interpolated from a coarse 83 BZ grid. The equilibrium vacancy concentration at temperature T in copper is estimated using48

$${C}_{{\rm{v}}}(T)={e}^{-(\Delta {H}_{{\rm{v}}}-T\Delta {S}_{{\rm{v}}})/{k}_{{\rm{B}}}T},$$
(11)

where ΔHv and ΔSv are the vacancy formation enthalpy and entropy, respectively, and kB the Boltzmann constant. In copper, ΔSv = 3.0 kB and ΔHv = 1.19 eV48.