Electron-phonon interactions theory

Electron-phonon interactions from perturbation theory

This theory section is meant to outline the most important equations and their physical meaning in the context of electron-phonon calculations using VASP. One can obtain these quantities rigorously from many-body perturbation theory and then apply approximations to make the calculations feasible at the DFT level. For a comprehensive theoretical description of electron-phonon interactions using many-body perturbation theory, we recommend to read Ref. or other contemporary literature.

Electron self-energy

The interaction between electrons and lattice vibrations (phonons) leads to a modification of the electronic energies, described by the electron self-energy, [math]\displaystyle{ \Sigma_{n\mathbf{k}}(E, T) }[/math], for the band [math]\displaystyle{ n }[/math] and wavevector [math]\displaystyle{ \mathbf{k} }[/math] at energy [math]\displaystyle{ E }[/math] and temperature [math]\displaystyle{ T }[/math] (often, the energy dependence is written as a frequency dependence, [math]\displaystyle{ \hbar \omega }[/math]). This complex quantity has a real part, [math]\displaystyle{ \text{Re}[\Sigma_{n\mathbf{k}}(E, T)] }[/math], which represents the shift or "renormalization" of the electron's energy. The imaginary part, [math]\displaystyle{ \text{Im}[\Sigma_{n\mathbf{k}}(E, T)] }[/math], describes the finite lifetime of the electronic state due to scattering with phonons.

Using second-order perturbation theory, the electron self-energy can be separated into two contributions: the Fan-Migdal (FM) and Debye-Waller (DW) terms:

[math]\displaystyle{ \Sigma_{n\mathbf{k}}(E, T) = \Sigma^{\text{FM}}_{n\mathbf{k}}(E, T) + \Sigma^{\text{DW}}_{n\mathbf{k}}(E, T). }[/math]

In practice, we often use the on-the-mass-shell approximation, where the self-energy is evaluated at the Kohn-Sham (KS) eigenvalue, [math]\displaystyle{ E = \varepsilon_{n\mathbf{k}} }[/math]. Let us highlight the FM contribution, which involves the electron-phonon matrix elements [math]\displaystyle{ g_{mn\mathbf{k}, \nu \mathbf{q}} }[/math]:

[math]\displaystyle{ \Sigma^{\text{FM}}_{n\mathbf{k}}(T) = \frac{1}{N_q} \sum_{m\nu\mathbf{q}} |g_{mn\mathbf{k}, \nu \mathbf{q}}|^2 \left[ \frac{n_{\nu\mathbf{q}}(T) + 1 - f_{m\mathbf{k}+\mathbf{q}}(T)}{\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k}+\mathbf{q}} - \omega_{\nu\mathbf{q}} + i\delta} + \frac{n_{\nu\mathbf{q}}(T) + f_{m\mathbf{k}+\mathbf{q}}(T)}{\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k}+\mathbf{q}} + \omega_{\nu\mathbf{q}} + i\delta} \right], }[/math]

where [math]\displaystyle{ n_{\nu\mathbf{q}}(T) }[/math] is the Bose-Einstein occupation of phonon mode [math]\displaystyle{ (\nu, \mathbf{q}) }[/math] with index [math]\displaystyle{ \nu }[/math] and phonon wavevector [math]\displaystyle{ \mathbf{q} }[/math], [math]\displaystyle{ \omega_{\nu\mathbf{q}} }[/math] is the phonon frequency, [math]\displaystyle{ N_q }[/math] is the number of [math]\displaystyle{ \mathbf{q} }[/math]-points, and [math]\displaystyle{ f_{m\mathbf{k}+\mathbf{q}}(T) }[/math] is the Fermi-Dirac occupation of the electronic state with band [math]\displaystyle{ m }[/math] and wavevector [math]\displaystyle{ \mathbf{k}+\mathbf{q} }[/math]. The infinitesimal imaginary shift [math]\displaystyle{ i \delta }[/math] in the denominator ensures the correct pole structure of the self-energy. In practice, [math]\displaystyle{ \delta }[/math] is often finite and used as a smearing or convergence parameter.

The DW contribution is a purely real shift and can be calculated using the phonon frequencies and eigenvectors, as well as the second derivatives of the KS potential with respect to atomic displacements:

[math]\displaystyle{ \Sigma^{\text{DW}}_{n\mathbf{k}}(T) = \frac{1}{N_q} \sum_{\nu\mathbf{q}} g^{2,\text{DW}}_{n \mathbf{k}, \nu\mathbf{q}} (2n_{\nu\mathbf{q}}(T) + 1), }[/math]

where [math]\displaystyle{ g^{2,\text{DW}}_{n \mathbf{k}, \nu\mathbf{q}} }[/math] represents the second-order change in the KS potential due to phonon mode [math]\displaystyle{ (\nu, \mathbf{q}) }[/math]. To evaluate this term in practice, we employ the rigid-ion approximation, which simplifies the calculation by assuming that the potential changes rigidly with atomic displacements. This approximation is implemented in VASP. It means that we do not need to compute the second derivatives of the KS potential explicitly, but can instead use the first derivatives (i.e., the electron-phonon matrix elements) to estimate the DW term.

The combination of [math]\displaystyle{ \Sigma^{\text{FM}}_{n\mathbf{k}}(T) + \Sigma^{\text{DW}}_{n\mathbf{k}}(T) }[/math] at this level of approximation is often called the Allen-Heine-Cardona (AHC) theory. To be precise, the expressions we showed here are have come to be know as nonadiabatic AHC theory, owing to the fact that the phonon frequency [math]\displaystyle{ \omega_{\nu \mathbf{q}} }[/math] appears in the denominator of the FM self-energy. Physically, this means that energy can be exchanged between the electrons and phonons. In contrast, there is also an adiabatic formula where the phonon frequency is neglected in the denominator, which assumes that the electrons respond instantaneously to the slower lattice vibrations:

[math]\displaystyle{ \Sigma^{\text{FM, ad}}_{n\mathbf{k}}(T) = \frac{1}{N_q} \sum_{m\nu\mathbf{q}} |g_{mn\mathbf{k}, \nu \mathbf{q}}|^2 \frac{2 n_{\nu\mathbf{q}}(T) + 1}{\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k}+\mathbf{q}} + i\delta}. }[/math]

The adiabatic approximation is sometimes used for calculating band-structure renormalizations in semiconductors and insulators, where the phonon frequencies are typically much smaller than the electronic energy differences. However, the adiabatic approximation breaks down in polar materials due to the long-range Fröhlich interaction, which can lead to divergences in the self-energy.

To compute both the FM and DW contributions, we need the electron-phonon matrix elements. In the context of KS DFT, the matrix elements can be defined using the KS orbitals, [math]\displaystyle{ \psi_{n\mathbf{k}} }[/math], and the perturbation of the KS potential, [math]\displaystyle{ V_{\text{KS}} }[/math], due to phonon displacements:

[math]\displaystyle{ g_{mn\mathbf{k}, \nu \mathbf{q}} = \left\langle \psi_{m\mathbf{k}+\mathbf{q}} \left| \partial_{\nu \mathbf{q}} V_{\text{KS}} \right| \psi_{n\mathbf{k}} \right\rangle, }[/math]

We usually call the object [math]\displaystyle{ \partial_{\nu \mathbf{q}} V_{\text{KS}} }[/math] the electron-phonon potential. The phonon displacement operator [math]\displaystyle{ \partial_{\nu \mathbf{q}} }[/math] can be expressed as a superposition of atomic displacements:

[math]\displaystyle{ \partial_{\nu \mathbf{q}} = \sum_{l\kappa\alpha} \sqrt{\frac{\hbar}{2 m_\kappa \omega_{\nu\mathbf{q}}}} e_{\kappa\alpha,\nu \mathbf{q}} \mathrm{e}^{\mathrm{i} \mathbf{q}\cdot\mathbf{R}_l} \frac{\partial}{\partial u_{l\kappa\alpha}}, }[/math]

where [math]\displaystyle{ \kappa }[/math] and [math]\displaystyle{ \alpha }[/math] index the atoms in the unit cell and Cartesian directions, respectively, [math]\displaystyle{ m_\kappa }[/math] is the mass of atom [math]\displaystyle{ \kappa }[/math], [math]\displaystyle{ e_{\kappa\alpha,\nu \mathbf{q}} }[/math] is the phonon eigenvector component, and [math]\displaystyle{ u_{l\kappa\alpha} }[/math] is the displacement of atom [math]\displaystyle{ \kappa }[/math] in cell [math]\displaystyle{ \mathbf{R}_l }[/math].

There is a set of instructions on how to compute the electron-phonon potential.

Band-structure renormalization

As mentioned above, the electron-phonon interactions lead to a renormalization of the electronic band structure. The renormalized electronic energies, [math]\displaystyle{ E_{n\mathbf{k}}(T) }[/math], can be obtained by adding the real part of the self-energy to the KS eigenvalues:

[math]\displaystyle{ E_{n\mathbf{k}}(T) = \varepsilon_{n\mathbf{k}} + \text{Re}[\Sigma_{n\mathbf{k}}(T)], }[/math]

where [math]\displaystyle{ \varepsilon_{n\mathbf{k}} }[/math] are the Kohn-Sham eigenvalues.

Of particular interest is the renormalization of the bandgap in semiconductors and insulators, which refers to the change in the energy difference between the conduction band minimum (CBM) and the valence band maximum (VBM) due to these electron-phonon interactions. In this case, we evaluate [math]\displaystyle{ \Sigma_{n\mathbf{k}}(T) }[/math] at the KS eigenvalues that correspond to the CBM and VBM to get the bandgap renormalization. The renormalization varies with temperature as phonon populations change. However, it also occurs at 0 K due to zero-point atomic motion where it is referred to as zero-point renormalization (ZPR).

For the computation of the band-structure renormalization, it is important to include both the FM and DW contributions to the self-energy. Neglecting the DW term can lead to significant errors in the predicted bandgap renormalization, as the DW term often partially cancels the FM contribution.

For information on how to compute the band-structure renormalization using VASP, please refer to this how-to page of the manual.

Electron-phonon interactions from statistical sampling

In contrast to the perturbative approach outlined above, this section shows how to describe electron-phonon interactions from a statistical point of view. To begin with, the probability distribution of finding an atom within the coordinates [math]\displaystyle{ \kappa+d\kappa }[/math] (where [math]\displaystyle{ \kappa }[/math] denotes the Cartesian coordinates as well as the atom number) at temperature [math]\displaystyle{ T }[/math] in the harmonic approximation is given by the following expression

[math]\displaystyle{ dW_{\nu}(\kappa,T)=\frac{1}{2\pi \langle u^{2}_{\nu \kappa}\rangle} e^{-\kappa^{2}/(2 \langle u^{2}_{\nu \kappa}\rangle)} d\kappa, }[/math]

where the mean-square displacement of the harmonic oscillator is given as

[math]\displaystyle{ \langle u^{2}_{\nu \kappa}\rangle = \frac{\hbar}{2 M_{\kappa} \omega_{\nu}} \coth{\frac{\hbar \omega_{\nu}}{2 k_{B}T}}. }[/math]

Here [math]\displaystyle{ M_{\kappa} }[/math], [math]\displaystyle{ \nu }[/math] and [math]\displaystyle{ \omega_{\nu} }[/math] denote the mass, phonon eigenmode and phonon eigenfrequency, respectively. The equation for [math]\displaystyle{ dW }[/math] is valid at any temperature and the high (Maxwell--Boltzmann distribution) and low temperature limits are easily regained. In order to obtain an observable [math]\displaystyle{ O(T) }[/math] at a given temperature [math]\displaystyle{ T }[/math], the average of the observable sampled at different coordinate sets [math]\displaystyle{ x_{T}^{\textrm{MC},i} }[/math] with sample size [math]\displaystyle{ n }[/math] is taken

[math]\displaystyle{ \langle O(T)\rangle = \frac{1}{n} \sum\limits_{i=1}^{n} O(x_{T}^{\textrm{MC,i}}). }[/math]

Full Monte-Carlo sampling

Each set [math]\displaystyle{ i }[/math] is obtained from the equilibrium atomic positions [math]\displaystyle{ x_{\textrm{eq}} }[/math] as

[math]\displaystyle{ x_{T}^{\textrm{MC,i}} = x_{\textrm{eq}} + \Delta \tau^{\textrm{MC,i}} }[/math]

with the displacement

[math]\displaystyle{ \Delta \tau^{\textrm{MC,i}} = \sqrt{\frac{1}{M_{\kappa}}} \sum\limits_{\nu}^{3(N-1)} \varepsilon_{\kappa,\nu} \mathcal{N}. }[/math]

Here [math]\displaystyle{ \varepsilon_{\kappa,\nu} }[/math] denotes the unit vector of eigenmode [math]\displaystyle{ \nu }[/math] on atom [math]\displaystyle{ \kappa }[/math]. The magnitude of the displacement in each Cartesian direction is obtained from the normal-distributed random variable [math]\displaystyle{ \mathcal{N} }[/math] with a probability distribution according to [math]\displaystyle{ dW }[/math].

ZG configuration (one-shot method)

Motivated by the empirical observation that for increasing super-cell sizes the number of required structures in the MC method can be decreased, M. Zacharias and F. Giustino proposed a one-shot method where only a single set of displacements is used

[math]\displaystyle{ \Delta \tau^{\textrm{OS}} = \sqrt{\frac{1}{M_{\kappa}}} \sum\limits_{\nu}^{3(N-1)} (-1)^{\nu-1} \varepsilon_{\kappa,\nu} \sigma_{\nu,T}, }[/math]

where the summation over the eigenmodes runs in an ascending order with respect to the values of the eigenfrequencies, and the magnitude of each displacement is given by

[math]\displaystyle{ \sigma_{\nu,T} = \sqrt{(2 n_{\nu,T}+1) \frac{\hbar}{2 \omega_{\nu}}}. }[/math]

Here [math]\displaystyle{ n_{\nu,T}=[\mathrm{exp}(\hbar \omega_{\nu} /k_{B}T)-1]^{-1} }[/math] denotes the Bose-Einstein occupation number. In this way, the sum for the observable [math]\displaystyle{ \langle O(T)\rangle }[/math] is reduced to a single calculation. In Ref. it was shown that for super-cell sizes [math]\displaystyle{ N\rightarrow \infty }[/math] the structural configuration obtained using the ZG configuration should lead to equivalent results as fully converged MC calculations. In practice, it was shown that already relatively small supercell sizes are sufficient to achieve good accuracy, but the convergence with respect to the cell size can vary between different materials. In Ref. we have also used a slightly modified approach, in which the signs of the displacements are chosen randomly instead of [math]\displaystyle{ \pm 1 }[/math]. This was only necessary when calculating volume-dependent ZPR, since the modes sometimes change the order as the volume changes. Using alternating signs for the displacement then causes small discontinuities in the ZPR volume curve of the order of 5 meV for carbon diamond. By averaging over many random phases, this problem can be eliminated. Nevertheless 5 meV difference between the ZG configuration and full MC is considered as accurate enough in most calculations.

For information on how to compute the band-structure renormalization using VASP, please refer to this how-to page of the manual.

---