Systematically improvable multi-scale solver for correlated electron systems
SSystematically improvable multi-scale solver for correlated electron systems
Alexei A. Kananenka, Emanuel Gull, and Dominika Zgid Department of Chemistry, University of Michigan, Ann Arbor, Michigan, 48109, USA Department of Physics, University of Michigan, Ann Arbor, Michigan, 48109, USA
The development of numerical methods capable of simulating realistic materials with stronglycorrelated electrons, with controllable errors, is a central challenge in quantum many-body physics.Here we describe a framework for a general multi-scale method based on embedding a self-energy ofa strongly correlated subsystem into a self-energy generated by a method able to treat large weaklycorrelated systems approximately. As an example, we present the embedding of an exact diagonal-ization self-energy into a self-energy generated from self-consistent second order perturbation theory.Using a quantum impurity model, generated from a cluster dynamical mean field approximation tothe 2D Hubbard model, as a benchmark, we illustrate that our method allows us to obtain accurateresults at a fraction of the cost of typical Monte Carlo calculations. We test the method in multipleregimes of interaction strengths and doping of the model. The general embedding framework wepresent avoids difficulties such as double counting corrections, frequency dependent interactions, orvertex functions. As it is solely formulated at the level of the single-particle Green’s function, itprovides a promising route for the simulation of realistic materials that are currently difficult tostudy with other methods.
PACS numbers: 71.15.-m, 71.20.-b, 71.30.+h,71.10.Fd
I. INTRODUCTION AND GENERALFRAMEWORK
The theoretical description of strongly correlated ma-terials has proven to be challenging, mainly becausemany of their interesting properties are caused by theinterplay of subtle electronic correlation effects on low en-ergy scales. Since simultaneous treatment of both strongand weak correlations is of major importance for thequantitative description of these systems, two main con-ceptual approaches are used: the reduction to a few ‘rel-evant’ degrees of freedom or essential orbitals around theFermi level and the subsequent construction of a modelsystem or, alternatively, the treatment of the entire sys-tem using methods which significantly approximate cor-relation effects.The first approach, with methods including exact di-agonalization (ED) and its variants, density ma-trix renormalization group (DMRG), dynamical meanfield theory (DMFT) and lattice quantum Monte Carlo(QMC) applied to model Hamiltonians, can yield veryprecise results for model systems. When applied to real-istic systems, its main uncertainties and possible sourcesof errors lie in the construction of the parameters of theeffective model.The second approach, which includes implementationsof the density functional theory (DFT), HartreeFock (HF), GW, the random phase approximation(RPA), Møller-Plesset second order perturbation the-ory (MP2), GF2 or QMC, avoids constructing aneffective model by treating the full Hamiltonian with allorbitals and interactions, but relies on potentially severeapproximations to the electronic correlations.Multi-scale methods for extended systems combiningthe best aspects of both approaches, e.g. by solving thesystem using DFT or GW and using the result to con- struct a model system, have been implemented as theGW+DMFT and DFT+DMFT method.Constructing a robust multi-scale method is aformidable problem and an active field of research. First,different energy scales have to be defined and a set ofstrongly correlated orbitals requiring a higher level treat-ment has to be chosen. Second, the non-local Coulombinteractions present in realistic materials have to be in-cluded by a suitable choice of ‘screened’ interactions.Third, correlations in the weakly correlated orbitalsshould not be completely neglected but rather be treatedperturbatively, if a quantitative material-dependent de-scription is desired.In this paper, we present a general framework for amulti-scale algorithm in which a self-energy describingstrongly correlated orbitals is self-consistently embeddedinto the a self-energy obtained from a method able totreat long-range interaction and correlation effects. Wecall this general framework ‘self-energy embedding the-ory’ (SEET). A mathematically rigorous procedure foridentifying the strongly correlated orbitals and system-atically increasing the accuracy of the treatment of theweakly correlated orbitals is an integral part of our proce-dure. The strongly and the weakly correlated subspacesare treated using different methods. As an example, wewill choose exact diagonalization (ED) to yield the self-energy for the strongly correlated orbitals and the self-consistent second order Green’s function method (GF2) for the weakly correlated part, and call the resulting al-gorithm SEET(ED-in-GF2). However,we note that ourscheme is general and independent of the algorithms usedto treat the weakly and strongly treated subspaces: EDcould be replaced by, e.g. , CT-QMC, while GF2 could bereplaced with FLEX, GW, or the Parquet method.Here we illustrate SEET(ED-in-GF2) and calibrateit using impurity problems, since a method capable of a r X i v : . [ c ond - m a t . s t r- e l ] A p r treating multi-scale problems such as realistic materialsshould also yield accurate results for model systems. Im-purity problems have a continuous dispersion but only afinite number of interaction terms. Nevertheless, they ex-hibit strongly and weakly correlated regimes and a num-ber of phases and phase transitions that are very wellunderstood and for which numerically exact comparisonalgorithms exist. The restriction to these models allowsus to provide stringent tests on the accuracy of SEET in awell controlled test environment. In the future, we aim toapply SEET to realistic materials, where it has the poten-tial to become a method complementary to GW+DMFTor DFT+DMFT.In order to generate a wide range of correlated phases,we generate our impurity parameters from a four-site dy-namical cluster approximation (DCA) to the 2D Hub-bard model, i.e. test our method as a ‘DCA impuritysolver’, so that our results can be compared against nu-merically exact continuous time CT-QMC data.
Weemphasize that while we have developed a numericallyefficient impurity solver, we do not envisage this as themain use of SEET, and we mainly resort to impuritymodels for the sake of generating comparisons to reliableknown results.SEET in the ED-in-GF2 variant is computationally af-fordable, as GF2 scales as O ( N ), where N is the numberof orbitals in a unit cell. Additionally, SEET is amenableto parallelization on large machines. The scaling of GF2can further be reduced to O ( N ) by employing densityfitted integrals, and the strongly correlated orbitals canbe treated by ED as pairs, further reducing the numer-ical cost. Consequently, large systems containing manyunit cells or k -points containing multiple orbitals can betreated simultaneously, providing non-local effects andmomentum dependence.In Section II, we introduce SEET(ED-in-GF2). Sec. IIIshows results for our test model, and Sec. IV containsconclusions of our work. II. THE SEET(ED-IN-GF2) METHOD APPLIEDTO AN IMPURITY MODEL
We consider an impurity problem with N impurity or-bitals a i coupled to an infinitely many bath orbitals c µ described by a general Hamiltonianˆ H = (cid:88) ij t ij a † i a j + (cid:88) ijkl U ijkl a † i a † j a l a k ++ (cid:88) iλ V iλ a † i c λ + (cid:88) λ (cid:15) λ c † λ c λ + h.c., (1)where t and U are material specific one- and two-bodyoperators, V is the hybridization strength, and (cid:15) λ is the c -electron dispersion. The single-particle properties of thisHamiltonian are described by a non-interacting Matsub-ara Green’s function for a -electrons G ( iω ) = [ iω + µ − t − ∆( iω )] − , (2) with ∆( iω ) encapsulating the properties of the c -electrons and µ being the chemical potential. InSEET(ED-in-GF2), we obtain the interacting Green’sfunction G GF of this N -orbital impurity problem it-eratively, starting from G GF = G , by self-consistentsecond order perturbation theory (GF2), [Σ GF ( τ )] ij = − (cid:88) klmnpq [ G GF ( τ )] kl [ G GF ( τ )] mn ×× [ G GF ( − τ )] pq U iqmk (cid:0) U lnpj − U nlpj (cid:1) (3)and the corresponding GF2 Green’s function G GF ( iω ) = (cid:2) [ G GF ( iω )] − − Σ GF ( iω ) (cid:3) − . (4)Note that GF2 can be solved self-consistently and in-cludes an exchange diagram important for describing sys-tems with a localized electronic density, but does not in-clude higher order RPA-like diagrams that are present, e.g. , in GW. We then evaluate the one-body densitymatrix using the converged GF2 Green’s function andchoose a set of n < N orbitals corresponding to eigenval-ues of the one-body density matrix which are significantlydifferent from 0 or 2. These n orbitals, which we will call‘strongly correlated’, are used to build an n -orbital im-purity problem which is then solved with a method moreaccurate than GF2 to compute a self-energy. Here, weuse ED to solve this impurity problem. The resultingED self-energy is used to modify the GF2 self-energy andto obtain the total self-energy in the natural orbital basisas [Σ] ij = [Σ EDstrong ] µν + [Σ GF ] ij − [Σ GF strong ] µν . (5)The indices i and j run over all N orbitals, while µ and ν run only over the n strongly correlated orbitals. Thetotal self-energy is schematically illustrated in Fig. 1.As the n correlated orbitals are chosen in the eigenbasisof the one-body density matrix, a transformation of theone-body and two-body integrals in this n -orbital sub-space to the eigenbasis is necessary. Note, that even forcases where model Hamiltonians with simplified (e.g. lo-cal or density-density) interaction structures are studied,this transformation generates general interactions U ijkl .This n -orbital impurity problem with non-local interac-tion U ijkl is then treated by the ED solver requiring anadditional bath discretization step which may introducefitting errors. We emphasize that these are small for thecases studied here and that, in principle, any solver suit-able to describe strong correlations and able to treat gen-eral multi-orbital interactions can be employed, includ-ing QMC solvers based on the hybridization expansion which do not require a bath discretization step.The ED-in-GF2 procedure is iterated, and the GF2calculation updates [Σ GF ] ij for the N orbitals since i, j = 1 , . . . , N where[Σ GF weak ] µν = [Σ GF ] µν − [Σ GF strong ] µν (6) + total frequency dependentself-energy accurate solver perturbative solver ⌃( i! ) = ⌃ GF ( i! ) ⌃ ED strong ( i! ) ⌃ ED strong ( i! ) ⌃ GF weak ( i! ) ⌃ GF weak ( i! ) ⌃ GF ( i! ) FIG. 1. The total self-energy in natural orbital basis pro-duced in the SEET with ED-in-GF2 scheme. is responsible for removal of diagrams later included atthe ED level. Subsequent ED calculation updates thestrongly part of self-energy [Σ
EDstrong ] µν . The iterative up-dates stop when the total self-energy in Eq. 5 is convergedto a predefined accuracy. We present a detailed algorith-mic description of this framework in the supplementarymaterial.The algorithm, both in its general form and in the ED-in-GF2 variant, is based on a diagrammatic formulationin which a ‘double counting’ problem does not appear.The single-particle formulation avoids vertex functions,which are often difficult to handle, and is based on static(or frequency-independent) interactions. III. RESULTS
We calibrate SEET(ED-in-GF2) for the 2 × and consequently the Hamiltonian from Eq. 1is defined for t describing nearest neighbor hopping onlyand U exclusively on-site interactions. DCA provides thenon-interacting Green’s function (in Eq. 2) which is thenemployed to obtain the GF2 self-energy from Eq. 3. Sub-sequently, we construct the one-body density matrix andchoose a pair of two-site impurities to be treated by ED.The occupations of the four site cluster in natural or-bitals are 2- x , 1, 1, x , where for most regimes x is nota small number, thus the orbitals with occupations 2- x ,and x are not any longer weakly correlated. This moti-vates us to choose two separate impurity problems withorbitals occupied as (1,1) and (2- x , x ) and treat them asa pair of two-site impurities embedded into the GF2 de-scription. This means that only the interactions betweenthese two-site impurities are treated at the GF2 level. Aschematic description of the DCA+ED-in-GF2 iterativescheme is shown in Fig. 2. SEET allows us to treat mul-tiple embedded impurities which is computationally ad-vantageous since in realistic cases the number of stronglycorrelated orbitals may be too large for current solverssuch as ED or the hybridization expansion.Testing ED-in-GF2 using SEET on the DCA approxi-mation to the 2D Hubbard model provides a worst casescenario for a multi-scale embedding scheme, since inmultiple regimes the 4-site cluster does not display a sep-aration of energy scales or any ‘weakly’ and ‘strongly’correlated orbitals as typically found in realistic mate- f r o m D C A k-space real-space G F SEETED-in-GF2 b a c k t o D C A G DCA ( i! n ) G GF ( i! ) = ⇥ [ G DCA ( i! )] ⌃ GF ( i! ) ⇤ natural orbitals density matrixAIM 1 with ED AIM 2 with ED G ( i! ) = ⇥ [ G ( i! )] (⌃ EDstrong +⌃ GF weak )] FIG. 2. Schematic view of the DCA+ED-in-GF2 procedureused for treating the 4-site cluster DCA approximation to the2D Hubbard model. Note that ‘strong’ denotes quantities inthe correlated subspace, and ‘weak’ quantities defined on theentire system. −0.28−0.24−0.20−0.16−0.12−0.08 0 5 10 15 20 I m Σ ( i ω n ) U/t=3CTQMCGF2ED−in−GF2 −0.50−0.40−0.30−0.20−0.10 0 5 10 15 20U/t=4CTQMCGF2ED−in−GF2−5.00−4.00−3.00−2.00−1.000.00 0 5 10 15 20 I m Σ ( i ω n ) ω n U/t=6CTQMCGF2ED−in−GF2 −12.00−10.00−8.00−6.00−4.00−2.000.00 0 5 10 15 20 ω n U/t=8CTQMCGF2ED−in−GF2
FIG. 3. The imaginary part of the on-site CT-QMC, GF2and ED-in-GF2 self-energy obtained for a 4-site cluster ofthe half-filled 2D Hubbard model for various regimes with β = 10 t . rials. Rather, in the Mott regime of the 2D Hubbardmodel, all orbitals are strongly correlated, providing astringent test of the SEET with ED-in-GF2 method.In Fig. 3, the imaginary part of the self-energy isplotted for the half-filled case. For weak coupling, i.e. U/t < , GF2 recovers the QMC results well. While for
U/t = 3 ED-in-GF2 corrects the GF2 result only slightly,for
U/t = 4, the improvement is more substantial. Inthis case, the ED-in-GF2 recovers QMC results and is aquantitative correction to the qualitatively correct GF2 −0.60−0.50−0.40−0.30−0.20−0.10 0 1 2 3 4 5 I m G ( i ω n ) U/t=3CTQMCGF2ED−in−GF2 −0.50−0.40−0.30−0.20−0.10 0 1 2 3 4 5U/t=4CTQMCGF2ED−in−GF2−0.45−0.40−0.35−0.30−0.25−0.20−0.15−0.10 0 1 2 3 4 5 I m G ( i ω n ) ω n U/t=5CTQMCGF2ED−in−GF2 −0.40−0.35−0.30−0.25−0.20−0.15−0.10 0 1 2 3 4 5 ω n U/t=6CTQMCGF2ED−in−GF2
FIG. 4. The imaginary part of the on-site CT-QMC, GF2and ED-in-GF2 Matsubara Green’s function obtained for a4-site cluster of the 2D Hubbard model at 10% doping, for
U/t = 3 , , , and 6 with β = 10 t . curve.As expected, in the Mott regime, U/t = 6 and 8, GF2fails to recover the self-energy even qualitatively. Notethat an IPT-like fitting of the large- U limit to the atomiclimit would be possible for this particular example butnot in general, as it requires the determination of thelocal physics at exponential (in n ) cost. In the Mottregime, ED-in-GF2 recovers to a decent quantitative ac-curacy the QMC self-energy for both U/t = 6 and 8.In Figs. 4 and 5, we examine several interesting regimesat 10% doping, where the system exhibits the behavior ofa strongly correlated Fermi liquid. In these cases, we re-port real and imaginary parts of Green’s functions ratherthan self-energies, since a slight difference in chemical po-tentials between different methods results in a shift of theHartree term in the large- ω limit.The imaginary part of Green’s function shows a goodquantitative agreement between CT-QMC and ED-in-GF2 for multiple U/t regimes. The real part of Green’sfunction shows more differences than the imaginary part.In the weak coupling regime illustrated in Fig. 5, for U/t = 3 and
U/t = 4, all the QMC, GF2 and ED-in-GF2real parts of Green’s functions are close. The
U/t = 5and
U/t = 6 regimes are more correlated and GF2 yieldsa qualitatively incorrect result. ED-in-GF2 corrects thisresult and provides a quantitative agreement with CT-QMC.
IV. CONCLUSIONS
We introduced a general self-energy embedding theory(SEET) for correlated systems and performed a compari- R e G ( i ω n ) U/t=3CTQMCGF2ED−in−GF2 −0.010.000.010.020.030.040.050.060.07 0 5 10 15 20 25 30U/t=4CTQMCGF2ED−in−GF20.000.010.020.030.04 0 5 10 15 20 25 30 R e G ( i ω n ) ω n U/t=5CTQMCGF2ED−in−GF2 −0.0050.0000.0050.0100.0150.020 0 5 10 15 20 25 30 ω n U/t=6CTQMCGF2ED−in−GF2
FIG. 5. The real part of the on-site CT-QMC, GF2 andED-in-GF2 Matsubara Green’s function obtained for a 4-sitecluster of the 2D Hubbard model at 10% doping, for
U/t =3 , , , and 6 with β = 10 t . son of SEET(ED-in-GF2) on a strongly correlated systemfor which the exact solution is known, the 4 site clusterDCA approximation to the 2D Hubbard model. Thismodel has a continuous dispersion and shows a rangeof correlated phases, thus providing us with a detailedassessment of strengths and weaknesses of our method.However, it does not illustrate the effect of non-local in-teractions. Since in multiple regimes a clear separationof energy scales is not present, this model provides a rig-orous test for a multi-scale method. We were able toshow that ED-in-GF2 provides accurate results for the4 site Hubbard model in the weakly correlated, interme-diately correlated, and strongly correlated regimes, atand away from half-filling. While the solution in thestrongly correlated embedded subset of orbitals has ex-ponential scaling in our case, the total self-energy forthe strongly correlated orbitals can be assembled usingsolutions of multiple small impurity problems. The cal-culation of the properties of the weakly coupled orbitalswith GF2 scales as O ( N ), making SEET(ED-in-GF2) anideal tool for the simulation of realistic materials. Exten-sions using other diagrammatic or correlated methods,such as SEET(QMC-in-GF2) or methods based on GWare straightforward.In real materials, the number of weakly correlated or-bitals in the unit cell is significantly larger than the num-ber of strongly correlated orbitals, thus providing an idealsituation where many orbitals can be treated cheaply byGF2 while the number of orbitals treated by ED remainssmall. Moreover, the SEET(ED-in-GF2) hybrid is easyto implement and, since it does not use frequency depen-dent effective interactions, can be trivially extended toemploy different solvers for the strongly correlated part,such as truncated CI variants with a suitably chosen ac-tive space or QMC hybridization expansions. Similarly,the weakly correlated part can be treated by different lev-els of perturbation theory or cheap truncated CI meth-ods, instead of GF2. Our ED-in-GF2 method can be ad-justed to yield more accurate results, either by increasingthe order of the perturbative treatment (e.g. by employ-ing FLEX or GW), or by increasing the number of or-bitals treated by ED. These limits therefore provide a rig-orous assessment of the convergence of the self-energies.Since a set of strongly correlated orbitals in SEET is cho-sen based on a unique mathematical criterion,the method has the potential for becoming a black box method forrealistic correlated materials calculations. ACKNOWLEDGMENTS
D.Z. and A.K. acknowledge support by DOE ER16391,E.G. support from the Simons collaboration on the many-electron problem. We thank P. Werner for helpful com-ments on the first draft of this manuscript. E. R. Davidson, Journal of Computational Physics , 87(1975). H. Lin and J. Gubernatis, Computers in Physics , 400(1993). A. Liebsch and H. Ishida, Journal of Physics: CondensedMatter , 053201 (2012). M. Capone, L. de’ Medici, and A. Georges, Phys. Rev. B , 245116 (2007). D. Zgid and G. K.-L. Chan, J. Chem. Phys. , 094115(2011). D. Zgid, E. Gull, and G. K.-L. Chan, Phys. Rev. B ,165128 (2012). S. R. White, Phys. Rev. Lett. , 2863 (1992). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev.Mod. Phys. , 1027 (2005). R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys.Rev. D , 2278 (1981). P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). R. Parr and W. Yang,
Density Functional Theory of Atomsand Molecules (Oxford Science, 1989). L. Hedin, Phys. Rev. , A796 (1965). D. Bohm and D. Pines, Phys. Rev. , 625 (1951). C. Møller and M. S. Plesset, Phys. Rev. , 618 (1934). J. J. Phillips and D. Zgid, J. Chem. Phys. , 241101(2014). N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. ,164102 (2005). W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, Rev. Mod. Phys. , 33 (2001). S. Biermann, F. Aryasetiawan, and A. Georges, Phys. Rev.Lett. , 086402 (2003). S. Biermann, F. Aryasetiawan, and A. Georges, in
Physicsof Spin in Solids: Materials, Methods and Applications ,NATO Science Series II: Mathematics, Physics and Chem-istry, Vol. 156, edited by S. Halilov (Springer Netherlands,2005) pp. 43–65. K. Karlsson, Journal of Physics: Condensed Matter ,7573 (2005). G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko,O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. ,865 (2006). K. Held, Adv. Phys. , 829 (2007). J. M. Tomczak, M. Casula, T. Miyake, F. Aryasetiawan,and S. Biermann, EPL (Europhysics Letters) , 67001 (2012). M. Casula, P. Werner, L. Vaugier, F. Aryasetiawan,T. Miyake, A. J. Millis, and S. Biermann, Phys. Rev.Lett. , 126408 (2012). M. Casula, A. Rubtsov, and S. Biermann, Phys. Rev. B , 035115 (2012). T. Ayral, P. Werner, and S. Biermann, Phys. Rev. Lett. , 226401 (2012). C. Taranto, M. Kaltak, N. Parragh, G. Sangiovanni,G. Kresse, A. Toschi, and K. Held, Phys. Rev. B ,165119 (2013). R. Sakuma, P. Werner, and F. Aryasetiawan, Phys. Rev.B , 235110 (2013). P. Hansmann, T. Ayral, L. Vaugier, P. Werner, andS. Biermann, Phys. Rev. Lett. , 166401 (2013). T. Ayral, S. Biermann, and P. Werner, Phys. Rev. B ,125149 (2013). S. Biermann, Journal of Physics: Condensed Matter ,173202 (2014). E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov,M. Troyer, and P. Werner, Rev. Mod. Phys. , 349(2011). M. H. Hettler, M. Mukherjee, M. Jarrell, and H. R. Kr-ishnamurthy, Phys. Rev. B , 12739 (2000). E. Gull, P. Werner, O. Parcollet, and M. Troyer, EPL(Europhysics Letters) , 57003 (2008). E. Gull, P. Staar, S. Fuchs, P. Nukala, M. S. Summers,T. Pruschke, T. C. Schulthess, and T. Maier, Phys. Rev.B , 075122 (2011). M. Caffarel and W. Krauth, Phys. Rev. Lett. , 1545(1994). P. Werner and A. J. Millis, Phys. Rev. B , 155107 (2006). M. Karolak, G. Ulm, T. Wehling, V. Mazurenko,A. Poteryaev, and A. Lichtenstein, Journal of ElectronSpectroscopy and Related Phenomena , 11 (2010),proceedings of International Workshop on Strong Corre-lations and Angle-Resolved Photoemission Spectroscopy2009. F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar,S. Biermann, and A. I. Lichtenstein, Phys. Rev. B ,195104 (2004). L. Vaugier, H. Jiang, and S. Biermann, Phys. Rev. B ,165105 (2012). R. Sakuma and F. Aryasetiawan, Phys. Rev. B87