Systematically improvable multi-scale solver for correlated electron systems

Systematically improvable multi-scale solver for correlated electron systems

Alexei A. Kananenka Department of Chemistry, University of Michigan, Ann Arbor, Michigan, 48109, USA    Emanuel Gull Department of Physics, University of Michigan, Ann Arbor, Michigan, 48109, USA    Dominika Zgid Department of Chemistry, University of Michigan, Ann Arbor, Michigan, 48109, USA
Abstract

The development of numerical methods capable of simulating realistic materials with strongly correlated 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 of a strongly correlated subsystem into a self-energy generated by a method able to treat large weakly correlated systems approximately. As an example, we present the embedding of an exact diagonalization 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 to the 2D Hubbard model, as a benchmark, we illustrate that our method allows us to obtain accurate results at a fraction of the cost of typical Monte Carlo calculations. We test the method in multiple regimes of interaction strengths and doping of the model. The general embedding framework we present avoids difficulties such as double counting corrections, frequency dependent interactions, or vertex functions. As it is solely formulated at the level of the single-particle Green’s function, it provides a promising route for the simulation of realistic materials that are currently difficult to study with other methods.

pacs:
71.15.-m, 71.20.-b, 71.30.+h,71.10.Fd

I Introduction and General Framework

The theoretical description of strongly correlated materials has proven to be challenging, mainly because many of their interesting properties are caused by the interplay of subtle electronic correlation effects on low energy scales. Since simultaneous treatment of both strong and weak correlations is of major importance for the quantitative description of these systems, two main conceptual approaches are used: the reduction to a few ‘relevant’ degrees of freedom or essential orbitals around the Fermi level and the subsequent construction of a model system or, alternatively, the treatment of the entire system using methods which significantly approximate correlation effects.

The first approach, with methods including exact diagonalization (ED) Davidson (1975); Lin and Gubernatis (1993); Liebsch and Ishida (2012); Capone et al. (2007) and its variants,Zgid and Chan (2011); Zgid et al. (2012) density matrix renormalization group (DMRG),White (1992) dynamical mean field theory (DMFT)Georges et al. (1996); Maier et al. (2005) and lattice quantum Monte Carlo (QMC)Blankenbecler et al. (1981) applied to model Hamiltonians, can yield very precise results for model systems. When applied to realistic systems, its main uncertainties and possible sources of errors lie in the construction of the parameters of the effective model.

The second approach, which includes implementations of the density functional theory (DFT),Hohenberg and Kohn (1964, 1964); Parr and Yang (1989) Hartree Fock (HF), GW,Hedin (1965) the random phase approximation (RPA),Bohm and Pines (1951) Møller-Plesset second order perturbation theory (MP2),Møller and Plesset (1934) GF2 Phillips and Zgid (2014); Dahlen and van Leeuwen (2005) or QMC,Foulkes et al. (2001) avoids constructing an effective model by treating the full Hamiltonian with all orbitals and interactions, but relies on potentially severe approximations to the electronic correlations.

Multi-scale methods for extended systems combining the best aspects of both approaches, e.g. by solving the system using DFT or GW and using the result to construct a model system, have been implemented as the GW+DMFTBiermann et al. (2003, 2005); Karlsson (2005); Kotliar et al. (2006); Held (2007); Tomczak et al. (2012); Casula et al. (2012a, b); Ayral et al. (2012); Taranto et al. (2013); Sakuma et al. (2013); Hansmann et al. (2013); Ayral et al. (2013); Biermann (2014) and DFT+DMFTKotliar et al. (2006); Held (2007) method.

Constructing a robust multi-scale method is a formidable problem and an active field of research. First, different energy scales have to be defined and a set of strongly correlated orbitals requiring a higher level treatment has to be chosen. Second, the non-local Coulomb interactions present in realistic materials have to be included by a suitable choice of ‘screened’ interactions. Third, correlations in the weakly correlated orbitals should not be completely neglected but rather be treated perturbatively, if a quantitative material-dependent description is desired.

In this paper, we present a general framework for a multi-scale algorithm in which a self-energy describing strongly correlated orbitals is self-consistently embedded into the a self-energy obtained from a method able to treat long-range interaction and correlation effects. We call this general framework ‘self-energy embedding theory’ (SEET). A mathematically rigorous procedure for identifying the strongly correlated orbitals and systematically increasing the accuracy of the treatment of the weakly correlated orbitals is an integral part of our procedure. The strongly and the weakly correlated subspaces are treated using different methods. As an example, we will 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)Phillips and Zgid (2014) for the weakly correlated part, and call the resulting algorithm SEET(ED-in-GF2). However,we note that our scheme is general and independent of the algorithms used to treat the weakly and strongly treated subspaces: ED could be replaced by, e.g., CT-QMC, while GF2 could be replaced with FLEX, GW, or the Parquet method.

Here we illustrate SEET(ED-in-GF2) and calibrate it using impurity problems, since a method capable of treating multi-scale problems such as realistic materials should also yield accurate results for model systems. Impurity problems have a continuous dispersion but only a finite number of interaction terms. Nevertheless, they exhibit strongly and weakly correlated regimes and a number of phases and phase transitions that are very well understood and for which numerically exact comparison algorithms exist.Gull et al. (2011a) The restriction to these models allows us to provide stringent tests on the accuracy of SEET in a well controlled test environment. In the future, we aim to apply SEET to realistic materials, where it has the potential to become a method complementary to GW+DMFT or DFT+DMFT.

In order to generate a wide range of correlated phases, we generate our impurity parameters from a four-site dynamical cluster approximation (DCA)Hettler et al. (2000); Maier et al. (2005) to the 2D Hubbard model, i.e. test our method as a ‘DCA impurity solver’, so that our results can be compared against numerically exact continuous time CT-QMC data.Gull et al. (2008, 2011b) We emphasize that while we have developed a numerically efficient impurity solver, we do not envisage this as the main use of SEET, and we mainly resort to impurity models for the sake of generating comparisons to reliable known results.

SEET in the ED-in-GF2 variant is computationally affordable, as GF2 scales as , where is the number of orbitals in a unit cell. Additionally, SEET is amenable to parallelization on large machines. The scaling of GF2 can further be reduced to by employing density fitted integrals,Phillips and Zgid (2014) and the strongly correlated orbitals can be treated by ED as pairs, further reducing the numerical cost. Consequently, large systems containing many unit cells or -points containing multiple orbitals can be treated simultaneously, providing non-local effects and momentum dependence.

In Section II, we introduce SEET(ED-in-GF2). Sec. III shows results for our test model, and Sec. IV contains conclusions of our work.

Ii The SEET(ED-in-GF2) method applied to an impurity model

We consider an impurity problem with impurity orbitals coupled to an infinitely many bath orbitals described by a general Hamiltonian

(1)

where and are material specific one- and two-body operators, is the hybridization strength, and is the -electron dispersion. The single-particle properties of this Hamiltonian are described by a non-interacting Matsubara Green’s function for -electrons

(2)

with encapsulating the properties of the -electrons and being the chemical potential. In SEET(ED-in-GF2), we obtain the interacting Green’s function of this -orbital impurity problem iteratively, starting from , by self-consistent second order perturbation theory (GF2),Phillips and Zgid (2014); Dahlen and van Leeuwen (2005)

(3)

and the corresponding GF2 Green’s function

(4)

Note that GF2 can be solved self-consistently and includes an exchange diagram important for describing systems with a localized electronic density, but does not include higher order RPA-like diagrams that are present, e.g., in GW. We then evaluate the one-body density matrix using the converged GF2 Green’s function and choose a set of orbitals corresponding to eigenvalues of the one-body density matrix which are significantly different from 0 or 2. These orbitals, which we will call ‘strongly correlated’, are used to build an -orbital impurity problem which is then solved with a method more accurate than GF2 to compute a self-energy. Here, we use EDCaffarel and Krauth (1994) to solve this impurity problem. The resulting ED self-energy is used to modify the GF2 self-energy and to obtain the total self-energy in the natural orbital basis as

(5)

The indices and run over all orbitals, while and run only over the strongly correlated orbitals. The total self-energy is schematically illustrated in Fig. 1.

As the correlated orbitals are chosen in the eigenbasis of the one-body density matrix, a transformation of the one-body and two-body integrals in this -orbital subspace to the eigenbasis is necessary. Note, that even for cases where model Hamiltonians with simplified (e.g. local or density-density) interaction structures are studied, this transformation generates general interactions . This -orbital impurity problem with non-local interaction is then treated by the ED solver requiring an additional bath discretization step which may introduce fitting errors. We emphasize that these are small for the cases studied here and that, in principle, any solver suitable to describe strong correlations and able to treat general multi-orbital interactions can be employed, including QMC solvers based on the hybridization expansion Werner and Millis (2006) which do not require a bath discretization step.

The ED-in-GF2 procedure is iterated, and the GF2 calculation updates for the orbitals since where

(6)

is responsible for removal of diagrams later included at the ED level. Subsequent ED calculation updates the strongly part of self-energy . The iterative updates stop when the total self-energy in Eq. 5 is converged to a predefined accuracy. We present a detailed algorithmic description of this framework in the supplementary material.

Figure 1: The total self-energy in natural orbital basis produced in the SEET with ED-in-GF2 scheme.

The algorithm, both in its general form and in the ED-in-GF2 variant, is based on a diagrammatic formulation in which a ‘double counting’Karolak et al. (2010) 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.Aryasetiawan et al. (2004); Vaugier et al. (2012); Sakuma and Aryasetiawan (2013)

Iii Results

We calibrate SEET(ED-in-GF2) for the dynamical cluster approximation (DCA) to the 2D Hubbard model,Hettler et al. (2000); Maier et al. (2005) and consequently the Hamiltonian from Eq. 1 is defined for describing nearest neighbor hopping only and exclusively on-site interactions. DCA provides the non-interacting Green’s function (in Eq. 2) which is then employed to obtain the GF2 self-energy from Eq. 3. Subsequently, we construct the one-body density matrix and choose a pair of two-site impurities to be treated by ED. The occupations of the four site cluster in natural orbitals are 2-, 1, 1, , where for most regimes is not a small number, thus the orbitals with occupations 2-, and are not any longer weakly correlated. This motivates us to choose two separate impurity problems with orbitals occupied as (1,1) and (2-,) and treat them as a pair of two-site impurities embedded into the GF2 description. This means that only the interactions between these two-site impurities are treated at the GF2 level. A schematic description of the DCA+ED-in-GF2 iterative scheme is shown in Fig. 2. SEET allows us to treat multiple embedded impurities which is computationally advantageous since in realistic cases the number of strongly correlated orbitals may be too large for current solvers such as ED or the hybridization expansion.

Figure 2: Schematic view of the DCA+ED-in-GF2 procedure used for treating the 4-site cluster DCA approximation to the 2D Hubbard model. Note that ‘strong’ denotes quantities in the correlated subspace, and ‘weak’ quantities defined on the entire system.

Testing ED-in-GF2 using SEET on the DCA approximation to the 2D Hubbard model provides a worst case scenario for a multi-scale embedding scheme, since in multiple regimes the 4-site cluster does not display a separation of energy scales or any ‘weakly’ and ‘strongly’ correlated orbitals as typically found in realistic materials. Rather, in the Mott regime of the 2D Hubbard model, all orbitals are strongly correlated, providing a stringent test of the SEET with ED-in-GF2 method.

Figure 3: The imaginary part of the on-site CT-QMC, GF2 and ED-in-GF2 self-energy obtained for a 4-site cluster of the half-filled 2D Hubbard model for various regimes with .

In Fig. 3, the imaginary part of the self-energy is plotted for the half-filled case. For weak coupling, i.e. GF2 recovers the QMC results well. While for ED-in-GF2 corrects the GF2 result only slightly, for , the improvement is more substantial. In this case, the ED-in-GF2 recovers QMC results and is a quantitative correction to the qualitatively correct GF2 curve.

As expected, in the Mott regime, and , GF2 fails to recover the self-energy even qualitatively. Note that an IPT-like fitting of the large- limit to the atomic limit would be possible for this particular example but not in general, as it requires the determination of the local physics at exponential (in ) cost. In the Mott regime, ED-in-GF2 recovers to a decent quantitative accuracy the QMC self-energy for both and .

In Figs. 4 and 5, we examine several interesting regimes at doping, where the system exhibits the behavior of a strongly correlated Fermi liquid. In these cases, we report real and imaginary parts of Green’s functions rather than self-energies, since a slight difference in chemical potentials between different methods results in a shift of the Hartree term in the large- limit.

The imaginary part of Green’s function shows a good quantitative agreement between CT-QMC and ED-in-GF2 for multiple U/t regimes.

Figure 4: The imaginary part of the on-site CT-QMC, GF2 and ED-in-GF2 Matsubara Green’s function obtained for a 4-site cluster of the 2D Hubbard model at doping, for and with .

The real part of Green’s function shows more differences than the imaginary part. In the weak coupling regime illustrated in Fig. 5, for and , all the QMC, GF2 and ED-in-GF2 real parts of Green’s functions are close. The and regimes are more correlated and GF2 yields a qualitatively incorrect result. ED-in-GF2 corrects this result and provides a quantitative agreement with CT-QMC.

Figure 5: The real part of the on-site CT-QMC, GF2 and ED-in-GF2 Matsubara Green’s function obtained for a 4-site cluster of the 2D Hubbard model at 10% doping, for and with .

Iv Conclusions

We introduced a general self-energy embedding theory (SEET) for correlated systems and performed a comparison of SEET(ED-in-GF2) on a strongly correlated system for which the exact solution is known, the 4 site cluster DCA approximation to the 2D Hubbard model. This model has a continuous dispersion and shows a range of correlated phases, thus providing us with a detailed assessment of strengths and weaknesses of our method. However, it does not illustrate the effect of non-local interactions. Since in multiple regimes a clear separation of energy scales is not present, this model provides a rigorous test for a multi-scale method. We were able to show that ED-in-GF2 provides accurate results for the 4 site Hubbard model in the weakly correlated, intermediately correlated, and strongly correlated regimes, at and away from half-filling. While the solution in the strongly correlated embedded subset of orbitals has exponential scaling in our case, the total self-energy for the strongly correlated orbitals can be assembled using solutions of multiple small impurity problems. The calculation of the properties of the weakly coupled orbitals with GF2 scales as , making SEET(ED-in-GF2) an ideal tool for the simulation of realistic materials. Extensions using other diagrammatic or correlated methods, such as SEET(QMC-in-GF2) or methods based on GW are straightforward.

In real materials, the number of weakly correlated orbitals in the unit cell is significantly larger than the number of strongly correlated orbitals, thus providing an ideal situation where many orbitals can be treated cheaply by GF2 while the number of orbitals treated by ED remains small. Moreover, the SEET(ED-in-GF2) hybrid is easy to implement and, since it does not use frequency dependent effective interactions, can be trivially extended to employ different solvers for the strongly correlated part, such as truncated CI variants with a suitably chosen active space or QMC hybridization expansions. Similarly, the weakly correlated part can be treated by different levels of perturbation theory or cheap truncated CI methods, instead of GF2. Our ED-in-GF2 method can be adjusted to yield more accurate results, either by increasing the order of the perturbative treatment (e.g. by employing FLEX or GW), or by increasing the number of orbitals treated by ED. These limits therefore provide a rigorous assessment of the convergence of the self-energies. Since a set of strongly correlated orbitals in SEET is chosen based on a unique mathematical criterion,the method has the potential for becoming a black box method for realistic correlated materials calculations.

Acknowledgements.
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 comments on the first draft of this manuscript.

References

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
16600
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description