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\). The first-order variation of the self-consistent Kohn-Sham potential with respect to the normal coordinate of the phonon mode \((\mathbf{q}, \nu)\) is

(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

Marios Zacharias 08/11/2026

The Zacharias-Giustino (ZG) special displacement method enables the construction of an optimal supercell configuration that accurately represents the properties of the system at thermal equilibrium [31], [32], [33]. The special displacement method incorporates efficiently the effect of electron-phonon coupling in a non-perturbative fashion. For notation and conventions please refer here.

Supercell formulation

The concept of using a single set of atomic displacements, later referred to as ZG displacements \(\Delta \tau^{\rm ZG}_{\kappa\alpha p}\), was originally introduced in Ref. [31] for the calculation of phonon-induced band gap renormalization and phonon-assisted optical absorption. Within a supercell formulation, the ZG displacements at temperature \(T\) are given by:

(93)\[\Delta \tau^{\rm ZG}_{\kappa\alpha} = (M_0/M_\kappa)^{1/2} \sum_{\nu} (-1)^{\nu-1} e_{\kappa \alpha,\nu} \, \sigma_{\nu,T}\]

where indices \(p\) and \({\mathbf q}\) are set to 0 and therefore they are omitted. \(\sigma^2_{\nu,T}\) is the mean square displacement of the atoms along a single phonon mode defined as:

(94)\[ \sigma^2_{\nu,T} = l^2_{\nu} \, (2n_{\nu}+1),\]

where the temperature dependence of \(\sigma_{\nu,T}\) (and hence \(\Delta \tau^{\rm ZG}_{\kappa\alpha}\)) enters through the Bose-Einstein occupation factor \(n_{\nu}\).

The special displacement method essentially replaces the Monte Carlo averaging of an observable [32] with its evaluation using a single ZG configuration, i.e.:

(95)\[ \left\langle O \right\rangle_T = \prod_\nu \int dx_\nu \frac{e^{-x_\nu^2/2\sigma^2_\nu}}{ \sigma_\nu \sqrt{2\pi}} \, O({\{x\}}) \simeq O({\{ \boldsymbol{\tau}_{\kappa}^{\rm ZG}\}}),\]

where \(x_\nu\) is the normal coordinate associated with the phonon mode \(\nu\) and \(\{\boldsymbol{\tau}_{\kappa}^{\rm ZG}\}\) represents the set of atomic coordinates generated using the displacements in Eq. (93). The above expression becomes exact in the thermodynamic limit, i.e. for a large supercell size, as demonstrated in Ref. [31]. For generating ZG displacements via Eq. (93), one needs to calculate the phonon modes at the \(\Gamma\)-point for a given supercell size and then enforce the same choice of gauge for each vibrational mode.

Reciprocal space formulation

The reciprocal-space formulation and formal theory of the special displacement method has been developed in Ref. [33]. The theory, implemented in ZG.x, enables the generation of ZG configurations for arbitrary supercell sizes, provided that the phonons have been computed on a given phonon \({\mathbf q}\)-grid. In other words, the ZG supercell does not need to be commensurate with the phonon grid used in the phonon calculation, since ZG.x exploits Fourier interpolation of both the long-range and short-range components of the dynamical matrix. Within the reciprocal space formulation, the ZG displacements are given by:

(96)\[ \Delta \boldsymbol{\tau}^{\rm ZG}_{\kappa p} = \left[\!\frac{M_0}{ N_p M_\kappa } \!\right]^{\!\frac{1}{2}} 2 \!\sum_{ {\mathbf q} \in \mathcal{B}, \nu}\! S_{{\mathbf q} \nu} \,{\rm Re} \Big[ e^{i {\mathbf q} \cdot {\mathbf R}_p} {\bf e}_{\kappa,\nu} ({\mathbf q} ) \Big] \sigma_{{\mathbf q} \nu, T} + \left[\!\frac{M_0}{ N_p M_\kappa } \!\right]^{\!\frac{1}{2}} \!\sum_{{\mathbf q} \in \mathcal{A}, \nu}\! S_{{\mathbf q} \nu} {\mathbf e}_{\kappa,\nu} ({\mathbf q} ) \sigma_{{\mathbf q} \nu, T}.\]

Here, set \(\mathcal{B}\) includes the wave vectors \({\mathbf q}\) whose time-reversal partners are not contained within the same set, while a set \(\mathcal{C}\) includes the corresponding time-reversed wave vectors obtained by reversing the sign of each wave vector in set \(\mathcal{B}\). Set \(\mathcal{A}\) includes phonon wave vectors that coincide with their time-reversal partners, modulo a reciprocal lattice vector. For the first term on the right hand side, the real part and factor of 2 arise from grouping together a phonon with wave vector \({\mathbf q}\) that belongs to set \(\mathcal{B}\) and its time-reversal partner at \(-{\mathbf q}\) that belongs to set \(\mathcal{C}\). The code can apply a smooth Berry connection between phonon eigenmodes by applying unitary transformations for adjacent \({\mathbf q}\)-points based on singular value decomposition. The smooth Berry connection accelerates the combinatorial search for an optimal set of signs \(S_{{\mathbf q} \nu}= \pm 1\), so that the following error function is minimized with respect to a predetermined threshold:

(97)\[\begin{split} E(\{S_{{\mathbf q}\nu}\},T) = \sum_{\substack{\kappa\alpha \\ \kappa'\alpha'}} \frac{\Big|\displaystyle\sum_{\substack{{\mathbf q} \in \mathcal{A},\mathcal{B}\\\nu<\nu'}} {\rm Re} [ e^*_{\kappa\alpha,\nu} ({\mathbf q} ) e_{\kappa'\alpha',\nu'} ({\mathbf q} ) ] \sigma_{{\mathbf q} \nu,T} \sigma_{{\mathbf q} \nu',T} S_{{\mathbf q}\nu}S_{{\mathbf q}\nu'}\Big|} {\Big|\displaystyle\sum_{\substack{{\mathbf q}\in\mathcal{A},\mathcal{B}\\\nu}} {\rm Re} [ e^*_{\kappa\alpha,\nu} ({\mathbf q} ) e_{\kappa'\alpha',\nu} ({\mathbf q} ) ] \sigma_{{\mathbf q} \nu,T}^2 \Big|}\quad.\end{split}\]

This minimization ensures that the ZG displacements closely reproduce the exact mean-square anisotropic displacement tensor at temperature \(T\) [33]:

(98)\[{\mathcal U}_{\kappa,\alpha\alpha'} (T) = \frac{M_0}{M_\kappa} \sum_\nu \int_{\rm BZ} \frac{d{\mathbf q}}{\Omega_{\rm BZ}} e_{\kappa\alpha,\nu} ({\mathbf q} ) e^*_{\kappa\alpha',\nu} ({\mathbf q} ) \sigma_{{\mathbf q} \nu,T}^2\]

and the diffuse (phonon-induced inelastic) scattering intensity [34]:

(99)\[ \left\langle I({\mathbf Q}) \right\rangle_T \simeq I_{\rm ZG}({\mathbf Q},T) = \Big| \sum_{\kappa p} f_\kappa({\mathbf Q}) e^{i {\mathbf Q} \cdot ({\mathbf R}_p + \boldsymbol{\tau}_\kappa + \Delta \boldsymbol{\tau}^{\rm ZG}_{\kappa p} )} \Big|^2,\]

where \({\mathbf Q}\) is the scattering wavevector and \(f_\kappa({\mathbf Q})\) is the atomic scattering factor.

The special displacement method for treating electron–phonon coupling offers both strengths and weaknesses compared with state-of-the-art perturbative approaches. The most important include:

Strengths

• As a non-perturbative approach, the method naturally accounts for higher-order electron-phonon coupling contributions in the evaluation of the observable.

• Includes the effect of electron-phonon coupling on the electronic structure simultaneously with the calculation of the target observable, rather than treating the renormalization as a separate correction.

• Can be straightforwardly combined with anharmonic phonons (see ASDM) to incorporate anharmonic electron-phonon coupling effects.

Weaknesses

• Relies on supercell calculations and therefore does not retain the computational efficiency and elegance of the unit-cell formulation employed in perturbative approaches, particularly with respect to dense sampling of the \({\mathbf q}\)-point and \({\mathbf k}\)-point grids.

• Relies on the adiabatic approximation and therefore misses non-adiabatic effects.

Anharmonic special displacement method

Marios Zacharias 08/11/2026

Equations (93) and (96) become ill-defined in the presence of imaginary phonon frequencies. These can originate from numerical inaccuracies (e.g., insufficient convergence of the k-mesh, energy cutoff, phonon grid, or acoustic sum rule used) or from a true dynamical instability indicating that the harmonic approximation is inadequate. In the latter case, anharmonic effects must be taken into account to obtain a physically meaningful description of the vibrational and electron-phonon properties. A convenient framework for this purpose is provided by the self-consistent phonon theory [35] in which anharmonic interactions are incorporated through a renormalization of the phonon frequencies and eigenvectors. Some of the ab initio implementations of the self-consistent phonon theory can be found in Refs. [36], [37], [38], [39], [40], and [41]. In the following, please refer here for notation and conventions.

Self-consistent phonon theory

Within the self-consistent phonon (SCP) theory, anharmonic effects are incorporated through a temperature-dependent effective harmonic Hamiltonian. The corresponding effective force constants are obtained variationally by minimizing a trial free energy based on the Gibbs–Bogoliubov inequality.

The Gibbs–Bogoliubov inequality states that:

(100)\[F_{\mathrm{exact}}(T) \leq F_{\mathrm{trial}}(T),\]

where \(F_{\mathrm{exact}}(T)\) is the exact anharmonic free energy and \(F_{\mathrm{trial}}(T)\) is the variational trial free energy associated with an effective harmonic Hamiltonian. The latter therefore provides an upper bound to the exact free energy.

The optimal effective harmonic system is obtained by minimizing this upper bound with respect to the temperature-dependent force constants \(C_{\kappa\alpha p,\kappa'\alpha' p'}(T)\).

The trial free energy can be written as:

(101)\[F_{\mathrm{trial}}(T) = \langle U \rangle_T - U_{\mathrm{h}}(T) + F_{\mathrm{vib}}(T),\]

where \(\langle U \rangle_T\) is the thermal average of the full anharmonic potential-energy surface evaluated over the effective harmonic distribution, \(U_{\mathrm{h}}(T)\) is the thermal average of the effective harmonic potential, and \(F_{\mathrm{vib}}(T)\) is the vibrational free energy of the effective harmonic system.

Using the normal-mode representation, the trial free energy becomes:

(102)\[F_{\mathrm{trial}}(T) = \langle U \rangle_T - \frac{M_0}{2} \sum_{\nu} \omega_{\nu}^{2}\sigma_{\nu,T}^{2} + \sum_{\nu} \left[ \frac{\hbar\omega_{\nu}}{2} - k_{\mathrm{B}}T \ln\left(1+n_{\nu}\right) \right].\]

The temperature-dependent effective force constants are determined by imposing the stationarity condition:

(103)\[\frac{\partial F_{\mathrm{trial}}(T)} {\partial C_{\kappa\alpha p,\kappa'\alpha' p'}} = 0.\]

The derivative of the trial free energy can be written as:

(104)\[\frac{\partial F_{\mathrm{trial}}(T)} {\partial C_{\kappa\alpha p, \kappa'\alpha' p'}} = \frac{\partial \langle U \rangle_T} {\partial C_{\kappa\alpha p,\kappa'\alpha' p'}} - \frac{\partial U_{\mathrm{h}}(T)} {\partial C_{\kappa\alpha p,\kappa'\alpha' p'}} + \frac{\partial F_{\mathrm{vib}}(T)} {\partial C_{\kappa\alpha p,\kappa'\alpha' p'}}.\]

Applying the stationarity condition yields the self-consistent relation ([40] and [41]):

(105)\[C_{\kappa\alpha p,\kappa'\alpha' p'}(T) = \left\langle \frac{\partial^{2} U(\{\boldsymbol{\tau}_{\kappa p}\})} {\partial \tau_{\kappa\alpha p} \partial \tau_{\kappa'\alpha' p'}} \right\rangle_T.\]

Therefore, the effective harmonic force constants correspond to the thermal average of the Hessian of the full anharmonic potential-energy surface.

Importantly, the thermal average on the right-hand side is evaluated using the Gaussian distribution defined by the effective harmonic Hamiltonian itself. The effective force constants therefore determine the phonon frequencies and eigenvectors entering the same thermal distribution used to evaluate the anharmonic average.

Eq. (105) must consequently be solved iteratively. Starting from an initial set of force constants, the corresponding phonon frequencies and eigenvectors are calculated and used to evaluate:

(106)\[\left\langle \frac{\partial^{2} U(\{\boldsymbol{\tau}_{\kappa p}\})} {\partial \tau_{\kappa\alpha p} \partial \tau_{\kappa'\alpha' p'}} \right\rangle_T,\]

The resulting temperature-dependent force constants define a new effective harmonic Hamiltonian and therefore a new vibrational distribution. This procedure is repeated until \(C_{\kappa\alpha p,\kappa'\alpha' p'}(T)\) is converged, yielding the self-consistent anharmonic phonons at temperature \(T\).

Within this formulation, the anharmonic potential-energy surface does not need to be explicitly truncated at a particular order. The SCP equation therefore contains all first-order anharmonic self-energy diagrams generated by the full potential-energy surface, with the quartic loop diagram providing the leading contribution. Higher even-order anharmonic terms generate the corresponding higher-order first-order diagrams. By contrast, cubic anharmonicity does not contribute at first order because odd moments of the harmonic Gaussian distribution vanish. The cubic contribution first appears at second order through the bubble diagram, which involves two third-order anharmonic vertices. This diagram is therefore not contained in the standard SCP self-consistency condition and must be included either as an additional perturbative self-energy correction or, in the static limit, through the second derivatives of the anharmonic free energy with respect to the centroid atomic positions [40].

Anharmonic ZG displacements

The anharmonic special displacement method (ASDM) combines the special displacement method with self-consistent phonon theory to efficiently calculate temperature-dependent anharmonic phonons and thus generate anharmonic ZG displacements. The central idea is to replace the thermal average entering the SCP equations by an evaluation at a single, optimally chosen ZG configuration.

Within SCP theory, the temperature-dependent effective interatomic force constants (IFCs) are given by Eq. (106). In the special displacement formulation, this thermal average is approximated by evaluating the curvature of the potential-energy surface at the ZG configuration:

\[C_{\kappa\alpha p,\kappa'\alpha' p'}(T) = \left\langle \frac{\partial^2 U(\{\boldsymbol{\tau}_{\kappa p}\})} {\partial \tau_{\kappa\alpha p} \partial \tau_{\kappa'\alpha' p'}} \right\rangle_T \simeq \left. \frac{\partial^2 U(\{\boldsymbol{\tau}_{\kappa p}\})} {\partial \tau_{\kappa\alpha p} \partial \tau_{\kappa'\alpha' p'}} \right|_{\boldsymbol{\tau}_{\kappa p}=\boldsymbol{\tau}^{\mathrm ZG}_{\kappa p}}.\]

Therefore, instead of explicitly performing the multidimensional thermal average over nuclear configurations, ASDM evaluates the effective IFCs from a single ZG configuration representative of the nuclear thermal distribution at temperature \(T\).

To this aim, ZG.x is combined with pw.x to generate the temperature-dependent ZG configurations and calculate the corresponding effective IFCs iteratively until self-consistent phonons are obtained.

The ASDM procedure consists of the following steps:

1. Choose the temperature and supercell

   Select the target temperature \(T\) and the supercell size to be used for the calculation. The supercell determines the real-space range of the effective IFCs and the corresponding sampling of vibrational modes.

2. Search for a stable reference structure

   Starting from an initial ZG configuration, perform a geometry optimization to identify a stable structure. In systems exhibiting strong structural instabilities, this relaxation can lower the symmetry and lead to a polymorphous structure.This step provides a stable reference configuration around which the initial vibrational problem is constructed.

3. Calculate the initial harmonic phonons

   Compute the harmonic IFCs \(C_{\kappa\alpha p,\kappa'\alpha' p'}\) of the relaxed structure using finite differences. The appropriate crystal symmetries and sum rules can then be imposed on the calculated IFCs. Diagonalization of the corresponding dynamical matrix provides the initial set of phonon frequencies and eigenvectors:

   \[\left\{ \omega_{\mathbf{q}\nu}, e_{\kappa\alpha,\nu}(\mathbf{q}) \right\}_{\rm Init}.\]

   These phonons define the initial nuclear thermal distribution used to construct the first finite-temperature ZG configuration.

4. Generate the finite-temperature ZG configuration

   Using the current phonon frequencies and eigenvectors, generate a new set of ZG displacements \(\Delta\tau^{\rm ZG}_{\kappa\alpha p}(T)\), according to Eq. (96). The atomic positions of the corresponding ZG configuration are therefore:

   \[\tau^{\rm ZG}_{\kappa\alpha p}(T) = \tau^{0}_{\kappa\alpha p} + \Delta\tau^{\rm ZG}_{\kappa\alpha p}(T).\]

   The ZG configuration provides a single configuration designed to represent the thermal fluctuations of the nuclei at the target temperature.

5. Calculate the temperature-dependent effective IFCs

   For the generated ZG configuration, calculate the effective IFCs using finite differences:

   \[C_{\kappa\alpha p,\kappa'\alpha' p'}(T) \simeq \left. \frac{\partial^2 U(\{\boldsymbol{\tau}_{\kappa p}\})} {\partial \tau_{\kappa\alpha p} \partial \tau_{\kappa'\alpha' p'}} \right|_{\boldsymbol{\tau}_{\kappa p}=\boldsymbol{\tau}^{\mathrm ZG}_{\kappa p}}.\]

   Crystal symmetries are subsequently applied to the calculated effective IFCs. Under a symmetry operation \(S\), the force constants transform according to:

   \[C_{PK\alpha,P'K'\alpha'} = \sum_{\beta\beta'} S_{\alpha\beta} S_{\alpha'\beta'} C_{\kappa\beta p,\kappa'\beta' p'}.\]

   Here, \((\kappa p,\kappa' p')\) and \((KP,K'P')\) denote pairs of atoms related by the corresponding crystal symmetry operation.

6. Calculate the anharmonic phonons

   Construct and diagonalize the dynamical matrix using the temperature-dependent effective IFCs. This yields a new set of temperature-dependent phonon frequencies and eigenvectors:

   \[\left\{ \omega_{\mathbf{q}\nu}(T), e_{\kappa\alpha,\nu}(\mathbf{q},T) \right\}.\]

   These phonons contain the anharmonic renormalization encoded in the temperature-dependent effective IFCs.

7. Iterate to self-consistency

   The updated phonon frequencies and eigenvectors are used to construct a new ZG displacement, \(\Delta\tau^{\rm ZG}_{\kappa\alpha p}(T)\), from which a new ZG configuration and a new set of effective IFCs are obtained. Steps 4–6 are repeated until the phonon frequencies, eigenvectors, or effective IFCs satisfy the desired convergence criterion. Iterative mixing of the IFCs can be employed to improve the numerical stability of the self-consistent cycle. Schematically, the ASDM cycle can be written as:

   \[\left\{ \omega_{\mathbf{q}\nu}(T), e_{\kappa\alpha,\nu}(\mathbf{q},T) \right\}^{(n)} \rightarrow \tau_{\kappa\alpha p}^{{\rm ZG},(n)} \rightarrow C_{\kappa\alpha p,\kappa'\alpha' p'}^{(n+1)}(T) \rightarrow \left\{ \omega_{\mathbf{q}\nu}(T), e_{\kappa\alpha,\nu}(\mathbf{q},T) \right\}^{(n+1)},\]

   and is repeated until self-consistency is reached.

8. Obtain the final anharmonic phonons and ZG configuration

   After convergence, the final temperature-dependent phonon frequencies and eigenvectors can be used to calculate anharmonic phonon dispersions and other vibrational properties at finite temperature. At the same time, the converged ZG displacements define an anharmonic ZG configuration representative of the nuclear thermal fluctuations of the self-consistent anharmonic system. This configuration can subsequently be used to evaluate temperature-dependent electronic and optical properties including anharmonic electron–phonon coupling effects [42] , [43].

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.

[31] M. Zacharias and F. Giustino, One-shot calculation of temperature-dependent optical spectra and phonon-induced band-gap renormalization, Phys. Rev. B 94, 075125 (2016), doi:10.1103/PhysRevB.94.075125.

[32] M. Zacharias, C. E. Patrick, and F. Giustino, Stochastic approach to phonon-assisted optical absorption, Phys. Rev. Lett. 115, 177401 (2015), doi:10.1103/PhysRevLett.115.177401.

[33] M. Zacharias and F. Giustino, Theory of the special displacement method for electronic structure calculations at finite temperature, Phys. Rev. Research 2, 013357 (2020), doi:10.1103/PhysRevResearch.2.013357.

[34] M. Zacharias, H. Seiler, F. Caruso, D. Zahn, F. Giustino, P. C. Kelires, and R. Ernstorfer, Multiphonon diffuse scattering in solids from first principles: Application to layered crystals and two-dimensional materials, Phys. Rev. B 104, 205109 (2021), doi:10.1103/PhysRevB.104.205109.

[35] D. Hooton, LI. A new treatment of anharmonicity in lattice thermodynamics: I, London Edinburgh Philos. Mag. J. Sci. 46, 422 (1955), doi:10.1080/14786440408520575.

[36] P. Souvatzis, O. Eriksson, M. I. Katsnelson, and S. P. Rudin, Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory, Phys. Rev. Lett. 100, 095901 (2008), doi:10.1103/PhysRevLett.100.095901.

[37] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011), doi:10.1103/10.1103/PhysRevB.84.180301.

[38] I. Errea, M. Calandra, and F. Mauri, First-principles theory of anharmonicity and the inverse isotope effect in superconducting palladium-hydride compounds, Phys. Rev. Lett. 111, 177002 (2013), doi:10.1103/PhysRevLett.111.177002.

[39] T. Tadano and S. Tsuneyuki, Self-consistent phonon calculations of lattice dynamical properties in cubic SrTiO3 with first-principles anharmonic force constants, Phys.Rev. B 92, 054301 (2015), doi:10.1103/PhysRevB.92.054301.

[40] R. Bianco, I. Errea, L. Paulatto, M. Calandra and F. Mauri, Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation, Phys. Rev. B 96, 014111 (2017), doi:10.1103/PhysRevB.96.014111.

[41] M. Zacharias, G. Volonakis, F. Giustino, and J. Even, Anharmonic lattice dynamics via the special displacement method, Phys. Rev. B 108, 035155 (2023), doi:10.1103/PhysRevB.108.035155.

[42] M. Zacharias, G. Volonakis, F. Giustino, and J. Even, Anharmonic electron-phonon coupling in ultrasoft and locally disordered perovskites, npj Comput. Mater. 9, 153 (2023), doi:10.1038/s41524-023-01089-2.

[43] M. Zacharias, G. Volonakis, L. Pedesseau, C. Katan, F. Giustino, and J. Even, Electron-phonon couplings in polymorphous crystals, Phys. Rev. B 113, L081104 (2026), doi:10.1103/n52n-g9nr.