First-Principles Modeling of Defects in Lead-Halide Perovskites: Best Practices and Open Issues

Lead-halide perovskites are outstanding materials for photovoltaics, showing long lifetimes of photo-generated carriers which induce high conversion efficiencies in solar cell and light-emitting devices. Native defects can severely limit the efficiency of optoelectronic devices by acting as carrier recombination centers. The study of defects in lead halide perovskites thus assumes a prominent role in further advancing the exploitation of this class of materials. The perovskites defect chemistry has been mainly investigated by computational methods based on Density Functional Theory. The complex electronic structure of perovskites, however, poses challenges to the accuracy of such calculations. In this work we review the state of the art of defects calculations in lead halide perovskites, discussing the major technical issues commonly encountered and what we believe to be the best practices. By keeping as a test case the prototype MAPbI 3 compound, we discuss the impact of exchange-correlation functional on the electronic structure and on the defect The semi-local are compared and with results the functionals (HSE06 and PBE0) as a function of the amount of exact exchange, with and without including SOC. This analysis shows that sizable deviations in the thermodynamics and electronic properties of defects in perovskites are obtained at different levels of theory. The large deviations in the electronic properties of MAPbI 3 obtained with semi-local and hybrid functionals also induce remarkable energies differences between the calculated lattice defects equilibrium geometries.

formation energies (DFEs), by comparing semi-local and hybrid functionals, with and without spinorbit coupling corrections. The convergence of calculated defects structures and their DFEs with respect to the simulation supercell size and the performance of commonly employed charge corrections is thus discussed. A summary of results concerning the defect chemistry of MAPbI 3 is provided, providing hints onto the so called defect tolerance of this material class and casting possible scenarios to be addressed by experimental investigations.

TOC graphics
Metal-halide perovskites are among the most promising materials for photovoltaics 1 due to their attractive optoelectronic properties, such as direct band gap, 2, 3 large absorption coefficient and low recombination rate of photo-generated charge carriers [4][5][6][7] coupled to a fairly high mobility. 8 The high efficiencies of these materials in solar cells, exceeding 22%, are mainly related to the very long lifetimes of photogenerated carriers, [8][9][10] even in the presence of moderate-high defects density. [11][12][13] As for any semiconductor, native defects assume a prominent role in the optoelectronic properties of perovskites due to their ability to trap photo-generated carriers and to promote their recombination on the defects sites, leading to solar cell efficiency loss, which is especially manifested as open circuit voltage losses. 14 According to the Schokley-Read-Hall theory, 15 the charge trapping rate R of carriers on a defect center is defined by = − 2 ( + 1 ) + ( + 1 ) (1) where n, p, and n i are the electron/hole carrier densities and intrinsic carriers density, respectively; τ p,n are the holes and electrons lifetimes; 1 = exp (− − ) and 1 = exp (− − ); and E T is the energy of the considered trap (E C and E V are the conduction and valence band energies, with related N C and N V density of states). Based on Eq. 1, trap states with ionization levels placed in the middle of the band gap show an increased trapping activity compared with shallow charge traps. 15 Despite the limited lattice stability of perovskites compared to traditional semiconductors, only low to moderate trap densities of ~10 11 -10 16 cm -3 for single crystals and polycrystalline samples have been reported, respectively. [16][17][18][19][20] Such densities are comparable to those usually reported in conventional semiconductors like polycrystalline Si (10 13 -10 14 cm -3 ), 21 CdTe (10 13 -10 15 cm -3 ) 22 and CIGS (10 13 cm -3 ). 23 Nevertheless, the efficiencies of perovskites in solar cells are only slightly affected. 17,24 Although the high dielectric screening of the organic/inorganic lattice can increase the lifetimes of charge carriers through the formation of large polarons, [25][26][27][28] the origin of the apparent defect tolerance in this class of materials is still under debate and several studies aimed to unveil the nature and the photophysics of charge traps in lead-halides perovskites have been performed. 12, 13, 16-19, 29-42, 46- showing that the Pb sublattice is also largely involved in the defect chemistry, by forming vacancies in the lattice, with potentially harmful impact on the overall PL properties. 32 Spatially resolved PL experiments performed on polycrystalline films of MAPbI 3 showed that I-rich surfaces are less emissive than Pb-rich surfaces, suggesting that most of the charge traps forms in I-rich portions of the perovskite. 33,34 This result is related to the study by Minns et al., who showed that high structural disorder in bulk MAPbI 3 is associated to iodine interstitials, through X-ray and neutron diffraction techniques. 35 The charge carriers dynamics has been also investigated in order to figure out the origin of the apparent defects tolerance. PL and transient absorption experiments performed by Leijtens et al. confirmed the existence of electron traps in polycrystalline MAPbI 3 , with higher densities compared to single crystals, i.e. ~10 16 cm -3 . 17 Further insights in the photophysics of defects have been also provided by the same authors, by showing that after the initial trapping of electrons on the defects centers, filled traps become essentially inactive for very long times, of the order of μs. This largely limits the non-radiative recombination on these centers and leads to a predominant hole conductivity. 17 Indirect hints concerning the nature of charge traps in lead-halide perovskites have been provided by investigating PL properties under different environment conditions. Exposure of MAPbI 3 to low pressures of gaseous I 2 was found to induce a quenching of the PL quantum yield (PLQY), 36 accelerating the degradation of the perovskite through the formation of metallic lead, as revealed by XPS experiments. 37 An irreversible p-doping of the perovskite associated to a decreasing of ionic conductivity was also reported. 38 MAPbI 3 and avoiding trapping on single charged defects, as discussed by Walsh and coworkers. 64 The large part of these studies, however, have been carried out by using semi-local exchange correlation functionals, such as the Perdew-Burke-Ernzherof (PBE) functional. 65 The use of such semi-local functionals, although usually providing a good approximation to defect geometrical structures and lattice energies, may lead to unsatisfactory results for the calculations of defects formation energies (DFEs) and thermodynamic ionization levels when both spin-orbit coupling (SOC) and self-interactions corrections are neglected. 66  reported, with connection to available experimental observables.

Defects formation energies: computational protocol
The stability of native defects in solids is commonly evaluated by calculating the associated DFE through DFT calculations in the supercell approach following the equation: [83][84][85] ( ) = ( ) − ( ) − ∑ + ( + + ∆ ) + where E(X q ) is the energy of the defect supercell of charge q, E(bulk) the energy of the pristine supercell, n and µ are the number and the chemical potentials of the species added or subtracted to the perfect bulk in order to form the defect. The fourth term represents the energy associated to the exchange of charges with the electrons reservoir (the Fermi level of the system E F ), referenced to the valence band maximum (VBM) of the pristine crystal. Adopting this convention, the energy of the VBM explicitly contributes to the DFE. The absolute VBM energy should also be corrected for the electrostatic potential shift ΔV induced by the defect in the supercell calculation compared to the bulk reference potential, from which the VBM is taken. Finally, the last term in Eq. 2, E q , is introduced to correct the electrostatic interaction between charged point defects in the periodically repeated supercells.
As expressed by Eq. 2, DFEs are function of the chemical potentials of the elements constituting the lattice µ i (i = MA, Pb, I) and of the Fermi level of the system E F . By fixing the chemical potentials, DFEs show a linear dependence on the Fermi level in the system which is allowed to vary within the material band gap, whose slope is determined by the charge of the defect.
Where DFE(D q , E F =0) are the formation energies of defects in state of charge q calculated for E F = 0, i.e. coinciding with the VB maximum (VBM). 83 The relative concentration c(D q ) of defects can be estimated by the calculated DFE by applying the Boltzmann law ( ) = exp (− ( , ) ), where N s is the density of the defects sites in the modelled supercell, T is the absolute temperature and k the Boltzmann constant. Furthermore, from the calculated defect density it is possible to predict the position of the native Fermi level in the system, thus whether the material will be n-or p-doped (or intrinsic) by solving the associated electro-neutrality equation: 83 where p and n are the intrinsic densities of electrons and holes.
Based on the above framework, three terms affect the accuracy of calculated DFEs and charge transitions: i) the relative energies of the bulk and defective systems; ii) the electronic term accounting for the exchange of charges between the defect site and the Fermi level, i.e. the absolute VBM value of the bulk reference phase; and iii) the electrostatic corrections need to account for the presence of charged defects. In the following, our analysis will focus on the computational methods applicable in order to accurately predict these terms.

Stability and electronic properties of bulk MAPbI 3
To simulate the defect chemistry of MAPbI 3 in different environments, appropriate values of the chemical potentials to the constituting species should be assigned. The energy ranges of chemical potentials are set in order to preserve the thermodynamic stability of the perovskite, by simulating halide rich (metal poor) or halide poor (metal rich) conditions. The field of thermodynamic stability of MAPbI 3 is defined by the following relations for the chemical potentials μ: These relations express the stability of the three relevant compounds in equilibrium, i.e. MAPbI 3 with its constituents (Eq. 5), MAI (Eq. 6) and PbI 2 (Eq. 7); they are commonly expressed by referring the chemical potentials of the species to their respective values in the standard state. While the definition of standard states of iodine and lead is straightforward, by using the chemical potentials of molecular iodine and of bulk metallic lead, a difficulty clearly arises in assigning a definite standard state to the charged MA molecule. 12 The definition of a MA standard state is however not essential in the simulation of I-rich or I-poor conditions and its value can be evaluated by solving Eq. 5-7. In Table 1 the expression of chemical potentials simulating MAPbI 3 in equilibrium with PbI 2 and MAI in I-rich and I-poor conditions are reported for clarity. with PbI 2 and MAI phases. Points A and C in the thermodynamic diagram are showed in Figure 1.

I-poor μ(I) = [μ(MAPbI 3 )-μ(MAI)-μ(Pb bulk )]/2 μ(Pb) = μ(Pb bulk ) μ(MA) = [3μ(MAI)-μ(MAPbI 3 )+ μ(Pb bulk )]/2
In Table 1 halide rich conditions are simulated by setting the chemical potential of I to that of the I 2 molecule (I 2 mol ), while halide poor conditions can be simulated by setting the chemical potential of Pb to that of the metallic lead, Pb bulk . Although the standard state of iodine at room temperature is the I 2 molecule in the orthorhombic solid phase, we consider hereafter as a standard state the iodine molecule in the gas phase. This assumption is justified by considering the high vapor pressure of iodine at the temperature usually employed in the synthesis and annealing of MAPbI 3 , i.e. above 100 °C.
By calculating chemical potentials according to expressions in Table 1, it is possible to define the thermodynamic field of stability of MAPbI 3 , as illustrated in Figure 1a. Notably, the perovskite is stable only in a limited area of I and Pb chemical potentials, in equilibrium with the PbI 2 and MAI phases, whose width is the enthalpy of formation of the perovskite starting from PbI 2 and MAI, i.e. -0.11 eV as calculated at the PBE level. Similarly important to the rational definition of the chemical potentials for the elements, a quantitative prediction of the band edges is essential in order to estimate accurately the electronic term in Eq. 2, since a possible inaccuracy in this term can have dramatic effects on the calculated DFEs.
Semi-local exchange correlation functionals, such as PBE 65 and PBEsol, 86 Figure 1b). However, as previously discussed, such matching is only fortuitous. 87,88 The inclusion of SOC strongly narrows the band gap to 0.55 eV (see Figure 1b). 89,90 Notably, this is the correct band-gap predicted by semi-local functionals, retrieving the typical band-gap underestimate known for most of semiconductors. By using GW calculations including SOC on tetragonal MAPbI 3 , a band gap of 1.64 eV was obtained, 87  providing a good description of the electronic properties, is currently limited to the calculation of quasi-particle energies at fixed geometries, although promising procedures aimed to combine the GW addition/removal electrons energies with lattice relaxation energies obtained by DFT methods have been recently proposed. 91 The high computational cost of the GW technique, however, currently limits its application to small/medium systems. A good compromise for an accurate treatment of both lattice and electronic energies is represented by hybrid exchange correlation functionals which employ a fraction of the exact Hartree-Fock exchange in their formulations in order to partially correct the self-interaction error of semi-local functionals. The advantage of using hybrid functionals, such as PBE0 92 and HSE06 68 is that, beside a quantitative estimation of the electronic band edges energies, they also provide accurate lattice energies, generally improving upon semi-local functionals. The fraction of exact exchange (α) included in the hybrid functional is a tunable parameter and is set to a value of 0.25 for both PBE0 and for the screened HSE06 functional, an optimal value for low-medium band gap semiconductors. 93 Table S1, Supporting Information, where a list of properties of MAPbI 3 , PbI 2 and I 2 have been compared as calculated by using two different values of exchange fractions. Notice that while in PBE0 the exact exchangemultiplied by α -is used for any inter-electronic distance, in HSE06 the same exchange terms is used within a radius of inter-electronic distances, to switch to the semi-local exchange term at long range. As we will show below, this has important consequences for the calculation accuracy.
In the calculations of solids it is customary to remove core (and semi-core) electrons from the actual electronic structure problem, for computational convenience. This is either afforded by replacing the core electronic states by a pseudo-potential (PP) 94,95 or by freezing such states at those of the isolated atom, as in the Projector Augmented Wave (PAW) approach. 96 Notably, the use of PPs can be more computationally effective, but we found that PPs in combination with hybrid functionals can lead to unsatisfactory results in hybrid calculations if partially incomplete valence shells are considered in the pseudopotentials. This is (unfortunately) the case of Pb when the semi-core n=5 valence shell is excluded in the pseudopotential construction, leaving only the 5d, 6s and 6p electrons in the valence, the most common pseudopotentials usually employed in calculations.
The exclusion of the semi-core 5s and 5p electrons in the calculation of exact exchange leads to an unphysical down shift of the CB band edge by increasing the fraction of exchange with large underestimations of the calculated band gaps, see Figure S1,

Effects of self-interaction and spin-orbit coupling on defect properties
The terms affecting the accuracy of calculated DFEs in Eq. 2 are the energies of the bulk and defective supercells and the electronic term accounting for the charges exchange with the Fermi level. In this section we discuss the influence of the self-interaction error and spin-orbit coupling on these terms.   Table 2.  By the analysis of results in  An overview of the thermodynamic ionization levels of native defects in MAPbI 3 at different levels of theory are reported in Figure S2, Supporting Information.
It is worth mentioning that HSE06 functional provides more accurate results than PBE0 for thermodynamic ionization levels of shallow defects. Inaccuracies of the PBE0 functional in the prediction of shallow transition levels have been reported in the literature. 99 We illustrate this point in Figure 3d, where the thermodynamic ionization levels of the pristine phase P*, iodine interstitial The impact of dispersion corrections on the calculated DFEs is also related to the choice of chemical potentials references. As a comparison, in Figure S3 of the SI we report the effects of dispersion corrections on the DFEs of MAPbI 3 in I-rich conditions by using as chemical potential reference for iodine the I 2 gas molecule and the I 2 in the orthorhombic solid phase. In the last case the impact on the DFEs is less dramatic than in the former and an overall trend to stabilize/destabilize interstitials/vacancies is observed.

Charge corrections and elastic effects
Within the supercell approach the simulation of charged defects is limited by periodic interactions of the charges in periodic boundary conditions (PBC). The calculation of charged systems in PBC is problematic due to the divergence of the G=0 term in the Hartree potential. This term is usually In the Makov-Payne approach the correction term is an extension of the Madelung energy for localized charge q in a periodically repeated cubic lattice of length L, whose expression is   Table 3 the DFEs of a selected list of neutral defects calculated at the PBE level in I-rich conditions have been reported as calculated in the 2x2x1 and 2x2x2 supercells. While most of the defects show converged values already in the 2x2x1 supercell, in the most critical case, the I i interstitial is stabilized by ~0.31 eV going from the 2x2x1 to 2x2x2 supercell. This suggests that in order to avoid elastic effects in the evaluation of DFEs, calculations should be performed at least in 2x2x2 supercells of the tetragonal phase or in comparable supercells for cubic and orthorhombic phases. Based on these considerations, in this last section we provide an overview of the main defects properties of the MAPbI 3 perovskite, as calculated in the 2x2x2 tetragonal supercell by using the HSE06-SOC functional (α=0.43) and by only applying the potential alignment term.

Defect chemistry of MAPbI
Dispersion corrections have been also added a-posteriori according to the DFT-D3 method of Grimme. 104 In Figure 5a-c we report the calculated DFEs in I-rich, I-medium and I-poor conditions; the thermodynamic ionization levels diagram is shown in Figure 5d; and the optimized structures of a selected list of defects are reported in Figure 5e.

Effects of the environment
The calculated DFEs and ionization levels also deepens our knowledge of the defects chemistry when the perovskite interacts with the environment, particularly regarding the activity of charge traps. As an example, in Figure 6, a schematic of two possible mechanisms of PL quenching (PLQ) and enhancement (PLE) induced by exposure of MAPbI 3 to I 2 and O 2 gas, respectively, is reported. 80,81 The exposure of MAPbI 3 to moderate pressure of I 2 gas induces: i) the formation of metallic conditions (Figure 6c), simulating MAPbI 3 in equilibrium with I 2 gas at room temperature, also confirms this picture by indicating that most stable defects in these conditions are I i + and V Pb The exposure of MAPbI 3 to moderate quantities of dry oxygen, on the other hand, leads to an enhancement of the PL (PLE), as reported by several experiments. [43][44][45][46][47][48][49] The reversibility of the process suggests a weak interaction of O 2 with the charge traps of the perovskite, likely leading to their passivation. The interaction with O 2 however is a complex process, and can also induce a degradation of the perovskite through the formation of superoxides species assisted by light. 82 The beneficial effect of oxygen on the PL can be investigated by simulating its interaction with charge traps, in particular negative iodine interstitials which can trap holes through the (0/-) transition, leading to radiative recombination processes. The oxidation of negative iodine interstitials I i mainly leads to the formation of typical hypoiodite, iodite and iodate compounds, i.e.
IO -, IO 2 and IO 3 -. The formation of these species is only slightly thermodynamically favored, with calculated PBE energies of 0.02, 0.13 -0.05 eV, respectively, highlighting that such reactions are partially reversible. 81 Furthermore, besides the oxidation of interstitials, oxidation of the lattice iodines can also take place, although this process is less favored than the former. 81 Among the aforementioned reactions, the formation of iodate, i.e. I i -+ 3/2 O 2 = IO 3 -(see Figure 6e) is the most favored reaction and likely IO 3 is the most stable defect product forming upon exposition of MAPbI 3 to O 2 -rich environment. The effect of the oxidation of interstitials defects on the PL properties can be investigated by the analysis of the deep (0/-) transition responsible of the charge trapping on the center. The oxidation of iodine interstitial to form IO n induces a shallowing of the (0/-) transition, that shifts down 0.25 eV toward the VB as reported for the IO 3 defect (see Figure   6f), by reducing the trapping activity of the defect. 81 By our model thus we predict a reversible passivation of iodine traps through the formation of iodate species in the lattice associated to an increase of the PLQY of the perovskite.

Concluding remarks
In summary, an analysis of the technical issues commonly encountered in the DFT modelling of defects in lead halide perovskites has been carried out. Our discussion, focused on the prototype We also discussed convergence issues in the calculation of the defects formation energies, by comparing the performance of commonly employed correction schemes. In this regard, the high dielectric screening of the lattice strongly reduces the periodic interactions of charged defects in the supercell approach and potential alignment procedures are sufficient in order to have accurate DFEs, provided that calculations are performed in medium/large supercells.
Finally, we discussed the defect properties of MAPbI 3 by using a converged computational strategy. The results show that the defect chemistry of MAPbI 3 is dominated by the typical iodine redox chemistry, as exemplified by the defect properties of interstitial iodine. As such, we expect quantitative differences by changing the halide in e.g. MAPbBr 3 , although the underlying chemistry should remain similar. The lack of strong directional preference for lead-halide bonds makes Pbbased defects less relevant for trapping.
We believe this study will serve as useful reference for further defect calculations on leadhalide perovskites.

Computational details
All defects calculations throughout the work have been carried out in the tetragonal phase of the Fock grid has been reduced to 80 Ryd, without sizable impact on the total energies and preserving the good convergence of calculations. 1x1x2 and Γ-only k-points grids in Brillouin zone were used for the 2x2x1 and 2x2x2 supercells, respectively. All calculations have been carried out by using the Quantum Espresso package. 106 Dispersion corrections have been applied to defects formation energies a-posteriori on the PBE geometries by using the DFT-D3 code of Grimme. 104 The 22 electrons norm-conserving pseudopotential used in hybrid calculations has been generated by using the ONCV code, 107 while FNV corrections to the DFEs have been applied by using