Correcting the corrections for charged defects in crystals

While the theory of imperfections in solids is firmly established, procedures for first-principles calculations of defect quantities continue to evolve. A plethora of ad hoc correction schemes is being replaced by sophisticated self-consistent procedures that will enable more quantitative predictions of the formation energies of defect species and their spectroscopic signatures.

The formation of defects in crystals is a natural consequence of chemical thermodynamics-the balance between the enthalpic cost of perturbing the atomic bonding environments and the entropic gain of introducing an ensemble of imperfections.
The key quantity that determines, in equilibrium, if a particular defect species will be abundant or rare is the free energy of formation, ΔG f . The concentration (n d ) can be decomposed into contributions from enthalpy and vibrational entropy where N site and g denote the number and degeneracy of available sites in the host crystal. The enthalpy change dominates under standard conditions, so the more burdensome vibrational term 1 is most often neglected. The first objective in defect modelling is to calculate ΔH f as accurately as possible, which can be achieved in a number of ways, as illustrated in Fig. 1. Due to the exponential factor in Eq. (1), the prediction of defect concentrations is sensitive to the quality of the underlying energy terms. My brief perspective on major developments in this field is given below. I apologise in advance to the many contributions that are overlooked in this synopsis (e.g. see ref. 2 for a more thorough overview).

EARLY 20TH CENTURY-THEORETICAL FOUNDATIONS
The theory of lattice dynamics developed by Born and others described the motion of atoms around their equilibrium crystallographic positions. However, what if atoms left its ideal position and wandered into an interstitial position? A classical numerical procedure for calculating the formation energy of such point defects in ionic solids was proposed by Frenkel in 1926 3 . The semiclassical (continuum) Mott-Littleton method for charged defects energetics was reported in 1938 4 , which was later expanded into a more general framework for multi-region embedded crystal calculations 5 . The extensive theoretical infrastructure subsequently developed for describing the electronic structure of donor and acceptor levels in semiconductors was reviewed by Pantelides in 1978 6 .

LATE 20TH CENTURY-PRACTICAL SOLUTIONS
The combination of efficient algorithms to solve the Kohn-Sham equations and massively parallel computing enabled the modern era of computational materials science. Two types of approaches were developed for modelling charged point defects in crystals.
Those based on embedding potentials (a dilute defect in a host matrix as shown in Fig. 1c) offer some advantages, but remain technically challenging to set up and analyse 7 . A self-consistent Green's function procedure for treating defect perturbations was already reported in 1979 8 and quickly adapted to describe the deep states associated with transition metal impurities in Si using a local density functional 9 .
The alternative, and more widely employed, supercell approach benefits from the robust infrastructure for calculation of crystals within periodic boundary conditions (Fig. 1d). The principal issue in describing charged defects is the electrostatic interaction between repeating centres. The long-range Coulomb interaction depends on the defect charge q, their spatial arrangement, and the dielectric response of the host ϵ(r). The standard solution is to introduce a homogeneous 'jellium' background charge to enforce charge neutrality in each repeat unit and ensure convergent Coulomb energy. This fix results in a shift in the average electrostatic potential and total energy of the defective supercell that must be corrected to become useful (Fig. 2).
Leslie and Gillan 10 showed how to properly account for a point charge q interacting with its periodic images through an isotropic dielectric medium. The approach was expanded by Makov and Payne 11 to include the quadrupole moment (Q). In practice, Q cannot be determined from the defect charge distribution. Defects come in many forms and their wavefunctions may be localised or delocalised and are invariably coupled to the screening host charge. In addition, the limitations of a given density functional can require further correction terms, such as valence and conduction band edge shifts, to reduce the divide between modelling and measurement 12 .
Exciting results and insights were gained using these methods, especially in the area of metal oxides. There were problems, however. Research groups performing similar simulations would obtain different conclusions depending on the flavour of corrections that were employed; the curious reader can look into the case of defects in ZnO. In my opinion, the absolute defect energetics reported during this era should be taken with a pinch of salt.

EARLY 21ST CENTURY-ROBUST PREDICTIONS
To tackle these issues, a new set of correction schemes emerged for supercell calculations. Freysoldt et al. 13 modelled the defect charge as a Gaussian distribution in an isotropic medium. There were other efforts to include anistropic dielectric screening 14,15 . Kumagai and Oba 16 refined these approaches with a more practical alignment procedure, based on atomic site potentials, that accounts for the full anisotropic low-frequency dielectric response.
More recent developments have moved away from a posteriori corrections of the total energy to direct modification of the underlying self-consistent calculation. The advantages are that the total energy and electronic eigenvalues can be corrected directly and more physical long-range dielectric screening can be incorporated.
A self-consistent potential correction was proposed by da Silva et al. in 2021 17 . The change in potential, from either a reference pristine supercell or neutral defect, is used to correct for the difference between the periodic and isolated charged defect. One limitation in the current formalism is the use of an isotropic dielectric constant; however, spatial variation is allowed in one direction for the case of slabs. In contrast, the image charge correction proposed by Suo et al. 18 avoids the input of a dielectric constant at all and is instead based on the self-consistent charge density difference between charged and neutral defects, which already contains the relevant screening information. Both approaches show promising behaviour for a series of test cases including NaCl and MgO.
An issue not addressed in these two recent works is how to incorporate the low-frequency dielectric response involving ion displacements. This is critical for accurately describing defect transitions, which may be excited optically or thermally 19 . Suo et al. 18 suggest introducing a Gaussian broadening scheme to capture ion displacements, but this remains to be implemented and tested. An entirely different approach has been taken by Xiao et al. 20 in a move away from the standard jellium model to consider charge compensation by realistic valence or conduction band edge states of the host crystal. Their self-consistent correction scheme, which avoids the artificial jellium background, is intuitive for traditional semiconductors where charged defects are compensated by electrons/ holes, but less so for cases where ionic compensation dominates 21,22 . The approach gives good agreement with conventional methods and can also describe defects in low dimensional structures.
This new wave of research and development in the field is exciting and will further increase the predictive power of firstprinciples calculations of materials. A family of reliable and general correction schemes for charged defects will allow us to confidently tackle important scientific challenges such as equilibrium defect distributions, defect vibrations, and non-equilibrium charge transitions. They will also support more robust workflows for defect automation that can be used to identify new behaviour and physical trends in imperfect crystals.