Diagrammatic analysis of the Hubbard model: Stationary property of the thermodynamic potential V. A. Moskalenko1,2,∗ L. A. Dohotaru3, and I. D. Cebotari1 1Institute of Applied Physics, Moldova Academy of Sciences, Chisinau 2028, Moldova 2BLTP, Joint Institute for Nuclear Research, 141980 Dubna, Russia and 3Technical University, Chisinau 2004, Moldova (Dated: January 12, 2010) Diagrammatic approach proposed many years ago for strong correlated Hubbardmodel is devel- opedforanalyzingofthethermodynamicpotentialproperties. Thenewexactrelationbetweensuch 0 renormalized quantities as thermodynamic potential, one-particle propagator and correlation func- 1 tionisestablished. Thisrelationcontainsadditionalintegrationoftheone-particlepropagatorbythe 0 auxiliary constant. The vacuum skeleton diagrams constructed from irreducible Green’s functions 2 andtunnelingpropagator linesaredetermined andspecial functionalisintroduced. Theproperties of such functional are investigated and its relation to the thermodynamic potential is established. n The stationary properties of this functional with respect to first order changing of the correlation a J function is demonstrated and as a consequence the stationary properties of the thermodynamic potential is proved. 1 1 PACSnumbers: 71.27.+a,71.10.Fd ] l r-e I. INTRODUCTION HereC−→+i σ(C−→i σ)arethe−→creation(annihilation)electron t operators with local site i and spin σ. Because in the s TheHubbardmodelisoneofthemostimportantmod- . thermodynamicperturbationtheoryweshallusethermal t els for the electron of solids which describes quantum a averagesinagrandcanonicalensemblewehaveaddedto m mechanical hopping of electron between lattice sites and the Hamiltonian (1) the term −µNˆ e their short ranged repulsive Coulomb interaction. - d This model was discussed by Hubbard [1] in order to Nˆ = n−→ (4) n describe a narrow-band system of transition metals and e −→X i σ o has been revised to investigate the properties of highly i σ [c correlated electron systems such as the copper oxide su- whereµisthechemicalpotentialandNˆeelectronnumber perconductors and others. operator. The quantities U and Nˆe are the fundamental 1 Hubbard model exhibits various phenomena includ- parameters of the model and because of large value of v ingmetal-insulatortransition,antiferromagnetism,ferro- the Coulomb repulsion it is taken into account in zero 7 4 magnetism and superconductivity. This model assumes approximation of our theory. The operator H′, which thateachatomofthecrystallatticehasonlyoneelectron 6 describes hopping of the electrons between sites of the orbitandthe correspondingorbitalstate is nondegener- 1 crystal lattice is regarded as a perturbation. . ate. For investigating this model new physical and mathe- 1 The Hamiltonian of Hubbard model is a sum of the 0 matical concepts and techniques have been elaborated. two terms 0 We can only enumerate some of them. There are 1 H = H0+H′, (1) many analytical approximations as Hubbard approxi- : mation, noncrossing approximation (NCA), slave-boson v i where H0 is the atomic contribution, which contains the method, Dynamical Mean Field Theory (DMFT), com- X Coulomb interaction term U and local electron energy ǫ posite and others projection operator methods and ap- r on the atom proaches. Also numerical simulations of thermodynamic a quantities and density of states have been performed H0 = H−→0 , i by various methods. Some exact results for one- and X−→ i two-dimensional Hubbard model are known. Each ap- H−→0 = ǫn−→ +Un−→ n−→ , (2) proximate method has its advantages and disadvan- i X i σ i ↑ i ↓ tages. A short and comprehensive reviews of the meth- σ ǫ = ǫ−µ, n−→ =C−→+ C−→ ods can be found in papers and books [2−7]. Besides i σ i σ i σ these approacheswe must enumerate also the special di- agramtechniqueselaboratedforstronglycorrelatedelec- and hopping Hamiltonian tron systems. Stasyuk [8], Zaitsev [9], Izyumov [10] and −→ −→ H′ = t( i − j )C−→+ C−→ , (3) coworkers have developed a diagram technique for Hub- −→X−→Xσ i σ j σ bard model, based on the disentangling of nondiagonal i j Hubbard operator χmn out of a time-ordered products −→ −→ −→ −→ t( i − j ) = t∗(j − i ), t(0)=0. of other such operators. Due to a more complicate al- 2 gebra of Hubbard transfer operator than that for Fermi II. PERTURBATION TREATMENT operators the essential features of this technique remain till now poorly developed. We shallusethe definitionofthe one-particleMatsub- The other diagrammatic approach around the atomic ara Green’s functions in interaction representation as in limit also has been proposedfor Hubbard model both in paper [11,12]: normal[11,12] andinsuperconductingstate[13]. Thisthe- ory introduces the Generalized Wick Theorem (GWT) G(x|x′) = − TC−→xσ(τ)C−→x′σ′(τ′)U(β) c0, (5) that uses cumulant expansion of the statistical average (cid:10) (cid:11) −→ c where x stands for (x,σ,τ) and index c for h...i means values for the products of the Fermion operators. The 0 GWT takes into account the fact that Hamiltonian H0 the connected part of the diagrams which appear in the right-hand of part of (5). We use the series expansion is nonquadraticin fermionoperatorsdue to Coulombin- fortheevolutionoperatorU(β)withsomegeneralization teraction. This last circumstance is responsible for the because we introduce the auxiliary constant λ and use appearance of the nonvanishing site cumulants called ir- λH′ instead of H′: reducibleGreen’sfunctions. ThesenewGreen’sfunctions takeintoaccountallthespin,chargeandpairingfluctua- β tions of the system. The perturbation formalism around U (β) = T exp(−λ H′(τ)dτ). (6) atomic limit has the advantagethat local(atomic)phys- λ Z ical properties can be evaluated exactly and that the 0 transferofthislocalinformationtoneighboringsitesdue In the presence of this constant we shall use index λ for tothekineticmobilityoftheconductionelectronscanbe all dynamical quantities as G (x|x′) and so on. At the λ handled perturbatively in powers of hopping integral. last stage of the calculations this constant will be put AGWT forchronologicalaveragesofproducts ofelec- equal to one. tronoperatorswas formulatedlater alsoby Metzner [14]. In zero order approximation the one-particle Green’s Metzner did not derive a Dyson-type equation for the function G(0)(x|x′) is local renormalized one-particle Green’s function and the role of our correlation function has not established. G(0)(x|x′) = δ−→x−→x′δσσ′G(0)(τ −τ′) (7) In spite of the similarity to our works, the diagrams and the Fourier representation of the last function is in his approach are quite different. The n - particle cu- mulant is represented a 2n - valent-point vertex with at- β tached n entering and n leaving lines, where as our such G(0)(iω ) = G(0)(τ)eiωnτdτ a cumulant or irreducible Green function is represented σ n Z σ by a rectangle with 2n vertices with n entering and n 0 leaving arrows. The rectangle contours the vertices with 1 e−βE0+e−βEσ e−βEσ +e−βE2 = { + }, the same site index but different time and spin labels. Z iω +E −E iω +E −E 0 n 0 σ n σ 2 InadditionMetzner investigatedthe limit ofhighlattice (8) dimensions. Z = e−βE0 +e−βEσ +e−βEσ +e−βE2, 0 In this paper we shall develop the diagrammatic the- ory proposed before for Hubbard model [11,12] with the E0 = 0, Eσ =Eσ =ǫ, E2 =2ǫ+U, aimtodemonstratethe existenceoftherelationbetween ωn = (2n+1)π/β. renormalizedquantities ofthermodynamic potential and Here E ,E and E are the energies of atomic sites. one-particle Green’s function and also to prove the sta- 0 σ 2 Ashasbeenprovedinpapers[11,12] propagator(5)has tionary properties of this potential. the diagrammatic representationdepicted on the Fig. 1. Such theorem was proved firstly by Luttinger and TheirreducibleGreen’sfunctionsofordernaredepicted Ward[15] foruncorrelatedsystemsbyusingthediagram- withrectangleswith2nvertices. Thearrowswhichenter matic technique of weak coupling field theory. inthevertexpointdepictannihilationelectronsandthat The strong coupling diagram theory used by us needs which go out the created electrons. new conceptions and new equations and they are used In paper [11] we have introduced the notion of correla- toprovestationarypropertyofthermodynamicpotential tion function Z (x|x′) which is the sum of strongly con- λ for strongly correlated systems. Such a proof has been necteddiagramscontainingirreducibleGreen’sfunctions already achieved for Anderson impurity model in paper and related to a more convenient function Λ (x|x′) = λ [16]. G0(x|x′)+Z (x|x′). λ The paper is organized in the following way. In sec- InFig. 1thefourthandseventhdiagramsoftheright- tionIIwedevelopthediagrammatictheoryinthestrong hand site belong to correlationfunction. Due to the fact coupling limit andthe skeletondiagramsare introduced. that irreducible functions arelocaland tunneling matrix −→ −→ InSectionIII weprovethe stationarytheoremforrenor- elements have the property t(x − x) = 0 in Fig. 1 malized thermodynamic potential. Section IV has the are omitted all the diagrams which contain self-locked conclusions. tunneling elements. 3 Gλ(x | x′) = +λ +λ2 − x x′ x x′ x x′ −λ2 − λ3 − λ3 − x x′ x x′ x x′ −λ3 + ... x x′ FIG.1: Firstthreeordersofperturbationtheoryforone-particlepropagator. Thethinsolidlinedepictszeroorderpropagator, thindashed linedepictsthetunnelingmatrix element. Therectangles depict irreducibletwo-particles Green’s functionsG(0)ir 2 and points are the vertices of diagrams. AsisseenfromFig. 1theprocessofpropagatorrenor- of this Figure are local and therefore their Fourier rep- malization is accompanied by the analogous process for resentation don’t depend on momentum. The last dia- tunnelingmatrixelementsrenormalizationandreplacing gramandotheroneswithmorenumberofrectanglesare oftheinstantquantityλt(x−x′)=λt(−→x −−→x′)δ(τ−τ′− non localized and their Fourier representation depends 0+)) by dynamical one T (x|x′) equal to on momentum. Only the first category of diagram are λ taken into account in Dynamical Mean Field Theory. e T (x|x′) = λt(−→x −−→x′)δ(τ −τ′−0+) As was proved in papers [11,12] the knowledge of this λ function permits to formulate the following Dyson-type e + λt(x−1)G (1|2)λt(2−x′) (9) λ equation for one-particle Green’s function X 12 Λ (k) which in Fourier representation Gλ(k)= λ−→ . (11) 1−λǫ(k)Λ (k) λ −→ 1 −→ −→−→ t(x) = ǫ(k)exp(−ik x), N −→ X−→ Here k stands for (k,iωn)with odd Matsubarafrequen- k cies. Equation (10) and (11) gives us the results 1 1 −→ G (x|x′) = G (k|iω ) λ λ n N β X−→k Xωn Tλ(k) = λTλ(k), −→ × exp[−i−→k(−→x −−→x′)−iω (τ −τ′)] e ǫ(k) n Tλ(k) = −→ . (12) 1−λǫ(k)Λ (k) λ has a form: −→ −→ Equation(12)hasaformofDysonequationfortunnelind Tλ(k|iωn) ≡ λTλ(k|iωn) Green’sfunction,andtheroleofmassoperatorΣλ iscar- −→ −→ −→ e = λǫ(k)(1+λǫ(k)Gλ(k|iωn)) (10) ried out by correlation function multiplied by auxiliary constant λ: TherenormalizedtunnelingmatrixelementT reallyis λ tunnelingGreen’sfunctionandwillbedepictedasdouble Σλ(k) = λΛλ(k). (13) dashedline. T isrepresentedbysuchdoubledashedline λ multiplied by λ. InHubbardIapproximationweneglectthecorrelation e Now we introduce the skeleton diagrams which con- functionZ (k)andconsiderfunctionΛ (k)equaltozero λ λ tainonlyirreducibleGreen’sfunctionsandsimpledashed order Green’s function G0(k). In this approximation we lines without any renormalization. In such skeleton di- have agrams the thin dashed lines are replaced by double dashed lines with realizing the complete renormalization G0(k) GI(k)= , of dynamical quantities. λ 1−λǫ(k)G0(k) The skeleton diagrams for correlation Λ function are λ depicted on Fig. 2. whichdescribestwoHubbardenergysubbandsdistanced As can bee seen from Fig. 2 the skeleton diagrams for between them for U 6= 0 and without the possibility to Λ(k) function are of two kinds. The first four diagrams describe the Mott-Hubbard transition. 4 2 2 1 1 λ2 Λ (x|x′)= − λ − 1 1 − λ x x′ x x′ 2 x x′ 3 3 2 2 1 1 λ3 2 2 − λ3 + ... − 2 6 1 1 x 3 3 x′ x x′ FIG.2: Theskeletondiagramsforcorrelation functionΛλ. DoubledashedlinesdepictthetunnelingGreen’sfunctionsTλ and therectangles depict theirreducible Green’s functions. III. THERMODYNAMIC POTENTIAL In expression (17) the coefficients 1 before each dia- n DIAGRAMS gram are absent, where n is the order of perturbation theory. These coefficients are present in Fig. 3. In order The thermodynamic potential of the system is deter- to restore these 1 coefficients in (17) and obtain the co- n c mined by the connected part of the mean value of the incidence with hUλ(β)i0 series it is enough to integrate evolution operator [11,12] by λ the expression (17) and obtain: F = F0− β1 hU(β)ic0. (14) −−→xX−→x′Xσ βZ dλt(−→x′−−→x)Gλσ(−→x −−→x′|−0+).(18) Let us consider a more general quantity first The expression (18) displayed in a diagrammatic repre- sentation coincides exactly with mean value of the evo- 1 lution operator: c F(λ) = F − hU (β)i , (15) 0 β λ 0 hU (β)ic = − βt(−→x′−−→x) λ 0 and put then λ=1. −→xX−→x′ Byusingtheperturbationtheorywehaveobtainedthe λ c firsItnoorrddeerrs toof doibatgarianmtshfeorbehUttλe(rβu)in0d,edrsetpainctdeidnginoFf itgh.es3e. × Z dλ′Gλ′σ(−→x −−→x′|−0+). (19) diagrammatic contributions we examine the expression 0 In Fourier representation we have Gλ(x|x′)λt(−→x′−−→x)δ(τ −τ′−0+)δσσ′, (16) λ Xxx′ hUλ(β)ic0 =−Z dλ′ ǫ(−→k)Gλ′σ(−→k|iωn)exp(iωn0+). −→X where double repeated indices suppose summation and 0 kσωn integration. Consequently (16) is equal to (20) − β G (−→x −−→x′|−0+)λt(−→x′−−→x) From (15) and (20) we obtain λσ −→xX−→x′Xσ λ = −λ ǫ(−→k)Gλσ(−→k|iωn)exp(iωn0+). (17) F(λ) = F0+Z dλ′ β1 ǫ(−→k) −→XkσXωn 0 −→Xkσ Xωn −→ × Gλ′σ(k|iωn)exp(iωn0+). (21) Here we have carried out the integration by time. Fromdiagrammaticpointofviewequation(16)implies Using the definition (14), the equation (21) can be writ- the procedureoflockingofthe externallines ofpropaga- ten in the form: tors Gλ diagrams depicted on the Fig. 1 with the tun- λ neling matrix element t(−→x′−−→x) and obtaining in such dλ′ 1 aonwesayforthheUd(iaβg)riacmdespwicittehdouotneFxitge.rn3a.lTlhineessestwimoilsaerriwesitohf F(λ) = F0+Z0 λ′ −→Xkσ β Xωn Tλ′(k) λ 0 diagrams differ by coefficients in front of them. × Σλ′(k)exp(iωn0+). (22) 5 1 4 1 3 1 2 1 4 hU (β)ic = −λ2 −λ3 1 3 −λ4 + λ 0 2 3 4 2 3 1 2 2 2 2 3 3 4 2 2 3 3 3 4 +λ4 +λ4 +... 2 8 4 4 1 2 1 1 1 2 FIG. 3: The first four orders of perturbation theory for hUλ(β)ic. 0 From (22) we have By using these definitions we can prove the equations F(λ) 1 δβY1(λ) λ dλ = −→Xkσ β Xωn Tλ(k)Σλ(k)exp(iωn0+) δδβTYλ′((kλ)) = −λΛλ(k)=−Σλ(k), = λΛ (k)=Σ (k), (27) 1 λ λ δT (k) = Tr(T Σ ). (23) λ λ λ β and as a result we obtain the stationary property: Inordertohaveafullsystemofequationsweaddto(22) δβY(λ) the definition of the chemical potential of the system =0. (28) δT (k) λ −→ N = 1 G (k|iω )exp(iω 0+), (24) e −→Pkσ β Pωn σ n n wBeycuasninrgewthreitedetfihneitfuionnct(i1o3n)alofYth(λe)minastsheopfoerramtor Σλ(k) 1 where N is the electron number. e 1 −→ Equation(21)establishestherelationbetweenthermo- Y1(λ) = − [ln(ǫ(k)Σλ(k)−1) β dynamic potential and renormalized one-particle propa- −→XkσXωn gator. This last quantity depends on the auxiliary pa- + T (k)Σ (k)]exp(iω 0+), (29) λ λ n rameter λ and (21) contains an additional integration over it and is awkwardbecause that. and prove the second form of stationary property We shall obtain more convenient equation for thermo- dynamic potential without such integration by λ. For δ Y(λ) =0, (30) that we shall introduce a special functional δΣ (k) λ Y(λ)=Y (λ)+Y′(λ), (25) if we take into accountthe Dysonequation for tunneling 1 function T (k). λ where Now we shall discuss the derivative by λ of the func- tional Y(λ). As was mentioned above there is twofold 1 −→ dependence of λ throughΣ (k) and explicit through the Y (λ) = − [ln(ǫ(k)λΛ (k)−1) λ 1 β λ factors λn before skeleton diagrams of Y′(λ) Fig. 4. −→X k,σ,ωn Due to the stationary property (30) we obtain + T (k)λΛ (k)]exp(iω 0+), (26) λ λ n dY(λ) δY(λ) δΣ (k) ∂Y(λ) λ = + | (31) andY′(λ)isconstructedfromskeletondiagramswithout dλ XδΣλ(k) δλ ∂λ Σλ k external lines and is depicted on Fig.4 ∂Y′(λ) ThedependenceonλinthefunctionalY′(λ)istwofold, = ∂λ |Σλ. justthroughthedependenceoftherenormalizedGreen’s functionsG ,T ,Λ andΣ andthroughtheexplicitfac- Here we have taken into account the equation (30) and λ λ λ λ tors λn in front of each diagram of Y′(λ). that that Y (λ) has not explicit dependence of λ. 1 6 Y′(λ) = −1 P P{λ − λ2 − λ3 − β −→k ωn 2 6 σ λ4 − λ4 +... } − 8 24 FIG. 4: The simplest skeleton diagrams for functional Y′(λ). Doubledashed lines are tunnelingfunctions Tλ(k). By using Fig. 4 for Y′(λ) and the definition of Λ of notion of renormalized tunneling Green’s function T in λ Fig. 2 it is easy to establish the property: addition to known before. We have defined correlation function Λ, the mass operator Σ for tunneling function ∂Y′(λ) 1 λ | = T (k)λΛ (k)exp(iω 0+) andestablishtheDysonequationforthem. Themassop- ∂λ Σλ β −→XkσXωn λ λ n eforratλor=w1a.sfound (13)to be equalto correlationfunction 1 = Tr(TλΛλ). (32) We have obtained the diagrammatic representation of β the correlation function through the skeleton diagrams Therefore from (31) and (32) we obtain: whichcontainthemany-particleirreducibleGreen’sfunc- tions G(0)ir with all possible values of n− order of per- dY(λ) 1 2n λ = T (k)Σ (k) turbationtheoryandtherenormalizedtunnelingGreen’s λ λ dλ β −→XkσXωn functions. 1 Wehaveestablishedtherelationbetweenrenormalized = Tr(TλΣλ). (33) quantities of thermodynamic potential and one-particle β Green’s function. This last function has additional de- From equation (23) and (33) we have pendence on auxiliary constant λ and has to be inte- grated over it. dF(λ) dY(λ) λ =λ , (34) We have proved the possibility to avoid such integra- dλ dλ tion by λ and to introduce the special functional Y(λ) and as the consequence we obtain constructed from skeleton diagrams. F(λ)=Y(λ)+const. (35) At the final part of the paper we have proved the sta- tionarypropertyofthis functionalandfound its relation Because for λ = 0 perturbation is absent and F(0) = to stationary of thermodynamic potential with respect F ,Y(0)=0 we have variation of mass operator or full tunneling function. 0 This theorem is the generalization of known Luttinger F(λ)=Y(λ)+F0. (36) and Word theorem [15] proved for weakly correlated sys- temsonthecaseofstronglycorrelatedsystemsdescribed Now we can put λ=1 and to obtain with the Hubbard model. F(1)=Y(1)+F . (37) 0 with the stationary property δF =0. (38) Acknowledgments δΣ IV. CONCLUSIONS Two of us (V.A.M. and L.A.D.) would like to thank We havedevelopedthe diagrammatictheoryproposed Professor N.M. Plakida and Dr. S. Cojocaru for a very for Hubbard model many years ago and introduced the helpful discussion. 7 ∗ Electronic address: [email protected] 10 Yu. A. Izyumov, F. L. Kassan-Ogly and Y. M. Skryabin, 1 J.Hubbard,Proc.Roy.Soc.A276,238(1963),A277,237 Theory of statistical mechanics of magnito-oriented sys- (1964), A281, 401 (1964). tems, Nauka, Moskow, (1987). 2 P. Fulde, Electron correlations in Molecules and Solids, 11 M.I.VladimirandV.A.Moskalenko, Teor.Mat. Fiz.82, Springer-Verlag, Berlin (1993). 428 (1990); [Theor. Math. Phys. 82, 301 (1990)]. 3 A. C. Hewson , The Kondo Problem to Heavy Fermions, 12 S. I. Vakaru, M. I. Vladimir and V. A. Moskalenko, Teor. Cambridge UniversityPress, Cambridge, England (1993). Mat. Fiz. 85, 248 (1990);[ Theor. Math. Phys. 85, 1185 4 A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg, (1990)]. Rev.Mod. Phys.8, 13 (1996). 13 N. N. Bogoliubov and V. A. Moskalenko, Teor. Mat. Fiz. 5 H. Matsumoto and F. Mancini, Phys. Rev. B55, 2095 86, 16 (1991);[ Theor. Math. Phys. 86, 10 (1991)]; Teor. (1997). Mat. Fiz. 92, 182 (1992); [Theor. Math. Phys. 92, 820 6 A. Avella, F. Mancini and R. Munzner, Phys. Rev. B63, (1992)]; DokladyANSSSR316,1107(1991);JINRRapid 245117 (2001). Communications 44, 5 (1990). 7 G.KotliarandD.Vollhardt,PhysicsToday57,53(2004). 14 W. Metzner, Phys. Rev.B43, 8549 (1991). 8 P. M. Slobodyan and I. V. Stasyuk, Teor. Mat. Fiz. 19, 15 J.M.LuttingerandJ.C.Ward,Phys.Rev.118,N5,1417 423 (1974) ”Diagram technique for Hubbard operators”, (1960). Preprint 73-136R (in Russian), Jnstitute of Theoretical 16 V. A.Moskalenko, P.Entel, L.A. Dohotaru Physics,UkranianSSRAcademyofSciences,Kiev(1973). and R. Citro, Teor. Mat. Fiz. 159, 154 (2009); 9 R. O. Zaitsev, Zh. Eksp. Teor. Fiz. 68, 207 (1974) ”Dia- [Theor. Math. Phys. 159, 550 (2009)]; Cond. mat. gram technique for Hubbard model”, Preprint IAE -2378 arXiv:0804.1651, 914 (2008). (in Russian),Jnstitute ofAtomicEnergy, Moskow (1974).