Systematically improvable multi-scale solver for correlated electron systems
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.
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
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
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)
and the corresponding GF2 Green’s function
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
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
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.
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)
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.
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.
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.
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.
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.
- Davidson (1975) E. R. Davidson, Journal of Computational Physics 17, 87 (1975).
- Lin and Gubernatis (1993) H. Lin and J. Gubernatis, Computers in Physics 7, 400 (1993).
- Liebsch and Ishida (2012) A. Liebsch and H. Ishida, Journal of Physics: Condensed Matter 24, 053201 (2012).
- Capone et al. (2007) M. Capone, L. de’ Medici, and A. Georges, Phys. Rev. B 76, 245116 (2007).
- Zgid and Chan (2011) D. Zgid and G. K.-L. Chan, J. Chem. Phys. 134, 094115 (2011).
- Zgid et al. (2012) D. Zgid, E. Gull, and G. K.-L. Chan, Phys. Rev. B 86, 165128 (2012).
- White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- Maier et al. (2005) T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
- Blankenbecler et al. (1981) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
- Parr and Yang (1989) R. Parr and W. Yang, Density Functional Theory of Atoms and Molecules (Oxford Science, 1989).
- Hedin (1965) L. Hedin, Phys. Rev. 139, A796 (1965).
- Bohm and Pines (1951) D. Bohm and D. Pines, Phys. Rev. 82, 625 (1951).
- Møller and Plesset (1934) C. Møller and M. S. Plesset, Phys. Rev. 46, 618 (1934).
- Phillips and Zgid (2014) J. J. Phillips and D. Zgid, J. Chem. Phys. 140, 241101 (2014).
- Dahlen and van Leeuwen (2005) N. E. Dahlen and R. van Leeuwen, J. Chem. Phys. 122, 164102 (2005).
- Foulkes et al. (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- Biermann et al. (2003) S. Biermann, F. Aryasetiawan, and A. Georges, Phys. Rev. Lett. 90, 086402 (2003).
- Biermann et al. (2005) S. Biermann, F. Aryasetiawan, and A. Georges, in Physics of Spin in Solids: Materials, Methods and Applications, NATO Science Series II: Mathematics, Physics and Chemistry, Vol. 156, edited by S. Halilov (Springer Netherlands, 2005) pp. 43–65.
- Karlsson (2005) K. Karlsson, Journal of Physics: Condensed Matter 17, 7573 (2005).
- Kotliar et al. (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
- Held (2007) K. Held, Adv. Phys. 56, 829 (2007).
- Tomczak et al. (2012) J. M. Tomczak, M. Casula, T. Miyake, F. Aryasetiawan, and S. Biermann, EPL (Europhysics Letters) 100, 67001 (2012).
- Casula et al. (2012a) M. Casula, P. Werner, L. Vaugier, F. Aryasetiawan, T. Miyake, A. J. Millis, and S. Biermann, Phys. Rev. Lett. 109, 126408 (2012a).
- Casula et al. (2012b) M. Casula, A. Rubtsov, and S. Biermann, Phys. Rev. B 85, 035115 (2012b).
- Ayral et al. (2012) T. Ayral, P. Werner, and S. Biermann, Phys. Rev. Lett. 109, 226401 (2012).
- Taranto et al. (2013) C. Taranto, M. Kaltak, N. Parragh, G. Sangiovanni, G. Kresse, A. Toschi, and K. Held, Phys. Rev. B 88, 165119 (2013).
- Sakuma et al. (2013) R. Sakuma, P. Werner, and F. Aryasetiawan, Phys. Rev. B 88, 235110 (2013).
- Hansmann et al. (2013) P. Hansmann, T. Ayral, L. Vaugier, P. Werner, and S. Biermann, Phys. Rev. Lett. 110, 166401 (2013).
- Ayral et al. (2013) T. Ayral, S. Biermann, and P. Werner, Phys. Rev. B 87, 125149 (2013).
- Biermann (2014) S. Biermann, Journal of Physics: Condensed Matter 26, 173202 (2014).
- Gull et al. (2011a) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011a).
- Hettler et al. (2000) M. H. Hettler, M. Mukherjee, M. Jarrell, and H. R. Krishnamurthy, Phys. Rev. B 61, 12739 (2000).
- Gull et al. (2008) E. Gull, P. Werner, O. Parcollet, and M. Troyer, EPL (Europhysics Letters) 82, 57003 (2008).
- Gull et al. (2011b) E. Gull, P. Staar, S. Fuchs, P. Nukala, M. S. Summers, T. Pruschke, T. C. Schulthess, and T. Maier, Phys. Rev. B 83, 075122 (2011b).
- Caffarel and Krauth (1994) M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (1994).
- Werner and Millis (2006) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006).
- Karolak et al. (2010) M. Karolak, G. Ulm, T. Wehling, V. Mazurenko, A. Poteryaev, and A. Lichtenstein, Journal of Electron Spectroscopy and Related Phenomena 181, 11 (2010), proceedings of International Workshop on Strong Correlations and Angle-Resolved Photoemission Spectroscopy 2009.
- Aryasetiawan et al. (2004) F. Aryasetiawan, M. Imada, A. Georges, G. Kotliar, S. Biermann, and A. I. Lichtenstein, Phys. Rev. B 70, 195104 (2004).
- Vaugier et al. (2012) L. Vaugier, H. Jiang, and S. Biermann, Phys. Rev. B 86, 165105 (2012).
- Sakuma and Aryasetiawan (2013) R. Sakuma and F. Aryasetiawan, Phys. Rev. B 87, 165118 (2013).