Theory

Notation and conventions

Jon Lafuente-Bartolomé 07/17/2026

The main definitions and conventions used throughout this page are collected below. We follow the notation of Ref. [1].

Crystal and reciprocal-space notation

The primitive lattice vectors are denoted by \(\mathbf{a}_i\), and the reciprocal lattice vectors \(\mathbf{b}_j\) satisfy \(\mathbf{a}_i\cdot\mathbf{b}_j=2\pi\delta_{ij}\). We describe the crystal using Born-von Kármán boundary conditions in a supercell containing \(N_p\) primitive cells. The lattice vector of cell \(p\) is denoted as \(\mathbf{R}_p\), and the equilibrium position of atom \(\kappa\) in that cell is \(\boldsymbol{\tau}_{\kappa p}^{0} =\mathbf{R}_p+\boldsymbol{\tau}_{\kappa}\). Its component along the Cartesian direction \(\alpha\) is \(\tau_{\kappa\alpha p}^{0}=R_{p\alpha}+\tau_{\kappa\alpha}\). Coordinates without the superscript \(0\) denote general, possibly displaced, atomic positions. Electronic wavevectors are denoted by \(\mathbf{k}\), phonon wavevectors by \(\mathbf{q}\), and reciprocal-lattice vectors by \(\mathbf{G}\).

Electronic states

Electronic bands are labelled by \(n\) and \(m\). A Kohn-Sham eigenstate is written in Bloch form as

\[\psi_{n\mathbf{k}}(\mathbf{r}) = N_p^{-1/2} u_{n\mathbf{k}}(\mathbf{r}) e^{i\mathbf{k}\cdot\mathbf{r}},\]

where \(u_{n\mathbf{k}}\) is lattice-periodic. The Bloch state \(\psi_{n\mathbf{k}}\) is normalized over the Born-von Kármán supercell, while \(u_{n\mathbf{k}}\) is normalized over one primitive cell. Its Kohn-Sham eigenvalue is denoted as \(\varepsilon_{n\mathbf{k}}\). Spin indices are suppressed unless they are explicitly needed.

At temperature \(T\), the Fermi-Dirac occupation factor of this state is

\[f_{n\mathbf{k}} = \frac{1}{e^{(\varepsilon_{n\mathbf{k}}-\mu)/(k_{\mathrm B}T)}+1},\]

where \(\mu\) is the chemical potential and \(k_{\mathrm B}\) is the Boltzmann constant.

Phonons and lattice dynamics

Let \(U(\{\boldsymbol{\tau}_{\kappa p}\})\) be the total potential energy of the electrons and nuclei, with the electrons in their ground state and the nuclei clamped at fixed coordinates. The interatomic force constants are the second derivatives of \(U\) with respect to the nuclear coordinates,

\[C_{\kappa\alpha p,\kappa'\alpha' p'} = \left. \frac{\partial^2 U} {\partial\tau_{\kappa\alpha p}\, \partial\tau_{\kappa'\alpha' p'}} \right|_0.\]

The derivatives are evaluated at the equilibrium atomic configuration. Using translational invariance and taking cell \(0\) as the reference cell, the mass-weighted dynamical matrix is their Fourier transform,

\[D_{\kappa\alpha,\kappa'\alpha'}(\mathbf{q}) = \frac{1}{\sqrt{M_\kappa M_{\kappa'}}} \sum_p e^{i\mathbf{q}\cdot\mathbf{R}_p} C_{\kappa\alpha 0,\kappa'\alpha' p},\]

where \(M_\kappa\) is the mass of atom \(\kappa\). The dynamical matrix can be computed using DFPT [2]. The phonon frequencies and polarization vectors are obtained by diagonalizing the dynamical matrix:

\[\sum_{\kappa'\alpha'} D_{\kappa\alpha,\kappa'\alpha'}(\mathbf{q}) e_{\kappa'\alpha',\nu}(\mathbf{q}) = \omega_{\mathbf{q}\nu}^{2} e_{\kappa\alpha,\nu}(\mathbf{q}).\]

A phonon mode is labelled by its wavevector \(\mathbf{q}\) and branch index \(\nu\); its frequency and polarization vector are \(\omega_{\mathbf{q}\nu}\) and \(e_{\kappa\alpha,\nu}(\mathbf{q})\), respectively. The polarization vectors are normalized according to

\[\sum_{\kappa\alpha} e^{*}_{\kappa\alpha,\nu}(\mathbf{q}) e_{\kappa\alpha,\nu'}(\mathbf{q}) = \delta_{\nu\nu'}.\]

Their phases are chosen such that

\[e_{\kappa\alpha,\nu}(-\mathbf{q}) =e^{*}_{\kappa\alpha,\nu}(\mathbf{q}).\]

We define the zero-point displacement amplitude as

\[l_{\mathbf{q}\nu} = \left(\frac{\hbar}{2M_0\omega_{\mathbf{q}\nu}}\right)^{1/2},\]

where \(\hbar\) is the reduced Planck constant and \(M_0\) is an arbitrary reference mass, conventionally the proton mass. The three zero-frequency translational modes at \(\mathbf{q}=0\) are omitted from expressions containing \(l_{\mathbf{q}\nu}\).

The Bose-Einstein occupation factor of a mode \((\mathbf{q}, \nu)\) is

\[n_{\mathbf{q}\nu} = \frac{1}{e^{\hbar\omega_{\mathbf{q}\nu}/(k_{\mathrm B}T)}-1}.\]

Electron-phonon matrix elements

Jon Lafuente-Bartolomé 07/17/2026

The electron-phonon matrix element is defined as

\[g_{mn\nu}(\mathbf{k},\mathbf{q}) = \left\langle u_{m\mathbf{k}+\mathbf{q}} \middle| \Delta_{\mathbf{q}\nu}v^{\mathrm{KS}} \middle| u_{n\mathbf{k}} \right\rangle_{\mathrm{uc}},\]

where the integral is evaluated over one primitive cell. This quantity gives the scattering amplitude for an electron to transition from the state \((n,\mathbf{k})\) to \((m,\mathbf{k}+\mathbf{q})\) through a phonon with wavevector \(\mathbf{q}\) and branch index \(\nu\). For this phonon, the full first-order variation of the self-consistent Kohn-Sham potential is written as

\[\Delta_{\mathbf{q}\nu}V^{\mathrm{KS}}(\mathbf{r}) = e^{i\mathbf{q}\cdot\mathbf{r}} \Delta_{\mathbf{q}\nu}v^{\mathrm{KS}}(\mathbf{r}),\]

where \(\Delta_{\mathbf{q}\nu}v^{\mathrm{KS}}\) is lattice-periodic. It can be computed using density-functional perturbation theory [2] and is defined as

\[\Delta_{\mathbf{q}\nu}v^{\mathrm{KS}} = l_{\mathbf{q}\nu} \sum_{\kappa\alpha} \left(\frac{M_0}{M_\kappa}\right)^{1/2} e_{\kappa\alpha,\nu}(\mathbf{q}) \partial_{\kappa\alpha,\mathbf{q}}v^{\mathrm{KS}},\]

with

\[\partial_{\kappa\alpha,\mathbf{q}}v^{\mathrm{KS}}(\mathbf{r}) = \sum_p e^{-i\mathbf{q}\cdot(\mathbf{r}-\mathbf{R}_p)} \left.\frac{\partial V^{\mathrm{KS}} (\mathbf{r};\{\boldsymbol{\tau}_{\kappa p}\})} {\partial\tau_{\kappa\alpha p}}\right|_0.\]

The derivative is evaluated at the equilibrium atomic configuration.

Wannier-Fourier interpolation

Calculations of observables related to electron-phonon coupling typically require matrix elements on much denser grids than can be computed directly with DFPT. In EPW they can be efficiently calculated using the Wannier-Fourier interpolation technique, as described in Ref. [3]. The Wannier functions are defined as

\[w_{mp}(\mathbf{r}) = \frac{1}{N_p}\sum_{n\mathbf{k}} e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{R}_p)} U_{nm\mathbf{k}}u_{n\mathbf{k}}(\mathbf{r}),\]

where \(U_{nm\mathbf{k}}\) is a unitary transformation within the chosen band manifold. The maximally localized Wannier functions are obtained by choosing these matrices so as to minimize their spatial spread [4]. For entangled bands, a smooth subspace is first selected using a disentanglement procedure [5], and \(U_{nm\mathbf{k}}\) then acts within this subspace. The inverse transformation is

\[u_{n\mathbf{k}}(\mathbf{r}) = \sum_{mp}e^{-i\mathbf{k}\cdot(\mathbf{r}-\mathbf{R}_p)} U^{\dagger}_{mn\mathbf{k}}w_{mp}(\mathbf{r}).\]

The electron-phonon matrix elements in the Wannier representation are defined as

\[g_{mn\kappa\alpha}(\mathbf{R}_p,\mathbf{R}_{p'}) = \left\langle w_{m0}(\mathbf{r})\middle| \frac{\partial V^{\mathrm{KS}}}{\partial\tau_{\kappa\alpha}} (\mathbf{r}-\mathbf{R}_{p'}) \middle|w_{n0}(\mathbf{r}-\mathbf{R}_p)\right\rangle_{\mathrm{sc}},\]

where the integral is evaluated over the Born-von Kármán supercell and the index \(0\) denotes the reference cell.

The matrix elements on the fine grids are obtained from their real-space counterparts through

\[\begin{split}\begin{aligned} g_{mn\nu}(\mathbf{k},\mathbf{q}) ={}& \sum_{pp'}e^{i(\mathbf{k}\cdot\mathbf{R}_p +\mathbf{q}\cdot\mathbf{R}_{p'})} \\ &\times\sum_{m'n'\kappa\alpha} U_{mm',\mathbf{k}+\mathbf{q}} g_{m'n'\kappa\alpha}(\mathbf{R}_p,\mathbf{R}_{p'}) U^{\dagger}_{n'n\mathbf{k}} u_{\kappa\alpha,\mathbf{q}\nu}, \end{aligned}\end{split}\]

where \(u_{\kappa\alpha,\mathbf{q}\nu}\) is the mass-scaled phonon eigenvector,

\[u_{\kappa\alpha,\mathbf{q}\nu} = \left(\frac{\hbar}{2M_\kappa\omega_{\mathbf{q}\nu}}\right)^{1/2} e_{\kappa\alpha,\nu}(\mathbf{q}).\]

Conversely, the coarse-grid matrix elements are transformed to real space using

\[\begin{split}\begin{aligned} g_{mn\kappa\alpha}(\mathbf{R}_p,\mathbf{R}_{p'}) ={}& \frac{1}{N_pN_{p'}} \sum_{\mathbf{k}\mathbf{q}} e^{-i(\mathbf{k}\cdot\mathbf{R}_p +\mathbf{q}\cdot\mathbf{R}_{p'})} \\ &\times\sum_{m'n'\nu} u^{-1}_{\kappa\alpha,\mathbf{q}\nu} U^{\dagger}_{mm',\mathbf{k}+\mathbf{q}} g_{m'n'\nu}(\mathbf{k},\mathbf{q}) U_{n'n\mathbf{k}}, \end{aligned}\end{split}\]

where \(u^{-1}_{\kappa\alpha,\mathbf{q}\nu}\) denotes the inverse phonon-mode transformation,

\[u^{-1}_{\kappa\alpha,\mathbf{q}\nu} = \left(\frac{2M_\kappa\omega_{\mathbf{q}\nu}}{\hbar}\right)^{1/2} e^{*}_{\kappa\alpha,\nu}(\mathbf{q}).\]

The symbols \(N_p\) and \(N_{p'}\) allow for different Born-von Kármán supercells for the electron and phonon grids. If \(g_{mn\kappa\alpha}(\mathbf{R}_p,\mathbf{R}_{p'})\) decays rapidly with \(|\mathbf{R}_p|\) and \(|\mathbf{R}_{p'}|\), the coarse-grid calculation contains all the information needed to reconstruct \(g_{mn\nu}(\mathbf{k},\mathbf{q})\) on arbitrarily fine grids. In practice, the Wannier-Fourier interpolation in EPW proceeds as follows:

1. Compute the electron-phonon matrix elements on coarse \(\mathbf{k}\)- and \(\mathbf{q}\)-point grids, the Wannier-Bloch rotation matrices on the coarse electronic grid, and the phonon frequencies and eigenvectors on the coarse phonon grid.

2. Transform the matrix elements from the coarse reciprocal-space grids to the real-space Wannier representation.

3. Exploit their spatial localization by neglecting matrix elements outside the corresponding Wigner-Seitz supercells.

4. Transform back to reciprocal space on the fine grids. The electronic Hamiltonian and phonon dynamical matrix are interpolated alongside the electron-phonon matrix elements to obtain the required Kohn-Sham energies, phonon frequencies and eigenvectors, and Wannier-Bloch rotation matrices.

Long-range contributions

In semiconductors and insulators, atomic displacements can generate long-range electrostatic potentials. The resulting matrix elements are not localized in the Wannier representation and must therefore be treated separately. The main idea is to decompose the electron-phonon matrix elements into

\[g_{mn\nu}(\mathbf{k},\mathbf{q}) = g^{\mathrm{S}}_{mn\nu}(\mathbf{k},\mathbf{q}) + g^{\mathrm{L,D}}_{mn\nu}(\mathbf{k},\mathbf{q}) + g^{\mathrm{L,Q}}_{mn\nu}(\mathbf{k},\mathbf{q}),\]

where \(g^{\mathrm{S}}\) is the short-range component, and \(g^{\mathrm{L,D}}\) and \(g^{\mathrm{L,Q}}\) are the long-range dipole and quadrupole terms, respectively. For three-dimensional bulk crystals, the dipole term contains the so-called Fröhlich interaction, which diverges as \(|\mathbf{q}|^{-1}\) for long-wavelength longitudinal-optical phonons in polar materials, and is determined by the Born effective charges and high-frequency dielectric tensor [6]. The quadrupole term is determined by the dynamical quadrupole tensors. It remains finite but can be direction-dependent as \(\mathbf{q}\rightarrow 0\), as discussed in Refs. [7] and [8]. The corresponding formulations for two-dimensional materials, where different electrostatic screening modifies the long-wavelength behavior, are presented in Refs. [9] and [10].

In practice, EPW treats the long-range contributions as follows:

1. Compute the complete electron-phonon matrix elements with DFPT on the coarse grids.

2. Evaluate and subtract the analytic dipole and quadrupole terms, leaving a localized short-range component.

3. Wannier-Fourier interpolate the short-range component from the coarse to the fine grids.

4. Evaluate the long-range terms on the fine grids and add them back to the interpolated short-range component. The corresponding long-range contributions to the phonon dynamical matrix are treated consistently.

The explicit expressions and implementation details are given in the Methods section of Ref. [11].

Self-energies and coupling strength

Jon Lafuente-Bartolomé 07/17/2026

Phonon linewidths and coupling strength

Neglecting the coupling between different phonon branches, the complex poles of the dressed phonon propagator satisfies

\[\widetilde{\Omega}_{\mathbf{q}\nu}^{2} = \omega_{\mathbf{q}\nu}^{2} + 2\omega_{\mathbf{q}\nu} \Pi^{\mathrm{NA}}_{\mathbf{q}\nu\nu} (\widetilde{\Omega}_{\mathbf{q}\nu}), \qquad \widetilde{\Omega}_{\mathbf{q}\nu} = \Omega_{\mathbf{q}\nu}-i\gamma_{\mathbf{q}\nu},\]

where \(\Pi^{\mathrm{NA}}\) is the nonadiabatic phonon self-energy, \(\Omega_{\mathbf{q}\nu}\) is the renormalized phonon frequency, and \(\gamma_{\mathbf{q}\nu}\) is the half-width. If the frequency shift obeys \(|\Omega_{\mathbf{q}\nu}-\omega_{\mathbf{q}\nu}|\ll\omega_{\mathbf{q}\nu}\) and \(\gamma_{\mathbf{q}\nu}\ll\omega_{\mathbf{q}\nu}\), linearizing the pole equation and evaluating the self-energy on shell gives

\[\Omega_{\mathbf{q}\nu} \simeq \omega_{\mathbf{q}\nu} +\operatorname{Re}\Pi^{\mathrm{NA}}_{\mathbf{q}\nu\nu} (\omega_{\mathbf{q}\nu}), \qquad \gamma_{\mathbf{q}\nu} \simeq -\operatorname{Im}\Pi^{\mathrm{NA}}_{\mathbf{q}\nu\nu} (\omega_{\mathbf{q}\nu}).\]

The real and imaginary parts of the phonon self-energy therefore describe the phonon frequency shift and broadening, respectively [1]. For a spin-degenerate metal, the low-temperature double-delta approximation introduced in Ref. [12] (and used in the interpolation tutorial) gives

\[\gamma_{\mathbf{q}\nu} = 2\pi\omega_{\mathbf{q}\nu} \sum_{mn}\int_{\mathrm{BZ}}\frac{d\mathbf{k}}{\Omega_{\mathrm{BZ}}} \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 \delta(\varepsilon_{n\mathbf{k}}-\varepsilon_{\mathrm F}) \delta(\varepsilon_{m\mathbf{k}+\mathbf{q}}-\varepsilon_{\mathrm F}),\]

where \(\Omega_{\mathrm{BZ}}\) is the Brillouin-zone volume and \(\varepsilon_{\mathrm F}\) is the Fermi energy. Here \(\gamma_{\mathbf{q}\nu}\) has units of angular frequency; EPW reports \(\hbar\gamma_{\mathbf{q}\nu}\) in energy units. The phonon lifetime is \(\tau^{\mathrm{ph}}_{\mathbf{q}\nu}=(2\gamma_{\mathbf{q}\nu})^{-1}\).

The mode-resolved electron-phonon coupling strength is [13]

\[\lambda_{\mathbf{q}\nu} = \frac{\gamma_{\mathbf{q}\nu}} {\pi\hbar N_{\mathrm F}\omega_{\mathbf{q}\nu}^{2}},\]

where \(N_{\mathrm F}\) is the electronic density of states at the Fermi energy. The isotropic Eliashberg spectral function is

\[\begin{split}\begin{aligned} \alpha^2F(\omega) ={}& \frac{1}{N_{\mathrm F}} \sum_{mn\nu}\int_{\mathrm{BZ}} \frac{d\mathbf{k}\,d\mathbf{q}}{\Omega_{\mathrm{BZ}}^2} \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 \\ &\times \delta(\varepsilon_{n\mathbf{k}}-\varepsilon_{\mathrm F}) \delta(\varepsilon_{m\mathbf{k}+\mathbf{q}}-\varepsilon_{\mathrm F}) \delta(\hbar\omega-\hbar\omega_{\mathbf{q}\nu}), \end{aligned}\end{split}\]

and the total isotropic coupling strength can be obtained from

\[\lambda = \sum_\nu\int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \lambda_{\mathbf{q}\nu} = 2\int_0^\infty \frac{\alpha^2F(\omega)}{\omega}\,d\omega.\]

Electron self-energy

Within the band-diagonal approximation, the finite-temperature Fan-Migdal self-energy is

\[\begin{split}\begin{aligned} \Sigma^{\mathrm{FM}}_{nn\mathbf{k}}(\omega,T) ={}& \frac{1}{\hbar} \sum_{m\nu}\int_{\mathrm{BZ}} \frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 \\ &\times\left[ \frac{1-f_{m\mathbf{k}+\mathbf{q}}+n_{\mathbf{q}\nu}} {\omega-\varepsilon_{m\mathbf{k}+\mathbf{q}}/\hbar -\omega_{\mathbf{q}\nu}+i\eta} + \frac{f_{m\mathbf{k}+\mathbf{q}}+n_{\mathbf{q}\nu}} {\omega-\varepsilon_{m\mathbf{k}+\mathbf{q}}/\hbar +\omega_{\mathbf{q}\nu}+i\eta} \right], \end{aligned}\end{split}\]

where \(\eta\) is a positive infinitesimal.

The static Debye-Waller contribution is

\[\Sigma^{\mathrm{DW}}_{nn\mathbf{k}}(T) = \sum_\nu\int_{\mathrm{BZ}} \frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} g^{\mathrm{DW}}_{nn\nu\nu} (\mathbf{k},\mathbf{q},-\mathbf{q}) \left(2n_{\mathbf{q}\nu}+1\right),\]

where \(g^{\mathrm{DW}}\) is the second-order electron-phonon matrix element; its detailed definition is given in Ref. [1].

Spectral functions and quasiparticles

In terms of the retarded Green’s function, the electron spectral function is [14]

\[A(\mathbf{k},\omega) = -\frac{1}{\pi}\sum_n \operatorname{Im}G^{\mathrm{ret}}_{nn\mathbf{k}}(\omega).\]

When the self-energy is diagonal in the band indices, this becomes

\[A(\mathbf{k},\omega) = \sum_n \frac{-\pi^{-1}\operatorname{Im} \Sigma_{nn\mathbf{k}}(\omega,T)} {\left[\hbar\omega-\varepsilon_{n\mathbf{k}} -\operatorname{Re}\Sigma_{nn\mathbf{k}}(\omega,T)\right]^2 +\left[\operatorname{Im} \Sigma_{nn\mathbf{k}}(\omega,T)\right]^2}.\]

Its maxima and widths describe the renormalized dispersion and excitation lifetimes, while additional structures can appear as kinks or satellites.

If the spectral function exhibits a well-defined quasiparticle peak, the Green’s function near the pole can be approximated as

\[G^{\mathrm{ret}}_{nn\mathbf{k}}(\omega) \simeq \frac{Z_{n\mathbf{k}}} {\hbar\omega-E_{n\mathbf{k}}(T)+i\Gamma_{n\mathbf{k}}}.\]

Here \(E_{n\mathbf{k}}\) is the quasiparticle energy and \(Z_{n\mathbf{k}}\) is the quasiparticle weight,

\[Z_{n\mathbf{k}} = \left[ 1-\frac{1}{\hbar} \left. \frac{\partial\operatorname{Re}\Sigma_{nn\mathbf{k}} (\omega,T)}{\partial\omega} \right|_{\omega=E_{n\mathbf{k}}/\hbar} \right]^{-1},\]

where \(\Gamma_{n\mathbf{k}}\) and \(\tau_{n\mathbf{k}}\) are the quasiparticle half-width in energy units and lifetime, respectively:

\[\Gamma_{n\mathbf{k}} = Z_{n\mathbf{k}} \left|\operatorname{Im}\Sigma_{nn\mathbf{k}} \left(E_{n\mathbf{k}}/\hbar,T\right)\right|, \qquad \tau_{n\mathbf{k}}=\frac{\hbar}{2\Gamma_{n\mathbf{k}}}.\]

The electron-state-resolved electron-phonon coupling strength is

\[\lambda_{n\mathbf{k}} = -\frac{1}{\hbar} \left. \frac{\partial\operatorname{Re}\Sigma^{\mathrm{FM}}_{nn\mathbf{k}} (\omega,T)}{\partial\omega} \right|_{\omega=E_{n\mathbf{k}}/\hbar},\]

and the quasiparticle weight can then be written as \(Z_{n\mathbf{k}}=(1+\lambda_{n\mathbf{k}})^{-1}\).

Band renormalization

Neglecting polaron self-trapping effects (see Sec. Ab initio polaron equations) and band off-diagonal self-energy matrix elements, the quasiparticle energy satisfies

\[E_{n\mathbf{k}}(T) = \varepsilon_{n\mathbf{k}} + \operatorname{Re}\Sigma^{\mathrm{FM}}_{nn\mathbf{k}} \left(E_{n\mathbf{k}}/\hbar,T\right) + \Sigma^{\mathrm{DW}}_{nn\mathbf{k}}(T),\]

which leads to the theory of temperature-dependent band renormalization developed by Allen and Heine [15]. In the on-shell approximation, \(E_{n\mathbf{k}}\) on the right-hand side is replaced by \(\varepsilon_{n\mathbf{k}}\).

Ab initio polaron equations

Jon Lafuente-Bartolomé 07/17/2026

The ground-state wave function \(\psi(\mathbf{r})\) and atomic displacements \(\Delta\tau_{\kappa\alpha p}=\tau_{\kappa\alpha p}-\tau_{\kappa\alpha p}^{0}\) forming a polaron can be found by minimizing the total DFT energy functional of an excess electron added to a crystal, as described in Refs. [16] and [17]. This yields the coupled equations

(1)\[\hat{H}_{\mathrm{KS}}^{0}\psi(\mathbf{r}) + \sum_{\kappa\alpha p} \left. \frac{\partial V^{\mathrm{KS}}} {\partial\tau_{\kappa\alpha p}} \right|_0 \Delta\tau_{\kappa\alpha p}\psi(\mathbf{r}) = \varepsilon\psi(\mathbf{r}),\]

(2)\[\Delta\tau_{\kappa\alpha p} = -\sum_{\kappa'\alpha'p'} (C^{-1})_{\kappa\alpha p,\kappa'\alpha'p'} \int d\mathbf{r}\, \left. \frac{\partial V^{\mathrm{KS}}} {\partial\tau_{\kappa'\alpha'p'}} \right|_0 |\psi(\mathbf{r})|^2.\]

Here \(\hat{H}_{\mathrm{KS}}^{0}\) is the Kohn-Sham Hamiltonian of the pristine system, and \(\varepsilon\) is the polaron eigenvalue. The integral is performed over a Born-von Kármán supercell containing \(N_p\) primitive cells.

Reciprocal-space formulation

The polaron wave function can be expanded in terms of the single-particle Kohn-Sham states,

(3)\[\psi(\mathbf{r}) = \frac{1}{\sqrt{N_p}} \sum_{n\mathbf{k}}A_{n\mathbf{k}}\psi_{n\mathbf{k}}(\mathbf{r}),\]

and the atomic displacements in terms of the phonon eigenmodes,

(4)\[\Delta\tau_{\kappa\alpha p} = -\frac{2}{N_p}\sum_{\mathbf{q}\nu} B_{\mathbf{q}\nu}^{*} \left(\frac{\hbar} {2M_\kappa\omega_{\mathbf{q}\nu}}\right)^{1/2} e_{\kappa\alpha,\nu}(\mathbf{q}) e^{i\mathbf{q}\cdot\mathbf{R}_p}.\]

The real-space equations can then be transformed into a coupled set of equations for the expansion coefficients:

(5)\[\frac{2}{N_p}\sum_{\mathbf{q}m\nu} B_{\mathbf{q}\nu} g_{mn\nu}^{*}(\mathbf{k},\mathbf{q}) A_{m\mathbf{k}+\mathbf{q}} = (\varepsilon_{n\mathbf{k}}-\varepsilon)A_{n\mathbf{k}},\]

(6)\[B_{\mathbf{q}\nu} = \frac{1}{N_p}\sum_{mn\mathbf{k}} A_{m\mathbf{k}+\mathbf{q}}^{*} \frac{g_{mn\nu}(\mathbf{k},\mathbf{q})} {\hbar\omega_{\mathbf{q}\nu}} A_{n\mathbf{k}}.\]

These expressions involve the Kohn-Sham eigenvalues, phonon frequencies, and the electron-phonon matrix elements. The coupled equations can be written as an effective eigenvalue problem,

(7)\[\sum_{n'\mathbf{k}'} H_{n\mathbf{k},n'\mathbf{k}'} A_{n'\mathbf{k}'} = \varepsilon A_{n\mathbf{k}},\]

where

(8)\[H_{n\mathbf{k},n'\mathbf{k}'} = \delta_{n\mathbf{k},n'\mathbf{k}'} \varepsilon_{n\mathbf{k}} - \frac{2}{N_p}\sum_\nu B_{\mathbf{k}-\mathbf{k}',\nu}^{*} g_{nn'\nu}(\mathbf{k}',\mathbf{k}-\mathbf{k}').\]

In practice, the equation for \(B_{\mathbf{q}\nu}\) and the effective eigenvalue problem are solved iteratively until self-consistency is reached. EPW first obtains the Kohn-Sham eigenvalues, phonon frequencies, and electron-phonon matrix elements on the required fine grids using Wannier-Fourier interpolation, and then performs an iterative self-consistent solution of these equations [11].

For visualization, the atomic displacements are recovered from \(B_{\mathbf{q}\nu}\), while the real-space wave function can be obtained from the maximally localized Wannier functions as

(9)\[\psi(\mathbf{r}) = \sum_{mp}A_m(\mathbf{R}_p) w_{mp}(\mathbf{r}),\]

with

(10)\[A_m(\mathbf{R}_p) = \frac{1}{N_p}\sum_{n\mathbf{k}} e^{i\mathbf{k}\cdot\mathbf{R}_p} U^\dagger_{mn\mathbf{k}}A_{n\mathbf{k}}.\]

Here \(U^\dagger_{mn\mathbf{k}}\) is the Wannier-Bloch rotation matrix introduced in the Wannier-Fourier interpolation section.

Non-diagonal supercells

A Born-von Kármán supercell is defined by an integer, nonsingular transformation matrix \(S\) relating its direct lattice vectors \(\mathbf{a}_{si}\) to the primitive-cell vectors \(\mathbf{a}_{pj}\):

(11)\[\mathbf{a}_{si}=\sum_j S_{ij}\mathbf{a}_{pj}.\]

The corresponding reciprocal lattice vectors transform according to

(12)\[\mathbf{b}_{si} =\sum_j\left[(S^{-1})^{T}\right]_{ij}\mathbf{b}_{pj}.\]

The supercell contains \(N_p=|\det S|\) primitive cells. Off-diagonal elements of \(S\) allow its shape to differ from that of the primitive cell. The compatible electronic and phonon grids are the points of the supercell reciprocal lattice that lie within the primitive-cell Brillouin zone. EPW can generate these grids directly from \(S\), enabling polaron calculations in supercells of general shape. Further details are given in Proc. Natl. Acad. Sci. USA 121, e2318151121 (2024).

Formation energy

The polaron formation energy \(\Delta E_{\mathrm f}\), defined as the energy required to trap a conduction-band state with eigenvalue \(\varepsilon_{\mathrm{CBM}}\) into a localized polaron state, can be written as

(13)\[\Delta E_{\mathrm f} = \frac{1}{N_p}\sum_{n\mathbf{k}} |A_{n\mathbf{k}}|^2 (\varepsilon_{n\mathbf{k}}-\varepsilon_{\mathrm{CBM}}) - \frac{1}{N_p}\sum_{\mathbf{q}\nu} |B_{\mathbf{q}\nu}|^2\hbar\omega_{\mathbf{q}\nu}.\]

The first and second terms on the right-hand side are the electron and phonon parts of the formation energy, respectively [17]. To obtain the formation energy of an isolated polaron, the equations are solved on progressively denser \(\mathbf{k}\)- and \(\mathbf{q}\)-point grids and the results are extrapolated to the infinite-supercell limit. Further practical details are given in the polaron tutorial.

Polaron energy landscapes

For a fixed configuration of atomic displacements \(\{\Delta\tau_{\kappa\alpha p}\}\), the energy can be minimized with respect to the polaron wave function only. The resulting atomic-configuration-dependent formation energy is

(14)\[\Delta\widetilde{E}_{\mathrm f} = \varepsilon-\varepsilon_{\mathrm{CBM}} + \frac{1}{N_p}\sum_{\mathbf{q}\nu} |B_{\mathbf{q}\nu}|^2\hbar\omega_{\mathbf{q}\nu}.\]

The tilde indicates that the energy is minimized with respect to the electronic degrees of freedom, while the forces on the atoms are generally nonzero. If this quantity is also minimized with respect to the atomic displacements, the formation energy \(\Delta E_{\mathrm f}\) is recovered.

To compute polaron energy landscapes in EPW, a displacement configuration is first transformed into the coefficients \(B_{\mathbf{q}\nu}\). The effective Hamiltonian is then diagonalized once, without a self-consistency cycle for the displacements, and the associated formation energy is evaluated. Repeating this procedure for a sequence of configurations maps the adiabatic polaron energy landscape. A practical example is given in the polaron energy-landscape exercise.

Boltzmann transport equation

This section is currently being updated.

Superconductivity

This section is currently being updated.

Optical absorption

This section is currently being updated.

Special-displacement method

This section is currently being updated.

References

[1] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. 89, 015003 (2017), doi:10.1103/RevModPhys.89.015003.

[2] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001), doi:10.1103/RevModPhys.73.515.

[3] F. Giustino, M. L. Cohen, and S. G. Louie, Electron-phonon interaction using Wannier functions, Phys. Rev. B 76, 165108 (2007), doi:10.1103/PhysRevB.76.165108.

[4] N. Marzari and D. Vanderbilt, Maximally localized generalized Wannier functions for composite energy bands, Phys. Rev. B 56, 12847 (1997), doi:10.1103/PhysRevB.56.12847.

[5] I. Souza, N. Marzari, and D. Vanderbilt, Maximally localized Wannier functions for entangled energy bands, Phys. Rev. B 65, 035109 (2001), doi:10.1103/PhysRevB.65.035109.

[6] C. Verdi and F. Giustino, Fröhlich electron-phonon vertex from first principles, Phys. Rev. Lett. 115, 176401 (2015), doi:10.1103/PhysRevLett.115.176401.

[7] G. Brunin et al., Phonon-limited electron mobility in Si, GaAs, and GaP with exact treatment of dynamical quadrupoles, Phys. Rev. B 102, 094308 (2020), doi:10.1103/PhysRevB.102.094308.

[8] J. Park, J.-J. Zhou, V. A. Jhalani, C. E. Dreyer, and M. Bernardi, Long-range quadrupole electron-phonon interaction from first principles, Phys. Rev. B 102, 125203 (2020), doi:10.1103/PhysRevB.102.125203.

[9] S. Poncé, M. Royo, M. Stengel, N. Marzari, and M. Gibertini, Long-range electrostatic contribution to electron-phonon couplings and mobilities of two-dimensional and bulk materials, Phys. Rev. B 107, 155424 (2023), doi:10.1103/PhysRevB.107.155424.

[10] W. H. Sio and F. Giustino, Unified ab initio description of Fröhlich electron-phonon interactions in two-dimensional and three-dimensional materials, Phys. Rev. B 105, 115414 (2022), doi:10.1103/PhysRevB.105.115414.

[11] H. Lee et al., Electron-phonon physics from first principles using the EPW code, npj Comput. Mater. 9, 156 (2023), doi:10.1038/s41524-023-01107-3.

[12] P. B. Allen, Neutron spectroscopy of superconductors, Phys. Rev. B 6, 2577 (1972), doi:10.1103/PhysRevB.6.2577.

[13] G. Grimvall, The Electron-Phonon Interaction in Metals (North-Holland, Amsterdam, 1981).

[14] G. D. Mahan, Many-Particle Physics, 3rd ed. (Kluwer Academic/Plenum, New York, 2000), doi:10.1007/978-1-4757-5714-9.

[15] P. B. Allen and V. Heine, Theory of the temperature dependence of electronic band structures, J. Phys. C: Solid State Phys. 9, 2305 (1976), doi:10.1088/0022-3719/9/12/013.

[16] W. H. Sio, C. Verdi, S. Poncé, and F. Giustino, Polarons from first principles, without supercells, Phys. Rev. Lett. 122, 246403 (2019), doi:10.1103/PhysRevLett.122.246403.

[17] W. H. Sio, C. Verdi, S. Poncé, and F. Giustino, Ab initio theory of polarons: Formalism and applications, Phys. Rev. B 99, 235139 (2019), doi:10.1103/PhysRevB.99.235139.