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

(1)\[\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

(2)\[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,

(3)\[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,

(4)\[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:

(5)\[\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

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

Their phases are chosen such that

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

We define the zero-point displacement amplitude as

(8)\[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

(9)\[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

(10)\[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

(11)\[\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

(12)\[\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

(13)\[\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

(14)\[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

(15)\[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

(16)\[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

(17)\[\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,

(18)\[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

(19)\[\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,

(20)\[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

(21)\[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

(22)\[\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

(23)\[\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

(24)\[\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]

(25)\[\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

(26)\[\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

(27)\[\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

(28)\[\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

(29)\[\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]

(30)\[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

(31)\[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

(32)\[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,

(33)\[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:

(34)\[\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

(35)\[\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

(36)\[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

(37)\[\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}),\]

(38)\[\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,

(39)\[\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,

(40)\[\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:

(41)\[\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}},\]

(42)\[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,

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

where

(44)\[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

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

with

(46)\[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}\):

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

The corresponding reciprocal lattice vectors transform according to

(48)\[\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

(49)\[\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

(50)\[\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.

Superconductivity

Shashi B. Mishra 08/05/2026

Phonon-mediated superconductivity is described in EPW within the Migdal-Eliashberg (ME) formalism [18-21], which extends the weak-coupling Bardeen-Cooper-Schrieffer theory by treating the retarded, frequency-dependent electron-phonon interaction on the same footing as the repulsive electron-electron interaction. The starting point is the electron self-energy in the superconducting state, written as the sum of a Fan-Migdal electron-phonon term and a Coulomb term. Applying the Eliashberg procedure [19] to the Dyson equation then reduces the problem to a small set of coupled non-linear equations on the imaginary frequency axis. The practical aspects of these calculations are covered in the superconductivity tutorial.

Nambu-Gor’kov formalism

The superconducting state is described with the two-component Nambu spinor [22]

(51)\[\begin{split}\hat{\Psi}_{n\mathbf{k}} = \begin{pmatrix} \hat{c}_{n\mathbf{k}\uparrow} \\ \hat{c}^{\dagger}_{n-\mathbf{k}\downarrow} \end{pmatrix}, \qquad \hat{\Psi}^{\dagger}_{n\mathbf{k}} = \begin{pmatrix} \hat{c}^{\dagger}_{n\mathbf{k}\uparrow} & \hat{c}_{n-\mathbf{k}\downarrow} \end{pmatrix},\end{split}\]

where \(\hat{c}^{\dagger}_{n\mathbf{k}\sigma}\) and \(\hat{c}_{n\mathbf{k}\sigma}\) create and annihilate an electron in the state \((n,\mathbf{k})\) with spin \(\sigma\). The corresponding \(2\times2\) matrix Green’s function in imaginary time is

(52)\[\hat{G}_{n\mathbf{k}}(\tau) = -\left\langle T_\tau\,\hat{\Psi}_{n\mathbf{k}}(\tau)\, \hat{\Psi}^{\dagger}_{n\mathbf{k}}(0) \right\rangle,\]

with \(T_\tau\) the imaginary-time ordering operator. Its diagonal elements are the normal Green’s functions describing single-particle electronic excitations, while the off-diagonal elements are the anomalous Green’s functions introduced by Gor’kov [23], which describe the Cooper-pair amplitude and become non-zero below \(T_{\mathrm c}\). A central feature of this representation is that the usual Feynman-Dyson rules of many-body perturbation theory remain valid.

Since \(\hat{G}_{n\mathbf{k}}(\tau)\) is antiperiodic in \(\tau\), it can be expanded in the fermionic Matsubara frequencies \(i\omega_j = i(2j+1)\pi k_{\mathrm B}T\) with \(j\) an integer,

(53)\[\hat{G}_{n\mathbf{k}}(\tau) = k_{\mathrm B}T\sum_{j} e^{-i\omega_j\tau}\hat{G}_{n\mathbf{k}}(i\omega_j).\]

In this representation the Green’s function obeys the Dyson equation

(54)\[\left[\hat{G}_{n\mathbf{k}}(i\omega_j)\right]^{-1} = \left[\hat{G}^{0}_{n\mathbf{k}}(i\omega_j)\right]^{-1} - \hat{\Sigma}_{n\mathbf{k}}(i\omega_j), \qquad \left[\hat{G}^{0}_{n\mathbf{k}}(i\omega_j)\right]^{-1} = i\omega_j\hat{\tau}_0 - (\varepsilon_{n\mathbf{k}}-\mu)\hat{\tau}_3,\]

where \(\hat{\tau}_i\) \((i=0,1,2,3)\) are the Pauli matrices and \(\mu\) is the chemical potential.

Electron self-energy: Fan-Migdal and Coulomb terms

The self-energy entering Eq. (54) contains an electron-phonon and a Coulomb contribution,

(55)\[\hat{\Sigma}_{n\mathbf{k}}(i\omega_j) = \hat{\Sigma}^{\mathrm{ep}}_{n\mathbf{k}}(i\omega_j) + \hat{\Sigma}^{\mathrm{c}}_{n\mathbf{k}}(i\omega_j).\]

Keeping only the first-order (Fan-Migdal) electron-phonon scattering diagram, in line with Migdal’s theorem [18] (see Sec. Migdal’s approximation and beyond), the electron-phonon term reads

(56)\[\hat{\Sigma}^{\mathrm{ep}}_{n\mathbf{k}}(i\omega_j) = -k_{\mathrm B}T\sum_{m\nu}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 D^{0}_{\mathbf{q}\nu}(i\omega_j-i\omega_{j'})\, \hat{\tau}_3\hat{G}_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})\hat{\tau}_3,\]

where \(g_{mn\nu}(\mathbf{k},\mathbf{q})\) are the electron-phonon matrix elements and

(57)\[D^{0}_{\mathbf{q}\nu}(i\omega_j-i\omega_{j'}) = \frac{2\omega_{\mathbf{q}\nu}} {(i\omega_j-i\omega_{j'})^2-\omega_{\mathbf{q}\nu}^2} = -\frac{2\omega_{\mathbf{q}\nu}} {(\omega_j-\omega_{j'})^2+\omega_{\mathbf{q}\nu}^2}\]

is the bare phonon propagator. The Coulomb term is evaluated within the GW approximation,

(58)\[\hat{\Sigma}^{\mathrm{c}}_{n\mathbf{k}}(i\omega_j) = -k_{\mathrm B}T\sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}\, \hat{\tau}_3 \hat{G}^{\mathrm{od}}_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'}) \hat{\tau}_3,\]

where \(V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}\) is the statically screened Coulomb interaction between electrons. Only the off-diagonal (anomalous) components \(\hat{G}^{\mathrm{od}}\) are retained, because the Kohn-Sham eigenvalue \(\varepsilon_{n\mathbf{k}}\) in \(\hat{G}^{0}_{n\mathbf{k}}\) already accounts for part of the Coulomb interaction [20].

It is convenient to absorb the phonon propagator into an anisotropic electron-phonon coupling strength,

(59)\[\begin{split}\begin{aligned} \lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega_j-\omega_{j'}) ={}& -N_{\mathrm F}\sum_\nu \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 D^{0}_{\mathbf{q}\nu}(i\omega_j-i\omega_{j'}) \\ ={}& N_{\mathrm F}\sum_\nu \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 \frac{2\omega_{\mathbf{q}\nu}} {(\omega_j-\omega_{j'})^2+\omega_{\mathbf{q}\nu}^2} \\ ={}& \int_0^\infty d\omega\, \frac{2\omega}{(\omega_j-\omega_{j'})^2+\omega^2}\, \alpha^2F_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega), \end{aligned}\end{split}\]

where the fully anisotropic Eliashberg spectral function is

(60)\[\alpha^2F_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega) = N_{\mathrm F}\sum_\nu \left|g_{mn\nu}(\mathbf{k},\mathbf{q})\right|^2 \delta(\omega-\omega_{\mathbf{q}\nu}).\]

This is the momentum-resolved counterpart of the isotropic \(\alpha^2F(\omega)\) defined in Sec. Self-energies and coupling strength.

Anisotropic Migdal-Eliashberg equations: full-bandwidth formulation

Following the standard Eliashberg procedure, the self-energy is expanded on the Pauli-matrix basis in terms of three scalar functions,

(61)\[\hat{\Sigma}_{n\mathbf{k}}(i\omega_j) = i\omega_j\left[1-Z_{n\mathbf{k}}(i\omega_j)\right]\hat{\tau}_0 + \chi_{n\mathbf{k}}(i\omega_j)\hat{\tau}_3 + \phi_{n\mathbf{k}}(i\omega_j)\hat{\tau}_1,\]

where \(Z_{n\mathbf{k}}\) is the mass renormalization function, \(\chi_{n\mathbf{k}}\) the energy shift, and \(\phi_{n\mathbf{k}}\) the order parameter. The superconducting gap function is \(\Delta_{n\mathbf{k}}(i\omega_j)=\phi_{n\mathbf{k}}(i\omega_j)/ Z_{n\mathbf{k}}(i\omega_j)\). \(Z_{n\mathbf{k}}\) controls the quasiparticle renormalization, \(\chi_{n\mathbf{k}}\) shifts the single-particle energies, and \(\Delta_{n\mathbf{k}}\) opens the superconducting gap. Inserting Eq. (61) into the Dyson equation and inverting the resulting \(2\times2\) matrix gives

(62)\[\hat{G}_{n\mathbf{k}}(i\omega_j) = -\frac{i\omega_jZ_{n\mathbf{k}}(i\omega_j)\hat{\tau}_0 +\left[\varepsilon_{n\mathbf{k}}-\mu +\chi_{n\mathbf{k}}(i\omega_j)\right]\hat{\tau}_3 +\phi_{n\mathbf{k}}(i\omega_j)\hat{\tau}_1} {\Theta_{n\mathbf{k}}(i\omega_j)},\]

with

(63)\[\Theta_{n\mathbf{k}}(i\omega_j) = \left[\omega_jZ_{n\mathbf{k}}(i\omega_j)\right]^2 + \left[\varepsilon_{n\mathbf{k}}-\mu +\chi_{n\mathbf{k}}(i\omega_j)\right]^2 + \left[\phi_{n\mathbf{k}}(i\omega_j)\right]^2.\]

Equating the coefficients of \(\hat{\tau}_0\), \(\hat{\tau}_3\), and \(\hat{\tau}_1\) yields the anisotropic full-bandwidth (FBW) Migdal-Eliashberg equations [24]:

(64)\[Z_{n\mathbf{k}}(i\omega_j) = 1 + \frac{k_{\mathrm B}T}{\omega_jN_{\mathrm F}} \sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \frac{\omega_{j'}Z_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} {\Theta_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} \lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega_j-\omega_{j'}),\]

(65)\[\chi_{n\mathbf{k}}(i\omega_j) = -\frac{k_{\mathrm B}T}{N_{\mathrm F}} \sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \frac{\varepsilon_{m\mathbf{k}+\mathbf{q}}-\mu +\chi_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} {\Theta_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} \lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega_j-\omega_{j'}),\]

(66)\[\begin{split}\begin{aligned} \phi_{n\mathbf{k}}(i\omega_j) ={}& \frac{k_{\mathrm B}T}{N_{\mathrm F}} \sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \frac{\phi_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} {\Theta_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} \\ &\times\left[ \lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega_j-\omega_{j'}) - N_{\mathrm F}V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} \right], \end{aligned}\end{split}\]

(67)\[N_{\mathrm e} = \sum_n\int_{\mathrm{BZ}}\frac{d\mathbf{k}}{\Omega_{\mathrm{BZ}}} \left[ 1 - 2k_{\mathrm B}T\sum_j \frac{\varepsilon_{n\mathbf{k}}-\mu+\chi_{n\mathbf{k}}(i\omega_j)} {\Theta_{n\mathbf{k}}(i\omega_j)} \right].\]

The first three equations are solved self-consistently, starting from \(Z_{n\mathbf{k}}=1\), \(\chi_{n\mathbf{k}}=0\), and a model gap, and iterating until \(\phi_{n\mathbf{k}}\) is converged. The transition temperature \(T_{\mathrm c}\) is the highest temperature at which the equations retain a non-trivial solution \(\phi_{n\mathbf{k}}\neq0\). The gap on the real axis is obtained by Padé continuation or by analytic continuation of the Matsubara-axis solution.

Equation (67) is a constraint that fixes the electron number per unit cell \(N_{\mathrm e}\) and thereby determines the chemical potential. When \(\mu\) is updated self-consistently at every temperature the scheme is denoted FBW+\(\mu\); when it is kept fixed at \(\varepsilon_{\mathrm F}\) it is simply called FBW. Because the electronic states are not restricted to the immediate vicinity of the Fermi surface, the FBW formulation retains the full energy dependence of the density of states and is essential for materials with narrow bands or with van Hove singularities close to \(\varepsilon_{\mathrm F}\), such as the superhydrides [24].

The Coulomb interaction enters only through Eq. (66). Evaluating \(V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}\) at the same level as the electron-phonon interaction is expensive, so it is commonly replaced by the semi-empirical Morel-Anderson pseudopotential [25]

(68)\[\mu_{\mathrm c} = N_{\mathrm F}\left\langle\left\langle V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} \right\rangle\right\rangle_{\mathrm{FS}}, \qquad \mu^{*}_{\mathrm c} = \frac{\mu_{\mathrm c}} {1+\mu_{\mathrm c}\ln(\omega_{\mathrm{el}}/\omega_{\mathrm{ph}})},\]

where \(\langle\langle\cdots\rangle\rangle_{\mathrm{FS}}\) denotes a double Fermi-surface average, \(\omega_{\mathrm{el}}\) is the electronic bandwidth, and \(\omega_{\mathrm{ph}}\) a characteristic phonon frequency.

Constant-DOS Fermi-surface-restricted approximation

The FBW equations can be simplified by restricting the energy range to a narrow window around \(\varepsilon_{\mathrm F}\). Inserting the identity \(\int_{-\infty}^{+\infty}d\varepsilon\, \delta(\varepsilon_{n\mathbf{k}}-\varepsilon)=1\) for each Green’s function and assuming that the density of states is constant within this window, \(N(\varepsilon)\rightarrow N_{\mathrm F}\), one may replace \(\delta(\varepsilon_{n\mathbf{k}}-\varepsilon)\) by \(\delta(\varepsilon_{n\mathbf{k}}-\varepsilon_{\mathrm F})\) and carry out the energy integrals analytically [20]. Using

(69)\[\begin{split}\int_{-\infty}^{+\infty}\!\!d\varepsilon\, \frac{\left\{\begin{array}{c} 1 \\ \varepsilon-\mu+\chi_{n\mathbf{k}}(i\omega_j) \end{array}\right\}} {\left[\omega_jZ_{n\mathbf{k}}(i\omega_j)\right]^2 +\left[\varepsilon-\mu+\chi_{n\mathbf{k}}(i\omega_j)\right]^2 +\left[\phi_{n\mathbf{k}}(i\omega_j)\right]^2} = \left\{\begin{array}{l} \pi/\widetilde{\Theta}_{n\mathbf{k}}(i\omega_j) \\ 0\quad(\text{odd in }\varepsilon) \end{array}\right\},\end{split}\]

with \(\widetilde{\Theta}_{n\mathbf{k}}(i\omega_j) =\sqrt{[\omega_jZ_{n\mathbf{k}}(i\omega_j)]^2 +[\phi_{n\mathbf{k}}(i\omega_j)]^2}\), the energy-shift equation vanishes identically and the particle-number constraint is automatically satisfied. Only two coupled equations remain, the anisotropic Fermi-surface-restricted (FSR) Migdal-Eliashberg equations [21]:

(70)\[\begin{split}\begin{aligned} Z_{n\mathbf{k}}(i\omega_j) ={}& 1 + \frac{\pi k_{\mathrm B}T}{\omega_jN_{\mathrm F}} \sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \frac{\omega_{j'}} {\sqrt{\omega_{j'}^2 +\Delta^2_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})}} \\ &\times \lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega_j-\omega_{j'})\, \delta(\varepsilon_{m\mathbf{k}+\mathbf{q}}-\varepsilon_{\mathrm F}), \end{aligned}\end{split}\]

(71)\[\begin{split}\begin{aligned} Z_{n\mathbf{k}}(i\omega_j)\Delta_{n\mathbf{k}}(i\omega_j) ={}& \frac{\pi k_{\mathrm B}T}{N_{\mathrm F}} \sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} \frac{\Delta_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})} {\sqrt{\omega_{j'}^2 +\Delta^2_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'})}} \\ &\times\left[ \lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\omega_j-\omega_{j'}) - N_{\mathrm F}V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} \right] \delta(\varepsilon_{m\mathbf{k}+\mathbf{q}}-\varepsilon_{\mathrm F}). \end{aligned}\end{split}\]

The Dirac delta functions project the equations onto the Fermi surface, so that only states within a chosen energy window contribute. The FSR equations remain accurate whenever the density of states varies slowly on the scale of the phonon energies.

Intermediate representation and sparse sampling

All quantities in the Migdal-Eliashberg equations are defined on Matsubara frequencies, whose spacing is proportional to \(T\). Covering a fixed energy range therefore requires a number of frequencies growing as \(T^{-1}\), which becomes prohibitive at low temperature, and especially in the FBW+\(\mu\) scheme, where a cutoff 6-12 times the Fermi energy window is needed to converge \(\mu\). Two complementary strategies are available in EPW to control this cost.

Intermediate representation. The FBW equations can be recast in terms of the normal and anomalous Green’s functions,

(72)\[G_{n\mathbf{k}}(i\omega_j) = -\frac{i\omega_jZ_{n\mathbf{k}}(i\omega_j) +\varepsilon_{n\mathbf{k}}-\mu+\chi_{n\mathbf{k}}(i\omega_j)} {\Theta_{n\mathbf{k}}(i\omega_j)}, \qquad F_{n\mathbf{k}}(i\omega_j) = -\frac{\phi_{n\mathbf{k}}(i\omega_j)} {\Theta_{n\mathbf{k}}(i\omega_j)},\]

so that the normal and pairing self-energies become plain convolutions in Matsubara frequency,

(73)\[\begin{split}\begin{aligned} \Sigma^{\mathrm N}_{n\mathbf{k}}(i\omega_j) ={}& -k_{\mathrm B}T\sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} G_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'}) V^{\mathrm{ep}}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} (i\omega_j-i\omega_{j'}), \\ \phi_{n\mathbf{k}}(i\omega_j) ={}& k_{\mathrm B}T\sum_{m}\sum_{j'} \int_{\mathrm{BZ}}\frac{d\mathbf{q}}{\Omega_{\mathrm{BZ}}} F_{m\mathbf{k}+\mathbf{q}}(i\omega_{j'}) \left[ V^{\mathrm{ep}}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} (i\omega_j-i\omega_{j'}) + V^{\mathrm c}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} \right], \end{aligned}\end{split}\]

where \(V^{\mathrm{ep}}_{n\mathbf{k},m\mathbf{k}+\mathbf{q}} =-\lambda_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}/N_{\mathrm F}\) is the electron-phonon kernel of Eq. (59). The renormalization and energy-shift functions follow from the odd and even parts of \(\Sigma^{\mathrm N}_{n\mathbf{k}}\),

(74)\[i\omega_j\left[1-Z_{n\mathbf{k}}(i\omega_j)\right] = \tfrac{1}{2}\left[ \Sigma^{\mathrm N}_{n\mathbf{k}}(i\omega_j) -\Sigma^{\mathrm N}_{n\mathbf{k}}(-i\omega_j)\right], \qquad \chi_{n\mathbf{k}}(i\omega_j) = \tfrac{1}{2}\left[ \Sigma^{\mathrm N}_{n\mathbf{k}}(i\omega_j) +\Sigma^{\mathrm N}_{n\mathbf{k}}(-i\omega_j)\right].\]

A convolution in frequency is a product in imaginary time,

(75)\[\sum_{j'}f(i\omega_j-i\omega_{j'})\,g(i\omega_{j'}) = \mathcal{F}^{-1} \left[\mathcal{F}(f)\cdot\mathcal{F}(g)\right],\]

where \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) denote the Fourier transformations between the Matsubara-frequency and imaginary-time domains. These transformations are performed compactly in the intermediate representation (IR) [26], in which any imaginary-time correlation function is expanded on a small set of material-independent basis functions,

(76)\[G(i\bar{\omega}_j) = \sum_l G_l\,\hat{U}_l(i\bar{\omega}_j), \qquad G(\bar{\tau}_i) = \sum_l G_l\,U_l(\bar{\tau}_i),\]

with the expansion coefficients obtained from either domain through the pseudo-inverse \(U^{+}=(U^{\dagger}U)^{-1}U^{\dagger}\),

(77)\[G_l = \sum_j\left[\hat{U}^{+}\right]_{lj}G(i\bar{\omega}_j) = \sum_i\left[U^{+}\right]_{li}G(\bar{\tau}_i).\]

The sampling points \(\{\bar{\tau}_i\}\) and \(\{\bar{\omega}_j\}\) and the basis functions are generated beforehand with the sparse-ir library [27] from two parameters, \(\Lambda=\omega_{\max}/T\) and a truncation threshold \(\varepsilon_{\mathrm{IR}}\). The number of IR coefficients grows only logarithmically with \(\Lambda\), so frequency ranges of several tens of eV are covered with roughly one hundred sampling points. This makes anisotropic FBW calculations tractable well below 1 K and allows the electron-phonon and Coulomb interactions to be treated consistently [28].

Sparse sampling of Matsubara frequencies. As a simpler alternative, EPW can prune the uniform Matsubara grid directly [24]. Denoting by \(N_j\) the integer index of the \(j\)-th frequency, \(i\omega_j=i(2N_j+1)\pi k_{\mathrm B}T\), the uniform and sparse grids are generated as

(78)\[N_j = j, \qquad j=0,\pm1,\pm2,\pm3,\ldots\]

and

(79)\[N_j = N_{j-1} + \mathrm{INT}\left[\exp\left(\frac{j}{N_{\max}W}\right)\right], \qquad N_0=0,\quad j=1,2,3,\ldots\]

where \(\mathrm{INT}[\;]\) rounds to the closest integer, \(W\) is an adjustable weight factor, and

(80)\[N_{\max} = \mathrm{INT}\left[ \frac{1}{2}\left(\frac{\omega_{\max}}{\pi k_{\mathrm B}T}-1\right) \right]\]

is the largest Matsubara index for a given cutoff \(\omega_{\max}\) and temperature \(T\). Negative indices follow from \(\omega_{-(j+1)}=-\omega_j\) with \(j>0\). The resulting mesh is uniform at low frequencies, which dominate the Matsubara sums, and becomes logarithmically sparser at high frequencies. Increasing (decreasing) \(W\) produces a denser (sparser) grid; with the default \(W=1.0\) about 30% fewer frequencies are needed than in the uniform scheme, reducing the computational cost by roughly 40% while reproducing the gap function of the uniform grid.

Isotropic approximation

For most materials, apart from a few strongly anisotropic layered systems, the momentum dependence of the gap is weak and an isotropic treatment is sufficient. The isotropic average of a quantity \(A_{n\mathbf{k}}\) replaces the sums over momenta and bands by an energy integral weighted by the density of states,

(81)\[A(\varepsilon) = \frac{1}{N(\varepsilon)} \sum_{n\mathbf{k}}A_{n\mathbf{k}}\, \delta(\varepsilon_{n\mathbf{k}}-\varepsilon), \qquad N(\varepsilon)=\sum_{n\mathbf{k}} \delta(\varepsilon_{n\mathbf{k}}-\varepsilon).\]

The energy dependence of the electron-phonon kernel is neglected and it is evaluated at the Fermi level [20], \(\lambda(\varepsilon,\varepsilon',i\omega_j-i\omega_{j'})\approx \lambda(i\omega_j-i\omega_{j'})\), which decouples the kernel from the Green’s function in the energy integrals. The isotropic kernel is the double Fermi-surface average of Eq. (59),

(82)\[\lambda(i\omega_j-i\omega_{j'}) = \int_0^\infty d\omega\, \frac{2\omega}{(\omega_j-\omega_{j'})^2+\omega^2}\, \alpha^2F(\omega),\]

with \(\alpha^2F(\omega)\) the isotropic Eliashberg spectral function defined in Sec. Self-energies and coupling strength. The kernel is even in the Matsubara frequency and reduces at \(\omega_j=\omega_{j'}\) to the isotropic coupling strength \(\lambda=2\int_0^\infty d\omega\,\alpha^2F(\omega)/\omega\).

Isotropic FBW equations. Retaining the full energy dependence of the density of states gives [11]

(83)\[Z(i\omega_j) = 1 + \frac{k_{\mathrm B}T}{N_{\mathrm F}\omega_j} \sum_{j'}\lambda(i\omega_j-i\omega_{j'}) \int_{-\infty}^{\infty}d\varepsilon'\,N(\varepsilon')\, \gamma^{Z}(\varepsilon',i\omega_{j'}),\]

(84)\[\chi(i\omega_j) = -\frac{k_{\mathrm B}T}{N_{\mathrm F}} \sum_{j'}\lambda(i\omega_j-i\omega_{j'}) \int_{-\infty}^{\infty}d\varepsilon'\,N(\varepsilon')\, \gamma^{\chi}(\varepsilon',i\omega_{j'}),\]

(85)\[\phi(i\omega_j) = \frac{k_{\mathrm B}T}{N_{\mathrm F}} \sum_{j'}\left[\lambda(i\omega_j-i\omega_{j'})-\mu^{*}_{\mathrm c}\right] \int_{-\infty}^{\infty}d\varepsilon'\,N(\varepsilon')\, \gamma^{\phi}(\varepsilon',i\omega_{j'}),\]

(86)\[N_{\mathrm e} = \int_{-\infty}^{\infty}d\varepsilon'\,N(\varepsilon') \left[1-2k_{\mathrm B}T\sum_{j'} \frac{\varepsilon'-\mu+\chi(i\omega_{j'})} {\Theta(\varepsilon',i\omega_{j'})}\right],\]

where

(87)\[\gamma^{Z}(\varepsilon,i\omega_j) = \frac{\omega_jZ(i\omega_j)}{\Theta(\varepsilon,i\omega_j)}, \qquad \gamma^{\chi}(\varepsilon,i\omega_j) = \frac{\varepsilon-\mu+\chi(i\omega_j)}{\Theta(\varepsilon,i\omega_j)}, \qquad \gamma^{\phi}(\varepsilon,i\omega_j) = \frac{\phi(i\omega_j)}{\Theta(\varepsilon,i\omega_j)},\]

and

(88)\[\Theta(\varepsilon,i\omega_j) = \left[\omega_jZ(i\omega_j)\right]^2 + \left[\varepsilon-\mu+\chi(i\omega_j)\right]^2 + \left[\phi(i\omega_j)\right]^2.\]

Isotropic FSR equations. Taking the constant-DOS limit \(N(\varepsilon)\rightarrow N_{\mathrm F}\) and performing the energy integrals analytically as in Eq. (69), the energy shift vanishes and only two equations survive:

(89)\[Z(i\omega_j) = 1 + \frac{\pi k_{\mathrm B}T}{\omega_j} \sum_{j'}\lambda(i\omega_j-i\omega_{j'}) \frac{\omega_{j'}}{\sqrt{\omega_{j'}^2+\Delta^2(i\omega_{j'})}},\]

(90)\[Z(i\omega_j)\Delta(i\omega_j) = \pi k_{\mathrm B}T\sum_{j'} \left[\lambda(i\omega_j-i\omega_{j'})-\mu^{*}_{\mathrm c}\right] \frac{\Delta(i\omega_{j'})} {\sqrt{\omega_{j'}^2+\Delta^2(i\omega_{j'})}}.\]

These are the isotropic Migdal-Eliashberg equations solved in EPW, and are the form most directly comparable with the McMillan and Allen-Dynes expressions for \(T_{\mathrm c}\).

Migdal’s approximation and beyond

All the equations above rest on Migdal’s theorem [18], which assumes that electrons and ions evolve on two well-separated energy scales, in accordance with the adiabatic Born-Oppenheimer approximation. Migdal showed that each additional order in the perturbative expansion of the electron-phonon vertex is suppressed by a factor of order

(91)\[\lambda\,\frac{\hbar\omega_0}{\varepsilon_{\mathrm F}} \sim \lambda\sqrt{\frac{m_{\mathrm e}}{M}},\]

where \(\omega_0\) is a characteristic phonon frequency and \(m_{\mathrm e}\) and \(M\) are the electron and ion masses. Since electronic energy scales are of order several eV while phonon energies are a fraction of an eV, \(\lambda\hbar\omega_0/\varepsilon_{\mathrm F}\ll1\) holds in most metals and the vertex corrections can safely be dropped.

This condition breaks down in systems with a small effective Fermi energy or very stiff phonons, such as superhydrides, fullerides, and low-carrier-density superconductors. In H₃S, for instance, a flat-band feature produces van Hove singularities near \(\varepsilon_{\mathrm F}\) whose width sets an effective Fermi energy of only a few tens of meV, while the characteristic phonon energy is \(\sim140\) meV, giving an adiabatic parameter \(\hbar\omega_0/\varepsilon_{\mathrm F}\) of order unity or larger [29], [30]. The natural generalization is a perturbative expansion in \(\lambda\hbar\omega_0/\varepsilon_{\mathrm F}\), in which the lowest-order (crossed) vertex-correction diagram is added to the Fan-Migdal self-energy. This diagram involves four electron-phonon matrix elements and two phonons, and is characterized by a vertex-corrected Eliashberg spectral function \(\alpha^2F^{\mathrm V}(\omega,\omega')\) and the associated coupling strength

(92)\[\lambda^{\mathrm V} = 4\int_0^\infty d\omega\int_0^\infty d\omega'\, \frac{\alpha^2F^{\mathrm V}(\omega,\omega')}{\omega\,\omega'}.\]

Implementation of the isotropic vertex-corrected equations in EPW is described in Ref. [29]. For H₃S the vertex-corrected coupling remains a large fraction of the adiabatic one (\(\lambda^{\mathrm V}=1.78\) versus \(\lambda=2.29\) with harmonic phonons), and combining vertex corrections with phonon anharmonicity and the energy dependence of the density of states brings the predicted \(T_{\mathrm c}\) into agreement with experiment for both H₃S and D₃S [30]. For elemental Pb, where the adiabatic limit holds, the vertex corrections are negligible and the standard Migdal-Eliashberg results are recovered. Vertex corrections are therefore worth the added cost precisely in the regime where they are expected to matter, and the choice between the adiabatic and non-adiabatic treatment can be guided by the size of \(\lambda\hbar\omega_0/\varepsilon_{\mathrm F}\) for the material at hand.

Boltzmann transport equation

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.

[18] A. B. Migdal, Interaction between electrons and lattice vibrations in a normal metal, Sov. Phys. JETP 7, 996 (1958).

[19] G. M. Eliashberg, Interactions between electrons and lattice vibrations in a superconductor, Sov. Phys. JETP 11, 696 (1960).

[20] P. B. Allen and B. Mitrović, Theory of superconducting \(T_c\), in Solid State Physics, Vol. 37 (Academic Press, New York, 1982), p. 1, doi:10.1016/S0081-1947(08)60665-7.

[21] E. R. Margine and F. Giustino, Anisotropic Migdal-Eliashberg theory using Wannier functions, Phys. Rev. B 87, 024505 (2013), doi:10.1103/PhysRevB.87.024505.

[22] Y. Nambu, Quasi-particles and gauge invariance in the theory of superconductivity, Phys. Rev. 117, 648 (1960), doi:10.1103/PhysRev.117.648.

[23] L. P. Gor’kov, On the energy spectrum of superconductors, Sov. Phys. JETP 7, 505 (1958).

[24] R. Lucrezi et al., Full-bandwidth anisotropic Migdal-Eliashberg theory and its application to superhydrides, Commun. Phys. 7, 33 (2024), doi:10.1038/s42005-024-01528-6.

[25] P. Morel and P. W. Anderson, Calculation of the superconducting state parameters with retarded electron-phonon interaction, Phys. Rev. 125, 1263 (1962), doi:10.1103/PhysRev.125.1263.

[26] H. Shinaoka, J. Otsuki, M. Ohzeki, and K. Yoshimi, Compressing Green’s function using intermediate representation between imaginary-time and real-frequency domains, Phys. Rev. B 96, 035147 (2017), doi:10.1103/PhysRevB.96.035147.

[27] M. Wallerberger et al., sparse-ir: Optimal compression and sparse sampling of many-body propagators, SoftwareX 21, 101266 (2023), doi:10.1016/j.softx.2022.101266.

[28] H. Mori, T. Nomoto, R. Arita, and E. R. Margine, Efficient anisotropic Migdal-Eliashberg calculations with an intermediate representation basis and Wannier interpolation, Phys. Rev. B 110, 064505 (2024), doi:10.1103/PhysRevB.110.064505.

[29] S. B. Mishra, H. Mori, and E. R. Margine, Electron-phonon vertex correction effect in superconducting H₃S, npj Comput. Mater. 11, 342 (2025), doi:10.1038/s41524-025-01818-9.

[30] S. B. Mishra and E. R. Margine, Nonadiabatic and anharmonic effects in high-pressure H₃S and D₃S superconductors, Ann. Phys. (Berlin) 538, e00553 (2026), doi:10.1002/andp.202500553.