Electron-phonon interaction using Wannier functions

We introduce a technique based on the spatial localization of electron and phonon Wannier functions to perform first-principles calculations of the electron-phonon interaction with an ultradense sampling of the Brillouin zone. After developing the basic theory, we describe the practical implementation within a density-functional framework. The proposed method is illustrated by considering a virtual crystal model of boron-doped diamond. For this test case, we first discuss the spatial localization of the electron-phonon matrix element in the Wannier representation. Then, we assess the accuracy of the Wannier-Fourier interpolation in momentum space. Finally, we study the convergence of the electron-phonon self-energies with the sampling of the Brillouin zone by calculating the electron and phonon linewidths, the Eliashberg spectral function, and the mass enhancement parameter of B-doped diamond. We show that more than ${10}^5$ points in the irreducible wedge of the Brillouin zone are needed to achieve convergence.

Figure 1

(Color online) First-order e-ph diagrams considered in this work (red). Left: the self-energy of a phonon with momentum $\mathbf{q}$ (black wiggly line) is given by the Fermion loop containing two electron lines (red lines) connected by the bare e-ph vertices (red disks) (Refs. 1, 12). Right: the self-energy of an electron with momentum $\mathbf{k}$ (black straight line) is given by the loop with one electron line (straight) and one phonon line (wiggly), connected by the bare e-ph vertices. Equations (1, 2) were obtained from these diagrams.

Figure 2

(Color online) Simplified scheme of the electron and phonon Wannier functions entering the three-point e-ph matrix element [Eq. (23)]. The square lattice indicates the unit cells of the crystal, the red lines the electron Wannier functions $\mid {\mathbf{0}}_e⟩$ and $\mid {\mathbf{R}}_e⟩$, and the blue line the phonon perturbation in the Wannier representation ${\partial }_{{\mathbf{R}}_p}V$. Whenever two of these functions are centered on distant unit cells, the e-ph matrix element in the Wannier representation $⟨{\mathbf{0}}_e\mid {\partial }_{{\mathbf{R}}_p}V\mid {\mathbf{R}}_e⟩$ vanishes.

Figure 3

(Color online) Comparison between the electronic band structure of pristine diamond (lines) and of B-doped diamond within a virtual crystal approximation (circles). The Fermi level of doped diamond with a B concentration of 1.85% is located $0.57\mathrm{eV}$ below the valence band top at $\Gamma$ and is indicated by a black solid line.

Figure 4

(Color online) Comparison between the phonon dispersions of pristine diamond (solid lines) and B-doped diamond within a virtual crystal approximation (dashed lines). The disks correspond to inelastic neutron scattering data from Ref. 79.

Figure 5

(Color online) Spatial decay of the electronic Hamiltonian in the Wannier representation $H_{{\mathbf{R}}_e,{\mathbf{R}}_e^{\prime }}^{\mathrm{el}}=⟨m{\mathbf{R}}_e^{\prime }\mid {\widehat{H}}^{\mathrm{el}}\mid n{\mathbf{R}}_e⟩$ [Eq. (26)] as a function of $R=\mid {\mathbf{R}}_e-{\mathbf{R}}_e^{\prime }\mid$. The data points correspond to the largest value taken over the Wannier functions indices and over the unit cells ${\mathbf{R}}_e$, ${\mathbf{R}}_e^{\prime }$ located at the same distance $R$: $\parallel H\left(R\right)\parallel ={\max }_{mn,\mid {\mathbf{R}}_e-{\mathbf{R}}_e^{\prime }\mid =R}\mid ⟨m{\mathbf{R}}_e^{\prime }\mid {\widehat{H}}^{\mathrm{el}}\mid n{\mathbf{R}}_e⟩\mid$. The data are normalized to their largest value. The inset shows the same quantity on a logarithmic scale $\left({\log }_{10}\right)$.

Figure 6

(Color online) Spatial decay of the dynamical matrix in the Wannier representation $D_{{\mathbf{R}}_p,{\mathbf{R}}_p^{\prime }}^{\mathrm{ph}}=⟨{\mathbf{R}}_p\mid {\widehat{D}}^{\mathrm{ph}}\mid {\mathbf{R}}_p^{\prime }⟩$ [Eq. (29)] as a function of the distance $R=\mid {\mathbf{R}}_p-{\mathbf{R}}_p^{\prime }\mid$. The data points correspond to the largest value taken over the ions in the unit cell, the Cartesian directions, and the unit cells ${\mathbf{R}}_p$, ${\mathbf{R}}_p^{\prime }$ located at the same distance $R$: $\parallel D\left(R\right)\parallel ={\max }_{\kappa {\kappa }^{\prime }\alpha {\alpha }^{\prime },\mid {\mathbf{R}}_p-{\mathbf{R}}_p^{\prime }\mid =R}\mid ⟨{\kappa }^{\prime }{\alpha }^{\prime }{\mathbf{R}}_p\mid {\widehat{D}}^{\mathrm{ph}}\mid \kappa \alpha {\mathbf{R}}_p^{\prime }⟩\mid$. The data are normalized to their largest value. The inset shows the same quantity on a logarithmic scale $\left({\log }_{10}\right)$.

Figure 7

(Color online) Spatial decay of the e-ph vertex in the joint electron-phonon Wannier representation $g_{mn,\nu }({\mathbf{R}}_e,{\mathbf{R}}_p)=⟨m{\mathbf{0}}_e\mid {\partial }_{\kappa \alpha ,{\mathbf{R}}_p}V\mid n{\mathbf{R}}_e⟩$ [Eq. (23)] as a function of ${\mathbf{R}}_p$ and ${\mathbf{R}}_e$: (a) the limiting case $g({\mathbf{R}}_e,{\mathbf{R}}_p=\mathbf{0})$ and (b) the limiting case $g({\mathbf{R}}_e=\mathbf{0},{\mathbf{R}}_p)$. The data points correspond to the largest value taken over the Wannier functions indices, the ions in the unit cell, the Cartesian directions, and the unit cells located at the same distance $R=\mid {\mathbf{R}}_e\mid$ [panel (a)] or $R=\mid {\mathbf{R}}_p\mid$ [panel (b)] from the origin of the reference frame: (a) $\parallel g(R,0)\parallel ={\max }_{mn,\kappa \alpha ,\mid {\mathbf{R}}_e\mid =R}\mid g_{mn,\kappa \alpha }({\mathbf{R}}_e,{\mathbf{0}}_p)\mid$ and (b) $\parallel g(0,R)\parallel ={\max }_{mn,\kappa \alpha ,\mid {\mathbf{R}}_p\mid =R}\mid g_{mn,\kappa \alpha }({\mathbf{0}}_e,{\mathbf{R}}_p)\mid$. The data are normalized to their largest value. The insets show the same quantities on a logarithmic scale $\left({\log }_{10}\right)$. When two Wannier functions are located on a C–C bond crossing a cell boundary, identical e-ph matrix elements appear in adjacent unit cells, resulting in the steplike behavior seen in (b).

Figure 8

(Color online) Comparison of the e-ph matrix elements $g_{mn,\nu }^{\mathrm{SE}}(\mathbf{k},\mathbf{q})$ [Eq. (5)] obtained by direct first-principles calculations [panel (c), disks] and those computed with the Wannier-Fourier interpolation method discussed in Sec. 4 [panel (c), lines]. The interpolated matrix elements are computed starting from an initial $4^3$ Brillouin-zone grid (dotted line), a $6^3$ grid (dashed line), or a $8^3$ grid (solid line). For illustration, we fixed the initial electronic state $\mid n\mathbf{k}⟩$ for the valence band top at $\Gamma$ $(\mathbf{k}=\mathbf{0})$; we let the final electronic state $\mid m\mathbf{k}+\mathbf{q}⟩$ span the ${\Lambda }_3$, ${\Delta }_5$, and ${\Sigma }_2$ bands as shown in panel (a) (dashed line), and we take the phonon perturbation corresponding to the highest optical branches as shown in panel (b) (dashed line).

Figure 9

(Color online) Calculated electron linewidths for the ${\Lambda }_3$ and ${\Delta }_5$ bands of B-doped diamond [Fig. 8a]. Plots (a)–(c) on the left were obtained by using ${10}^3$ irreducible $\mathbf{q}$ points in the momentum integration of Eq. (3), while plots (d)–(f) on the right were obtained with ${50}^3$ irreducible points in the Brillouin zone. We report the results for three broadening parameters $\eta$: $10\mathrm{meV}$ [panels (a) and (d)], $50\mathrm{meV}$ [panels (b) and (e)], and $100\mathrm{meV}$ [panels (c) and (f)]. The curves are cut off at half the Brillouin-zone size for clarity.

Figure 10

(Color online) Calculated phonon linewidths for the highest optical mode of B-doped diamond [Fig. 8b]. Plots (a)–(c) on the left were obtained with ${10}^3$ irreducible $\mathbf{q}$ points in the momentum integration of Eq. (4), while plots (d)–(f) on the right were obtained with ${50}^3$ irreducible points. The results for three broadening parameters 10, 50, and $100\mathrm{meV}$ are shown in panels (a) and (d), panels (b) and (e), and panels (c) and (f), respectively. Note the different vertical scales in the plots corresponding to different broadening parameters.

Figure 11

(Color online) Eliashberg function ${\alpha }^{2}F\left(\omega \right)$ calculated for B-doped diamond in the virtual crystal approximation with (a) ${10}^3\mathbf{k}$ and $\mathbf{q}$ points in the irreducible part of the Brillouin zone and (b) ${50}^3$ points. For each case, we report the results corresponding to the smearing parameters $10\mathrm{meV}$ (dotted lines), $50\mathrm{meV}$ (dashed lines), and $100\mathrm{meV}$ (solid lines). The Dirac delta function in Eq. (7) was replaced by a Gaussian with a standard deviation of $0.5\mathrm{meV}$. Note the different energy scales and scaling factors for the acoustic and the optical frequency ranges.