Ab initio electron-defect interactions using Wannier functions

Computing electron-defect (e-d) interactions from first principles has remained impractical due to computational cost. Here we develop an interpolation scheme based on maximally localized Wannier functions (WFs) to efficiently compute e-d interaction matrix elements. The interpolated matrix elements can accurately reproduce those computed directly without interpolation, and the approach can significantly speed up calculations of e-d relaxation times and defect-limited charge transport. We show example calculations of vacancy defects in silicon and copper, for which we compute the e-d relaxation times on fine uniform and random Brillouin zone grids (and for copper, directly on the Fermi surface) as well as the defect-limited resistivity at low temperature. Our interpolation approach opens doors for atomistic calculations of charge carrier dynamics in the presence of defects.


INTRODUCTION
The interactions between electrons and defects control charge and spin transport at low temperature, and give rise to a range of quantum transport phenomena [1][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 method 6 , which is based on multiple-scattering theory. Several properties have been computed with this approach, including the magnetoresistance due to impurity scattering 7 , the residual resistivity of metals and alloys 8 , and spin relaxation [9][10][11][12] and the spin Hall effect [13][14][15] . In contrast, ab initio e-d calculations using pseudopotentials or projector augmented waves [16][17][18] have progressed more slowly, mainly due to the high computational cost of obtaining the e-d interaction matrix elements needed for perturbative calculations 19 .
We recently developed an ab initio method 18 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 functions 16,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 set 21 .
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 AOs 22 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 calculations [25][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 properties 25,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.

Theory
The perturbation potential ΔV e−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 where nk j i 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 M mn of size N 2 b (N b is the number of bands) for each pair of crystal momenta k 0 and k in the BZ. Within DFT, the perturbation potential ΔV e−d can be computed as the difference between the Kohn-Sham potential of a defectcontaining supercell and that of a pristine supercell with no defect 18 .
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 where jR j i 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 WFs 29 where N k is the number of k-points in the BZ. The e-d matrix element M ij ðR 0 ; RÞ between two WFs centered at the lattice vectors R 0 and R in the Wannier representation is defined as If the center of the perturbation potential lies in the unit cell at the origin, the absolute value of the e-d matrix elements, jM ij ðR 0 ; RÞj, decays rapidly (within a few lattice constants) for increasing values of the lattice vectors R 0 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 k c : 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 k 0 f and k f , using Although U kc 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 k f in Eq. (6), is obtained by diagonalizing the fine-grid Hamiltonian, where H(R) is the electronic Hamiltonian in the Wannier basis. These equations are analogous to those used for interpolating the e-ph matrix elements [22][23][24] , except that here the lattice vectors R 0 and R are both associated with electronic states. The lattice vectors R 0 and R in the WS supercell are determined-through the periodic boundary conditions-by the k 0 c and k c 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 defectcontaining supercells.
Similar to the e-ph case 22 , 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 R 0 and R at which MðR 0 ; RÞ needs to be computed. In particular, computing Mðk 0 f ; k f Þ at small k 0 f and k f vectors would in principle require summing the Fourier transform in Eq. (6) up to correspondingly large lattice vectors of length jR 0 j ¼ 2π=jk 0 f j and |R| = 2π∕|k f |; 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 k f 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 k c 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 k f ; (iv) interpolate the e-d matrix elements to any desired pair of fine-grid points k 0 f and k f using the matrix elements in the Wannier representation and the fine-grid unitary matrices [see Eq. (6)].
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 T d , whereas previous work and experiment find a D 2d symmetry 31 . 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 D 2d vacancy structure. Although our calculations use a vacancy with T d rather than D 2d 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 4 3 and 10 3 points k c (here and below, we denote an N × N × N uniform grid as N 3 ). The interpolated results can qualitatively reproduce the direct computation for both coarse grids, but the results from the 10 3 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 k cgrid (see Table 1). This trend implies that the e-d perturbation potential decays to a negligible value over >4 but <10 lattice constants.
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 R 0 and R the maximum absolute value of the e-d matrix elements as jjMðR 0 ; RÞjj ¼ max ij jM ij ðR 0 ; RÞj. Figure 2b, c show the spatial behavior of jjMðR 0 ; RÞjj for a vacancy defect in silicon as a function of jR 0 j while keeping R = 0 and as a function of |R| while keeping R 0 ¼ 0, respectively. Note that jMðR 0 ; 0Þj 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 approximation 3,18 : where ℏ is the reduced Planck constant, n at the number of atoms in a primitive cell, C d the (dimensionless) defect concentration (defined as the number of defects divided by the number of atoms), N k 0 the number of k-points used in the summation, and ε nk the unperturbed energy of the Bloch state nk j i in the primitive cell. The delta function is implemented as a normalized Gaussian with a small broadening η, δ η ðxÞ ¼ e Àx 2 =2η 2 = ffiffiffiffiffi ffi 2π p η. The e-d RTs are proportional to the defect concentration because our approach assumes that the scattering events are independent and uncorrelated 18 .
For defect-limited carrier transport, we first compute the conductivity tensor σ(T) at temperature T using 25 where e is the electron charge and E the electron energy, f(T, E) the Fermi-Dirac distribution, and Σ E ð Þ the transport distribution function (TDF) at energy E, defined as where α and β are Cartesian directions, and Ω uc the volume of the primitive cell. The TDF is computed with a tetrahedron integration method 25 , using our calculated e-d RTs and Wannier-interpolated band velocities v nk 32,33 . The mobility is obtained as μ = σ∕n c e, where n c 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 10 6 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 40 3 to 600 3 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 200 3 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 elements 18 .
The interpolation method allows us to use extremely dense BZ grids with up to 600 3 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 semiconductors 18,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 k f in the entire fine BZ grid. In the direct computation, one computes the e-d matrix elements Mðk 0 f ; k f Þ between all crystal momentum pairs and thus a number of matrix elements of order ðαN k f Þ 2 . In the interpolation approach, the most time-consuming step is directly computing the N 2 kc 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 onlỹ 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 $ ðαN k f Þ 2 =N 2 kc . The typical value of N k f is around 10 3 -10 4 N kc . For our silicon calculations, the value of α is Fig. 3 Carrier relaxation times and the hole mobility for vacancy defects in silicon. 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. The data are for a neutral vacancy in silicon, with e-d matrix elements computed for the four lowest valence bands and along high-symmetry BZ lines L-Γ-X-K-Γ shown in Fig. 2.
I.-T. Lu et al. 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 kpoints on the Fermi surface of copper, using interpolated e-d matrix elements obtained from a moderate size (8 × 8 × 8) coarse k c -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 method 9,10 . These results show that our firstprinciples approach can access microscopic details of the e-d scattering processes. 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 results 34,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 result [36][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 exists 41 for RRR in metals. As at room temperature the resistivity is usually phonon-limited, combined with e-ph calculations 25,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 interactions 25,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 formalism 20,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 10 8 -10 9 kpoints. 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.

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 pseudopotentials 45 with the QUANTUM ESPRESSO code 46 . 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 Å Fig. 4 Relaxation times and the defect-limited resistivity for vacancy defects in copper. 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 . 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 grid 47 for the primitive cells of both materials to converge the charge density and total energy, and interpolate their band structures using maximally localized WFs 29 with the Wannier90 code 32,33 . Coarse k c -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 k f -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 sp 3 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 ( 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 300 3 k f -points, whereas the interpolated matrix elements are calculated on the same dense grid starting from a coarse grid with 10 3 k c -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 mobility 25 ; 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 300 3 points for the RTs; for the mobility, we use a 1 meV broadening and e-d matrix elements interpolated from a 10 3 coarse BZ grid. In copper, we compute the RTs and resistivity on a fine BZ grid with 240 3 points, using a 1 meV broadening and e-d matrix elements interpolated from a coarse 8 3 BZ grid. The equilibrium vacancy concentration at temperature T in copper is estimated using 48 where ΔH v and ΔS v are the vacancy formation enthalpy and entropy, respectively, and k B the Boltzmann constant. In copper, ΔS v = 3.0 k B and ΔH v = 1.19 eV 48 .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.