First-principles calculation of intrinsic defect chemistry and self-doping in PbTe

Semiconductor dopability is inherently limited by intrinsic defect chemistry. In many thermoelectric materials, narrow band gaps due to strong spin–orbit interactions make accurate atomic level predictions of intrinsic defect chemistry and self-doping computationally challenging. Here we use different levels of theory to model point defects in PbTe, and compare and contrast the results against each other and a large body of experimental data. We find that to accurately reproduce the intrinsic defect chemistry and known self-doping behavior of PbTe, it is essential to (a) go beyond the semi-local GGA approximation to density functional theory, (b) include spin–orbit coupling, and (c) utilize many-body GW theory to describe the positions of individual band edges. The hybrid HSE functional with spin–orbit coupling included, in combination with the band edge shifts from G0W0 is the only approach that accurately captures both the intrinsic conductivity type of PbTe as function of synthesis conditions as well as the measured charge carrier concentrations, without the need for experimental inputs. Our results reaffirm the critical role of the position of individual band edges in defect calculations, and demonstrate that dopability can be accurately predicted in such challenging narrow band gap materials. Accurate modeling of defects in PbTe lies on the combination of hybrid functionals with spin-orbit coupling and Green function theory. Dopability of thermoelectric materials (that can convert thermal into electrical energy and vice versa) is crucial for their improvement; however most modeling approaches fail on an atomic level, especially if spin-orbit coupling effects (that modify the band position) are present. In this work, Vladan Stevanovic and coauthors manage to accurately reproduce the electronic structure of PbTe with first-principles based density functional theory. They prove that the only approach that can model the intrinsic defect chemistry, in agreement with experimental results, is the combination of hybrid functionals (based on a screened Coulomb potential) with spin-orbit coupling and the band edge shifts calculated through the G0W0 approximation (which calculates the Green function theory that models excited-state properties of extended systems). Such level of understanding is necessary to predict the dopability of PbTe.


INTRODUCTION
The dopability of semiconductor materials plays a decisive role in device performance. Dopability refers to the carrier concentration limits achievable in a semiconductor material. These limits are set by the compensating intrinsic (or native) defects and the solubility of extrinsic dopants. The computational prediction of dopability has a multi-decade history in microelectronic and optoelectronic materials (e.g. III-V compounds, 1 transparent conducting oxides 2,3 ). Successful computational prediction of dopability has been enabled by the accurate description of native defect chemistry, their formation energies and the absolute position of band edges in these materials. 4,5 In thermoelectric materials, computationally guided search for new materials continues to be a strong research focus. [6][7][8][9][10][11] In order to realize the full potential of these materials and to guide the search, accurate predictions of dopability is crucial. 12,13 However, atomic level understanding of native defects in many thermoelectric materials is very challenging due to the presence of heavy elements with strong spin-orbit coupling (SOC). Spin-orbit coupling shifts the band edge positions and significantly reduces the magnitude of the band gap. Describing defect chemistry and dopability in narrow band gap systems such as PbTe is challenging, because small uncertainties in the defect formation energies as well as in the position of the band edges can have strong implication towards predictions of the intrinsic defect conductivity type (n or p), and defect and carrier concentrations. Figure 1 shows a sketch of the formation energy of donor and acceptor defects for a model system with a narrow band gap of 0.2 eV. Error of about 0.3-0.4 eV in the band edge positions (gray boxes) leads to qualitatively distinct conclusions concerning the dopability. In case (a), the system is hard to extrinsically dope p-type due to compensation from native donor defects ('killer' donors). In contrast, predicting case (b) would indicate a highly extrinsically p or n-type dopable system, and case (c) suggests a system with 'killer' native acceptors that limits n-type doping.
Consistent with the challenge presented in Fig. 1, the literature in thermoelectric materials often shows qualitatively different results. [14][15][16][17][18] Here, we consider PbTe a well studied material with strong spin-orbit coupling, to demonstrate the challenges posed in first-principles based density functional theory (DFT) calculations aimed at bulk properties, [19][20][21][22][23][24][25] and intrinsic defects. 14,15,[26][27][28] The DFT calculations employing the local density approximation (LDA) 29 or the generalized gradient approximation (GGA) 30 for exchange correlation provide only a qualitative description of the electronic structure of PbTe. For example, the calculated band gap of PbTe with LDA and GGA is~0.8 eV, 22,23 compared to the experimental value of 0.19 eV at 4.2 K. 31 On the other hand, inclusion of SOC reduces the band gap to nearly zero (−0.01-0.09 eV). 22,23 Issues with the PbTe electronic structure have been addressed 22,24 by using SOC on higher accuracy methods that go beyond the semilocal approximations, such as hybrid functionals, 32 or the GW approach. 33 However, in the context of defect calculations, point defects in PbTe have only been modeled using GGA 14,26 or GGA + SOC. 15,27,28 Bajaj et al., 14 calculated the formation energies of native point defects using GGA, and estimated the equilibrium position of the Fermi energy and the resulting concentrations of free carriers by scaling the GGA band gap (0.82 eV) to match the experimental value. Interestingly, the authors 14 found good agreement between their predicted and experimentallymeasured carrier concentrations even without taking SOC into account in their defect calculations. In another study, Wang et al., 15 performed defect calculations with GGA functional, including SOC effects and found PbTe to be intrinsically n-type under both Te-poor and Te-rich conditions, and found p-type behavior under intermediate Te-rich conditions. This is in contrast to experimentally-observed p-type conductivity under Te-rich conditions and n-type under Te-poor conditions. [34][35][36][37] Given the contested accuracy of defect calculations in compounds with strong spin-orbit coupling, this work focuses on PbTe to establish best practices in defect and dopability calculations. By using hybrid functionals (HSE stands for HSE06 38,39 ) along with spin-orbit coupling to perform defect calculations and quasi-particle GW approach (at G 0 W 0 level) to describe the valence and conduction band edges, we obtain (1) a quantitatively accurate description of the defect chemistry and associated free carrier concentrations in PbTe, and (2) bipolar doping behavior, in agreement with experiments. We also show that all other levels of theory (GGA, GGA + SOC, GGA + SOC + GW, HSE, HSE + SOC) qualitatively fail to describe the experimentally observed self-doping behavior in PbTe. We find that the primary difference between different levels of theory is in the position of the valence and conduction band edges relative to each other and relative to vacuum. Our findings highlight (1) the importance of accurate band edge energies in predicting dopability for systems with strong SOC and (2) that quantitatively accurate dopability predictions are achievable in such challenging material systems.

RESULTS AND DISCUSSION
Defect calculations A realistic description of point defects and intrinsic and extrinsic doping behavior in semiconductors require the knowledge of point defect formation energies, and their resulting concentrations. In this work we employ the standard supercell approach to calculate formation energies of native point defects in PbTe using the following equation: where ΔH D,q represents the formation energy of a defect D in charge state q; E D ,q and E H are the total energies of the supercells with and without the defect, respectively; {μ i } are the chemical potentials of different atomic species describing exchange of particles with the respective reservoirs; E F is the Fermi energy, which in semiconductors ranges between the valence band maximum (VBM) and conduction band minimum (CBM); and E corr is the correction term that accounts for the finite-size corrections within the supercell approach. 40 Next, we describe how each of these terms are computed. Chemical potentials. Numerical values of the chemical potentials (μ i ¼ μ 0 i þ Δμ i ) depend on reference energies, μ 0 i , and enthalpy of formation, ΔH f . The reference energies are obtained from the DFT total energies of bulk Pb and Te metals using respective functionals. Limits to the respective elemental chemical potentials are determined by the thermodynamic stability condition, Δμ Pb + Δμ Te = ΔH f (PbTe), as there are no other competing phases. The enthalpies of formation of PbTe, computed as the difference of the compound total energy and the sum of total energies of Pb and Te in their reference phases, are summarized in Table 1 from each computational method and compared to experiment. Having accurate ΔH f values helps establishing correct equilibrium conditions of different phases, which in turn, allows calculating the limits for the elemental chemical potentials. It is important to note here that the errors in ΔH f of different computational methods fall within the expected errors the first-principles methods typically make. 41 Somewhat surprising is the influence of spin-orbit interactions, which bring ΔH f closer to experiments for both GGA and HSE. It has been shown that the spin-orbit contributions can be seen to a good approximation, as atomic effects that should cancel when calculating total energy differences. 41 Apparently, this does not apply to PbTe and adding spin-orbit interactions helps both in establishing better band edge energies as well as more accurate enthalpies of formation.
Fermi energy. Formation energy of charged defects (q ≠ 0) is a linear function of the Fermi energy as given in Eq. 1. Therefore, the defect formation energies depend on the magnitude of the band gap as well as the energies of individual band edges. The computed PbTe band gap from GGA calculation is 0.79 eV,   The computed band gap with HSE is 1.1 eV, and with HSE + SOC calculation it is 0.29 eV, which is in much better agreement with the experimental value at 300 K. Moreover, HSE + SOC calculations also correctly reproduces the known unusual order of the band gap 0.40 eV E g (PbS) > 0.30 eV E g (PbTe) > 0.27 eV E g (PbSe) within the lead chalcogenides series. 42 As discussed previously, 43 accurate calculation of band gap and band edge energies require GW quasiparticle energy calculations. 33 The band edge shifts from GW calculations are computed relative to the respective underlying exchange-correlation functional. The addition of GW band edge shifts on top of GGA + SOC calculations, opens the band gap to 0.16 eV. For HSE + SOC calculations only a single step G 0 W 0 calculation is performed, resulting in a band gap of 0.194 eV, which is in excellent agreement with the experimental band gap of 0.19 eV at 4.2 K. We do not apply the GW method to defect calculations, as it is computationally too demanding. The GW band edge shits (ΔE VBM , ΔE CBM ) are included in the defect calculations as part of the E corr term through the so called band gap corrections. 43,44 Once the band edge shifts ΔE VBM and ΔE CBM between GGA or HSE calculations and GW are determined, the formation energies are corrected by −z h ΔE VBM (z e ΔE CBM ) for all shallow acceptor (donor) defects occupied by z h holes (z e electrons). This correction accounts for the holes or electrons introduced in the defect calculations occupying VBM or CBM energy levels that need to be shifted to their GW values. 44 To allow comparison across different levels of theory, we reference the respective band edge energies relative to the absolute zero, the vacuum level (Fig. 2). These calculations are performed following the three-step approach 45 for referencing bulk electronic energies to the vacuum, 46,47 involving band edge shifts from GW calculations in combination with DFT calculations of the potential steps normal to the (100) surface of rocksalt PbTe. The computed absolute valence band maxima and conduction band minima are shown in Fig. 2 together with the experimentally measured value of electron affinity in PbTe of 4.6 ± 0.3 eV. 48 Because of the relatively large experimental uncertainty it is difficult to gauge which of the two methods GGA + SOC + GW or HSE + SOC + G 0 W 0 is better in reproducing the measured electron affinity. Certainly, for GW calculations, the specific choice of functional (GGA or HSE) and the level of self-consistency (single step G 0 W 0 or fully self-consistent GW) influences the band edge position. 49,50 Based on our results on native defects chemistry and self-doping behavior in PbTe, we find HSE + SOC + G 0 W 0 band edge energies to be more appropriate than GGA + SOC + GW, thus necessitating the need for more precise experimental measurements.
Finite-size corrections. There are many approaches to compute finite-size corrections in defect calculations. 51 In this work we follow the approach proposed by Lany and Zunger, 40,44 as implemented in our computational framework. 52 The corrections include (1) potential alignment, which restores the relative position of the host VBM in the calculations of charged defect (affecting the Fermi energy), (2) image-charge correction that accounts for the spurious electrostatic interactions of the charged defect with its periodic images, and (3) band filling correction that takes into account the Moss-Burstein-type band filling effects to shallow defects that appear due to high defect concentrations in a typical finite-size supercell calculations. 40 The magnitude of the respective corrections is found to vary between 0.01 to 0.45 eV, depending on the type and charge state of the defect. The computed static dielectric constant (electronic + ionic) needed for image-charge correction is given in Table 1 and is found to be smaller compared to the experimental value, which is in part a consequence of the larger band gaps obtained by the GGA and HSE functionals.
Charge transition levels. Charge transition level or thermodynamic transition level (q 1 /q 2 ) is defined as the Fermi energy for which the formation energies of two charge states q 1 and q 2 of a defect are equal, and is given as where ΔH D,q (E F = 0) is the formation energy of the defect D in the charge state q when the Fermi energy is at the VBM (E F = 0). These levels corresponds to Fermi energy where transition from one defect charge state to another occurs, and are often used as the basis for experimental detection of the defects. 51 The defect formation energies and charge transition levels resulting from the described calculations are shown in Fig. 4 and Fig. 5.
Defect and charge carrier concentrations. The defect formation energy allow calculation of the defect and carrier concentrations in the dilute limit. In this section, we establish set of equations that can be solved self-consistently 3,53 to yield defect concentrations C D,q , and the equilibrium position of the Fermi energy E F . Key to extracting these quantities is the charge neutrality conditions, which sets the position of the Fermi energy and therefore, the corresponding defect and free carrier concentrations. The charge neutrality condition is given as X where q is the charge state, and n and p are the electron and hole concentrations, respectively. The concentration of a defect is obtained by where N is the concentration of the corresponding lattice sites, and k B is the Boltzmann constant. n and p are computed from the density of states (DOS) and the Fermi-Dirac distribution function f (ε) as where D C (ε) and D V (ε) are the DOS at the bottom of the conduction band and top of the valence band, respectively, and using the single parabolic band approximation it is given as For intrinsic defects in PbTe and under non-degenerate doping such that (ε−ε F ) is at least 2 or 3 times k B T, the Fermi-Dirac distribution can be replaced with the Maxwell-Boltzmann distribution function. The carrier concentration integrals in Eqs. 5a and 5b can then be analytically be approximated as The Fermi energy E F is a variable and is determined selfconsistently by solving the charge neutrality condition in Eq. 3 under the convergence criteria of 10 12 cm −3 for concentration for charge balance.
Defects formation energies and charge transition levels In this section, we compare the formation energies and transition levels of native defects in PbTe that result from different levels of theory. Formation energies of vacancy, interstitial and anti-site defects in PbTe under Te-rich conditions are presented in Fig. 4. Calculated with both HSE and GGA with and without spin-orbit coupling, as well as with GW band edge shifts added on top of spin-orbit calculations. We focus on Te-rich conditions (Fig. 4), as they exhibits the largest disagreements between different levels of theory.
The defect plots in the top panel of Fig. 4a and b, obtained using only HSE and GGA functionals are very similar to each other. Both levels of theory overestimate the band gap and, serendipitously, correctly predict the p-type nature of PbTe under Te-rich conditions. The Fermi energy is pinned near the intersection of negatively charged Pb vacancy (q = 2−) and positively charged Te Pb anti-site defect (q = 2+) at about 0.36 eV above the valence band maximum. This will correctly imply p-type intrinsic conductivity, but with relatively low concentration of free holes (Eq. 7b), and more importantly, only a limited extrinsic dopability, which is in disagreement with the experiments. 35,54 Indeed, any attempt to lower the E F below the crossing point using an external acceptor will result in lowering the formation energy of donor Te þ2 Pb , which will compensate extrinsic acceptor and prevent production of free holes. In their GGA study, Bajaj et al. 14 remedied this problem by scaling the defect plots to the experimental band gap of 0.2 eV. As a result, the pining of the Fermi level occurs closer to the VBM (less than 0.1 eV), which, in turn, produces higher hole concentrations. It is important to note that HSE and GGA calculations also predict the correct n-type behavior under Pb-rich conditions, but with low concentration of free electrons (Eq. 7a) compared to experiments.
After including spin-orbit coupling, the band gap reduces significantly and shifts of about 0.1-0.4 eV in the position of charge transition levels occurs between Fig. 4a and c and Fig. 4b and d. Across Fig. 4a-f the major difference is in the position of the band edges and the band gap, and not so much (within 0.2 eV) in the charge transition levels for the majority of the defects. In the HSE + SOC calculations (Fig. 4c), the E F pining occurs below the crossing between V 2À Pb and Te 2þ Pb and is marginally above the mid gap, resulting in low concentration of electrons and a weak n-type behavior. However, in the GGA + SOC calculations (Fig. 4d), the positively charged Te þ2 Pb donor is the dominating defect within the band gap (only 0.02 eV), which pins the E F close to the conduction band edge and results in a high concentration of free electrons under Te-rich conditions. Therefore, neither HSE nor GGA, with or without spin-orbit coupling, is able to satisfactorily reproduce the experimentallyobserved self-doping behavior in PbTe. The resolution of this apparent contradiction is in the observation that the defect formation energies are very similar between different levels of theory if the differences in band edge positions are ignored. To a good approximation, the differences result from different positions of the band edges. As already discussed in Ref. Peng et al. 43 one way of obtaining more accurate positions of valence and conduction band edges is to employ GW calculation, and correct the defect formation energies using the GW band edge shifts, which we discuss next.
The main effect of the G 0 W 0 on the HSE + SOC band structure of PbTe is the shift of the valence band edge by approximately 0.1 eV. This brings the band gap in close agreement with the experiments and affects the defect formation energies by changing the range of allowed E F values. V 2À Pb and Te 2þ Pb remain the low energy defects under Te-rich conditions, but pinning of the E F now occurs below the mid gap, leading to the p-type nature of PbTe. However, addition of GW band edge shifts on top of GGA + SOC calculations in Te-rich conditions, incorrectly results in solely n-type behavior due to Te 2þ Pb , which pin the E F close to the conduction band.
Experimentally PbTe is well-known to be p-type due to Pb vacancies under Te-rich conditions. 34,37 However, in our results we find that only HSE + SOC + G 0 W 0 results correctly reproduce the experiments, without any experimental inputs. The n-type intrinsic behavior in GGA + SOC and GGA + SOC + GW is much more prominent than HSE + SOC, due to singularly dominating donor type Te 2þ Pb defect resulting in the equilibrium Fermi level close to the conduction band. Their failure to predict p-type intrinsic dopability in PbTe under Te-rich conditions, originates mainly from the lower position of the valence band energy, which favors formation of positively charged donor defects. Figure 5, shows defect plots for both Te-poor and Te-rich conditions calculated using HSE + SOC + G 0 W 0 and GGA + SOC + GW. The differences are obvious. Due to lower absolute positions of the band edges in GGA + SOC + GW the donor defects dominate, regardless of the conditions, rendering the systems always n-type behavior. As plots in Fig. 5 assume the most extreme Te-poor and Te-rich conditions (equilibrium with pure Pb and Te, respectively), the system will be n-type for any other set of chemical potentials. The HSE + SOC + G 0 W 0 results are qualitatively different. Under Te-poor conditions the lowest energy defect is donor V 2þ Te for any position of the Fermi energy, rendering the system n-type. Under Te-rich conditions, as already discussed, the equilibrium position of the Fermi energy is determined by the charge neutrality between V 2À Pb and Te 2þ Pb , and the system is correctly predicted to be p-type.
The p-type and n-type intrinsic conductivity in PbTe has been attributed to Pb and Te vacancies, respectively, via experiments. 34,55 Proposition of Pb interstitial as dominating defect over Te vacancies under Te-poor conditions has been made by Schenk et al., 37 and in an earlier work by Brebrick and Allgaier. 56 However, Brebrick and Grubner 34 later corrected their previous conclusions 56 by citing the presence of Cu impurities as donor in Te-poor conditions and not Pb interstitials. In Te-poor conditions, we find Pb 2þ Te , leading to a concentration that is two to three orders of magnitude lower than Te vacancies.
Extrinsic doping is outside the scope of this study, but based on the HSE + SOC + G 0 W 0 results (Fig. 5a), PbTe can be extrinsically doped both n-type and p-type. This is because the energy of native defects is relatively high at the band edges (>0.5 eV), which gives adequate room for formation of extrinsic defects, without having to form compensating defects. On the other hand, the GGA + SOC + GW (Fig. 5b) incorrectly suggests that PbTe will be difficult to ever dope p-type due to the formation of the compensating V 2þ Te and Te 2þ Pb defects, in Te-poor and Te-rich conditions, respectively. 35 Overall, the main differences in predicting self-doping based on defect chemistry between different levels of theory comes from variations of the band edge energies. Our results support similar conclusions made for defect calculations 57-59 comparing defect formation energies and charge transitions levels between standard GGA and hybird HSE functional calculations. We find relatively small differences in atomic relaxations around the defect structures between GGA and HSE, which likely affect the charge transitions levels by about 0.2 eV within the two functionals. Differences of this magnitude may not be notable for a large band gap system but for a narrow band gap system these yield qualitatively different results. However, as discussed by Lyon and Van de Walle, 59 even for large band gap system such as GaN (3.51 eV) differences in atomic relaxations within GGA and HSE can undermine accurate description of optical transition levels.
At this point it is important to discuss the relation of HSE + SOC + G 0 W 0 to the approach frequently used in defect calculations with HSE functionals. The amount of exact exchange, i.e., the α parameter, is often adjusted to match the experimental band gap, and then the same α is used for subsequent defect calculations. 51 This would, in principle, alleviate the need for G 0 W 0 and the associated band edge corrections. Indeed, our HSE + SOC calculations show that for α = 0.18, which gives the band gap of 0.2 eV instead 0.3 eV with the default α = 0.25, the absolute band edges come very close to those of HSE + SOC + G 0 W 0 (see Fig. 2 in the supporting information). Based on our previous discussion this would imply that HSE(α = 0.18) + SOC could be an alternative approach that would offer similarly accurate description of defects in PbTe. While this is likely, one weakness is the need for experimental band gap to tune the value of α. As our primary goal is to construct an approach that is independent of experimental data and applicable to systems not characterized well, we argue that in those cases the choice should be HSE(α = 0.25) + SOC + G 0 W 0 .
Free carrier concentration as a function of temperature In this section, we show that calculations with G 0 W 0 band edge shifts on top of HSE + SOC provides quantitatively correct estimates of the free carrier concentrations in PbTe, whereas other levels of theory fail. Defect and carrier concentrations are computed using the approach already described in Sec. 2.1. To account for the effect of temperature on defect formation energy (Eq. 1), the temperature dependence of band gap (Fig. 6a) as well as band edge energies (Fig. 6b) are taken into account based on the available experimental data. 35,60 The experimental data in Fig.  6c and d is from Hall measurements [34][35][36][37] done between 77-100 K on quenched samples that are annealed at a much higher temperature range, between 400 and 1100 K. Therefore carrier concentrations from the charge neutrality condition (Eq. 3) are computed at 100 K, using defect concentrations corresponding to the annealing temperature of the experiment. Optical experiment measurements 31,42,60,61 have shown that the band gap increases linearly with temperature (at constant pressure) up to 400 K (Fig. 6a), and approaches a constant value of 0.36 eV at higher temperature. The positive sign of the temperature coefficient is interesting since most semiconductors have a negative temperature coefficient. This unusual behavior of the band gap with temperature is explained based on the downward shift of the valence band edge at the L point in the Brillouin zone 31,60 to about 400 K with negative temperature coefficient between 2 and 3 × 10 −4 eV K −1 (Fig. 6b). The conduction band edge moves upward with increasing temperature but its temperature dependence (1.6 × 10 −4 eV K −1 ) is weaker compared to the valence band edge. The temperature dependence of both the conduction band and the valence band edge are adopted from Ref. 35.
The free carrier concentrations computed with the HSE + SOC + G 0 W 0 calculations are not only quantitatively accurate, but also predict the correct bipolar intrinsic conductivity in PbTe, i.e., ptype behavior in Te-rich conditions and n-type behavior in Te-poor conditions. The computed free carrier concentrations are within half to an order magnitude of the experimental values 34,36,37 depending on the temperature, as shown in Fig. 6c and d.
HSE + SOC calculations without the G 0 W 0 band edge shifts, predict the correct n-type conductivity under Te-poor conditions, but fail to predict the p-type conductivity under Te-rich conditions at low temperatures. We find that intrinsic conductivity changes from n-type to p-type at temperatures above 700 K in HSE + SOC calculations under Te-rich conditions. This crossover is due to higher concentration of thermally activated holes compared to that of electrons at high temperatures. The higher effective mass of holes compared to that of electrons in PbTe creates this difference in their concentrations, and because of the narrow band gap of PbTe, easier activation of holes and electron begins to control the conductivity at high temperatures. However, even at high temperatures the computed free carrier concentrations from HSE + SOC are about two order of magnitude lower than experimental values under Te-rich conditions.
Lastly, we raise the question of the possible influence of vibrations to the predictions of defect and carrier concentrations. In narrow gap systems such as PbTe relatively small deviations in defect formation energies and/or band edge positions could come from vibrations and could have a profound effect on the predictions. This is certainly a valid point that warrants further investigations of the phonon contributions to defect formation energies, especially in relatively soft systems like PbTe in which the occupation of phonon modes can be significant even at moderate temperatures, and which also exhibit a significant amount of anharmonicity. 23,25 While the defect calculations in the presence of phonons are beyond the scope of the present study, our results clearly show that by neglecting phonon contributions to defect formation energies defined in Eq. 1 it is possible to obtain qualitatively correct and quantitatively reasonably accurate description of defect chemistry and intrinsic carrier concentration in PbTe. This implies that the phonon contributions amount to differences between the predicted and experimentally measured carrier concentrations, that are estimated to be roughly of an order of magnitude. Despite the challenges, our results show that it is possible to accurately predict the intrinsic dopability in PbTe from firstprinciples point defect calculations and thermodynamic simulationss. These results have important implication towards the other thermoelectric materials such as PbSe, PbS, and Bi 2 X 3 (X = S, Se, Te). These systems have narrow band gaps due to strong spin-orbit interactions. However, an accurate description of point defects in these systems is still missing. Intrinsic dopability controls the effectiveness of external dopants, as native point defects with low formation energy can act as compensating 'killer' defects. Therefore, using the computational methodology described in this paper will not only enhance the understanding of point defects, but also help in making accurate prediction about dopability in these materials.
In conclusion, we have systematically reviewed first-principles calculations of native point defects in PbTe using different levels of theory, with the goal of achieving qualitatively correct and quantitively accurate description of its intrinsic defect chemistry and self doping behavior. Similar to other narrow band gap systems, achieving accurate atomic-level description of intrinsic defect chemistry is challenging because small uncertainties in the calculations can lead to large deviations in the predictions. We showed here that accurate description of the experimentally observed bipolar doping behavior as a function of the synthesis conditions and measured charge carrier concentrations, requires defect calculations that combine hybrid functionals with spin-orbit coupling included and quasi-particle G 0 W 0 description of the position of the individual band edges. The primary difference between different levels of theory considered here (DFT and Hybrid functionals, with and without SOC and GW band edge shifts) is actually in the position of the valence and conduction band edges relative to vacuum. The correct description of intrinsic defect chemistry and self-doping emerges only if total energies are calculated at the HSE + SOC level and the results are corrected using the G 0 W 0 band-edge shifts. Our results reaffirm the importance of the band edge positions in defect calculations and help formulate reliable first-principles procedure for predicting dopability in PbTe and other narrow band gap systems.

METHODS
DFT calculations are performed using the projector augmented wave (PAW) method 62 as implemented in VASP. 63 The plane wave energy cutoff of 340 eV, and a Monkhorst-Pack k-point sampling 64 is used. The bulk properties (Table 1) are calculated using a 12 × 12 × 12 and 8 × 8 × 8 kpoint mesh for the primitive cell with GGA and HSE calculations, respectively. The static total (electronic + ionic) dielectric constant with GGA functional is calculated via the density functional perturbation theory 65,66 (DFPT). With HSE functional, the electronic component to the dielectric tensor is calculated from the self-consistent response of the system to a finite electric field, 67 and the ionic part is taken same as that of the GGA functional from DFPT calculation. Defect calculations are performed on 64 atom bulk supercell with a Γ-centered 2 × 2 × 2 k-point mesh. Static self-consistent SOC calculations are performed on GGA and HSE relaxed defect structures. In the present VASP implementation of the GW method, 63 the inclusion of SOC is not possible. Therefore, the GW band edge shifts are estimated relative to the GGA or HSE computed band edges, and subsequently added to the SOC results.
To test the convergence of our results with supercell size we have performed additional calculations using DFT-PBE on a 216-atom supercell. The defect formation energies are converged within 0.1-0.2 eV depending on the type and charge state of the defect, and the net carrier concentrations are converged within an order of magnitude or less depending on the temperature. The overall description of the defect chemistry remains essentially the same, and the position of the equilibrium Fermi energy is well converged within 0.02 eV between the 64 and 216atom supercell calculations. Results from 216-atom supercell calculations are summarized in Fig. 1 in the Supporting information.

Data availability statement
The authors declare that the data supporting the findings of this study are available within the paper and in its supplementary information file.