A finite electric-field approach to evaluate the vertex correction for the screened Coulomb interaction in the quasiparticle self-consistent GW method
aa r X i v : . [ c ond - m a t . s t r- e l ] M a r A finite electric-field approach to evaluate the vertex correction for the screenedCoulomb interaction in the quasiparticle self-consistent GW method
Hirofumi Sakakibara and Takao Kotani
Department of Applied Mathematics and Physics, Tottori university, Tottori 680-8552, Japan
Masao Obata and Tatsuki Oda
Department of Computational Science, Institute of Science and Engineering,Kanazawa University, Kakuma, Kanazawa 920-1192, Japan (Dated: March 6, 2020)We apply the quasiparticle self-consistent GW method (QSGW) to slab models of ionic mate-rials, LiF, KF, NaCl, MgO, and CaO, under electric field. Then we obtain the optical dielectricconstants ǫ ∞ (Slab) from the differences of the slopes of the electrostatic potential in the bulk andvacuum regions. Calculated ǫ ∞ (Slab) show very good agreements with experiments. For example,we have ǫ ∞ (Slab)=2.91 for MgO, in agreement with the experimental value ǫ ∞ (Experiment)=2.96.This is in contrast to ǫ ∞ (RPA)=2.37, which is calculated in the random-phase approximation forthe bulk MgO in QSGW. After we explain the difference between the quasiparticle-based pertur-bation theory and the Green’s function based perturbation theory, we interpret the large difference ǫ ∞ (Slab) − ǫ ∞ (RPA) = 2 . − .
37 as the contribution from the vertex correction of the properpolarization which determines the screened Coulomb interaction W . Our result encourages the the-oretical development of self-consistent G W approximation along the line of QSGW self-consistency,as was performed by Shishkin, Marsman and Kresse [Phys. Rev. Lett. , 246403(2007)]. PACS numbers: 71.10.-w, 71.15.-m, 71.15.Dx
I. INTRODUCTION
The quasiparticle self-consistent GW (QSGW) is oneof the most reliable method to determine the one-particleeffective Hamiltonian which describes the independent-particle picture, or the quasiparticle (QP) picture, fortreating electric excitations of materials [1–3]. Othercompetitive methods such as HSE [4] and Tran-Blaha-09 functional [5] may work well in many systems, al-though we may need to use material-dependent parame-ters [6]. In contrast, QSGW is virtually parameter freeand gives reliable descriptions for a wide range of ma-terials, not only metals and semiconductors, but alsotransition-metal oxides, type-II superlattice, and 4 f sys-tems [7–10]. Since heterogeneous mixtures of materialsare used in current technologies, QSGW is worth to bedeveloped more as a tool to treat electronic structures ofsuch materials, where methods including such material-dependent parameters are hardly applicable.However, QSGW as it is has a shortcoming that it givesa systematic overestimation of the exchange effects. Thisresults in a little larger band gaps in QSGW for mate-rials. In fact, Faleev, van Schilfgaarge, and Kotani [1–3]made a suggestion that the overestimation is removed ifwe perform improved QSGW calculations taking into ac-count the enhancement of the screening effect due to theelectron-hole correlation in the evaluation of the screenedCoulomb interaction W . This is based on the theoreti-cal consideration combined with the observation that thecalculated optical dielectric constant ǫ ∞ in the random-phase approximation (RPA) in QSGW gives ∼
20 per-cent smaller ǫ ∞ than experiments for kinds of materials[3, 11]. Such an improved calculation was performed byShishkin, Marsman and Kresse, where they include theenhancement of the screening effect [12]. The enhance-ment is via the vertex correction for the proper polar-ization P , which determines W = v/ (1 − vP ), where v denoted the Coulomb interaction. They approximatelyinclude the lowest-order vertex correction due to theelectron-hole correlation; see Eq.(15) around in Ref.13.Their results are theoretically quite satisfactory in thesense that both band gaps and ǫ ∞ , which are calcu-lated simultaneously and self-consistently without pa-rameters as in HSE, are in agreement with experiments.For example, calculated values ǫ ∞ = 2 .
96 and band gap E G = 8 . GW (e-h) inTable I in Ref.12. Furthermore, based on these theoreti-cal analyses, we can introduce QSGW80 to avoid the veryexpensive computational costs of the method by Shishkin et al. . QSGW80 is just a simple hybridization, 80 %QSGW+ 20 % GGA to include such enhancement of thescreening effectively. The hybridization is very differentfrom the hybridization in HSE, where the mixing ratio α between GGA and the Hartree-Fock strongly affectsto final results. QSGW80 works well to describe experi-mental band gaps [14]. The performance of QSGW80 issystematically examined in Ref.8 by Deguchi et al. , wherewe see both the calculated band gaps and effective massesare in good agreements with experiments. QSGW80 issuccessfully used for practical applications, for example,to the type-II superlattice of InAs/GaSb [15, 16].In this paper, we evaluate ǫ ∞ not in bulk calculationswith such approximations used in Ref.12, but by slabmodels with finite electric bias voltage. We treat fiveionic materials, LiF, KF, NaCl, MgO and CaO. We puta slab in the middle of vacuum region in a supercell. Theelectric field is applied by the effective screening medium(ESM) method given by Otani and Sugino [17]. We ob-tain ǫ ∞ from the ratio of slopes of the electrostatic fieldsat slab region and at vacuum region. Our approach isbased on the self-consistent method, thus we do not needto utilize approximations as was used in Ref.12. Sincewe explicitly treat the response to the bias, our methodincludes higher-order effects in a self-consistent manner.Our findings are that the calculated ǫ ∞ in QSGW forthe slab models are very close to experimental values.This is in contrast to the fact that ǫ ∞ in RPA of QSGWare generally ∼
20 percent smaller than experimental val-ues. This indicates that the vertex correction at the levelof derivative of the QSGW self-energy should make W bein agreements with experiments. Our results is consistentwith Table II in Ref.12.We can interpret the enhancement of screening, rep-resented by the enlargement of ǫ ∞ , as the size of thevertex correction for the proper polarization P . Notethat the vertex correction we evaluate is not what is de-fined in the Hedin’s equation [18]. In the equation, we see P = − iGG Γ, that is, the vertex function Γ is for the cor-rection to P = − iGG , where G denotes the Green’s func-tion. Instead, we rather evaluate Γ for P = − iG G Γ,where G is the bare Green’s function. To clarify theabove theoretical point on Γ, we give an extensive discus-sion in Sec.II. We explain role of Γ in the two kinds of per-turbation theories. In Sec.III, we explain QSGW+ESM,an implementation of QSGW combined with ESM. TheQSGW+ESM for slab models should be very useful notonly for our purpose here, but also for others where usualGGA+ESM have difficulties. In Sec.IV, we show our re-sults of ǫ ∞ . Then they are interpreted as the vertex cor-rection. In Sec. IV B, we give a rationale of QSGW80,followed by a summary. II. QP-BASED PERTURBATION VS. G -BASEDPERTURBATION To make our motivation in this paper clarified, we haveto clarify the difference between the quasiparticle-basedperturbation (QbP) and the Green’s function-based per-turbation (GbP). QbP is based on the Landau-Silin’s QPtheory, while GbP is on the Hedin’s one. To illustrate thedifference between QbP and GbP, we give a narrow-bandmodel as follows. The model represent situations wherewe have good QP picture (= independent-particle pic-ture).In advance, remind that we mainly have two kinds ofexcitations in the paramagnetic electronic systems. Thatis, the multi-particle excitations, and the collective exci-tation such as plasmons. The former is described by QPsinteracting each other. Note that plasmons is located athigh energy because of the long-range Coulomb inter-action [19]. These excitations can be hybridized. For example, we know pseudo plasmon in Silver, where one-particle excitations of 3 d electrons are hybridized withplasmons. A. narrow-band model to explain the QP-basedperturbation
The QPs based on the Landau-Silin’s Fermi liquid the-ory is originally for metals [19]. However, the idea ofQPs are rather easily applicable to insulators. We canconsider QbP based on the QPs. To illustrate this, letus consider a narrow-band model, a paramagnetic casegiven by a Hamiltonian H , which has an one-body termrepresented by finite numbers of Wannier functions in theprimitive cell, and the Coulomb-like interaction e ǫ ′ | r − r ′ | ,where ǫ ′ is a constant. We consider a case that it givesthe QPs shown in Fig.1 given by the one-particle Hamil-tonian H . FIG. 1. Band structure of a narrow-band model to illustratethe quasiparticle based perturbation. CBM and VBM arethe acronyms of conduction-band minimum and valence-bandmaximum. The width of bands B C and B V are smaller thanthe band gap E G . As in the text, we expect well-defined QPsin this model. Based on the perturbation, we expect that the lifetimesof all the electrons are infinite because the band gap E G is large enough to forbid all the electrons decaying intolower energy electrons accompanying electron-hole pairs,holes as well. In other words, all bands are within thethreshold of impact ionization [20]. Thus QPs describedby G = 1 / ( ω − H ) should give well-defined one-particleexcitations of the narrow-band model.This conclusion should be essentially kept even whenwe fully turn on the interaction as long as the follow-ing conditions are well satisfied. First, excitons, bindingstates of electron-hole pairs, should be only slightly atlower energies than E G . Second, plasmons should be lo-cated at high enough energies so that the plasmon arehardly hybridized with the one-particle excitations. Un-der these assumptions, we have well-defined QPs. Wecan consider a path of adiabatic connection with keepingthe well-defined QPs given by H since the QP spectrumis clearly separated from the other excitations. It can bewritten as H λ = H + ( H λ − H ) for λ is from zero tounity, where H = H λ =1 . Note that ( H λ − H ) need tocontain λ -dependent one-body term.In this narrow-band model, QbP should be suitable;QPs are interacting each other by ( H λ =1 − H ). Asfor the proper polarization function, we can take allthe non-interacting two-body QP excitations correctlyby P = − iG G . Based of this QbP, G W approxi-mation, Σ = iG W where W = v/ (1 − vP ), is phys-ically justified. That is, it describes how the motion ofQPs in G is perturbed by the dynamical self-interactiongiven by W in RPA. Most of all first-principles calcu-lations in literatures are implicitly based on this QbP,while some latest literatures [21, 22] are based on GbPexplained in Sec.II C. QSGW is a method to determine H self-consistently in QbP. B. vertex function Γ QbP in the QP-basedperturbation
To improve G W , we may include electron-hole cor-relation. The corrections replace P = − iG G with P where we include the correlation via the Bethe-Salpeterequation (ladder diagrams). That is, we include two-body spectrum in the proper polarization accurately.This lets us to use W = 1 / (1 − vP ) instead of W , re-sulting the G W approximation as Σ = iG W . In thispaper, we concentrate on the G W approximation as inthe case of sc GW (e-h) in Ref.12. In QbP, we thus definethe vertex function Γ QbP for W as P = − iG G Γ QbP .Roughly speaking, P = P Γ QbP . We give the numeri-cal evaluation for this Γ
QbP via the evaluation of ǫ ∞ asshown in Sec.IV A.To go beyond Σ = iG W approximation here, we needto take into account three-particle intermediate states.However, it is not theoretically straightforward becauseof a double counting problem that Σ = iG W alreadypartially takes into account such states. We will need toconstruct theories of three-particle problem without dou-ble counting along the line of first-principles calculations.We do not treat this problem in this paper. C. G -based perturbation and vertex function Γ GbP
Let us consider how we can apply GbP to the narrow-band model. In contrast to G in QbP, the one-bodyGreen function G has complex meanings. We have imag-inary part of G at high energies, representing QPs hy-bridized with plasmons (plasmarons). Because of sumrule for Im[ G ], the QP parts are suppressed by Z -factor.Thus there is a problem that P = − iGG do not contain the two-body non-interacting excitations with the correctweight, in contrast to the case P = − iG G .In principle, this problem is corrected by including thevertex function Γ GbP in the Hedin’s equation to deter-mine the one-particle Green’s function G (1 ,
2) [18]. Be-cause Hedin’s equation is theoretically rigorous, we ex-pect − iGG Γ GbP ≈ P = − iG G in the model, underthe discussion of Sec.II A that P gives good approxi-mation for the model. That is, contributions related tothe collective excitations and renormalization factors Z in − iGG should be virtually taken away by the factorΓ GbP . However, such numerical calculations should becomputationally very demanding [23]. Similar discussionof Z -factor cancellation is also seen when we multiply G to W . That is, we should have G W ≈ GW Γ GbP since QbP correctly treats the model. Our analysis hereis consistent with Takada’s analysis based on the Wardidentity [24].QbP should be generally superior to GbP even in realmaterials. In contrast to GbP, QbP is quite simple andphysically convincing. We should not be confused withthe similarity of QbP and GbP. In this paper, we evaluateΓ
QbP . In the following, we calculate the enhancement ofthe screening effect. Then we evaluate the size of theratio
P/P from the comparison between calculated ǫ ∞ in RPA and ǫ ∞ in the slab models. This ratio gives thesize of Γ QbP . III. QSGW COMBINED WITH EFFECTIVESCREENING MEDIUM METHOD
To calculate ǫ ∞ , we put a slab in the middle of vacuumregion in a supercell. ǫ ∞ is calculated from the differenceof slopes in the vacuum region and in the slab region un-der small bias voltage. The supercell we use are detailedat the beginning of Sec.IV. In such calculations, we canobtain ǫ ∞ beyond the bulk calculation in RPA as we ex-plain in the next paragraph. That is, we can obtain ǫ ∞ including the effect of the vertex correction.To explain how the effect is included, let us considerslab calculations in the case of GGA at first. We first per-form self-consistent calculation under zero bias. Then weperform self-consistent calculation under the finite bias(theoretically, it should be infinitesimally small). Thenwe have the difference of the electron density δn ( r ) be-tween the two calculations. Simultaneously, we havea corresponding response of the one-particle potentialgiven as δV ( r ) = R d r ′ v ( r − r ′ ) δn ( r ′ ) + ∂V GGAxc ∂n ( r ) δn ( r ).The last term is the difference in the exchange-correlation(xc) potential caused by δn ( r ) self-consistently. That is,the derivative δV ( r ) δn ( r ′ ) contains the contribution of the xckernel f xc = ∂V GGAxc ∂n ( r ) . Under the bias, we can obtain ǫ ∞ ,from the ratio of slopes of the electrostatic potential inthe vacuum region and in the slab region. It should con-tain the contribution from f xc , which is identified as thevertex correction in GGA.This is essentially the same in QSGW. Recall thatthe self-energy in QSGW denoted as V QSGWxc ( r , r ′ ) is astatic non-local potential, replacing V GGAxc . The deriva-tive of the one-particle potential is given as δV ( r , r ′ ) = R d r ′′ v ( r − r ′′ ) δn ( r ′′ ) + δV QSGWxc ( r , r ′ ), where the lastterm play a role of ∂V GGAxc ∂n ( r ) δn ( r ). Note that δV QSGWxc ( r , r ′ )is determined self-consistently, although it is not so sim-ply given as ∂V QSGWxc ∂n ( r ) δn ( r ). Our calculations include thecontribution of δV QSGWxc ( r , r ′ ) self-consistently as in thecase of GGA. Our method is on the same spirit of solvingthe Bethe-Salpeter equation in Ref.25.Our QSGW+ESM is implemented in a first-principlespackage ecalj [8, 26] which is based on a mixed-basismethod, the augmented plane wave (APW) and Muffin-tin (MT) orbital method (the PMT method) [27–30].The PMT method is an all-electron full potential methodwhich uses not only the APW basis in the LAPWmethod, but also the MT orbitals in the LMTO methodsimultaneously in the expansion of eigenfunctions. Italso use the local orbital basis [31]. On top of the PMTmethod, we had implemented the QSGW method [8, 29].In PMT, we use very localized untuned MTOs which con-tains damping factor exp( − κr ), where κ are fixed to be1/bohr and/or 2/bohr, together with low-cutoff APWs( ≤ GW methods which requires theWannier-interpolation technique to make band plots inthe whole Brillouin zone, we can make band plots easilywithout resorting to the technique [8]. In the following,we show how to implement ESM in the PMT method,after an explanation of general theory of ESM. A. The electrostatic potential in the Effectivescreening medium method
We apply the ESM method [17] to slab models un-der an external electric field. We treat a supercell withperiodic boundary condition where we have a slab withperiodicity in the xy -plane. The slab is at the middle ofsupercell. Position in the cell is specified by r = ( r // , z ).Planes at z = − z and at z = z are the left and rightends of the supercells. The electrostatic potential is cal-culated from the charge density in the supercell assum-ing two electrodes are at z = ± z (we set z = z inFig.1 of Ref.17.) for applying voltage to the supercell.As we summarize as follows, the ESM in DFT is formu-lated from the total energy minimization, however, it isnot true in QSGW since QSGW itself is not formulatedfrom the total energy minimization. After we obtain thefollowing key equation Eq. (2) to determine electrostaticpotential. We use it even in QSGW.Let us start from the energy functional of DFT in the ESM. It is written as E [ n ] = E kin [ n ] + E xc [ n ] + E es [ n ] + E app [ n ] . (1)Here, we have kinetic energy E kin [ n ], xc energy E xc [ n ],and the electrostatic energy E es [ n ] terms. In addition,the last term is the applied electrostatic term E app [ n ] = R d rV app ( r )( n ( r )+ n N ( r )), where n ( r ) and n N ( r ) are theelectron density and the the charge density of nuclei, re-spectively; V app ( r ) is a linear function of z , representingthe external field.In ESM, we enforce the periodicity in the supercell forthe electrostatic potential. Thus we use V app ( r ) s ( r ) in-stead of V app ( r ), where we introduce a support function s ( r ) which is unity for most of all regions, but is goingto be zero at z = − z and z = z . It is different fromunity only near the boundaries, z ≈ − z or z ≈ z . Thusthe potential V app ( r ) s ( r ) recover the periodicity of thesupercell. A constant can be added to V app ( r ) so that itkeep smooth periodicity over z = ± z . As long as we uselarge enough vacuum region, we have little electrons nearthe boundaries. Thus the choice of s ( r ) is irrelevant.A key in ESM is that we use the Green function ¯ v ( r , r ′ )for the electrostatic energy E es [ n ] instead of the Coulombinteraction v ( r − r ′ ) in usual the GGA calculations. Asin Ref.17, ¯ v ( r , r ′ ) contains not only the Coulomb inter-action v ( r − r ′ ) but also the effects due to the polariza-tion of virtual electrodes, which are at z = ± z (we use z = z in our calculations here). Polarization of the slaboccurs with keeping the electrostatic potential being con-stants at electrodes. Corresponding to V app ( r ) s ( r ), weuse s ( r )¯ v ( r , r ′ ) s ( r ′ ) instead of ¯ v ( r , r ′ ) in practice. Thenwe have well-defined Kohn-Sham total energy with keep-ing the periodic boundary condition for given V app ( r ).The minimization of E [ n ] with respect to n ( r ) givesthe Kohn-Sham potential V ( r ) as V ( r ) = Z d r ′ ¯ v ( r , r ′ )( n ( r ′ ) + n N ( r ′ )) + V app ( r ) + V xc ( r ) , (2)Hereafter, we skip s ( r ) for simplicity.In QSGW [29], we cannot derive its fundamentalequation from the energy minimization. Thus theformulation of QSGW+ESM is not exactly along theline above. However, we can use the one-particlepotential of Eq. (2) in the self-consistent cycle, where V xc ( r ) is replaced by a static version of the self-energy[1]. Thus, in principle, it is straightforward to performQSGW+ESM. B. ESM in the PMT method
In the PMT method, electron density (and also thecharge density) is represented by the three componentformalism described in Ref.30, originally introduced bySoler and Williams [33–35]. At first, space is dividedinto MT regions and interstitial regions. Then elec-tron density is represented by three components as n = { n ( r ) , { n ,a ( r ) } , { n ,a ( r ) }} where a is the index ofatomic sites in the primitive cell. Following Ref.30, thisis simply expressed as n = n ⊕ n ⊖ n . The 0thcomponent n ( r ) is the spatially smooth functions, ex-panded in analytic functions, that is, planewaves, Gaus-sians, and smooth Hankel functions [29]. The 1st com-ponents n ,a ( r ) is the true electron density within MTat R a . The 2nd components n ,a ( r ) is the counter part,that is, the projection of n ( r ) into the MT at R a . n ( r )and n ,a ( r ) are identical within MT at R a up to givenangular momentum cutoff in their spherical harmonicsexpansion.We can get all charge density n Zcv = n Zcv0 ⊕ n Zcv1 ⊖ n Zcv2 by adding the ion-core contribution to n .Then we apply the multipole transformation clearly de-fined in Ref. 30, resulting ¯ n Zcv0 ⊕ ¯ n Zcv1 ⊖ ¯ n Zcv2 as shown inEq.(28-30) in Ref.30. The transformation makes ¯ n Zcv0 ( r ),¯ n Zcv1 ,a ( r ), and ¯ n Zcv2 ,a ( r ) have the same multipole in eachMT site at R a , although physically observable densityunchanged. The 1st components ¯ n Zcv1 ,a ( r ), unchanged bythe transformation, are the sum of ion-core charge den-sity and n ,a ( r ).From the smooth density ¯ n Zcv0 ( r ), we cancalculate electrostatic potential as V es0 ( r ) = R d r ′ ¯ v ( r , r ′ )¯ n Zcv0 ( r ′ ) + V app ( r ). This gives correctinterstitial part of the potential V es0 ( r ) calculated fromall charge density. The values of V es0 ( r ) at MT bound-aries are used to determine the electrostatic potentialwithin MTs.We can use usual procedure to determine the electro-static potential within MTs. In each MT, we have 1stand 2nd components ¯ n Zcv1 ,a ( r ) and ¯ n Zcv2 ,a ( r ), which havethe same multipole. With the condition that the elec-trostatic potential is zero at the MT boundary, we cancalculate the potential generated by the difference of the1st and 2nd components.Thus we finally have the electrostatic potential V es ( r )represented in the three component formalism. With thispotential, we can perform self-consistent calculations forslab models. IV. RESULTSA. Optical dielectric constants via the slab model
In Fig.2, we illustrate our treatments in the slab mod-els for five NaCl-structure ionic materials, where we use ± z = ±
30 a.u. We use slabs made of nine layers, 18atoms in the supercell. We use experimental lattice con-stants of bulk materials, without relaxation of atomicpositions. The electrostatic potential V es z ( z, E ) are theaverage of V es0 ( r ) in the xy plane under the bias voltage E . We plot the cases of E = 0 . E = 0 . V es z ( z ) = V es z ( z, . − V es z ( z, . V es z ( z, E ) changes linearly as afunction of z in the vacuum regions and in the slab re-gions. From the ratio of two slopes of ∆ V es z in the slab FIG. 2. A slab (18 atoms per cell) is placed in in the middleof a supercell (60 a.u. width along the z axis which is per-pendicular to the slab), with electrodes at the left and rightends. At the top panel, we show V es z ( z, E ) for E = 0 . E = 0 . V es z ( z ). From the ratio of two slopes of ∆ V es z in theslab region(green) and in the vacuum region(violet), we ob-tain ǫ ∞ (Slab). We have better numerical accuracy by using∆ V es z ( z ) instead of V es z ( z, . region and in the vacuum region, we obtain ǫ ∞ (Slab).Our main results are ǫ ∞ calculated from slab mod-els in QSGW, ǫ ∞ (QSGW , Slab), in Table I. Note that ǫ ∞ (Slab) contains the effect of vertex corrections basedon QbP (See Sec.II), because changes of the self-energycaused by the bias E are self-consistently taken into ac-count. Numerical reliability of our calculations are es-timated to be . ǫ ∞ (RPA). To obtain them, we first performself-consistent calculations in QSGW for bulk materials.Then we calculate ǫ ∞ in the random-phase approxima-tion (RPA) with/without local field correction (LFC).We also show ǫ ∞ in GGA together.The QSGW values are in good agreements withexperiments. For example, ǫ ∞ (QSGW , Slab)=1.94for LiF gives surprisingly good agreement with ǫ ∞ (Experiment)=1.96. In contrast, ǫ ∞ (QSGW , RPA)=1.67 is very smaller than ǫ ∞ (QSGW , Slab)=1.94. These are generallytrue in all other materials. We see that ratios
TABLE I. Calculated optical dielectric constant ǫ ∞ . ‘RPA’are in bulk calculations with local field correction (LFC).‘RPA(noLFC)’ are without LFC. ‘Slab’ are calculated fromthe slab models in the setting of Fig.2. Ratios η = ǫ ∞ (RPA) ǫ ∞ (Slab) and γ = ǫ ∞ (Slab) − ǫ ∞ (RPA) − are calculated just simply from the valuesof ǫ (QSGW , RRA) and ǫ (QSGW , Slab).RPA RPA Slab η γ
Experiments(noLFC) [36–38]LiF GGA 2.04 1.95 2.01 1.96QSGW 1.73 1.67 1.94 0.86 1.40KF GGA 2.16 1.96 1.94 1.85QSGW 1.79 1.68 1.86 0.90 1.26NaCl GGA 2.70 2.33 2.42 2.34QSGW 2.13 1.92 2.31 0.83 1.42MgO GGA 3.17 2.96 3.09 2.96QSGW 2.50 2.37 2.91 0.81 1.39CaO GGA 3.94 3.59 3.68 3.33QSGW 2.88 2.68 3.31 0.81 1.38 η = ǫ ∞ (QSGW , RPA) /ǫ ∞ (QSGW , Slab) in Table Iare ∼ ǫ ∞ forZnO, Cu O, MnO, and NiO are presented. From apoint of view to estimate the enhancement factors ( ≈ vertex Γ ) of the proper polarization, we may considerratios γ = ǫ ∞ (Slab) − ǫ ∞ (RPA) − . As shown in Table I, γ ∼ ǫ ∞ (QSGW , Slab) gives very good agreements with ǫ ∞ (experiment), we can say that the vertex correctionfor bulk should give the difference between ǫ ∞ (RPA) and ǫ ∞ (experiment) very well, where the vertex correction iscalculated at the level of the functional derivative of theself-energy in QSGW; See Sec.III.This is in contrast to the case of GGA. For example,look into the case of LiF. The difference ǫ ∞ (GGA , Slab) − ǫ ∞ (GGA , RPA) = 2 . − .
95 = 0 .
06 is very small. Thedifference is originated from the xc kernel f xc in the den-sity functional perturbation theory. This is consistentwith results in Ref.40 where they explicitly evaluate f xc in GGA for bulk materials. Note that ǫ ∞ (GGA , Slab) =2 .
01 is a little larger than ǫ ∞ (experiment) = 1 .
96; thisis true for all other materials. We see that the contribu-tions of vertex corrections f xc do not necessarily improveagreements; ǫ ∞ (GGA , Slab) give poorer agreement with ǫ ∞ (experiment) than ǫ ∞ (GGA , RPA).
B. Rationale for QSGW80
Our result in Sec.IV A shows that the vertex correc-tion should be included in the proper polarization P toobtain ǫ ∞ in agreement with experiments. We have touse such P in the QSGW self-consistent cycle. Such im-proved QSGW self-consistency can be identified as a self-consistent method in the G W approximation on the ba-sis of QbP. Ref.12 by Shishkin et al. gives a method onthis idea.However, their method is too expensive for computa- TABLE II. Calculated band gaps (eV) of bulk materials. InQSGW80, we show self-consistent results with the hybrid xcpotential, 80 % QSGW+ 20 % GGA. QSGW80nosc specifiesone-shot calculations with the hybrid potentials after QSGW100% self-consistent calculations. QSGW80nosc is slightlylarger because it is not fully self-consistent under such xc po-tential. See Ref.8.Experiments QSGW QSGW80 QSGW80nosc GGA[41–44]LiF 13.6 16.04 14.53 14.85 9.52KF 10.9 11.78 10.53 10.82 6.43NaCl 8.6 9.51 8.55 8.76 5.37MgO 7.77 8.86 7.91 8.10 4.86CaO 7.1 7.45 6.57 6.74 3.69 tional efforts to apply wide-range of materials. In fact,although their method was applied to calculate ionizationpotentials in Ref.45, it was not so satisfactory becausecalculations are performed in the combination of simplematerials (bulk calculations) with the supercell calcula-tions in GGA. Furthermore, no papers available to treattransition-metal oxides such as LaMnO in their method.We have to develop such an improved QSGW method ap-plicable to wide range of materials. Two requirements arethe computational efficiency and the theoretical validity.As a possibility to respect the efficiency, we can con-sider a hybridization method between QSGW and thedensity functional xc [14]. In QSGW80, a simple hy-bridization, 80 % QSGW+ 20 % GGA, we can see thatit works well for wide range of materials. Our present re-sults support the method of QSGW80, which takes onlythe 80 percent of QSGW self-energy. We can identifyQSGW80 as a simplification of the method in Ref.12.Ref.46 also presents one another approximation at thelevel of QSGW80 for the vertex correction in QSGW,resulting in similar good agreement with experiments.Let us examine how QSGW80 is justified for mate-rials calculated here. This is by the fact that η = ǫ ∞ (QSGW , RPA) /ǫ ∞ (QSGW , Slab) in Table I are ap-proximately 80 %, and show little material-dependency .Thus we expect that QSGW80 can mimic QSGW withthe vertex corrections; too large screened-exchange effectis reduced by the factor 0.8, with adding 0.2 GGA termso as to keep the total size of the xc term. In Table II,we show band gaps in QSGW and QSGW80 for materialstreated here. The band gaps are systematically too largein QSGW in comparison with experimental values [8],while QSGW80 gives rather better agreements with ex-perimental values. In Ref.8, we checked the performanceof the QSGW80 for ranges of materials. As in the case ofRef.12, QSGW80 is theoretically reasonable in the sensethat the band gaps are improved by using the corrected W . To go beyond QSGW80, we have to develop meth-ods to take the vertex correction into W as was donein Ref.12 in a simple manner. Considering the fact thatQSGW80 works well as shown in Ref.8, we may expectsimple methods to represent the vertex correction by ascalar factor or by limited number of parameters. As longas we know, the vertex correction can be relatively insen-sitive to materials, thus we expect some simple methodmight be available. V. SUMMARY
To clarify the importance of quasiparticle self-consistency in QSGW, we have explained the quasiparti-cle based perturbation in Sec.II. Then we emphasize theimportance of the self-consistency in the G W approxi-mation. Then we obtain the quasiparticles (independent-particle) given by H , and the interaction between thequasiparticles given by W . The vertex correction in QbPis introduced.We have performed QSGW calculations for slab mod-els under electric field by means of the ESM method. Thecalculated ǫ ∞ are in good agreements with experimentalvalues. Compared with ǫ ∞ in bulk calculation in RPA,we evaluated the size of vertex corrections as the func-tional derivative of the static self-energy in QSGW. Ourresults on ǫ ∞ give a support to the method by Shishkin, Marsman and Kresse [12]. As a simplified substitution oftheir method, we examined the performance of QSGW80[8] for materials treated here. The method QSGW+ESMdeveloped for the calculations should be useful even forother purposes such as bias-dependent spin susceptibilityin material theory, as well as practical device applicationsand materials designs.To go beyond usual QSGW, we should develop im-proved QSGW method in the G W approximation,whereas we should use accurate W by including the ver-tex correction. Then we have virtually best division of H = H + ( H − H ) where H gives the optimum inde-pendent particle picture. ACKNOWLEDGMENTS
T. K. thanks to supporting by JSPS KAKENHI GrantNumber 17K05499. We also thank the computing timeprovided by Research Institute for Informat ion Technol-ogy (Kyushu University). H. S. thanks to the comput-ing resources provided by the supercomputer system inRIKEN (HOKUSAI) and the supercomputer system inISSP (sekirei). [1] S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys.Rev. Lett. , 126406 (2004).[2] M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys.Rev. Lett. , 226402 (2006).[3] T. Kotani, M. van Schilfgaarde, and S. V. Faleev,Phys. Rev. B , 165106 (2007).[4] J. Heyd, G. E. Scuseria, and M. Ernzerhof,The Journal of Chemical Physics , 219906 (2006).[5] F. Tran and P. Blaha,Phys. Rev. Lett. , 226401 (2009).[6] J. He and C. Franchini,Physical Review B , 235117 (2012), arXiv: 1209.0486.[7] F. Bruneval and M. Gatti, in First Principles Approaches to Spectroscopic Properties of Complex Materials ,Vol. 347, edited by C. Di Valentin, S. Botti, and M. Co-coccioni (Springer Berlin Heidelberg, Berlin, Heidelberg,2014) pp. 99–135.[8] D. Deguchi, K. Sato, H. Kino, and T. Kotani,Jpn. J. Appl. Phys. , 051201 (2016).[9] H. Okumura, K. Sato, and T. Kotani,Physical Review B (2019), 10.1103/PhysRevB.100.054419.[10] Y. Lee, T. Kotan, and L. Ke,arXiv:2002.12417 [cond-mat] (2020), arXiv:2002.12417.[11] C. Bhandari, M. van Schilgaarde,T. Kotani, and W. R. L. Lambrecht,Phys. Rev. Materials , 013807 (2018).[12] M. Shishkin, M. Marsman, and G. Kresse,Phys. Rev. Lett. , 246403 (2007).[13] F. Bruneval, F. Sottile, F. Olevano, R. D. Sole, andL. Reining, Phys. Rev. Lett. , 186402 (2005).[14] A. N. Chantis, M. van Schilfgaarde, and T. Kotani, Phys.Rev. Lett. , 086405 (2006). [15] J. Otsuka, T. Kato, H. Sakakibara, and T. Kotani,Jpn. J. Appl. Phys. , 021201 (2017).[16] A. Sawamura, J. Otsuka, T. Kato, and T. Kotani,Journal of Applied Physics , 235704 (2017).[17] M. Otani and O. Sugino,Phys. Rev. B (2006), 10.1103/PhysRevB.73.115407.[18] L. Hedin and S. Lundqvist, Effects of Electron-Electronand Electron-Phonon interactions on the One-ElectronStates of Solids , Vol. 12 (Oxford university press NewYork, 1969).[19] D. Pines and P. Nozieres,
The Theory of Quantum Liquid.Vol I (W.A.Benjamin Inc., New York, 1966).[20] T. Kotani and M. van Schilfgaarde,Physical Review B , 125201 (2010),wOS:000276248900074.[21] A. L. Kutepov, Physical Review B , 155101 (2016).[22] A. L. Kutepov, Physical Review B , 195120 (2017).[23] F. Bechstedt, K. Tenelsen, B. Adolph, and R. D. Sole,Phys. Rev. Lett. , 1528 (1997).[24] Y. Takada, Molecular Physics , 1041 (2016).[25] N. L. Nguyen, H. Ma, M. Govoni, F. Gygi, and G. Galli,Phys. Rev. Lett. , 237402 (2019).[26] The ecalj package is available fromhttps://github.com/tkotani/ecalj/.[27] T. Kotani and M. van Schilfgaarde,Phys. Rev. B , 125117 (2010).[28] T. Kotani and H. Kino,J. Phys. Soc. Jpn. , 124714 [8 Pages] (2013),WOS:000327350700043.[29] T. Kotani, J. Phys. Soc. Jpn. , 094711 [11 Pages] (2014),WOS:000340822100029.[30] T. Kotani, H. Kino, and H. Akai,J. Phys. Soc. Jpn. , 034702 (2015). [31] D. Singh, Phys. Rev. B , 6388 (1991).[32] P. E. Blochl, Phys. Rev. B , 17953 (1994).[33] J. M. Soler and A. R. Williams,Phys. Rev. B , 1560 (1989).[34] J. M. Soler and A. R. Williams, Phys. Rev. B , 9728(1990).[35] J. M. Soler and A. R. Williams, Phys. Rev. B , 6784(1993).[36] N. Ashcroft, N. Mermin, and D. Wei, Solid State Physics:Revised Edition (Cengage Learning Asia, 2016).[37] M. J. L. Sangster and A. M. Stone-ham, Philosophical Magazine B , 597 (1981),https://doi.org/10.1080/01418638108222162.[38] J. Jacobson and E. Nixon,J. Phys. Chem. Solids , 967 (1968).[39] See a supplementary material, https://github.com/hskkbr/supplement4esmpaper.[40] S. Botti, F. Sottile, N. Vast, V. Ole-vano, L. Reining, H.-C. Weissker, A. Rubio,G. Onida, R. Del Sole, and R. W. Godby,Phys. Rev. B (2004), 10.1103/PhysRevB.69.155112.[41] D. Roessler and W. Walker,J. Phys. Chem. Solids , 1507 (1967).[42] J. E. Eby, K. J. Teegarden, and D. B. Dutton,Phys. Rev. , 1099 (1959).[43] D. M. Roessler and W. C. Walker,Phys. Rev. , 733 (1967).[44] R. Whited, C. J. Flaten, and W. Walker,Solid State Commun. , 1903 (1973).[45] A. Gr¨uneis, G. Kresse, Y. Hinuma, and F. Oba,Phys. Rev. Lett. , 096401 (2014).[46] W. Chen and A. Pasquarello,Phys. Rev. B92