Ab initio self-consistent many-body theory of polarons at all couplings

We present a theoretical framework to describe polarons from first principles within a many-body Green's function formalism. Starting from a general electron-phonon Hamiltonian, we derive a self-consistent Dyson equation in which the phonon-mediated self-energy is composed by two distinct terms. One term is the Fan-Migdal self-energy and describes dynamic electron-phonon processes, the other term is a contribution to the self-energy originating from the static displacements of the atomic nuclei in the polaronic ground state. The lowest-order approximation to the present theory yields the standard many-body perturbation theory approach to electron-phonon interactions in the limit of large polarons, and the ab initio polaron equations introduced [Sio et al., Phys. Rev. B 99, 235139 (2019); Phys. Rev. Lett. 122, 246403 (2019)] in the limit of small polarons. A practical recipe to implement the present unifying formalism in first-principles calculations is outlined. We apply our method to the Fröhlich model, and obtain remarkably accurate polaron energies at all couplings, in line with Feynman's polaron theory and diagrammatic Monte Carlo calculations. We also recover the well-known results of Fröhlich and Pekar at weak and strong coupling, respectively. The present approach enables predictive many-body calculations of polarons in real materials at all couplings.

Figure 1

Diagrammatic representation of the self-consistent Green's function theory of polarons. Legends for the different parts of the diagrams are given in the lower left corner. (a) Fan-Migdal self-energy in Eq. (18). (b) Self-consistent definition of the vertex function in Eq. (24). (c) Polaronic self-energy in Eq. (33). (d) Dyson equation (17). (e) Schematic representation of the self-consistent solution of Eqs. (17), (18), (24), and (33).

Figure 2

Schematic representation of the different energies involved in the polaron formation process. The solid gray line represents the original periodic $N$-electron system, and the solid blue line represents the ($N+1$)-electron system for which a distorted polaron configuration is the ground state in the presence of a pinning potential (cf. Sec. 6). The dashed gray line represents a fictitious ($N+1$)-electron system in which an electron with no phonon-mediated renormalization has been added to the conduction band minimum (${\varepsilon }_{\mathrm{CBM}}^{\mathrm{QP}}$). The excitation energy from the ($N+1$)-electron ground state to the distorted $N$-electron state is the electron addition or removal energy ${\varepsilon }_s$. The energy released by the relaxation of the lattice back to the periodic configuration is represented by $E_{\mathrm{lattice}}$. The polaron formation energy, that is the energy gained by the system when a delocalized electronic state becomes localized in a polaronic state, is represented by $\mathrm{\Delta }E$.

Figure 3

Schematic representation of the self-consistent procedure required to solve the many-body polaron equations (48) and (49).

Figure 4

Variational ansatz for the electronic component of the polaron wave function in the Fröhlich model. (a) Wave-function plot for a few values of the variational parameter $r_p$, which corresponds to the polaron radius, from Eq. (72). $a_0$ is the Bohr radius. (b) Reciprocal-space coefficients of the wave functions shown in (a), from Eq. (73).

Figure 5

Ground-state polaron energy as a function of the polaron radius $r_p$, within the Fröhlich model for $\alpha =4.94$. (a) $r_p$ dependence of the different terms contributing to the total energy in Eq. (78), namely, the kinetic energy (dotted), the Coulomb energy composed by the Landau-Pekar and the lattice energy (dotted-dashed), and the Fan-Migdal self-energy contribution (dashed). (b) Total energy as a function of the polaron radius $r_p$. The red solid line represents the Landau-Pekar result [35], where only the kinetic and the Coulomb energies are considered. The blue solid line represents the total energy including the FM contribution, as in Eq. (78). The energy minima are highlighted by the solid crosses.

Figure 6

Total energy of the ground state of the Fröhlich polaron as a function of the coupling strength $\alpha$. (a) Fan-Migdal solution in Rayleigh-Schrödinger perturbation theory (FM-RS, blue line), Landau-Pekar solution (LP, red line), Feynman's variational path-integral solution (gray line) [26, 91, 92], and diagrammatic Monte Carlo results (DMC, black circles) [28, 29, 93]. (b) Relative deviation between FM-RS and LP energies with respect to Feynman's result, following the same color scheme as in (a). (c) Ground-state energy of the Fröhlich polaron evaluated using the present Green's function approach. The light green line represents the self-consistent solution, the dark-green line represents the perturbative calculation. The relative errors of each approximation with respect to Feynman's result are shown in (d).

Figure 7

Comparison between the hole polaron wave function in LiF obtained by (a) direct DFT (HSE) supercell calculation with a removed electron and (b) the solution of the polaron equations in Eqs. (59) and (62).

Figure 8

Analysis of the various contributions to the polaron formation energy in the lowest-order approximation, Eqs. (58), (62), and (63). (a) Polaron eigenvalue (light blue triangles) and lattice energy (dark blue crosses) as a function of the supercell size for the electron polaron in LiF. The supercell size is given as $L^{-1}$, where $L^3$ is the supercell volume. The numbers next to the data points indicate the number of unit cells in each direction of the homogeneous Born–von Karman supercell or, equivalently, the number of $\mathbf{k}$ points in each direction of the homogeneous mesh of the Brillouin zone. The dashed lines correspond to the extrapolation of each contribution to the infinite supercell size. (b) Same as (a) but for the hole polaron in LiF. (c) Expectation value of the $\mathrm{FM}+\mathrm{DW}$ self-energy over the polaron quasiparticle amplitudes for the electron polaron in LiF, as a function of the number of $\mathbf{k}$ points in each direction of the homogeneous mesh of the Brillouin zone. (d) Same as (c) but for the hole polaron in LiF.