GWΓ + Bethe-Salpeter equation approach for photoabsorption spectra: Importance of self-consistent GWΓ calculations in small atomic systems
aa r X i v : . [ c ond - m a t . o t h e r] S e p Physical Review B, Rapid Communications (2016), to appear. GW Γ + Bethe–Salpeter equation approach for photoabsorption spectra:Importance of self-consistent GW Γ calculations in small atomic systems Riichi Kuwahara,
1, 2
Yoshifumi Noguchi, and Kaoru Ohno ∗ Department of Physics, Yokohama National University,79-5 Tokiwadai, Hodogaya-ku, Yokohama 240-8501, Japan Dassault Syst`emes BIOVIA K.K., ThinkPark Tower,2-1-1 Osaki, Shinagawa-ku, Tokyo 140-6020, Japan Institute for Solid State Physics, The University of Tokyo,5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8581, Japan (Received 30 July 2015; Revised 29 June 2016)
Abstract
The self-consistent GW Γ method satisfies the Ward–Takahashi identity (i.e., the gauge invariance orthe local charge continuity) for arbitrary energy ( ω ) and momentum ( q ) transfers. Its self-consistent first-principles treatment of the vertex Γ = Γ v or Γ W is possible to first order in the bare ( v ) or dynamically-screened ( W ) Coulomb interaction. It is developed within a linearized scheme and combined with theBethe–Salpeter equation (BSE) to accurately calculate photoabsorption spectra (PAS) and photoemission(or inverse photoemission) spectra (PES) simultaneously. The method greatly improves the PAS of Na,Na , B , and C H calculated using the standard one-shot G W + BSE method that results in significantlyredshifted PAS by 0.8-3.1 eV, although the PES are well reproduced already in G W . v , and solve the QP equation,which is equivalent to the Dyson equation, as a self-consistent (SC) eigenvalue problem. TheHartree–Fock (HF) approach provides the first-order approximation. In Hedin’s set of equations[1] known as the GW Γ approach, the exchange-correlation part of the self-energy is expressed as Σ xc σ = iG σ W Γ σ , where G σ and Γ σ are the one-particle Green’s function and the vertex function ( σ isthe spin index), respectively, and W = (1 − vP ) − v represents the dynamically screened Coulombinteraction ( P = − i P σ G σ G σ Γ σ is the polarization function). The simplest approximation is toassume Γ σ =
1, which is called the GW approximation.It is well known that the SC GW method usually overestimates the energy gap [2, 3], whilethe one-shot GW approach ( G W ) using the Kohn–Sham (KS) wave functions and eigenvalues[4] results in a better energy gap. However, quite recently, it has been pointed out that the pho-toabsorption spectra (PAS) for small molecules obtained by solving the Bethe–Salpeter equation(BSE) [5, 6] using G W are often significantly redshifted by about 1 eV [7, 8]. The use of theHeyd-Scuseria-Ernzerhof (HSE) functional or the SC GW calculation (hereafter referred to asGW) improves the results, but they are not perfect [8, 9]. For a spin-polarized sodium atom (Na)and trimer (Na ), G W + BSE is extremely bad, although the G W QP energies are reasonablygood [10]. The calculated and experimental [11] optical gaps for Na are 1.32 eV and 2.10 eV,respectively, and the calculated and experimental [12] PAS for Na are shown in Fig. 1. Thesecalculated results are far o ff from the experimental data [13].Here, we develop a GW Γ method, which involves a SC treatment of the vertex Γ = Γ v or Γ W and satisfies the Ward–Takahashi identity [14–16] to first order in v or W , and show that itremarkably improves the QP energies and the optical gaps of spin-polarized Na, Na , B , andclosed-shell C H . In this method, the SC one-particle Green’s function, i.e., SC QP energies andwave functions are obtained in the GW Γ scheme. We use the all-electron mixed basis approach,in which single particle wave functions are expanded with both plane waves (PWs) and atomicorbitals (AOs) [10, 17]. This Rapid Communication reports a first-principles SC GW Γ calculationand its application to the BSE, which has never been performed so far except for some recentreports of non-SC GW calculations including the second-order screened exchange by Ren et al .2 IG. 1: Photoabsorption spectra of Na calculated using G W , GW, LGW, LGW Γ v , and LGW Γ W . Experi-mental data are taken from Ref. [12]. [18] and the GW Γ method (i.e., GW TC − TC + single-shot vertex correction for the self-energy withthe static approximation) by Gr¨uneis et al . [19]. All these authors used the KS, HF, or HSE wavefunctions throughout the calculations.In the present SC GW Γ +
BSE calculations, we will show the following: (1) Highly reliablePES and PAS are simultaneously obtained for every system. (2) All calculated results deviateby 0.1 eV at most from the available experimental data. (3) The failure of the G W + BSEcalculations for the PAS is caused by the use of localized KS wave functions above the vacuumlevel (VL), and hence accurate QP wave functions are required.Except for the G W and GW calculations, we use our recently developed technique [17] tolinearize the energy dependence of the self-energy Σ σ ( ǫ n ) to avoid the non-Hermitian problemcaused by the energy dependence and to perform fully SC calculations. The important point of thistechnique is that Λ σ = lim ( ω, q ) → Γ σ ( r , r , q ; µ + ω, µ ) = − ∂ Σ σ ( ω ) /∂ω | ω = µ is the vertex functionin the limit ( ω, q ) →
0. This is the Ward–Takahashi identity in the same limit. [Here, µ = ( ǫ HOMO + ǫ LUMO ) / H σ | n σ i = ǫ n σ Λ σ | n σ i with H σ = ˆ T + ˆ v nuc +Σ σ ( µ ) + µ ( Λ σ − L σ in the Choleski decomposition [20] Λ σ = L σ L † σ ,the renormalized QP states are given by | f n σ i = L † σ | n σ i , which satisfy e H σ | f n σ i = ǫ n σ | f n σ i e H σ = L − σ H σ L − † σ as well as the orthogonality and completeness conditions. Moreover, therenormalized Green’s function is given by e G σ ( ω ) = L † σ G σ ( ω ) L σ ; see Ref. [17] for more details. FIG. 2: Skeleton diagrams of the first-order vertex part (a) Γ v σ and (b) Γ W σ ; (c), (d), and (e) are thepolarization part P ; (f), (g), and (h) are the exchange-correlation part of the self-energy Σ xc σ ; (i), (j), and (k)are the interaction kernel e I σ σ ′ of the BSE. (c), (f), and (i) are usual diagrams without vertex correction;(a), (d), (g), (j), and (k) involve the first-order vertex in v (dotted line); (b), (e), and (h) involve the first-ordervertex in W (wavy line). (Theorem 1) In this linearized formulation, we can additionally introduce the vertex part Γ v σ ( r , r , q ; ǫ + ω, ǫ ) to first order in v [Fig. 2(a)], which we call the LGW Γ v method, or Γ W σ ( r , r , q ; ǫ + ω, ǫ ) to first order in W [Fig. 2(b)], which we call the LGW Γ W method. Thesevertex parts depend fully on the energy and momentum transfers ω and q , respectively, at thecenter (cross in those figures). See the Supplemental Material (SM) [21] for the proof of thistheorem.Then, the polarization function and the self-energy include the skeleton diagrams as shown inFigs. 2(c)-2(h). Figures 2(c) and 2(f) represent the diagrams without a vertex; Figs. 2(d) and 2(g)and Figs. 2(e) and 2(h) are the corresponding vertex corrections to first order in v ( Γ = Γ v ) and W ( Γ = Γ W ), respectively. Figure 3 illustrates the flow chart of the SC LGW Γ W method. The formsof the polarization function (Fig. 2(e)) and self-energy (Fig. 2(h)) are given in the SM.(Theorem 2) The present LGW Γ v and LGW Γ W methods satisfy the generalized Ward–Takahashi identity for arbitrary ω and q , which is equivalent to the gauge invariance (continuityequation for the electron density) [14–16], up to first order in v and W , respectively. The proof is4iven in the SM. FIG. 3: Flow chart of the SC LGW Γ W method. W and Γ W are replaced by v and Γ v in the SC LGW Γ v method. Recently, the BSE has been solved in the one-shot second-order approach [22]. Inwhat follows, we formulate the BSE for the LGW Γ v approach to spin-polarized sys-tems. In the linearized formulation, we use the renormalized two-particle Green’s function e S σ σ σ ′ σ ′ ( x , x ′ ; x , x ′ ) = L σ ′ L † σ S σ σ σ ′ σ ′ ( x , x ′ ; x , x ′ ) L σ L † σ ′ instead of S σ σ σ ′ σ ′ ( x , x ′ ; x , x ′ ), defined bysubtracting δ σ σ ′ δ σ σ ′ G ( x , x ′ ) G ( x , x ′ ) from the original two-particle Green’s function. It satisfiesthe BSE e S σ σ σ ′ σ ′ ( x , x ′ ; x , x ′ ) = e G σ ( x , x ) e G σ ′ ( x ′ , x ′ ) δ σ ′ σ ′ δ σ σ + X σ σ ′ Z e G σ ( x , x ) δ Σ σ σ ′ ( x , x ′ ) δ e G σ σ ′ ( x , x ′ ) e G σ ′ ( x ′ , x ′ ) ( L σ ′ L † σ ) − × e S σ σ σ ′ σ ′ ( x , x ′ ; x , x ′ ) d x d x ′ d x d x ′ , (1)where we used the fact that the original kernel Ξ σ σ σ ′ σ ′ ( x , x ′ ; x , x ′ ) is given by Ξ σ σ σ ′ σ ′ = δ Σ σ σ ′ δ G σ σ ′ = δ e G σ σ ′ δ G σ σ ′ δ Σ σ σ ′ δ e G σ σ ′ = L σ δ Σ σ σ ′ δ e G σ σ ′ L † σ ′ . (2)The functional derivative δ Σ σ σ ′ ( x , x ′ ) /δ e G σ σ ′ ( x , x ′ ) is given by − i δ σ σ ′ δ σ σ ′ δ ( x − x ′ ) δ ( x − x ′ ) v ( r − r ) + δ σ σ δ σ ′ σ ′ e I σ σ ′ ( x , x ′ ; x , x ′ ) . Ignoring all terms having functional derivatives of W [ e G ] with respect to e G as usual [23], we have e I σ σ ′ ( x , x ′ ; x , x ′ ) expressed as i { δ ( x − x ) δ ( x ′ − x ′ ) W ( x , x ′ ) + δ ( x − x )[ W Γ v σ ′ ]( x ′ , x ′ ; x ) + δ ( x ′ − x ′ )[ W Γ v σ ]( x , x ; x ′ ) } , which is representedby the skeleton diagrams of Figs. 2(i), 2(j), and 2(k). Here the last two terms (Figs. 2(j) and2(k)) come from vertex correction to first order in v . From these equations, we find that Λ σ = L σ L † σ should be multiplied to the polarization function as e P σ GG ′ = P σ GG ′ Λ σ [17]. Then, putting5 σ σ νµ dc = P G h g νσ | e i G · r | g µσ ih d σ | e − i G · r ′ | c σ i v ( G ) and using the expression for e I σ σ ′ νµ dc ( ω ) givenin the SM, we obtain the matrix eigenvalue equation of the BSE [23]( ǫ c ′ σ − ǫ d ′ σ ′ − Ω r ) A rd ′ σ ′ , c ′ σ = − occ X d emp X c δ σ σ ′ X σ V σ σ c ′ d ′ dc A rd σ , c σ − e I σ σ ′ c ′ d ′ dc ( Ω r ) A rd σ ′ , c σ (3)in the Tamm–Danco ff approximation [23]. We also use this formulation in the LGW Γ W approach,because the resulting error is on the order of 0.01 eV.For spin-polarized systems, we have to generally solve the eigenvalue equation (3) in the ↑↑ - ↓↓ subspace, h ↑↑ + v ↑↑ v ↑↓ v ↓↑ h ↓↓ + v ↓↓ A r ↑↑ A r ↓↓ = , (4)where we put h σ σ ′ = ( ǫ c ′ σ − ǫ d ′ σ ′ − Ω r ) δ cc ′ δ dd ′ − e I σ σ ′ c ′ d ′ dc and v σ σ = V σ σ c ′ d ′ dc .We used a face-centered cubic unit cell with edge length of 14 ˚A for Na and B , 15 ˚A forC H , and 18 ˚A for Na . All of the core and (truncated) valence numerical AOs are used togetherwith the PWs. The atomic geometries are optimized with DMol [24, 25]. The bond lengthsare 3.23 ˚A, 3.23 ˚A, and 5.01 ˚A for Na , 1.61 ˚A for B , 1.20 ˚A for C ≡ C, and 1.06 ˚A for C-H atthe Becke three-parameter Lee-Yang-Parr (B3LYP) functional level. We used 3.61 (50.76) Ry,1.23 (30.7) Ry, 6.82 (38.1) Ry, and 11.1 (44.2) Ry cuto ff energies for PWs (for Σ x σ ), respectively,for Na, Na , B , and C H . The cuto ff energy for P and Σ c σ is the same as that for PWs for Na andNa , and is set at 4.57 Ry for B and 3.98 Ry for C H . We used the full ω -integration [26] andthe projection operator for the GW-related calculations, but used the plasmon-pole model [27] and600 empty states for the Γ -related calculations as well as for solving the BSE in order to save thecomputational cost.The resulting ionization potential (IP), electron a ffi nity (EA), and optical gap E opt g (correspond-ing to the first dipole-allowed transition) of Na, Na , B , and C H calculated using the G W ,GW, LGW, and LGW Γ W methods are listed in Tables I and II, together with the results of previousmultireference single and double configuration interaction (MRDCI) calculations [28–34], config-uration interaction single and double (CID) calculations [35], and the corresponding experimentalvalues [11, 12, 36–45]. For Na and Na , the results of LGW Γ v are also listed in Table I. Notethat EA of C H is negative and not shown in Table II. Let us first compare the results of IP andEA with the experimental values. G W results in reasonable IP and EA (IP of C H is similarto those obtained in Ref. [46]) while GW has a tendency to overestimate IP and underestimate6A, although the experimental error bar is large for B . LGW improves GW [17], but is not per-fect. In contrast, LGW Γ v and LGW Γ W almost perfectly improve both IP and EA. For LGW Γ W ,the deviation from the experimental values is 0.03 eV for Na, 0.07 eV for Na , and 0.01 eV forC H . Compared with previous MRDCI calculations [28–32, 34], our results are closer to theexperimental IP and EA for almost all cases. TABLE I: Ionization potential (IP), electron a ffi nity (EA), and optical gap E opt g (corresponding to S → Pand B → A transitions) of Na and Na (in units of eV).Na Na IP EA E opt g IP EA E opt g G W Γ v Γ W a b b c / b b Expt. 5.14 d e f g / h i a Reference [28]. b Reference [29]. c Reference [30]. d Reference [36]. e Reference [37]. f Reference [11]. g Reference [38]. h Reference [39]. i Reference [12].
Next we compare the results of the optical gap E opt g with experiments. G W significantly under-estimates the experimental E opt g for all systems and GW overestimates the experimental E opt g . Thedeviation from the experimental values is 0.8-3.1 eV for G W and 0.13-0.26 eV for GW. LGWimproves the results except for Na ; the deviation from the experimental values is 0.08 eV for Na,0.27 eV for Na , 0.04 eV for B , and 0.07 eV for C H . In contrast, LGW Γ v and LGW Γ W give ex-cellent E opt g for all systems. For LGW Γ W , the di ff erence between the theoretical and experimentalvalues is less than 0.06 eV for Na and Na , 0.05 eV for B , and 0.09 eV for C H . Compared with7 ABLE II: IP, EA, and E opt g (corresponding to the Σ − g → Σ − u transition) of B , and IP and E opt g (corre-sponding to the Σ + g → Π u transition) of C H (in units of eV).B C H IP EA E opt g IP E opt g G W Γ W a b c d (8.06) e Expt. 10.3 ± . f ± . g h i j a Reference [31]. b Reference [32]. c CASSCF / MRDCI: Reference [33]. d Reference [34]. e CID: Reference [35]. f Reference [40]. g Reference [41]. h Reference [42, 43]. i Reference [44]. j Reference [45]. the experimental values, our E opt g is better than the previous MRDCI results for Na [29], and CIDresults for C H [35], or comparable to (di ff ers only by 0.01 eV from) previous MRDCI results forNa [29], and complete-active-space self-consistent-field (CASSCF) / MRDCI results for B [33].The LGW Γ W + BSE photoabsorption peak, i.e., the exciton wave function mainly consists of thefollowing QP hole and electron level pair(s): 6( s ↑ ) → p ↑ ) for Na, 17( s ↑ ) → p ↑ ) for Na ,4( σ ↑ ) → π ↑ ) and 3( σ ↓ ) → π ↓ ) for B , and 7( π ) → σ ) for C H . Figure 1 shows the PASof Na calculated using G W , GW, LGW, LGW Γ v , and LGW Γ W and the experimental spectra[12]. The overall spectral shapes are similar in all these methods except for G W , although thepeak positions are almost constantly shifted by an amount indicated by the di ff erence betweenthe calculated and experimental E opt g ’s in Table I, and the peak heights somewhat change afterthe inclusion of the vertex correction. Obviously, GW and LGW overestimate the peak positions,while LGW Γ v and LGW Γ W give good peak positions except for P and P . (LGW Γ v has a small8eak at 2.9 eV, which may correspond to P .) The remaining discrepancy between the theory andexperiment in the case of Na may be mainly attributed to the neglect of isomers and the atomicvibration e ff ects.Figure 4 shows the QP (or KS) energy spectra calculated using the local density approximation(LDA), G W , and LGW Γ W . The experimental IP and EA are indicated by IP and EA on the rightvertical border line. Now we discuss the reason why the PAS calculated using G W + BSE areso poor. For Na , the number of up-spin (down-spin) levels below the VL is 26 (26) for LDA,20 (19) for GW, 20 (20) for LGW, and 22 (20) for LGW Γ v and LGW Γ W . We confirmed that theKS and QP wave functions very much resemble each other for the first 20 levels below the VL.However, they are quite di ff erent for the QP levels above the VL. For example, the spin-up andspin-down KS wave functions at the 21st down-spin level and the 23rd up-spin level below the VLare depicted in Fig. 4. They are localized. However, the corresponding G W QP energies are bothabove the VL and the full QP wave functions are not localized. In our G W + BSE calculationof Na , the first small photoabsorption peak around 0.35 eV (see the top panel in Fig. 1) mainlyconsists of the QP hole and electron level pairs between 16 →
17 (19.8%), 21 (8.7%) for down-spin, and 17 →
18 (21.8%), 20 (31.4%), 21 (9.0%), 23 (2.3%) for up-spin. The unphysically boundKS wave functions of the 21 ↓ and 23 ↑ levels contribute to the BSE matrix elements, leading tounphysically large electron-hole interactions and in turn, to the optical transitions with very smallphotoabsorption energies. This gives the wrong spectra for G W in Fig. 1. It has already beenknown for more than 50 years [47] that the BSE should be solved with the fully SC Green’sfunction in order to satisfy the conservation laws as well as the longitudinal f -sum rule. However,the QP gap and the optical gap obtained using the GW method are blueshifted because they do notsatisfy the generalized Ward–Takahashi identity. To improve the result, it is necessary to use theGW Γ method.In this Rapid Communication, we presented the G W , GW, LGW, and LGW Γ W (LGW Γ v )calculations for Na, Na , B , and C H . If the G W QP energies are used together with the KSwave functions, there is inconsistency between the QP energies and wave functions at some levelsabove the VL. Moreover, the GW and LGW methods are not su ffi cient because they overestimateboth the QP energy gap and optical gap. To obtain better gap estimates, it is necessary to performthe GW Γ calculation. We showed that the LGW Γ W method produces consistent and the best PESand PAS among all of the methods used in this study. The self-consistent treatment of Γ is requiredto obtain consistently good results for both PES and PAS, and its computational cost scales as9 IG. 4: QP (or KS) energy spectra of Na calculated using the LDA, G W , and LGW Γ W . Unphysicallybound KS wave functions at the 21st spin-down level and the 23rd spin-up level are also depicted. Red ballsare Na atoms, while yellow and blue colors indicate the positive and negative regions of the wave function,respectively. O( N M ), where N and M are the numbers of basis functions and empty states, respectively, if weuse the plasmon-pole model. The present method is applicable to vertical transitions but cannothandle relaxation processes.This work was supported by the Grant-in-Aid for Scientific Research B (Grant No. 25289218)from JSPS and also by the Grant-in-Aid for Scientific Research on Innovative Areas (Grant No.25104713) from MEXT. We are also indebted to the HPCI promoted by MEXT for the use ofthe supercomputer SR16000 at Hokkaido University and at IMR, Tohoku University (Project IDs.hp140214, hp150231, hp160072, and hp160234). ∗ [email protected] [1] L. Hedin, Phys. Rev. , A796 (1965).[2] T. Kotani, M. van Schilfgaarde, and S. V. Faleev, Phys. Rev. B , 165106 (2007).[3] M. Shishkin, M. Marsman, and G. Kresse, Phys. Rev. Lett. , 246403 (2007).[4] M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. , 1418 (1985).[5] G. Onida, L. Reining, R. W. Godby, R. Del Sole, and W. Andreoni, Phys. Rev. Lett. , 818 (1995).[6] M. Rohlfing and S. G. Louie, Phys. Rev. Lett. , 3320 (1998).
7] D. Hirose, Y. Noguchi, and O. Sugino, Phys. Rev. B , 205111 (2015).[8] D. Jacquemin, I. Duchemin, X. Blase, J. Chem. Theory Comput. , 3290 (2015).[9] F. Bruneval, S. M. Hamed, and J. B. Neaton, J. Chem. Phys. , 244101 (2015).[10] Y. Noguchi, S. Ishii, K. Ohno, and T. Sasaki, J. Chem. Phys. , 104104 (2008).[11] J. E. Sansonetti, J. Phys. Chem. Ref. Data , 1659 (2008).[12] C. R. C. Wang, S. Pollack, T. A. Dahlseid, G. M. Koretsky, and M. M. Kappes, J. Chem. Phys. ,7931 (1992).[13] For the even number clusters of sodium, the G W + BSE result is not so bad; see Ref. [5] and: Y.Noguchi and K. Ohno, Phys. Rev. A , 045201 (2010).[14] J. C. Ward, Phys. Rev. , 182 (1950);[15] Y. Takahashi, Il Nuovo Cimento , 371 (1957).[16] J. R. Schrie ff er, Theory of Superconductivity (Advanced Book Classic) (Westview Press, 1999) Chap.8-5.[17] R. Kuwahara and K. Ohno, Phys. Rev. A. , 032506 (2014).[18] X. Ren, P. Rinke, G. E. Scuseria, and M. Sche ffl er, Phys. Rev. B , 035120 (2013).[19] A. Gr¨uneis, G. Kresse, Y. Hinuma, and F. Oba, Phys. Rev. Lett. , 096401 (2014).[20] It is also possible to use L = Λ / instead of the Choleski decomposition in the linearization procedureas in Ref. [3].[21] See Supplemental Material at http: // link.aps.org / supplemental / / PhysRevB.xx.xxxxxx for theproof of Theorem 1 and Theorem 2.[22] D. Zhang, S. N. Steinmann, and W. Yang, J. Chem. Phys. , 154109 (2013).[23] G. Strinati, Phys. Rev. Lett. , 1519 (1982).[24] B. Delley, J. Chem. Phys. , 508 (1990).[25] B. Delley, J. Chem. Phys. , 7756 (2000).[26] M. Zhang, S. Ono, N. Nagatsuka, and K. Ohno, Phys. Rev. B , 155116 (2016).[27] W. von der Linden and P. Horsch, Phys. Rev. B , 8351 (1988).[28] R. J. Buenker and S. Krebs, in Recent Advances in Multireference Methods , edited by K. Hirao, RecentAdvances in Computational Chemistry Vol. 4 (World Scientific, Singapore, 1999) pp.1-29.[29] V. Bonaˇci´c-Kouteck´y, P. Fantucci, and J. Kouteck´y, J. Chem. Phys. , 3802 (1990).[30] V. Bonaˇci´c-Kouteck´y, I. Boustani, M. Guest, and J. Kouteck´y, J. Chem. Phys. , 4861 (1988).[31] P. J. Bruna and J. S. Wright, J. Phys. Chem. , 1774 (1990).
32] P. J. Bruna and J. S. Wright, J. Phys. B , 2197S (1990).[33] S. R. Langho ff and C. W. Bauschlicher, Jr., J. Chem. Phys. , 5882 (1991).[34] W. P. Kraemer and W. Koch, Chem. Phys. Lett. , 631 (1993).[35] W. E. Kammer, Chem. Phys. , 408 (1974).[36] A. Herrmann, E. Schumacher, and L. W¨oste, J Chem. Phys. , 2327 (1978).[37] H. Hotop and W. C. Lineberger, J. Phys. Chem. Ref. Data , 731 (1985).[38] A. Herrmann, S. Leutwyler, E. Schumacher, and L. Woste, Helv. Chim. Acta , 453 (1978).[39] K. M. McHugh, J. G. Eaton, G. H. Lee, H. W. Sarkas, L. H. Kidder, J. T. Snodgrass, M. R. Manaa,and K. H. Bowen, J. Chem Phys. , 3792 (1989).[40] L. Hanley, J. L. Whitten, and S. L. Anderson, J. Phys. Chem. , 5803 (1988).[41] C. J. Reid, Int. J. Mass Spect. Ion Proces. , 147 (1993).[42] W. R. M. Graham and W. Weltner Jr., J. Chem. Phys. , 1516 (1976).[43] K. P. Huber and G. Herzberg, Molecular Spectra and Molecular Structure: IV. Constants of DiatomicMolecules (Van Nostrand Reinhold Company, New York, 1979).[44] G. Bieri and L. ˚Asbrink, J. Elec. Spec. Rel. Phenom. , 149 (1980).[45] W. C. Price, Phys. Rev. , 444 (1935).[46] M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren, M. Sche ffl er, F. Liu, J. Lischner, L. Lin, J. R.Deslippe, S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Evers, and P. Rinke, J. Chem. TheoryComput. , 5665 (2015).[47] G. Baym and L. P. Kadano ff , Phys. Rev. , 287 (1961). upplemental Material * * * * * * GW Γ + Bethe–Salpeter equation approach for photoabsorption spectra:Importance of self-consistent GW Γ calculations in small atomic systems * * * * * * *Riichi Kuwahara, , Yoshifumi Noguchi, and Kaoru Ohno , ∗ DepartmentofPhysics,YokohamaNationalUniversity,79-5 Tokiwadai,Hodogaya-ku,Yokohama240-8501,Japan DassaultSyst`emesBIOVIA K.K.,ThinkParkTower,2-1-1 Osaki,Shinagawa-ku,Tokyo140-6020,Japan InstituteforSolidState Physics,TheUniversityofTokyo,5-1-5 Kashiwanoha,Kashiwa, Chiba277-8581,JapanIn this Supplemental Material, we prove Theorem 1 and Theorem 2, and present explicit formsof the polarization function, the self-energy, and the interaction kernel e I σ σ ′ νµ dc ( ω ). Figure and refer-ence numbers refer to the figures and references in the Letter except for Ref [S1].Theorem 1. In the linearized formulation, one can additionally introduce the vertex part Γ v σ ( r , r , q ; ǫ + ω, ǫ ) at the first order in the bare Coulomb interaction v (Fig. 2(a)), which we will call theLGW Γ v method, or Γ W σ ( r , r , q ; ǫ + ω, ǫ ) at the first order in the dynamically screenedCoulomb interaction W (Fig. 2(b)), which we will call the LGW Γ W method.Theorem 2. The LGW Γ v and LGW Γ W methods satisfy the Ward–Takahashi identity to the first order inthe bare Coulomb interaction v and in the dynamically screened Coulomb interaction W ,respectively. 13 ROOF OF THEOREM 1
The statement of Theorem 1 holds because using e G σ ( ω ) in place of G σ ( ω ) in the linearizedformulation introduces the ω = q = Λ σ just at the other side of the interactionline where we introduced the first-order vertex part Γ v σ ( r , r , q ; ǫ + ω, ǫ ) or Γ W σ ( r , r , q ; ǫ + ω, ǫ )for arbitrary ω and q . Therefore there is no double counting in the vertex correction up to the firstorder in v or W . In this way, the LGW Γ v and LGW Γ W methods rigorously treat the vertex parts tothe first order in v (Fig. 2(a)) and in W (Fig. 2(b)), respectively. These vertex parts depend fullyon the energy and momentum transfers ω and q at the center (cross) as well as the frequencies andthe coordinates at both ends; see Figs. 2(a) and (b).Moreover, in the LGW Γ v method, it is possible to show that there is no interference betweenthe first-order vertex part of Fig. 2(a) (i.e., the vertex part at the first order in v ) and Λ σ in the ω = q = ω = q = − ω derivative of the ω -independent Fock exchange self-energy Σ x σ and hence exactly equals to zero, i.e., Γ v σ ( r , r , q = ǫ, ǫ ) = − ∂ Σ x σ /∂ω =
0. Therefore, the full vertex part Λ σ in the ω = q = Γ W method, however, the first-order vertex part of Fig. 2(b) (i.e., the vertex part at the first order in W ) may interfere with the full vertex part Λ σ in the ω = q = PROOF OF THEOREM 2
Ward–Takahashi identity [16-18] is given by δ ( x − x ) G − ( x , x ) − G − ( x , x ) δ ( x − x ) = i ∇ · Γ ( x , x , x ) + i ∂∂ t Γ ( x , x , x ) . (S.5)If we Fourier transform this equation with respect to t − t and t − t into ǫ + ω and ǫ , respectively,(i.e., if we multiply exp[ i ( ǫ + ω )( t − t ) + i ǫ ( t − t )] and integrate with respect to t − t and t − t ),this equation becomes δ ( r − r ) G − ( r , r ; ǫ + ω ) − G − ( r , r ; ǫ ) δ ( r − r ) = i ∇ · Γ ( r , r , r ; ǫ + ω, ǫ ) + ω Γ ( r , r , r ; ǫ + ω, ǫ ) . (S.6)14t the lowest order, the first term in the right hand side is approximately given by i ∇ · Γ ( r , r , r ; ǫ + ω, ǫ ) ∼ ∇ · ( ∇ + ∇ ) δ ( r − r ) δ ( r − r ) = − δ ( r − r ) 12 ∇ δ ( r − r ) + ∇ δ ( r − r ) δ ( r − r ) − i Z ∞−∞ d ω ′ π G ( r , r ; ǫ − ω ′ ) − ∇ ! G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π − ∇ ! G ( r , r ; ǫ − ω ′ ) G ( r , r ; ǫ + ω − ω ′ ) W ( r , r , ω ′ ) . (S.7)Substituting this into (S.6), we have δ ( r − r ) G − ( r , r ; ǫ + ω ) − G − ( r , r ; ǫ ) δ ( r − r ) = − δ ( r − r ) 12 ∇ δ ( r − r ) + ∇ δ ( r − r ) δ ( r − r ) − i Z ∞−∞ d ω ′ π G ( r , r ; ǫ − ω ′ ) − ∇ ! G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π − ∇ ! G ( r , r ; ǫ − ω ′ ) G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + ω Γ ( r , r , r ; ǫ + ω, ǫ ) , (S.8)which is eqivalent to Σ ( r , r ; ǫ ) δ ( r − r ) − δ ( r − r ) Σ ( r , r ; ǫ + ω ) = − i Z ∞−∞ d ω ′ π G ( r , r ; ǫ − ω ′ ) − ∇ ! G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π − ∇ ! G ( r , r ; ǫ − ω ′ ) G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + ω [ Γ ( r , r , r ; ǫ + ω, ǫ ) − δ ( r − r ) δ ( r − r )] . (S.9)In the linarized formulation, this equation can be rewritten as Σ ( r , r ; ǫ ) δ ( r − r ) − δ ( r − r ) Σ ( r , r ; ǫ + ω ) = − i Z ∞−∞ d ω ′ π e G ( r , r ; ǫ − ω ′ ) − ∇ ! e G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π − ∇ ! e G ( r , r ; ǫ − ω ′ ) e G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + ω [ Γ ( r , r , r ; ǫ + ω, ǫ ) − δ ( r − r ) δ ( r − r )] . (S.10)15n order to prove this equality, we calculate the left hand side of (S.10) within the linearizedformulation as follows: i Z ∞−∞ d ω ′ π [ e G ( r , r ; ǫ − ω ′ ) δ ( r − r ) − δ ( r − r ) e G ( r , r ; ǫ + ω − ω ′ )] W ( r , r ; ω ′ ) = i Z ∞−∞ d ω ′ π h r | ǫ − ω ′ − e H − i δ e H | r ih r | r i − h r | r ih r | ǫ + ω − ω ′ − e H − i δ e H | r i W ( r , r ; ω ′ ) = i Z ∞−∞ d ω ′ π h r | ǫ − ω ′ − e H − i δ e H | r ih r | ǫ + ω − ω ′ − e H − i δ e H ǫ + ω − ω ′ − e H − i δ e H | r i W ( r , r ; ω ′ ) − i Z ∞−∞ d ω ′ π h r | ǫ − ω ′ − e H − i δ e H ǫ − ω ′ − e H − i δ e H | r ih r | ǫ + ω − ω ′ − e H − i δ e H | r i W ( r , r ; ω ′ ) = i Z ∞−∞ d ω ′ π h r | ǫ − ω ′ − e H − i δ e H | r ih r | ωǫ + ω − ω ′ − e H − i δ e H | r i W ( r , r ; ω ′ ) − i Z ∞−∞ d ω ′ π h r | ǫ − ω ′ − e H − i δ e H | r ih r | e H ǫ + ω − ω ′ − e H − i δ e H | r i W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π h r | e H ǫ − ω ′ − e H − i δ e H | r ih r | ǫ + ω − ω ′ − e H − i δ e H | r i W ( r , r ; ω ′ ) = i ω Z ∞−∞ d ω ′ π e G ( r , r ; ǫ − ω ′ ) e G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) − i Z ∞−∞ d ω ′ π e G ( r , r ; ǫ − ω ′ ) − ∇ ! e G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π − ∇ ! e G ( r , r ; ǫ − ω ′ ) e G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) − i Z ∞−∞ d ω ′ π Z d r ′ e G ( r , r ; ǫ − ω ′ ) Σ xc ( r , r ′ ) e G ( r ′ , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) + i Z ∞−∞ d ω ′ π Z d r ′ Σ xc ( r , r ′ ) e G ( r ′ , r ; ǫ − ω ′ ) e G ( r , r ; ǫ + ω − ω ′ ) W ( r , r ; ω ′ ) , (S.11)where we used the fact that the electron-nucleus potential in e H commutes with the operator | r ih r | .The first term of Eq. (S.11) is exactly equal to the vertex part at the lowest order in the dynamicallyscreened Coulomb interaction W shown in Fig. 2(b) except for the prefactor ω , and thus equalsto the last two terms of the right hand side of Eq. (S.10). The second and the third terms equalto the first and second terms of the right hand side of Eq. (S.10). The fourth and fifth terms areat least one order higher in v compared to the other terms that are lowest order in W , and can beignored. Therefore, the LGW Γ W method satisfies the Ward–Takahashi identity to the lowest orderin W . This discussion holds also in the case we expand to the first order in the electron-electronCoulomb interaction v of Fig. 2(a), since W ( r , r ; ω ′ ) can be replaced by v ( r − r ) in this case.Therefore, the LGW Γ v method also satisfies the Ward–Takahashi identity to the lowest order in v .16 OLARIZATION FUNCTION AND SELF-ENERGY
The contributions to the polarization function and self-energy coming from the first-order ver-tex correction, FIGs. 2(e) and (h) are explicitly given by e P (1) GG ′ ( ω ) = Ω Z d ω ′ π e i ω ′ η Z d ω ′′ π e i ω ′′ η X σ X i jkl X G ′′ G ′′′ W G ′′ G ′′′ ( ω ′′ ) × h e i σ | e − i G · r | f j σ ih f j σ | e − i G ′′ · r ′′ | f k σ ih f k σ | e i G ′ · r ′ | e l σ ih e l σ | e i G ′′′ · r ′′′ | e i σ i ( ω ′ − ǫ i σ − i η i σ )( ω + ω ′ − ǫ j σ − i η j σ )( ω ′ + ω ′′ + ω − ǫ k σ − i η k σ )( ω ′ + ω ′′ − ǫ l σ − i η l σ ) , (S.12)and ( α | Σ (2) σ ( ω ) | β ) = Z d ω ′ π e i ω ′ η Z d ω ′′ π e i ω ′′ η X σ X i jk X GG ′ X G ′′ G ′′′ W GG ′ ( ω ′ ) W G ′′ G ′′′ ( ω ′′ ) × ( α | e i G · r | e i σ ih e i σ | e i G ′ · r ′ | f j σ ih f j σ | e − i G ′′ · r ′′ | f k σ ih f k σ | e − i G ′′′ · r ′′′ | β )( ω + ω ′ − ǫ i σ − i η i σ )( ω + ω ′ + ω ′′ − ǫ j σ − i η j σ )( ω + ω ′′ − ǫ k σ − i η k σ ) , (S.13)respectively, where W GG ′ ( ω ) is the dynamically screened Coulomb interaction, η is a positiveinfinitesimal ( η i σ is η for occupied i σ and − η for empty i σ ), Ω is the volume of the unit cell, and( α | , | β ) denote basis functions. INTERACTION KERNEL e I σ σ ′ νµ dc ( ω ) As written in the Letter, the interaction kernel e I σ σ ′ ( x , x ′ ; x , x ′ ) is expressed as i { δ ( x − x ) δ ( x ′ − x ′ ) W ( x , x ′ ) + δ ( x − x )[ W Γ v σ ′ ]( x ′ , x ′ ; x ) + δ ( x ′ − x ′ )[ W Γ v σ ]( x , x ; x ′ ) } , which arediagrammatically represented by Figs. 2(i), (j), and (k), respectively. In the third term correspond-ing to Fig. 2(k), for example, the two points connected by the dotted line, i.e., the bare Coulombinteraction, are the same time t = t . On the other hand, one end of the wavy line, which is in-volved in this vertex part Γ v σ is another time, say t ′′ , and the other end of the wavy line, which isnot involved in the vertex part is obviously a unique time t ′ = t ′ . Therefore, this interaction kernel e I σ σ ′ ( x , x ′ ; x , x ′ ) has only the time dependence through τ = t − t ′ . Moreover, the vertex part Γ v σ in this diagram has only the time dependence through τ ′ = t − t ′′ , and the wavy line W has onlythe time dependence through τ ′′ = t ′′ − t ′ = τ − τ ′ . Therefore, the interaction kernel e I σ σ ′ ( τ ), i.e.,the dynamically screened Coulomb interaction W ( τ − τ ′ ) times the vertex correction Γ v σ ( τ ′ ), has17 convolution type in time, so that its Fourier transformation from τ to ω ′ is written as e I σ σ ′ ( ω ′ ),which can be written as the product of the Fourier transforms of W ( τ ′′ ) and Γ v σ ( τ ′ ), W ( ω ′ ) and Γ v σ ( ω ′ ); thus we have e I σ σ ′ ( ω ′ ) = W ( ω ′ ) Γ v σ ( ω ′ ). Since the ω ′ -dependence in e I σ σ ′ ( ω ′ ) can betreated in a standard way given first by Strinati [S1] when we solve the Bethe–Salpeter equation,it is allowed to replace e I σ σ ′ ( ω ′ ) with W ( ω ′ ) Γ v σ ( ω ′ ) in the final equation (2.18) of Ref. [S1], andthus we can derive e I σ σ ′ νµ dc ( ω ) = i X GG ′ Z d ω ′ π e − i ω ′ η W GG ′ ( ω ′ ) ( h g νσ | e i G · r | c σ ih d σ ′ | e − i G ′ · r ′ | g µσ ′ i + X αβ (cid:0) δ α occ δ β emp − δ α emp δ β occ (cid:1) h g νσ | e i G · r | c σ i V σ ′ σ ′ d αβµ h g ασ ′ | e − i G ′ · r ′ | g βσ ′ i ǫ ασ ′ − ǫ βσ ′ + ω ′ + i η ασ ′ + V σ σ ναβ c h g ασ | e i G · r | g βσ i ǫ ασ − ǫ βσ − ω ′ + i η ασ h d σ ′ | e − i G ′ · r ′ | g µσ ′ i !) × ( ω − ω ′ − ǫ c σ + ǫ µσ ′ + i η + ω + ω ′ − ǫ νσ + ǫ d σ ′ + i η ) (S.14)with V σ σ νµ dc = P G h g νσ | e i G · r | g µσ ih d σ | e − i G · r ′ | c σ i v ( G ) (note that all states are renormalizedQP states except for d and c ).Reference:[S1] G. Strinati, Phys. Rev. B29