Defect Formation Energies without the Band-Gap Problem: Combining DFT and GW for the Silicon Self-Interstitial

We present an improved method to calculate defect formation energies that overcomes the band-gap problem of Kohn-Sham density-functional theory (DFT) and reduces the self-interaction error of the local-density approximation (LDA) to DFT. We demonstrate for the silicon self-interstitial that combining LDA with quasiparticle energy calculations in the G0W0 approach increases the defect formation energy of the neutral charge state by ~1.1 eV, which is in good agreement with diffusion Monte Carlo calculations (E. R. Batista et al. Phys. Rev. B 74, 121102(R) (2006), W.-K. Leung et al. Phys. Rev. Lett. 83, 2351 (1999)). Moreover, the G0W0-corrected charge transition levels agree well with recent measurements.

Defects often noticeably influence the electrical and optical properties of a material by introducing defect states into the band gap. Reaching a microscopic understanding of the physical and chemical properties of defects in solids has long been a goal of first-principles electronic structure methods. Probably the most widespread theoretical method in this realm today is density functional theory (DFT) in the local-density (LDA) and generalized gradient approximation (GGA), but certain intrinsic deficiencies limit their predictive power. Artificial self-interaction and the absence of the derivative discontinuity in the exchange-correlation potential [1] present the most notable deficiencies in this context. They give, amongst other things, rise to the band-gap problem -the fact that the band gap in LDA and GGA underestimates the quasiparticle gap [1,2]. In this Letter we show that the band-gap problem in LDA/GGA not only affects the reliable computation of defect levels, but in certain cases (e.g. filled defect states in the band gap) also that of formation energies. We present a formalism for calculating formation energies of defects in solids that combines LDA with quasiparticle energy calculations in the G 0 W 0 approximation [3] to reduce the self-interaction error and to overcome the band-gap problem. In some cases a heuristic "scissor operator" approach may approximately correct the problem. However, particularly when the experimental answer is unknown, a more accurate method is needed.
We illustrate our approach with the example of a selfintersitital in silicon (Si i ), a defect of high technological relevance [4,5,6]. In the neutral charge state the Si i has several stable and metastable atomic configurations [7,8] (see Fig. 1), in all of which two electrons occupy a defect level in the band gap. The LDA formation energies of all these configurations are underestimated severely (by ∼1.5 eV) compared to diffusion Monte Carlo (DMC) calculations [9,10]. However, no insight into this dis- crepancy is provided by the DMC calculations.
In our formalism the formation of the neutral defect is expressed as successive charging of its 2+ charge state, for which the defect level is unoccupied. This allows us to decompose the formation energy into that of the 2+ state (E f (2+)), a lattice and an electron addition part. This decomposition is not only insightful for analyzing the underestimation of the LDA formation energy, but also permits us to apply the most suitable method for each type of contribution [11]. For the lattice part we retain the LDA and argue that the relaxation energies and E f (2+) are not as strongly affected by the deficiencies of the LDA as in the positive and the neutral case since the defect level in the band gap is unoccupied. For the electron affinities, on the other hand, we employ Hedin's GW method [3]. Since self-consistency in GW is still discussed controversially [2] we obtain the quasiparticle corrections arXiv:0812.2492v1 [cond-mat.mtrl-sci] 12 Dec 2008 to the LDA Kohn-Sham energies from first order perturbation theory (G 0 W 0 ), which is currently the method of choice for computing quasiparticle band structures in solids [2,12]. While not completely self-interaction free [13] G 0 W 0 significantly reduces the self-interaction error. With this combined approach the formation energy in the neutral charge state increases by ∼1.1 eV compared to the LDA. Recent DFT calculations employing the HSE hybrid functional, which also significantly reduces the self-interaction error, yield a similar improvement [10] and lend further substance to this notion. Moreover, the G 0 W 0 -corrected charge transition levels agree well with recent experimental measurements [4].
We will present our combined DFT+G 0 W 0 approach for the example of the Si i , but it can easily be generalized to defects in other materials. An additional silicon atom in an interstitial site can adopt different configurations with similar formation energies (cf. Fig. 1). In the tetrahedral (tet) geometry the extra silicon atom gives rise to an a 1 and a threefold degenerate t 2 state. The latter is empty and in resonance with the conduction band. The partial occupation of the degenerate t 2 state triggers a Jahn-Teller distortion along the <111>-axis into a geometry with C 3v symmetry, also referred to as "displaced hexagonal structure" in previous studies [14]. The addition of a second electron displaces the atom further in the <111>-direction stabilizing the neutral charge state. This sequence is illustrated in Fig. 2. Moving the interstitial atom along the <111>-direction into the center of a six-membered ring (hexagonal (hex) configuration) lowers the neutral state further in energy. It reaches its lowest position in the split<110> configuration, where the added atom and a host atom share an interstitial site oriented in the <110>-direction.
For the 2+ charge state the tetrahedral is the most stable configuration [8] (see also Table I) and we will use it as a starting point for building our scheme. The positive charge state is then formed by adding one electron as depicted in steps 1 and 2 in Fig. 1. Mathematically this can be expressed by starting from the expression for the formation energy in the positive charge state is the total energy in charge state q and atomic positions R D q of defect configuration D in charge state q . E ref is a suitably chosen reference system, here bulk silicon, and F the Fermi level of the electrons referenced to the top of the valence band. Adding and substracting first E(+, R tet 2+ ) and then E(2+, R tet 2+ ) leads to The energy difference E(+, R tet 2+ ) − E(2+, R tet 2+ ) defines the vertical electron affinity A(2+, R tet 2+ ) of the 2+ state (in its tetrahedral configuration), step 1 in Fig. 1, referenced to the top of the valence band, whereas gives the subsequent relaxation energy ∆(+, R D + , R tet 2+ ) in the positive charge state (step 2).
Similarly the neutral charge state emerges from the positive one by addition of an electron. Mathematically we again achieve this by adding and substracting first E(+, R D 0 ) and then E(+, R D + ) to and from the expression for the neutral formation energy E f Again A(+, R D 0 ) denotes a vertical electron affinity ) the relaxation energy from the neutral to the positive geometry in the positive charge state (step 3). An expression for the negative charge state can be obtained completely analogously once E f D (0, F = 0) has been computed.
The decomposition in Eq. (2) and (3) is not only appealing from an intuitive point of view, but also groups the required total-energy differences into two categories: lattice contributions in a fixed charge state and electron addition energies at fixed geometry. This permits us to go beyond a pure DFT description in an easy fashion by employing the most suitable method for each type of contribution [11]. Since we expect relaxation energies in the same charge state to be given reliably by LDA we retain DFT for the lattice part. For the electron addition energies, i.e., changes in charge state, which are typically problematic in LDA, we instead resort to the G 0 W 0 method.
The last remaining quantity to be assigned is E f tet (2+, F = 0), which we compute in the LDA. Unlike for the neutral state, the absence of DMC reference data unfortunately does not permit an assessment of the LDA error in this case. However, since the conductionband-derived defect levels are unoccupied the effect of the self-interaction and the band-gap error on the formation energy should be small. We therefore expect LDA to be more reliable for the tetrahedral 2+ state than for the neutral or the positive states.
The LDA calculations in the present work have been performed with the plane-wave, pseudopotential code S/PHI/nX [15]. 64-atom supercells were used throughout, unless otherwise noted. To remove the contributions arising from the homogeneous compensation charge density that is added to charged supercell calculations we have performed calculations for supercells with 64, 216 and 512 atoms. In these the interstitial atom was placed in the tetrahedral (2+) and the C 3v (+) position of a perfect (undistored) Si lattice. Fitting the formation energies up to cubic order in the inverse cell length and   extrapolating to infinite length we obtain corrections to the 64-atom cell of 0.17 eV and 0.04 eV for the 2+ and + state, respectively. Our extrapolated formation energy for the unrelaxed tetrahedral 2+ state of 3.19 eV agrees well the 3.31 eV obtained by Wright and Modine for a slightly larger lattice constant [16]. With this correction the formation energy of the relaxed tetrahedral 2+ configuration amounts to 2.65 eV. Applying a recently developed improved correction scheme [17] yields a corrected value of 2.66 eV, in excellent agreement with our extrapolated value.
For the G 0 W 0 calculations [18] we have employed the G 0 W 0 space-time code gwst [19,20,21]. For computational convenience we calculate the electron affinity of positive charge states (A(+, R D 0 )) by their inverse process, the electron removal from the neutral state, since no spin polarization or partially filled defect states are encountered then. Separate G 0 W 0 calculations for bulk silicon yield a band gap of 1.27 eV in good agreement with previous pseudopotential G 0 W 0 calculations [12].
The computed vertical electron affinities are shown in Figure 2. For comparison the LDA affinities calculated as Slater transition states [22] at half occupation have been included. The G 0 W 0 corrections for the +/0 state are similar for the three configurations and relatively small (∼0.2 eV). For states that in the LDA are in resonance with the conduction band, however, the G 0 W 0 corrections are much more pronounced. Since these states have a contribution from delocalized conduction-band states the delocalization error of the LDA [23] leads to a breakdown of Slater transition state theory. The resulting severe underestimation of the vertical affinities is akin to the band-gap problem. In LDA the band gap E g = I − A, where I is the ionization potential and A the electron affinity, is underestimated regardless of whether I and A are calculated as total energy differences or by Kohn-Sham eigenvalues, because the exchange-correlation functional is a continuous function of the electron density and therefore does not exhibit a derivative discontinuity. Many-body perturbation theoy in the GW approach, on the other hand, does not suffer from this problem.
Having identified the relevant electron affinities we can now return to the formation energies in Eqs. (2) and (3). Table I shows that already upon adding the first electron to the 2+ state we observe a large correction (∼0.9 eV) for the formation energy of the positive state. This error subsequently carries over to the neutral charge state, and adding the second electron incurs a further increase. The G 0 W 0 -corrected formation energies are now on average 1.1 eV larger than in the LDA. Since the quasiparticle shift of the empty defect state in the split<110> configuration is smaller than the band-gap opening the state is moved into the band gap (A(0, R split 0 )=1.1 eV). As a result the negative charge state becomes stable in G 0 W 0 , which is not the case in LDA, and has a formation energy of 5.53 eV.
For the neutral charge state our G 0 W 0 corrected formation energies compare well with recent DMC calculations that find an average increase of ∼1.5 eV (with a statitistical error bar of ±0.09 eV) and DFT HSE hybrid functional calculations that significantly reduce the selfinteraction error and yield an average increase of ∼1.2 eV [10]. Earlier DMC calculations give a larger average increase of 1.7 eV compared to the LDA, but also a much larger statistical error bar (±0.48 eV) [9]. Assuming a migration barrier of ∼0.2 eV [9] our computed activation enthalpy (formation energy + migration barrier) of ∼4.7 eV for the neutral split<110> interstitial is also in very good agreement with the experimentally determined value of 4.95 eV [24].
Finally we address the stability of the different defect configurations when the Fermi energy is varied throughout the band gap (cf Fig. 3). For clarity this is shown only for the split<110> configuration (lower panels) and the configuration with lowest energy at a given Fermi level (upper panels). The situation for the hex and C 3v configurations is qualitatively similar to that of the split<110>. If each configuration is considered separately, the formation energy diagram looks startlingly different in LDA and G 0 W 0 . Since LDA underestimates the formation energies of the + and 0 state relative to the 2+ it does not exhibit the negative-U behavior } for all Fermi energies) observed in G 0 W 0 . In addition G 0 W 0 stabilizes the negative charge state for the split<110> interstitial. If, on the other hand, the configuration with the lowest energy is considered, LDA and G 0 W 0 superficially give a more similar picture: the tetrahedral 2+ state is stable for 60-70% of the respective band gaps. While LDA then gives preference to the neutral split<110> for larger Fermi levels, the G 0 W 0 corrections marginally stabilize the neutral hex configuration, in agreement with the earlier DMC calculations [9]. The actual energies and transition levels between LDA and G 0 W 0 , however, differ appreciably.
Every point at which two lines in Fig. 3 cross corresponds to a charge-state transition level ε q/q . Bracht et al. have recently determined these for the silicon self-interstitial in high temperature diffusion experiments [4]. They identified two levels, at ≈ 0.1-0.2 eV and at ≈ 0.4 eV above the valence-band maximum, that they ascribed to ε 0/+ and ε +/2+ , respectively. These would most closely correspond to the G 0 W 0 -corrected chargestate transition levels ε 0/+ =0.09 eV and ε +/2+ =0.58 eV for the hexagonal configuration or ε 0/+ =0.05 eV and ε +/2+ =0.50 eV for the split<110>, while those of the C 3v configuration are noticeably different (0.62 eV and 1.24 eV). Although lowest in formation energy and therefore highest in concentration, the C 3v 2+ configuration (which is identical to tet 2+) most likely plays a negligible role in the diffusion experiments, since its diffusion would