All-electron Exact Exchange Treatment of Semiconductors: Effect of Core-valence Interaction on Band-gap and d-band Position 5 0 S. Sharma,∗ J. K. Dewhurst, and C. Ambrosch-Draxl 0 2 Institut fu¨r Physik, Karl–Franzens–Universit¨at Graz, n Universit¨atsplatz 5, A–8010 Graz, Austria a J (Dated: February 2, 2008) 7 1 Abstract ] i c Exact exchange (EXX) Kohn-Sham calculations within an all-electron full-potential method are s - l performedonarangeofsemiconductorsandinsulators(Ge, GaAs,CdS,Si,ZnS,C,BN,Ne, Ar,Kr r t m andXe). Wefindthattheband-gapsarenotasclosetoexperimentasthoseobtainedfromprevious . t a pseudopotential EXX calculations. Full-potential band-gaps are also not significantly better for m - sp semiconductors than for insulators, as had been found for pseudopotentials. The locations of d n d-band states, determined using the full-potential EXX method, are in excellent agreement with o c experiment, irrespective of whether these states are core, semi-core or valence. We conclude that [ 1 the inclusion of the core-valence interaction is necessary for accurate determination of EXX Kohn- v 3 Sham band structures, indicating a possible deficiency in pseudopotential calculations. 5 3 1 PACS numbers: 71.15.Mb,71.10.-w,71.22.+i 0 5 0 / t a m - d n o c : v i X r a 1 Exact treatment of exchange within the Kohn-Sham (KS) formulation of density func- tional theory (DFT) has been one of the most interesting developments in the last few years [1, 2, 3, 4, 5]. The KS band-gaps (EKS) for sp semiconductors, determined using exact g exchange (EXX) and the pseudopotential (PP) method, are found to be very close to ex- perimental gaps [3, 6], even in the absence of any correlation. This is surprising because the relation between EKS and the fundamental band-gapE of solids is given by E = EKS+∆ g g g g xc where, ∆ is the discontinuity in the exchange-correlation potential. This discontinuity is xc lacking in the local density approximation (LDA), resulting in the equality of the KS and fundamental band-gap [5] which leads to the famous band-gap problem. But for the EXX potential even in the absence of exact correlation, the discontinuity does not vanish and hence one cannot expect that the EXX KS gap alone can be identified directly with ex- perimental values. In fact, St¨adele et al. [3] determined the discontinuities in the exchange potentialexplicitly forSi, GeandGaAs, andfoundthemtobeabout3-5timestheband-gap. Furthermore, Magyar et al. [5] noted that the EXX KS gap for wide band-gap noble gas solids is not as close to the experimental gap (either optical or fundamental) as it is for sp semiconductors. One of the major conclusions of this work was that the EXX method does notprovide aKSbandstructure thatagreesequally wellwithexperiment forsemiconductors and insulators. All the above work has been done using the PP approach. This approach neglects both the exact core-valence exchange interaction as well as the relaxation of the core states in the solid-state environment. There do exist all-electron EXX calculations of the band-gaps in C, Si and Ge [4]. In this work, the linear combination of muffin-tin orbitals (LMTO) within the atomic sphere approximation (ASA) used indicates that the band-gaps are an improvement over LDA but do not agree with experiment as well as PP-EXX. These calculations however suffer from the shape approximation of the potential which ignores non-spherical contributions. In our recent article, the importance of this non-spherical part of the self-consistent EXX potential was demostrated[7]. At this point it is relevant to note that another method used to improve upon the single particle KS band-gaps is the GW approximation. Here a partially or fully self-consistent so- lution of the Dyson equation is used to determine the self-energy which provides a correction to the KS-LDA band-gap. It has been inferred from GW calculations performed using the PP approach that full self-consistency, even though it is essential for charge conservation, worsens the agreement of quasi-particle band-gaps with experimental data [8, 9, 10]. Very 2 recently it was noted by Ku and Eguiluz [11] that the success of the partially self-consistent GW approach is due to some kind of error cancellation which is an artifact of the PP and also that the deep-core states have an important effect on the band-gaps. Such an effect may be of similar importance for EXX calculations and requires thorough investigation. Another important contributory factor to the gap error is the cation d-band position in sp semiconductors, especially when these levels are relatively shallow (e.g. GaAs, ZnS, CdS) [12]. In these cases, use of the LDA leads to a misplacement of the d-bands, although this can be improved upon by the GW approximation or the self-interaction correction (SIC) to the potential [12, 13, 14]. Previous PP-EXX results indicate that even though the EXX potential is self-interaction free and the band-gaps are an improvement over those of LDA, the d-band positions are still not correct [15]. In this Letter, we present results of all-electron full-potential EXX calculations applied to representative semiconductors and insulators. This is to ascertain the effect of core-valence interaction and asymmetry in the potential on band-gaps and d-band positions, and also to find the reason for the inconsistency in EXX band-gaps of semiconductors and insulators. We have implemented EXX potential using the all-electron full-potential linearized augmented-plane wave (FP-LAPW) method [16] within the EXCITING code [17]. The start- ing point for these calculations is the exchange energy Ex[n] = −1 occ wkwk′ dr dr′φ∗ik(r)φ∗jk′(r′)φjk′(r)φik(r′), 2 |r−r′| ik,jk′ Z X where φ is the ith Kohn-Sham orbital of k-point k, w is the k-point weight, n(r) ≡ ik k occw φ∗ (r)φ (r) is the density, and i and j run over both core and valence states, in ik k ik ik cPontrast to PPs where the indices run only over the valence states. Proceeding in the same manner as G¨orling et al. [1, 2], the functional derivative chain rule is employed to obtain the exchange potential δE [n] x v [n](r) ≡ x δn(r) occ δE δφ (r′′) δE δφ∗ (r′′) δv (r′) = w dr′ dr′′ x ik + x ik s k δφ (r′′) δv (r′) δφ∗ (r′′) δv (r′) δn(r) ik Z (cid:20) ik s ik s (cid:21) X occ unocc φ∗ (r′)φ (r′) δv (r′) = dr′ w hφ |vˆNL|φ i jk ik +c.c. s , (1) k ik x jk ε −ε δn(r) " ik jk # Z ik j X X 3 where ε are the Kohn-Sham eigenvalues and ik hφik|vˆxNL|φjki = occ wk′ dr dr′φ∗ik(r)φ∗lk′(|rr′−)φrlk′′|(r)φjk(r′). (2) lk′ Z X In order to obtain the functional derivative δv (r′)/δn(r), the linear-response operator χ, s defined by the relation δn(r) = dr′ χ(r,r′)δv (r′), s Z must be inverted. This is not possible directly as χ has a zero eigenvalue corresponding to a constant eigenfunction. However, in a basis which excludes such functions the inversion of χ can be performed. Application of elementary perturbation theory to n leads to an explicit expression for the response occ unocc φ∗ (r)φ (r)φ∗ (r′)φ (r′) χ(r,r′) = w ik jk jk ik +c.c. k ε −ε ik jk ik j X X Implementations of EXX for pseudopotentials have used the plane wave basis for expanding χ. This is not feasible in a FP-LAPW approach because of the strongly varying wave functions inside the muffin-tin. Instead, we propose a new basis as follows. The overlap densitiesρ (r) ≡ φ∗ (r)φ (r)andtheircomplexconjugatesareusedasaspanningset,where α ik jk α ≡ (ik,jk) with i and j labelling occupied and unoccupied states, respectively. These are the states of some representative k-point which is chosen in advance. This spanning set is then reduced to a basis set as follows. First the overlap matrix of the spanning set elements O ≡ dr ρ∗(r)ρ (r) αβ α β Z is obtained and diagonalized. Note that O, which is positive semidefinite, has nonnegative eigenvalues. Next, all the eigenvectors of O corresponding to eigenvalues smaller than a certain tolerance, ǫ, are discarded. Finally, denoting the remaining eigenvectors as vβ, γ where β labels the vector and γ its coordinate, we seek a transformation matrix C such that if ρ˜ (r) = Cα vβρ (r) α β γ γ β γ X X then dr ρ˜∗(r)ρ˜ (r) = δ . α β αβ Z Straight-forward algebra shows that C should satisfy CC† = v†Ov −1 (3) (cid:0) (cid:1) 4 where v is the matrix of eigenvectors vβ in column-wise form. The matrix C is obtained by γ performing a Cholesky decomposition on the right hand side of Eq. 3. By virtue of their construction and the fact that dr ρ (r) = 0, the set of functions {ρ˜ } form an optimal α α basis for the expansion and inveRrsion of χ, as well as for the term in square brackets in Eq. 1. We should also point out that this basis may be useful for FP-LAPW time-dependent DFT response and GW methods. Specialattention should now bedrawn to the calculation ofthe non-localmatrix elements (Eq. 2)whichmaybedeterminedbyawell-establishedmethodforsolvingPoisson’sequation inaFP-LAPWenvironment [16]. Thedifferencesherearethatthedensitiesarenowcomplex and there is a long-range term arising when q ≡ k−k′ is close to zero. This is treated by considering theso-calledpseudo-chargedensity[16], whichischosentobesufficiently smooth within the muffin-tins so that it may be expanded in terms of plane waves, and yet has the same multipole expansion as the real density. If we restrict k and k′ such that q is in the first Brillouin zone (BZ), then the long-range (LR) contribution to the matrix elements from the pseudo-charge is occ 4πΩ hφ |vˆNL|φ i = w ρ∗(q)ρ (q), ik x jk LR q q2 il lj lq X where ρ (q) and ρ (q) are the pseudo-charge densities in reciprocal space. This sum suffers il lj from poor convergence with respect to the number of q-points. We therefore approximate it by an integral over a sphere of volume equivalent to that of the BZ. The final expression is 6Ω5 1/3 occ hφ |vˆNL|φ i ≃ 2 w ρ∗(q)ρ (q). ik x jk LR π q il lj (cid:18) (cid:19) lq X All the calculations in the present work are performed on a mesh of 27 special k-points in the irreducible BZ. As the response and its inverse are very sensitive to the number of empty states, we have taken special care to use a sufficient number to converge the band- gaps to within 0.01 eV. We find that almost all the materials studied here have converged d-band positions and band-gaps with around 25 empty states. All the results presented are generated using 30 empty states. One final practical note is that the eigenvalue cut-off, ǫ, discussed above, is taken to be 10−5 times the largest eigenvalue of O. For example for Si this results in 73 elements in the full basis set. In Fig. 1 we present the difference between the experimental fundamental band-gap and the KS band-gap for various semiconductors and insulators. The EXX potential is self- 5 FP-EXX 2 PP-EXX 0 LMTO-ASA-EXX FP-EXX-ncv FP-EXX-sph 1 FP-LDA ) V e -2 ( S K g E 0 - g -4 E -1 -6 -2 Ge GaAs CdS Si ZnS C BN Ar Ne Kr Xe FIG. 1: Difference between theoretical KS gaps and experimental fundamental gaps in eV. Full- potential exact exchange (FP-EXX), FP-EXX without core-valence interaction (FP-EXX-ncv), FP-EXX with only spherical exchange (FP-EXX-sph) and LDA (FP-LDA) results are from this work. The PP-EXX results are from Refs. 5, 6, 15. The LMTO-ASA-EXX results are from Ref. 4 and the experimental data is taken from Refs. 18, 19, 20. Note that LDA gap difference for Ne is -10.25 eV, and is off the scale of the plot. interaction free and hence an enhancement of the KS band-gap over that of the LDA is expected. This is indeed the case here. We find that for semiconductors the FP-EXX band- gaps are overestimated with respect to experiments by up to 60%. The band-gaps in all the insulators except Xe are underestimated with respect to the experimental fundamental gap up to about 24%. For Kr and Xe the KS band-gap is overestimated with respect to the experimental optical gap by 3.8% and 24.6% respectively, while for other insulators it is underestimated up to 6% . 6 Unlike for PP-EXX, there is not the excellent agreement of the KS gaps of sp semicon- ductors with experiment. In contradistinction to the PP-EXX results we find that the EXX is not inconsistent in its treatment of semiconductors and insulators, in fact agreement with experiment of the KS gaps for insulators is not any worse than for the semiconductors. In order to determine the reason for this behaviour we performed ground state calculations for these materials by ignoring the core-valence interaction term in the EXX potential. We emphasise here that this term was removed only from the exchange potential and not from the rest of the potential. Even this has substantial effect on the KS band-gaps. The results are marked as FP-EXX-ncv in Fig. 1. Compared to FP-EXX, the FP-EXX-ncv band-gaps for all the semiconductors are closer to the PP-EXX and the experimental values. On the other hand, the agreement of the KS band-gaps with experiments for insulators now wors- ens. This trend is concomitant with the findings of the PP-EXX calculations [5, 6] where the performance of the EXX potential in determining the KS gaps is excellent for the semi- conductors and lacking for the insulators. Previous all electron EXX calculations performed within the ASA show a similar overestimation of the gaps for semiconductors [4]. All this suggests that the lack of core-valence interaction is responsible for the spurious agreement of the PP-EXX KS band-gaps with the experiments for sp semiconductors. In our work on magnetic metals [7] we found that inclusion of the non-spherical contribu- tions to the exchange potential is crucial to obtain the correct ground-state. The results for the KS band-gaps obtained using only the spherical contribution to the EXX potential are marked as FP-EXX-sph in Fig. 1. The effect of the shape approximation to the potential in all the materials other than noble gas solids is substantial. Which is expected since the noble gas solids are composed of almost independent atoms. Another property which is known to lead to errors in the band-gaps[12]and for which the LDA performs badly is the position of the semi-core d-states. The LDA under-binds leading to displacement of d-bands much above the experimental value. In the past, a quasi-particle bandstructure was determined using the self energy calculated with the GW approximation to correct the d-band position [13]. Also the SIC to the potential was found to lead to good agreement between the theoretical and experimental position of these d-bands. Since the EXX is self-interaction free in nature, it should improve the d-band positions. Counter intutively Rinke et al. [15] found that even though the PP-EXX calculations result in a downward motion of the d-bands with respect to LDA, still GW calculations are needed to 7 Ge GaAs InP ZnS CdS ) V 6 e ( t n e 4 m ri e 2 p x e o t 0 e v ti a -2 el FP-EXX r PP-EXX n o -4 GW-FPLMTO ti GW-LMTO-ASA si o GW-PP-EXX p -6 HF-FPLMTO d n PP-SIC a FP-EXX-ncv b - -8 FP-LDA d FIG. 2: Difference between theoretical d-band eigenvalues and experimental data in eV. Full- potential exact exchange (FP-EXX), FP-EXXwithoutcore-valence interaction (FP-EXX-ncv) and LDA (FP-LDA) results are from this work. The results obtained using the GW calculation on top of the PP-EXX ground-state (GW-PP-EXX) are from Ref. 15, the self-interaction corrected potential within PP approach (PP-EXX) are from Ref. 14. GW results calculated using a FP- LMTO method (GW-FPLMTO) and Hartree-Fock FP-LMTO (HF-FPLMTO) are from Ref. 12. The GW calculation on top of the LMTO-ASA-LDA ground-state (GW-LMTO-ASA) are from Ref. 13. fully correct the positions. The FP-EXX potential, in which the core and valence states are treatedonthesame footing,istrulyself-interaction freeinnatureandit wouldbeinteresting to see how this effects these semi-core/core-like d-bands. In the present work, we compare FP-EXX d-band eigenvalues with those of previous PP-EXX, experiment and various GW and SIC calculations in Fig. 2. The compounds are 8 chosen so that d-bands are in various energy regimes to get an overall picture: core-like (Ge: experimentally -29.6 eV from the top of the valence band), semi-core (GaAs: -18.8 eV, InP: -16.8 eV) and valence (ZnS: -8.7 eV, CdS: -9.2 eV). As can be seen, the eigenvalues of d-states obtained using the FP-EXX method are in excellent agreement with experiments. In fact, in several cases it is better than that of the FP-GW results. Here again the agreement between the FP-EXX-ncv and PP-EXX values of Rinke et al. [15] is quite good. This indicates that the inclusion of core-valence interaction is crucial for determination of the correct d-band eigenvalues. The main aim of our paper is to establish the effects of the inclusion of core-valence interaction in full-potential EXX calculations. Prior to this investigation, it appeared that EXX might be the answer to the band-gap problem at least for sp semiconductors. We find that this is not the case if core-valence interaction is included in the calculation of the EXX potential. We also find that EXX is consistent in its treatment of semiconductors and insulators. Ourfindings indicateapossible shortcoming inthePP-EXXmethodarising from thelackofexplicitcore-valenceinteractionsinthePP-EXXpotential,asmentionedearlierby Ku and Eguiluz [11]. Furthermore, we show that the d-band eigenvalues determined using the FP-EXX method are in excellent agreement with experiment irrespective of whether these states are in the core, semi-core or valence. In most cases the agreement is better than that achieved by the full-potential GW calculations. WeacknowledgetheAustrianScienceFund(projectP16227)andtheEXCITINGnetwork funded by the EU (Contract HPRN-CT-2002-00317) for financial support. ∗ Electronic address: [email protected] [1] A. Go¨rling and M. Levy, Phys. Rev. A 50, 196 (1994). [2] A. Go¨rling, Phys. Rev. B 53, 7024 (1996). [3] M. Sta¨dele, J. A. Majewski, P. Vogl, and A. Go¨rling, Phys. Rev. Lett. 79, 2089 (1997). [4] T. Kotani, Phys. Rev. Lett. 74, 2989 (1995). [5] R. J. Magyar, A. Fleszar, and E. K. U. Gross, Phys. Rev. B 69, 045111 (2004). [6] M. Sta¨dele, M. Moukara, J. A. Majewski, P. Vogl, and A. Go¨rling, Phys. Rev. B 59, 10031 (1999). 9 [7] S. Sharma, J. K. Dewhurst, C. Ambrosch-Draxl, and C. Persson (2005), URL http://arxiv.org/abs/cond-mat/0501258. [8] B. Holm and U. von Barth, Phys. Rev. B 57, 2108 (1998). [9] W.-D. Scho¨ne and A. G. Eguiluz, Phys. Rev. Lett. 81, 1662 (1998). [10] B. Holm, Phys. Rev. Lett. 83, 788 (1999). [11] W. Ku and A. G. Eguiluz, Phys. Rev. Lett. 89, 126401 (2002). [12] T. Kotani and M. V. Schilfgaarde, Solid State Comm. 121, 461 (2002). [13] F. Aryasetiawan and O. Gunnarsson, Phys. Rev. B 54, 17564 (1996). [14] D. Vogel, P. Kru¨ger, and J. Pollmann, Phys. Rev. B 54, 5495 (1996). [15] P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, and M. Scheffler, unpublished (2004). [16] D. J. Singh, Planewaves pseudopotentials and the lapw method (1994). [17] J. K. Dewhurst, S. Sharma, and C. Ambrosch-Draxl, EXCITING code developed under the Research and Training Network EXCITING funded by the EU, contract No. HPRN-CT-2002- 00317 (2004), URL http://www.exciting.physics.at. [18] K.-H.Hellwege, O.Madelung,M.schulz,andH.Weiss,eds.,semiconductors. Physicsof Group IV Elements and III-IV Compounds., vol. 17 (Springer Berlin, 1982). [19] A. Mainwood, Properties and growth of diamond (1994). [20] D. R. T. Zahn, G. Kudlek, U. Rossow, A. Hoffman, I. Broser, and W. Richter, Adv. Mater. Opt. Electron. 3, 11 (1994). 10