# A multi-orbital iterated perturbation theory for model Hamiltonians and real material-specific calculations of correlated systems

## Abstract

Perturbative schemes utilizing a spectral moment expansion are well known and extensively used for investigating the physics of model Hamiltonians and real material systems (in combination with density functional theory). However, such methods are not always reliable in various parameter regimes such as in the proximity of phase transitions or for strong couplings. Nevertheless, the advantages they offer, in terms of being computationally inexpensive, with real frequency output at zero and finite temperatures, compensate for their deficiencies and offer a quick, qualitative analysis of the system behavior. In this work, we have developed such a method, that can be classified as a multi-orbital iterative perturbation theory (MO-IPT) to study N-fold degenerate and non degenerate Anderson impurity models. As applications of the solver, we have combined the method with dynamical mean field theory to explore lattice models like the single orbital Hubbard model, covalent band insulator and the multi-orbital Hubbard model for density-density type interactions in different parameter regimes. The Hund’s coupling effects in case of multiple orbitals is also studied. The limitations and quality of results are gauged through extensive comparison with data from the numerically exact continuous time quantum Monte Carlo method (CTQMC). In the case of the single orbital Hubbard model, covalent band insulators and non degenerate multi-orbital Hubbard models, we obtained an excellent agreement between the Matsubara self-energies of MO-IPT and CTQMC. But for the degenerate multi-orbital Hubbard model, we observe that the agreement with CTQMC results gets better as we move away from particle-hole symmetry. We have integrated MO-IPT with density functional theory based electronic structure methods to study real material systems. As a test case, we have studied the classic, strongly correlated electronic material, SrVO. A comparison of density of states and photo emission spectrum (PES) with results obtained from different impurity solvers and experiments yields good agreement.

###### pacs:

71.27.+a, 71.10.Fd, 71.10.-w, 71.30.+h, 71.15.Mb## I Introduction

The development of efficient methods to solve quantum impurity problems, especially those involving multiple orbitals, has been a significant research direction in the field of theoretical condensed matter physics. Subsequent to the development of the dynamical mean field theory (DMFT)RevModPhys.68.13, which is exact in the limit of infinite dimensions and an excellent local approximation in finite dimensions, the importance of obtaining reliable solutions to general quantum impurity problems has increased further.

Within the DMFT framework, a lattice model may be mapped onto a quantum impurity embedded in a self-consistently determined host. The impurity problem may then be solved by a variety of techniques including– numerically exact methods like quantum Monte Carlo (QMC), numerical renormalization group (NRG), exact diagonalization (ED) and density matrix renormalization group (DMRG) or semi-analytical methods like iterated perturbation theory (IPT), local moment approach (LMA), non-crossing approximation (NCA) and fluctuation exchange approximation (FLEX). Each method has its own advantages as well as pitfalls. For example, QMCGull is a numerically exact method, but is computationally expensive. It yields data on the Matsubara axis (or imaginary time) so to obtain dynamical quantities such as the density of states and transport quantities, analytic continuation of the data to real frequencies is essentialm_jarrell_96a, which is a mathematically ill-posed problem. Additionally, it is very difficult to access the low temperature region where statistical errors become important. As a real frequency method, NRGRevModPhys.80.395 can avoid the difficulties that arise from the need to carry out analytic continuation. However, the method becomes extremely cumbersome for more than one impurity or channel. NRG is better suited for low temperature studies. Recently,Pruschke NRG was applied to study degenerate multi-orbital lattice problems, but the non-degenerate case remains unexplored. EDPhysRevLett.72.1545 is also a real-frequency method, but one considers only a finite number of bath states, so the resulting energy spectrum is discrete, and the broadening procedure for obtaining continuous spectra is not free of ambiguities. Moreover, large systems or multi-orbitals are not accessible. DMRGPhysRevLett.93.246403 for the single site case has some numerical artifacts and its accuracy as an impurity solver is not entirely clearGull.

The semi-analytical methods are perturbation theory based solvers that attempt to capture the essential physics by constructing an ansatz for the single-particle quantities. The ansatz is based on satisfying various limits or conservation laws, and comprises diagrams up to a certain order or sums a specific class of processes to infinite order. The main advantages of these methods are that they are computationally less expensive than the numerically exact methods listed above, while also yielding real frequency data. However, semi-analytical methods are, by definition, approximate and need to be benchmarked against exact results to gauge their range of validity. For example, although NCAAndreas gives qualitatively correct results for temperatures higher than the Kondo temperature, spurious non-analyticity at the Fermi energy develops at lower temperaturesPhysRevB.53.1850. To recover the correct Fermi liquid behavior at low temperatures, one needs to consider a larger class of diagramsPhysRevB.64.155111. The FLEX approximation is conserving in the Baym-Kadanoff sense, but it does not have the correct strong coupling behavior. So when it is employed for the half-filled Hubbard model, strong coupling physics like the Mott transition is not capturedMorita2002547; RevModPhys.78.865. The FLEXHaule has been extended to study degenerate multi-orbital problems but the issues plaguing single-orbital problems remain. The LMA is a highlyLMA1; LMA2 accurate technique, that has been benchmarked extensively LMA3 with NRG, but the method has not been used to study lattice problems except the periodic Anderson modelPramod. Moreover, extensions to symmetry broken phases or multiple orbital problems remain to be carried out.

The IPT is a simple, second order perturbation theory based method and it has been used widely to solve impurityMartin2; Yosida1 and lattice problemsGeorge1 at zero as well as at finite temperature. In the IPT, a self-energy ansatz is constructed that interpolates between known limits (i.e., weak coupling, atomic and high frequency limits) which is why it is also called an interpolative approach. It is clear that even the single-orbital IPT is not free of ambiguities so different constraints or limits to construct the ansatz yield different results. Hence, an IPT for multi-orbital problems has been ‘synthesized’ in many different ways by various groupsKajueter2, and we discuss these variations next.

The IPT ansatz for the self-energy is based on a rational or continued fraction expansion of a specific subset of diagrams, and consists of a small number of free variables that are fixed by satisfying various limits, such as atomic and high frequency limits and conservation laws such as the Luttinger’s theorem. Such an interpolative approach was first initiated by Martín RoderoMartin1; Martin2 for the single impurity and periodic Anderson models. The approach used the second order self-energy as a building block and the pseudo-chemical potential , was fixed by assuming that the occupation of the non-interacting part of the Anderson impurity problem is equal to the lattice occupation . Soon after the development of DMFT, the single band Hubbard model was studied by Muller-Hartmannmuller using self-consistent perturbation theory . The self-consistent perturbation method is able to produce a coherence peak in the single particle spectral function; however, it fails to reproduce the high energy Hubbard bands. For the single impurity Anderson model (SIAM), Yosida and YamadaYosida1; Yosida2 demonstrated that perturbation theory in is quite well behaved for the symmetric case when expanding around the Hartree-Fock solution. Based on these findings, Georges and KotliarGeorge1; George2 introduced an impurity solver called iterative perturbation theory (IPT) to solve the single band Hubbard model within DMFT which is based on mean-field and is able to capture coherent and incoherent features of single particle spectrum quite successfully. This is one of the biggest advantages of theories based on rather than the fully dressed G. Another advantage they offer is that, these theories naturally avoid any form of two particle divergence and are therefore able to provide a reliable description also “beyond” the perturbative regimePhysRevLett.110.246405. This is very important, because the twoÂparticle divergencies, which might induce ambiguities in the numerical determination of the DMFT self-energy () within standard perturbation theory schemes, occur in a rather large portion of the phase diagram, including the metallic regime much before the Mott transitionPhysRevLett.114.156402.

Subsequently, Kajueter and KotliarKajueter1; PhysRevB.53.16214 proposed a modification to the IPT called modified iterative perturbation theory (MIPT). In addition to the usual constraints of IPT, the MIPT constrains the zero frequency behavior of the self energy by adding a pseudo chemical potential to the Hartree corrected bath propagators. This pseudo-chemical potential, , can be obtained in different ways so there is an ambiguity in the method. KajueterKajueter1 fixed this free parameter by satisfying the Friedel’s sum rule (equivalently Luttinger theorem), hence his method is called IPT-L. The Luttinger theorem and Friedel’s sum rule are valid only at zero temperature, hence for finite temperature calculations, KajueterKajueter2 used the same that was obtained at zero temperature.

To study spontaneous magnetism in the single band Hubbard model, Potthoff, Wegner and NoltingNolting1; Nolting2 improved MIPT further by taking into account the spectral moments up to third order and instead of fixing by using Luttinger theorem, they fixed it by the constraint. This method may be called IPT-. They also considered the simpler option, where lattice chemical potential is equal to the pseudo chemical potential . This is called IPT- and they bench-marked IPT-L and IPT- with IPT-. Recently, Arsenault, Sémon, and TremblayTremblay bench-marked IPT- with CTQMC and found the pathology in IPT- that, in the strong coupling regime, the method does not recover a Fermi-liquid for filling close to . They suggested a new method (IPT-D) to fix the through a double occupancy constraint. The range of schemes originating from the inherent ambiguities at the single-orbital level give an idea of the far larger range of approximations that can be built at the multi-orbital level. These schemes will be described next.

KajueterKajueter2 extended his single orbital perturbative scheme to the degenerate multi-orbital case. He used the coherent potential approximation (CPA) to calculate higher order correlation functions in the self energy. He showed, by benchmarking against ED, that the scheme provides reasonable results only if the total particle density per site is less than one. For fillings greater than one, his scheme produced a false double peaked structure at the Fermi level instead of a single resonance. The reason for such a spurious structure is that the high frequency tails in the continued fraction expansion can be systematically improved by considering poles involving higher-order correlations functions in the self-energy, but this in turn seriously degrades the low frequency behavior when the Luttinger’s theorem is attempted to be satisfied. To study quantum transport in mesoscopic systems such as multi-level quantum dots, Yeyati et al. Martin3 introduced an interpolative scheme based on IPT-. LiebschLiebsch applied an extension of IPT to study the orbital selective Mott-transition, using which he showed that inter-orbital Coulomb interactions gives rise to a single first-order transition rather than a sequence of orbital selective transitions. In Liebsch’s extension of IPT for the multi-orbital case, he chose the self energy to be the combination of Hartree term and second order pair-bubble diagram with interaction vertices between electrons in different orbitals on the impurity. Laad et al.Laad constructed an interpolative scheme for multi-orbitals that was used extensively to study real materials through the LDA+DMFT framework. In a similar context, Fujiwara et al.Fujiwara developed an interpolative approach for degenerate multi-orbitals. The novelty of their method was that they used ligand field theory in the atomic limit to find the higher-order correlation functions.

Although there exist a large range of schemes for extending IPT to the multi-orbital case, extensive benchmarking of any single method has not been carried out. Recently Savrasov et al. savrasov, and Oudovenko et al.Oudovenko developed an interpolative approach for degenerate multi-orbitals based on a simple rational form of the self-energy, where the unknown coefficients in the self-energy are determined using slave boson mean-field and Hubbard I approximations. In their Hirsch-Fye-QMC work on the SU(4) Hubbard model, they have observed a good agreement in the particle-hole asymmetric cases.

In the present work, we build upon the previous knowledge to develop an interpolative scheme for solving a general multi-orbital quantum impurity problem. Our scheme is also based on the second-order self-energy as a building block and we use the generic name for the method as simply multi-orbital iterative perturbation theory (MO-IPT). Our method has a single pseudo-chemical potential , that is found by satisfying the Luttinger’s theorem. We impose the correct high frequency and atomic limits to get the unknown coefficients in the self-energy ansatz. In the single orbital case, we find that MO-IPT recovers the usual MIPT self energy expression and for the degenerate multi-orbital case, our MO-IPT self-energy expression reduces to that of KajueterKajueter2. The main novelty lies in handling the high frequency poles in a systematic way. The method is general enough that it can be applied to study symmetry broken phases, Hund’s coupling (density-density type) and crystal field effects.

Since MO-IPT is a semi analytical method it needs to be bench-marked. Subsequent to the description of the method, we embark upon an extensive benchmarking of MO-IPT with numerically exact, hybridization expansion continuous time quantum Monte Carlo method (S-CTQMC)Gull2 as implemented in the ALPSBauer libraries and our own implementation of interaction expansion CTQMC (W-CTQMC). Our main conclusion is that the MO-IPT method works very well when used away from integer-fillings, even at reasonably strong coupling. Using MO-IPT, we have addressed issues disputed in the current literature of doped Mott insulatorsTremblay and covalent band insulatorswerner2. In the multi-orbital degenerate case, the method proposed by KajueterKajueter2 shows spurious features which restricted its applicability to fillings smaller than one. However, our approach circumvents all the above issues and moreover captures the filling dependent effect of the Hund’s couplingJanus in the low energy scale. Our study of crystal field effects (the non-degenerate case) is the first attempt to extend the ansatz beyond the degenerate case. In addition, we have shown that our method produces a good default model for the analytic continuation of CTQMC data using the maximum entropy method Mark3. We have also integrated the MO-IPT with material-specific, density functional theory based calculations, and thus tested it for a prototypical example of strongly correlated electronic system, SrVO. A rather good agreement is obtained when the MO-IPT photo-emission spectra (PES) is compared with experiments.

We have organized the paper as follows. In section II, we outline the formalism for MO-IPT. In section III, we discuss results when the MO-IPT is applied in the DMFT context for lattice problems. In section IV, we present our conclusions along with future directions and possible improvements.

## Ii Model and Formalism

The multi-orbital Hubbard model for density-density type interactions and for cubic environment in standard second quantization notation is given by

(1) |

where creates the electron at lattice site , in orbital with spin and annihilates the electron at site , in orbital with spin . We are mainly interested in the local single particle electron dynamics, which is given by the momentum sum of the lattice Green’s function

(2) |

Where = +i and is the convergence factor. Here comprises intra-unit-cell hybridization and inter-unit-cell hopping, namely

(3) | ||||

(4) | ||||

(5) |

where are orbital energies, T are intra-unit cell hybridization matrix elements, and is the dispersion of the lattice, that depends on its geometry. For example, in the case of a simple cubic lattice, assumes the form, .

Within DMFT, one can map the multi-orbital Hubbard model on to an auxiliary impurity problem with a self consistently determined bath. The Hamiltonian of the corresponding single impurity multi-orbital Anderson model, is expressed in standard notation as:

(6) |

Here and are impurity orbital indices including spin. The first term in the above equation represents the orbital energy; the second term is the hybridization between the impurity and the host conduction electrons, the third term represents the host kinetic energy and the final term is the local Coulomb repulsion between electrons at the impurity. The corresponding impurity Green’s function is given by,

(7) |

where . Here is the hybridization matrix or equivalently the self-consistently determined bath and is the impurity self-energy obtained by solving the impurity problem. The set of equations is closed by noting that, within DMFT, the lattice self-energy is momentum-independent and is the same as the impurity self-energy, i.e . The local Green’s function obtained in Eq. (2) is used for defining a new hybridization as

(8) |

Obtaining the self-energy however is the most challenging step, and we employ multi-orbital iterated perturbation theory to solve the multi-orbital Anderson model. The starting point, as usual, is an ansatz for the impurity self-energy, given byKajueter1

(9) |

The self-energy is thus restricted to being diagonal in the orbital basis. In the above ansatz, the first term is simply the Hartree energy and the second term contains the second order pair-bubble diagram of matrix size NN, where N is the number of orbitals. The second order pair-bubble diagram on the real frequency axis is given by

(10) |

where and is the Hartree corrected bath propagator, which is obtained from a Dyson like equation, and is given by

(11) |

The pseudo chemical potential, , is found at by satisfying the Luttinger’s theorem,

(12) |

At finite temperature, an ambiguity exists in the determination of the pseudo-chemical potential. We choose to use the determined at zero temperature for all finite temperatures. The chemical potential, , is found by fixing the total occupancy from the local Green’s function, , to be equal to the desired filling,

(13) |

where the trace is over the spin and orbital indices. The unknown coefficients from Eq. (9) are obtained in the standard way by satisfying the high frequency limit and the atomic limit respectively. The detailed procedure to derive and their expressions are discussed in the Appendix. These coefficients contain higher order correlation functions. The order of the correlation functions depends on number of poles in the self energy. For example a pole of order n involves (n+1) order correlation functions. For a two pole ansatz and involve two and three particle correlation functions. We calculate the two particle correlation function Himadri using the equation of motion method to obtain Fetter

(14) |

This single equation is not sufficient to find all the two-particle correlators. Hence as an approximation, we use the following:

(15) |

We calculate the three particle correlation function encountered in B approximately by decoupling it in terms of two and single particle correlation functions. In this work, we have ignored the three particle correlation function.

## Iii Results and Discussion

The formalism developed in the previous section is applied to a wide variety of correlated systems. We begin with a discussion of the well studied paramagnetic Mott transition in the half-filled single-band Hubbard model. Then we examine the doped Mott insulator. The covalent insulator is considered next, followed by the two-orbital Hubbard model. For the latter, We investigate the effect of filling, Hund’s coupling and crystal field splitting. Finally, we move to real material calculations considering specifically the case of SrVO. As mentioned earlier we bench-mark our results with those from numerically exact CTQMCGull2 methods. The CTQMC formalism yields results on the Matsubara frequency axis so to get the real frequency data, analytical continuation is required. We avoid analytic continuation by transforming the real frequency data obtained from MO-IPT to imaginary frequencies using the following spectral representations:

(16) |

and

(17) |

where and . In order to quantify the efficiency of the method, the imaginary part of the self energy needs to be bench-marked rather than the Green’s function. This is because the former is far more sensitive than the latter and moreover, the low energy scale of the system depends on the imaginary part of the self energy.

### iii.1 Single band Hubbard model: Half-filled case

The Hamiltonian for the single band Hubbard model is given by

(18) |

We study the above model within DMFT for a semi-elliptical density of states, given by

(19) |

Here W is the full-band width. In our calculations, we choose the energy unit to be = 1.

The half-filled Hubbard model exhibits an interaction-driven metal-insulator Mott transition at a critical . Terletska et al. PhysRevLett.107.026401 found that the critical exponents and scaling functions obtained by IPT are identical to those from CTQMC. Here, we revisit this case and benchmark the quasiparticle weight, double occupancy, spectra and imaginary part of the self-energy. The MO-IPT method reduces to the second order perturbation theory in terms of Hartree-corrected propagators. In Figure 1(a) we compare the quasi-particle weight obtained from different impurity solvers and several values of the Coulomb interaction. The values of obtained from S-CTQMC match well with those from NRGBulla1 for all values of U/W except close to the Mott-transition. This is most likely because we have done CTQMC calculations at , while NRG is at zero temperature. The critical interaction strength, obtained from both the methodsBulla2 agrees very well. The obtained from MO-IPT at matches quantitatively with CTQMC and NRG in the weak coupling limit and only qualitatively in the proximity of the transition. On the other hand, the results of the self energy functional approach (SFA)SFA agree with MO-IPT in the strong coupling limit rather than in the weak coupling limit. The MO-IPT yields the critical value of = 1.42, which is in good agreement with the critical value = 1.45 obtained from SFASFA at zero temperature. The double occupancy obtained from MO-IPT and S-CTQMC (shown in panel (b) of Fig. 1) also match, except very close to the transition. A detailed comparison of spectra from S-CTQMC and W-CTQMC with the same from MO-IPT (transformed to imaginary frequencies) is shown in Figure 2. The left panels show the imaginary part of the Green’s function at (top panel) and (bottom panel), while the right panels show the imaginary part of the corresponding self-energies. The excellent agreement between the three methods is clearly evident.

### iii.2 Single band Hubbard model: Doped Mott insulator case

The single band Hubbard model has gained a lot of interest, because the doped Mott insulator regime is believed to capture the essential physics of high T superconductorsRevModPhys.78.17. This regime is, in reality, highly complex, because many different factors such as proximity to the antiferromagnetic Mott insulator, disorder, d-wave superconducting fluctuations and pseudogap physics have to be treated on an equal footing. Hence, investigations of the doped Mott insulator in all its glory represents one of the toughest challenges in condensed matter. Here, we take a simplistic approach to the problem, and investigate the performance of MO-IPT in the paramagnetic doped Mott insulator in infinite dimensions. Our MO-IPT reduces basically to the IPT-L in this regime.

A comparison of quasi-particle weight at obtained from MO-IPT and S-CTQMC as a function of filling (Fig. 3) yields, surprisingly, an excellent agreement. We observe that as we decrease the filling (from 1) for a given , the Mott insulator turns into a strongly correlated metal and finally ends up as a simple metal. In the strong coupling limit, for filling close to , the IPT-n method gives an insulating solution, while the IPT-L correctly predicts a metal in agreement with exact methods. Kajueter and KotliarKajueter2 have benchmarked the real-frequency spectral functions obtained from IPT-L with exact diagonalization calculations and had found good agreement.

We find that the imaginary part of the Green’s function and self-energy obtained from IPT, when transformed to the Matsubara frequency axis using Eqs. (16), (17) are almost identical to those obtained from the strong coupling and weak-coupling variants of CTQMC (see Fig. 4). The slope of the as is , and the good agreement of shown in Fig. 3 is simply a reflection of the detailed agreement for all frequencies. Such an excellent agreement is truly surprising because IPT is a perturbative method by construction and the strongly correlated, doped Mott insulator regime should not, in general, be amenable to perturbative methods.

### iii.3 Covalent Insulator:

The discovery of topological insulators RevModPhys.82.3045 has led to a renewed interest in the role of e-e correlations in band insulators (BI)kunes. The prime examples of such materials would be FeSiFesi and FeSbFesb, since experimental measurements indicate a small optical gap and large thermopower (at low ). Increasing temperature leads to closing of the gap, and concomitantly a insulator-metal crossover in the resistivity. Such large scale spectral weight transfers are highly indicative of strong correlations. Specific heat measurements also seem to validate this observation. The band gap in these systems is a simple consequence of the structure of the hopping matrix and not of completely filled electronic shells werner2. Hence these materials are called covalent insulatorskunes; werner2. A Hamiltonian that describes the covalent insulator is given bywerner2

(20) |

where = a and b are two sub-lattices with semi-elliptic bands and having dispersion and - respectively. The two sub-lattices are coupled by a independent hybridization . While the unit of energy is chosen to be W = 2 throughout, for this subsection W=4 has been chosen in order to benchmark with earlier resultswerner2. This is the first two-band model we have studied in this work, since the previous cases pertained to the single-band Hubbard model. Hence this will be the first real test of the ‘multi-orbital’ part of MO-IPT. Since this is still the half-filled case, Luttinger’s theorem does not have to be satisfied explicitly. The and for all orbitals. Thus, the MO-IPT used for the covalent insulator case is equivalent to that employed by Liebsch Liebsch for studying the Mott transition in the two-band Hubbard model.

Before presenting the numerical results, we will present a simple analytical argument about the absence of any intermittent metallic phase in the CBI at zero temperature unlike the interaction-induced metallic phase in the Ionic Hubbard modelPhysRevLett.97.046403. The structure of in the CBI is given by,

(21) |

and the corresponding impurity Greens function for sublattice ’a’ is given by,

(22) |

where

and . In the half-filling case, the Hamiltonian has mirror type symmetry between sublattices, which reflects in the impurity Green’s function and self-energy in the following way,

(23) | ||||

(24) |

By using above self-energy symmetry relation, it is easy to see that

(25) |

then Eq. (22) can be written as,

(26) |

With the assumption that in the band insulator and metallic phases, a Fermi-liquid expansion of self-energy holds, namely that . Then, the value of imaginary part of self-energy at zero frequency (and zero temperature) is , and the corresponding density of states (DOS) is given by,

(27) |

where =. In the limit ,

Since the argument of the Dirac delta function is positive definite for any , the density of states at the chemical potential, will necessarily vanish for any , and hence the system will be gapped. Thus for the CBI, interactions do not close the gap, no matter how strong they are. This implies a clear absence of metallicity in these insulators.

The quasi-particle weights (Fig. 5(a)) and double occupancy (Fig. 5(b)) of sublattice ’a’ obtained from MO-IPT and S-CTQMC (shown as black circles and red squares respectively) are in close agreement except in the proximity of the transition of the correlated band insulator to a Mott insulator. Unlike the ionic Hubbard model casePhysRevB.89.165117, we do not see any intervening metallic phase between the correlated band insulator and the Mott insulator. This is also consistent with the S-CTQMC and analytical results. At high temperatures, the correlated band insulator should be gapless, and must develop the gap with decreasing temperature. Precisely this behavior is seen in the real frequency spectra (left panels, Fig. 6), which arises from the spectral weight transfer in the self-energy as a function of temperature. The high reliability of these spectra and self-energies computed through MO-IPT is apparent in the excellent agreement with the same obtained through strong coupling CTQMC (on the Matsubara axis, in Fig. 7). The crossover of the band-insulator to Mott insulator is also visible in the increasing (negative) slope of the imaginary part of self-energy with increasing . Similar kind of correlation induced transitions have also observed between topologically trivial and non-trivial band insulatorsTsui2015330; PhysRevLett.114.185701.

### iii.4 Two orbital Hubbard model:

Encouraged by the excellent benchmarking of MO-IPT with CTQMC for the two-band covalent insulator system, we now move on to the two-orbital Hubbard modelAntipov; Bierman. The Hamiltonian, in standard notation, for a cubic environment and for unbroken spin symmetry, is described in Eq. (1). Throughout the paper, we have considered local interactions of density-density type which are obtained by neglecting spin flip and pair-hopping terms that must be present for a rotationally invariant Hund’s coupling. The hopping is taken to be diagonal in orbital indices for simplicity.

#### (a) Half-filling: J = 0

We begin by considering the half-filled case (total occupancy is two) with . The Hamiltonian (Eq. (1)) has SU(4) symmetry in this situation. We have employed a semi-elliptic non-interacting density of states of full-band width W = 2 for the MO-IPT-DMFT calculations.

In Fig. 8a, we plot the quasi-particle weight () obtained from different impurity solvers for the particle-hole symmetric case. The results from strong coupling CTQMC, EDMedici and SFASFA, including the critical U, where the system transitions from metal to Mott-insulator, are in good agreement. The critical value U obtained in the multi-orbital case is greater than the value obtained in the single band case. The Mott transition is absent in the FLEX result Haule. The MO-IPT is seen to underestimate the quasiparticle weight as compared to the other methods (except MO-IPTnn=0; see below). However, the critical agrees reasonably well with that from hybridization expansion CTQMC. The green diamonds are from a variant of MO-IPT (used e.g. by Fujiwara et al. Fujiwara) where the two-particle correlation function is simply decoupled into two single-particle terms (=). The neglect of two particle correlations leads to a much worse comparison than MO-IPT. In contrast to the not-so-good agreement with exact methods for the quasiparticle weight, the average double occupancy obtained from MO-IPT shows excellent agreement with CTQMC (see Fig. 8b). Since the total energy of the system depends on single particle and two particle correlation functions, we expect that thermodynamic quantities like total energy or specific heat computed through MO-IPT might be reliable. One more important observation is that the double occupancy remains finite and almost constant even beyond the Mott transition, unlike the the single band case. We also compare the single-particle Green’s function and self-energy on the Matsubara frequency axis (Fig. 9). At high frequencies, the agreement between MO-IPT and S-CTQMC is seen to be excellent, while the agreement worsens at low frequencies, especially with increasing .

#### (b) Half-filling: Effect of Hund’s coupling ()

The interplay of Hund’s coupling, , and local interaction, , has been investigated by several groups. The main consensus is that strong correlation effects can be affected significantly through Medici; Janus. For example, in the half-filled case, the for Mott transition is lowered by , where is the number of orbitals, while the critical is enhanced by in the non-half-filled (but integral occupancy) case Medici.

It is important to know the extent to which the interplay between and is captured by the MO-IPT method. In Fig. 10(a) and (b), the quasi-particle weight for different values of J obtained from S-CTQMC and ED Medici is shown. Indeed, with increasing , the at which decreases sharply, as expected from the atomic limit. Also, for each , the quasiparticle weight decreases monotonically with increasing interaction strength. Although the latter trend is qualitatively captured by the MO-IPT result (shown in Fig. 10c) for larger , there is a disagreement with the exact results at lower values. The MO-IPT yields a that is a non-monotonic function of . The insets of panels a and b in Fig. 10 zoom in on the low interaction () part of the main panels. Unlike for , where increasing leads to a monotonic reduction of , a rise and fall of is observed for . Although such a trend is achieved by MO-IPT as well, the non-monotonicity sustains even for larger . A frozen local-moment phase is seen in the S-CTQMC calculations for any given in the strong coupling limit, while such a phase is not observed either by EDMedici or MO-IPT calculations. It must be mentioned here that the CTQMC calculations employ a density-density type Hund’s coupling, while the ED employs a fully rotationally invariant . Although the quasiparticle weight dependence on and is not accurately captured by MO-IPT, the single-particle dynamics on all scales is in qualitative agreement with S-CTQMC calculations (as seen in Fig. 11).

#### (c) Away from Half-filling: Effect of

The MO-IPT method works best away from half-filling, which is consistent with the results of comparisons carried out previously by other groupskotliar. In order to illustrate this, here we study the two orbital Hubbard model for a total occupation of . The imaginary part of the Matsubara self-energy obtained from S-CTQMC matches well with that from MO-IPT (Fig. 12, panels (a) and (b)), hence the latter does well in this regime. This observation is reinforced by the panels (c)-(e), which show a comparison of the quasiparticle weights as a function of for three values of , namely and . The results of MO-IPT are seen to agree very well with those from CTQMC. For most real material calculations, the regime considered in this subsection is perhaps the most relevant. Hence, accurate results from MO-IPT prove its efficacy for integration into first-principles approaches.

The Hund’s coupling and Coulomb interaction have a synergistic effect at half-filling, while in the doped case, the reverse occursJanus. This is shown in panel (f) of Fig. 12, where an increase of is seen with increasing Hund’s coupling at a fixed interaction strength.

It is quite instructive to study the real frequency spectral functions and self-energies as obtained from MO-IPT. These are shown in Fig. 13 for various values of interaction strength and Hund’s coupling, . In the absence of Hund’s coupling, the spectrum (shown in the left panels of Fig. 13) exhibits spectral weight transfers characteristic of increasing correlation strength: a central resonance that becomes sharper, and Hubbard bands that grow in prominence with increasing . However, at a fixed , increasing Hund’s coupling leads to a reversal of the aforementioned trend, i.e, a broadening of the resonance and a melting of the Hubbard band (see e.g. left panel bottom figure). In this parameter regime, a previous formulation of the multi-orbital iterated perturbation theorykotliar found a double peak structure at the chemical potential. Such a feature was shown by the authorsKajueter1 to be spurious by comparison to results from exact diagonalization. The reason we do not observe such a spurious feature is that we have considered only two poles in the self-energy, in contrast to the formulation of Ref. Kajueter1, where they have retained all the eight poles (for a two-orbital model). Although our ansatz seems like an ad-hoc truncation scheme, the justification for such a scheme lies in its excellent agreement with CTQMC results (shown in Fig. 14) and the absence of spurious features. In Fig. 14, the imaginary part of Matsubara Green’s functions and self energies obtained from MO-IPT are compared with those from CTQMC for three values of at and 64. For all values of the Hund’s coupling, an excellent agreement is obtained.

### iii.5 Two orbital Hubbard model: Crystal field splitting and Hund’s coupling

We now proceed to the case of a non-degenerate two-orbital model with crystal field splittingNondeg in the presence of Hund’s coupling. In most materials, the crystalline environment lifts the orbital degeneracypavarini1. For example in transition metal oxides, due to crystal field effects, the five fold degenerate -level splits into triply degenerate t and doubly degenerate e levels and the corresponding energy gap is 1-2 eV. The degeneracy of each of these levels (t, e) is further lifted by distortions such as the ones occurring in GdFeO, or arising through the Jahn-Teller effect or spin-orbit coupling. The energy cost for such distortion induced splitting is a few meV. Recently, Pavarini et al.pavarini2 studied crystal field effects in d type perovskites such as SrVO, CaVO, LaTiO and YTiO. It was found that crystal field effects and cation-covalency (GdFeO -type distortion) lift the orbital degeneracy and reduce the orbital fluctuations. Thus, investigating crystal field effects in model HamiltoniansPhysRevB.78.045115; PhysRevB.88.195116 is highly relevant for understanding of real materials.

We have investigated the Hamiltonian in Eq. (1) by considering two orbitals with energies and , which corresponds to a crystal field splitting of . The results from MO-IPT, for a fixed total filling of , are compared with those from strong coupling CTQMC at the corresponding orbital occupancies. In Fig. 15, we compare the quasi particle weights of the two orbitals obtained from MO-IPT with that of CTQMC. We observe a better agreement of for orbital-1 than for orbital-2. This must be expected, since orbital-1 is further away from particle-hole symmetry than orbital-2. The corresponding orbital occupancies as a function of increasing interaction (and hence ) are shown in Fig. 15. The deviation between results from the two methods increases with increasing and , which indicates that MO-IPT is almost exact for .

Next, we benchmark the single-particle dynamics in the presence of crystal field splitting. In Fig. 16, we show the imaginary part of the Matsubara frequency self energies obtained from MO-IPT and CTQMC for orbitals-1 and 2 (left and right panels respectively). The agreement between the results is quite evident, this suggesting that the MO-IPT should serve as a good method to study interacting, real material systems with finite crystal field effects and Hund’s couplings. This is especially true if the material in question has a large number of bands, which would make it prohibitively expensive to treat with CTQMC, while MO-IPT would be able to handle it with ease. We now demonstrate the efficacy of MO-IPT when applied to a well studied, real material system, namely SrVO.

### iii.6 Application to real materials: SrVO

Over the past decade or so, the combination of density functional theory (DFT) with dynamical mean field theory, such as LDA+DMFTHeld, has emerged as one of the most powerful methods for electronic structure calculations of strongly correlated electronic systems. Although the DFT results contain rich, material specific information, being a single particle theory, it works well only for weakly correlated systems where the ratio of Coulomb interaction () to bandwidth () is small i.e., . If we consider the opposite limit of , we have successful methods like the Hubbard-I and Hubbard-III approximations or the LDA+U method for predicting the ground state of the system. But these also have limitations, such as the neglect of dynamical fluctuations in the LDA+U method. In nature, there are many materials, for example transition metal oxides, which lie in between these two limits. It has been established in the context of model Hamiltonians that the DMFT can handle both limits quite efficiently. Hence a natural combination of LDA with DMFT is expected to bring predictive capabilities in the theory of strongly correlated electronic systems. Nevertheless, LDA+DMFT is not without its own bottlenecks.

One of the central issues of the LDA+DMFT method is the correct definition of a correlated subspace. The basic idea of a correlated subspace is to make an appropriate choice of energy window around the Fermi level and fit the band structure to a few-orbital tight-binding model. Many techniques have been proposed to construct such a material specific ‘non-interacting’ Hamiltonian. The two major techniques for this purpose are down-folding Andersen19951573 and projection based Wannier function technique PhysRevB.56.12847. In general, bands which are crossing the Fermi level are considered in the desired energy window for Hamiltonian construction. For example in transition metal compounds bands having d-orbital character normally cross the Fermi level. This process becomes simple if there is no hybridization in the system and these bands with d-orbital character are well separated from other bands like the p-bands. As Dang et al.PhysRevB.90.125114 pointed out, a mixing of these d orbital bands with p orbital bands can create several complications.

After getting the ’non-interacting’ Hamiltonian , one can add various types of interactions terms to obtain a full material-specific multi-orbital model. The solution of such a Hamiltonian is however a major challenge and this is where the MO-IPT can be most useful, since it scales only algebraically with increasing number of bands, while yielding real frequency quantities directly. In contrast, impurity solvers like CTQMC and ED scale exponentially with increasing number of orbitals and are very expensive, especially for investigations of real materials. As a test case, we study SrVO which is considered a prototypical example of a strongly correlated electronic system.

#### (a). Computational Details

We perform our density functional theory (DFT) calculations with linearized augmented plane wave (LAPW) based method as implemented in the all-electron package WIEN2KBlaha. The experimentally determined structureRey1990101 of cubic SrVO in a non-magnetic phase was used for the calculations (neglecting spin-orbit coupling). The product of the plane-wave cut off and the smallest atomic sphere radius was chosen as for controlling the basis set. The radii of the muffin-tin spheres were chosen to be larger than the corresponding atomic radii. Thus, the values used for were for Sr, for V and for O. With these parameters, charge leakage was absent and our DFT results agree well with DFT calculation using other basis sets 0953-8984-20-13-135227. We utilize the generalized gradient approximation (GGA) of Perdew, Burke and Ernzerhofperdew for the exchange and correlation functional. In this calculation, we consider 512 k-points in the irreducible part of the Brillouin zone. After getting the Bloch-eigen states, all the necessary inputs for constructing the maximally localized Wannier functions (MLWFs) are prepared by the WIEN2WANNIER codekune1. Finally, the Hamiltonian is constructed in the maximally localized Wannier basis by taking a projection of the three orbitals within the energy window of -1.0 eV to 1.8 eV with the standard procedure implemented in Wannier90Mosti. We begin by discussing the DFT results.

#### (b). GGA+DMFT: Results and discussion

Our computed band structure and density of states (DOS) are presented in Fig. 17 and Fig. 18. The three bands, crossing the Fermi level, are highlighted in cyan, violet and grey colors. These bands originate from the states, and are located between -1.1eV and 1.5eV. The states lie at higher energies, between 1.1eV to 5.8eV (see the projected density of states in Fig. 18). The band structure agrees well with previous results by Ishida et. al.PhysRevB.73.245421 obtained in the LAPW basis. When compared with results from the linear muffin-tin orbital (LMTO) calculations of Nekrasov et al. PhysRevB.72.155106, the position of bands agrees well but the position of states differs by about 0.3 eV. This discrepancy is, most likely, due to the difference in basis sets used in the two calculations. A significant computational simplification results from ignoring the hybridization between and orbitals, since the low energy correlated subspace comprises just the three orbitals.

Thus, the DFT results yield a ‘non-interacting’ Hamiltonian , which in this case is a matrix for each . Thus, the full Hamiltonian is given by

(28) |

where is the interaction term is given by

(29) |

In the above expression, stands for V sites and is the orbital index with spin . , and are the local, intra orbital and inter orbital Coulomb repulsion respectively and is the Hund’s exchange. The local, non-interacting lattice Green’s function, in the orbital basis, , can be obtained from the DFT calculated by the following equation as

(30) | ||||

(31) |

where is the chemical potential and is the hybridization. In the DFT approach electronic correlations are partially entered through the LDA/GGA exchange-correlation potential. This part of the interaction () has to be subtracted in the LDA+DMFT approach to avoid double-counting. This is not an important issue when the low energy effective Hamiltonian contains only the d-manifold because we can absorb it into the chemical potential. However it is an important issue when the low energy effective Hamiltonian contains O-2p orbitals also. Various schemes for finding the double-counting correction exist, each with a different physical motivation. Details about such schemes may be found in the work by Lechermann et al.PhysRevB.74.125120 and Nicolaus ParraghPhysRevB.86.155158; Parragh. In general we can construct the modified host Green’s function for the orbital as

(32) |

We find the pseudo-chemical potential using the same procedure as in the model calculations. The self-energy can be found, e.g. through the MO-IPT method outlined in Section II. The second-order self-energy in Eq. (9) is a functional of the modified host Green’s functions, . The full local Green’s function for the lattice Hamiltonian (Eq. (28)) is given by

(33) |

The above Green’s function may be used to obtain a new host Green’s function through the Dyson’s equation:

(34) |

In general, the chemical potential, is found by fixing the total occupancy from the full Green’s function, to be equal to the value found from DFT,

(35) |

where the trace is over spin and orbital indices.

Thus the full solution of the problem proceeds as follows. Given the , we guess an initial self-energy, as well as the and ; and use these to find the local and the host Green’s functions through Eqs. (33) and (34). The host Green’s functions are then used to find the self-energy, and Eqs.(33) and (35) are used to find the chemical potential. For a fixed , these equations are then iterated, until the self-energy converges. With the chosen pseudo-chemical potential, the Luttinger’s integral, Eq. (12) is computed using the converged self-energy and local Green’s functions. If the Luttinger’s theorem is satisfied within a numerical tolerance, the solution is considered to be obtained, else the is tuned, and the DMFT equations are iterated, until the Luttinger’s theorem is satisfied.

The DFT predicted occupancy per spin on the three correlated V-t orbitals in SrVO is 0.166, which implies SrVO is a d system. For the DMFT calculations, we employ interaction parameters eV and eV, that were obtained by Taranto et al. Taranto through the random phase approximation (RPA). For SrVO, we have not introduced explicit double counting correction because we choose the correlated subspace that is identical with the set of Wannier bands. We absorb the double counting correction and orbital energies in the lattice chemical potential, which we find by using Eq. (35).

Our computed GGA+DMFT spectrum for SrVO is shown in Fig. 19 and compared with results obtained from other impurity solvers. The GGA result (shown in blue) has no signatures of correlation, while each of the DMFT calculations exhibit a three peak structure. The CTQMC results from GW+DMFT (black) agree qualitatively with those from LDA+DMFT. However, the details do differ. Namely, the positions and weights of the resonance at the Fermi level and of the Hubbard bands differ to a significant extent. This difference, naturally, can be attributed to the different starting points, namely GW vs LDA, of the CTQMC calculations. Results from the MO-IPT solver agree with those from CTQMC in the neighborhood of the chemical potential as well as in the proximity of the lower Hubbard band. The upper Hubbard band is clearly in disagreement with the CTQMC results.

As a final benchmark of the GGA+DMFT(MO-IPT) calculation, we compare our result with the experimentally measured photo emission spectrum (PES) which is shown in Fig. 20. A Hubbard satellite at eV is seen in the experimental PES spectrum. Our GGA+DMFT(MO-IPT) calculation predicts the Hubbard satellite at -1.25 eV. Results from other approaches, namely LDA, LDA+DMFT(CTQMC) and GW+DMFT(CTQMC) are also reproduced. Surprisingly, the closest match with the experiment is achieved by the GGA+DMFT(MO-IPT) in terms of the position and width of the resonance at the Fermi level and of the lower Hubbard band. Thus, we infer that the MO-IPT method outlined in this work may be used as an efficient tool to study the electronic structure of real material systems.

## Iv Conclusions

The development of iterated perturbation theory as an impurity solver for single band models and for multi-band models dates back to almost two decades. Although a few comparisons with numerically exact methods have been made, being a perturbative approach, the method has suffered from reliability issues, especially for multi-orbital systems. Nevertheless, several multi-orbital extensions of IPT have been proposed and used to investigate model Hamiltonians and even real material systems. In this work, we have outlined a multi-orbital extension of IPT, and benchmarked it extensively against continuous time quantum Monte Carlo results. Our work is the first systematic study of the multi-orbital Hubbard model using MO-IPT as a solver and varying parameters such as filling, Hund’s coupling, and Coulomb repulsion, as well as including crystal field effects and application to real materials. One of the main bottlenecks in methods based on spectral moment expansions is the evaluation of high-order correlation functions. We find that including such correlations that are beyond two-particle type through approximate methods such as CPA or lower order decomposition, can lead to spurious features at the chemical potential. We find the best benchmarks simply by neglecting correlations beyond two-particle. We conjecture that evaluation of the higher-order correlations through exact methods such as ligand field theory might be able to circumvent the issues mentioned abovePhysRevB.77.195124; PhysRevB.57.6884. We are presently implementing such a procedure. This procedure will also enable us to treat the Hund’s coupling term in the rotationally invariant form rather than the simpler and approximate density-density type treated in the present work. We are also planning to extend our method to incorporate the off-diagonal hybridization, which at present is very difficult to handle by numerical exact methods like CTQMC. Apart from the benchmarks for model Hamiltonians in various parameter regimes, we have also carried out a GGA+DMFT(MO-IPT) study of the perovskite SrVO, and compared our predicted photoemission spectrum with experiments and results from other methods. The agreement with experiments was found to be excellent. A full scale implementation of the method outlined here, with detailed instructions for installation and use may be found at http://www.institute.loni.org/lasigma/package/mo-ipt/.

## V Acknowledgments

We thank CSIR and DST (India) for research funding. Our simulations used an open source implementationHafer of the hybridization expansion continuous-time quantum Monte Carlo algorithmComanac and the ALPS Bauer libraries. This work is supported by NSF DMR-1237565 and NSF EPSCoR Cooperative Agreement EPS-1003897. Supercomputer support is provided by the Louisiana Optical Network Initiative (LONI) and HPC@LSU. NSV acknowledges an international travel fellowship award from IUSSTF. Dasari acknowledges the hospitality of the department of Physics & Astronomy and the Center for Computation Technology, at Louisiana State University.

## Vi Appendix

In this appendix, we provide the derivations of the unknown parameters appearing in the MO-IPT ansatz for the self energy (Eq. (9)).

Derivation for A:

The spectral representation of the -orbital Green’s function is given by

(36) |

This can be Taylor expanded to obtain the Green function in terms of spectral moments,

(37) |

where ’s are the spectral moments. We can also represent the Green function in terms of a continued fraction expansion and this is given by

(38) |

By comparing Eq. (38) with Eq. (37), we obtain the continued fraction expansion coefficients in terms of spectral moments. Now we can calculate the spectral moments exactly up to any order by using the following expressionsNolting1:

The relation between the first few spectral moments and the continued fraction expansion coefficients is given by,

(39) |

(40) |

(41) |

(42) |

(43) |

For sufficiently large values of , one can truncate the continued fraction expansion of the Green’s function (Eq. (38)) at the appropriate level and take the limit . Up to the second order moment

(44) |

After substituting the continued fraction expansion coefficients in Eq. (44), we find the self energy contribution to the Green’s function in the high frequency limit as

(45) |

In the high frequency limit the self energy ansatz reduces to the following form:

(46) |

It is easy to show that in the limit of high frequencies, has the following formMartin3,

(47) |

Here n is the Hartree-corrected impurity charge because the propagators used in the second order pair bubble diagram are Hartree-corrected propagators. We obtain the expression for by substituting Eq. (47) in Eq. (46) and comparing with Eq. (45) as

Note that a two-particle correlation function is needed to find .

Derivation for B:

The relation between the impurity Green’s function and the self energy in the atomic limit is,

(48) |

where the self-energy, may be represented as a continued fraction:

(49) |

As a simple case we consider only two poles in the self energy. In principle we can keep all the poles of the self energy but the difficulty is that a pole of order involves the order correlation function. These functions are very hard to calculate without making approximations. With the two pole approximation for the self energy, Eq. (49) reduces to the following formOudovenko

(50) |

(51) |

(52) |

(53) |

In the atomic limit (V 0), the second order pair bubble diagram reduces to the following formKajueter1; Martin1,

(54) |

Here is the pseudo-chemical potential. As mentioned earlier, we find this quantity by satisfying the Luttinger’s theorem or equivalently the Friedel’s sum rule. Now, the self energy ansatz becomes

(55) | |||||

By comparing Eq. (55) with Eq. (49) we find the expression for in terms of spectral moments as,

(56) |

After substituting the spectral moments in Eq. (56) B becomes,