First-Principles Supercell Calculations of Small Polarons with Proper Account for Long-Range Polarization Effects

We present a density functional theory (DFT) based supercell approach for modeling small polarons with proper account for the long-range elastic response of the material. Our analysis of the supercell dependence of the polaron properties (e.g., atomic structure, binding energy, and the polaron level) reveals long-range electrostatic effects and the electron-phonon interaction as the two main contributors. We develop a correction scheme for DFT polaron calculations that significantly reduces the dependence of polaron properties on the DFT exchange-correlation functional and the size of the supercell in the limit of strong electron-phonon coupling. Using our correction approach, we present accurate all-electron full-potential DFT results for small polarons in rocksalt MgO and rutile TiO$_2$.


Introduction
The electron-phonon (el-ph) interaction is fundamental to materials. It mediates, for example, the excitation of phonons in response to electronic excitations, which is especially pronounced in polar materials. These phonon excitations can stabilize a lattice distortion around a single excess charge (electron or hole). The excess charge and its accompanying lattice distortion are then referred to as a quasiparticle or more specifically as polaron. The formation and migration of polarons determine the properties of functional materials, such as their catalytic [1,2] and photovoltaic [3] behavior. The direct observation of polarons in experiments, e.g. with electron-paramagnetic resonance [4], UV/IR spectroscopy [5], or scanning tunneling microscopy or spectroscopy [6] is difficult, and computational studies are required to interpret the experimental data correctly. In this work, we develop a new method that addresses challenges faced in computational modeling of small polarons in materials with strong electronphonon coupling, in particular in oxides, with density-functional theory (DFT).
Polarons can be classified by their size as quantified by the extent of their total wave function (electrons and ions). Large polarons are delocalized over several unit cells and usually appear, if the el-ph interaction is weak. Such polarons were first investigated by Fröhlich [7]. In contrast, small polarons are mainly localized on one atomic site and form when the el-ph interaction is strong. Intermediate polarons [5] cover the size range in between. Pioneering work on small polarons was performed by Holstein [8]. He took only shortrange interactions into account and identified the Fröhlich coupling constant α Fröhlich [7] as good indicator for the el-ph interaction strength. Oxides fall into the intermediate to strong coupling regime, i.e., α Fröhlich > 1. For instance, for MgO α Fröhlich is 4.4 and for rutile TiO 2 2.2. We therefore expect small polaron formation in both of these oxides. However, since MgO and TiO 2 are strongly ionic, the distortion of the ionic lattice can be long-ranged in violation of the polaron classification scheme. Such polarons in which the excess charge is localized, but the lattice distortion delocalized are referred to as small Fröhlich polarons.
To describe small Fröhlich polarons accurately in computational materials modeling, both long-and short-range interactions have to be treated appropriately. How to accomplish this task in DFT calculations that employ supercells, whose extend is typically smaller than the ionic lattice distortions, is the subject of this paper. Since small polarons can be regarded as a special type of a point defect, our study is also useful for point defect calculations of this type, which have so far eluded a reliable theoretical treatment.
The paper is structure as follows. In Section 2 we derive the electrostatic and the electron-phonon contributions to the elastic long-range response of a material to a localized excess charge. We will then use this derivation to develop a correction scheme that removes artificial interactions from the supercell approach to obtain polaron properties in the dilute limit. In light of our new understanding, we analyze shortcomings of previous polaron approaches in Section 3. In Section 4, we demonstrate the efficiency of our approach for hole polarons in MgO and electron polarons in rutile TiO 2 .
2 Elastic long-range behavior DFT in combination with the supercell approach has become the method of choice for the ab inito calculation of point defects in solids. However, the supercell approach suffers from finite-size effects, especially for charged defects. These finite-size effects include the interaction of the excess charge with its periodic images, with the compensating constant background charge introduced to keep the unit cell neutral, and with the periodic constraint on the atomic relaxation. To overcome these finite-size limitations, two strategies are commonly used: (a) extending the supercell and extrapolating to the dilute limit based on a scaling law, or (b) applying an a posteriori correction. For (a) only general knowledge about the size dependence is necessary. For example, the formation energy of a charged defect in the bulk as a function of the supercell size L (L = Ω 1/3 , where Ω is the supercell volume) can be written as an inverse powerlaw: where E(∞) is the formation energy in the dilute limit. This scaling law was derived by Makov and Payne [9]. The disadvantage of this procedure is that at least three supercell calculations of increasing size are needed to fit E(∞) in Eq. (1), which is computationally very demanding, especially if atomic relaxations are included 1 . Conversely, approach (b) requires an appropriate physical model for the long-range interactions in the solid. If only the electronic response to the excess charge is considered, its long-range contribution to the energy is described by a term proportional to 1/ǫ ∞ r (e.g. Ref. [11]). However, if the ionic response cannot be neglected, the problem becomes challenging, and so far this case has not been solved. It has been suggested that the long-range elastic contribution is similar to the electronic one, but with the high-frequency dielectric constant ǫ ∞ replaced by the static one ǫ 0 , i.e., the long-range potential behaves classically like 1/ǫ 0 r. However, corrections based on this assumption generally overestimate E(∞), especially for vacancies. This overestimation has two reasons. First, the aforementioned long-range behavior is a crude approximation, neglecting all details of the underlying phonon structure. Second, the short-range screening of the excess charge depends on the el-ph coupling, which determines the interplay of the electronic and ionic responses. Screening can then be more efficient than in the static limit. In the following we analyze the screening effects in detail and show that only in the strong-coupling limit of the el-ph interaction the substitution of ǫ ∞ with ǫ 0 is a good approximation.
We start by splitting the long-range elastic potential V lr into the electronphonon interaction V lr el-ph and electrostatic potential V lr el-st : V lr el-st is generated by the charge density ρ d (r) of the localized excess charge. The Fourier transform of V lr el-st is then given by: where ρ d (k) is the Fourier transform of ρ d (r).
To obtain a corresponding expression for V el-ph we first have to introduce additional assumptions. First, we will focus on ionic crystals. Second, we only consider the interaction of an electron with a single phonon at a time, neglecting higher-order contributions. Third, we assume that the adiabatic approximation (factorization of the electron and phonon wave functions) and strong electronphonon coupling limit (α Fröhlich 4) are applicable. With these assumptions the long-range part of V el-ph reduces to [12]: The potential in Eq. (4) is attractive, lending further stabilization to the polaron. An analytic expression for g lr el-ph was recently derived by Verdi and Guistino [13]: where ν labels the phonon mode, ω kν is the corresponding phonon frequency of ion κ with mass M κ . Z * κ is the Born effective charge tensor and e κν (k) are the phonon eigenvectors of the dynamical matrix.
Equations (4) and (5) describe the scattering of all phonon modes with ρ d . Thus, the long-range behavior of the el-ph interaction depends on the phonon structure across the entire phonon Brillouin zone, and the elastic behavior is not captured by the classical 1/ǫ 0 r limit. If we only consider the interaction of ρ d with a single dispersion-less longitudinal-optical mode ω LO , we recover the limit of the Fröhlich electron-phonon interaction in the strong-coupling limit, which was first investigated by Pekar [14]. With the Fröhlich matrix element (for the anisotropic case we refer to [15]): we obtain the potential: Upon substituting Eq. (7) and Eq. (3) into Eq. (2), we finally arrive at the classical limit of a screened potential for a localized charged distribution in an anisotropic medium: The el-ph potential given by Eq. (7) is an upper bound and, consequently, Eq. (8) is also an upper bound. This explains why any correction based on Eq. (8) overestimates the actual limit. We find that, despite the approximations we made, V lr (k) in Eq. (8)  Based on the knowledge of the long-range behavior, the errors due to finte size of the supercell can be corrected using a posteriori methods, such as the method of Freysoldt et al. [11]. For technical details we refer to [11,16]. Generalizing the Freysoldt method to an arbitrary interaction potential V (r) and anisotropic media (in the standard approach of Freysoldt et al., V (r) = 1/ǫ ∞ r), the correction for the interaction energy is obtained as the difference between the energy of the artificial lattice of charged defects, E latt , and the energy of an isolated defect, E iso : where V can be V lr el-ph , V lr el-st or the sum of both, and q d (k) is the Fourier transform of the excess charge distribution, and q is the total charge. Taking into account the alignment terms q∆V [11,16] A summary of the Freysoldt et al. correction scheme including the meaning of the alignment terms can be found in the Appendix C, the corrected energy is obtained as: Having derived the correction for the elastic contribution, we can apply it to the polaron problem and investigate the effects of the two parts in Eq. (2) separately.
3 The Polaron in a supercell

The charged supercell
An important property of a polaron is its binding energy where the energies have not been corrected for finite-size effects, yet. The plus in E ± bind corresponds to electron removal (hole polaron), while the minus sign corresponds to electron addition (electron polaron). E polaron is the total energy of the distorted system (polaron geometry), E perf the total energy of the undistorted system. The number of electrons in the system are given in parenthesis, with N corresponding always to the neutral system. A negative E ± bind indicates an energy gain and a stable (self-trapped) polaron.
In the following we focus on the hole polaron for brevity, since only small adjustments of the formalism are needed for the electron polaron case. The simplest way to calculate the polaron binding energy is straightforward: in Eq. (11) E polaron (N ∓ 1) is computed with DFT and full structure relaxation in the charged supercell. To ease the system out of possible high symmetry configurations an initial symmetry-breaking distortion might have to be applied. Finite-size effects are expected to be small, since the elastic long-range interaction falls off with 1/ǫ 0 r and the static dielectric constant ǫ 0 is usually large ( 10) for ionic crystals (however, as demonstrated and explained below, the dependence of the polaron binding energy defined by Eq. (11) on the approximations in the exchange-correlation functional is strong). The supercell dependence of E + bind for MgO is shown in Fig. 2, panel a, where we used HSE06 hybrid functional [17,18] with the fraction of exact exchange α = 1 [denoted HSE06(α = 1); see Section 4 for more computational details]. We find a small hole polaron mainly localized at the central oxygen atom. The displacements of the nearest neighbors are of the order of 0.1Å and decaying fast away from the center. The shape of the excess charge density distribution is p-like. For sufficiently large supercells, when the long-range regime is valid, the dependence of the binding energy on the supercell size L becomes 1/ǫ 0 L. From the slope of E + bind (1/L) at 1/L = 0 we obtain ǫ 0 = 10.32, in good agreement with the experimental static dielectric constant for MgO ǫ 0 = 9.8.
Next, we calculate the correction for the artificial electrostatic interaction due to the periodic arrangement of the holes and their interaction with the constant background, using Eq. (9) with the potential1/ǫ ∞ r. To model the excess charge density ρ d (r) needed here and for following finite-size corrections, we fit the envelope of the KS eigenstate density (decays exponentially for a localized state) with an exponential function ρ model = A exp(−|r − r 0 |/γ), where A is a normalization constant, r 0 is the center of the polaron, and γ the fitting parameter corresponding to the polaron radius (cf. App. A). Additionally, we calculated the alignment term ∆V in Eq. (10) between the charged, neutral, and model (i.e., including the model excess charge density compensated by a constant background charge) systems following the approach outlined in Ref. [16]. After this correction, according to Eq. (2) the remaining contribution is due to the long-range electron-phonon interaction. This contribution is shown by the blue line in Fig. 2, panel a. The line is almost perfectly straight, and the slope is equal to ǫ −1 ∞ − ǫ −1 0 , where ǫ 0 = 10.32 is taken from the fit of E + bind presented above, and ǫ ∞ = 2.4 is obtained from an independent calculation 2 . This analysis explains the role of different long-range interactions in Eq. (2) in the supercell dependence of polaron properties. Thus, the approximations in Eq. (7) work well for MgO, which is expected since it has only one longitudinal optical phonon mode, strong el-ph coupling, and is an isotropic material. However, we find that the polaron binding energy defined by Eq. (11) is extremely sensitive to the approximations in the exchangecorrelation functional. Figure 2, panel b, shows the dependence of the binding energy on the fraction of exact exchange α in the HSE06 functional. Within a small range ±0.05 of α around the standard value (0.25) the binding energy changes by about 0.5 eV. This leads to a qualitative change in small polaron stability, from a stable self-trapped polaron (negative binding energy) to an unstable small polaron (positive binding energy). This strong functional dependence makes even a qualitative assessment of the existence of a self-trapped polarons impossible. Several approaches have been suggested in the literature for determining the correct or at least optimal value of α [19,20,21,22,23,24,25]. Here we focus on restoring the IP theorem [21] as a consistent DFT-based solution of the problem.
In (exact) DFT within the scope of Kohn-Sham (KS) scheme the vertical ionization potential IP should be equal to the negative of the highest occupied KS state energy ε ho in the system: where E(N − 1) and E(N ) are total energies of the ionized and neutral system, respectively. In this work we refer to this relation as IP-theorem, but it is also know as HOMO-I condition [19] or Generalized Koopmans' Theorem [21], and is directly related to the straight-line dependence of the total energy on occupation of the highest-occupied state [26] or the fact that the position of ε ho is independent on its occupation. Equation (12) is always correct for any extended (delocalized) state, as was already pointed out by Janak (1978) and extended to the case of the generalized KS scheme by Perdew et al. [27]. However, for a given density-functional approximation (DFA) Eq. (12) does not necessarily hold, if the orbital is localized, unless the satisfaction of the straight-line condition is explicitly included in the design of the functional. The deviation from the straight line ∆ XC (α) is described by two contributions to Eq. 12: with the self-interaction error Π causing a convex curvature of the total energy as a function of occupation, and the orbital relaxation Σ a concave curvature. The optimal α = α opt minimizing the self-interaction error [20,28] is then determined from the condition ∆(α opt ) = 0. The straight-line theorem (Eq. (12)) was originally proven for finite systems, and transferring it to a solid with periodic boundary conditions needs special care. For any finite supercell with volume Ω the energy of the artificial electrostatic interactions due to the periodic arrangement [E el-st corr (Ω), obtained using Eq. (9) with potential from Eq. (3)] has to be removed from E(N − 1): since it would only vanish in the limit of an infinite supercell. Combining Eq. (14) and Eq. (11), we get: where E el-st corr stands for the artificial electrostatic interaction energy for the distorted geometry, since for the perfect geometry it is zero. The quantity is calculated using only neutral unit cells, with the energy of distortion from perfect to polaronic geometry ∆E polaron = E polaron (N ) − E perf (N ) and the polaron level energy with respect to the VBM E 0 = ε polaron ho (N )−ε polaron VBM (N ). According to Eq. (15), when ∆ XC (α) is zero, E 0 bind represents the polaron binding energy corrected for the artificial electrostatic interaction.
E 0 bind is shown in Fig. 2, panel a, as red line (top-most line), where the optimized polaron geometry of the charged supercell is used. Note that, despite including only quantities calculated using neutral unit cells, E 0 bind has a strong dependence on the unit cell size. As discussed below, this dependence is due to the interaction of the ionic relaxations in different unit cells. Taking the difference between the blue and the red lines, we find that the exchange-correlation error ∆ XC is practically independent on the unit cell size (green line in Fig. 2, panel a), starting from the smallest supercell with 64 atoms we have considered. This implies that ∆ XC (α) could be calculated even in the smallest supercell and then removed in any larger supercell. In order to obtain the optimized α = α opt we have to remove E el-st (Ω) from the binding energy E + bind and determine the intersection with E 0 bind . The result is shown in Fig. 2, panel b, and we obtain α opt = 0.48. Since the dependence on α is not linear, at least three different values of α have to be calculated to estimate α opt . Additionally for each value of α the dielectric constant ǫ ∞ has to be calculated. Thus, the simulation of the polaron in a charged supercell is computationally demanding, since it is extremely sensitive to the underlying functional. In the next subsection we demonstrate an approach to overcome this problem.

The neutral supercell
As mentioned above, E 0 bind in Eq. (15) is equal to the polaron binding energy corrected for the artificial electrostatic interaction, only when ∆ XC (α) vanishes. However, similar to previous work [29,30] we find that E 0 bind is far less sensitive to the underlying functional than E + bind , as can be seen for MgO in Fig. 2, panel b. The same is true for TiO 2 , but the remainig dependence is larger than for MgO (see Fig. 3). This has an interesting implication: E 0 bind is the polaron binding energy with most of the exchange-correlation error removed (because for an exact functional E 0 bind is equal to the polaron binding energy). The reason for the insesitivity of E 0 bind on the functional remains unclear [29]. As a consequence, even with PBE we find a stable self-trapped hole polaron in MgO, which is not the case when charged supercells are used. Also, we find that the polaron level with respect to the band edge (E 0 ), calculated using a neutral supercell, is insensitive to the functional, as can be seen in Fig. 4. A stronger functional dependence of E 0 is expected when the character of the polaronic state or states of the band edges are sensitive to the functional.
Using E 0 bind for calculating polaron binidng energies has been first implicitly introduced by Zawadski et al. [30]. The independence of E 0 bind on the functional has been discussed by Sadigh et al. [29]. In their work, Sadigh et al. have also suggested a way to obtain forces for a polaronic distortion directly using E 0 bind potential-energy surface (PES). This facilitates the calculation of accurate elastic response to the excess charge at the level of a hybrid functional, but at the cost of a PBE calculation.
Thus, using the E 0 bind PES instead of charged supercells allows us to significantly reduce the functional dependence. Naively, one may expect that the supercell dependence is also reduced, since only neutral supercell calculations are performed. However, this is not the case. As can be seen from Fig. 2, panel a, the dependence of E 0 bind on the supercell size is much stronger than in the case of charged suercells. This dependence is due to the artificial interaction between ionic relaxation fields in different supercells. Indeed, the E + bind − E el-st corr and E 0 bind supercell dependence are practically identical and correspond to the long-range part of the electron-phonon interaction potential given by Eq. (7) in the strong el-ph coupling limit. This understanding allows us to introduce an a posteriori correction E corr removing the dependence of E 0 bind on the supercell size. To remove the artificial interaction terms, we use the approach of Freysoldt et al. [11,16], but for a different long-range potential, namely the one given by Eq. (7). This new correction scheme relies on the assumption of a strong el-ph coupling, but, as demonstrated below, works reasonably well also for intermediate coupling regimes.
The polaron level E 0 also depends on the supercell size. Because of special properties of the small polaron in the adiabatic strong-coupling limit, it is possible to relate the polaron binding energy to the polaron level, in accordance with Pekar's 1:2:3:4 theorem [31] . It follows from the theorem that (see details in App. B): Thus, the correction to E 0 (Ω) in a finite supercell is expected to be about twice as large as for the polaron binding energy calculated using neutral supercells. Indeed, this is what we observe for MgO (see Fig. 3, panel a), where the absolute value of the E 0 (L) slope is almost exactly twice of the absolute value of the E 0 bind slope. For TiO 2 , the relation between the E 0 (L) and E 0 bind dependencies deviates from the one derived from Pekar's model (see Fig. 3, panel b) due to a weaker electron-phonon coupling, as discussed in detail in Section 4.
In summary, we find that in this approach the dependence on the exchangecorrelation approaximation is drastically reduced, but the finite-size effects are significantly more pronounced. However, these effects, caused by the electronphonon long-range potential (eq. (7)), can be corrected using the approach of Freysoldt et al., but with the potential V lr el-ph . This makes possible using moderately sized supercells and semi-local functionals to predict polaron properties, as demonstrated in the next section.

Polarons in rocksalt MgO and rutile TiO 2
Building on the findings and understanding obtained in the previous sections, we formulate our approach for a reliable calculation of polaron properties: 1. We obtain the atomic structure of the polaron using the PBE functional [corresponding to HSE06(α = 0)] and the approach of Sadigh et al. [29] for each supercell size.
2. HSE06(α = 1) calculations (as a limiting case) are performed for the fixed geometries obtained with PBE. This allows the estimation of the functional dependence for the systems.
3. The polaron binding energies are calculated using Eq. (15). The finitesize correction for the binding energy is calculated using Eq. (9) with the potential from Eq. (7). The correction for the polaron level is calculated as twice the correction for the binding energy.
The different sign of the correction for the hole polaron versus the electron polaron (compare panels a and b in Fig. 3) is explained by the fact that the equation for the electron affinity has to be used for the electron instead of the ionization potential for the hole. We use the hybrid-functional implementation [32] in the all-electron fullpotential electronic-structure package FHI-aims [33,34,35]. The evaluation of forces and total energies are computed with FHI-aims using the default light settings, to obtain consistent results for all unit cell sizes. As is shown in the Supplementary Information (SI), using default tight settings, which are the recommended settings for well-converged calculations, does not affect the results for the smallest supercell. As a demonstration, we apply our new approach to polarons in MgO and rutile TiO 2 . For the cubic 8-atom MgO unit cell we use a lattice constant of a = 4.211Å obtained with HSE06 (α = 0.25), and a Γcentered 8 × 8 × 8 k-grid. The number of k-points for each direction is scaled down linearly for larger supercell sizes. For the tetragonal 6-atom TiO 2 unit cell we use a = 4.64Å and c = 2.97Å obtained with the PBE functional and a 9 × 9 × 15 k-grid. Due to one more degree of freedom the positions of the atoms are optimized, too, using the PBE functional (for details see SI).
The results for a hole polaron in MgO and an electron polaron in rutile TiO 2 are shown in Fig. 3. For every supercell size we allow all atoms to relax to obtain the full elastic contribution within the cell. The corrected E 0 bind values for each supercell are shown for PBE. Clearly, the supercell-size dependence of E 0 bind for both MgO and TiO 2 agrees very well with the behavior corresponding to the electron-phonon long-range contribution described by Eq. (7). As mentioned above, the Frölich coupling constant α Fröhlich is equal to 4.4 for MgO and 2.2 for TiO 2 . Thus, MgO is better described by Pekar's model Eq. (7), and the size-corrected binding energy practically coincides with the binding energy obtained from a linear extrapolation to the dilute limit. For TiO 2 , the corrected energy deviates (surprisingly only slightly) from the extrapolated one (within 0.05 eV), reflecting approximations in Pekar's model. Also, the functional dependence of the energies is stronger for TiO 2 , indicating a larger contribution of the short-range effects to the binding energy. Additionally, we observe that the atomic structure is sensitive to the functional as well demonstrating limitations of obtaining polaron atomic geometries with only the PBE functional, even on the PES corresponding to E 0 bind , which is much less sensitive to the approximations in the functional than the PES of a charged supercell. This sensitivity is connected to delocalization errors and missing static correlation originated in the d-orbitals. However, the changes in the geometry as a function of α are still small, and we use the configurations of the perfect system obtained with PBE. We find final polaron binding energies in the dilute limit -0.38. . . These results remain both qualitatively and quantitatively consistent across a broad range of functionals generated by varying the fraction of exact exchange. This consistency is remarkable when compared to previous theoretical studies, especially for TiO 2 , since it was either shown that the small polaron formation is expected only for a certain range of a parameter, e.g. for DFT+U [6,36] or HSE(α) [37], or it was demonstrated only for a specific value of a parameter, e.g. for HSE(α=0.25) [38].
To make a connection to experimentally accessible quantities, in particular photoluminescence (PL) measurements, accurately predicting the position of the polaron level is important. Since the quantities obtained with the neutral PES E 0 bind are weakly dependent on the underlying functional, the fraction of exact exchange α can be used to tune the gap E gap to recover the experimental band gap. The main PL peak due to the small polaron formation can be expected at: For MgO the experimental band gap was measured to 7.8-7.9 eV [39], which can be simulated by a fraction α = 0.4. Based on our HSE06(α = 0.4) calculations, the PL peak should be at 6.3 ± 0.1eV. Unfortunately we could not find any experimental reference for this region of PL. For TiO 2 a fraction α = 0.2 is needed in order to reproduce the experimental band gap of 3.0 eV [40], and the corresponding photoluminescence peak is predicted to be at 2.1 ± 0.1 eV. This is in good agreement with experimental findings of P L = 2.34 eV for rutile powders [41] or direct measurements of the polaron level E 0 = 0.7 ± 0.1 eV with scanning tunneling spectroscopy [6]. We note that the results provided here only represent an upper limit for the polaron level or lower limit for the PL peak, since neither finite-temperature nor non-adiabatic effects are taken into account.

Conclusions
In this work, we developed a new approach for first-principles modelling of small polarons in materials using DFT supercell calculations. Because on the one hand, the standard charged supercell approach allows us to obtain polaron properties in the dilute limit (for moderately large finite supercells and values of ǫ 0 finite-size errors can be even neglected), but the results strongly depend on the underlying exchange-correlation functional. On the other hand, the approach of Sadigh et al. [29] significantly reduces the dependence on the functional, but, as we demonstrate, introduces a strong dependence on the supercell size. We show that the large finite-size errors in the latter approach are due to constraints imposed on the elastic response to the excess charge by the periodic boundary conditions, and suggest a way to correct the errors for finite supercells. The correction relies on the validity of Pekar's model [14] for the long-range response, based on approximations corresponding to the adiabatic strong (in Frölich's sense) electron-phonon coupling limit. As a result, our approach allows us to obtain polaron properties in the dilute limit and at the same time reduce the exchange-correlation errors, so that even semi-local functionals can be used to reliably estimate polaron level, binding energy, and atomic structure. For more accurate modelling of polaron effects on photoluminescence in materials, the use of hybrid-density functionals or methods beyond DFT, such as the GW approach [42,43,44], is still necessary. We apply the developed approach to small polarons in MgO and rutile TiO 2 . We find that the hole polaron in MgO indeed behaves as Pekar's polaron at the long range, as expected based on the large value of Frölich's constant. For electron polarons in TiO 2 , our approach also works surprisingly well, considering the weaker electron-phonon coupling in this material. Our all-electron fullpotential results support the existence of a small electron polaron in rutile TiO 2 in agreement with previous work [37,6].

A Pekar's Polaron and its relation to KS eigenstates
The objective of the appendix is to show how analytical polaron models are connected to the actual many-body problem treated with DFT. Especially, the relation of the polaron wave function to the highest occupied (ho) or lowest unoccupied (lu) KS state is discussed below.
Pekar's polaron model [14] evolves from the Fröhlich Hamiltonian [7] in the strong coupling limit, as was shown for example by Devreese [12] (cf. also citation in it for original works). In this limit, assuming adiabatic separation of ionic and electronic degrees of freedom, the electron-phonon interaction has the form: which is the classical response of a polar dielectric to an extended charge distribution. The inverse dielectric constant κ −1 = ǫ −1 ∞ − ǫ −1 0 describes the polarization of the rigid ions in the medium by the electron or hole. For simplicity, here we assume an isotropic medium (the dielectric response is described by a single constant). Let us regard Eq. 19 as a perturbation of the perfect system H perf -i.e. the single-electron Hamiltonian, where the electron has been placed at the bottom of the conduction band minimum (CBm) φ CBm with energy ε CBm of the non-interacting system (this is the scenario for the electron polaron): Following the Kohn-Luttinger perturbation theory [45] the solution of: in first order is given by: where E 0 and Φ are obtained from the solution of the effective Hamiltonian of the charge with an effective mass m * [45] (without taking into account polaronic effects): with the effective mass m * from the CBm. With Eq. (23) we recover the original problem of Pekar's polaron and E 0 is the energy of the bound (polaron) state relative to the conduction-band edge for the case of an electron polaron. Eq. (23) does not contain microscopic details. However, it can be regarded as describing asymptotic el-ph interaction far away from the localized part of the excess electron charge distribution and, thus, Φ(r) is the asymptotic solution of Eq. 23. According to Eq. (22), Φ(r) represents the envelop of the original electronic state φ CBm and is expected to decay exponentially with distance. The electron KS eigenstate φ DFT lu corresponding to ε lu (N ) in the DFT calculation at the distorted (polaron) geometry is the polaron wave function Ψ. Thus, the envelop of φ DFT lu shows the localization of ρ d (r) needed for the correction scheme Eq. (9) in order to fit ρ m (r). An example of ρ d (r) calculated with DFT is shown in Fig. 4(b).

B Pekar's 1:2:3:4 theorem
For arbitrary coupling constants in the Fröhlich Hamiltonian: with the Hamiltonian of the phonons H ph , it has been shown [31] that there exist fixed ratios of the effective kinetic energy E kin,eff , lattice distortion (phonon field) energy ∆E polaron , the polaron state energy E 0 , and the electron-phonon interaction energy E el-ph : where η depends on the value of Frölich coupling constant α F . In the limit of strong electron-phonon coupling (α F → ∞), the polaron energy is dictated by the polarization of the lattice, and η approaches 2. From this it follows: Eqs. (26) and (27) clearly show the dependence of the binding energy and the polaron level on the energy of the el-ph interaction. The latter energy is the one that remains to be corrected for the artificial supercell interactions, and thus the correction for the polaron level has to be twice of the correction for the binding energy, which leads to Eq. (17). However, these ratios Eq. (25) are only based on an effective single-particle model (Eq. 23). In our microscopic (DFT) model, additional (short-range) contributions to the energy components and deviations of η ≤ 2 lead to violation of the above ratios. In particular for TiO 2 the ratios are not preserved. However, for MgO, where the Frölich constant is 4.4, indicating indeed a strong electronphonon coupling, the ratios are close to the ones found by Pekar, and the polaron level and binding energies calculated from the model are close to the ones from DFT calculations, as described in the text.

C Freysoldt et al. correction scheme for finitesize effects in a nutshell
A repeating point in this paper is the correction of finite-size effects for supercell calculations. For completeness we present the main ideas of the correction scheme proposed by Freysoldt et al. [11,16]. Starting point is the simulation of a charged point defect in an otherwise pristine crystal causing a localized excess-charge distribution ρ d . It is assumed that for a sufficient large supercell the quantum nature of the defect is simulated properly and only long-ranged interactions do affect the defect potential in neighboring cells. If the longrange potential possess a Fourier-transformation, e.g. as shown here for V lr el-st (eq. (3)) and V lr el-ph (eq. (7)), then, it is possible to correct the biased energies a posteriori. For this, to have simple evaluable sums and integrals Freysoldt et al suggest to model ρ d with a simple isotropic function ρ m , such as an exponential or a Gaussian (the fitting of ρ d by ρ m is demonstrated in App.A). The actual detailed excess-charge distribution is not necessary to know and would change the correction only negligibly. (As Freysoldt et al. in their original paper note it is not even important to imitate the proper localization of ρ d as long as the distribution is well-localized within the supercell.) With this, it is possible to evaluate the lattice sum of the long-range potential (i.e. the potential energy due to their periodic arrangement): (for the detailed nomenclature see main text), where V (G) is the Fouriertransform of the long-range potential, and the sum runs over all reciprocal lattice vectors |G| < G cut . The cut-off G cut has to be chosen carefully to ensure convergence of the sum. Eq. 28 is the artificial energy, which has to be removed from the regarded energy (e.g. the polaron binding energy or level). What is missing is the long-range energy of the isolated defect. This is easily calculated by: (29) and the total correction is given by E corr = E latt − E iso . To obtain the desired energy in its dilute E ∞ limit the correction E corr (Ω) has to be removed from the energy E(Ω) calculated in the supercell of size Ω: The last term q∆V is the so-called alignment term and has to be considered for the following reasons: First, usually E(Ω) is calculated with respect to a reference system, often the pristine bulk system. Due to defect or the charge there might be difference in the potentials for the defect system and the pristine system even far away from the defect center. This difference can be obtained by aligning the electrostatic potentials (or Hartree potentials). Second, the absolute position of the long-range potential calculated from ρ m might not be equal to the one from the original ρ d . This difference must be aligned, too. Hence, in general the term q∆V should include these two contributions.