Table Of Content

Exploring approximations to the GW self-energy ionic gradients C. Faber1,2,3, P. Boulanger1,2, C. Attaccalite1,2, E. Cannuccia4,5, I. Duchemin2,3, T. Deutsch2,3, X. Blase1,2 1CNRS, Institut Néel, 38042 Grenoble, France, 2Université Grenoble Alpes, 38042 Grenoble, France, 3INAC, SP2M/L sim, CEA cedex 09, 38054 Grenoble, France. 4Institut Laue Langevin, BP 156, 38042 Grenoble, France, 5Laboratoire de Physique des Interactions Ioniques et Moléculaires, Aix-Marseille Université/CNRS, Campus de Saint Jéroˆme, 13397 Marseille, France. 5 (Dated: March 12, 2015) 1 0 Theaccuracyofthemany-bodyperturbationtheoryGW formalismtocalculateelectron-phonon 2 coupling matrix elements has been recently demonstrated in the case of a few important systems. n However, the related computational costs are high and thus represent strong limitations to its a widespread application. In the present study, we explore two less demanding alternatives for the J calculation of electron-phonon coupling matrix elements on the many-body perturbation theory 8 level. Namely, we test the accuracy of the static Coulomb-hole plus screened-exchange (COHSEX) 2 approximation and further of the constant screening approach, where variations of the screened CoulombpotentialW uponsmallchangesoftheatomicpositionsalongthevibrationaleigenmodes ] areneglected. Wefindthislatterapproximationtobethemostreliable,whereasthestaticCOHSEX h ansatz leads to substantial errors. Our conclusions are validated in a few paradigmatic cases: p diamond, graphene and the C fullerene. These findings open the way for combining the present - 60 p many-body perturbation approach with efficient linear-response theories. m PACSnumbers: 71.15.Qe,71.38.-k o c . s I. INTRODUCTION are not yet available within the framework of MBPT. c Existing GW calculations of EPC matrix elements have i s Electron-phononcoupling(EPC)occupiesaprominent beenthereforebasedonthefrozen-phononapproach,ne- y cessitating a stepwise displacement of the atoms along roleinvariousfieldsofcondensedmatterphysics,includ- h thephononmodeswithanexplicitevaluationoftheelec- p ing phonon-mediated superconductivity, photoemission tronic structure on the GW level for each step. This de- [ bandgaprenormalization,currentcarriersinelasticscat- manding approach cannot be reconciled with very large teringorthelifetimeofhotelectrons. Concerningthecal- 1 unit cells and is thus only feasible for zone-center or culationofEPCmatrixelementsontheabinitiolevel,up v tonow, mainlydensityfunctionaltheory(DFT)1 andits zone-boundaryphononmodes. Goingbeyondthefrozen- 8 5 perturbative linear response extensions (DFPT)2,3 have phonon approach, in order to access EPC matrix elements on the GW level for a large number of electron 0 beenapplied,providingimportantinformationatthemi- 7 croscopic level. and phonon wave vectors, remains a considerable chal- 0 lengetotheab initio community. Significantworkisstill Recently,severalstudiesquestionedtheaccuracyofthe . neededtocomeupwithaschemethatmayallowtogen- 1 DFT-based EPC matrix elements when obtained with eralizethepioneeringabove-mentionedstudiestoamuch 0 (semi)local functionals such as LDA or PBE.4 By way 5 larger set of systems. of example, the electron-phonon coupling involving spe- 1 cific electronic and phonon modes in graphene5,6 and In the present study, we explore the merits : v graphite,7 the value of the electron-phonon coupling po- of the static coulomb-hole plus screened exchange Xi tential to states at the Femi level in the electron-doped (COHSEX)14,15,19–27 approximation to the GW self- fullerene,8,9superconductingbismuthatesandtransition- energy for the calculation of the electron-phonon cou- r a metalchloronitrides,10 ortherenormalizationofthepho- pling matrix elements in the isolated fullerene molecule, toemission band structure of pentacene11 and diamond in diamond and in graphene. Our definitions of the rele- crystals,12 were shown to be affected by a significant un- vantEPCmatrixelementsfollowpreviousstudiesrelated derestimation of the EPC matrix elements when calcu- to the superconducting transition in the fullerides,8,9 lated within DFT and (semi)local functionals. the zero-point motion renormalization of the gap in As a cure to such problems, many-body perturba- diamond12 andthegapopeningthroughelectron-phonon tiontheory(MBPT)techniqueswithintheso-calledGW coupling in graphene,5 respectively. For the same set of approximation13–18 showed a clear improvement when systems, we also explore the accuracy of the constant- compared to available experimental data.5,7,9–12 Unfor- screening approximation, namely the assumption that tunately, the kind of theories that are available within thescreenedCoulombpotentialcanbeconsidered,tolin- DFT,andinparticularthepowerfulDFPTformalism,2,3 earorder,asaconstantuponsmallchangesoftheatomic 2 positions. WhereasCOHSEXleadstonon-negligibledis- crepancies,thislatterapproximationwillbeshowntobe occ robust and accurate. The present results offer impor- ΣSEX(r,r(cid:48);E) = (cid:88)W(r,r(cid:48);E ε )φ (r)φ∗(r(cid:48)) tantperspectivesinmakingtheGW formalismamenable − − n n n n to the study of condensed matter phenomena involving (cid:88) (cid:90) ∞ B(r,r(cid:48);ω(cid:48)) electron-phonon coupling. ΣCOH(r,r(cid:48);E) = φ (r)φ∗(r(cid:48)) dω(cid:48) , n n P E ε ω(cid:48) n 0 − n− II. METHODOLOGY where indicates the principal value. ΣSEX, which in- P volves a summation over the occupied states only, originates from the poles of G. It is called the screened ex- Webrieflyintroducethemany-bodyperturbationthe- change interaction in analogy to the bare exchange term ory Green’s function formalism, providing a solid frame- that can be obtained by replacing W with the energy- work for the calculation of quasiparticle energies E. In independentbareCoulombpotentialv. ΣCOH originates such an approach, the one-body quasiparticle eigenvalue from the poles of W and represents the Coulomb-hole equation reads: contribution, since it can be shown to be related to the interaction of an electron with its related adiabatically (cid:18) 2 (cid:19) built correlation hole. −∇ +Vion(r) φ(r)+VH(r)φ(r) 2 The static approximation to the exact COHSEX de- (cid:90) compositionassumesthat(E εn) 0forall(n),leading + dr(cid:48) Σ(r,r(cid:48);E)φ(r(cid:48))=Eφ(r), to simplified static screened e−xchan(cid:39)ge and Coulomb-hole expressions, namely: where Vion and VH are the ionic and Hartree potential, respectively. The self-energy Σ(r,r(cid:48);E) replaces occ (cid:88) the well-known exchange-correlation potential of den- ΣSEX(r,r(cid:48);0) = W(r,r(cid:48);0)φ (r)φ∗(r(cid:48)) static − n n sity functional theory or the exchange operator in the n Hartree-Fock formalism. In general, it is non-local, 1(cid:88) energy-dependent and non-Hermitian. Derived within ΣCstOatHic(r,r(cid:48);0) = 2 φn(r)φ∗n(r(cid:48))W(cid:102)(r,r(cid:48);0) Schwinger’s functional derivative approach to perturba- n 1 tion theory,13 the GW approximation to the self-energy = δ(r r(cid:48))W(cid:102)(r,r(cid:48);0). leads to: 2 − Here, W(cid:102) =(W v)isthedifferencebetweenthescreened i (cid:90) and bare Coulo−mb potential. As a result, besides being ΣGW(r,r(cid:48);E)= dωeiδωW(r,r(cid:48);ω)G(r,r(cid:48);E+ω), 2π a static approximation, this Coulomb-hole term is also local in space. Such a static COHSEX approximation whereGandW arethetime-orderedone-particleGreen’s (labeled COHSEX here below) was shown to yield too function and the dynamically screened Coulomb poten- large gaps in the case of semiconductors.15 By way of tial, respectively, and δ = 0+ a small positive infinitesi- example, in the present case of the C molecule, the 60 mal. COHSEX gap is found to be 5.3 eV, i.e. 0.4 eV larger TheCoulomb-holeplusscreened-exchange(COHSEX) than the 4.9 eV experimental gap.67 Nevertheless, this representationoftheGW self-energywasoriginallyintro- has to be∼compared to the starting 1.6 eV DFT-LDA ducedusingatimerepresentationofG, W andΣ.14,20,21 Kohn-Sham gap which is dramatically too small. We follow Hybertsen and Louie15 by using the following While it cannot be claimed that the static COHSEX spectral representations of G and W: approach is a good approximation to absolute quasiparticle energies, we emphasize that we are interested in quasiparticle energy differences upon small (infinitesi- G(r,r(cid:48);E+ω) = (cid:88) φn(r)φ∗n(r(cid:48)) mal) atomic lattice motions. The main assumption on E+ω ε iδsgn(ε µ) n − n− n− which we rely to calculate the electron-phonon coupling (cid:90) ∞ 2ω(cid:48)B(r,r(cid:48);ω(cid:48)) within the COHSEX approximation is that the varia- W(r,r(cid:48);ω) = v(r,r(cid:48))+ dω(cid:48) , ω2 (ω(cid:48) iδ)2 tions of the dynamical contribution to the self-energy 0 − − can be neglected. This can be a priori rationalized by where G has been written in terms of one-body eigen- emphasizing that dynamical interactions are driven by states φ and eigenenergies ε , typically starting DFT the plasmons dynamics, collective excitations less sensi- n n Kohn-Sham or Hartree-Fock solutions, and W in terms tive to small atomic displacements than single-particle of its spectral function B. These expressions allow to excitation energies and wave functions. It remains, how- obtain the pole structure of both G and W. From the ever, that besides the static approximation, the spatially residue theorem, one then rapidly obtains an exact de- local character of the static COH term is at odds with composition of Σ: the nonlocality of the full GW self-energy. 3 Inordertofurtherjustifythefollowingresultsconcern- the superconducting transition temperature through the ing the constant screening approach, namely the second dimensionless parameter λ = N(E )Vep, where N(E ) F F approximation we explore in this study, it is instructive is the density of states at the Fermi level. Namely, Vep to consider the GW plus Bethe-Salpeter formalism.28–30 reads in the so-called molecular limit:44,45 Thelatterisamany-bodyperturbationtheoryapproach concerned with describing the linear response of a sys- tteiomnwanitdhtrheuspsetchtetoBaetthiem-Se-adlpepeteenrdeeqnutaetxitoenrn(aBlSpEe)rtiusrtbhae- Vep = 9M1 (cid:88)ω12 (cid:88)3 (cid:12)(cid:12)(cid:12)(cid:12)∂∂(cid:15)umν(cid:12)(cid:12)(cid:12)(cid:12)2. (1) MBPT analogue to time-dependent DFT. At the heart ν ν m=1 of the GW/BSE approach lies the variation (∂GW/∂λ), Here, M is the mass of the carbon ions, ω is the fre- wherethe“perturbation”λistheone-bodyGreen’sfunc- ν quency of the vibrational mode with index (ν) and uν is tion G. The most common approximation that has been the vibrational polarization vector. As already discussed shown to be remarkably accurate31–33 is to replace the before, the electronic states with energy ((cid:15) ) are lim- GW self-energy by its static COHSEX approximation, m ited to the 3-fold degenerate LUMO level. From group and to consider further that (∂W/∂λ) = 0. It is such a theory analysis it follows that only the H and A vi- simplifiedschemeweaimtoexploreinthepresentstudy, g g brational modes can couple to these states. A schematic differing in the fact that the perturbation λ is now in- representation of the t level splitting as a function of duced by a vibrational distortion of the system. 1u the deformation amplitude along a H mode is provided g in Fig. 1(a). The EPC matrix elements are consequently relatedtotheslopes(∂(cid:15) /∂uν),showingastrongdepen- m III. TECHNICAL DETAILS denceontheformalismadoptedtocalculatetheenergies (cid:15) . These energy derivatives are calculated within the m The many-body GW and COHSEX calculations for frozen phonon approach, using a symmetric five points theisolatedC fullereneareperformedusingtheFiesta finite-difference formula. 60 code, an implementation of the GW formalism within Such an effective electron-phonon potential was calcu- a Gaussian basis.34–36 We start from DFT-LDA eigen- latedforC inRefs.8,9,bothwithinDFTusingdiverse 60 statescalculatedwiththeSiestapackage37 andatriple- (semi)local and hybrid exchange-correlation functionals zeta plus polarization (TZP) basis38 for the description and GW. While DFT-LDA values were shown to signif- of the valence orbitals combined with standard norm- icantly underestimate the coupling energy, DFT calcula- conserving pseudopotentials.39 We use the resolution of tionsperformedwithglobalhybridscontaining20%-30% the identity technique (RI-SVS) with an even-tempered ofexactexchangeandtheGW approachwerefound8,9 to auxiliary basis set formed of four Gaussians per angular compare favorably to gas phase experimental data.46,47 channel (up to d character) with exponents between 0.2 As a drawback of DFT calculations with hybrids, it was and 3.2 Bohr−2. With such a basis, our G W @LDA evidenced thatthe resulting Vep potential wouldquickly 0 0 ionization potential and HOMO-LUMO gap are found increasewiththeamountofexactexchange, yieldingthe to be 7.33 eV and 4.39 eV (B3LYP geometry), respec- usual question of the proper choice of the functional pa- tively, starting from a 1.66 eV DFT-LDA gap. This is in rameters for a given system. closeagreementwiththe7.45eVand4.40eVplanewave In the present study, we are concerned with explor- G W @LDA values of Ref. 40 (LDA geometry). The ing simplified GW schemes, i.e approximations reducing 0 0 G W @LDA gap remains smaller than the 4.9 eV ex- the needed computational effort. As such, quasiparti- 0 0 ∼ perimental value, an issue that will be addressed by the cle energies are evaluated at a single-shot COHSEX and partiallyself-consistentschemeusedinthepresentstudy, a partially self-consistent evCOHSEX level, where the as discussed below. quasiparticle energies are reinjected self-consistently in Following Ref. 8, we use the relaxed structure and the construction of G and W, while keeping the start- phonon eigenmodes generated within a DFT-B3LYP ap- ing Kohn-Sham LDA wave functions unchanged. Our proach using a 6-311G(d) basis. This was shown to pro- reference point are partially self-consistent evGW calcu- vide vibrational frequencies in excellent agreement with lations with self-consistency on the eigenvalues as pre- Raman measurements.8,41 Concerning the calculation of sented in Ref. 9 for the calculation of the EPC matrix the EPC matrix elements, our present study follows elements in C . Such a simplified approach to full self- 60 the approach described in Refs. 8,9, where the electron- consistency,justifiedbythedramaticallytoosmallstart- phonon coupling strength associated with the (electron- ing point DFT-LDA gap for small isolated molecules, doped) lowest unoccupied 3-fold degenerate t molecu- has been shown in the case of gas phase organic systems 1u lar orbital (LUMO) has been explored within DFT and to produce more accurate ionization potentials and elec- GW. The specific choice of the LUMO level is dictated tronicaffinities,34,35,48–51 togetherwithimprovedoptical by the physics of the phonon-mediated superconducting excitation energies at the GW/BSE level,36,52–54 when transition in the fullerides.42,43 We adopt the standard LDAorPBEKohn-Shameigenstatesarechosenasstart- definition of the effective electron-phonon coupling po- ing points. Consistently, when starting from DFT-LDA, tential Vep that enters e.g. the McMillan formula for EPC matrix elements in C were also found to be in 60 4 FIG.1: Schematicrepresentationofchangesintheelectronicstructureunderphonondistortion: a)inC ,the3-foldLUMO 60 levels (blue) split under a distortion (red) along a H phonon mode. b) In graphene, an energy gap ∆E at the Dirac point is g opened(red)whendistortingtheequilibriumstructure(blue)alongtheΓ−A(cid:48) orK−E phononmodes. c)Indiamond,the 1 2g lowestlyingconductionbandsatΓ/X (blue)arerenormalisedbytheircouplingtoazone-boundaryphononmode. Forclarity, only the states showing the largest shift are plotted (red dot, see text). closer agreement with experiment within evGW as com- tive (non self-consistent) GW approach to obtain quasi- pared to G W .9 particle energies. 0 0 For the periodic diamond and graphene systems, we For diamond, following Antonius and coworkers,12 we use the Yambo code which implements many-body per- haveusedarelaxedDFT-LDAgeometry(latticeparame- turbation theory within a plane wave formalism.55 The tera=3.591˚A)withaplanewavecutoffenergyof80Ry startingDFTKohn-Shameigenstatesaregeneratedwith and norm-conserving pseudopotentials. All calculations the Abinit code.56 In the case of graphene, the running were performed with a 8x4x4 uniform Monkhorst-Pack parameters are identical to those used in a previous GW grid to describe the supercell necessary for the frozen- study of the electron-phonon coupling in this material,5 phonon simulation of the (q = X) phonon. We have namely we use a = 2.46 ˚A as lattice parameter, 15 a.u. used 160 bands and a cutoff of 8 Ha for the construc- distance between the layers,57 36 36 1 k-points for tion of the dielectric function and the GW self-energy. × × the GW calculations in the Brillouin zone (BZ) corre- Following Refs. 12 and 58 where the diamond band-gap sponding to the primitive (2-atoms) cell and the nearest renormalization by the zero-point motion was studied, equivalent k-grid in the BZ corresponding to the super- we define the relevant EPC matrix elements as the sec- cell needed to describe zone-boundary (q=K) phonons ond derivative of the electronic eigenvalues with respect inthefrozen-phononapproach. Theenergycutoffonthe to the phonon displacements at equilibrium: plane wave basis was set to 60 Ry. All conduction bands within an energy range of 45 eV above the Fermi level wereusedtobuildGandW,andacutoffof4Haforthe ∂(cid:15) ¯h ∂2 mk = (cid:15) , (3) constructionofthedielectricfunction. TheGodby-Needs ∂n 2Mω ∂u2 mk νq νq νq plasmon pole model was employed for the dynamical dielectric constant entering in W. where n is the occupation of the phonon state (ν) of νq As schematically described in Fig. 1(b), we focus on wave vector q. The EPC matrix elements are now re- the splitting ∆E of the degenerate occupied and unoc- lated to the temperature dependent renormalization of cupied levels at the Dirac point caused by the coupling the electronic bands. We evaluate the second-derivative to the qν =Γ-E optical mode at zone center together 2g of the eigenvalues using a symmetric five-point finite- with the zone-boundary qν = K-A(cid:48) phonon, for which difference formula with displacements of 0.01, 0.02 and 1 the strongest renormalization of the coupling constant 0.03 ˚A. In Ref. 59, it was already shown that stan- have been observed upon replacing the DFT-LDA ap- dardDFTandDFPTmethodswith(semi)localfunction- proach by the G0W0 formalism.5,6 The EPC matrix ele- als strongly underestimate (by ca. 30%) the zero-point mentsofinterestarethushererelatedtothederivatives: renormalization of the direct gap of diamond compared (∂∆E/∂uqν), namely (see Ref. 5): to experiment, while Ref. 12 demonstrated that G0W0 calculationsrestoredthegoodagreementbetweentheory andexperiment. Wefocusherebelowontherenormaliza- 1 (cid:18)∂∆E(cid:19)2 D2 = , (2) tion of the lowest conduction band edge at (mk = X1c) (cid:104) qν(cid:105) Aqν ∂uνq and (mk = Γ15c) by the coupling to the (νq = X4) phonon modes at the zone-boundary (see Fig. 1c). It expressed in (eV/˚A)2, where Aqν equals 16 and 8 for is for such a phonon wave vector that the EPC matrix qν = Γ-E and qν = K-A(cid:48), respectively. We per- elements were found to be large and that also the GW 2g 1 formed“single-shot”G W andCOHSEXcalculationson correction was shown to be the most significant.12 Fol- 0 0 top of DFT-LDA results, namely the standard perturba- lowing Ref. 12, the provided EPC matrix elements are 5 LDA G W G W (W) COHSEX COHSEX(W) 0 0 0 0 Graphene Γ-E 44.3 68.1 (65.2) 72.0 81.5 83.2 2g (cid:48) K-A 88.6 201(187) 207 260 256 1 Diamond Γ 1.281 0.611 0.657 1.311 1.445 15c X -2.009 -1.031 -1.066 -2.124 -2.273 1c TABLE I: Calculated EPC-related quantities (see text) for graphene and diamond. In the case of graphene, these correspond to the splitting of the degenerate highest valence and lowest conduction band at the Dirac point under the influence of a (qν = Γ-E ) and (qν = K-A(cid:48)) phonon mode and are expressed in (eV/˚A)2. Number in parenthesis are obtained with a 2g 1 plasmonpoleenergyof7eV,closetothegrapheneπ-plasmonenergy,insteadofthe27eVoriginallyusedinRef.5(seetext). Fordiamond,thesearetheresultsofequation(3)fortheX andΓ conductionstatesandfortheX zone-boundaryphonon 1c 15c 4 mode only (values in eV). This corresponds to twice the zero-point renormalization arising from this given phonon mode (see Ref. 12). averaged over the degenerate electronic manifold consid- tronic state cannot be made. Our findings are, however, ered (namely the 3 degenerate levels at zone-center and very consistent with what is observed for the Γ state, 15c the 2 degenerate levels at the X zone boundary). namelyadramaticdecreaseofthecouplingstrengthwith 1c the X phonon mode which, again, is the dominant cou- 4 pling mode. The coupling energy is indeed found to be IV. RESULTS reduced by 52% and 49% from LDA to GW for the Γ 15c andX states,respectively. SincetheX stateiscloser 1c 1c A. Reference GW calculations to the true conduction band minimum, this is a strong indication that zero-point motion renormalization of the indirect gap will certainly also be strongly affected by We first reproduce for validation and for reference the the GW correction. The full study of such an effect is previously published GW calculations related to EPC beyond the scope of the present paper. matrix elements in graphene, diamond and C . Our re- 60 For C , the total evGW coupling potential is within sults are compiled in Table I for diamond and graphene 60 8% of that found by Ref. 9 as a result of the larger TZP andinTableIIforC ,respectively. OurLDAandG W 60 0 0 basis used here.38 Our evGW Vep potential is found to values for graphene are close to those found in Ref. 5, namely e.g. 201 eV/˚A (G W value, present study) for increase by 44% with respect to the corresponding LDA 0 0 the largest matrix element with the (K-A(cid:48)) phonon, to value, to be compared with the 48% increase obtained 1 in Ref. 9 with a smaller basis. This indicates the good be compared to 193 eV/˚A in the previous study.57 The convergenceoftheGW correctiontotheLDAvalue. We difference can be explained by an increased five points, note that while DFT and GW calculations are known instead of only two, finite-difference formula, and by the toconvergeverydifferentlywithrespecttobasissize, we increaseofthedielectricmatrixsizeenergycutoff,from2 compareinthepresentstudyapproacheswhicharemuch to4Hacutoff. WealsoverifiedthatchangingtheGodby- moresimilar,namely“standardandapproximated”GW Needs plasmon model input finite frequency, from the calculations, indicating certainly an even better conver- (default) 27 eV value in Ref. 5 to a 7 eV value closer to gence of the differences we are interested in. the π-plasmon resonance in graphene, does not change significantly the calculated EPC strength (see numbers in parenthesis in the Table I). Such differences are neg- ligible with respect to the more than 100% increase as B. The COHSEX approximation compared to the LDA value. For diamond, the EPC contribution associated with We now explore approximations to the full GW cal- the mk = Γ conduction states are in good agree- culations performed in Refs. 5,9,12 and reproduced here 15c ment with the LDA and G W results of Antonius and above. Westartbythecaseofdiamond. Ascomparedto 0 0 coworkers,12 i.e. we find a 11 meV difference with re- the G W calculations, the static COHSEX approxima- 0 0 spect to their EPC in the case of LDA, while we find a tion is shown to induce errors larger than 100% for the difference of 8 meV for the G W results.68 These small relevantelectron-phononcouplingenergytotheΓ and 0 0 15c variations can be ascribed to the different pseudopoten- X states. In both cases, the static COHSEX approach 1c tials (Ref. 59 showed that pseudopotentials can lead to strongly overestimates the electron-phonon coupling. A errors up to 50 meV), and to the different convergence representation of the errors as compared to the G W 0 0 parameters. calculations is provided in Fig. 2. SinceRef.12focusedontherenormalizationofthedi- The case of graphene is of specific interest since, in rect band gap at the zone-center, comparison with the greatcontrasttodiamondandC ,thefundamentalgap 60 present results for the phonon coupling to the X elec- reduces to zero in the equilibrium geometry and opens 1c 6 Theory Experiment Mode LDA evGW COHSEX evCOHSEX evCOHSEX(W) Fullerene H (1) 4.7 5.85 4.99 6.3 6.3 g H (2) 9.3 10.1 9.99 10.3 11.9 g H (3) 8.8 12.95 11.96 13.7 15.3 g H (4) 4.2 5.5 4.9 5.4 5.7 g H (5) 3.98 4.9 4.5 5.3 5.5 g H (6) 1.8 1.99 2.05 2.2 2.4 g H (7) 15.8 25.4 23.3 27.9 26.5 g H (8) 13.1 18.3 17.1 19.8 19.5 g A (1) 1.2 1.9 1.85 1.8 1.2 g A (2) 7.2 13.98 13.2 15.4 12.2 g Total A 8.4 15.9 15.0 17.2 13.4 g Total H 61.7 84.96 78.8 90.9 93.1 96.2b,96.5c g Total 70.2 100.9 (108.6a) 93.8 108.1 106.6 106.8b TABLEII:Mode-resolvedVep couplingpotentialassociatedwiththe3-folddegenerateC LUMOandthecorrespondingVep ν 60 totals (in meV). For group symmetry reasons, only the H and A vibrational modes couple to this electronic state. g g a Ref. 9, b Ref. 46, Table V, c Ref. 47. Diamond Graphene 350 300 Γ15cstate W) X1cstate W) 40 Γ(E2g)mode K(A1’)mode %err112205050000 LDA W) COHSEX COHSEX( LDA W) COHSEX COHSEX( %err−22000 LDA GW00 W(W)0 HSEX EX(W) LDA GW00 (W) HSEX EX(W) 50 GW00 GW(00 GW00 GW(00 −40 G0 CO COHS GW00 CO COHS 0 60 − FIG. 2: Relative error of the different approximations with FIG. 3: Relative error of the different approximations with respect to the G0W0 results for the “thermal renormaliza- respecttotheG0W0 resultsfortheelectron-phononcoupling tion” EPC energies (Eqn. 3) for the lowest conduction band in graphene. of diamond. is observed again to systematically increase the coupling upon phonon distortions. For the optical phonon mode strength as compared to evGW.62 This result is consis- at the zone-center, the COHSEX approximation leads tent with what was observed for diamond and graphene, again to an increase as compared to the G W value, 0 0 even though clearly the error appears to vary from one eventhoughsmallerthaninthediamondcasewitha20% system to another. error. Such an error becomes larger for the K-A(cid:48) mode 1 witha29%increaseoftheEPC-relatedgap-openingrate. In the present graphene case, the error induced by the COHSEX approach remains smaller than the 34% and C. The constant screening approximation 56% reduction observed upon using DFT-LDA instead of G W . However, the close to 30% error observed for Following the above-mentioned analogy with the 0 0 the COHSEX calculation of D2 is clearly significant Bethe-Salpeter formalism, we finally test another ap- (cid:104) K(cid:105) (see Fig. 3), questioning again the applicability of the proximation, that is the constant screening approach. static COHSEX formalism for EPC matrix elements cal- Namely, we explicitly calculate the screened Coulomb culations. potential for the undistorted structure and upon chang- We now consider the C case, where we compare ing the positions of the ions along the (R ) vibrational 60 qν evGW and evCOHSEX calculations. A global increase eigenmodes, we keep W frozen to the undistorted value. of about 7% from evGW to its static approximation is Formally,thisamountstoassumingthat: ∂GW/∂uqν observed. Besides the H (4) and A (1) modes showing (∂G/∂uqν)W. SuchanapproachwillbelabeledGW(W≈) g g verysmallcouplings,thestaticCOHSEXapproximation or COHSEX(W) in the following. 7 Fullerene Comparingthetotalcoupling,evCOHSEX(W)andev- 30 COHSEX agree within 1.5 %. 20 Hgmode We now explore the constant-screening approximation 10 for periodic diamond. Results are again compiled in Ta- 0 ble I. For diamond, we observe a 3.4% and 7.5% dif- %err−−2100 evGW COHSEX HSEX(W) cEXP dEXP fcΓeo1rn5ecsntcasetnatbtseecstrw,eeereennsinpGgec0atWipvpe0rl(yoW.xi)mAaatntidotnhGey0iWCelOd0sHfeoSrrErotXrhseloeXfvet1hlc,eatonhrde- 30 v O der of 7.0% and 10.2%, respectively. In summary, while − A e C −40 LD ev the COHSEX approximation was failing for diamond, the constant-W approximation leads again to very rea- 50 − sonnable results, with errors well within 10% at the GW level. As illustrated in Fig. 2, this error is much smaller FIG. 4: Relative error of the different approximations with respect to the evGW results for the H -related Vep of equa- than the error induced by the LDA approximation. g tion (1). For convenience, a comparison to the experimental We finally address the case of graphene. Concern- data is provided. ing the coupling with the zone-center optical mode, the constant-W approach leads to a 5.7% and 2.1% error respectively when comparing G W (W) to G W and 0 0 0 0 Inthecaseofthefullerene,wherecalculationsareper- COHSEX(W) to COHSEX. In the case of coupling with formed within a Gaussian basis, special care must be the zone-corner K-A(cid:48) phonon mode, the discrepancy is 1 takenintheimplementationofsuchaconstant-screening of 3% and 1.5% applying the constant-screening approx- approximation. Intheauxiliarybasis,orresolutionofthe imation to GW and COHSEX approaches, respectively. identityapproach,thebareandscreenedCoulombpoten- Notice that the results for the constant-screening of the tialsareexpressedintermsofanatom-centeredauxiliary K-A(cid:48) phonon are affected by a numerical error of about 1 basis with the following relations: 4%, therefore the discrepancy could be slightly larger. Anywaythisdifferenceismuchsmallerthantheerrorin- (cid:90) (cid:90) duced by the standard DFT-LDA calculations, but also [W]β,β(cid:48) = drdr(cid:48)β(r)W(r,r(cid:48))β(cid:48)(r(cid:48)), (4) much smaller than the error induced by the static COH- SEX formalism as compared to the reference GW calcu- W(r,r(cid:48))=(cid:88)β(r)(cid:0)S−1[W]S−1(cid:1) β(cid:48)(r(cid:48)), (5) lation. β,β(cid:48) β,β(cid:48) where S is the overlap matrix in the auxiliary basis. Us- V. DISCUSSION ingnowthenotation: W andβ forthescreenedCoulomb potential and the auxiliary basis for the slightly dis- torted system, the assumption: W(r,r(cid:48)) W(r,r(cid:48)) Clearly, among the two approximations explored here leads straightforwardly to the condition: (cid:39) above, the constant-screening stands as a much better approach than the static COHSEX for the calculation of electron-phonon coupling matrix elements within the [W] S [W] S , (6) GW formalism. Overall, the largest error induced by β,β(cid:48) ββ β,β(cid:48) β(cid:48)β(cid:48) (cid:39) the static COHSEX approximation is larger than 100% where S =< β β > is an overlap matrix between the in the diamond case, while it is reduced to 7% in the ββ | auxiliary bases for the perturbed and unperturbed sys- caseoftheconstant-screeningapproximationtofullGW tems, respectively. calculations. For C again, we test the constant screening approx- The static COHSEX approximation is known to al- 60 imation at the COHSEX level only, namely comparing ways overestimate band gaps. As proposed in Ref. 8, a evCOHSEXand evCOHSEX(W)calculationsof theVep too large band gap should lead to an underscreening of energies. While the plasmon-pole approach used for dia- the electron-phonon interaction and consequently to en- mond and graphene is based on a fixed finite frequency hancedEPCmatrixelements. Thisisconsistentwiththe pole value independent of the geometry, the real axis fact that, within Hartree-Fock or with increasing exact poles contribution to the correlation energy in the con- exchangeinhybridfunctionals,theEPCmatrixelements tourdeformationapproach(seeRefs.19,34)implemented are found to steeply increase. The Hartree-Fock D2 intheFiestapackagechangesfromonestructuretoan- and D2 coupling constants in graphene were foun(cid:104)dΓto(cid:105) (cid:104) K(cid:105) other,leadingtodifficultieswhentryingtoimplementthe be about 5 and 30 times larger, respectively, as com- constant-W approach within the frozen-phonon scheme. pared to their GW analog,5 while in C increasing the 60 The results of the constant-screening approximation percentage of exact exchange from 20% to 30% in hy- are compiled in Tables I and II. For C , an excel- brid functionals was found to enhance the Vep energy 60 lentagreementisobtainedwithinevCOHSEX(W)com- by about 15%.8 Such an interpretation matches the ob- pared to the corresponding evCOHSEX calculations. servation that in graphene and C , the COHSEX EPC 60 8 coupling constants are larger than their corresponding VI. CONCLUSIONS GW reference. In the case of diamond, where we are concernedwiththesecond-orderderivatives,theeffectof assuming a static approximation may be more difficult We have explored two approximations for calculating to interpret. self-energy gradients of interest for electron-phonon cou- AnimportantobservationfurtheristhattheCoulomb- pling, namely the two simplifications commonly used in hole term is a local potential in the static COHSEX ap- the GW/Bethe-Salpeter calculations, that is the static proximation, washing out spatial local-field effects, a ob- COHSEX and the constant screening approximations. servation that may potentially allow to understand the We explored these approaches in the case of diamond, failure of the COHSEX approximation to reproduce ac- graphene and C . Our findings suggest that the COH- 60 curatelytheevolutionofthebandstructurewithphonon SEXapproximationcannotbetrustedtoimproveonthe deformation. In any case, the large variations of the in- DFT-LDA values as clearly illustrated in particular in duced errors from diamond to graphene and C60 is still the case of diamond. On the contrary, the constant to be understood, and we will just stand here on the ob- screening hypothesis, namely assuming that W remains servationthatthestaticCOHSEXapproximationcannot to first order constant with respect to small ionic dis- be trusted to improve on DFT-LDA calculations. placements, seems to be reliable, with a discrepancy no Concerning the neglect of the gradient of W with larger than 10%, even in the difficult case of graphene respect to the ionic positions, such an approximation where the phonon perturbation dramatically affects the was already tested by Ismail-Beigi and Louie in the Diracconeandthe(semi)metallicnatureofthegraphene context of excited state ionic forces within the Bethe- sheet. Even though a deeper understanding of the ori- Salpeter formalism,63 showing good agreement with ex- gin of the specific difficulties uncountered by the static plicit finite-difference ”exact” BSE calculations for small COHSEX approximation would allow to better rational- CO and NH3 molecules. It is commonly assumed that ize the validity and limits of the tested approximations, the GW/BSE approach is much more resistant to ap- the present results offer promising perspectives to carry proximations on W as compared to the GW approach on such many-body evaluations of the electron-phonon for (charged) excitations. This is due to cancellations coupling gradients with much reduced computer cost on of errors between the electron-electron and electron-hole realisticsystems,includingthestudyofperiodicsystems interactions. Namely, any error introduced in W is ex- witharbitrarywavevectorperturbation. Furtherstudies pected to affect excitonic interactions and quasiparticle arehoweverrequiredonalargersetofsystemsandphysi- gapsinoppositeways. Clearly, thepresentGW studyof calobservablesinordertobetterassesstheinterest,with thevariationsofagivenquasiparticleenergywithrespect respect to common DFT calculations, of using the GW toionicpositionscannotbenefitfromsuchcancellationof approach,anditsvariousapproximations,forcalculating errors. Still, the constant-screening approximation turns self-energy gradients with respect to ionic positions. to be a reliable approach to save on computational cost. An important consequence of the present findings is Acknowledgments. P.B. warmly acknowledges G. that once the screened Coulomb potential W(r,r(cid:48);ω) is Antonius for providing the diamond phonon eigenmodes built for the equilibrium geometry, the calculation of used in this study and for guiding our comparison with the variations of the quasiparticle energies with respect Ref.12. C.A.acknowledgesDanieleVarsanoforthecor- to the perturbation (λ) only requires the evaluation of rectionofsomebugsintheYambocode. C.F.isindebted the variations of the Green’s function G with respect to to the French CNRS and CEA for Ph.D funding. P.B. the perturbation. This can be performed within stan- acknowledges a postdoctoral fellowship from the French dard DFPT techniques, at least in the case of non-self- National Research Agency under Contract No. ANR- consistentG W calculationswheretheGreen’sfunction 2012-BS04PANELS.Computingtimehasbeenprovided 0 0 assumes an explicit form in function of the input DFT by the “Curie” national GENCI-IDRIS supercomputing eigenstates. This may invite, for systems such as C , to center under contract No. i2012096655 and a PRACE 60 useG W calculationsstartingfromDFTeigenstatesob- European project under contract No. 2012071258. X.B. 0 0 tainedwithhybridfunctionals,whichhavebeenshownto acknowledges illuminating discussions with Steven G. be a better starting point for finite-size systems.54,64–66 Louie and G. Antonius. 1 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864-B871 5 M. Lazzeri, C. Attaccalite, L. Wirtz, and F. Mauri, Phys. (1964). Rev. B 78, 081406 (2008). 2 S.Baroni,S.deGironcoli,A.DalCorso,andP.Giannozzi, 6 D. Basko and I. L. Aleiner, Phys. Rev. B, 77, 041409, Rev. Mod. Phys. 73, 515 (2001). (2008). 3 X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997). 7 A. Gru¨neis, J. Serrano, A. Bosak, M. Lazzeri, S. L. 4 J.P.Perdew,K.BurkeandM.Ernzerhof,Phys.Rev.Lett. Molodtsov, L. Wirtz, C. Attaccalite, M. Krisch, A. Ru- 77, 3865-3868 (1996). bio, F. Mauri, and T. Pichler, Phys. Rev. B 80, 085423 9 (2009). descriptionoftheeigenstates.Thisbasisislargerthanthe 8 J. Laflamme Janssen, M. Coˆté, S.G. Louie, and M.L. Co- double-zeta plus polarization (DZP) basis used in Ref. 9, hen, Phys. Rev. B 81, 073106 (2010). withtheinclusioninparticularofdiffuses andp orbitals. 9 C. Faber, J. Laflamme Janssen, M. Coˆté, E. Runge, and 39 N. Troullier and J.L. Martins, Phys. Rev. B 43, 1993 X. Blase Phys. Rev. B 84, 155104 (2011). (1991). 10 Z.P. Yin, A. Kutepov, and G. Kotliar, Phys. Rev. X 3, 40 T. Anh Pham, H.-V. Nguyen, D. Rocca, and G. Galli, 021011 (2013). Phys. Rev. B 87, 155148 (2013). 11 S.Ciuchi,R.C.Hatch,H.Ho¨chst,C.Faber,X.Blase,and 41 Onehastonotehoweverthattheimpactofthefunctional S. Fratini, Phys. Rev. Lett. 108, 256401 (2012). on the vibrational eigenmodes used to calculate the EPC 12 G. Antonius, S. Poncé, P. Boulanger, M. Coˆté, and X. matrixelementswasfoundtobemarginal,themaineffect Gonze, Phys. Rev. Lett. 112, 215501 (2014). comingfromtheframeworkwithinwhichthequasiparticle 13 P.C. Martin and J. Schwinger, Phys. Rev. 115, 1342 energies are calculated. (1959). 42 A.F. Hebard et al., Nature (London) 350, 600 (1991). 14 L. Hedin, Phys. Rev. 139, A796 (1965). 43 O. Gunnarsson, Rev. Mod. Phys. 69, 575 (1997). 15 M.S.Hybertsen,S.G.Louie,Phys.Rev.B34,5390(1986). 44 V.P. Antropov, O. Gunnarsson, A.I. Liechtenstein, Phys. 16 R.W. Godby, M. Schlu¨ter, and L.J. Sham, Phys. Rev. B Rev. B 48, 7651 (1993). 37, 10159 (1988). 45 Coˆté,M.andGrossman,J.C.andCohen,M.L.andLouie, 17 G.Onida, L.Reining, A.Rubio, Rev.Mod.Phys.74, 601 S. G., Phys. Rev. Lett. 81, 697 (1998). (2002). 46 N.Iwahara,T.Sato,K.Tanaka,andL.F.Chibotaru,Phys. 18 F.Aryasetiawan,O.Gunnarsson,Rep.Prog.Phys.61,237 Rev. B 82, 245409 (2010). For an interpretation of the (1998). experimental data, see Note 21 of Ref. 9. 19 B. Farid, R. Daling, D. Lenstra, and W. van Haeringen, 47 I.D. Hands, J.L. Dunn, L. Janette, C.A. Bates, A. Colin, Phys. Rev. B 38, 7530 (1988); B. Farid, in Electron Cor- M.J. Hope, J. Michael, S.R. Meech, D.L. Andrews, Phys. relation in the Solid State, edited by N. H. March (World Rev. B 77, 115445 (2008). For an interpretation of the Scientific,Singapore,1999),p.217,andreferencestherein. experimental data, see Note 22 of Ref. 9. 20 F.Bruneval,Ph.D.thesis,EcolePolytechnique,Palaiseau, 48 P.H.Hahn,W.G.Schmidt,andF.Bechstedt,Phys.Rev. France, 2005. B72,245425(2005).Inthisstudy,theimportanceofcor- 21 F.Bruneval,N.Vast,L.Reining,Phys.Rev.B74,045102 rectingthestartingeigenvaluesintheconstructionofWis (2006). emphasized for molecular systems, bringing some light in 22 M.Gatti,F.Bruneval,V.Olevano,L.Reining,Phys.Rev. thediscussionsconcerningtheupdateoftheeigenvaluesin Lett. 99, 266402 (2007). theGreen’sfunctiononly(theGW scheme)orinbothG 0 23 F.Bruneval,N.Vast,L.Reining,M.Izquierdo,F.Sirotti, and W. N. Barrett, Phys. Rev. Lett. 97, 267601 (2006). 49 L.Tiago,P.R.C.Kent,R.Q.Hood,andF.A.Reboredo, 24 F. Trani, J. Vidal, S. Botti, and M.A.L. Marques Phys. J. Chem. Phys. 129, 084311 (2008). Rev. B 82, 085115 (2010). 50 Sahar Sharifzadeh, Ariel Biller, Leeor Kronik, and Jeffrey 25 J. Vidal, S. Botti, P. Olsson, J.-F. Guillemoles, and L. B. Neaton, Phys. Rev. B 85, 125307 (2012). Reining, Phys. Rev. Lett. 104, 056401 (2010). 51 Conor Hogan, Maurizia Palummo, Johannes Gierschner, 26 M.J. Lucero, I. Aguilera, C. Diaconu, P. Palacios, P. and Angel Rubio, J. Chem. Phys. 138, 024312 (2013). Wahnoń, and G. E. Scuseria, Phys. Rev. B 83, 205128 52 C. Faber, P. Boulanger, I. Duchemin, C. Attaccalite, and (2011). X. Blase, J. Chem. Phys. 139, 194308 (2013). 27 W. Kang, M.S. Hybertsen, Phys. Rev. B 82, 195108 53 B. Baumeier, D. Andrienko, Y. Ma, and M. Rohlfing, J. (2010). Chem. Theory Comput. 8, 997 (2012). 28 L.J. Sham and T.M. Rice, Phys. Rev. 144, 708 (1966). 54 P. Boulanger, D. Jacquemin, I. Duchemin, X. Blase, J. 29 W.HankeandL.J.Sham,Phys.Rev.Lett.43,387(1979). Chem. Theory Comput. 10, 1212 (2014). The effect of 30 G. Strinati, La Rivista del Nuovo Cimento 11 (12), 1-86 partial self-consistency on the eigenvalues, with various (1988). starting-point wave functions, is discussed, in comparison 31 M. Rohlfing and S.G. Louie, Phys. Rev. Lett. 80, 3320 with high-level coupled-cluster data. (1998). 55 A.Marini,C.Hogan,M.Gruu¨ning,andD.Varsano,Com- 32 L.X. Benedict, E.L. Shirley, and R.B. Bohn, Phys. Rev. put. Phys. Commun. 180, 1392 (2009). Lett. 80, 4514 (1998). 56 X. Gonze et al., Comp. Mat. Science, 25, 478 (2002). 33 S.Albrecht,L.Reining,R.DelSole,G.Onida,Phys.Rev. 57 Notice that in Ref. 5 the interlayer distance between the Lett. 80, 4510 (1998). graphene sheets is 15 a.u. and not 20 a.u. as it is wrongly 34 X.Blase,C.Attaccalite,V.Olevano,Phys.Rev.B83,536 reported in the paper. 115103 (2011). 58 P.B.AllenandV.Heine,J.Phys.C9,2305(1976);P.B. 35 C. Faber, C. Attaccalite, V. Olevano, E. Runge, X. Blase, Allen and M. Cardona, Phys. Rev. B 23, 1495 (1981). Phys. Rev. B 83, 115123 (2011). 59 S. Poncé, G. Antonius, P. Boulanger, et al; Comput. 36 X.BlaseandC.Attaccalite,Appl.Phys.Lett.99, 171909 Mater. Sci. 83, 341 (2014). (2011). 60 For sake of indication, at the LDA level, the two lowest 37 J.M.Soler,E.Artacho,J.D.Gale,A.Garc´ıa,J.Junquera, X and X vibrational modes lead to a (negative) -0.14 1 2 P.Ordejoń,D.Sańchez-Portal,J.Phys.: Condens.Matter eV coupling energy, one order of magnitude smaller than 14, 2745 (2002). the X mode coupling strength. 4 38 Thistriple-zetapluspolarization(TZP)basisconsistsof3 61 P. Trevisanutto, C. Giorgetti, L. Reining, M. Ladisa and shells of s and p orbitals, complemented by one polariza- V. Olevano, Phys. Rev. Lett. 101, 226405, (2008) tion d shell, that is 19 orbitals per carbon atom, for the 62 We provide for information in Table II the results asso- 10 ciated with the (single-shot) COHSEX approach, starting Marques, S. Botti, J. Chem. Theory Comput. 10, 3934 fromDFT-LDAeigenstates.Clearly,theEPCcouplingen- (2014). ergiesapproachfrombelowtheevGW valuewithanerror 67 TheexperimentalbandgapistakenfromtheNISTchem- notlargerthanthatinducedbytheevCOHSEXtechnique. istry webbook at [http://webbook.nist.gov/chemistry/]. ThisstemsfromacancellationoferrorsbetweentheCOH- Thebandgapiscalculatedasthedifferenceoftheexperi- SEXapproximationandthelackofself-consistencywhich mental ionization energy and electronic affinity. Weadopt leads to overscreening when building W directly from the the most recent experimental data. DFT-LDA spectrum. 68 DirectcomparisonwithFigure2ofRef.12isnotpossible, 63 S. Ismail-Beigi and S. G. Louie, Phys. Rev. Lett. 90, sincethefigurepresentsthetotalelectron-phononcoupling 076401 (2003). atwavevectorq=X,i.e.allphononmodesareconsidered 64 Thomas Ko¨rzdo¨rfer and Noa Marom, Phys. Rev. B 86, andnotjustX asinthisstudy.Thecomparisonwasmade 4 041110(R) (2012). directly with the data used in Ref. 12 via personal corre- 65 F. Bruneval, M.A.L. Marques, J. Chem. Theory Comput. spondence. 9, 324 (2013). 66 S. Korbel, P. Boulanger, I. Duchemin, X. Blase, M.A.L.