Electronic orders in multi-orbital Hubbard models with lifted orbital degeneracy
aa r X i v : . [ c ond - m a t . s t r- e l ] F e b Electronic orders in multi-orbital Hubbard models with lifted orbital degeneracy
Shintaro Hoshino and Philipp Werner Department of Basic Science, The University of Tokyo, Meguro, Tokyo 153-8902, Japan Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland (Dated: October 15, 2018)We study the symmetry-broken phases in two- and three-orbital Hubbard models with liftedorbital degeneracy using dynamical mean field theory. On the technical level, we explain howsymmetry relations can be exploited to measure the four-point correlation functions needed forthe calculation of the lattice susceptibilities. In the half-filled two-orbital model with crystal fieldsplitting, we find an instability of the metallic phase to spin-orbital order with neither spin nororbital moment. This ordered phase is shown to be related to the recently discovered fluctuating-moment induced spin-triplet superconducting state in the orbitally degenerate model with shiftedchemical potential. In the three-orbital case, we consider the effect of a crystal field splitting onthe spin-triplet superconducting state in the model with positive Hund coupling, and the spin-singlet superconducting state in the case of negative Hund coupling. It is demonstrated that forcertain crystal field splittings the higher energy orbitals instead of the lower ones are relevant forsuperconductivity, and that T c can be enhanced by the crystal field effect. We comment on theimplications of our results for the superconductivity in strontium ruthenates, and for the recentlyreported light-enhanced superconducting state in alkali-doped fullerides. I. INTRODUCTION
The local Slater-Kanamori interaction [1, 2] originatingfrom Coulomb repulsion leads to highly nontrivial phasediagrams and crossover phenomena in multi-orbital Hub-bard models. Depending on the filling and the energysplittings between the orbitals one finds antiferromag-netic or ferromagnetic order [3, 4], orbital order [3], high-spin/low-spin transitions [5, 6], staggered high-spin/low-spin order [7], excitonic insulating phases [8, 9], or intra-orbital spin-triplet (equal-spin) pairing [4]. Some of theseinstabilities are linked to the spin-freezing crossover [10],which occurs in models with nonzero (positive) Hundcoupling, and underlies the unusual finite-temperatureproperties of Hund metals [11]. In models with nega-tive Hund coupling, an intra-orbital spin-singlet super-conducting state appears [12, 13], but alternative orderedstates have not yet been systematically explored.An unbiased way of mapping out the electronic insta-bilities to long-range orders is to compute the correspond-ing susceptibilities and to look for divergences as a func-tion of the model parameters. While such calculationscannot be easily performed with numerically exact latticemethods, they become computationally tractable withinthe dynamical mean field theory (DMFT) approximation[14]. This theory assumes a spatially local self-energy andvertex, and produces qualitatively correct solutions forhigh-dimensional lattice models. The use of a featurelessdensity of states in DMFT ensures that the results do notdepend on subtle bandstructure effects. In this work, weemploy the DMFT formalism and a semi-circular den-sity of states to map out the instabilities to uniform andstaggered long-range ordered phases in two and three or-bital systems, focusing on the physically most interestingintermediate coupling regime.One main purpose of this study is to show how sym-metry relations can be exploited in the numerical sim- ulations based on continuous-time Monte Carlo impu-rity solvers [15], and in the analysis of the ordering in-stabilities. In particular, we will demonstrate that thespin-orbital ordered phase (also called excitonic insula-tor phase [8, 9]) appearing in half-filled two-orbital mod-els with crystal field (CF) splitting is related by symme-try operations to the orbital-singlet spin-triplet super-conducting state found in the orbitally degenerate modelaway from half-filling [16–21]. We also use symmetry re-lations to analyze and explain the Cooper pair formationin three orbital systems with lifted orbital degeneracy.For positive (ferromagnetic) Hund coupling, the Slater-Kanamori interaction describes t -based three orbitalsystems, as realized for example in Sr RuO and SrRuO [11]. These materials indeed exhibit spin-triplet super-conductivity or ferromagnetism, as well as the non-Fermiliquid properties associated with spin-freezing. Recentlythe authors have demonstrated that the spin-triplet pair-ing is induced by the fluctuating local moments which ap-pear at the border of the spin-frozen regime and eventu-ally order ferromagnetically or antiferromagnetically [4].Similar multi-orbital physics is also expected to be rele-vant for other types of spin-triplet pairing, such as real-ized in U-based ferromagnetic superconductors [22]. Inthis paper we will examine the effect of CF splittings onthe spin-triplet pairing. We will show that the Cooperpairs can be formed in the higher-energy lifted orbitalsinstead of lower one. This somewhat counter-intuitive re-sult can be explained by exploiting a particle-hole trans-formation.The half-filled three-orbital model with negative Hundcoupling exhibits an intra-orbital pairing state which isrelevant for fulleride superconductors [12, 23, 24]. Theeffect of CF splittings on this superconducting state isof interest in connection with the recently reported light-enhanced superconductivity in K C [25]. We show thatin equilibrium, the lifting of the orbital degeneracy canslightly stabilize the pairing, but the effect on T c is small.Hence, it is unlikely that the experimentally observeddramatic enhancement of T c can be explained by quasi-static distortions of the C molecules within an equilib-rium picture.The outline of the paper is as follows. In Sec. II, wepresent the model and explain its basic properties, inparticular the symmetries. We discuss in Sec. III sometechnical details concerning the calculation of lattice sus-ceptibilities within DMFT. Section IV shows the resultsfor the two-orbital model, and explains the connectionbetween spin-orbital order and spin-triplet superconduc-tivity, while Sec. V presents results for the three-orbitalmodel with lifted orbital degeneracy. Section VI is a briefsummary and conclusion. II. MODEL AND METHOD
We consider the Hubbard model with M degenerateorbitals given by the general Hamiltonian H = X k γσ ( ε k − µ ) c † k γσ c k γσ + U X iγ n iγ ↑ n iγ ↓ + U ′ X iσ,γ<γ ′ n iγσ n iγ ′ ¯ σ + ( U ′ − J ) X iσ,γ<γ ′ n iγσ n iγ ′ σ − αJ X i,γ<γ ′ ( c † iγ ↑ c iγ ↓ c † iγ ′ ↓ c iγ ′ ↑ + c † iγ ↑ c † iγ ↓ c iγ ′ ↑ c iγ ′ ↓ + H . c . ) . (1)The operator c iγσ annihilates the electron with orbital γ and spin σ at site i . The Fourier transformation isdefined by c k γσ = N − / P i c iγσ e i k · R i , where N is thenumber of sites and R i the spatial coordinate at site i .The number operator is defined by n iγσ = c † iγσ c iγσ . InEq. (1), the parameter α is introduced to describe theeffect of spin anisotropy and orbital anisotropy. Physi-cally this anisotropy originates from the spin-orbit cou-pling. We note that for repulsively interacting systems,the pair hopping term is irrelevant. This is because thestrongest intra-orbital repulsion U disfavors the doublyoccupied orbital state and as a result the pair hoppingprocess between orbitals rarely occurs. If the Hund cou-pling J is negative, as in models discussed in connectionwith fulleride superconductors [12, 24], the pair hoppingbecomes relevant, while spin-flip processes are effectivelysuppressed. With an anisotropy parameter α <
1, the ef-fect of these spin-flip and pair-hopping process becomesweaker.Let us first consider the isotropic case with α = 1. Inrotationally invariant systems, where the relation U = U ′ + 2 J holds, the interaction part of the Hamiltonian can be rewritten in the form H int = X iγγ ′ σσ ′ (cid:20) U ′ c † iγσ c iγσ c † iγ ′ σ ′ c iγ ′ σ ′ + J c † iγσ c iγ ′ σ c † iγ ′ σ ′ c iγσ ′ + J c † iγσ c iγ ′ σ c † iγσ ′ c iγ ′ σ ′ (cid:21) , (2)which is suitable for discussing the symmetries of thismodel. Using this expression, we can easily verify thatthe transformation S c iγσ S − = e i θ X σ ′ U σσ ′ X γ ′ V γγ ′ c iγ ′ σ ′ (3)does not change H int , if U and V are 2 × M × M orthogonal matrices, respectively. This invarianceimplies an U(1) × SU(2) × SO( M ) symmetry.On the other hand, in the anisotropic case with α = 0,there are only density-density type interactions. Thenthe symmetry is described by [ U (1)] M where 2 M is thenumber of spin/orbital indices. The relevant symme-try operation is given by S c iγσ S − = e i θ γσ c iγσ for any( γ, σ ).When we consider the degenerate-orbital model in-cluding kinetic terms, the total Hamiltonian is also un-changed by this local transformation. In this case wehave the relation h O i = h S OS − i (4)for arbitrary operators O . This relation will be used inSec. IIIC to derive measurement formulas for the two-particle Green functions.We analyze the multi-orbital Hubbard model withinthe framework of DMFT [14]. This theory becomes ex-act in the limit of infinite dimensions, where only localcorrelations are relevant. Because of the local self-energy,the lattice problem can be mapped onto a multi-orbitalimpurity problem with local interaction terms identi-cal to those of Eq. (1) and noninteracting baths whoseproperties can be encoded by the hybridization functions∆ γσ . We consider the infinite-dimensional Bethe lat-tice, whose non-interacting density of states has a semi-circular shape: ρ ( ε ) = (1 / πt ) √ t − ε with t = 1 the(rescaled) hopping integral. For this lattice, the DMFTself-consistency condition simplifies to ∆ γσ = t G γσ ,with G γσ the impurity Green function. Since the wavevector in the Bethe lattice is ill-defined, it is natural towork in a real space representation. On the other handwe can also consider an alternative description in terms of“pseudo-wave-vectors” (see Appendix A), which enablesa discussion analogous to that for ordinary periodic lat-tices. III. EVALUATION OF SUSCEPTIBILITIESA. Definition
In this paper we discuss the instabilities toward long-ranged orders. We consider ordered states correspondingto operators O ( q ) = X i O i e − i q · R i , (5)which can be detected by the divergence of the suscepti-bilities defined by χ O ( q ) = 1 N Z β h O ( q , τ ) O ( − q ) i d τ. (6)Here O ( τ ) = e τ H O e − τ H is the Heisenberg picture withimaginary time τ . In DMFT, these susceptibilities canbe calculated from the local vertex parts extracted fromthe effective impurity problem [14].Let us consider the specific forms of O by taking thetwo-orbital model as an example. For diagonal orders,the operators O i are given by n i = X γσ c † iγσ c iγσ , (7) s µi = X γσσ ′ c † iγσ σ µσσ ′ c iγσ ′ , (8) τ νi = X γγ ′ σ c † iγσ σ νγγ ′ c iγ ′ σ , (9) o νµi = X γγ ′ σσ ′ c † iγσ σ νγγ ′ σ µσσ ′ c iγ ′ σ ′ , (10)corresponding to charge, spin, orbital and spin-orbitalmoments, respectively. Here µ, ν = x, y, z . For offdiago-nal orders, we have p s µi = X γγ ′ σσ ′ c † iγσ ǫ γγ ′ ( σ µ ǫ ) σσ ′ c † iγ ′ σ ′ + H . c ., (11) p ν s i = X γγ ′ σσ ′ c † iγσ ( σ ν ǫ ) γγ ′ ǫ σσ ′ c † iγ ′ σ ′ + H . c ., (12)which correspond to the orbital-singlet-spin-triplet andorbital-triplet-spin-singlet pairing amplitudes, respec-tively. The anti-symmetric unit tensor is defined by ǫ = i σ y . All of these operators are hermitian.To further classify the diagonal operators given above,we consider the time-reversal operation defined by T =exp( − i π P i s yi / K with complex conjugation operator K . This transforms the electron operator as T c iγσ T − = X σ ′ ǫ σσ ′ c iγσ ′ . (13)Note that T is an antiunitary operator: T z T − = z ∗ for a complex number z . This operation changes themomentum ( k → − k ) and flips the spin state ( ↑→↓ or ↓→↑ ). As expected, the charge is time-reversal even( T n i T − = n i ) and spin is odd ( T s µi T − = − s µi ). Theorbital moments τ xi and τ zi are time-reversal even, while τ yi is odd. The difference between the orbital momentsarises due to the presence of the imaginary unit in the identifier operator time-reversal(i) n even(ii) s x , s y , s z odd(iii) τ x , τ z even(iv) τ y odd(v) o xx , o xy , o xz , o zx , o zy , o zz odd(vi) o yx , o yy , o yz even(vii) p sx , p sy , p sz even(viii) p xs , p zs even(ix) p ys evenTABLE I: Classification of operators relevant to diagonalorders for the rotationally invariant case ( α = 1). Pauli matrix σ y , which gives an additional minus signunder the time-reversal operation.The results are summarized in Tab. I. While in thispaper we focus on a model study, a physical interpre-tation in terms of the doubly degenerate e g orbitals isgiven in Appendix B, where the operator τ y is shown tobe a (generalized) magnetic moment. These argumentscan also be applied to the three orbital model, where thePauli matrix for orbital is replaced by the Gell-Mann ma-trix. In this case we have eight kinds of orbital orders,three of which are time-reversal odd operators.For the pairing state characterized by Eqs. (11) and(12), the time-reversal symmetry is not broken. Morespecifically, while the time-reversal operation T canchange the sign of these operators, this can be absorbedby simultaneously performing a global U(1) gauge trans-formation. On the other hand, if two or more Cooperpairs become finite with different phases, then the result-ing state has broken time-reversal symmetry [26]. In thiscase the sign cannot be absorbed by the gauge transfor-mation, and thus the time-reversal symmetry is broken. B. Two-particle Green function and susceptibilities
Here we discuss the calculation of the susceptibility inthe framework of DMFT [14]. For diagonal and offdiago-nal orders, the relevant two-particle Green functions aregiven by χ diag ij,a a a a ( τ , τ , τ , τ ) = h T τ c † ia ( τ ) c ia ( τ ) c † ja ( τ ) c ja ( τ ) i− h T τ c † ia ( τ ) c ia ( τ ) ih T τ c † ja ( τ ) c ja ( τ ) i , (14) χ offd ij,a a a a ( τ , τ , τ , τ ) = h T τ c † ia ( τ ) c † ia ( τ ) c ja ( τ ) c ja ( τ ) i , (15)respectively. The index a denotes the spin and orbitalindices ( γ, σ ). Here T τ represents the imaginary time-ordering operator, and we subtract the disconnected partin Eq. (14). We also define the Fourier transformationwith respect to imaginary time by χ ξij,a a a a (i ε n , i ε n ′ ) =1 β Z β d τ d τ d τ d τ χ ξij,a a a a ( τ , τ , τ , τ ) × e i ε n ( τ − τ ) e i ε n ′ ( τ − τ ) , (16)where ξ means ‘diag’ or ‘offd’. Since we are interestedonly in static susceptibilities, we set the bosonic fre-quency to zero. The susceptibilities defined in Eq. (6)can be derived by summing up the fermionic Matsubarafrequencies in the two-particle Green function, which isexplicitly written in the form χ ξη ( q ) = TN X nn ′ X ij X a a a a A ηa a A ηa a e − i q · ( R i − R j ) × χ ξij,a a a a (i ε n , i ε n ′ ) , (17)where the A ηaa ′ are form factors originating from Eqs. (7–10) for ξ =‘diag’, and from Eqs. (11,12) for ξ =‘offd’.Now we employ the Bethe-Salpeter equation which re-lates the two-particle Green function and vertex. Sincethe vertex part Γ is local in DMFT [14], it can be evalu-ated from the local two-particle Green function as χ ξii,a a a a (i ε n , i ε n ′ ) = χ ξ ii,a a a a (i ε n , i ε n ′ )+ X n n X aa ′ a ′′ a ′′′ χ ξ ii,a a aa ′ (i ε n , i ε n )Γ i,a ′ aa ′′′ a ′′ (i ε n , i ε n ) × χ ξii,a ′′ a ′′′ a a (i ε n , i ε n ′ ) . (18)Here the two-particle Green function without vertexparts is written as χ ξ . The site index for the vertexpart can be neglected if the system is uniform.The local vertex extracted above is inserted into thenon-local Bethe-Salpeter equation χ ξij,a a a a (i ε n , i ε n ′ ) = χ ξ ij,a a a a (i ε n , i ε n ′ )+ X ℓn n X aa ′ a ′′ a ′′′ χ ξ iℓ,a a aa ′ (i ε n , i ε n )Γ ℓ,a ′ aa ′′′ a ′′ (i ε n , i ε n ) × χ ξℓj,a ′′ a ′′′ a a (i ε n , i ε n ′ ) . (19)In practice, it is convenient to perform the Fourier trans-formation with respect to the site index before solvingthe matrix equation. Thus we obtain the two-particlelattice Green functions, which contain the information ofthe susceptibilities. C. Measurement trick for local two-particle Greenfunctions
As explained in the previous section, the local verticescan be extracted from the local two-particle Green func-tions. In DMFT calculations of the multi-orbital Hub-bard model (1), we must consider the following two-particle Green functions of the corresponding impurity model: χ a γγ ′ σσ ′ = h T τ c † γσ ( τ ) c γσ ( τ ) c † γ ′ σ ′ ( τ ) c γ ′ σ ′ ( τ ) i , (20) χ b γγ ′ σσ ′ = h T τ c † γσ ( τ ) c γσ ′ ( τ ) c † γ ′ σ ′ ( τ ) c γ ′ σ ( τ ) i , (21) χ c γγ ′ σσ ′ = h T τ c † γσ ( τ ) c γ ′ σ ( τ ) c † γσ ′ ( τ ) c γ ′ σ ′ ( τ ) i . (22)All of the diagonal and offdiagonal two-particle Greenfunctions that appear in the previous subsection can becalculated from χ a , b , c . A powerful method to solve im-purity models is the continuous-time Monte Carlo tech-nique [15]. Here, we employ the hybridization expansionmethod in the matrix formulation [27]. In this approach,the function χ a can be calculated by the standard pro-cedure which involves the removal of two hybridizationfunctions [15]. However, the other functions χ b , c , whichare in general finite, cannot be obtained by the standardtechnique unless the impurity model has offdiagonal hy-bridizations.If we have the continuous SU(2) × SO( M ) symmetrythe two-particle Green functions χ b and χ c can be ex-pressed in terms of χ a using Eq. (4). To see this, weconsider the function χ b12 ↑↓ in the two-orbital model( M = 2). The SU(2) and SO(2) matrices in the spinand orbital spaces are given by U = z − w ∗ w z ∗ ! , V = x − yy x ! , (23)where | z | + | w | = 1 ( z, w ∈ C ) and x + y = 1( x, y ∈ R ). The SO(2) matrix in orbital space repre-sents a rotation around the τ y axis. The quantity χ b12 ↑↓ is transformed by U as χ b12 ↑↓ −→ ( | z | + | w | ) χ b12 ↑↓ + 2 | z | | w | ( χ a12 ↑↑ − χ a12 ↑↓ ) , (24)where we have used the equivalence between γ = 1 , σ = ↑ , ↓ . If the transformation does not changethe whole Hamiltonian, which implies α = 1, the left- andright-hand sides of Eq. (24) can be regarded as identical.Thus we obtain the relation χ b12 ↑↓ = χ a12 ↑↑ − χ a12 ↑↓ . (25)Note that this result is not dependent on the values of z and w . This relation has been previously used in thestudy of a two-channel Kondo lattice [28]. In a similarmanner, if the system has orbital SO(2) symmetry, itsatisfies the additional relation χ c12 ↑↓ = χ a11 ↑↓ − χ a12 ↑↓ + ˜ χ b12 ↑↓ , (26)where we have defined the time-sorted quantity˜ χ ( τ , τ , τ , τ ) = − χ ( τ , τ , τ , τ ) . (27)By an analogous trick we can also obtain the relation χ c12 ↑↑ = χ a11 ↑↑ − χ a12 ↑↑ + ˜ χ a12 ↑↑ . (28) -2 0 2 4 6 8 10 12 14 0.6 0.8 1 1.2 1.4 1.6 1.8 2 i nv e r s e s u s ce p ti b ilit y [ a r b . un it ] n (ii)(iii)(iv)(v)(vi) FIG. 1: (Color online) Filling dependence of inverse stag-gered susceptibilities for diagonal orderings in the isotropiccase ( α = 1). The orders (ii)-(vi) are defined in Tab. I. Wechoose the parameters U = 2, J/U = 1 / T = 0 .
02. Thestaggered charge susceptibility [(i) in Tab. I] is not shown herebecause of low numerical accuracy.
The other types of χ b and χ c can be represented by thequantities obtained above. Thus, for the calculation ofthe two-particle Green functions χ a , b , c , it is sufficient tomeasure χ a .A key prerequisite for this technique is the SU(2) orSO(2) symmetry in the Hamiltonian. Hence, these for-mulae do not hold when the symmetry is lowered by, e.g.,a CF splitting. The generation of new terms on the right-hand side of Eq. (24) is also an important ingredient forthe present trick. For example the U(1) transformationonly changes the phase factor, which is not enough toderive expressions for χ b and χ c . Thus, at least a two-dimensional representation of the symmetry operation isnecessary.Let us also add a comment on the phase transfor-mation. While usually it changes the phase of thefermion operators uniformly ( c iγσ → c iγσ e i θ ), here wecan slightly generalize it as c iγσ → c iγσ e i( θ σ + θ γ ) becauseof the SU(2) × SO(2) symmetry. For the invariance ofthe Hamiltonian, the phases θ ↑ and θ ↓ are arbitrary, but θ − θ = πm ( m ∈ Z ) must be satisfied due to the pair-hopping term. These are related to conservation laws ofspin/orbital indices. Combining this transformation withEq. (4), it can be explicitly shown that expectation valuessuch as h c † ↑ c ↓ c † ↑ c ↓ i and h c † ↑ c ↑ c † ↑ c ↑ i are identicallyzero.The above techniques can also be used in models witha general number M of orbitals. Namely, if we choosetwo of the M orbitals, we can utilize the partial SO(2) symmetry within this subspace, so that the procedureremains unchanged. IV. TWO-ORBITAL MODELA. Degenerate orbitals
In this section, we study the ordering instabilities inthe two-orbital model with
J > Q ” (see Appendix A) in the spin-isotropic systemwith α = 1. Here the symmetry tricks explained in theprevious section are used for the evaluation of the two-particle Green functions. The uniform susceptibilities fordiagonal and offdiagonal orders do not diverge for thechosen parameters, and therefore are not shown.Near half filling ( n . s ( Q ) appears as demonstrated by the sign changeof the corresponding inverse susceptibility in Fig. 1. Theother susceptibilities remain positive and hence no fur-ther instabilities are present. For the low-filling casewith n .
1, all of the susceptibilities behave in a similarmanner. This is because the interaction effect is weakerand an approximate SU(4) symmetry appears. For thepresent choice of parameters no other symmetry break-ing, including superconductivity, occurs.In our study, although the survey is limited due to thehigh numerical cost for α = 0, we did not find any inter-esting ordered states for the spin isotropic case ( α = 1)with J >
0. In the following, we therefore focus on themodel with anisotropic interaction ( α = 0). The corre-sponding phase diagram is richer, since superconductiv-ity appears in addition to trivial magnetic orders. Be-cause the spin-flip and pair-hopping terms in Eq. (2) aredropped in the spin-anisotropic case with α = 0, we donot have to consider χ b , c and the calculation of the two-particle Green functions simplifies substantially. The in-dependent order parameters are also changed from theprevious section, and are listed in Tab. II. We will con- identifier operator(i’) n (ii’) s z (iii’) s x , s y , o zx , o zy (iv’) τ z (v’) τ x , τ y , o xz , o yz (vi’) o zz (vii’) o xx , o xy , o yx , o yy (viii’) p sx , p sy (ix’) p zs , p sz (x’) p xs , p ys TABLE II: Classification of independent operators for theIsing-anisotropic case ( α = 0) without CF splitting. J / U U FL MI
FIG. 2: (Color online) Phase boundaries determined by thedivergent points of susceptibilities at n = 2 and T = 0 . J/U = 1 / J/U = 1 / J/U = 2 / sider only uniform ( q = ) and staggered ( q = Q ) order-ing vectors.In the following of this section, we focus on s z , o xx , and p s x , since the corresponding symmetry breaking appearsin regions which are not dominated by other ordering,for the parameters considered here. Figure 2 shows thephase diagram in the plane of U and J/U at half filling( n = 2) and T = 0 .
02. For J = 0, because of the highsymmetry SU(4), all the diagonal orders except for thecharge order are degenerate. On the other hand, a finite J substantially stabilizes the magnetic order with s z ( Q ).At the same time, the critical value of the metal-insulatortransition is reduced by the Hund coupling. There ap-pears staggered spin-orbital order ( o xx ( Q )) and orbital-singlet/spin-triplet pairing ( p s x ( )), although these or-dered regions are covered by the s z ( Q ) order. The spin-orbital order realized at J = 0 is destabilized by a smallHund coupling J . This is because the c † ↑ c ↓ operatoris relevant for the spin-orbital order with o xx , and thisspin-orbital flipping process is suppressed by a positive J which favors equal-spin states. However, it is stabilizednear the Mott transition due to local spin and chargefluctuations. As will be shown later, these orders becomemost stable away from half-filling or in the presence of aCF splitting.It is notable that at J/U = 0 . o xx ( Q ) U n FL spin-freezingcrossover NFL
FIG. 3: (Color online) Filling ( n ) versus interaction ( U ) phasediagram of the two-orbital model with J/U = 1 / T =0 .
02. The Fermi liquid (FL) and non-Fermi liquid (NFL)regimes are separated by the spin-freezing crossover line. and p s x ( ) have the same transition point. This can beunderstood by symmetry considerations. We first definethe particle-hole transformation for orbital 2, P c i σ P − = X σ ′ σ xσσ ′ c † i σ ′ e i Q · R i , (29)which leaves c i σ unchanged. This is similar to the trans-formation from repulsive to attractive interactions in thehalf-filled single-band Hubbard model [29]. One can showthat this transformation does not change the Hamilto-nian with Ising anisotropy ( α = 0) at J/U = 2 / n = 2. (In fact, this invariance holds also for the casewhere the spin-flip term is included, but not in the pres-ence of the pair-hopping term.) On the other hand, P transforms the offdiagonal order into a diagonal order: P p s x ( ) P − = o xx ( Q ) . (30)Thus, these two orders are degenerate at this particularvalue of J/U .Figure 3 plots the phase diagram of the orbitally de-generate model away from half-filling, for T = 0 .
02 and
J/U = 1 /
4. Near half-filling, an antiferromagneticallyordered phase is found, while in the large U regime, aferromagnetic phase extends over a wide filling range.Between these two magnetically ordered regions, we findan instability toward an orbital-singlet spin-triplet su-perconducting state with p s x ( ). Here, the quasi-localCooper pair results from purely repulsive interactions, bya mechanism that has been proposed in the cold atomcontext [30]. This phase extends from the end-pointof the half-filled Mott insulator along the so-called spinfreezing line [10] into the metallic regime, in analogy tothe three orbital results discussed in Ref. 4. The spin-freezing phenomenon is well characterized by the localquantity∆ χ loc = Z β dτ [ h s zi ( τ ) s zi i − h s zi ( β/ s zi i ] . (31)The first term is the local spin susceptibility. Thelong-time spin correlator, which appears as the secondterm, reflects the presence of frozen moments. Thedifference ∆ χ loc thus quantifies the local spin fluctua-tions, whose maxima in parameter space define the spin-freezing crossover line [4]. This line is closely relatedto the spin-triplet superconductivity as shown in Fig. 3.While the main topic of this paper is the effect of CFsplitting, Figs. 2 and 3 are important for a deeper under-standing of the resulting ordered states under the CF, asdiscussed in the following. B. Split orbitals
We next consider the half-filled two-orbital model withan additional crystal field (CF) splitting given by H CF = ∆ X iσ ( n i σ − n i σ ) . (32)For α = 1, the metal-insulator phase diagram of thismodel without ordering has been studied in Ref. 5. Aqualitatively similar diagram is obtained also for the α = 0 case [31]. At large U and small ∆, there is ahigh-spin (Mott) insulating phase with small orbital po-larization. For sufficiently large ∆, the system is in alow-spin insulating state with large orbital polarization,which for U → U and finite T , there is atransition from the high-spin to the low-spin insulator at∆ ≃ J/ U , there exists an intermediatemetallic phase, which for large J/U and low T extendsalong the ∆ ≃ J/ U phase diagram for theanisotropic case ( α = 0) at half filling. While thehigh-spin insulator is antiferromagnetically ordered, themetallic nose separating the high-spin and low-spin insu-lators is unstable to the staggered spin-orbital orderingwith o xx ( Q ). This result is consistent with the excitonicinsulator state found in this parameter regime in Ref. 8.The low-spin insulating region is not susceptible to long-range order. We also find an instability toward the stag-gered high-spin/low-spin order reported in Ref. 7. Thisinstability, which can be detected by the orbital suscep-tibility with τ z ( Q ), appears for U &
6. However, we donot show this phase, since it is covered by the o xx ( Q )order at least for U ≤ U ∆ FL LSIHSMI spin-freezingcrossover
FIG. 4: (Color online) Phase diagram of the half-filled two-orbital model with
J/U = 1 / U at temperature T = 0 .
02. In the disordered case,we have Fermi liquid (FL), high-spin Mott insulator (HSMI),and low-spin insulator (LSI) phases.
We now discuss the connections between the two phasediagrams in Figs. 3 and 4, which follow from the sym-metry relation (30). Away from half-filling, in themodel with CF splitting, the spin-triplet superconductiv-ity p s x ( ) is stabilized along the spin-freezing crossoverline, which indicates that local spin-fluctuations inducethe pairing state [4]. If we perform the particle-holetransformation, this pairing state changes into the or-dered state with o xx ( Q ) as discussed above. At the sametime, the chemical potential term H chem = − µ P iγσ n iγσ is transformed into a CF splitting term: PH chem P − = − H CF (33)with ∆ = µ . In other words, the spin-triplet supercon-ductivity away from half filling corresponds to the spin-orbital order under the CF splitting . Furthermore the lo-cal spin moment is unchanged: P s zi P − = s zi . Hence, inboth cases the local spin-fluctuations play an importantrole in stabilizing the ordered state. As shown in Fig. 4,there is indeed a spin-freezing line extending along themetallic nose in the region where spin-orbital order ap-pears. While this mapping is exact only for J/U = 2 / J/U = 1 / ∆0 J > 0 J < 0 γ = 3γ = 1,2 FIG. 5: (Color online) Schematic illustration of atomic con-figurations as a function of ∆ at half-filling ( n = 3). Thecircles with dotted lines show the orbital which is responsiblefor superconductivity. For the case with n = 2, one of threeelectrons should be removed from figure. V. THREE-ORBITAL MODELA. General remarks
We are interested in the stability of the inter-orbital-spin-triplet (
J >
0) and intra-orbital-spin-singlet (
J <
0) superconducting phases to CF splittings of the “2/1”type H CF = − ∆ X iσ ( n i σ + n i σ − n i σ ) . (34)This term breaks the orbital SO(3) symmetry, but theSO(2) symmetry within the γ = 1 , α = 1.At half-filling and without CF splitting, the Hamilto-nian is invariant under the simple particle-hole transfor-mation defined by P c iγσ P − = c † iγσ e i Q · R i , (35)i.e., we have P H P − = H . On the other hand,the CF Hamiltonian is transformed as P H CF P − = − H CF . Thus the sign of the CF splitting ∆ is reversed.For some interaction and filling parameters, as seen be-low, the present system shows superconductivity withinter-orbital spin-triplet pairs p γγ ′ t ( ) = X i c † iγ ↑ c † iγ ′ ↑ + H . c . (36)for J > γ = γ ′ , or intra-orbital spin-singlet pairs p γ s ( ) = X i c † iγ ↑ c † iγ ↓ + H . c . (37)for J <
0. Here we choose these expressions for thepair amplitudes, although we can classify them usinga similar method as given in Eqs. (11) and (12) usingPauli and Gell-Mann matrices. We have the symmetryrelations P p γγ ′ t P − = − p γγ ′ t and P p γ s P − = − p γ s .While the signs are reversed after the transformation, i.e., the phases of the pair amplitudes are rotated by π ,their forms are unchanged after the particle-hole trans-formation.To understand the implications of the above symmetryrelations, let us first consider the J > p , p and p are degenerate. For ∆ >
0, asillustrated in the top part of Fig. 5, one expects that p is more stable than p and p because the orbitals with γ = 1 , γ = 3. In a similarmanner one may also speculate that for ∆ < p and p are more stable than p . However this is in-correct . The pair p is the most stable even for ∆ < > <
0, according to the above particle-holesymmetry argument, even though the Cooper pairs areformed in the higher-energy orbital (broken circle in theupper panel of Fig. 5) when ∆ <
0. This at first sightcounter-intuitive result can be rationalized by consider-ing the degeneracy: the p and p pairs are destabi-lized by fluctuations among these energetically degener-ate states, which may diminish the pairing compared tothe non-degenerate case. On the other hand, the p pairs are not subject to such fluctuations and are there-fore more stable.One can also think of other effects which help explainthis peculiar behavior. Whereas the lower-energy orbitalis easily occupied and pairs between low-energy and high-energy electrons can be formed, the mobility of thesepairs, which is also important for realizing superconduc-tivity, is reduced because of the high occupancy of thelow-energy orbitals. On the other hand, for higher-energyorbitals the mobility of the pairs is high, although thepairs themselves are harder to form. As a result, thesuperconducting state resulting from pairs in the highenergy orbitals is actually more robust, which is also sup-ported by the numerical results.The same argument can also be applied to the caseof J <
0. At ∆ = 0, the three pairs p , p and p aredegenerate. For ∆ <
0, as illustrated in the lower panelof Fig. 5, one expects that p is more stable than p and p . On the other hand, for ∆ > p and p are more stable than p , but this is againincorrect. Although the level of the orbital γ = 3 ishigher than the others, the most stable Cooper pair isformed in this high-energy orbital.Note that in the above discussion, we do not invokeany subtleties of the system, only the particle-hole sym-metry of the original Hamiltonian without CF splitting.Keeping these facts in mind, we can readily understandthe result in the following subsections. B. J > case In the remainder of this section, we show results forthe three-orbital model with CF splitting, focusing again ∆ FL spin-freezingcrossover (a) NFL U T (b) U =3.2 FIG. 6: (Color online) (a) Phase boundaries determined bythe divergence of susceptibilities for the parameters n = 2, J/U = 1 / T = 0 .
02. The limit ∆ → + ∞ correspondsto the two-orbital model at half filling. (b) Phase diagram inthe temperature ( T )-CF splitting (∆) space for U = 3 . on the case with Ising anisotropy ( α = 0). For J >
0, thepairing amplitude is given by Eq. (36), and the magneticmoments for the three-orbital model are defined by s z ( q ) = X iγσ c † iγσ σ zσσ c iγσ e − i q · R i . (38)Other ordered states are not considered here since theseare covered by the orders p γγ ′ t ( ) and s z ( q ) in the pa-rameter range considered in this paper.At ∆ = 0 the phase diagram for the three-orbitalmodel [4] has properties similar to those of the two-orbital model shown in Fig. 3. Here we discuss the effectof the CF splitting on the three-orbital model with n = 2,because it is relevant to Sr RuO as discussed later. Fig-ure 6(a) shows the U -∆ phase diagram at n = 2 and T = 0 .
02. Note that we do not have a symmetry between∆ > < s z ( Q )) dominates for suf-ficiently large ∆ because in the limit ∆ → ∞ the systembecomes the half-filled two-orbital model discussed in theprevious section. In the large- U region the ferromagneticorder ( s z ( )) appears.The inter-orbital spin-triplet superconductivity is re-alized at intermediate values of U and for a large range of ∆. The spin-freezing crossover defined by the maxi-mum of Eq. (31) is also plotted, and it is obvious thatthe superconductivity is stabilized near this line. For∆ = 0, the degeneracy of the three pairs is lifted, andthe most stable one is p regardless of the sign of ∆, asexplained in the previous subsection. The other pairings p , (filled circles in Fig. 6) are quickly suppressed.The CF splitting dependence of the transition tem-perature is shown in Fig. 6(b). It is notable that thetransition temperature is enhanced by the CF splittingfor ∆ >
0. This behavior might be due to a reductionof fluctuations among orbitals and an enhanced proba-bility of pair formations by the splitting of the degen-erate orbitals. For ∆ < p inter-orbital spin-triplet superconductivity can beobserved in a wide region of | ∆ | at low temperatures.Our simulations are relevant for strontium ruthenatecompounds. In Sr RuO , the filling is n = 4, which isidentical to n = 2 due to particle-hole symmetry, andthe interaction is estimated as U ≃ . orbitals hasthe form of Eq. (34) due to the tetragonal symmetry ofthe crystal. According to Figs. 6, the spin-triplet super-conductivity occurs at low temperatures even for finiteCF splittings in the physically reasonable range. Thepairing is most robust for same-spin electrons in the de-generate orbitals. We thus obtain consistent results forthis strontium ruthenate compound, although for a moredetailed comparison with experiments we should considerthe realistic band structure.Here we briefly comment on the effect of α which re-stores the isotropy. For the present spin-triplet super-conductivity with J >
0, the effect of spin-flip terms isrelevant and suppresses this pairing due to the fluctua-tions among the three spin-triplet states [4]. Indeed wecould not find the spin-triplet superconductivity down tothe lowest accessible temperatures for α = 1. Thus theanisotropy in spin space is important for its realization.Despite this fact, the spin-triplet superconducting statecan be realized in Sr RuO , since only a small anisotropyis necessary for its stabilization [4]. C. J < case Because of the relevance for fulleride compounds, it isalso interesting to consider the half-filled ( n = 3) threeorbital model with negative Hund coupling [24]. In addi-tion to the magnetic moment defined in the last subsec-tion, we consider two kinds of orbital moments definedby τ , ( q ) = X iγσ c † iγσ λ , γγ c iγσ e − i q · R i , (39)0 ∆ FLMI (a)
HOUMLOOM U (b) U =4.5 T FIG. 7: (Color online) (a) Phase boundaries determined bythe divergence of susceptibilities for n = 3, J/U = − /
10 and T = 0 .
02. The ∆ < > > <
0. Open(filled) pentagons correspond to p s ( p , s ) pairing. (b) Phasediagram in the CF splitting (∆)-temperature ( T ) space for U = 4 . where λ = diag(1 , − ,
0) and λ = diag(1 , , − / √ τ ( Q ) is more stablethan the other, so we plot only the phase boundaries for τ ( Q ).Figure 7(a) shows the U -∆ phase diagram at J/U = − /
10 and T = 0 .
02. Although this ratio of
J/U ismuch larger than the estimate for alkali-doped fullerides[23], it is suitable for clarifying the qualitative behaviorsoriginating from the negative Hund coupling. Withoutany long-range order, we have three kinds of states de-pending on the parameters. In the small and large U regions, metallic and Mott insulating states are realized,respectively. For sufficiently large CF splitting, a newstate appears in-between, which is a metallic state witha completely empty (occupied) orbital for ∆ > < τ ( Q ) moments is thedominant phase near ∆ = 0. For large U , antiferromag-netism is also stabilized, which is explained by the stillremaining active spin degrees of freedom in the half-filledsystem.As shown in Fig. 7, the intra-orbital spin-singlet pair-ing with p is realized in the intermediate U region. Thispairing is also caused by the negative Hund coupling,which favors doubly occupied orbitals. The other pair-ings with p , are less stable (see filled pentagons in thefigure) as discussed before. For small ∆, this supercon-ducting state is covered by the orbital order. With in-creasing CF splitting, we find a region between the or-bital ordered phase and the occupied/unoccupied orbitalmetal, where the p pairing state is most stable. Thetransition temperature as a function of ∆ is plotted inFig. 7(b). It is notable that the transition temperature isslightly enhanced by the CF splitting. This enhancementcan be intuitively interpreted again as resulting from thesuppression of fluctuations among orbitals and increas-ing probability of pair formation. At the same time, for∆ < T c is almost independent of ∆. For ∆ >
0, on the otherhand, the probability of pair formation instead of the mo-bility of pairs is decreased, to give the exactly same resultas in the ∆ < p pairing state disappearswhen orbital γ = 3 becomes fully occupied or empty.Finally we discuss the effect of the parameter α . Inthe J < α >
0, which is in contrast to the
J > J is tiny ( J/U ∼ − .
025 for fullerides[23]), the pairing from negative Hund coupling has thisadvantage and has a chance to be realized in systemswith anisotropic interaction. The effect of the increasingorbital fluctuations on the intra-orbital spin-singlet su-perconducting state for α >
0, however, remains to beinvestigated.
VI. SUMMARY AND DISCUSSION
We have presented a systematic analysis of the elec-tronic ordering instabilities in two- and tree-orbital Hub-bard models. We have used DMFT in combination witha numerically exact hybridization expansion impuritysolver and a semi-circular density of states. Our re-sults thus represent the generic phase diagrams of high-dimensional lattice models, irrespective of details of thebandstructure. Only uniform and staggered order pa-rameters have been considered, and we focused on theintermediate coupling regime, which is relevant for mostunconventional multi-band superconductors and cannotbe accessed reliably by approximate methods such as the1fluctuation-exchange approximation.In the case with positive Hund coupling, we found in-stabilities to antiferromagnetic and ferromagnetic orders,spin-orbital order and orbital-singlet spin-triplet super-conductivity. In the model with negative Hund couplingwe identified antiferro orbital order, antiferromagnetism,and intra-orbital spin-singlet pairing. We have shownhow symmetry relations can be used to connect someof these ordered phases, and to understand the Cooperpairings in three-orbital models with CF splittings.In particular, we showed that in the two-orbital modelwith
J/U = 2 /
5, the orbitally degenerate system awayfrom half-filling can be mapped exactly onto the half-filled system with CF splitting, and that the orbital-singlet spin-triplet superconductivity in the former sys-tem corresponds to spin-orbital order in the latter. Thequalitative correspondence between these at first sightdifferent physical situations can be expected to hold ina wider parameter regime. Since the fluctuating localmoments at the border of the spin-frozen metal regimeplay a crucial role in stabilizing the spin-triplet super-conducting state, and the spin moments are unaffectedby the mapping, this symmetry argument also impliesthat fluctuating moments drive the instability to spin-orbital order. Thus the diagonal and offdiagonal orderscan be understood in a unified manner.Crystal field splittings of the “2/1” type in the threeorbital model can have a surprising effect on the super-conducting state. In the spin-triplet state with filling n = 2 and J >
0, the most stable pairs are always formedin the two degenerate orbitals, even if the energy of theseorbitals is increased by the CF. On the other hand, in thehalf-filled spin-singlet state with
J >
0, the most stablepairs are always formed in the non-degenerate orbital,irrespective of the sign of the CF. We explained this be-havior using a particle-hole transformation, and also gavean intuitive interpretation based on fluctuations of thecondensate and mobility of Cooper pairs.The use of symmetry relations is also interesting froma computational point of view, since the rotational invari-ance of the Slater-Kanamori interaction allows to expresscertain four-point correlation functions in terms of otherswhich can be easily measured in standard hybridizationexpansion Monte Carlo simulations [15]. This trick how-ever fails as soon as the rotational invariance is broken,e.g. by the presence of CF splittings. Whether or notit is possible to extend this technique to more generalsituations by exploiting additional symmetries is an in-teresting open problem. An alternative strategy is todirectly measure all types of four-point correlation func-tions using a worm sampling algorithm [34].Our three-orbital calculations with filling n = 2 (equiv-alent to n = 4) can be regarded as toy model simulationsof Sr RuO . This material is a quasi-2D system withtetragonal symmetry and has a “2/1” type CF splitting.While the CF splitting quickly suppresses the supercon-ductivity in the non-degenerate orbital, the pairing be-tween the degenerate orbitals is robust. Our calculations are consistent with a fluctuating-moment induced spin-triplet superconductivity in this compound, which has apartial orbital degeneracy ( zx and yz orbitals). More ac-curate simulations would need to take into account therealistic bandstructure.The half-filled three-orbital model with negative Hundcoupling is of interest in connection with fulleride com-pounds, where the Jahn-Teller screening of the small J leads to a stabilization of the low-spin states [12]. The ef-fect of CF splitting is to further stabilize the spin-singletpairing in the non-degenerate orbital, and to suppress acompeting orbital ordered phase. While the J parame-ters used in our study are larger than the ab-initio esti-mates [23], this finding has some implications on the re-cently reported light-induced superconductivity in K C [25]. One possible explanation for the enhanced T c , pro-posed in Ref. [25], is that the driving of a phonon of T u symmetry leads to an essentially static distortion of theC molecules, and hence a splitting of the initially de-generate molecular orbitals. Our result shows that thissplitting, if studied within an equilibrium formalism, cancontribute to a slightly enhanced pairing. The relatedeffect of orbital differentiation in the interaction param-eters still needs to be systematically explored. Given thelarge enhancement of T c observed in the experiments,one may speculate that the J parameter is effectivelymodified in the driven state (enhanced dynamical Jahn-Teller effect). This, and other non-equilibrium phenom-ena should be investigated within a Floquet formalism orusing the non-equilibrium extension of DMFT [35]. Acknowledgments
S.H. acknowledges financial support from JSPS KAK-ENHI Grant No. 13J07701 and P.W. support from FP7ERC starting grant No. 278023. The authors benefitedfrom the Japan-Swiss Young Researcher Exchange Pro-gram 2014 coordinated by JSPS and SERI. The numeri-cal calculations have been performed on the BEO04 clus-ter at the University of Fribourg and the supercomputerat ISSP (University of Tokyo).
Appendix A: Bethe Lattice in Infinite Dimensions
Here we consider some properties of the Bethe lattice ininfinite dimension. Since the wave vectors are ill-definedin this special lattice, we have to deal with real space.However, we will show that the concept of wave vectorscan be partially applied, and that a treatment similar tothat of ordinary lattices is possible.2
1. Single-Particle Green Function
Let us here consider the tight-binding model for non-interacting spinless fermions on the Bethe lattice: H = − t X h ij i c † i c j , (A1)where the summation with respect to the site index isover the nearest neighbor pairs. We begin with the equa-tions of motion given by − ∂ τ c i = − t X δ c i + δ , (A2) − ∂ τ c i + δ = − tc i − t X δ + δ ′ =0 c i + δ + δ ′ , (A3)where − ∂ τ O = [ O, H ]. The site index i + δ denotes thenearest neighbor sites of the site i . The Fourier transfor-mation of the equations of motion for the Green functiondefined by G ij ( τ ) = −h T τ c i ( τ ) c † j i is given by zG ( z ) = 1 − tdG ( z ) , (A4) zG k ( z ) = − tG k − ( z ) − tdG k +1 ( z ) , (A5)for k = 1 , , · · · , which represents the number of sitesbetween the positions at i and j . Here z is a complexfrequency, d is a connectivity, and G is a local Greenfunction. We have taken the limit d → ∞ . Defining α k = G k /G k − , one can show the following relation α k = − tz + tdα k +1 = − tz + td − tz + td − tz + ··· = − tz + tdα k . (A6)Hence the ratio α k is independent of k , and we can writeit as α k = α . This quantity is explicitly derived by solv-ing the quadratic equation and we obtain the Green func-tion as G k ( z ) = − α ( z ) k +1 /t. (A7)One of the two solutions for α is chosen so that it behavesas G ( z ) → /z when | z | → ∞ . The form of the Greenfunction is the same as the one obtained in Ref. 36. Forinteracting systems, the local self-energy Σ is included bythe replacement z −→ z − Σ.
2. Two-Particle Green Function
We define the two-particle Green functions withoutvertex parts by χ ( z ) = − N X ij G ij ( z ) G ji ( z ) , (A8) χ ( z ) = − N X ij λ i λ j G ij ( z ) G ji ( z ) , (A9) which are relevant to uniform and staggered susceptibil-ities, respectively. Here λ i = ± χ = − ( α/t ) − dα , χ = − ( α/t ) dα . (A10)This leads to the simpler expressions χ = d G d z , χ = − G z . (A11)Thus we derive the uniform and staggered componentsof two-particle Green functions.We now consider another representation of the abovetwo-particle Green functions, and introduce a density ofstates and a “wave-vector k summation” by G ( z ) = Z ρ ( ε ) z − ε d ε ≡ N X k g k ( z ) , (A12)where ρ ( − ε ) = ρ ( ε ), g k ( z ) = 1 / ( z − ε k ) and N = P k ε k ,and only assume the existence of a vector Q defined by ε k + Q = − ε k . We also introduce the q -dependent two-particle Green function by˜ χ q ( z ) = − N X k g k ( z ) g k + q ( z ) . (A13)Using Eq. (A12), one can show the relations ˜ χ = χ and ˜ χ Q = χ . Thus the pseudo wave vectors q = and q = Q correspond to the uniform and staggeredcomponents, respectively. Namely, the Bethe lattice canbe partly handled as if it were a simple cubic lattice whichhas the same relation ε k + Q = − ε k with the staggeredordering vector Q = ( π, π, π ). Note that this analogy isvalid only for the uniform and staggered components. Appendix B: Physical Interpretation of OrderParameters
Let us physically interpret the orbital and spin-orbitalmoments in terms of doubly degenerate e g orbitals for d electrons under the cubic symmetry. The wave functionsare written as | γ = 1 i = | i , (B1) | γ = 2 i = ( | i + | − i ) / √ , (B2)where | m i in the right-hand side is the eigenstate of theangular momentum ℓ z . With use of this expression, wecan show the relations τ z ∝ ℓ z − ℓ , (B3) τ x ∝ ℓ x − ℓ y , (B4) τ y ∝ ℓ x ℓ y ℓ z , (B5)3where the overline symmetrizes the product of operatorsas e.g. ℓ x ℓ y = ( ℓ x ℓ y + ℓ y ℓ x ) / τ z and τ x are rank 2 operators,while τ y is a rank 3 operator. For the spin-orbital mo-ment, we further put the spin moment on the orbitalmoment τ µ , and hence the rank of the operator is in- creased by one. Namely, o zµ and o xµ are rank 3 op-erators, and o yµ is a rank 4 operator. In the contextof f -electron systems, rank 0,1,2,3,4 operators are called‘scalar (monopole)’, ‘dipole’, ‘quadrupole’, ‘octupole’ and‘hexadecapole’, respectively [37]. The odd (even) ranktensor is time-reversal odd (even). [1] J.C. Slater, Quantum Theory of Atomic Structure (McGraw-Hill, New York, 1960).[2] J. Kanamori, Prog. Theor. Phys. , 275 (1963).[3] C.-K. Chan, P. Werner, and A. J. Millis, Phys. Rev. B , 235114 (2009).[4] S. Hoshino and P. Werner, Phys. Rev. Lett. , 247001(2015).[5] P. Werner and A. J. Millis, Phys. Rev. Lett. , 126405(2007).[6] R. Suzuki, T. Watanabe, and S. Ishihara, Phys. Rev. B80, 054410 (2009).[7] J. Kuneˇs and P. Augustinsky, Phys. Rev. B , 115134(2014).[8] J. Kuneˇs, Phys. Rev. B , 235140 (2014).[9] J. Kuneˇs, J. Phys.: Condens. Matter , 333201 (2015).[10] P. Werner, E. Gull, M. Troyer, and A. J. Millis, Phys.Rev. Lett. , 166405 (2008).[11] A. Georges, L. d. Medici, and J. Mravlje, Annual Reviewof Condensed Matter Physics , 137 (2013).[12] M. Capone, M. Fabrizio, C. Castellani, and E. Tosatti,Rev. Mod. Phys. , 943 (2009).[13] A. Koga and P. Werner, Phys. Rev. B , 085108 (2015).[14] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996).[15] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M.Troyer, and P. Werner, Rev. Mod. Phys. , 349 (2011).[16] A. Klejnberg and J. Spalek, J. Phys. Condens. Matter , 6553 (1999).[17] J. Spalek, Phys. Rev. B , 104513 (2001).[18] J.E. Han, Phys. Rev. B , 054513 (2004).[19] S. Sakai, R. Arita, and H. Aoki, Phys. Rev. B , 172504(2004); Physica B - , 554 (2005).[20] K. Kubo, Phys. Rev. B , 224509 (2007).[21] M. Zegrodnik and J. Spalek, New J. Phys. , 033001(2014). [22] D. Aoki and J. Flouquet, J. Phys. Soc. Jpn. , 011003(2012).[23] Y. Nomura, K. Nakamura, and R. Arita, Phys. Rev. B , 155452 (2012).[24] Y. Nomura, S. Sakai, M. Capone, and R. Arita, ScienceAdvances , e1500568 (2015).[25] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser,A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M.Ricco, A. Subedi, S. R. Clark, D. Jaksch, A. Cavalleri,arXiv:1505.04529 (2015).[26] M. Sigrist and K. Ueda, Rev. Mod. Phys. , 239 (1991).[27] P. Werner and A. J. Millis, Phys. Rev. B , 155107(2006).[28] S. Hoshino and Y. Kuramoto, Phys. Rev. Lett. ,167204 (2014).[29] H. Shiba, Prog. Theor. Phys. , 2171 (1972).[30] K. Inaba and S. Suga, Phys. Rev. Lett. , 255301(2012).[31] H. Hafermann, K.R. Patton, and P. Werner, Phys. Rev.B (2012).[32] L. de’ Medici, J. Mravlje, and A. Georges, Phys. Rev.Lett. , 256401 (2011).[33] A similar orbital occupied/unoccupied state also appearsin the J > , 155102(2015).[35] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, andP. Werner, Rev. Mod. Phys. , 779 (2014).[36] M. Eckstein, M. Kollar, K. Byczuk, and D. Vollhardt,Phys. Rev. B , 235119 (2005).[37] For a review see, Y. Kuramoto, H. Kusunose and A. Kiss,J. Phys. Soc. Jpn.78