Development of a two-particle self-consistent method for multi-orbital systems and its application to unconventional superconductors Hideyuki Miyahara1, Ryotaro Arita1,2, and Hiroaki Ikeda3 1Department of Applied Physics, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan 2JST PRESTO, Kawaguchi, Saitama, 332-0012, Japan and 3Department of Physics, Kyoto University, Kyoto 606-8502, Japan (Dated: January 16, 2013) 3 1 We extend the two-particle self-consistent method proposed by Vilk and Tremblay (J. Phys. I 0 France7, 1309-1368 (1997)) to studysuperconductivityin multi-orbital systems. Startingwith the 2 sumrulesforthespinandchargesusceptibilities,wederiveself-consistentequationstodeterminethe n renormalized effectiveinteractions. We apply this method to the two-orbital dx2−y2-d3z2−r2 model for La CuO and the five-orbital d-model for LaFeAsO. Comparing the results with those of the a 2 4 J randomphaseapproximationorthefluctuationexchangeapproximationinwhichvertexcorrections areignored,wediscusshowthevertexcorrectionsaffectthepairinginstabilityofLa CuO andthe 5 2 4 dominant pairing symmetry of LaFeAsO. 1 ] n I. INTRODUCTION In this method, vertex corrections in the charge and o spin channel are assumed to be momentum and fre- c SincetheseminalstudiesbySuhl1 andKondo,2 super- quency independent, and they are determined in such a - r conductivityinmulti-orbitalsystemshasbeenoneofthe way that the correlation functions meet their sum rules. p major topics in condensed matter physics. So far, many Withthisnumericallyinexpensivetreatment,ithasbeen u s kinds of multi-orbital superconductors such as MgB2 demonstrated for the single-band Hubbard model that . (Ref. 3), alkali-doped C (Ref. 4), Na CoO · yH O TPSCshowsgoodagreementwithquantumMonteCarlo t 60 x 2 2 a (Ref.5),Sr RuO (Ref.6),iron-basedsuperconductors,7 (QMC) calculations. 2 4 m and heavy fermion superconductors8 have been discov- Inthispaper,weformulateTPSCforthemulti-orbital - ered and studied extensively. Theoretically, a variety Hubbard model. First, we derive a series of equa- d of exotic unconventional pairing mechanisms going be- tions to determine the vertex corrections in the spin n yondtheMigdal-Eliashbergtheory9 havebeenproposed. and charge channel, and then apply this method to a o c Forexample,ithas beenconsideredfor the cobaltatesu- two-orbital model for La2−x(Sr/Ba)xCuO4 and a five- [ perconductor that the Hund’s coupling (which of course orbital model for F-doped LaFeAsO. Recently, the two- 1 does not exist for single-orbital systems) induces triplet orbital model (which we call the dx2−y2-d3z2−r2 model) superconductivity,10 and it has recently become an is- was studied by FLEX18 to give an insight into the ma- v 4 sue of hot debates whether orbital fluctuations mediate terial dependence of superconducting transition temper- 5 superconductivity in the iron-based superconductors.11 ature (Tc). While FLEX successfully describes the dif- 2 To investigate these fascinating possibilities, accurate ference between La2−x(Sr/Ba)xCuO4 (Tc ∼ 40 K) and 3 calculations of superconductivity in correlated multi- HgBa2CuO4+δ (Tc ∼ 90 K), it underestimates the pair- 1. orbitalmodels areindispensable. Amongmany available ing instability for La2−x(Sr/Ba)xCuO4 and Tc is much 0 approaches,fromtheweakcouplingside,oneoftenstarts lowerthanthe experimentalvalue. We show that, in the 3 with the random phase approximation (RPA). Since the present multi-orbitalTPSC calculation, the inter-orbital 1 pioneering work for the single-band Hubbard model by scattering enhances the d-wave instability and reason- v: Scalapino et al.,12 RPA has been successfully applied to ablevalueofTc isobtainedfortheintermediatecoupling i various multi-orbital systems. The fluctuation exchange regime. Forthefive-orbitalmodel,ithasbeenextensively X approximation (FLEX) developed by Bickers et al.,13 studied by RPA19,20,22 and FLEX.21 There, strong spin ar whichincludestheself-energycorrectionself-consistently, fluctuation has been shown to mediate the s-wave su- has beenalsowidely used. Here, the self-energyis calcu- perconductivitywithsignchanges(theso-calleds±-wave latedinthemannerofBaymandKadanoff,14andconser- pairing). Ontheotherhand,recently,ithasbeenpointed vation laws for one-particle quantities such as the total out that vertex corrections can enhance orbital fluctua- energy and momentum are satisfied. However, due to tions, which mediate s-wave superconductivity without the absence of vertex corrections,FLEX violates conser- sign changes(the s++-wave pairing).16 In this paper, we vation laws for two-particle quantities. show that orbital fluctuations are enhanced in TPSC, Recently, several diagrammatic methods, which take while the dominant pairing symmetry is still s± when into account some vertex corrections, have been the system resides in the weak coupling regime. proposed.15,16 Among them, the two-particle self- This paper is organized as follows. In Sec. II, we for- consistent method (TPSC) proposed by Vilk and mulate multi-orbital TPSC for the Hubbard model. We Tremblay17 is a promising approach in that it is com- discuss how we calculate the charge (orbital) and spin patible with conservation laws in the two-particle level. correlationfunctions. InSec. III,weshowtheresultsfor 2 the two- and five-orbital Hubbard model and the sum- Note that the RPA violates the Pauli principles, and mary of the present study is given in Sec. IV. it does not fulfill the following two sum rules, T χsp(q)=h(n (r)−n (r))(n (r)−n (r))i ↑ ↓ ↑ ↓ N II. METHOD Xq =n−2hn n i, (2a) ↑ ↓ A. Model T χch(q)=h(n (r)+n (r))(n (r)+n (r))i−n2 ↑ ↓ ↑ ↓ N Xq The Hamiltonian of the multi-orbital Hubbard model =n+2hn n i−n2, (2b) is given by ↑ ↓ which are exact relations derived via the Pauli princi- H = ǫµnµσ(r)+ tµrrν′c†µσ(r)cνσ(r′) ples, hnσ(r)2i = hnσ(r)i (see Appendix A). Here, n is Xrµσ rrX′µνσ the particle number per site and for nonmagnetic states, + U nµ↑(r)nµ↓(r)+U′ nµσ(r)nνσ′(r) hn↑(r)i = hn↓(r)i = n/2. Note that the double occu- Xr h Xµ µX>νXσσ′ pancy,hn↑(r)n↓(r)i≡hn↑n↓iisalsotranslationinvariant and does not depend on site r. −J Sµ(r)·Sν(r)+J′ c†µ↑(r)c†µ↓(r)cν↓(r)cν↑(r) , In TPSC, to meet the above conditions [Eqs. (2)], we µX=6 ν µX6=ν i introduce two independent effective interactions, Usp for the spin channel and Uch for the charge channel. Then where c† (r) is a creation operator of an electron the full susceptibilities of Eq. (1) are replaced with µσ with spin σ and orbital µ at site r, and n (r) = µσ c†µσ(r)cµσ(r), Sµ(r) = (c†µ↑(r),c†µ↓(r))σ(cµ↑(r),cµ↓(r))T χsp(q)= 2χ0(q) , χch(q)= 2χ0(q) . (3) with the Pauli matrices σ. The on-site Coulomb inter- 1−Uspχ0(q) 1+Uchχ0(q) actions, U,U′,J,J′ denote the intra-orbital,inter-orbital Finally, we put the following ansatz: Coulomb repulsions, the Hund’s exchange, and the pair- hopping term, respectively. hn n i Usp = ↑ ↓ U, (4) hn ihn i ↑ ↓ which is compatible with the equations of motion (see B. Two-particle self-consistent method for the single-orbital Hubbard model Appendix B). Equations (2), (3), and (4) provide a set of self-consistent equations in TPSC. Namely, Uch, Usp, and hn n i are self-consistently determined for given n Let us start with a review of TPSC for the ↑ ↓ andU. Theone-particleGreen’sfunctionandself-energy single-orbital Hubbard model formulated by Vilk and are calculated by Tremblay.17 The central quantities in this method are the spin and charge correlation functions. In the non- G(k)=G0(k)+G0(k)Σ(k)G(k), magnetic state, the system holds SU(2) symmetry, and 1 T the spin-spin correlationfunctions do not depend on the Σ(k)= Uspχsp(q)U +Uchχch(q)U G(k−q). 4N spin directions. Thus we consider the z component of Xq h i the spin operator Sz(r) = n (r)−n (r) and the charge ↑ ↓ operator n(r) = n (r)+n (r). In RPA, the spin and For the single-band Hubbard model, it has been demon- ↑ ↓ charge correlationfunctions are evaluated as follows: strated that TPSC agrees well with QMC.17 2χ0(q) 2χ0(q) χsp (q)= , χch (q)= , (1) C. Extension to multi-orbital systems RPA 1−Uχ0(q) RPA 1+Uχ0(q) LetushereformulateTPSCforthemulti-orbitalHub- with the irreducible susceptibility, bard model. Hereafter, we follow the matrix form em- T ployed in Refs. 10 and 23. The irreducible susceptibility χ0(q)=− G0(k)G0(k+q), is defined as N Xk T χ0 (q)=− G0 (k)G0 (k+q), (5) where T and N are temperature and number of sites in λµνξ N νλ µξ the system, and G0(k) = 1/(iǫ +µ−ǫ(k)) is the bare Xk n Green’s function with chemical potential µ and energy which can be considered as a matrix element with a row dispersion ǫ(k). Here, we have introduced the abbrevi- λµ and a column νξ of a matrix χ0(q). In the non- ations k = (k,iǫ ) and q = (q,iν ) being the fermionic magnetic state, the system is invariant for spin rota- n n and bosonic Matsubara frequencies, respectively. tion, and then 2χspz(q)=χsp±(q) holds, where χsp±(q) 3 is the in-plane correlation function between S± (r) = ceptibilities and the double occupancy: λµ Sx (r) ± iSy (r) = c† (r)c (r) with σ =↑ or ↓. In λµ λµ λσ µσ¯ T TPSC, similar to the single-orbital case, the spin and χch (q)=h(n +n )(n +n )i−hn ihn i, N µµµµ µ↑ µ↓ µ↑ µ↓ µ µ charge susceptibilities are given by Xq =n +2hn n i−n2 (9a) µ µ↑ µ↓ µ χsp(q)=(1−χ0(q)Usp)−12χ0(q), (6a) T χch (q)=h(c† c +c† c )(c† c +c† c )i χch(q)=(1+χ0(q)Uch)−12χ0(q), (6b) N Xq µνµν µ↑ ν↑ µ↓ ν↓ ν↑ µ↑ ν↓ µ↓ =hn (1−n )i+hn (1−n )i µ↑ ν↑ µ↓ ν↓ where Usp(ch) is the renormalized effective interaction +hc† c c† c i+hc† c c† c i µ↓ ν↓ ν↑ µ↑ µ↑ ν↑ ν↓ µ↓ matrixforthespin(charge)channel.10,23 Fortwo-orbital T systems, for instance, these are represented as = χspz (q)+2χspz (q) . (9b) N(cid:18) µνµν µµνν (cid:19) Xq Usp Jsp 0 0 Finally, as in the single-band case, we introduce the 1111 Usp = Jsp U2s2p22 0 0 , (7a) following ansatz between the two-particle quantities and 0 0 Usp Jsp the interaction parameters (see Appendix C); 1212  0 0 Jsp Usp   2121 hn n i Uch 2Uch −Jch 0 0 Usp = σµ σ¯µ U, (10a) 2Uch11−11Jch 11U22ch 0 0  µµµµ hnσµihnσ¯µi Uch = 2211 2222 ,  0 0 −U1c2h12+2Jch 0  Usp = hnσµnσ¯νi U′, (10b)  0 0 0 −Uch +2Jch µνµν hn ihn i  2121  σµ σ¯ν (7b) hn n i Usp −Jsp= σµ σν (U′−J). (10c) µνµν hn ihn i σµ σν where Usp (Uch ) is the intra-orbital Coulomb in- µµµµ µµµµ Equations. (6)-(10) are a set of self-consistent equations teraction; Usp (Uch ) with µ 6= ν, the inter-orbital µνµν µνµν in the multi-orbital case. Coulombinteraction;Jsp(Jch),theHund’scoupling.10,23 Inthe presentstudy,forsimplicity,weignoretheHund’s coupling in the charge channel, namely, Jch = 0.24 In D. Eliashberg equation RPA, one employs the unperturbed bare vertex as fol- lows, Uµspµ(µcµh) = U, Uµspµ(νcνh) =Uµspν(µcνh) =U′, and Jsp(ch) = Superconductivity has been studied by the following J. linearized Eliashberg equation, Nextletus considerthe sumruleformulti-orbitalsys- tems(seeAppendixA). Inthen-orbitalHubbardmodel, λ∆ll′(k)= Vlm1m4l′(k,k′)Gm1m2(k′) there are n4 sum rules for χsp(q) and χch(q). Among kX′,mi (11) them, we use the following equations to determine Usp: ×∆ (k′)G (−k′). m2m3 m4m3 T Eigenstate ∆ll′(k) with the largest eigenvalue λ was nu- χspz (q)=2hn i−2hn n i, (8a) mericallyevaluatedbythepowermethod. Thesupercon- N µµµµ µ↑ µ↑ µ↓ Xq ducting transitionoccursatthe temperatureforwhichλ T becomesunity. Here,Gll′(k)isthe dressedGreen’sfunc- N χsµpν±µν(q)=2hc†µ↑cν↓c†ν↓cµ↑i tion, Xq =2hnµ↑i−2hnµ↑nν↓i, (8b) Gll′(k)=G0ll′(k)+G0lm(k)Σmm′(k)Gm′l′(k), (12) T N χzµ±µνν(q)=2hnµ↑nν↑i−2hnµ↑nµ↓i. (8c) and the self-energy Σll′(k) is given by Xq 1 T Σll′(k)= 4N Uspχsp(q)Us0p Note that the intra-orbital component of the sum rule Xq h (13) hmaosdtehl.eFsoarmtehefoirnmtera-sorthbaittalocfotmheposinnegnltes-,orwbeituasleHχusbpb±a(rqd) +Uchχch(q)Uc0hilml′m′Gmm′(k−q). ratherthanχspz(q),sincetheycanbeexpressedinterms In the present study, we omit the Hartree-Fock term, of the density operators. since a part of its contribution is already considered in For Uch, we use the followingsum rules for the charge the one-body part of the Hamiltonian, which is derived susceptibilitieswhichcanbe representedbythespinsus- from density functional calculation. 4 The effective interaction Vll′mm′(k,k′) for the spin- La2CuO4 singlet pairing can be expressed in a matrix form as fol- 2 lows: ] V V(k)=−3Uspχsp(k) Usp+ 1Uchχch(k)Uch [e 1 2 0 2 0 y (14) g − 1Usp− 1Uch, ner 0 2 0 2 0 E where Usp and Uch are the bare vertex in the spin and -1 0 0 charge channel, respectively.10,23 -2 Γ N X Γ III. RESULTS FIG. 1: (Color online) Band structure of the two-orbital Let us move on to the application of the multi-orbital model for La2CuO4. The model consists of the dx2−y2 or- TPSC method to the effective models for La2CuO4 and bital and the d3z2−r2 orbital. The Fermi level is set at 0 eV. LaFeAsO. Using the technique of the maximally local- ized Wannier functions,25 these models are derived from 1 first-principles calculations. In the density-functional TPSC calculations,weemployedthe exchangecorrelationfunc- eq. FLEX tional proposed by Perdew et al.,26 and the augmented erg 0.8 b planewaveandlocalorbital(APW+lo)methodasimple- ash 0.6 mentedintheWIEN2Kprogram.27Wethenconstructed Eli the Wannier functions for the d bands around the Fermi e of 0.4 level, using the WIEN2Wannier (Ref. 28) and the wan- alu nier90 (Ref. 29) codes. nv 0.2 e g Ei 0 A. La CuO 0 0.02 0.04 0.06 2 4 Temperature [eV] Recently, the two-orbital dx2−y2-d3z2−r2 Hubbard model for the cuprates were studied to understand the FIG. 2: (Color online) Temperature dependence of the max- materialdependenceofTc byFLEX.18 There,theenergy imum eigenvalue of the linearized Eliashberg equation ob- difference betweenthe dx2−y2 orbitalandthe d3z2−r2 or- tainedbyTPSC(redsolidline)andFLEX(bluedottedline). bital was found to be a key parameter to characterize U and J are 2.0 eV and 0.2 eV, respectively, and n is set to La CuO andHgBa CuO . Namely,inthe FLEXcalcu- be2.85. 2 4 2 4 lationforthetwo-orbitalmodel,thepairinginstabilityis stronger in the latter. On the other hand, in the former, the eigenvalue of the Eliashberg equation within FLEX orbital and strong commensurate inter-orbital fluctua- does not reach unity down to T ∼ 40 K. The purpose of tions. We here stress that there is no large peak in the this subsection is to examine how the vertex corrections charge susceptibilities in the RPA and FLEX calcula- in TPSC affect the superconductivity in La CuO . tions, so that these enhanced inter-orbital charge fluc- 2 4 The band structure of the effective two-orbital model tuations are purely due to the effects of vertex correc- for La CuO is shown in Fig. 1. We set U = 2.0 eV, tions. It should be noted that χch (Q) has a negative 2 4 1221 U′ = 1.6 eV, J = 0.2 eV, and n = 2.85. We employ peak around Q=(π,π), which works as attractive force 64×64k-pointmeshesand2048Matsubarafrequencies. betweenthetwoorbitalsford-wavepairingjustlikeanti- Hereafter,orbitals1and2denotethedx2−y2 andd3z2−r2 ferromagnetic spin fluctuations. [Note that the spin and orbitals, respectively.30 charge sectors in Eq. (14) have opposite signs.] As we In Fig. 2, we plot temperature dependence of λ, the will see below, there is a close correlation between the maximum eigenvalue of the Eliashberg equation. While characteristic enhancement of λ and the charge fluctua- λ does not show appreciable temperature dependence tions. in FLEX, λ is drastically enhanced at low temperature Figure 4 depicts temperature dependence of the max- <0.02eVinTPSC.Thecharacteristicenhancementofλ imum value in the spin susceptibility, the orbital suscep- in TPSC is attributed to the low-temperature behaviors tibility, and its inverse. We can see that the enhance- of the spin and charge susceptibility. In Fig. 3, we plot ment of the peaks in χch and −χch dominate over 1212 1221 χsp (q,ω = 0), χch (q,ω = 0), and −χch (q,ω = 0) that of χsp for T < 0.02 eV. This behavior comes 1111 1212 1221 1111 at T = 0.020 eV, which indicate that the system has from the fact that while Usp is renormalized substan- a strong incommensurate spin correlation in the dx2−y2 tially (U1s1p11 ∼ 1.16 eV at T = 0.020 eV), Uc3h3 (Uc4h4) 5 (a)χsp (b)χch (c) -χch 1111 1212 1221 9 12 12 6 6 6 3 0 (π, π) 0 (π, π) 0 (π, π) q q q (0, 0) q y (0, 0) q y (0, 0) q y x (π, 0) x (π, 0) x (π, 0) FIG.3: (Coloronline)(a)χsp (q,ω=0),(b)χch (q,ω=0),and(c)−χch (q,ω=0)atT =0.020eV,whereorbitals1and 1111 1212 1221 2 denote thedx2−y2 and d3z2−r2 orbitals, respectively. (a) (b) (c) 20 10 1.6 π-δ) 68 |π,π) 1126 -χχc1c1h2h21221 |-1π,π) 1.2 -χχc1c1h2h21221 spπ,χ(1111 24 χ(1212(1221) 48 chχ(1212(1221) 00..48 0 | 0 | 0 0 0.02 0.04 0.06 0 0.02 0.04 0.06 0 0.02 0.04 0.06 Temperature [eV] Temperature [eV] Temperature [eV] FIG. 4: (Color online) Temperature dependenceof themaximum value in (a) χsp , (b) |χch |, and (c) |χch −1|. 1111 1212(1221) 1212(1221) becomes ∼ 1.55 eV, which is even larger than the bare ture T = 0.018 eV. We see that characteristic structure value 1.2eVusedinRPAandFLEX.Infact, similaren- emerges at (π,0) and (0,π) in the d3z2−r2 gap function, hancement of charge channelis also observedin the case due to the inter-orbital effective interaction, dominantly of the single-orbitalmodel.17 mediated by the orbital fluctuation χch (q), which con- 1221 Figure5showsthegapfunctionsforthedx2−y2 orbital nects ∆11 and ∆22 [Eq. (11)]. Since χc1h221(q) takes a and the d3z2−r2 orbital (∆11 and ∆22) at T =0.022 eV. largenegativevaluearound(π,π),itcooperateswithan- They have the d-wave symmetry, which is mediated by tiferromagnetic spin fluctuation to enhance the d-wave the dominant spin fluctuation, denoted by the black ar- pairing instability. row. This is the conventionalsituationwhere the orbital fluctuations remains small. (a) dx2-y2 (b) dz2 (a) dx2-y2 (b) dz2 π π π π q 0 q 0 y y q 0 q 0 y y -π -π -π 0 π -π 0 π -π -π q q -π 0 π -π 0 π x x q q x x -0.003 0 0.003 -0.0015 0 0.0015 -0.006 0 0.006 -0.0008 0 0.0008 FIG.6: (Coloronline)PlotssimilartoFig.5forT =0.018eV. Thedashedblackarrowdenotesthepairscatteringmediated FIG. 5: (Color online) Gap function for the (a) dx2−y2 or- byorbital fluctuations. bital and the (b) d3z2−r2 orbital at T = 0.022 eV. The gap functionshaved-wavesymmetryandtheblackarrowdenotes the pair scattering mediated by antiferromagnetic spin fluc- This situationchangesinthe strongercouplingregime tuations. (U ∼ 2.5 eV), where the system goes away from a su- perconducting instability. This is because the dominant In Fig. 6, we plot the gap functions at lower tempera- spin/orbital fluctuations make the quasi-particle damp- 6 (a) (b) (c) 100 50 q. 1 π-δ) 6800 |π,π() 3400 -χχc1c1h2h21221 Eliashberg e 00..68 spχπ,(1111 2400 chχ1212(1221) 1200 nvalue of 00..24 0 | 0 ge 0 0 0.02 0.04 0.06 0 0.02 0.04 0.06 Ei 0 0.02 0.04 0.06 Temperature [eV] Temperature [eV] Temperature [eV] FIG. 7: (Color online) Temperature dependence of the maximum value in (a) χsp , (b) |χch |, and (c) the maximum 1111 1212(1221) eigenvalue of the linearized Eliashberg equation obtained by TPSC for U =2.5 eV, U′ =2.0 eV,and J =0.25 eV. ing around (0,π) and (π,0) significant, as is observed in LaFeAsO the previous FLEX calculation.18 Temperature depen- 2 dence of the maximum value in the spin, orbital suscep- tibility and the maximum eigenvalue of the Eliashberg ] V 1 equationλfor U =2.5eV, U′ =2.0 eVand J =0.25eV e [ areshowninFigs. 7(a),7(b),and7(c),respectively. The y g maximum value in χsp always dominates over that of er 0 1111 n χch and−χch ,andλdoesnotshowanyenhancement. E 1212 1212 -1 -2 B. Iron-based superconductor: LaFeAsO (0, 0) (π, 0) (π, π) (0, 0) Letusnowapplymulti-orbitalTPSCtotheiron-based superconductor, LaFeAsO. The recent discovery of high FIG. 8: (Color online) Band structure of the five-orbital d- Tc superconductivity in F-doped LaFeAsO32 has stimu- model for LaFeAsO.The Fermi level is set at 0 eV. latedarenewedinterestinmulti-orbitalsuperconductors. Asforthepairingmechanismoftheiron-basedsupercon- ductors, several scenarios have been proposed. Among than that of χsp (q,ω = 0) for T > 0.01 eV. While 2222 them, the possibility of the sign-reversing s -wave su- χch (q,ω = 0) has a broad maximum peak around ± 2424 perconductivity mediated by spin fluctuations22,33 have T = 0.01 eV, χsp (q,ω = 0) grows monotonously as 2222 been extensively studied. While the s -wave solution temperature lowers. ± has been obtained in the RPA19,20,22 or FLEX21 calcu- The enhancement in the orbital susceptibility comes lations for the five-orbital d-model, recently, it has been from the vertex correction in the charge susceptibility. proposedthatvertexcorrectionscanenhanceorbitalfluc- Tomakethispointclear,inFig. 11,weplottemperature tuations, and the s -pairing without sign reversingbe- dependenceofthemaximumvalueinχch (q,ω =0)and ++ 2424 comes dominant.11 In this subsection, we discuss how χsp (q,ω = 0) obtained by RPA, for which the bare 2222 vertex corrections in TPSC affects superconductivity in coupling constants are set to be U = 1.2 eV, U′ = 0.96 the five-orbital d-model for LaFeAsO. eV, and J = 0.12 eV. We see that while χsp (q,ω = 0) 2222 The band structure of the d-model is shown in Fig. 8. diverges around T = 0.01 eV, χch (q,ω = 0) has no 2424 The bare coupling constants are set to be U = 1.5 eV, significant temperature dependence. U′ =1.2eV,andJ =0.15eV.34 Thenumberofelectrons In Fig. 12, we show temperature dependence of the n is 6.1. We employ 64 × 64 k-point meshes and 2048 maximum eigenvalue of the Eliashberg equation. We see Matsubara frequencies. Hereafter orbitals 1, 2, 3, 4, and thatthesystemhasasuperconductingtransitionaround 5 denote the d3z2−r2, dxz, dyz, dx2−y2, and dxy orbitals, T ∼0.005 eV. respectively. The associated eigenfunctions of the Eliashberg equa- To see that TPSC can give enhanced orbital fluctua- tion at T = 0.015 eV (the gap functions) are shown tions,weplotχsp (q,ω =0)andχch (q,ω =0)atT = in Fig. 13 for the three bands crossing the Fermi level. 2222 2424 0.015 eV in Figs. 9(a) and 9(b), respectively. Clearly, We see that these gap functions have the s symme- ± both susceptibilities havepeaks around(π,0) and (0,π), try, indicating that the spin fluctuation is the primary and the peak in the orbital susceptibility is higher than glue of superconductivity. However, there is a notable that of the spin susceptibility. Temperature dependence difference between TPSC and RPA results in the am- of these peaks are shown in Fig. 10. We see that the plitudes of the gap functions on the Fermi surface. As peak of χch (q,ω = 0) is more drastically enhanced we can see in Fig. 14, the gap amplitude is larger for 2424 7 (a) χsp (b) χch 2222 2424 1.6 4.5 3 0.8 1.5 (π, π) (π, π) 0 0 q q (0, 0) y (0, 0) y q q x (π, 0) x (π, 0) FIG.9: (Color online) (a) χs2p222(q,ω=0) and (b)χc2h424(q,ω=0) at T =0.015 eV,whereorbitals 2and 4denotethedxz and dx2−y2 orbitals, respectively. ) 1 δ, 0 5 χs2p222 Eigenvalue of χspπ-(2222 34 χc2h424 hberg Eq. 00..68 Eliashberg Eq. d as an Eli π-δ, 0) 12 alue of 0.4 ( v 0.2 χch2424 0 0 0.01 0.02 0.03 0.04 0.05 Eigen 0 Temperature [eV] 0 0.01 0.02 0.03 0.04 0.05 Temperature [eV] FIG.10: (Coloronline)Temperaturedependenceofthemax- imum value in χsp (red solid line) and χch (blue dotted FIG.12: (Coloronline)Temperaturedependenceofthemax- 2222 2424 line) obtained byTPSC. imum eigenvalue of the linearized Eliashberg equation. The bluedotted line is a guide tothe eye. ) 0 10 χsp-RPAπ-δ, (2222 68 χχs2c2p2h4-22-RR24PPAA Vtthh(eeknF)ew=rem−fiins32udUrftsahpcaχetsapts(hkien)gUtahps0epRaimnPAptlhirteeusdEuellitab[seshecbeoemFrgeigse.lq1au4rag(dtei)o]on.n, d n a 4 ) 0 δ, IV. SUMMARY - 2 π (A χch-RP2424 0 0 0.01 0.02 0.03 0.04 0.05 conTsoisstuemntmmaertihzeo,dw(eThPaSvCe)dfeovretlohpeemdutlhtei-otwrboi-tpaalrHtiucblebaserldf Temperature [eV] model. Wederivedself-consistentequationstodetermine vertexcorrectionsinthespinandcharge(orbital)suscep- FIG.11: (Coloronline)Temperaturedependenceofthemax- tibilities. We applied this method to the effective mod- imum value in χs2p222 (red solid line) and χc2h424 (blue dotted els for La2CuO4 and LaFeAsO. We solved the linearized line) obtained by RPA. The bare coupling constants are set Eliashberg equation and found that vertex corrections to beU =1.2 eV, U′ =0.96 eV,and J =0.12 eV. play a crucial role in the multi-orbital superconductors. In the two-orbital dx2−y2-d3z2−r2 model for La2CuO4, while FLEX shows much lower T than its experimental c RPA than TPSC. This indicates that there is a frus- value ∼ 40 K, the present TPSC can increase Tc dra- tration between the orbital-fluctuation-mediated pair- maticallyduetoenhancedorbitalfluctuationsviavertex ing and the spin-fluctuation-mediated pairing. Indeed, corrections for intermediate U ∼ 2.0 eV. in TPSC, if we drop the contribution of the charge In the iron-based superconductor LaFeAsO, we have channel in the pairing interaction, namely consider only studiedwhetherorbitalfluctuationscanbeenhancedand 8 (a) Band 2 (b) Band 3 (c) Band 4 π π π q q q 0 0 0 y y y -π -π -π -π 0 π -π 0 π -π 0 π q q q x x x 0 0.001 0.002 -0.001 0 0.001 0.002 -0.0012 -0.0008 -0.0004 0 FIG. 13: (Color online) Gap functions obtained by TPSC for the bands with the (a) second, (b) third, and (c) fourth Kohn- Sham energy at T = 0.015 eV. The black line and dotted green line represent the Fermi surface and nodes of gap functions, respectively. (b) (c) (d) (a) ] π -30 6 FS1 6 6 FS4 1 FS2 [ FS3 FS1 FS3 s ∆ 4 FS4 4 4 q 0 n y FS2 ctio 2 2 2 n u -π p f 0 0 0 -π 0 π Ga q -2 -2 -2 x 0 π/2 0 π/2 0 π/2 Angle θ FIG.14: (Coloronline)Gapfunctiononthe(a)Fermisurfacesby(b)TPSC,(c)RPA,and(d)TPSCwithoutchargefluctuation, where θ is therotation angle from theky axis. induce the s -wave pairing within TPSC. Indeed we corrections in multi-orbital systems for cooperative and ++ have found that some kinds of orbital fluctuations are competitivephenomenabetweenspinandorbitaldegrees enhanced by considering vertex correction and become of freedoms. even stronger than spin fluctuations. However, their or- bitalfluctuationsarenotstrongenoughtocausethes - ++ wave pairing, and the obtained gap function has the s ± symmetry,althoughthe gapmagnitude is relativelysup- pressedduetoafrustrationbetweentwokindsofpairing interactions mediated by spin and orbital fluctuations. Acknowledgments It is an interesting problem in future research whether the pairing symmetry changes for larger interaction pa- rameters. Anotherimportantfutureissueisasystematic WethankH.KontaniandS.Onariforstimulatingdis- comparisonbetweenthepresentmulti-orbitalTPSCand cussions. This work was supported by Grants-in-Aid for other(diagrammatic)methodswhichconsiderthevertex ScientificResearch(No.23340095)fromMEXTandJST- corrections. Finally, we stress the importance of vertex PRESTO, Japan. 9 Appendix A: Definition of correlation functions (δ >0). Inthe multi-orbitalcase,we considercorrelationfunc- Inthesingle-orbitalcase,correlationfunctionsforspin tions for S (1) = (c† (1),c† (1))σ(c (1),c (1))T, µν µ↑ ν↓ µ↑ ν↓ Sz(r) and charge n(r) are defined as and n (1) = c† (1)c (1)+c† (1)c (1). These corre- µν µ↑ ν↑ µ↓ ν↓ lation functions for spin and charge channels are defined χsp(1,2)=hT Sz(1)Sz(2)i, (A1a) τ as χch(1,2)=hT n(1)n(2)i−hn(1)ihn(2)i, (A1b) τ χspz (1,2)=hT Sz (1)Sz (2)i, (A2a) where an abbreviation 1=(r ,τ ) denotes a position r λµνξ τ λµ ξν 1 1 1 and an imaginary time τ1, and Tτ is the time ordering χsλpµ±νξ(1,2)=hTτS+λµ(1)S−ξν(2)i, (A2b) operator. Time dependence of a generic operator Q(r) χch (1,2)=hT n (1)n (2)i−hn (1)ihn (2)i. is defined as Q(r,τ) = e−τHQ(r)eτH. χsp(ch)(q) in the λµνξ τ λµ ξν λµ ξν (A2c) maintextistheFouriertransformoftheabovereal-space representation. Thetwosumrules(Eqs.(2)),whichplay a central role in TPSC, originate from the definition at The sum rules of Eqs. (8) and (9) come from the fol- equal time, that is, χsp(ch)(1,1+) with 1+ = (r ,τ +δ) lowing definitions at equal time, 1 1 χspz (1,1+)=h(c† (1)c (1)−c† (1)c (1))(c† (1+)c (1+)−c† (1+)c (1+))i λµνξ λ↑ µ↑ λ↓ µ↓ ξ↑ ν↑ ξ↓ ν↓ (A3a) −h(c† (1)c (1)−c† (1)c (1))ih(c† (1+)c (1+)−c† (1+)c (1+))i λ↑ µ↑ λ↓ µ↓ ξ↑ ν↑ ξ↓ ν↓ χsp± (1,1+)=2hc† (1)c (1)c† (1+)c (1+)i (A3b) λµνξ λ↑ µ↓ ξ↓ ν↑ χch (1,1+)=h(c† (1)c (1)+c† (1)c )(c† (1+)c (1+)+c† (1+)c (1+))i λµνξ λ↑ µ↑ λ↓ µ↓ ξ↑ ν↑ ξ↓ ν↓ (A3c) −h(c† (1)c (1)+c† (1)c )(c† (1+)c (1+)+c† (1+)c (1+))i. λ↑ µ↑ λ↓ µ↓ ξ↑ ν↑ ξ↓ ν↓ Appendix B: Ansatz for effective interactions in with σ¯ = −σ. Here a bar over a number means the in- single-orbital case tegral over position and imaginary time. The four-point correlation function can be approximated by the local Following Vilk and Tremblay,17 let us derive the correlationfunction and the Green’s function as follows, ansatz, Eq. (4), used in TPSC calculations. The four- point vertex function, Γσσ′, between electrons with spin σ and σ′ is given by Γσσ′δ(1−3)δ(2−4)δ(2−1+)= δΣσ(1,2), (B1) −UhTτ[cσ†¯(1++)cσ¯(1+)cσ(1)c†σ(2)]i δGσ′(3,4) hn (1)n (1)i (B3) ∼U ↑ ↓ G (1,1+)G (1,2). σ¯ σ hn (1)ihn (1)i where G (1,2) and Σ (1,2) is the dressed Green’s func- ↑ ↓ σ σ tionandselfenergywithspinσ. Theequationofmotion and the Dyson equation leads to the relation, Σ (1,¯1)G (¯1,2) σ σ (B2) By substituting Eq. (B3) into Eq. (B1), we can obtain =−UhT [c†(1++)c (1+)c (1)c†(2)]i, τ σ¯ σ¯ σ σ δ U hn↑n↓i G (1,1+)δ(1−2) Γσσ′δ(1−3)δ(2−4)δ(2−1+)= δΣσ(1,2) = h hn↑ihn↓i σ¯ i δGσ′(3,4) δGσ′(3,4) (B4) = δhUhnhn↑↑inhn↓↓iiiG (1,1+)δ(1−2)+U hn↑n↓i δGσ¯(1,1+)δ(1−2). σ¯ δGσ′(3,4) hn↑ihn↓i δGσ′(3,4) The lastterm ofthis equationis proportionalto δσ¯σ′ via and then contributes to the spin channel, Usp = Γσσ¯ − δG (1,1+) Γσσ. This just provides Eq. (4) for Usp. σ¯ =δσ¯σ′δ(1−3)δ(4−1+), (B5) δGσ′(3,4) 10 Appendix C: Ansatz in multi-orbital systems pointvertexfunctionhasorbitalindices,λ,µ,ν,ξbesides spin index, σ. We here consider Γ , which can (µµσ)(ννσ′) In this section, let us extend the above-mentioned be written by only orbital-diagonalcomponents. ansatz into the multi-orbital case. In this case, the four- Γ δ(1−3)δ(2−4)δ(2−1+)= δΣµµσ(1,2) = δ Σµµσ(1,¯5)[G(¯5,¯6)G−1(¯6,2)]µµσ (µµσ)(ννσ′) δGννσ′(3,4) (cid:2) δGννσ′(3,4) (cid:3) δ hn n i hn n i hn n i ∼ −U µσ µσ¯ Gµµσ¯(1,1+)− U′ µσ ξσ¯ Gξξσ¯(1,1+)− (U′−J) µσ ξσ Gξξσ(1,1+) δGννσ′(3,4) hnµσihnµσ¯i ξX6=µ hnµσihnξσ¯i ξX6=µ hnµσihnξσi   hn n i δG (1,1+) hn n i δG (1,1+) hn n i δG (1,1+) ∼−U µσ µσ¯ µµσ¯ − U′ µσ ξσ¯ ξξσ¯ − (U′−J) µσ ξσ ξξσ . (C1) hnµσihnµσ¯i δGννσ′(3,4) ξX6=µ hnµσihnξσ¯i δGννσ′(3,4) ξX6=µ hnµσihnξσi δGννσ′(3,4) Here, following the single-orbital case, we have intro- Eq. (10a) for Usp can be obtained from the first term µµµµ duced the following approximations, of Eq. (C1), since the intra-orbital Coulomb interaction for the orbital µ is proportional to δ δ . In the same µµ σσ¯ Σ G G−1 ∼−U hnµσnµσ¯i G (1,1+), way, Eq. (10b) for Uµspνµν and Eq. (10c) for Uµspνµν −Jsp µµσ µµσ µµσ hn ihn i µµσ¯ with µ 6= ν can be obtained from the second and the µσ µσ¯ hn n i third terms in Eq. (C1), respectively. Σ G G−1 ∼−U′ µσ ξσ¯ G (1,1+), µµσ ξξσ ξξσ hn ihn i ξξσ¯ µσ ξσ¯ hn n i Σ G G−1 ∼−U′ µσ ξσ¯ G (1,1+). µµσ ξξσ¯ ξξσ¯ hn ihn i ξξσ µσ ξσ¯ 1 H. Suhl, B.T. Matthias, L.R. Walker, Phys. Rev. Lett., 3 14 G. Baym and L. P. Kadanoff: Phys. Rev. 124 (1961) 287, 552 (1959) G. Baym: ibid. 127 (1962) 1391. 2 J. Kondo,Progr. Theor. Phys., 29 1 (1963) 15 H. Kusunose,J. Phys.Soc. Jpn.79, 094707 (2010). 3 J. Nagamatsu, N. Nakagawa, T. Muranaka, Y. Zenitani 16 S. Onari and H. Kontani, Phys. Rev. Lett. 109, 137001 and J. Akimitsu,Nature (London) 410 63 (2001). (2012). 4 Forreview,see,e.g., O.Gunnarsson,Rev.Mod.Phys.69, 17 Y.M. Vilk, A.-M.S. Tremblay, J. Phys. I France 7, 1309 575(1997), M.Capone,M.Fabrizio, C.Castellani, andE. (1997). Tosatti, ibid. 81, 943 (2009). 18 H.Sakakibara, H.Usui, K.Kuroki,R.Arita and H.Aoki, 5 K. Takada et al.,Nature(London) 422, 53 (2003). Phys. Rev. Lett. 105, 057003 (2010), Phys. Rev. B 85 6 For review, see, e.g., A.P. Mackenzie and Y. Maeno, Rev. 064501 (2012). Mod. Phys. 75, 657 (2003). Y. Maeno, S. Kittaka, T. No- 19 S. Graser, T. A. Maier, P. J. Hirschfeld, and D. J. mura, S. Yonezawa, and K. Ishida, J. Phys. Soc. Jpn, 81, Scalapino, New J. Phys. 11, 025016 (2009). 011009 (2009). 20 K.Kuroki,S.Onari,R.Arita,H.UsuiandH.Aoki,Phys. 7 For review, see, e.g., J. Paglione and R. L. Greene, Nat. Rev.B 79,224511 (2009). Phys. 6 645 (2010). 21 H. Ikeda, J. Phys. Soc. Jpn. 77, 123707 (2008), R. Arita 8 For review, see, e.g., C. Pfleiderer, Rev. Mod. Phys., 81 and H. Ikeda, ibid. 78, 113707 (2009), H. Ikeda, R. Arita 1551 (2009). and J. Kuneˇs, Phys. Rev. B 82 024508 (2010); 81 054502 9 J. R. Schrieffer, Theory of Superconductivity (Westview (2010). Press, Colorado, 1971). 22 K. Kuroki, S. Onari, R. Arita, H. Usui, H. Kontani, Y. 10 M. Mochizuki, Y.Yanase and M. Ogata, Phys.Rev.Lett. TanakaandH.Aoki,Phys.Rev.Lett.101,087004(2008). 94, 147005 (2005). 23 T. Takimoto, T. Hotta, and K. Ueda, Phys. Rev. B 69, 11 H. Kontani and S. Onari, Phys. Rev. Lett. 104, 157001 104504 (2004). (2010). 24 TheHund’scouplinginthechargechannelisalwaysrenor- 12 D. J. Scalapino, E. Loh Jr., and J. E. Hirsch, Phys. Rev. malized tobesmall inourcase, andsometimescausesun- B 34, 8190 (1986). expected problematic behavior. 13 N. E. Bickers and D.J. Scalapino: Ann.Phys.(N.Y.) 193 25 N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1989) 206. (1997), I. Souza, N. Marzari and D. Vanderbilt, ibid. 65,