Full self-consistency versus quasiparticle self-consistency in diagrammatic approaches: Exactly solvable two-site Hubbard Model
aa r X i v : . [ c ond - m a t . s t r- e l ] M a r Full self-consistency versus quasiparticle self-consistency in diagrammatic approaches:Exactly solvable two-site Hubbard Model
A.L. Kutepov
Ames Laboratory USDOE, Ames, IA 50011
Self-consistent solutions of Hedin’s diagrammatic theory equations (HE) for the two-site Hub-bard Model (HM) have been studied. They have been found for three-point vertices of increasingcomplexity (Γ = 1 (GW approximation), Γ from the first order perturbation theory, and exactvertex Γ E ). The comparison is being made when an additional quasiparticle (QP) approximationfor the Green function is applied during the self-consistent iterative solving of HE and when QPapproximation is not applied. The results obtained with the exact vertex are directly related to thepresently open question - which approximation is more advantageous for future implementations -GW+DMFT or QPGW+DMFT. It is shown that in the regime of strong correlations only orig-inally proposed GW+DMFT scheme is able to provide reliable results. Vertex corrections basedon Perturbation Theory (PT) systematically improve the GW results when full self-consistency isapplied. The application of the QP self-consistency combined with PT vertex corrections showssimilar problems to the case when the exact vertex is applied combined with QP sc. The analysisof Ward Identity violation is performed for all studied in this work approximations and its relationto the general accuracy of the schemes used is provided. PACS numbers: 71.10.Fd, 71.27.+a
I. INTRODUCTION
One of the challenges for the computational theoristsworking in the solid state electronic structure field isthe robust implementation of the so called GW+DMFTmethod (combination of GW approximation (G- Green’sfunction, W - screened interaction) and DynamicalMean-Field Theory). The scheme was originally pro-posed by Sun and Kotliar and in slightly different (butprobably more commonly known) form by Biermann etal. The basic idea of the approach is to separate all ac-tive space of the basis set into ”weakly correlated” partfor which GW approximation is supposed to work welland ”strongly correlated” part for which one sums upmany diagrams (to infinite order if one uses the non-perturbative DMFT solver as for example the Quan-tum Monte Carlo - QMC). In order for the method tobe well defined, everything should be done till full self-consistency (sc), including the iterations of the GW partitself and also the ”internal” iterations of DMFT part toensure that the solution of the impurity problem repro-duces the same G and W in strongly correlated subspaceas the ones in the same subspace projected from the GWpart.Despite its obvious appeal GW+DMFT has made onlyslow progress during more than decade since its first ap-pearance in the above mentioned papers. The reason wasnot just because scGW is quite demanding computation-ally but mostly because one has to satisfy the impurity sccondition not only for G (as in the LDA+DMFT method- a combination of Local Density Approximation in Den-sity Functional Theory and DMFT) but also for the W,which seems to be not an easy task. To the best of myknowledge there were only ”one-shot” type calculationsfor real materials where GW iterations were neglectedaltogether and DMFT self-consistency was imposed only on G, whereas W was fixed at LDA level and correspond-ingly the U was considered as an external parameter(calculated in constrained random phase approximation- cRPA). Such an implementation, of course, is a definitestep towards the full GW+DMFT scheme, but still onecannot say that there was the summation of all ”corre-lated” graphs in it which would require that the W in GWpart and W imp in the impurity (DMFT) part were thesame. Instead they were totally decoupled which makesit unclear of what kind of diagrams from DMFT partwere actually added to GW part. Nevertheless, togetherwith development of QMC solvers capable to handle dy-namical interactions, the hope is growing that the fullGW+DMFT scheme will eventually be implemented.There are some subtleties about the GW part as well.In its full sc implementation the method is very time con-suming which in part has prevented its applications forthe real materials. But recently a very efficient imple-mentation of it was published where the most computa-tionally demanding parts (calculation of the polarizabil-ity P and the self energy Σ) are performed in real spaceand Matsubara’s time. As a result, it became possibleto successfully apply scGW to the actinides Pu and Amand (earlier) to simple sp-materials. However for the ma-jority of solids, scGW produces worse spectra than fast”one-shot” GW and this is the other reason why scGWis not popular. The reason for the failure of scGW withspectra may be traced as an extremely non symmetri-cal ”dressing” of the ”initial” Green’s function with selfenergy insertions of GW-only form in the course of self-consistency iterations and neglecting by vertex correc-tions.The origin of the problem with the scGW methodcan also be formulated in terms of the absence of theZ-factor cancelation, which again happens because weneglect by the vertex corrections. To resolve this prob-lem Kotani and Schilfgaarde devised a beautiful trickof doing yet another approximation. They used quasi-particle (QP) part of Green’s function G QP (instead offull Green’s function G calculated from Dyson’s equation(DE)) to calculate P and Σ on every iteration till self-consistency. The trick is, that the errors from the abovetwo approximations (using G QP instead of full G and ne-glecting by vertex corrections) mostly cancel each otherout and as a result the QPscGW (self-consistent quasi-particle GW) method usually gives much better spectrathan full scGW. The important fact is that QPscGW notjust slightly improves the one-shot GW description of sp-materials (which are good enough already in the one-shotGW), but often gives reasonable results for the materialswith d- or f-electrons too, and the method doesn’t relyon a particular starting point. It is totally self-consistent.Quite naturally, the success of QPscGW with spec-tral properties (as compared to the full scGW) hasignited the ideas to formulate another approach -QPscGW+DMFT , where one supposedly addsDMFT corrections to P and Σ as in GW+DMFT butuses QPscGW for the ”big” iterations in the ”weaklycorrelated” part.As it appears we have now two schemes proposed:GW+DMFT and QPGW+DMFT. In this work I amdoing an attempt to ”estimate” what to expect fromthe future implementation of the schemes. The analy-sis strongly depends on the ”separability” of weakly andstrongly correlated parts. I will assume here that they areperfectly separable and the correlations in the GW partare really weak. In this case, Z-factor in the GW partis close to 1, so that this part is equally well describedin both the GW and the QPGW approximations. Thedifference correspondingly comes only from the DMFTpart. For this part I assume that we are able to solveit exactly which means in particular that both sc condi-tions (for G and W) are satisfied. For GW+DMFT, itmeans G = G imp , W = W imp , and for QPGW+DMFT,it means G QP = G imp , W QP = W imp , where followingthe arguments of work [8] we have G QP ≈ GZ . The exactsolution of the impurity problem also means that corre-sponding self energies can be written in their exact dia-grammatic forms: Σ = − G Γ W and Σ QP = − G QP Γ W QP with Γ being the three-point vertex function which is in-cluded exactly. But following again the arguments inRef.[8], it becomes clear that in QPGW+DMFT case, i.e.when we include exact vertex and continue to apply QPapproximation for G, we will have a problem, because thefactor 1 /Z appears from the vertex and it doesn’t can-cel with the Z-factor in G as it happens in GW+DMFT.Basically it means the violation of Ward Identity (WI)in the QPGW+DMFT scheme. Thus, from this pointof view, QPGW+DMFT is highly problematic and itsdifficulties should grow up when the correlation strengthgrows up, because the 1 /Z -factor increases.The above simple argument against using theQPGW+DMFT scheme for strongly correlated materi-als needs some numerical support and I will provide it here using the exactly solvable two-site Hubbard model.For this model I calculate exactly the three-point ver-tex function and use it to calculate self-consistently theGreen functions G and G QP following two slightly differ-ent sc schemes (the scheme on the left side is basicallythe Hedin’s sc equations but with the 3-point vertexprecalculated exactly) P = G Γ Exact
G P QP = G QP Γ Exact G QP W = U + U P W W QP = U + U P QP W QP Σ = − G Γ Exact W Σ QP = − G QP Γ Exact W QP G = G + G Σ G G QP = Z G − G Σ QP , (1)where U is the bare interaction in the Hubbard modeland G if Green’s function in Hartree approximation. In(1) for QP case I formally represented the quasiparticleapproximation for G by simply dividing the full Greenfunction from Dyson’s equation by the Z-factor. Thisis for brevity. In fact I use the algorithm described inRef.[6] to construct G QP .It is obvious that at self-consistency the scheme on theleft side (I will call it G Γ E W ) is equivalent to the ex-act solution of GW+DMFT equations whereas the righthand scheme (I will call it QP G Γ E W ) is equivalent tothe exact solution of the QPGW+DMFT equations. Itis also obvious that the G Γ E W scheme is exact by con-struction. I use it in this work to check the numericalaccuracy of 3-point vertex evaluation. The QP G Γ E W scheme is approximate and below I will explore its accu-racy in different regimes of correlation strength for the 2-site Hubbard model. Also I will directly relate the prob-lems of the QP G Γ E W scheme with the degree of the WIviolation.Another goal of the present work was to explore thepossibility of combining the GW and the QPGW meth-ods with perturbative calculation of the 3-point vertexfunction. To this end I will use again the schemes similarto Eq.(1) but with Γ expanded to the first order in W (Γ )instead of exact Γ E . I also will show how the two corre-sponding perturbation theory based schemes ( G Γ W and QP G Γ W ) behave in the different regimes of parametersof the Hubbard model.This paper begins with the formal presentation of thetwo-site Hubbard model and the formulae used in thecalculations (Section II). In Section III the results arepresented and discussed. Finally in Section IV the con-clusions are given and the future plans are outlined. II. TWO-SITE HUBBARD MODEL
The Hamiltonian of two-site Hubbard model as it isused in this work is the following: H = − t X i = j,σ c + iσ c jσ + U X i X σσ ′ c + iσ c + iσ ′ c iσ ′ c iσ , (2)where t and U are the standard parameters of the Hub-bard model, c and c + are the destruction and creationoperators correspondingly, indexes i and j belong to thesites (1 or 2), and σ, σ ′ are the spin indexes.In this section, all equations which are used in the anal-ysis of the 2-site Hubbard model are collected for refer-ences. First I provide the energies and eigen vectors ofthe exact many-body states of the model for different oc-cupancies. Then, the exact expressions for Green func-tion and density-density correlation function are given.From them the exact self energy, polarizability, dielec-tric function, and screened interaction can be calculatedusing standard formulae. Next, two subsections providethe formulae for the exact 3-point vertex function andhow it is used to evaluate the corrections to the polar-izability and to the self energy in sc scheme (1). Per-turbation theory based equations for the 3-point vertexfunction are given next. Then I provide full and reduced(long-wave and long-wave+static limits) expressions forthe Ward Identity which are used later in this paper.Finally the formulae relevant to the evaluation of cur-rent three-point vertex function are given. For brevity,the derivations of the formulae are not provided, or onlysketch of derivation is given. In this work, the finite tem-perature framework is used so that in all equations be-low the time-frequency arguments are Matsubara’s time( τ ) and Matsubara’s frequencies ( ω for fermion frequencyand ν for boson frequency). A. Many-body states
In order to represent the many-body states of themodel it is convenient to introduce basis vectors | abcd i where all entries are equal to 0 or 1 in accordance withthe occupancies of corresponding one-electron states. Inthis work first two one-electron occupancies (a and b)correspond to spin-up and spin-down one-electron statesof the first site, and the third and fourth (c and d) corre-spond to the second site. Many-body energies and statesbelow have two indexes: the upper one corresponds tothe full occupancy (0,1,2,3, or 4) of the system, and thelower one distinguishes the states within the same occu-pancy.For the full occupancy N equal zero, we have corre-spondingly: E = 0; Ψ = | i , (3)for N=1 E = − t ; Ψ = √ ( | i + | i ) E = − t ; Ψ = √ ( | i + | i ) E = t ; Ψ = √ ( | i − | i ) E = t ; Ψ = √ ( | i − | i ) , (4) for N=2 E = U − c ; Ψ = 4 t | i−| i √ a ( c − U ) + | i + | i a E = 0; Ψ = | i E = 0; Ψ = | i E = 0; Ψ = √ ( | i + | i ) E = U ; Ψ = √ ( | i − | i ) E = U + c ; Ψ = 4 t | i−| i √ b ( c + U ) − | i + | i b , (5)with a = q t ( U − c ) , b = q t ( U + c ) , and c = √ t + U .For N=3 E = U − t ; Ψ = √ ( | i − | i ) E = U − t ; Ψ = √ ( | i − | i ) E = U + t ; Ψ = √ ( | i + | i ) E = U + t ; Ψ = √ ( | i + | i ) , (6)and for N=4 E = 2 U ; Ψ = | i . (7) B. The Partition function, Green’s function andself energy
For convenience, first I introduce ”shifted” many-bodyenergies E ′ Nn = E Nn − µN , with µ being the chemicalpotential. Then I renormalize them, defying the minimal E ′ min among them and subtracting it E ′′ Nn = E ′ Nn − E ′ min .This also factorizes the partition function: Z ( µ ) = X nN e − βE ′ Nn = e − βE ′ min X nN e − βE ′′ Nn = e − βE ′ min Z ′ ( µ ) , (8)where β is the inverse temperature.It is clear that now in every Gibbs average one can use E ′′ ; Z ′ instead of E ′ ; Z which is numerically more stable(big numbers have been subtracted).For exact Green’s function, one obtains through thestandard spectral decomposition: G σij ( ω ) = 1 Z ′ X N X n ∈ N X m ∈ N +1 e − βE ′′ N +1 m + e − βE ′′ Nn iω − E ′′ mN +1 + E ′′ nN × h Ψ Nn | c iσ | Ψ N +1 m ih Ψ N +1 m | c + jσ | Ψ Nn i . (9)The exact self energy is obtained by inversion of theDyson’s equation:Σ σij ( ω ) = G − σ ,ij ( ω ) − G − σij ( ω ) , (10)where Green’s function in Hartree approximation G − σ ,ij ( ω ) = ( iω + µ − U ρ i ) δ ij + t (1 − δ ij ) is used ( ρ i is the occupancy (”density”) of the site i ). C. Response function, polarizability, dielectricfunction, and W
The exact two-point density-density correlation func-tion is also obtained through the spectral decomposition χ ddij ( ν ) = 1 Z ′ X N X n ∈ N X m ∈ N e − βE ′′ Nm − e − βE ′′ Nn iν + E ′′ nN − E ′′ mN × h Ψ Nn | ˆ ρ i | Ψ Nm ih Ψ Nm | ˆ ρ j | Ψ Nn i , (11)with density operator ˆ ρ i = P σ c + iσ c iσ . It is convenient todefine also the density-density response function R ddij ( ν ) = βδ ν ρ i ρ j − χ ddij ( ν ) . (12)After that one can find the density-density dielectricfunction ǫ − ddij ( ν ) = δ ij + U R ddij ( ν ) , (13)the density-density polarizability P ddij ( ν ) = X k R ddik ( ν ) ǫ ddkj ( ν ) , (14) and the screened interaction W ij ( ν ) = U δ ij + U R ddij ( ν ) . (15)The response function, the dielectric function, the po-larizability, and the screened interaction calculated usingthe formulae (12)-(15) from the exact correlation functionare by construction exact and can be compared with thecorresponding quantities obtained using the PT. D. Exact 3-point vertex function in density channel
To find the exact 3-point density vertex function I firstcalculate the following three-point correlation function: χ σ,dijk ( τ ; τ ′ ) = h c iσ ( τ ) c + jσ ( τ ′ )ˆ ρ k (0) i with τ, τ ′ being Mat-subara’s times and ’d’ meaning ’density’. In the site-frequency domain its spectral decomposition reads as thefollowing: χ σ,dijk ( ω ; ν ) = 1 Z ′ X N n − X p ∈ N X n ∈ N +1 h Ψ N +1 n | c + jσ | Ψ Np i i ( ω − ν ) + E ′′ Np − E ′′ N +1 n X m ∈ N +1 h Ψ Np | c iσ | Ψ N +1 m ih Ψ N +1 m | ˆ ρ k | Ψ N +1 n i e − βE ′′ N +1 m − e − βE ′′ N +1 n iν + E ′′ N +1 n − E ′′ N +1 m + X p ∈ N +1 X m ∈ N h Ψ N +1 p | c + jσ | Ψ Nm i− i ( ω − ν ) + E ′′ N +1 p − E ′′ Nm X n ∈ N h Ψ Nn | c iσ | Ψ N +1 p ih Ψ Nm | ˆ ρ k | Ψ Nn i e − βE ′′ Nm − e − βE ′′ Nn iν + E ′′ Nn − E ′′ Nm + X p ∈ N X n ∈ N +1 h Ψ N +1 n | c + jσ | Ψ Np i i ( ω − ν ) + E ′′ Np − E ′′ N +1 n X m ∈ N +1 h Ψ Np | c iσ | Ψ N +1 m ih Ψ N +1 m | ˆ ρ k | Ψ N +1 n i e − βE ′′ N +1 m + e − βE ′′ Np iω + E ′′ Np − E ′′ N +1 m + X p ∈ N +1 X m ∈ N h Ψ N +1 p | c + jσ | Ψ Nm i− i ( ω − ν ) + E ′′ N +1 p − E ′′ Nm X n ∈ N h Ψ Nn | c iσ | Ψ N +1 p ih Ψ Nm | ˆ ρ k | Ψ Nn i e − βE ′′ N +1 p + e − βE ′′ Nn iω + E ′′ Nn − E ′′ N +1 p o . (16)After that the three-point density response function iscalculated R σ,dijk ( ω ; ν ) = δG σij ( ω ) δφ k ( ν ) = βδ ν G σij ( ω ) ρ k + χ σ,dijk ( ω ; ν ) , (17)where I have indicated that the three-point responsefunction is defined as the functional derivative of Green’sfunction with respect to the external perturbing field φ k ( ν ). The ”screened” 3-point density vertex function γ σ,dijk ( ω ; ν ) is defined as the functional derivative of theinverse Green’s function with respect to the external per-turbing field φ k ( ν ) and correspondingly can be related tothe above defined 3-point response function γ σ,dijk ( ω ; ν ) = − δG − σij ( ω ) δφ k ( ν )= X lt G − σil ( ω ) R σ,dltk ( ω ; ν ) G − σtj ( ω − ν ) . (18)Finally the three-point vertex function entering theEq.(1) (”bare” three-point vertex) is defined as the func-tional derivative of the inverse Green’s function with re-spect to the total field Φ k ( ν ) (perturbing external plusinduced internal) and is related to the ”screened” ver-tex through the density-density dielectric matrix (thisequation is the density-density part of the more generalequation (41))Γ σ,dijk ( ω ; ν ) = − δG − σij ( ω ) δ Φ k ( ν )= X l γ σ,dijl ( ω ; ν ) ǫ ddlk ( ν ) . (19) E. Vertex corrected polarizability and self energy
The vertex-corrected density-density polarizability andthe self energy entering the equations (1) are calculatedas the following P ddij ( ν ) = 1 β X ω X σ X kl G σik ( ω )Γ σ,dklj ( ω ; ν ) G σli ( ω − ν ) , (20)andΣ σij ( ω ) = − β X ν X kl G σik ( ω + ν )Γ σ,dkjl ( ω + ν ; ν ) W li ( ν ) . (21)The formulae (20) and (21) are used also when thethree-point vertex is obtained within the perturbationtheory. F. 3-point vertex function from PerturbationTheory in density channel
The first order (in W) term of the perturbation theoryfor the three-point density vertex function isΓ σ,dijk ( τ ; τ ′ ) = − W ji ( τ ′ − τ ) G σik ( τ ) G σkj ( − τ ′ ) . (22)In the frequency representation it can be convenientlyevaluated as the followingΓ σ,dijk ( ω ; ν ) = − Z dτ W ji ( τ ) 1 β X ω e − iωτ G σik ( ω ) G σkj ( ω − ν ) . (23) G. Ward Identities
In order to write down the Ward Identities one needs tospecify the current operator. It is convenient to introduceit through the substitution for the kinetic part of theHamiltonian H ′ = − t X iσ c + iσ e i ( A ˜ i − A i ) c ˜ iσ + U X i X σσ ′ c + iσ c + iσ ′ c iσ ′ c iσ , (24)where A i is the vector potential and the following con-vention for the sites was adopted: ˜1 = 2; ˜2 = 1. Withthe above definition the current operator for the modelreads as the followingˆ J i = − δH ′ δA i = it X σ (cid:8) c + iσ c ˜ iσ − c +˜ iσ c iσ (cid:9) , (25)with ˆ J ˜ i = − ˆ J i . From the equation of motion for thedensity operator one calculates ∂ ˆ ρ i ∂τ = − i ˆ J i , (26)which is the continuity equation for the model. Relat-ing the ”screened” 3-point current vertex function withthe corresponding 3-point response function (similarly tothe density case and in the space-time coordinates forbrevity; ’c’ goes for ’current’)) G σ (14) γ σ,c (453) G σ (52) = G σ (12) J (3)+ h c σ (1) c + σ (2) ˆ J (3) i , (27)one explicitly calculates the time derivatives and usingalso the continuity equation (26) one obtains G σ (14) n ∂∂τ γ σ,d (453) + iγ σ,c (453) o G σ (52)= ∂∂τ h c σ (1) c + σ (2)ˆ ρ (3) i + i h c σ (1) c + σ (2) ˆ J (3) i = G σ (12) n δ (13) − δ (23) o , (28)which is equivalent to ∂∂τ γ σ,d (123) + iγ σ,c (123) = G − σ (12) n δ (23) − δ (13) o . (29)In the site-frequency representation (29) reads as thefollowing i n νγ σ,dijk ( ω ; ν ) + γ σ,cijk ( ω ; ν ) o = G − σij ( ω ) δ jk − G − σij ( ω − ν ) δ ik . (30)The equation (30) can be simplified by removing theHartree-Fock contribution on the both sides of it. For thevertices, one obtains in the Hartree-Fock approximation γ σ,dijk ( ω ; ν ) = δ ik δ jk , (31)and γ σ,cijk ( ω ; ν ) = it n δ ik − δ jk o . (32)For the schemes with full sc (without the QP approx-imation) the removing of the Hartree-Fock contributionalso on the right side of (30) through Dyson’s equationgives i n ν △ γ σ,dijk ( ω ; ν ) + △ γ σ,cijk ( ω ; ν ) o = Σ c,σij ( ω − ν ) δ ik − Σ c,σij ( ω ) δ jk , (33)where △ γ means the ”screened” vertex part beyond theHartree-Fock approximation, and the Σ c is the correla-tion (frequency dependent) part of the self energy. Forthe QP-based schemes the Dyson equation is not satisfiedand instead of (33) one has i n ν △ γ σ,dijk ( ω ; ν ) + △ γ σ,cijk ( ω ; ν ) o = n H HFij − H QPij o [ δ ik − δ jk ] , (34)where the static effective Hamiltonians H HFij and H QPij (correspondingly in the Hartree-Fock and in the QP ap-proximation) were introduced: H HFij = − t (1 − δ ij ) + δ ij ( V Hi + Σ xi ) , (35)with V Hi and Σ xi being the Hartree potential and theexchange part of the self energy correspondingly, and H QPij = µ [ δ ij − Z ij ] + X kl Z / ik [ H HFkl + Σ ckl (0)] Z / lj , (36)where Z ij is the renormalization factor, and the Σ ckl (0)stands for the correlation part of the self energy at zerofrequency. When the effect of correlations is neglectedthe renormalization factor becomes equal to 1, and thecorrelation part of the self energy becomes zero, whichmeans that in this case H QP = H HF .The equations (33) and (34) are used in this work toevaluate the deviations from the full Ward Identity fordifferent approximate methods. I also use two reducedforms of the WI in the present study: the long-wave limitand the long-wave+static limit of the WI. The long-wavelimit of the WI for the two-site Hubbard model consists inthe summation over the index k in the equations (33) and(34) which correspondingly become (the current vertexdisappears after the summation) iν X k △ γ σ,dijk ( ω ; ν ) = Σ c,σij ( ω − ν ) − Σ c,σij ( ω ) , (37)and iν X k △ γ σ,dijk ( ω ; ν ) = 0 . (38)From (38) one can see that in order to satisfy the long-wave limit of the WI in the QP-based approximationsone has to neglect by vertex corrections altogether.The long-wave+static limit of the WI consists in takingthe limit ( ν →
0) in the equations (37) and (38): X k △ γ σ,dijk ( ω ; ν = 0) = lim ν → Σ c,σij ( ω − ν ) − Σ c,σij ( ω ) iν , (39)and X k △ γ σ,dijk ( ω ; ν = 0) = 0 . (40)The arguments supporting the quasiparticle approxi-mation in Ref.[8] are based on the long wave limit of theWard Identity. It is clear from the above considerationthat in the QPGW approximation (without the vertexcorrections) the corresponding limit of the WI is satis-fied exactly. H. 3-point current vertex function
In order to apply the full WI one needs the ”screened”current vertex function γ σ,cijk ( ω ; ν ). In this work both theexact and the PT-based vertex functions are used. Theexact one is evaluated following the formulae (16)-(18)with the replacement (ˆ ρ → ˆ J ) in the Eq.(16). The PT-based vertices are calculated following the scheme out-lined below.The ”bare” vertex functions are related to the”screened” ones through the full (density-current) dielec-tric function via the following equation which is the gen-eralization of Eq.(19) γ σ,Iijk ( ω ; ν ) = X J X l Γ σ,Jijl ( ω ; ν ) ǫ − JIlk ( ν ) , (41)where both I and J now run over the indices d (density)and c (current). The non-perturbed Hamiltonian has novector potentials so that the full dielectric matrix has theform (considering the first index as the density, and thesecond as the current) ǫ = (cid:18) ǫ dd ǫ dc (cid:19) , (42)and correspondingly its inverse ǫ − = (cid:18) ǫ − dd − ǫ − dd ǫ dc (cid:19) . (43)Thus, the bare current vertex can be evaluated as thefollowing γ σ,cijk ( ω ; ν ) = Γ σ,cijl ( ω ; ν ) − X lm Γ σ,dijl ( ω ; ν ) ǫ − ddlm ( ν ) ǫ dcmk ( ν ) . (44)The subtraction of the Hartree-Fock contribution (32)leaves us with the expression △ γ σ,cijk ( ω ; ν ) = △ Γ σ,cijl ( ω ; ν ) − X lm [ δ il δ jl + △ Γ σ,dijl ( ω ; ν )] ǫ − ddlm ( ν ) ǫ dcmk ( ν ) . (45)In the above expression the missing components are thenon-trivial part of the current vertex △ Γ σ,cijl ( ω ; ν ) and thedensity-current dielectric matrix ǫ dcmk ( ν ). The first oneis evaluated similar to the equation (23) for the densityvertex: △ Γ σ,cijk ( ω ; ν ) = − it Z dτ W ji ( τ ) 1 β X ω e − iωτ × n G σik ( ω ) G σ ˜ kj ( ω − ν ) − G σi ˜ k ( ω ) G σkj ( ω − ν ) o . (46)The second one is defined by the density-current po-larizability P dc ǫ dcij ( ν ) = − U P dcij ( ν ) , (47)which in its turn is evaluated using the current vertex(46): P dcij ( ν )= it Z dτ e iντ X σ n G σij ( τ ) G σ ˜ ji ( − τ ) − G σi ˜ j ( τ ) G σji ( − τ ) o + 1 β X ω X σ X kl G σik ( ω ) △ Γ σ,cklj ( ω ; ν ) G σli ( ω − ν ) . (48) I. Internal energy
The exact internal energy is evaluated directly as theaverage value of the Hamiltonian. In the spectral repre-sentation it reads -0.16-0.14-0.12-0.1-0.08-0.06-0.04-0.02 0 0 5 10 15 20 Im Σ Frequency
FIG. 1: (Color online) Imaginary part of the on-site self en-ergy as a function of frequency for different average occupan-cies for U=1. E = 1 Z ′ X nN E ′′ Nn e − βE ′′ Nn . (49)To evaluate the internal energy in the perturbation the-ory based methods the following formula is used E = t X i = j,σ G σji ( τ = β ) + U X i X σσ ′ ρ iσ ρ iσ ′ − X i X σ Σ x,σii G σii ( τ = β )+ 12 β X iσ X ω Σ c,σij ( ω ) G σji ( ω ) , (50)which is based on the Galitskii-Migdal expression forthe exchange-correlation energy. III. RESULTS
The two-site Hubbard model is studied here with thevalue of parameter t fixed and equal to 1. The tempera-ture was also fixed at T=0.05 t . Thus, only the parameterU was changing. The half-filling ( h N i = 2) case has beenconsidered. Such a choice for the occupancy has beenguided mostly by the fact that when one steps aside fromthe half-filling the correlations in the model quickly be-come unmanageable for the QP-based approaches. It isseen in the Fig.1 where the imaginary part of the on-site self energy is plotted as a function of Matsubara’sfrequency for the average occupancies 1, 2, and 3. |G-G Exact | Frequency G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 2: (Color online) Errors in Green’s function for U=0.5. |G-G
Exact | Frequency G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 3: (Color online) Errors in Green’s function for U=1. |G-G
Exact | Frequency G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 4: (Color online) Errors in Green’s function for U=2. |G-G
Exact | Frequency G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 5: (Color online) Errors in Green’s function for U=3. |G-G
Exact | Frequency G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 6: (Color online) Errors in Green’s function for U=4. |G-G
Exact | Frequency G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 7: (Color online) Errors in Green’s function for U=5.
The strong downturn in the function for the occupan-cies 1 and 3 at the small frequencies makes the lineariza-tion in the self energy pertinent to the QP approxima-tions highly inaccurate. At the half filling however, themodel shows slow increasing in the correlation strengthwith the increasing of U and, correspondingly, is conve-nient for the comparative studies. The calculations havebeen performed for different values of U which were in-creased till the methods based on the QP approximationbegan to fail seriously.I compare the exact results with the results obtainedwith the GW, the GΓ W, the GΓ E W, the QPGW,QPGΓ W, and the QPGΓ E W methods, where Γ andΓ E stand for the first order (in W) 3-point vertex andthe exact 3-point vertex correspondingly. At the full self-consistency GΓ E W reproduces the exact results, so thisapproach was used basically only as a mean to adjust(by comparing with the exact result) the calculationalparameters (such as the number of Matsubara’s frequen-cies included in the internal summations and the densityof mesh on the interval [0 : β ] for the τ -integrations.In addition to the Green function which serves in thiswork as basic representation quantity for the compar-isons, the analysis of the Ward Identities fulfillment hasbeen performed and was used as an indicator of the accu-racy of the approximations. Besides, the internal energyhas been evaluated and its accuracy was related to theerrors in the calculated G and to the deviations from theWard Identities.In the figures 2-7 the absolute error in the calculatedGreen’s function is shown as a function of the Matsub-ara’s frequency for the different values of U. The com-parison is being made between the exact result and theresults obtained with the approximate methods. The in-ternal energy as a function of U is presented in the Fig.8.The figure 9 shows the relative violation of the full WI(Eqs.33 and 34) in different approximate methods as afunction of U. The following form of the average devia-tion from the identity has been used P ω,ν P ijk (cid:13)(cid:13)(cid:13) ( LHS ) ijk ( ω ; ν ) − ( RHS ) ijk ( ω ; ν ) (cid:13)(cid:13)(cid:13)P ω,ν P ijk , (51)where ( LHS ) and (
RHS ) are the left- and the right-handsides in (33,34) correspondingly. The 3-point vertex per-tinent to the specific approximation (i.e. for instance thevertex △ γ is zero in the GW and the QPGW) was usedto evaluate the deviation from the WI in all approximateschemes, and the exact vertex was used in GΓ E W andQPGΓ E W. The summations over the frequencies ω and ν in (51) were performed for | ω | <
20 and 0 ≤ ν < k -summation.Finally in the figure 11 the deviation from the long-wave+static limit of the WI (39,40) is shown. In thiscase, the summation over ν has not been included. -2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4 0 1 2 3 4 5 Internal Energy
UExact, G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 8: (Color online) Internal energy as a function of Uparameter.
One can do a few observations from the results pre-sented.Most important for the present work is an observationfollowing from the errors in Green’s function, that theQPGΓ E W method doesn’t show any noticeable improve-ment as compared to the other approximate methodsstudied. Opposite to that - with U increasing it quicklybecomes the worst of the methods considered (especiallyfor the higher frequencies). The observation perfectlycorrelates with the results obtained for the internal en-ergy and for the Ward Identities (especially with the re-sults for the full WI). It means that putting the exactvertex together with the QP self-consistency is not com-patible.Grouping the methods based on the QP-approximationin one group and the rest of the methods in another(they satisfy the DE), one can do another conclusion.
Deviation from WI
UExact, G Γ E WGWG Γ WQPGWQPG Γ WQPG Γ E W FIG. 9: (Color online) Full Ward Identity average violationas a function of U parameter. → Γ → Γ E ) consistently improves the accuracy of G, finally mak-ing it exact when the exact vertex is being applied. Nat-urally, the same tendency can be noticed looking at theinternal energy graph and at all WI graphs.Opposite to that, in the QP-based methods the situ-ation is more complicated. At lowest U (0.5 and 1) theQPGW is the best and the QPGΓ E W is the worst at lowfrequencies with the reversed tendency at high frequen-cies. For U equal 2, 3, and 4 one can notice that nowthe QPGΓ E W is the best and the QPGW is the worst atlow frequencies with the reversed tendency at high fre-quencies. For U equal 5 the accuracy of the QPGΓ E Wapproach deteriorates also at low frequencies and its fail-ing at higher frequencies becomes severe. The importantpoint here is, that the above high-frequency tendenciesin the QP-based methods correlate well with the tenden-cies in the internal energy as one can easily see. Namely,for U ≤ E W gives the best (among theQP-based schemes) internal energy and for U > Deviation from WI
UExact, G Γ E W, QPGWGWG Γ WQPG Γ WQPG Γ E W FIG. 10: (Color online) q → Deviation from WI
UExact, G Γ E W, QPGWGWG Γ WQPG Γ WQPG Γ E W FIG. 11: (Color online) q → , ν → U=1. As one can see there is a strict separation in theaccuracy of G obtained within the QP-group and withinthe non QP-group. The reasons for this feature have notbeen studied in this work but it could be related to theabove mentioned swapping between different tendenciesin the QP-based methods at this U.The GW-based (without the QP approximation) meth-ods give numerically exact internal energy till U ≈ ) and in the weakly correlatedregime should produce accurate total energies. Generallyit is not a good idea to use the QP-based methods to eval-uate the total energy because even if they originate fromthe same Ψ-functional as the GW-based methods theDyson equation (DE) is not satisfied in them anymore.As it was indicated above, the deviations from the exactGreen’s function in the QP-based approaches are espe-cially noticeable at high frequencies which are importantfor the total energies evaluation. The results shown inFig.8 clearly support this point of view.Internal energy obtained in QP-based approximationsdeviates significantly from the exact one even at smallU. It is interesting that in the QPGW scheme the er-ror remains almost unchanged till U = 5 which prob-ably is accidental because at U ≥ W improves the internal energy ascompared to the GW, whereas applying the vertex ofimproved accuracy (1 → Γ → Γ E ) in the sequenceQPGW → QPGΓ W → QPGΓ E W makes the results worseand worse which also tells us that one should not applythe vertex correction combined with the QP approxima-tion.As it is clear, the degree of violation of the WI corre-lates well with the errors in the Green functions at av-erage and high frequencies, which are responsible for theaccuracy of the calculated total energies. And indeed,1the comparison of the Figures 8 and 9 tells us, that theaccuracy of the total energy and the accuracy of the fullWI fulfillment are closely related. Thus, the full WI canbe useful as a measure of the accuracy of the calculatedtotal energies.There was a hope, that the deviation from the long-wave+static limit of the WI correlates well with the er-rors in G at low frequencies. But for the two-site Hub-bard model this seems to be not the case (besides theabove mentioned success of the QPGW approach at thesmallest values of U). However, in real materials wherethe spectra obtained with the QPGW are generally no-ticeably better than the spectra obtained with the scGWthe situation might be different.
IV. CONCLUSIONS
In this work the Hedin’s equations for the two-siteHubbard model have been solved self-consistently withand without applying the quasiparticle approximation forthe Green function. The study has been performed bothwhen the exact three-point vertex function was used as aninput and when the perturbative theory (in its zero andfirst orders in W) was used to evaluate the correspond-ing vertices self-consistently. The results of this work ob-tained with the exact vertex have direct impact on whatone can expect from the future implementation of the scGW+DMFT and the sc QPGW+DMFT schemes. Asit has been shown here, only the GW+DMFT approachcan be considered as useful approximation. However, asit was said in the Introduction, this work deals with anideal situation, when the GW part and the DMFT partare perfectly separable and the subspace where GW isused is very weakly correlated (so that GW and QPGWgive identical results for the weakly correlated subspace).In practice, it is not always the case. The correlations in the ”GW”-subspace might be noticeable. In such situa-tion, the QPGW might be superior (with respect to theGW) for the subspace not included in the DMFT part.The conclusions about the DMFT part obviously remainas before - one should not impose the QP approximationon G in the DMFT part. As it seems, in such circum-stances the preference should be given to the approach(GW-based or QPGW-based) depending on which sub-space (the ”weakly” or the ”strongly” correlated) is moreimportant for the problem under consideration. However,on the fundamental level, such a situation should be re-solved either by the increasing of the subspace coveredby the DMFT part or (which seems to be easier practi-cally) by including more diagrams beyond GW for the”weakly” correlated subspace.It has been shown, that the methods with the PT-based vertices (when they are applicable) reveal similartendencies (for example if one chooses between the QPself-consistency and the full sc) as the methods basedon the exact vertices. Namely, when the correlationstrength increases, both the PT-based and the exactvertices-based schemes begin to fail if one uses the QPself-consistency. Also of practical importance is the find-ing that the violation of the WI in any particular methodcorrelates well with the general applicability of the givenmethod. This can be an useful information if one seesto apply the schemes with vertex corrections to the realmaterials where the exact solutions are not available toserve as a judgement.The results of this work are of the methodological im-portance. Of course the two-site HM doesn’t cover allpossible regimes of correlations which may happen in re-alistic materials. In order to cover a little bit more of thepossible regimes of correlations the similar work on thehomogeneous electron gas is now being performed (withthe three-point vertex functions calculated within the PTonly). P. Sun and G. Kotliar, Phys. Rev.B , 085120 (2002). S. Biermann, F. Aryasetiawan, and A. Georges,Phys. Rev. Lett. , 086402 (2003). J. M. Tomczak, M. Casula, T. Miyake, and S. Biermann,Phys. Rev. B , 165138 (2014). A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein,Phys. Rev. B , 035122 (2005). P. Werner, A. Comanac, L. de Medici, M. Troyer, andA. J. Millis, Phys. Rev. Lett. , 076405 (2006). A. Kutepov, K. Haule, S. Y. Savrasov, and G. Kotliar,Phys. Rev. B , 155129 (2012). A. Kutepov, S. Y. Savrasov, and G. Kotliar, Phys. Rev. B , 041103 (2009). T. Kotani and M. van Schilfgaarde, S. V. Faleev,Phys. Rev.B , 165106 (2007). M. van Schilfgaarde, T. Kotani, and S. Faleev,Phys. Rev. Lett. , 226402 (2006). J. M. Tomczak, M. van Schilfgaarde, and G. Kotliar,Phys. Rev. Lett. , 237110 (2012). J. Tomczak, arXiv.cond.mat.:1411.5180 (2014). L. Hedin, Phys. Rev. , A796 (1965). G. Baym, and L. P. Kadanoff, Phys. Rev. , 287 (1961). C.-O. Almbladh, U. von Barth and R. van Leeuwen, Int. J.of Mod.Phys. B13