Through the use of the piecewise-linearity condition of the total energy, we correct the self-interaction for the study of polarons by constructing nonempirical functionals at the semilocal level of theory. We consider two functionals, the $\gamma $DFT and $\mu $DFT functionals, both of which are based on the addition of a weak local potential to the semilocal Hamiltonian to enforce the piecewise-linearity condition. We show that the resulting polaron properties are in good agreement with reference hybrid functional calculations. This supports the use of semilocal functionals for calculating polaron properties.

## I. INTRODUCTION

The electron self-interaction is a spurious interaction that arises in density functional theory (DFT) due to inherent approximations in the exchange-correlation functional. For a many-electron system, two descriptions of the self-interaction can be distinguished. In an orbital-dependent formulation, the one-body self-interaction results from the self-terms of the Coulomb interaction and is found to vanish in Hartree-Fock theory, where the Hartree and exchange self-terms exactly compensate. The many-body self-interaction corresponds to the deviation of the energy functional from the piecewise-linearity condition upon electron addition.^{1–5} The piecewise linearity is a property of the exact functional, as first emphasized by Perdew *et al.*^{1} and by Yang *et al.*,^{4} with two complementary interpretations based on ensemble states and on replicas, respectively. Standard energy functionals deviate from the piecewise linearity, lacking derivative discontinuities at integer electron occupations.^{1,6} Suppressing the many-body self-interaction coincides with restoring the piecewise linearity of the total energy and thus with achieving derivative discontinuities of the energy at integer occupations.^{1,6}

The piecewise-linearity condition holds implications for systems involving integer number of electrons. For instance, in the hydrogen cation system, a piecewise-linear functional yields equilibrium bond length and dissociation energy closer to experiment compared to standard semilocal functionals.^{2} Another significant aspect is the ability of a functional to localize charges, which is often not achieved by standard semilocal functionals. This is due to the electron self-interaction, which causes charge delocalization. A prototypical example of charge delocalization is represented by polarons, which are quasiparticles consisting of a localized charge coupled with self-induced lattice distortions.^{7} The critical role of the self-interaction in determining the charge localization of polaronic states has led to extensive studies in the literature.^{8–14} In semiconductors or insulators, the formation of a polaron results in an energy level within the bandgap. In particular, an electron polaron is associated with an occupied energy level below the conduction band, while a hole polaron is associated with an unoccupied energy level above the valence band. The electron self-interaction reduces the energy separation between the polaron state and the respective delocalized band edge, thus favoring charge delocalization. When the energy cost of polaronic lattice distortions overcomes the energy gain from charge localization, the polaron state destabilizes. Consequently, the polaron charge delocalizes and the polaronic distortions vanish. These aspects are crucial when modeling small polarons, which are localized over short length scales, comparable to the extent of lattice bonds, and are associated with deep energy levels in the bandgap. Suppressing the self-interaction, thus, becomes essential in polaron calculations.

Various approaches have been proposed to address the self-interaction in DFT. In 1981, Perdew and Zunger introduced a method to suppress the self-interaction for each electronic state.^{15} Subsequently, approaches that specifically address the self-interaction of excess charges were proposed.^{16,17} More recently, Sio *et al*. introduced an *ab initio* formulation for polarons and derived therefrom a self-interaction-corrected functional.^{11,12} However, it remains unclear to which extent these approaches guarantee the piecewise linearity of the energy functional. Meta-GGA functionals partially tackle the absence of the derivative discountinuity but may struggle with localizing polarons.^{18} Hybrid functionals, which admix a fraction of Fock exchange to the semilocal exchange, can be tuned to suppress the self-interaction leading to polaron localization.^{6,8,10,19–35} However, the computational cost of hybrid functional calculations significantly exceeds that of semilocal calculations, making them impractical for large systems. Hubbard-corrected DFT+ $U$ functionals can suppress the many-body self-interaction by introducing on-site Coulomb repulsion for the occupation of specific orbitals.^{9,23,24,36–46} In the context of polarons,^{9,23,24,46} these orbitals are selected in consideration of their contribution to the polaron state. In another method, the piecewise-linearity condition is enforced on each electron state.^{47,48} Among previous approaches to address the self-interaction at the semilocal level, it is also worth mentioning a scheme relying on the electron chemical potential as slope of the total energy.^{27}

The $\gamma $DFT semilocal functional^{49,50} corrects for the self-interaction error through the addition of a weak local potential solely dependent on the polaron density rather than on specific orbitals. The polaron properties obtained with hybrid functionals, DFT+ $U$, and $\gamma $DFT have been extensively compared.^{46,49–51} Notably, it has been observed that enforcing the piecewise-linearity condition offers the advantage of yielding ground state and transport properties of polarons that remain robust irrespective of the adopted functional.^{46,49–51} Moreover, this robustness becomes pivotal for the reliability of theoretical predictions, particularly in the case of polarons, which often lack experimental data.

When addressing the self-interaction, it is important to distinguish between one-electron and many-electron systems. For a one-electron system, the exact functional is free from one-electron self-interaction.^{6,15} For a many-electron system, one-body and many-body forms of self-interaction can be defined. When using hybrid functionals or semilocal functionals, the description of the many-body self-interaction is in competition with the description of the one-body self-interaction.^{49} To illustrate this, we consider the global hybrid functional PBE0( $\alpha $),^{52} which admixes a fraction $\alpha $ of Fock exchange with a fraction $1\u2212\alpha $ of Perdew–Burke–Ernzerhof (PBE) exchange.^{53} For $\alpha =1$, the PBE0( $\alpha $) functional reduces to a Hartree–Fock-like functional, thus suppressing the one-body self-interaction. At variance, for a fraction $\alpha = \alpha k$ of Fock exchange, the total energy becomes piecewise linear, thereby suppressing the many-body self-interaction.^{49,50} Within this context, it has been shown that for a many-electron system, one should correct for the many-body self-interaction, as this allows one to include electron screening effects, which are not accounted for when correcting for the one-body self-interaction.^{49,50} Specifically, for the PBE0( $\alpha $) hybrid functional, the many-body self-interaction energy correction can be expressed as the one-body self-interaction energy correction divided by the high-frequency dielectric constant.^{49,50}

In this work, we construct nonempirical semilocal functionals to model polarons by suppressing the many-body self-interaction in the polaron state. We here introduce the $\mu $DFT functional, which is conceptually similar to the $\gamma $DFT functional developed previously.^{49,50} In both functionals, a weak local potential, which is linearly dependent on the polaron density, is added to the semilocal Hamiltonian. In the $\gamma $DFT functional, the localizing potential is constructed as a derivative of the exchange-correlation functional with respect to the polaron charge, modulated by a strength $\gamma $. In the $\mu $DFT functional, the localizing potential is taken to be proportional to the polaron density modulated by a strength $\mu $. The parameters $\gamma $ and $\mu $ are nonempirically fixed by enforcing the piecewise-linearity condition of the total energy. This is achieved by imposing the condition of constant polaron level upon charge addition, after the inclusion of proper finite-size corrections in the energetics. We show that both functionals yield polaron properties in agreement with reference hybrid functional calculations, despite their semilocal functional form. This supports the use of nonempirical semilocal functionals to study polarons in DFT calculations.

## II. PIECEWISE-LINEAR FUNCTIONALS

We illustrate our schemes for localizing polarons at the semilocal level of theory. The fundamental concept behind these schemes is that polaron localization can be achieved by introducing a weak local potential to the semilocal Hamiltonian, promoting localization over delocalization. This stems from the observation that localized and delocalized states are often in close energetic competition.^{54,55}

The notion that a localized potential could effectively address the self-interaction by enforcing the piecewise-linearity condition has previously been employed in the DFT+ $U$ scheme^{45} and in the approach proposed by Lany and Zunger.^{9} In both these schemes, the $U$ correction is applied to specific atomic orbitals. In contrast, the $\gamma $DFT and $\mu $DFT functionals are based on adding to the semilocal Hamiltonian a weak local potential that self-consistently originates from the polaron density and affects all states of the system, rather than being limited to selected atomic orbitals.

^{49,50}the weak local potential is essentially the first derivative of the exchange-correlation potential with respect to charge addition, modulated by a strength $\gamma $. This results in the following expression:

^{53}

^{,}$ V x \sigma $ is the PBE exchange-correlation potential, $q$ is the polaron charge, and $\sigma $ is the spin. The total energy corresponding to the Hamiltonian in Eq. (1) is

^{51}

The $\gamma $DFT and $\mu $DFT functionals share several properties. First, electron and hole polarons are treated analogously since the localizing potential acts on all states in the same way. Moreover, for $q=0$, the localizing potential vanishes and the electronic structure coincides with that of a regular PBE calculation. In Eqs. (2) and (4), the response of the total electron density upon the addition of the polaron charge is calculated by finite differences, which requires performing two successive calculations. The $\gamma $DFT and $\mu $DFT schemes also have several practical advantages. First, the same Hamiltonian is applied to all the electronic orbitals of the system, which ensures that the orthogonality of the wave functions can be achieved through the use of standard diagonalization algorithms. Second, the localized nature of the additional potential facilitates a smooth convergence of the Kohn–Sham equations. Third, it is rather straightforward to implement a local potential in existing DFT codes.

^{56}this corresponds to a polaron energy level $ \u03f5 \, p \xi $ being constant with respect to fractional electron occupation, namely,

^{50}the condition in Eq. (5) is achieved when the neutral and charged polaron levels coincide, namely, $ \u03f5 p \xi k(q)= \u03f5 p \xi k(0)$. We here illustrate the workflow to find $ \xi k$ in the case of an electron polaron, which corresponds to the addition of an electron ( $q=\u22121$). The case of a hole polaron is analogous and corresponds to the removal of an electron ( $q=+1$). First, one considers an initial structure in which the symmetry of the lattice is broken to accomodate the polaron lattice distortions. This is a crucial step since localized and delocalized polaronic states can be found in close energetic competition. For instance, for an electron polaron, the energy level of the extra electron coincides with the minimum of the conduction band if the symmetry of the lattice is not broken. In this case, the energy gain due to electron localization vanishes, and a structural relaxation would yield the delocalized state. The distortion can be imagined on the basis of how the lattice would respond to polaron charge localization at a given site. Once the structure is distorted, one makes an initial guess of the parameter $\xi $ to localize the polaron. For the chosen value of the parameter $\xi $, one performs structural and electronic relaxations. We note that the evaluation of the energy in Eq. (2) is not required in the search of the polaronic structure. Hence, the costs of $\gamma $DFT and $\mu $DFT structural relaxations essentially coincide with that of PBE calculations. For a sufficiently large value of $\xi $, the relaxed structure is characterized by polaronic structural distortions, which we denote $ R \u2212 1 \xi $. Given the structure $ R \u2212 1 \xi $, one evaluates the extent by which the piecewise-linearity condition is satisfied. The polaron level $ \u03f5 p \xi (\u22121)$ is found as the Kohn–Sham level of the last occupied state in a calculation with $q=\u22121$ and structure $ R \u2212 1 \xi $. The polaron level $ \u03f5 p \xi (0)$ is found as the Kohn–Sham level of the first unoccupied state in a calculation with $q=0$ and structure $ R \u2212 1 \xi $. We note that the neutral energy level $ \u03f5 p \xi (0)$ is equal to the PBE level $ \u03f5 p 0(0)$ because of the prefactor $q$ in Eqs. (1) and (3). Then, one repeats the same procedure for another value of $\xi $ that leads to polaronic lattice distortions. Assuming the linearity of polaron levels with $\xi $,

^{50}one can find an initial estimate of $ \xi k$, such that $ \u03f5 p \xi k(\u22121)= \u03f5 p \xi k(0)$. The structure is then optimized for such value of $ \xi k$, and the procedure is repeated until $ \xi k$ converges.

^{51}which is added to the Hamiltonian to shift the energy levels of the band resonating with the polaron state. Such a scissor operator can be constructed as a sum of projectors defined with the wave functions $ \psi i \sigma \xi $ obtained during the self-consistent optimization of the Kohn–Sham equations, namely,

^{57}

^{,}

^{50}

^{,}

^{50}

^{58}These effects are present for both the charged and neutral polaronic states.

^{59}To overcome these issues, finite-size electrostatic corrections need to be applied.

^{57–60}This can be achieved using the scheme introduced by Falletta, Wiktor, and Pasquarello,

^{59}which extends the finite-size correction scheme of Freysoldt, Neugebauer, and Van de Walle

^{60}to the case of defects involving built-in lattice distortions. The correction for the total energy of a calculation with polaronic charge $ q \u2032$ is

^{59}

^{,}

^{58,60}$ \epsilon \u221e$ is the high-frequency dielectric constant, and $ \epsilon 0$ is the static dielectric constant, and $ q pol=\u2212q(1\u2212 \epsilon \u221e/ \epsilon 0)$ is the ionic polarization charge associated to the polaron lattice distortions. The charge $q$ is set to $\u22121$ for electron polarons and to $+1$ for hole polarons. The charge $ q \u2032$ is set to 0 for the neutral polaronic state, to $\u22121$ for the electron polaron state, and to $+1$ for the hole polaron state. The correction for the polaron level is obtained through Janak’s theorem and is given by

^{59}

^{,}

^{61,62}Without the inclusion of such finite-size corrections, the polaron formation energies would be noticeably underestimated.

^{46,50}For simplicity, we here include these corrections in the calculated energetics, without explicitly showing them in the equations.

## III. RESULTS

The calculations are performed with the quantum espresso suite.^{63} The implementations of the $\gamma $DFT functional and of the self-consistent scissor operator are available in version 7.2, while the implementation of the $\mu $DFT functional will be made available for incorporation in the next official release of the code. The core–valence interactions are described by normconserving pseudopotentials.^{64} We consider the electron polaron in BiVO $ 4$ and the hole polaron in MgO. We model BiVO $ 4$ with a 96-atom orthorhombic supercell ( $a=10.34$ Å, $b=10.34$ Å, $c=11.79$ Å) and MgO with a 64-atom cubic supercell ( $a=8.45$ Å). The energy cutoff is set to 100 Ry in all calculations. The lattice parameters are determined at the PBE level of theory for the pristine system. The Brillouin zone is sampled at the $\Gamma $ point. Convergence tests involving larger simulations cells and denser **k**-point samplings have shown that the polaron formation energies obtained with the present computational setup are accurate within 0.06 eV.^{50} Through the application of finite electric fields^{65} at the PBE level of theory, we determine the high-frequency and static dielectric constants, which are used for the finite-size corrections. For BiVO $ 4$, our calculations give $ \epsilon \u221e=5.83$ and $ \epsilon 0=64.95$. For MgO, we obtain $ \epsilon \u221e=2.77$ and $ \epsilon 0=10.73$. In the following, we perform polaron calculations with the $\gamma $DFT and $\mu $DFT semilocal functionals and compare them with reference PBE0( $\alpha $) hybrid functional calculations.^{49,50}

First, we determine the parameters of the functionals such that the piecewise-linearity condition in Eq. (5) is satisfied. With this aim, we initialize the structure with polaronic distortions that break the symmetry of the lattice. For the electron polaron in BiVO $ 4$, we elongate the V–O bonds around a V atom. For the hole polaron in MgO, we elongate the O–Mg bonds around an O atom. Next, we perform electronic and structural relaxations to find the parameter $ \xi k$ and the corresponding structure $ R q \xi k$, as described in Sec. II. The parameters found with this procedure are given in Table I. As an example, we show in Fig. 1 the charged and neutral polaron levels obtained with the $\mu $DFT functional at the geometry $ R q \xi k$. We remark that the neutral level $ \u03f5 p \mu (0)$ is essentially independent of $\mu $ since the Hamiltonian $ H \sigma \mu (0)$ depends on $\mu $ only through the structure $ R q \mu $, which varies only marginally with $\mu $.

System . | γ_{k}
. | μ_{k}
. | α_{k}
. |
---|---|---|---|

BiVO_{4} | 1.80 | 2.85 | 0.14 |

MgO | 1.96 | 10.90 | 0.34 |

System . | γ_{k}
. | μ_{k}
. | α_{k}
. |
---|---|---|---|

BiVO_{4} | 1.80 | 2.85 | 0.14 |

MgO | 1.96 | 10.90 | 0.34 |

System . | γDFT
. | μDFT
. | PBE0 (α_{k})
. |
---|---|---|---|

BiVO_{4} | 1.82 | 1.82 | 1.80 |

MgO | 2.23 | 2.23 | 2.20 |

System . | γDFT
. | μDFT
. | PBE0 (α_{k})
. |
---|---|---|---|

BiVO_{4} | 1.82 | 1.82 | 1.80 |

MgO | 2.23 | 2.23 | 2.20 |

Finally, we compare the formation energies obtained with the piecewise-linear functionals using the expression in Eq. (9). We find that the formation energies obtained with $\gamma $DFT and $\mu $DFT are very close, with deviations of at most 0.03 eV (cf. Table III). A good agreement is also found with the PBE0( $ \alpha k$) results. This further corroborates the finding that upon enforcing the piecewise-linearity condition the polaron properties turn out to weakly depend on the adopted functional.^{46,49–51} This is particularly relevant in the case of the hole polaron in MgO, for which the experimental characterization remains inconclusive^{66,67} and the theoretical description contentious.^{8,10,50,59,68}

## IV. CONCLUSIONS

We investigate in this work the performance of suitably constructed semilocal density functionals in addressing polaron properties. As case studies, we consider both electron and hole polarons, for which our approach works in an equivalent manner. The two considered semilocal functionals, denoted $\gamma $DFT and $\mu $DFT, depend on a single parameter, which we fix in a nonempirical fashion by enforcing the piecewise-linearity condition and thus suppressing the many-body self-interaction. Both functionals are successful in achieving polaron localization. We find that the achieved results are in excellent agreement with those of reference hybrid-functional calculations. The present analysis encompasses both electronic and structural polaron properties, such as polaron densities, polaronic bond distortions, and polaron formation energies. This corroborates the finding that these polaron properties are robust and that they only weakly depend on the functional adopted, provided the piecewise-linearity condition is satisfied. Overall, our study shows that polaron properties can be accurately described through the use of semilocal functionals, without resorting to computationally more expensive hybrid-functional calculations.

## ACKNOWLEDGMENTS

We thank Michele Pavanello for useful discussions. The calculations have been performed at the Swiss National Supercomputing Centre (CSCS) (Grant under Projects ID No. s1122) and at SCITAS-EPFL.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors declare no conflict of interest.

### Author Contributions

**Stefano Falletta:** Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Alfredo Pasquarello:** Conceptualization (equal); Data curation (supporting); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (lead); Supervision (lead); Validation (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in Materials Cloud, Ref. 69.

## REFERENCES

*Ab Initio*finite-size corrections for charged-defect supercell calculations

*Ab initio*molecular dynamics in a finite homogeneous electric field