Polarons from first principles, without supercells

TL;DR

This paper introduces a first-principles computational method to study polarons in insulators and semiconductors without the need for large supercells, enabling accurate analysis of polaron properties.

Contribution

The authors develop a novel formalism that uses density-functional perturbation theory to study polarons, avoiding supercell limitations and unifying the description of large and small polarons.

Findings

Successfully calculated polaron wavefunctions in LiF and Li2O2.

Accurately determined polaron formation energies.

Demonstrated seamless description of large and small polarons.

Abstract

We develop a formalism and a computational method to study polarons in insulators and semi-conductors from first principles. Unlike in standard calculations requiring large supercells, we solve a secular equation involving phonons and electron-phonon matrix elements from density-functional perturbation theory, in a spirit similar to the Bethe-Salpeter equation for excitons. We show that our approach describes seamlessly large and small polarons, and we illustrate its capability by calculating wavefunctions, formation energies, and spectral decomposition of polarons in LiF and Li2O2.

Equations

$$E=\int d\mathbf{r}{\psi }^{\ast }{\hat{H}}_{\mathrm{KS}}^0\psi +\int d\mathbf{r}\frac{\partial V_{\mathrm{KS}}^0}{\partial {\tau }_s}|\psi {|}^2\Delta {\tau }_s+\frac{1}{2}{C}_{s{s}^{\prime }}^0\Delta {\tau }_s\Delta {\tau }_{s^{\prime }}.$$

$${\hat{H}}_{\mathrm{KS}}^0\psi +\frac{\partial V_{\mathrm{KS}}^0}{\partial {\tau }_s}\psi \Delta {\tau }_s=\epsilon \psi ,$$

$$\Delta {\tau }_s=-(C^0{)}_{s{s}^{\prime }}^{-1}\int d\mathbf{r}\frac{\partial V_{\mathrm{KS}}^0}{\partial {\tau }_{s^{\prime }}}|\psi {|}^2,$$

$$\frac{2}{N_p}\sum _{\mathbf{q}m\nu }{B}_{\mathbf{q}\nu }{g}_{mn\nu }^{\ast }(\mathbf{k},\mathbf{q})A_{m\mathbf{k}+\mathbf{q}}=({\epsilon }_{n\mathbf{k}}-\epsilon )A_{n\mathbf{k}},$$

$$B_{\mathbf{q}\nu }=\frac{1}{N_p}\sum _{mn\mathbf{k}}{A}_{m\mathbf{k}+\mathbf{q}}^{\ast }\frac{g_{mn\nu }(\mathbf{k},\mathbf{q})}{\hslash {\omega }_{\mathbf{q}\nu }}{A}_{n\mathbf{k}},$$

$$\Delta E_f=\frac{1}{N_p}\sum _{n\mathbf{k}}|A_{n\mathbf{k}}{|}^2({\epsilon }_{n\mathbf{k}}-{\epsilon }_{\mathrm{CBM}})-\frac{1}{N_p}\sum _{\mathbf{q}\nu }|B_{\mathbf{q}\nu }{|}^2\hslash {\omega }_{\mathbf{q}\nu },$$

Full text

Polarons have been attracting unrelenting attention ever since the polaron concept was formulated by Landau a century ago (Landau, 1933). For example polarons inspired the search for high-temperature superconducting oxides (Bednorz and Müller, 1986), they are considered one of the hallmarks of emergent behavior in quantum matter (Verdi et al., 2017; Wan, 2016; Cancellieri et al., 2016; Chen et al., 2015; Nie et al., 2015), and they have been linked to the extraordinary defect-tolerance of the new metal-halide perovskites (Zhu et al., 2016). At a more fundamental level, the quest for a satisfactory quantum-mechanical description of polarons stimulated much progress in diverse areas of theoretical physics. For example the solution of the Fröhlich polaron problem by Feynman was a landmark in the development of the path integral formulation of quantum mechanics (Feynman, 1955), and the Pekar polaron problem (Pekar, 1946) found applications in general relativity (Bahrami et al., 2014) and quantum state reduction (Penrose, 2014).

In the simplest picture, a self-trapped polaron forms when an excess electron or hole deforms a crystal lattice so as to create a potential well from which it cannot escape. Microscopic models of this effect have been developed and investigated in many seminal contributions of the last century (Pekar, 1946; Fröhlich et al., 1950; Lee et al., 1953; Feynman, 1955; Holstein, 1959), and subsequently refined to address increasingly more realistic scenarios, such as multiple electron bands, dispersive phonons, and transport properties (Lépine and Frongillo, 1992; Lépine, 1994; Alexandrov and Kornilovitch, 1999; Alexandrov, 2000; Alexandrov and Yavidov, 2004; Perroni et al., 2004; Rösch and Gunnarsson, 2004; Cruzeiro-Hansson et al., 2000; Alexandrov, 2003; Ortmann et al., 2009). More recently, significant progress in the theory of polarons has been achieved with the development of numerical many-body techniques, such as exact diagonalization (Fehske and Trugman, 2007), renormalization group (Jeckelmann and White, 1998; Grusdt, 2016), continuous-time quantum Monte Carlo (Kornilovitch, 2007), and diagrammatic Monte Carlo methods (Mishchenko et al., 2000; Prokof’ev and Svistunov, 2008; Hahn et al., 2018). Comprehensive reviews of the field can be found in Refs. Alexandrov, 2007; Devreese and Alexandrov, 2009; Alexandrov and Devreese, 2010; Emin, 2013; Devreese, .

The common denominator to most theoretical studies on polarons is that they focus on idealized mathematical models, for example the Fröhlich Hamiltonian (Fröhlich et al., 1950) and the Holstein Hamiltonian (Holstein, 1959), which describe a free or a tightly-bound electron interacting with a dispersionless optical phonon, respectively. These models offer an ideal testbed for methodological development, and shaped our current understanding of polarons. However, they are not suitable for studying real materials, as they lack essential features such as complex unit cells, band structures, phonon dispersion relations, and realistic electron-phonon coupling matrix elements. Furthermore, such models are not transferable to complex systems such as surfaces, interfaces, low-dimensional materials and heterostructures. Therefore, there is a need for supplementing correlated methods for polarons with more realistic materials parameters, as emphasized by authoritative reviews (Devreese and Alexandrov, 2009).

At the other end of the spectrum, ab initio calculations based on density functional theory (DFT) are ideally positioned to address the complexity of real materials. Indeed, studies of polarons under realistic conditions have begun to emerge during the past decade (Franchini et al., 2009; Setvin et al., 2014; Himmetoglu et al., 2014; Bondarenko et al., 2015; Reticcioli et al., 2017, 2018; Morbec and Galli, 2016; Kokott et al., 2018). However, also DFT faces important limitations: the calculations necessitate large supercells (Spreafico and VandeVondele, 2014; Erhart et al., 2014; Kokott et al., 2018), hence they are prohibitive for intermediate and large polarons which require several thousand atoms; local exchange-correlation functionals tend to suppress polaron self-trapping; calculations using Hubbard-corrected or hybrid functionals suffer from the sensitivity to the Hubbard parameter or the exchange fraction (Kokott et al., 2018). More fundamentally, the relation between DFT calculations of polarons and the vast literature on model Hamiltonians remains unclear.

In the present work we wish to overcome these limitations by filling the gap between model Hamiltonians and atomistic calculations of polarons. To this aim, we reformulate the direct calculation of polarons with DFT into a nonlinear eigenvalue problem. The ingredients of this formalism are the band structures, phonons, and electron-phonon matrix elements calculated in the crystal unit cell from density functional perturbation theory. The solution of this eigenvalue problem yields the formation energy of the polaron, its excitation energy, the electronic wavefunctions and atomic displacements, as well as the spectral decomposition of the polaron in terms of the underlying electron-phonon coupling mechanisms. We validate this methodology by studying two limiting cases, the large electron polaron in LiF and the small electron polaron in Li2O2. Complete derivations and extensive benchmarks are reported in the companion manuscript, Ref. Sio et al., .

We start by considering the DFT total energy of an excess electron added to a crystal with a finite band gap. The same reasoning holds for holes (Sio et al., ). The ground state is spin-unpolarized, and we make the working assumption that the perturbation to the valence Kohn-Sham (KS) states due to the extra electron can be neglected (to be validated a posteriori). By expanding the total energy in powers of the atomic displacements from their equilibrium positions in the ground state, at the lowest order which admits non-trivial solutions we obtain:

[TABLE]

Here the energy $E$ is referred to the ground state, $\psi({\bf r})$ is the wavefunction of the excess electron and $\Delta\tau_{s}$ are the atomic displacements; $s=(\kappa{\alpha}p)$ is a composite index denoting the Cartesian coordinate ${\alpha}$ of atom $\kappa$ in the unit cell $p$, and the Einstein summation convention is implied. The integrals are over a suitably large Born-von Kárman supercell. ${\hat{H}}_{\rm KS}^{0}$, $\partial V_{\rm KS}^{0}/\partial\tau_{s}$, and $C^{0}_{ss^{\prime}}$ represent the KS Hamiltonian, the variation of the KS potential resulting from an atomic displacement, and the matrix of interatomic force constants, respectively. The superscript ‘0’ indicates that these quantities are evaluated in the ground state, without excess electron. In Eq. (1) the spurious Hartree and exchange-correlation interactions of the polaron with itself and its periodic images are carefully eliminated via a self-interaction correction, as discussed in Ref. Sio et al., .

The total energy $E$ in Eq. (1) can be regarded as a functional of $\psi$ and $\Delta\tau_{s}$. By minimizing this functional subject to the constraint that $\psi$ be normalized, we obtain the nonlinear eigenvalue problem:

[TABLE]

where $\varepsilon$ is the polaron eigenvalue. In principle these equations could be solved in real space, but in most practical applications this is prohibitive, since the supercell must be large enough to accommodate the wavefunction $\psi$. To overcome this obstacle we proceed as in the calculation of excitons via the Bethe-Salpeter equation (Onida et al., 1995; Rohlfing and Louie, 1998), that is we expand the solution in terms of unperturbed KS states and phonons eigenmodes. To this aim we define $\psi=N_{p}^{-1/2}\sum_{n{\bf k}}A_{n{\bf k}}\,\psi_{n{\bf k}}$ and $\Delta\tau_{s}=-2N_{p}^{-1}\sum_{{\bf q}\nu}B^{*}_{{\bf q}\nu}\,(\hbar/2M_{\kappa}\omega_{{\bf q}\nu})^{1/2}e_{\kappa{\alpha},{\bf q}\nu}\,e^{i{\bf q}\cdot{\bf R}_{p}}$. Here $N_{p}$ is the number of unit cells in the supercell, $\psi_{n{\bf k}}$ is an unoccupied eigenstate of ${\hat{H}}_{\rm KS}^{0}$ for the band $n$ and wavevector ${\bf k}$ with energy $\varepsilon_{n{\bf k}}$, and $e_{\kappa{\alpha},{\bf q}\nu}$ is the vibrational mode with branch $\nu$, wavevector ${\bf q}$, and frequency $\omega_{{\bf q}\nu}$, obtained by diagonalizing $C^{0}_{ss^{\prime}}$. ${\bf R}_{p}$ is a vector of the direct lattice, and $M_{\kappa}$ is the mass of atom $\kappa$. Using these definitions in Eqs. (2)-(3) we obtain a nonlinear eigenvalue problem for the generalized Fourier amplitudes $A_{n{\bf k}}$ and $B_{{\bf q}\nu}$:

[TABLE]

where $g_{mn\nu}({\bf k},{\bf q})$ is the electron-phonon matrix element for the scattering between the electronic states $|n{\bf k}\rangle$ and $|m{\bf k}+{\bf q}\rangle$ via the phonon ${\bf q}\nu$ (Giustino, 2017). Equations (4)-(5) only require the band structures, phonon dispersions, and matrix elements calculated in the unit cell. In this representation the formation energy $\Delta E_{f}$ of the polaron, that is the total energy of the self-trapped polaron minus the energy of the undistorted crystal with an extra electron at the conduction band bottom, reads (Sio et al., ):

[TABLE]

where the KS eigenvalue is referred to the conduction band minimum. This expression shows that the polaron formation energy consists of a positive-definite electronic contribution and a negative-definite vibrational contribution. This spectral representation allows us to identify $|B_{{\bf q}\nu}|^{2}$ with the number of phonons taking part in the polaron, as discussed in greater detail in Ref. Sio et al., .

We now illustrate this approach using LiF and Li2O2 as case studies. LiF is a simple salt that crystallizes in the rocksalt structure. It is a paradigmatic wide gap polar insulator, and is known to host large electron polarons (Iadonisi, 1984). Li2O2 is a prototypical cathode for lithium-air batteries, and hosts small electron polarons (Feng et al., 2013). The structure consists of two-dimensional LiO2 layers intercalated by Li planes. We perform calculations within the PBE generalized gradient approximation to DFT (Perdew et al., 1996), using planewaves and pseudopotentials as implemented in the Quantum Espresso suite (Giannozzi et al., 2017). We employ optimized norm-conserving Vanderbilt pseudopotentials (Hamann, 2013) and planewaves kinetic energy cutoffs of 150 Ry and 105 Ry for LiF and Li2O2, respectively. We calculate phonons and electron-phonon matrix elements within density functional perturbation theory (Baroni et al., 2001; Giustino, 2017), and we perform Wannier interpolation of all properties using the Wannier90 (Marzari et al., 2012; Mostofi et al., 2014) and EPW (Giustino et al., 2007; Poncé et al., 2016) codes. We employ uniform Brillouin-zone grids in Eqs. (4)-(5), and we solve the linear system via a parallel steepest descent algorithm. As the initial seed for $A_{n{\bf k}}$ in Eq. (5) we use a simple Gaussian lineshape (Sio et al., ).

Figure 1 shows our results for LiF. In Fig. 1(a) we see that the electron wavefunction of the polaron extends over $\sim$10 unit cells, therefore we are in the presence of a large polaron. A cross-sectional view in the [010] direction of the same wavefunction is shown in Fig. 1(b). Here we recognize an envelope function of approximately Gaussian shape, which modulates the atomic Li and F $2s$ orbitals. From this plot we can quantify the spatial extent of the polaron using the full-width at half maximum, $2r_{\rm p}=9.0$ Å. This value is consistent with an earlier semiempirical estimate of 9.3 Å (Iadonisi, 1984). In Fig. 1(c) we show the atomic displacements along a line that cuts near the center of the polaron. As expected, also the displacements follow an approximately Gaussian profile.

In Fig. 1(d)-(e) we analyze the composition of the polaron in terms of the Fourier amplitudes $A_{n{\bf k}}$ and $B_{{\bf q}\nu}$. Panel (d) shows that the electron wavefunction draws weight primarily from KS states at the bottom of the conduction band. This localization in reciprocal space is consistent with the highly delocalized nature of the polaron; analogous structures are observed in the related physics of Wannier excitons (Bokdam et al., 2016). In panel (e) we see the phonon eigenmodes participating in the polaron. There is a strong contribution from the longitudinal optical (LO) phonons at the zone center, around 77 meV; this is an indication of Fröhlich-type electron-phonon coupling. Our approach also reveals a non-negligible contribution from longitudinal acoustic (LA) phonons, up to 40 meV. By integrating $|B_{{\bf q}\nu}|^{2}$ across the Brillouin zone and summing over the phonon branches, we find that the electron polaron in LiF involves $\sim$5 LO phonons and $\sim$3 LA phonons, respectively.

Figure 1(f) shows the polaron eigenvalue $\varepsilon$ and formation energy $\Delta E_{f}$ as a function of supercell size. For supercells smaller than 12$\times$12$\times$12 unit cells the nonlinear eigenproblem in Eqs. (4)-(5) does not admit localized solutions. This can be understood as a manifestation of the Mott transition (Mott, 1968) at a critical density of $4\cdot 10^{19}$ cm*-3*. Below this critical density we find localized solutions of the type shown in Fig. 1(a), with an energy that scales as a constant plus a term proportional to $L^{-1}$, where $L$ is the supercell size. This trend is understood as the Madelung energy of a superlattice of polarons. If we extrapolate to $L\rightarrow\infty$ we obtain the energy of one isolated polaron (Madelung, 1918; Makov and Payne, 1995). In this dilute limit the polaron formation energy is $-0.23$ eV, and the polaron eigenvalue is $-0.80$ eV with respect to the conduction band bottom. The ratio between these values follows approximately the 1/3 scaling law that is expected for the Pekar polaron, which considers exclusively the Fröhlich coupling (Sio et al., ; Alexandrov and Devreese, 2010).

To validate our approach for LiF, we performed self-interaction corrected DFT calculations up to supercells of size 7$\times$7$\times$7 unit cells, containing up to 686 atoms (Sio et al., ). In agreement with the results of Eqs. (4)-(5), these direct DFT calculations did not yield any localized solutions. Calculations for supercells large enough to have localized solutions would be prohibitively expensive, as they involve $>3400$ atoms. Therefore, in order to validate our theory in the dilute limit, we follow an alternative route and we compare with the prediction of the continuum Pekar polaron model (Pekar, 1946). To this aim we repeat all calculations in Fig. 1(f) after replacing the band structure by a parabolic model with the DFT effective mass, the phonon dispersions by a dispersionless LO mode, and the electron-phonon matrix elements by the long-range Fröhlich interaction following Ref. Verdi and Giustino, 2015. In Fig. 1(g) we show that our theory reproduces exactly the energetics of the Pekar model in the continuum limit.

Now we move to Li2O2 in Fig. 2. In this case we find a small polaron, as seen from the electron wavefunction in Fig. 2(a). The polaron is as small as two adjacent O-$2p$ atomic orbitals [Fig. 2(g)]; from the cross-sectional view in Fig. 2(b) we deduce a size $2r_{\rm p}=$ 1.3 Å in the (100) plane. Correspondingly the atomic displacements are highly localized, and only the first shell of atoms around the polaron center exhibits non-negligible distortions [Fig. 2(c)]. The electronic and vibrational Fourier amplitudes $A_{n{\bf k}}$ and $B_{{\bf q}\nu}$ of the polarons look very different from the case of LiF in Fig. 1. In fact in Fig. 2(d) we see that all states of the lowest conduction bands contribute uniformly to the polaron wavefunction; similarly phonons from the entire Brillouin zone and from the non-polar and polar branches in Fig. 2(e) contribute to the atomic displacements. These signatures are reminiscent of Holstein-type electron-phonon coupling (Holstein, 1959), albeit with multiple electron bands and phonon branches involved. From the square amplitudes $|B_{{\bf q}\nu}|^{2}$ we infer that the small polaron in Li2O2 involves $\sim$13  polar optical phonons centered around 72 meV and $\sim$46  non-polar optical phonons at energies near 96 meV.

In Fig. 2(f) we show the polaron eigenvalue and formation energy as a function of supercell size. In this case we observe the formation of polarons already for 3$\times$3$\times$3 supercells, and the corresponding critical density for the Mott transition is in the range $\sim 10^{21}$ cm*-3* (Sio et al., ). The small electron polaron in Li2O2 is very energetic, exhibiting $\Delta E_{f}=-4.9$ eV and $\varepsilon=-11.0$ eV (with respect to the conduction band bottom) in the dilute limit. Our calculated formation energy is in good agreement with our explicit supercell calculations, the deviation being of 3% for the smallest supercell [Fig. 2(f) and (g)].

Since Li2O2 exhibits localized polarons in fairly small supercells, in Fig. 2(f) we compare the wavefunction obtained within our method and that obtained via an explicit DFT calculation using the self-interaction correction scheme described in Ref. Sio et al., . Apart from the slight asymmetry in the wavefunction obtained via the DFT calculation, the two results are in very good agreement. Furthermore, our result is essentially identical to previous findings based on hybrid functional calculations (Kang et al., 2012). The agreement with explicit DFT calculations validates a posteriori our initial assumption leading to Eq. (1). Taken together, the results in Fig. 1 and 2 indicate that our theory is able to describe both large and small polarons on the same footing. A more in-depth analysis of the formalism and additional tests are provided in Ref. Sio et al., .

In summary, we developed a theoretical and computational approach that allows us to investigate, for the first time, polaron energies and wavefunctions across the length scales, without resorting to supercell calculations. Our work opens up many new directions in polaron physics. For example the spectral decomposition encoded in the Fourier amplitudes $A_{n{\bf k}}$ and $B_{{\bf q}\nu}$ could be used to construct model Hamiltonians with realistic materials parameters. This additional step will make it possible to supplement DFT with path-integrals or diagrammatic Monte Carlo calculations (Feynman, 1955; Prokof’ev and Svistunov, 2008), and ultimately open the way to predictive ab initio many-body calculations of polarons in real materials.