# Multi-orbital simplified parquet equations for strongly correlated electrons

###### Abstract

We extend an approximation earlier developed by us for the single-impurity Anderson model to a full-size impurity solver for models of interacting electrons with multiple orbitals. The approximation is based on parquet equations simplified by separating small and large energy fluctuations justified in the critical region of a pole in the two-particle vertex. We show that an -orbital model with most general interaction is described within this approximation by matrices and is Fermi liquid in the metallic phase. We explicitly calculate properties of a paramagnetic solution of a two-orbital Hubbard model with a Hund exchange and orbital splitting within the dynamical mean-field approximation. We trace the genesis of a metal-insulator transition induced by a crystal field and vanishing of the Kondo quasiparticle peak in strongly correlated orbitally asymmetric systems.

###### pacs:

72.15.Qm, 75.20.Hr## I Introduction

Understanding the effects of electron correlations is important not only for comprehending the behavior of a wide class of materials. It also allows one to get to fundamental principles of quantum theory. Quantum dynamical effects caused by electron correlations are still in the center of interest of condensed-matter physicists. Although a lot of aspects of electron correlations have been successfully explained, there are phenomena the picture of which is incomplete or lacks substantial pieces. The main obstacle in reaching a complete knowledge of the effects of electron correlations is the difficulty of a reliable description of dynamical many-body phenomena.

The major effort in grasping strong electron correlations is focused on single-impurity Anderson model (SIAM).Anderson61 () This model was conceived to provide a simple picture of formation of local magnetic moments. It appeared rather soon that formation of local moments is connected with the Kondo effect.Kondo64 () Since then SIAM has been attracting standing attention of condensed-matter theorists and experimentalists. The first wave of interest in SIAM was crowned by finding an exact ground state of the model with a flat band of conduction electrons for all interaction strengths.Tsvelik83 () The Kondo effect in SIAM was then explained and quantitatively described. The exact expression for the Kondo asymptotics in SIAM has become a hallmark and reference for dynamical effects of strong electron correlations.

Impurity models suit for a description of quantum dots and nano-particlesBorda03 (); Choi05 () but are unable to describe bulk, macroscopic effects of electron correlations. It is, however, generally believed that a local electron interaction in the tight-binding scheme offers enough space for the description of most of the dynamical effects of electron correlations. Success of the dynamical mean-field theory (DMFT) of electron correlations has proved this.Georges96 () Dynamical mean-field approach, that is local dynamical approximations, renewed interest in SIAM, since the basic ingredient of DMFT is an impurity solver. The latter essentially is an impurity model with a self-consistent condition on the local one-electron propagators. The exact thermodynamic solution of SIAM, however, cannot be used, since it does not produce the necessary dynamical renormalization of the energy spectrum. A new wave of attempts to find approximate dynamical solutions of SIAM in the strong-coupling regime arose.

Presently, the most comprehensive quantitative approaches producing dynamical properties of the impurity models are quantum Monte Carlo (QMC)Fye88 () with its continuous-time extension,Rubtsov05 (); Werner06a () exact diagonalization (ED) with a discretized bath,Caffarel94 () and the numerical renormalization group (NRG).Wilson75 (); Krishnamurthy80 () The first two approaches aim at a precise thermodynamics of the impurity problem, while the last one best reproduces low-temperature and low-energy scales of impurity models. The first two approaches work naturally in the Matsubara formalism and need a numerical analytic continuation to determine spectral properties, while NRG can be straightforwardly extended to dynamical and spectral quantities.Costi94 (); Bulla98 () The numerical renormalization group represents now-a-days the most accurate quantitative low-temperature and low-energy dynamical solution of impurity problems.Hewson93 ()

Apart from pure numerical methods a series of analytically controlled approximations have been proposed. An approximation based on low-frequency renormalizations of Fermi-liquid parameters is able to reproduce the Kondo scale in SIAM and fits for heavy Fermion systems.Hewson93a () It is not, however, suitable for the description of critical instabilities of the Fermi liquid state. Strong-coupling expansions based on the infinite-interaction model Keiter70 (); Pruschke89 () reproduce as well the Kondo scale but fail to reproduce the Fermi liquid regime in the weak-coupling limit.

Standard truncated weak-coupling expansions fail to reproduce correctly the strong-coupling regime. An improvement can be reached if the expansion is applied on vertex functions. Expansions for two-particle vertices are mostly based on multiple scatterings of electrons on electrons (T-matrix, TMA) or electrons on holes (random phase, RPA). The non-renormalized multiple electron-hole scatterings lead to an unphysical pole that must be removed by a self-consistent treatment. Standardly introduced one-electron self-consistency (fluctuation-exchange approximation) removes the pole but fails to reproduce the Kondo scale.Bickers89 () A static one-electron self-consistency of the Hartree type was introduced in an effective approximate scheme, called local-moment approach.Logan98 () The expansion for the polarization operator (RPA vertex) uses one-electron propagators with a broken spin symmetry in this approximation. The spin symmetry of the equilibrium thermodynamic state is then recovered via a construction of the self-energy calculated from the polarization operator. In this way the RPA pole is not crossed and the strong-coupling limit keeps the polarization operator in the critical region of the RPA pole. The spectral function displays the correct Kondo asymptotics of the quasiparticle peak. This theory, however, contains an artificial symmetry breaking and an auxiliary fitting parameter, the length of the local magnetic moment, in the strong-coupling regime.

An alternative way of introducing a self-consistency at the two-particle level was proposed by us.Janis07a (); Janis08 () We suggested to use parquet equations constructed from multiple electron-electron and electron-hole scatterings and use the criticality of the RPA pole for a simplification leading to a manageable theory. This approximation allows us to control analytically the RPA pole and the genesis of the Kondo asymptotics. The derived Kondo behavior reproduces universal features of the exact Bethe-ansatz solution.

Real materials are more complex and cannot be described by elementary impurity models. One has to extend the methods to multi-orbital Hamiltonians. Not all the methods used as impurity solvers for single-orbital models work also in the multi-orbital case. The exact solution can be extended to models only in the limit .Ogievetski83 () There are natural extensions of QMCBonca93 (); Kotliar06 (); Werner06b (), ED,Liebsch05 (); Inaba05 () and NRGGalpin05 (); Bulla08 () methods to multiple-orbitals. All the numerically exact approaches are, however, computationally rather expensive for realistic atomic structures. There is hence still a need for simpler, computationally less expensive semi-analytic approaches that could qualitatively correctly describe the transition from the weak to the strong coupling regimes in multi-orbital impurity and lattice models.

Some of the existing analytically controlled approaches to strong electron correlations have recently been extended to multi-orbital models, mostly Anderson models. Various spectral and thermodynamic quantities have recently been calculated within the renormalized Fermi-liquid expansion,Nishikawa10 () fluctuation exchange,Drchal05 () or local-moment approach.Galpin09 (); Kauch09 () The aim of this paper is to present a generalization of a two-particle approach based on simplified parquet equations to multi-orbital models and create a full-scale impurity solver for realistic models of correlated electrons. Unlike the most of the other multi-orbital approaches, we formulate the method so that it can directly be applied not only to impurity models but also to lattice models within the dynamical mean-field theory. The approximation is formulated in real frequencies and hence the spectral properties are directly available. The low and high frequency behavior of the one-electron spectral function, its three-peak structure, is qualitatively well reproduced and the method allows one to control analytically the critical behavior of Bethe-Salpeter equations for the two-particle vertex.

The paper is organized as follows. We introduce a general multi-orbital perturbation expansion for two-particle vertices in Sec. II. In Sec. III we derive simplified parquet equations with electron-hole and electron-electron multiple scatterings and show how the self-energy is calculated from the full two-particle vertex. We present explicit equations to be solved in a two-orbital case in Sec. IV. The derived approximation is used to obtain a numerical solution of a two-orbital model with Hund’s coupling and orbital splitting in Sec. V. We also compare our results with those from other approaches there. Reliability and domains of applicability of the derived approximation are discussed in the final Section VI.

## Ii Perturbation expansion for multi-orbital models of correlated electrons

We start with a model Hamiltonian with multiple orbitals, where the band structure is resolved and kinetic energy is diagonal in orbital indices. The unperturbed Hamiltonian then reads

(1) |

where , , and are momentum, spin and orbital index of valence electrons. The interaction multi-orbital Hamiltonian can generally be represented in the direct space as

(2) |

where all possible inter-orbital transitions are included. Orbital-momentum conservation poses a restriction on the selection of orbital indices of the interaction matrix being . We use a tight-binding approximation taking the coordinate from the lattice space. We will assume the interaction to be local in the lattice space, that is, Hubbard-like. A graphical representation of the interaction matrix is plotted in Fig. 1, where we used latin indices for all the dynamical degrees of freedom characterizing the electron states. They are generally four-momenta or Matsubara frequencies in local approximations (dynamical mean-field theory) together with orbital indices. We use this abbreviated notation to simplify summations over intermediate states in the multi-orbital perturbation expansion.

It is necessary to classify all possible interconnections of interaction vertices via fermion lines for building up the complete multi-orbital perturbation expansion. It is not an easy task and we will sum only selected classes of two-particle diagrams. We use the graphical representation of the interaction vertex from Fig. 1 and choose the upper part of it so that spin (not necessarily orbital index) is conserved along the horizontal fermion lines. We choose electron propagation from left to right and hole propagation from right to left. We have three possibilities how to connect two interaction vertices by two fermionic lines. They correspond to three two-particle scattering channels. Janis99b ()

The first possibility to connect interaction vertices is via propagation of two electrons, that is, by parallel fermionic lines, cf. Fig. 2. Mathematically this process can be represented via a matrix multiplication

(3) |

where the latin indices stand for all dynamical variables of the electron (momentum, frequency and orbital index). They obey a conservation law .

Second connection of interaction vertices is via propagation of an electron-hole pair, that is, by antiparallel fermionic lines. Here we must distinguish two options: electron-hole scatterings and polarization bubbles. The former process is plotted in Fig. 3 and is represented via another matrix multiplication

(4) |

where again the orbital-monentum conservation law holds.

A polarization bubble, shielding the interaction, is represented via a matrix multiplication

(5) |

and is plotted in Fig. 4. Notice that spin in this last process is not conserved, but orbital-momentum conservation restricts the intermediate indices of the dynamical variables to .

We use these elementary scattering processes to construct the full two-particle vertex from a two-particle perturbation expansion. Once we have this vertex we construct the one-electron self-energy as the fundamental quantity from which we derive all thermodynamic properties. The self-energy is calculated from the two-particle vertex via the Schwinger-Dyson equation, diagrammatically represented in Fig. 5. Mathematically it reads

(6) |

Again the conservation law for intermediate variables holds. The one-electron propagators in the perturbation expansion can either be chosen as fully dynamically renormalized as suggested by Baym and Kadanoff,Baym62 () or only statically renormalized via the Hartree approximation securing that the particle densities do not change during the summation of diagrams. We discussed earlier that the latter option is more appropriate in the strong-coupling regime, since we can better control the two-particle singularity in the electron-hole correlation function and consequently both low and high energy scales are qualitatively better reproduced than in approximate schemes with the fully renormalized propagators.Janis07a ()

## Iii Simplified parquet equations: matrix formulation

The generic function for our perturbation expansion is the two-particle vertex . Its approximate expression together with the Schwinger-Dyson equation (6) define our approximate scheme. We hence start with simple approximations on vertex . We resort our analysis only to local approximations, that is, we stay within the dynamical mean-field theory with only local fluctuations and Matsubara frequencies, orbital index and spin as dynamical variables. We will treat these degrees of freedom separately.

The simplest approximations for vertex , beyond the bare interaction, are multiple scatterings of the same form as defined by elementary processes from the preceding section. Sums of all multiple scatterings of the same type are solutions to Bethe-Salpeter equations. They can be represented with the appropriate matrix multiplications. For electron-hole and electron-electron multiple scatterings we obtain

(7) | |||

(8) |

respectively. The kernels of these equations are composed from the bare interaction and either the electron-hole bubble or the two-electron propagator. The kernel in the electron-hole channel explicitly reads

(9) |

with the electron-hole bubble

(10) |

Analogously we construct the kernel of the Bethe-Salpeter equation (8) with multiple electron-electron scatterings. The two-electron propagator then is

(11) |

We denote fermionic and bosonic Matsubara frequencies and , respectively (Boltzmann constant ). Notice that the transfer frequency , conserved in electron-hole scatterings, is different from that conserved in the electron-electron scattering channel, . If the incoming and outgoing frequencies of one electron are and , respectively, we have .

We do not take into account the third channel with triplet polarization bubbles. There are two reasons why we can afford to do this. First, the effects of the scatterings in the vertical channel are similar to those from the electron-hole channel, at least qualitatively. We hence do not lose track of any relevant effect in the strong-coupling regime. Second, taking into account only two scattering channels simplifies significantly approximate schemes for the two-particle vertex.Janis09 ()

From the solutions of the Bethe-Salpeter equations for and we obtain the full two-particle vertex

(12) |

Such a vertex defines the self-energy of the so-called FLEX approximation if the one-electron propagators in the Bethe-Salpeter equations (7) and (8) are fully renormalized, that is, they contain the self-energy determined by the approximate vertex from Eq. (12). Such an approximation is known to be unreliable in the strong-coupling regime.Janis99b ()

Intermediate electron coupling is marked by proximity of a pole in the electron-hole vertex . If the transfer frequency then we can define a small dimensionless scale

(13) |

approaching zero in the random-phase approximation with , a solution of Eq. (7), when . We used to denote the set of eigenvalues of matrix . This pole is unphysical and indicates that the perturbation expansion for must be renormalized. The critical interaction from RPA defines an upper limit on the weak-coupling regime where perturbation theory with selected classes of generic diagrams can safely be applied. We speak about strong coupling if , where only self-consistent approximations may be reliable.

One, and most often used option how to suppress the unphysical pole, or better to shift it to infinite interaction strength is to employ one-particle dynamical self-consistency. We use the self-energy calculated from the Schwinger-Dyson equation (6) with vertex from Eq. (12) to renormalize the one-electron propagators in Bethe-Salpeter equations (7) and (8). This one-particle self-consistency succeeds in shifting the unphysical pole to infinity, but the Kondo regime defined as is not properly reproduced. Moreover, FLEX-type approximations completely smear out the large-frequency structures. There are no satellite Hubbard bands.Janis99b ()

A better systematic alternative how to move the RPA critical interaction to infinity is to introduce a two-particle self-consistency. The most straightforward way to do so is to use the parquet approach. Its basic idea is to replace the bare interaction in Bethe-Salpeter equations (7) and (8) by irreducible vertex functions and , respectively. In the exact theory then and equation (12) transforms to . When we combine the two Bethe-Salpeter equations with the parquet equation (12) modified as above we obtain a set of non-linear integral equations for the irreducible vertices and finally also for the full vertex . It is, however, impossible to solve the parquet equations in the strong-coupling regime with a vanishing Kondo scale . The problem lies in that all the vertex functions, solutions of the full parquet equations, generally depend on three frequencies and there is presently no method, either analytic or numerical, that would predict the full analytic structure of the two-particle vertices. We hence have to approximate the full parquet equations.

We proposed a simplification of the parquet equations in the Kondo regime of the single-impurity Anderson model.Janis07a (); Janis08 () This simplification utilizes a partial separation of large and small frequency fluctuations in the Kondo regime when the solution of the Bethe-Salpeter equation in the electron-hole channel is almost singular. In the non-renormalized perturbation expansion it is the RPA pole. The principal idea of our simplification is that we separate large fluctuations diverging at the critical point from those remaining finite. In this way we describe correctly universal quantities connected with the critical point in the electron-hole Bethe-Salpeter equation. We neglect all finite fluctuations and keep only the relevant ones diverging at the critical point. We then replace all finite frequency-dependent functions by constants and keep only the frequency dependence in the variables controlling large fluctuations. This simplification is kind of a mean-field (simplest) approximation for the dynamics in the strong-coupling regime.

We know from the analysis of multiple electron-electron and electron-hole scatterings that only the latter can cause a singularity in the vertex function. It means that the irreducible vertex in the electron-hole channel remains finite. We hence can replace it by an effective interaction that we denote . The renormalized Bethe-Salpeter equation in the electron-hole channel is RPA-like with the following replacement

(14) |

The effective interaction will be determined from the Bethe-Salpeter equation in the electron-electron channel by using the parquet self-consistency. Within the scope of our approximation the Bethe-Salpeter equation in the electron-electron scattering channel holds only in an averaged form smearing out finite frequency differences. We proposed in Ref. Janis08, a decoupling of the Bethe-Salpeter equation in the electron-electron scattering channel so that an integral equation is converted to an algebraic, solvable one. In this approximate simplification an averaged integral kernel is assumed to diagonalize the integral equation. Such a decoupling is not uniquely defined and should be chosen so that the effective interaction is a real self-adjoint (positive definite) matrix. An equation for the effective interaction derived from a decoupling with the right-averaged integral kernel used in the single-orbital case reads

(15) |

where we denoted

(16) |

This one-sided averaging in the multi-orbital case can lead to a matrix of the effective interaction that is not self-adjoint. To avoid this spurious behavior we introduce also a left-averaged integral kernel in the electron-electron Bethe-Salpeter equation. We then define another equation for the effective interaction

(17) |

with a left-averaged integral kernel

(18) |

We demand that the matrix of the effective interaction be self-adjoint. It is always the case if we use a sum of equations (15) and (17) to determine the effective interaction. That is, we choose an equation symmetric with respect to right and left averaging of the integral kernel of the Bethe-Salpeter equation in the nonsingular electron-electron channel.

The dominant contribution to the integral kernel in the electron-electron scattering channel comes from the singular part of the vertex from the electron-hole scattering channel, being . Since the averaged integral kernels and depend on a fermionic Matsubara frequency, we assume that Eqs. (15), (17) are valid only for the lowest frequency, being zero at zero temperature. At finite temperatures we use the lowest or a few low frequencies, symmetric around zero, and add all equations for these frequencies. The effective interaction is then determined form an average over Eq. (15) and Eq. (17) for a few lowest-lying fermionic frequencies. The low-temperature asymptotics is in this way well reproduced. The higher the temperature the more Matsubara frequencies we must take into account to maintain accuracy of the zero-temperature solution. We know, however, that thermal fluctuation suppress the Kondo resonance and smear out the RPA singularity in the electron-hole correlation function.Janis08 ()

Equations (7), (15)-(18) are our simplified parquet equations determining the two-particle vertices and . The full two-particle vertex is determined from equation (12). Within our approximation it reduces to the renormalized RPA vertex . We use it to determine the self-energy. The respective Schwinger-Dyson equation reads

(19) |

where conservation of the orbital index holds.

The approximation is almost complete. The last thing we have to do is to specify the one-electron propagators in the parquet equations for the two-particle vertices. We know that the full one-particle self-consistency where the one-electron propagators are renormalized by the self-energy from Eq. (19) does not lead to the Kondo scale as and also smears out the satellite Hubbard bands.Janis08 () We hence use the Hartree one-electron propagators that in the dynamical mean-field theory is

(20) |

and for the impurity Anderson model

(21) |

We denoted , where is the center of the ’s orbital and is the chemical potential fixing the total charge density. The Hartree propagators are statically self-consistent, that is, the particle densities are determined from the fully renormalized propagator

(22) |

We denoted the dynamical part of the self-energy, a correction to the Hartree term. The chemical potential is then adjusted so that the total particle density is fixed. The approximate scheme based on the simplified parquet equations is now complete. We demonstrate explicitly applicability of this approximation on a two-orbital model in the next sections.

## Iv Two-orbital model

We use a general two-orbital interacting Hamiltonian comprising also Hund’s exchange couplingImada98 ()

(23) |

where . A graphical representation of single terms contributing to this two-orbital Hamiltonian is presented in Fig. 6. Each input parameter of this Hamiltonian will be renormalized independently.

Only one of the bare interacting terms denoted is spin triplet. This term does not mix with the other scattering processes in the Hamiltonian from Eq. (23). The simplified parquet equations for this term are identical with those from the single-orbital case. The effective interaction is

(24) |

with the averaged integral kernel from the Bethe-Salpeter equation in the electron-electron channel

(25) |

The dynamical part of the (reducible) vertex from the electron-hole scattering channel is

(26) |

The other terms in the interaction Hamiltonian from Eq. (23) mix in the parquet equations and must be treated simultaneously. Since only specific inter-orbital processes are allowed we can simplify the representation of the interaction term from the preceding section and reduce super-matrices (tensors) with four indices to regular matrices with two indices. A relation between matrix indices and super-indices can be chosen as follows , , , . If we go over from super-matrices to matrices we have to distinguish two representations corresponding to two matrix multiplications used in constructing the Bethe-Salpeter equations. The interaction matrix suitable for summation of multiple electron-electron scatterings is

(27) |

The interaction matrix for the electron-hole multiple scatterings interchanges the two Hund couplings and . Hence we have

(28) |

Each matrix representation is chosen so that the scattering processes in the respective channel are simple matrix multiplications.

The matrix integral kernel for the Bethe-Salpeter equation in the electron-hole channel is

(29) |

while the right-averaged kernel in the Bethe-Salpeter equation in the electron-electron channel reads

(30) |

with

(31) |

and . Analogously we construct the left-averaged kernel , if necessary. The matrix form of the simplified parquet equations then is in the electron-hole channel

(32) |

and in the electron-electron channel with right and left-averaged integral kernels

(33) |

All matrices appearing in the parquet equations (32) and (33) are block diagonal. We denote

(34) |

and solve the parquet equations for each block separately. The solution in the electron-hole scattering channel has the following representation

(35) |

and

(36) |

Notice that the Hubbard interaction splits into two renormalized values and if the bubbles in different orbitals are different, . Note that we suppressed the frequency index at the electron-hole bubble leading to a frequency dependence of vertex .

We now choose the singular part of the renormalized electron-hole vertex and use it to obtain a solution of the Bethe-Salpeter equation in the electron-electron channel with an averaged (right, left) kernel that determines the matrix of the effective (renormalized) interaction. Its upper matrix block reads

(37) |

The lower block analogously is

(38) |

The electron-electron bubble and vertices are taken here (zero temperature) only at zero frequency. If we must symmetrize Eqs. (37) and (38) by adding the Bethe-Salpeter equation with the left averaged kernel so that to keep the matrix of the effective interaction self-adjoint.

We use the electron-hole dynamical vertex resulting from the simplified parquet equations and calculate the self-energy correction to the Hartree term from Schwinger-Dyson equation (19). It is in the super-matrix representation

(39) |

The dynamical part of the self-energy is then used to determine the particle densities from Eq. (22) so that the static one-electron self-consistency is fulfilled.

## V Results

We analyze a general two-band Hubbard model within the dynamical mean-field theory. We set the model parameters , and , which corresponds to a spherically symmetric situation.Imada98 () We choose the Hund exchange so that all three electron interactions remain positive. We use the one-electron propagator from Eq. (20) with the semi-elliptic density of states, the bandwidth of which is . We investigate only the spin symmetric case where the one-electron propagators are spin-independent. The critical interaction in the random phase approximation, separating the weak and strong coupling regimes, then is . It does not depend on the number of orbitals if the Hund coupling . We perform all calculations at zero temperature.

We first investigate the simplest situation of the charge-symmetric model at half filling, that is , with no orbital splitting and with vanishing Hund coupling, . This situation mimics the single-orbital model. The emergence of the Kondo scale is plotted in Fig. 7. The quasiparticle peak starts to develop from the bare density of states for interactions around the critical one from RPA, . The satellite Hubbard bands are well formed in the strong-coupling regime and lie slightly beyond the atomic values . The Kondo scale remains exponentially small but non-zero even for large interaction strengths. The Kondo asymptotics for the semi-elliptic density of states is , where is the density of states at the Fermi energy. This result is in discrepancy with general confidence that the Hubbard model at half-filling undergoes in the paramagnetic phase a Mott-Hubbard metal insulator transition at a finite interaction strength. The existence of a metal-insulator transition in DMFT is concluded from extrapolations of low-temperature Monte-CarloGeorges96 () or NRGBulla01 (); Bulla08 () data. Although there is no doubt about the existence of an insulating solution for very large interaction strengths matching the atomic limit, numerical methods cannot presently exclude the existence of a metallic solution with a very narrow (below numerical resolution) quasiparticle peak. The existence of a sharp metal-insulator transition has not yet been rigorously proved and there is at present neither a consistent solution at the critical point of a transition nor there is an analytic formula for the vanishing width of the quasiparticle peak. If there is a metal-insulator transition, the metallic side of the transition cannot be a Fermi liquid. The low-frequency asymptotics of the electron-hole correlation function in Fermi liquid is due to sum rules of order . Such a behavior is then incompatible with Bethe-Salpeter equations that have to be fulfilled self-consistently in the electron-hole and electron-electron scattering channels by an exact solution. Only integrable singularities can emerge in an exact solution and Fermi liquid hence does not allow for a divergence in the electron-hole correlation function expected at the Mott-Hubbard metal-insulator transition. In this respect the simplified parquet approximation is consistent with the exact two-particle behavior.

Nonexistence of a sharp metal-insulator transition in the simplified parquet equations is not caused by approximations made in this approach but rather it is a consequence of a Fermi-liquid solution fulfilling self-consistently Bethe-Salpeter equations simultaneously in the electron-hole and electron-electron scattering channels. There is an insulating paramagnetic state in parquet-based theories for (critical interaction of the RPA pole) being a superposition of Hartree polarized solutions, but it coexists with the metallic one in the strong-coupling regime and has a higher total energy. A metal-insulator transition in the Fermi liquid can either be connected with a breaking of a global symmetry and/or can be noncritical as we discuss later for orbital splitting in a crystal field.

We compared the spectral function calculated from the simplified parquet equations with other diagrammatic approximations based on summations of Feynman diagrams for multiple two-particle scatterings. We chose a value of the interaction strength so that to be in the strong-coupling regime where renormalizations of the perturbation expansion are necessary. We used the spectral function for the self-consistent T-matrix, iterated perturbation theory in second order (IPT) and for T-matrix with IPT (topological) self-consistency, that is, the local renormalized propagator of the Hubbard model is defined as . The results of the comparison are plotted in Fig. 8. We can see that except for the fully self-consistent T-matrix, the chosen three approximations display the expected three-peak structure. They differ, however, in the way the satellite bands are attached to the central quasiparticle peak. The simplified parquet approximation keeps a small finite density of states between the central and satellite bands. It means that this approximation allows for an energy transfer between the central and the satellite bands. Approximations based on the IPT self-consistency develop a quasi-gap between the low and high energy states and no energy exchange occurs. The satellite bands are less spread in the parquet approximation than in the IPT-based ones. The detailed shape of the central quasi-particle peak in these approximations is magnified in Fig. 9. The most narrow is the central peak for IPT suggesting vicinity of a metal-insulator transition with zero width of the quasiparticle peak. The parquet-based approximations, on the other hand, slow down the tendency toward a metal-insulator transition. They actually screen it, since there may be no critical metal-insulator transition in these theories as discussed above.

The simplified parquet equations show best agreement with the exact solution of the single-impurity Anderson model at half filling, where the density of states at the Fermi surface is fixed. The same holds also for other fillings at intermediate and strong coupling calculated within the dynamical mean-field approximation. Apart from the charge-symmetric situation the simplified parquet equations reproduce qualitatively well the low and high frequency features of the spectral function, although the width of the quasiparticle peak and the Kondo temperature derived from it do not follow as closely the exact solution for the flat-band model as for the half-filled case.Janis08 () The spectral function of a strong-coupling solution changes when receding from half filling so that the central quasiparticle peak broadens and slightly moves below the Fermi surface. The lower satellite band moves toward the Fermi energy. The lower band and the central peak eventually merge, the quasiparticle peak is absorbed by the satellite one, when the electron-hole asymmetry becomes prominent. This feature seems to be universal for all theories and is demonstrated in our approach in Fig. 10.

We hitherto have shown the results for the two-orbital model without Hund exchange and orbital splitting, that is for and . Such a situation is close to the single-orbital model as we demonstrated. The most important change in the model with nonzero Hund coupling is a splitting of the interaction strength, that is, and . It means that the interaction strength controlling the approach to the critical point in the electron-hole channel is the Hubbard between different spins on the same orbital. The Hund exchange decreases the critical interaction at which the electron-hole vertex diverges in the bare perturbation theory. It would then narrow the width of the Kondo resonant peak at a fixed renormalized Coulomb repulsion . In the real situation we, however, do not fix the renormalized but rather the bare Coulomb repulsion. We find that though larger Hund exchanges lead to higher values of the averaged kernels , as we can see from Eqs. (37) and (38), they simultaneously cause decrease in the renormalized repulsions at a fixed bare Coulomb interaction. Consequently, the weight of the central quasiparticle peak in the self-energy decreases. We plotted the impact of the Hund coupling on the spectral function in Fig. 11. We observe that the width of the Kondo peak shrinks but only on a narrower frequency interval near the Fermi energy. The Hund coupling at a fixed bare Coulomb repulsion then leads to filling of the states between the central peak and the satellite bands. Further, the satellite bands are positioned closer to the central one and become less pronounced with the increasing Hund exchange.

The most interesting on the two-orbital model are changes in the spectral function due to a broken orbital symmetry induced by a crystal field. That is, when the one-electron propagators are different on different orbitals. We achieve this here by splitting the centers of the two bands. We set , but keep the same band structure, . The effect of the orbital splitting depends on whether we have a finite or infinite bandwidth of macroscopically occupied energies. Orbital splitting leads in the former case to a metal-insulator transition unlike the latter one. The genesis of a metal-insulator transition in the model with the semi-elliptic density of states is plotted in Fig. 12. We chose a moderate interaction strength, , so that to be able to follow the process more closely. When we increase the orbital splitting parameter there is almost no observable change at the Fermi energy with a pronounced quasiparticle peak. It slightly splits into two but stays at the Fermi energy. More apparent is the change caused by orbital splitting in the bulk of the energy band where the states are rearranged so that the orbital polarization is enhanced. Such a response on the band splitting effectively pushes the system out from the critical region of the RPA pole. With the increasing band splitting the dynamical renormalization of the band structure, apart from the vicinity of the Fermi energy, becomes less important and the solution approaches the Hartree one. The process continues up to the point when one of the orbital saturates and the other is emptied. The insulating solution is then that of the Hartree approximation with no dynamical correlations. Changes in the quasiparticle peak, split and a broadening, are observed only close to the transition. They vanish practically at the transition point. There are, however, no singular two-particle functions at the transition point, the metal-insulator transition is hence noncritical. It is worth noting that the Kondo quasiparticle peaks exist in a pronounced way only in the strong-coupling regime, cf. Fig. 13.

We further studied the effect of the Hund exchange on the orbital metal-insulator transition. We set the parameters of the model so that we could compare our results with the existing Monte-Carlo simulations of Ref. Werner07, . We chose three values of the bare orbital splitting . For each value of the band splitting we fixed four ratios of the Hund coupling and plotted the filling (per spin) of the upper orbital as a function of the bare interaction at half filling in Fig. 14. We observe that small Hund couplings do not prevent the system from falling into the orbitally polarized state as demonstrated in Fig. 12. For large the occupation of the upper orbital decreases practically monotonically with the increasing interaction. For smaller splittings the filling increases and may get to the Kondo regime at intermediate interaction strengths where both orbitals are almost equally occupied (). Unlike the Monte-Carlo data, we do not reach, however, the Mott insulating phase from the Fermi-liquid solution. If the system does not stay in the Kondo regime for ever. Strong electron repulsion increases orbital polarization due to the Hartree shift and the system eventually goes over into a polarized insulator . For ratios the Hartree term starts to act against the crystal field and diminishes orbital polarization. The filling then asymptotically approaches half filling. Apart from nonexistence of the Mott transition these effects of the Hund coupling on the occupation of orbitals in a crystal field are described by the simplified parquet equations in a good agreement with the Monte-Carlo simulations. It is worth mentioning that the Kondo resonance for survives to a certain value of the interaction strength beyond which it is rather abruptly suppressed, similarly to the scenario from Fig. 12. Vanishing of the Kondo resonance is followed by a sharp saturation of orbital polarization. The numerical solution then becomes unstable, since the metallic phase in the extreme Kondo regime and the orbitally polarized insulating solution coexist with almost the same energies.

The scenario of orbital splitting changes if the bare energy bands are infinite. When an orbital symmetry breaking is switched on, the system with infinite band-widths behaves for small splittings analogously to the model with finite energy bands. The changes take place in the bulk, away from the Fermi energy so that the orbital splitting is enhanced, cf. Fig. 15. Similarly there are almost no observable changes on the quasiparticle peak. The system is slowly pushed away from the critical region of the RPA transition by increasing the orbital splitting and approaches orbital saturation. But due to the infinite band width, there is no transition to the Hartree insulator and the density of states at the Fermi energy remains nonzero for all interaction strengths. The Kondo peaks do not survive to infinite interaction. A sudden macroscopic cleft breaks in between the central quasiparticle peaks of different orbitals at a finite interaction strength. The split Kondo peaks are then soon absorbed by the bulk energy bands. The solution approaches a weak-coupling (Hartree-Fock) state with almost separated energy bands centered around energies , where is the orbital moment.

## Vi Discussion and conclusions

Dynamical effects of strong electron correlations in impurity models manifest themselves in a three-peak structure of the spectral function with a narrow central quasiparticle peak near the Fermi energy and two satellite bands centered around eigenvalues of the local part of the interacting Hamiltonian. There is no exact solution proving such a picture, but a number of approximate numerical and analytic-numerical solutions confirm this behavior of the spectral function. There are a number of such approximate schemes for the simplest single-orbital impurity and DMFT models. Extensions of approximations showing the expected three-peak spectral function in strong coupling to multi-orbital models face a considerable increase of complexity in the strong-coupling regime. Monte-Carlo simulations as well as the numerical renormalization group applied on multi-orbital models need to resort to additional approximations to reach sufficient accuracy of dynamical and spectral properties.

Most of the existing approximations on dynamical properties of strongly-correlated electron systems make directly accessible only one-electron properties. Only few of them, even in single-orbital models, are able to approach two-particle quantities and in particular, to control two-particle divergencies that develop in Bethe-Salpeter equations. We extended in this paper an approximation based on parquet equations controlling the RPA singularity in the electron-hole scattering channel to a general multi-orbital model of correlated electrons. Approximations used to simplify complexity of the full parquet approach resulting from a separation of large and small dynamical fluctuations in the critical region of the RPA pole allow us to reproduce the expected three peak structure of the spectral function of SIAM in the Kondo regime. An extension of this approximation from single to multi-orbital models proved to be manageable with only a moderate increase in numerical requirements. An -orbital model is in the most general case described by matrices and the numerical expenses increase polynomially with the number of orbitals. In real situations symmetries and conservation laws further reduce complexity of the simplified parquet equations.

The impurity solver for the dynamical mean-field theory we built up within the parquet approximation suits best for describing the dynamical and spectral properties of the metallic (Fermi liquid) phase from weak up to moderately strong electron correlations. No self-consistency is needed in the weak-coupling regime below the critical interaction strength at which the RPA vertex diverges. The parquet self-consistency is indispensable in intermediate () and strong () coupling. The numerical procedure for the solution of the simplified parquet equations gets unstable for large interaction strengths. We were able to reach a stable metallic solution up to . Unlike other approaches such NRG, QMC or IPT this approximation does not allow for a critical metal-insulator transition without a symmetry breaking in Fermi liquid. This is a general feature of all impurity solutions where Bethe-Salpeter equations in the electron-hole and electron-electron channels should be obeyed simultaneously in a self-consistent manner. Parquet equations are hence unable to describe a Mott insulator without a symmetry breaking. The insulating solution obtained within the local parquet approach is a superposition of fluctuations-free Hartree saturated solutions and the transition to such a state is noncritical.

The simplified parquet equations presented in this paper are not quantitatively as accurate as numerically exact solutions, but they offer a numerically manageable scheme with an analytic control of two-particle singularities in Bethe-Salpeter equations. The self-energy and the spectral function calculated within this approximation obey Fermi-liquid rules in the metallic phase and produce the typical three-peak structure in the density of states in the strong-coupling regime. The width of the central quasiparticle peak is exponentially small for large interaction strengths. This approximation may not reproduce the exact asymptotics of the Kondo scale and temperature in the single-impurity model as well as e. g. the local moment approach, but unlike the latter it does not show any unphysical behavior and does not use auxiliary fitting parameters. The approximation based on the simplified parquet equations is formulated in such a way that it can directly be applied not only to impurity models but to lattice models within the dynamical mean-field approach as an impurity solver.

The impurity solver based on the simplified parquet equations is a semi-analytic theory for local spectral properties formulated in real frequencies with a full control of analytic properties of the approximation. From this reason the approximation best suits to the calculation of spectral properties and the dynamical self-energy directly accessible from a perturbation theory. Its unique feature is that it offers analytic control of two-particle singularities in the Bethe-Salpeter equations. The simplified parquet approximation is well suited to problems involving two-particle dynamical, and in particular critical, functions. The approximation is less suited to global thermodynamic quantities being averages of dynamical functions over frequencies. An extension to finite temperatures does not represent a problem and the effective interactions can be calculated also at nonzero temperatures as discussed in this paper. The quantitative reliability of the simplified parquet equations decreases with the reduced influence of the critical behavior of the vertex function in the vicinity of the RPA pole (Kondo regime). Only in the critical region of a pole in the two-particle vertex we can rely on a separation of large and small dynamical fluctuations, that is, a decoupling of singular and non-singular functions.

To our knowledge, the simplified parquet equations is the only scheme that uses a two-particle self-consistency to control the Kondo asymptotics and criticality of the RPA pole. There are attempts to use parquet equations for renormalizations of the perturbation expansion, but they are based on the full one- and two-particle self-consistency.Bickers04 (); Yang09 () It is, however, important to realize that the full dynamical one-particle self-consistency impedes control of two-particle divergencies and the Kondo scale in the spectral function is no longer correctly reproduced. Moreover, the full dynamical one-electron self-consistency ultimately leads to smearing of the satellite bands.Janis07a () It is hence of importance that the parquet equations are solved with the Hartree propagators so that to reproduce the three-peak spectral function with the correctly scaled central quasiparticle peak.

The extension of the simplified parquet equations was used in this paper as an impurity solver for the two-band Hubbard model with Hund coupling and orbital splitting. We showed emergence of the Kondo quasiparticle peak with the increasing interaction strength at half filling. The quasiparticle peak survives only close to the electron-hole symmetry where multiple electron-hole scatterings are dominant. When deviations from this symmetry are significant the quasiparticle peak is absorbed by one of the satellite bands. We also demonstrated how a metal-insulator transition induced by a crystal field takes place in a model with finite energy bands. The parquet equations do not allow for a critical transition in the local Fermi liquid. We showed that the orbitally induced metal-insulator transition is noncritical, but with a Kondo behavior almost up to the transition point, if the interaction is sufficiently strong.

Last but not least, the approximate impurity solver we built up on parquet equations uses the formalism of Green functions and can conveniently be combined with ab-initio computational schemes to assess correlation effects in materials with a nontrivial atomic structure. It should be suitable in particular for metallic systems where we expect Kondo-like behavior or strong valence fluctuations.

###### Acknowledgements.

We thank V. Drchal and A. B. Shick for valuable discussions and helpful comments. This research was carried out within project AV0Z10100520 of the Academy of Sciences of the Czech Republic and supported in part by Grant No. 202/07/J047 of the Czech Science Foundation.## References

- (1) P. W. Anderson, Phys. Rev. 124, 41 (196).
- (2) J. Kondo, Prog. Theor. Phys. 32, 37 (1964).
- (3) A. M. Tsvelik and P. B. Wiegmann, Adv. Phys. 32, 453 (1983).
- (4) L. Borda, G. Zarand, W. Hofstetter, B. I. Halperin, and J. von Delft, Phys. Rev. Lett. 90, 026602 (2003).
- (5) M. S. Choi, R. López, and R. Aguado, Phys. Rev. Lett. 95, 067204 (2005).
- (6) A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- (7) R. M. Fye and J. E. Hirsch, Phys. Rev. B 38, 433 (1988).
- (8) A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Phys. Rev. B 72, 035122 /2005).
- (9) P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
- (10) M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (1994).
- (11) K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
- (12) H. R. Krishna-murthy, J. W. Wilkins, and K. G. Wilson, Phys. Rev. B 21, 1003, 1044 (1980).
- (13) T. A. Costi, A. C. Hewson, and V. Zlatić, J. Phys.: Condens. Matter 6, 2519 (1994).
- (14) R. Bulla, A. C. Hewson, and T. Pruschke, J. Phys.: Condens. Matter 10, 8365 (1998).
- (15) A. C. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, Cambridge 1993).
- (16) A. C. Hewson, Phys. Rev. Lett. 70, 4007 (1993).
- (17) H. Keiter and J. C. Kimball, Phys. Rev. Lett. 25, 672 (1970).
- (18) T. Pruschke and N. Grewe, Z. Physik B 74, 439 (1989).
- (19) N. E. Bickers and D. J. Scalapino, Ann. Phys. (N. Y.) 193, 206 (19989).
- (20) D. E. Logan, M. P. Eastwood, and M. A. Tusch, J. Phys.: Condens. Matter 10, 2673 (1998).
- (21) V. Janiš and P. Augustinský, Phys. Rev. B 75, 165108 (2007).
- (22) V. Janiš and P. Augustinský, Phys. Rev. B 77, 085106 (2007).
- (23) E. Ogievetski, A,. M. Tsvelik, and P. B. Wiegmann, J. Phys. C: Solid State Phys. 16, L797 (1983).
- (24) J. Bonča and J. E. Gubernatis, Phys. Rev. B 47, 13137 (1993).
- (25) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, and O. Parcollet, Rev. Mod. Phys. 78 865 (2006).
- (26) P. Werner and A. J. Millis, Phys. Rev. B 74, 155107 (2006).
- (27) A. Liebsch, Phys. Rev. Lett. 95, 116402 (2005).
- (28) K. Inaba, A. Koga, S, Suga, and N. Kawakami, J. Phys. Soc. Jpn. 74, 2393 (2005).
- (29) M. R. Galpin, D. E. Logan, and H. R. Krishnamurthy, Phys. Rev. Lett. 94, 186406 (2005).
- (30) R. Bulla, T. Costi, and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
- (31) Y. Nishikawa, D. J. G. Crow, and A. C. Hewson, e-print arXiv:1005.5113.
- (32) V. Drchal, V. Janiš, J. Kudrnovský, V. S. Oudovenko, X. Dai, K. Haule, and G. Kotliar, J. Phys.: Condens. Matter 17, 61 (2005).
- (33) M. R. Galpin, A. B. Gilbert, and D. E. Logan, J. Phys.: Condens. Matter 21, 375602 (2009).
- (34) A. Kauch and K. Byczuk, e-print arXiv:0912.4278.
- (35) V. Janiš, Phys. Rev. B 60, 11345 (1999).
- (36) G. Baym, Phys. Rev. 127, 1391 (1962).
- (37) V. Janiš and M. Ringel, Acta Phys. Polon. A 115, 30 (2009).
- (38) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
- (39) R. Bulla, T. A. Costi, and D. Vollhardt, Phys. Rev. B 64, 045103 (2001).
- (40) P. Werner and A. J. Millis, Phys. Rev. Lett. 99, 126405 (2007).
- (41) N. E. Bickers, in Theoretical Methods for Strongly Correlated Electrons, edited by D. Senechal, A. Tremblay, and C. Bourbonnais (Springer-Verlag, New York, 2004), p. 237.
- (42) S. X. Yang at al, Phys. Rev. E 80, 046706 (2009).