The two orbital Hubbard model in a square lattice: a DMFT + DMRG approach

We develop a precise and reliable numerical method for the calculation of electronic properties of the two orbital Hubbard model in a square lattice at half filling, based on the Dynamical Mean Field Theory (DMFT). We use the Density Matrix Renormalization Group (DMRG) as the impurity solver for the DMFT's selfconsistent equations to obtain accurate values for the Green's functions on the real axis. This way, reliable densities of states are obtained that do not need to resort to analytical continuation methods as those using quantum Monte Carlo techniques. Large system sizes can be achieved with increasing accuracy.


Introduction
The Dynamical Mean Field Theory (DMFT) has become one of the basic methods to calculate realistic electronic band structure in strongly correlated systems [1,2]. A key point of the DMFT method is the solution of an associated quantum impurity model where the environment of the impurity has to be determined self-consistently. Therefore the possibility to obtain reliable DMFT solutions of lattice model Hamiltonians relies directly on the ability to solve quantum impurity models. Among the available numerical algorithms there is the Hirsch-Fye [3,4] and Continues Time [5] Quantum Monte Carlo [3,4] method and Wilson's Numerical Renormalization Group (NRG) [6,7,8]. While the former, a finite-temperature method, is very stable and accurate at the Matsubara frequencies, its main drawback is the access to real frequency quantities for the calculation of spectral functions which requires less controlled techniques for the analytic continuation of the Green functions. It also leads to sign problems for the case of two or more bands with interband hopping. The second method, the NRG, can be formulated both at T =0 and finite (small) T and provides extremely accurate results at very small frequencies, offering a less accurate description of higher energy features.
In order to overcome the difficulties encountered by these other methods, several other techniques have been proposed. We have shown that the Density Matrix Renormalization Group (DMRG) [9,10,11,12] can be used very reliably to solve the related impurity problem within DMFT [13,14]. By using the DMRG to solve the related impurity problem, no a priori approximations are made and the method provides equally reliable solutions for both gapless and gapfull phases. More significantly, it provides accurate estimates for the distributions of spectral intensities of high frequency features such as the Hubbard bands, that are of main relevance for analysis of x-ray photoemission and optical conductivity experiments. Other calculations using alternative methods for the calculation of dynamical properties of the impurity have been proposed [15,16,17] and more recenty methods using the time evolution DMRG algorithm results for the one and two orbital model where shown on [18] and a restricted active basis set based on natural orbitals [19] 2. The Model In this paper we extend the method described above (that is, the DMRG used as the impuritysolver of the DMFT self-consistent equations), to the two-orbital Hubbard model on a square lattice. The model we are considering is [20]: where i, j are the sites of a square lattice and brackets indicate nearest neighbors, m indicates each of the two orbitals and σ is the spin of the electron, whose creation and destruction operators are c † and c, respectively.
Defining n i and s i as the on-site charge and spin operators respectively, the rotationally invariant on-site Hamiltonian is: For two orbitals per site, we have: So the local Hamiltonian is: Considering a ferromagnetic J (J< 0) and a triplet ground state which leads to µ = 2U , the total Hamiltonian reads: (1) Applying DMFT to this model leads to a mapping of the original lattice model onto an associated quantum impurity problem in a self-consistent bath. In the particular case of the two orbital Hubbard model, the associated impurity problem is the single impurity Anderson model (SIAM) with two levels, where the hybridization function ∆(ω), which in the usual SIAM is a flat density of states of the conduction electrons, is now to be determined self-consistently.
Starting from a guessed hybridization ∆(ω +iη) for the impurity, its Green function G(ω +iη) is obtained using the correction vector method of DMRG [21,22]. From this we can compute the self energy Σ(ω + iη) = G −1 − g −1 0 where g 0 is the non interacting Green function corresponding to ∆(ω + iη). The self energy allows us to compute the Green function on a lattice, in this case on a square lattice (SL): where ω sl (k x , k y ) = t(cos(k x ) + cos(k y )) with t = 1/2 to have a band of half-width D = 1. All energies are given in units of D. The lattice Green function G SL defines a new non interacting Green function g −1 0 = G −1 SL + Σ which in turn defines a new hybridization t 2 ∆(w) = ω + iη − g −1 0 (ω + iη) which is the seed to restart the cicle. The procedure is repeated until converged lattice or impurity Green functions are obtained (typically between 5 to 10 iterations). [23] To implement the algorithm we consider [24,25] a general representation of the hybridization function in terms of continued fractions that define a parametrization of ∆(ω + iη) in terms of a set of real and positive coefficients. Since it is essentially a Green's function, ∆(z) can be decomposed into "particle" and "hole" contributions as ∆(z) = ∆ > (z) c|gs for a given non interacting Hamiltonian, H with ground-state energy E 0 . By standard Lanczos technique, H can be in principle tri-diagonalized and the functions ∆ > (z) and ∆ < (z) can be expressed in terms of respective continued fractions [26]. As first implemented in Ref. [24,25], each continued fraction can be represented by a chain of auxiliary atomic sites whose energies and hopping amplitudes are given by the continued fraction diagonal and off-diagonal coefficients respectively. As a result of the self-consistency condition, the two chains representing the hybridization, are "attached" to the right and left of each of both levels of the impurity to obtain a new SIAM Hamiltonian, H. In fact G(ω + iη) constitutes the local Green's function of the site plus the chain system.
The SIAM Hamiltonian therefore reads with a two-level impurity, each of which is coupled to 2N C sites (N C to the right and similarly to the left). The set of parameters {a α , b α } are directly obtained from the coefficients of the continued fraction representations of ∆(z) by the procedure just described.

Results
In Figs.1, 2 and 3 we show the DMFT+DMRG results for the density of states (DOS) and the real and imaginary parts of the self-energy, respectively, for two values of the interaction U and a finite ferromagnetic J. We have kept around 700 states in the DMRG decimation procedure and our results were benchmarked with exact diagonalization calculations for smaller systems. These figures show the insulating (large U ) and the metallic (small U ) regimes and the lower and upper Hubbard bands in both cases. We are currently analyzing the whole parameter space to obtain the coexistence regime and its dependence with the Hund coupling J.
As a conclusion, we have developed an accurate real-frequency axis method for calculating electronic structure properties of the two-level Hubbard model. It relies on using the Density Matrix Renormalization Group as the impurity-solver of the Dynamical Mean Field Theory selfconsistent equations. With this method, large systems can be analyzed, thus reducing finite-size effects. It also deals with all energy scales on equal footing which allows to find interesting substructure in the Hubbard bands of the correlated metallic state with the corresponding effect on other experimentally accessible quantities such as the optical conductivity and high resolution photoemission spectroscopies. In addition, realistic band-structure calculations of systems with a larger number of degrees of freedom can be handled.