Solving the multi-site and multi-orbital Dynamical Mean Field Theory using Density Matrix Renormalization

We implement an efficient numerical method to calculate response functions of complex impurities based on the Density Matrix Renormalization Group (DMRG) and use it as the impurity-solver of the Dynamical Mean Field Theory (DMFT). This method uses the correction vector to obtain precise Green's functions on the real frequency axis at zero temperature. By using a self-consistent bath configuration with very low entanglement, we take full advantage of the DMRG to calculate dynamical response functions paving the way to treat large effective impurities such as those corresponding to multi-orbital interacting models and multi-site or multi-momenta clusters. This method leads to reliable calculations of non-local self energies at arbitrary dopings and interactions and at any energy scale.


INTRODUCTION
Among the most intriguing problems in physics is the behaviour of strongly correlated materials which present emergent behavior such as high temperature superconductivity, ferroelectricity, magnetism and metal-insulator transitions. These systems have triggered a great deal of research and are still far from being understood. However, a complete theoretical understanding is still lacking due to the presence of strongly interacting local orbitals in these materials. Methods to calculate electronic structure of weakly correlated materials, such as the Density Functional Theory (DFT) (Hohenberg and Kohn, 1964) which use the local density approximation (LDA) (Jones and Gunnarsson, 1989) and other generalizations, are unable to describe accurately the strong electronic correlation case. Non-perturbative numerical methods are, thus, the only reliable approach.
To include correlations, the Dynamical Mean Field Theory (DMFT) was developed more than twenty years ago. Together with its sucessive improvements (Kotliar and Vollhardt, 2004;Georges et al., 1996;Kotliar et al., 2001a;Maier et al., 2005;Hettler et al., 1998;Sénéchal et al., 2000), these methods have led to more reliable results. The combination of the DMFT with LDA has allowed for band structure calculations of a large variety of correlated materials (for reviews see Refs. (Imada and Miyake, 2010;Held, 2007)), where the DMFT accounts mainly for local interactions (Anisimov et al., 1997;Lichtenstein and Katsnelson, 1998). A recent alternative proposal, the Density Matrix Embedding Theory, DMET, was developed, which relies on the embedding of the wave functions of a local cluster fragment (instead of the local Green functions) in a self-consistent finite environment (Knizia and Chan, 2012;Bulik et al., 2014).
The DMFT requires the calculation of an interacting quantum impurity for which the fermionic environment has to be determined self-consistently until convergence of the local Green functions and the local self-energies is reached. Therefore, the success and scope of the DMFT will depend on the existence of accurate methods to solve correlated and complex quantum impurities. This approach is exact for the infinitely coordinated system (infinite dimensions), the non-interacting model and in the atomic limit.
Several quantum impurity solvers have been proposed since the development of the DMFT, among which we can mention the iterative perturbation theory (IPT) (Georges and Kotliar, 1992;Rozenberg et al., 1994), exact diagonalization (ED) (Caffarel and Krauth, 1994;Lu et al., 2014), the Hirsch-Fye quantum Monte Carlo (HFQMC) (Hirsch and Fye, 1986), the continuous time quantum Monte Carlo (CTQMC) (Rubtsov et al., 2005;Werner et al., 2006;Park et al., 2008;Gull et al., 2011;Nomura et al., 2015), non-crossing approximations (NCA) (Pruschke et al., 1993), the numerical renormalization group (NRG) (Wilson, 1975;Bulla, 1999;Stadler et al., 2015Stadler et al., , 2016, the rotationally invariant slave-boson mean-field theory (RISB) (Lechermann et al., 2007;Isidori and Capone, 2009;Ferrero et al., 2009) and quantum chemistry-based techniques (Zgid and Chan, 2011). Although these methods allow for the calculation of relevant properties such as the metal-insulator transition and other low-lying energy properties, they present some problems. Among them, one can mention the sign problem and the difficulty in reaching low temperatures in the QMC-based algorithms, the failure of the NCA in obtaining a reliable solution for the metallic state, the limitation to few lattice sites, far from the thermodynamic limit of the ED and the reduced high-energy resolution of the NRG technique.
To overcome some of these difficulties an impurity solver based on the Densit Matrix Renormalization Group (DMRG) technique (White, 1992;I. Peschel and Hallberg, 1999;Schollwöck, 2005;Hallberg, 2006;Nishimoto et al., 2004) was proposed (García et al., 2004(García et al., , 2007Fernández et al., 2014;Karski et al., 2005). Subsequent improvements to this were introduced, such as those using the time evolution DMRG algorithm (Ganahl et al., 2015;Bauernfeind et al., 2017), dynamical calculations using the Kernel Polynomial Method (Chebyshev expansion for Green functions) (Weiße et al., 2006;Holzner et al., 2011;Wolf et al., 2014a;Ganahl et al., 2014) and the application to non-equilibrium DMFT using MPS (Wolf et al., 2014b). It was recently realized that converging the DMFT loop on the the imaginary-frequency axis rather than the real-frequency axis reduces computational costs by orders of magnitude because the bath can be represented in a controlled way with far fewer bath sites and, crucially, the imaginary-time evolution does not create entanglement. The imaginary time setup can therefore treat much more sophisticated model Hamiltonians, opening the possibility of studying more complicated and realistic models and performing cluster dynamical mean field calculations for multiorbital situations. The price to be paid, however, is a reduced resolution on the real-frequency axis (Wolf et al., 2015a).
In spite of these developments, several difficulties still remain which hinder the calculation of reliable spectral densities for complex multi-band and multi-orbital correlated systems (Hallberg et al., 2015). In this paper we present a novel technique based on the DMRG which includes important improvements and complements previous methods. It is based on an efficient selection of the relevant states due to low entanglement bath configurations and on the targetting of the correction vector for small real energy windows. This method, thus, provides detailed spectral functions for complex Hamiltonians at zero temperature and for any doping and correlations. In the following sections we describe the method and show some applications and potential uses.

GENERAL FORMULATION
In order to present a unified treatment of multi-site (or cluster) and multi-orbital Hamiltonians on the lattice, we start by interpreting the lattice as a superlattice such that: 1. The interactionV is local and completely contained in the unit cell:V = iVi , where i is the cell index.
2. The non-interacting HamiltonianĤ 0 is characterized by its local Green function matrix G 0 (ω 1 − T ); being T = (t IJ ) the coefficients of the local partĥ 0 creates an electron in cell i and local "orbital" I = 1, 2, .., N c with spin σ =↑, ↓.
These two points completely define our problem through the parametersV i , G 0 , T . Notice that G 0 and T are typically well known one-particle quantities for a given lattice problem.
The key idea of the DMFT is to neglect the self-energy between different cells i and j in the lattice, that is, to consider only the local self-energy: Σ ij (ω) ≈ Σ(ω)δ ij . In this way, we are neglecting spatial correlations up to a certain degree while a good treatment of the local dynamical correlations is made. The relevant point is that the problem becomes tractable, as we will see in the following. Note that G 0 , T , and Σ are N c × N c matrices for the spin-symmetric solution, and 2N c × 2N c matrices in the general case. Spatial correlations or the momentum dependence of Σ can be obtain by periodization (Stanescu and Kotliar, 2006).
The local Green function is now given by (dmf, 2014) which defines the self-consistency condition for the N c × N c matrices G and Σ. The lattice problem can now be mapped onto an auxiliar impurity problem that has the same local magnitudes G(ω) and Σ(ω). This impurity problem should be determined iteratively. The impurity Hamiltonian can be written: where the non-interacting part H b represents the bath: b † Iqσ corresponds to the creation operator for the bath-site q, associated to the "orbital" I and spin σ (see Fig. 1).
The self-consistent iterations can be summarized as follows: (i) Start with Σ(ω) = 0, (ii) Calculate the Green's function: (iii) Obtain the hybridization (iv) Find a Hamiltonian representation H imp with hybridizationΓ(ω) to approximate Γ(ω). The hybridizationΓ(ω) is characterized by the parameters Υ q = υ IJ q and Λ q = λ IJ q of H b through: (v) Calculate the impurity Green's function matrix G imp (ω) of the Hamiltonian H imp using DMRG. (vi) Obtain the self-energy Return to (ii) until convergence.
Step (iv) is a fitting problem for Υ q and Λ q , where we can use the general symmetries of the hybridization function matrix. If Γ can be diagonalized using the same unitary rotation R for all ω then we obtain (at most) N c independent fittings. This can be seen from Eq. (6) after applying R: where the superscript D is used to stress that these matrices are diagonal, and In this new basis (the so-called molecular-orbital basis), we have to fit Γ D 11 (ω) using the expression (9) forΓ D 11 (ω) which depends on the parameters Υ D q 11 and Λ D q 11 , and similarly for Γ D 22 (ω), etc. Once these independent fittings are done, we bring the parameters back to our original basis through M = R · M D · R † .
In general, symmetries can be expoited for a better performance and stability. For example, at half-filling we could also have the electron-hole symmetry, giving a conection between G(−ω) and G(ω), implying the same structure for the hybridization Γ(ω).
The most resource-demanding part of the algorithm is carried out at step (v), where the dynamics of a complex many-body problem (see Figure 1) is calculated. Here we use the correction-vector method together with the DMRG essentially following (Kühner and White, 1999;Ramasesha et al., 1997), although other methods to calculate dynamical response functions withing the DMRG can also be used (Hallberg, 1995;Jeckelmann, 2002). The one-dimensional representation of the problem (needed for a DMRG calculation) is shown in the figure, where we are also duplicating the graph when considering spin degress of freedom (not shown for clarity). In this configuration (star geometry), in spite of the high connectivity of the Hamiltonian, the DMRG shows a much better performance (Holzner et al., 2010;Wolf et al., 2015b) The correction vector method is implemented in DMRG by targeting not only the ground state |E 0 of the system but also the correction vector |CV r associated to the applied operator at frequency ω r (and its neighborhood). For example, to obtain the single-particle density of states (DOS), the correction vector reads: (ω r + iη − H imp − E 0 ) |CV r = c † 0Iσ |E 0 , where a Lorentzian broadering η was introduced to deal with the poles of a finite-length impurity model. In this way a suitable renormalized representation of the operators is obtained to calculate the properties (2) corresponding to the impurity problem for the one, two and four-site cellular DMFT. The circles (squares) represent the non-interacting (interacting/impurity) sites. The red lines correspond to the λ IJ q parameters between bath sites q related to impurities I and J (they are the only hybridization between the baths related to different impuritues). The blue lines are the υ IJ q with I = J while the black lines are the υ II q . In the bottom scheme we omit some obvious connections for clarity. of the excitations around ω r , particularly the Green's function, for instance, Here ω r with r = 1, 2, ..., N ω is a grid covering the frequencies of interest, typically N ω =40-50 and are treated independently. Thus each DMFT iteration uses around 30 cores totalling less than three hours for all cases considered in this work, considering system sizes of up to 36 sites.

ONE-SITE DMFT
As we remark, only three parameters should be defined in order to apply the DMFT algorithm:V i , G 0 , T . We study the paramagnetic solution of the DMFT in the square (and Bethe) lattice using the following: where (k) = −2t (cos k x + cos k y ) − 4t cos k x cosk y , with k = (k x , k y ) the Fourier space of the square lattice with N sites, N → ∞, and t ( t ) denotes the (next-)nearest-neighbor hopping integral (Sakai et al., 2010).

TWO-BAND BETHE LATTICE
We consider the interaction: where J > 0 is the Hund exchange, U (U 2 ) is the intra (inter)-orbital Coulomb repulsion, and I = 1, 2 are the orbitals. The on-site non-interacting coefficients are T = −µ t 12 t 12 −µ , and the local Green's function: where B = t 1 0 0 t 2 , and t 1 , t 2 are the nearest-neighbor hoppings for each orbital.
Concerning step (iv), if t 12 = 0 then all our 2 × 2 matrices are diagonal and we have only to calculate two Green's functions and do two independent fittings, one for each orbital. On the other hand, if t 1 = t 2 but t 12 = 0 then we can introduce the rotation R = 1 √ 2 1 1 1 −1 to diagonalize the hybridization and we do again only two independent fittings. In the general case, a non-diagonal matrix fitting should be done to obtain a bath representation of the given hybridization Γ(ω), that is, to find the parameters Υ q and Λ q which minimize ω Γ(ω) −Γ(ω) 2 using, for instance, the matrix norm M 2 = T r M T · M .
In Fig. 2 we present the results for this model where, by analyzing the DOS for the different bandwidth case, the orbital-selective Mott transition can be clearly observed for a finite Hund's coupling J. This phase is robust for a certain range of interband hybridization. Previous calculations (Koga et al., 2004(Koga et al., , 2005Winograd and de' Medici, 2014) either resorted to approximate analytic continuation methods to obtain the DOS or used exact diagonalization for small baths. The method presented here leads to much more precise results and has the potentiality of treating even larger clusters or more orbitals. For example, it was crucial to find the in-gap holon-doublon quasiparticle peaks in the DOS of the asymmetric Hubbard model (Núñez-Fernández et al., 2017), which would have been hindered using QMC or NRG solvers or which would have lacked a proper finite size analysis had an ED method been used.

CELLULAR DMFT ON THE SQUARE LATTICE
We consider the same physical problem of section 3 on the square lattice, but interpreted now in a superlattice of unit cell of size N c = 2 (or 4) corresponding to the two-site (c2) or four-site (c4) cellular DMFT (Liebsch et al., 2008;Kotliar et al., 2001b). This case is illustrated in Fig. (1). The next-nearestneighbor-hoppings t for the c4-DMFT connect the opposite vertices of the 4 impurity square depicted at the bottom of this figure. Our three parametersV i , T , G 0 are now: respectively. Here, T is the non-interacting intracluster matrix and˜ (k) is the intercluster hopping on the superlattice Fourier spacek, which is connected to the one-site lattice through˜ (k) IJ = 1 Nc K exp i(K +k) · R IJ (K +k) with K the intracluster Fourier-space vectors, see eq. (23) of (Maier et al., 2005).˜ for the c2-DMFT and for the c4-DMFT. Finally, the hybridization matrix Γ has the following form, see (Liebsch et al., 2008): which can be diagonalized using the corresponding unitary rotation R, obtaining (at most) N c independent fittings.
In Fig. 3 we show the DOS for the Hubbard Hamiltonian on the square lattice with nearest (t = 0.25) and next nearest-neighbor hopping (t = 0) for two values of U . Larger clusters lead to a smaller critical U 's (Imada and Miyake, 2010). To illustrate the results with finite doping, in Fig. 4 we show the DOS for the Hubbard Hamiltonian on the square lattice with nearest (t) and a finite next nearest-neighbor hopping (t ), together with the non-local Green's functions.

CONCLUSIONS
We have presented here an efficient and reliable numerical method to calculate dynamical properties of complex impurities based on the DMRG. This technique uses the correction vector to obtain precise Green's functions on the real frequency axis directly thus avoiding ill-posed analytic continuation methods from the Matsubara frequencies and fermionic sign problems present in quantum Monte Carlo-based techniques, allowing also for zero temperature calculations. When used as the impurity-solver of the DMFT algorithm it leads to highly reliable spectral functions by using a self-consistent bath with low entanglement for which the density matrix renormalization works best.
To illustrate the versatility of the method, we have shown examples of densities of state and response functions within the DMFT framework for two paradigmatic models such as the Hubbard model at half filling on the square lattice on the one, two and four-site effective impurity models and at finite doping on the four-site case and also for the two-band Kanamori-Hubbard model on the Bethe lattice in the presence of Hund's coupling and interband hybridization.
This method leads to reliable results for non-local self energies at arbitrary dopings, hybridizations and interactions, at any energy scale. It also paves the way to treating large effective impurities not only within the framework of the DMFT to study multi-band interacting models and multi-site or multi-momenta clusters, but also for complex impurity problems such as adsorbed atoms, cold atoms and interacting nanoscopic systems like quantum dot arrays among others.
There is room to include additional improvements such as the consideration of symmetries, finite temperature, and more realistic systems by taking into account configurations given by ab-initio methods.