A Dynamical Quantum Cluster Approach to Two-Particle Correlation Functions in the Hubbard Model S. Hochkeppel,1 F. F. Assaad,1 and W. Hanke1,2 1 Institute for Theoretical Physics, University of Wu¨rzburg, 97074 Wu¨rzburg, Germany 2 Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 93106, USA We investigate the charge- and spin dynamical structure factors for the 2D one-band Hubbard model in the strong coupling regime within an extension of the Dynamical Cluster Approximation (DCA) to two-particle response functions. The full irreducible two-particle vertex with three mo- menta and frequencies is approximated by an effective vertex dependent on the momentum and 8 frequency of the spin/charge excitation. In the spirit of the DCA, the effective vertex is calcu- 0 lated with quantum Monte Carlo methods on a finite cluster. On the basis of a comparison with 0 hightemperatureauxiliaryfieldquantumMonteCarlodataweshowthatnearandbeyondoptimal 2 doping, our results provide a consistent overall picture of the interplay between charge, spin and n single-particle excitations. a J PACSnumbers: 71.10.-w,71.10.Fd,71.30.+h 2 2 I. INTRODUCTION χ = ] l e Two-particlecorrelationfunctions,suchasthedynam- - r ical spin- and charge correlation functions, determine a st varietyofcrucialpropertiesofmany-bodysystems. Their + Γ χ t. polesasafunctionoffrequencyandmomentumdescribe a the elementary excitations, i.e. electron-hole excitations m and collective modes, such as spin- and charge-density FIG. 1: Bethe-Salpeter equation for the two-particle propa- d- waves. Furthermore, an effective way to identify contin- gator. uousphasetransitionsistosearchfordivergenciesofsus- n o ceptibilities, i.e. two-particle correlation functions. Yet, c compared to studies of single-particle Green’s functions through the Bethe-Salpeter equation, [ and their spectral properties, where a good overall ac- 1 cord between theoretical models (Hubbard type-models) χk,p(q)=χ0k,p(q)+ χ0k,k′(q)Γk′,k′′(q)χk′′,p(q), (3) v andexperiment(ARPES)hasbeenestablished(seeRefs. kX′,k′′ 7 1,2,3,4),thesituationisusuallynotsosatisfyingfortwo- 7 particle Green’s functions. This is especially so for the which is diagrammatically depicted in Fig. 1. 3 case of correlated electron systems such as high-T su- Within the Dynamical Cluster Approximation c 3 (DCA)5,6, which we consider in this paper, one can perconductors (HTSC). The primary reason for this is . consistently define the two-particle Green’s functions, 1 that calculations of these Green’s functions are, from a 0 numerical point of view, much more involved. by extracting the irreducible vertex function from the 8 cluster. This approximation maps the original lattice To expose the problem let us consider the spin- 0 problem to a cluster of size L = l l embedded : susceptibility which is given by: in a self-consistent host. Thusc, corrcel×atiocns up to a v i 1 range ξ < lc are treated accurately, while the physics at X χ(q)= χk,p(q) ,with (1) longer length-scales is described at the mean-field level. βL r Xk,p The DCA is conveniently formulated in momentum a space, i.e. via a patching of the BZ. Let K denote χk,p(q)=hc†k, ck+q, c†p, cp q, i such a patch, and k the original lattice momentum. ↑ ↓ ↓ − ↑ The approximationboils down to restricting momentum Here, L corresponds to the lattice size, β is the inverse conservation only to the patch wave vectors K. This temperature and q (q,Ω ), q being the momentum ≡ m approximation is justified if k-space correlation func- and Ω a (bosonic) Matsubara frequency. To simplify m tions are rather structureless and, thus, in real space the notation, we have adopted a path integral coherent short-ranged. From the technical point of view, the state notation with Grassman variables: approximations are implemented via the Laue function ∆(k ,k ,k ,k ), which guarantees momentum conser- 1 β 1 2 3 4 ck,σ ≡ck,ωm,σ = √βLXr Z0 dτei(ωmτ−kr)cr,σ(τ) (2) vLaatuioenfuunpcttoioanriescirperpolcaacleldatbtiyce∆veDcCtoAr(.KIn1,tKhe2,DKC3A,K, t4h)e, thereby insuring momentum conservation only between The two-particle irreducible vertex, Γk′,k′′(q), is defined the cluster momenta. It is understood that ki belongs 2 to the patch K . lattice susceptibility reads: i To define uniquely the DCA approximation, in par- χ (q) ticular in view of two-particle quantities, it is useful to 0 χ(q)= . (6) start with the Luttinger-Ward functional Φ, which is 1 U (Q) χ (q) eff 0 − · computed using the DCA Laue function. Hence, Φ DCA is a functional of a coarse-grained Green’s function, where χ0(q) corresponds to the bubble of the dressed Gs¯e(lfK-en,ieωrgmy), ≡anGd¯(tKhe).twIror-epdaurctiibcllee qiruraendtuictiibeslesuvcehrteaxs tahree latItnicethGerfeoelnlo’swfinugncsteiocntiso:nG, w(ke)d=esGcr−0i1b(ek)1−soΣm(Ke)a.spects of calculated on the cluster and correspond, respectively, the explicit implementation of the DCA which is based to the first- and second-order functional derivatives of onthe Hirsch-FyeQMCalgorithm11. Ournew approach Φ with respect to G¯. Using the cluster irreducible DCA requiressubstantialtesting. InSec. IIIAwecomparethe self-energy, Σ(K), and two-particle vertex, ΓK′,K′′(Q), N´eel temperature as obtained within the DCA without one can then compute the lattice single-particle and lat- any further approximations to the result obtained with tice two-particle correlation functions using the Dyson our new approach. In Sec. IIIB, we present results for and Bethe-Salpeter equations. This construction of two- the temperature and doping dependence of the spin and particle quantities has the appealing property that they chargedynamicalstructurefactorsandcomparethehigh are thermodynamically consistent7,8. Hence, the spin temperature data with auxiliary field QMC simulations. susceptibility as calculated using the particle-hole cor- Finally, Sec. IV draws conclusions. relationfunctions correspondspreciselyto the derivative of the magnetization with respect to an applied uniform static magnetic field. The technical aspects of the above II. IMPLEMENTATION OF THE DCA program are readily carried out for single-particle prop- erties. However a full calculation of the irreducible two- We considerthe standardmodelofstronglycorrelated particlevertex—evenwithintheDCA—isprohibitively electron systems, the single-band Hubbard model12,13,14 expensive9 and, thus, has never been carried out. Inthepresentwork,wewouldliketoovercomethissit- uation by suggesting a scheme where the K′− and K′′− H =−Xijσ tijc†iσcjσ +UXi ni↑ni↓, (7) dependencies of the irreducible vertex are neglected. At low temperatures, this amounts to the assumption that with hopping between nearest neighbors t and Hubbard in a energy and momentum window around the Fermi interaction U. The energy scale is set by t and through- surface, the irreducible vertex depends weakly on K′ out the paper we consider U =8t. − and K′′. Following this assumption, an effective two- Our goal is to compute the spin, S(q,ω) and charge particle vertex in terms of an average over the K′ and C(q,ω) dynamical structure factors. They are given, re- − K′′ dependencies of ΓK′,K′′(Q) is introduced: spectively, by: 1 1 βLUeff(Q)=hΓK′,K′′(Q)i. (4) hSz(q,τ)Sz(−q,0)i= π Z dw e−τω S(q,ω) (8a) 1 As shown in an earlier Quantum Monte Carlo (QMC) N(q,τ)N( q,0) = dw e τω C(q,ω). (8b) − h − i π Z study by Bulut etal. (Ref. 10)for asingle QMC cluster, this is reasonable for the 2D Hubbard model (on this Here, Sz(q) = 1 eiqj(n n ) and N(q) = QMC cluster of size 8 8 with U =8t). The authors of √L j j,↑− j,↓ Ref. 10 have also calcu×lated with Ueff(Q) the effective √1L jeiqj(nj,↑+njP,↓). Thelefthandsideoftheabove electron-electron interaction for the 8 8 single QMC equaPtions are obtained from the corresponding suscepti- × cluster. Both the momentum and frequency dependence bility as calculated from Eq. (6). Finally, a stochastic werein rathergoodagreementwith the QMC results for version of the Maximum Entropy method15,16 is used to the effective electron-electroninteraction. extract the dynamical quantities. Replacing the irreducible vertex by 1 U (Q) in the In order to cross check our results, we slightly modify βL eff clusterversionoftheBethe-SalpeterEq. (3)andcarrying Eq. (6) to: out the summations to obtain the cluster susceptibility χ (q) gives: χ(q)= 0 . (9) 1 α U (q) χ (q) eff 0 1 1 − · · U (Q)= , (5) eff χ¯ (Q) − χ(Q) Here,wehaveintroducedanadditional”controlling”pa- 0 rameterαinthesusceptibilitydenominator,whichiscal- whereχcorrespondstothefullyinteractingsusceptibility culatedina self-consistentmanner. Itassures,forexam- on the DCA cluster in the particle-hole channel and χ¯ pleinthecaseofthelongitudinalspinresponse,thatχ(q) 0 is the correspondingbubble asobtainedfromthe coarse- obeys the following sum rule (a similar idea, to use sum- grained Green’s functions. Finally, our estimate of the rulesforconstructingacontrolledlocalapproximationfor 3 111111222222 (0,π) (π,π) δ=0% 111111000000 LLLLLLcccccc======888888,,,,,,ββββββtttttt======666666 δ=5.5% ~ay δ=11.8% magnetic BZ 0)0)0)0)0)0) 888888 δ=15.1% c d (π,0) ====== δ=21.0% (0,0) mmmmmm QQQQQQ,iΩ,iΩ,iΩ,iΩ,iΩ,iΩ 666666 ~ax A.F. unit cell (((((( UUUUUUeffeffeffeffeffeff 444444 FIG.3: SU(2)symmetrybrokenDCAcalculation. Left: A.F. 222222 unitcellwithnewbasisvectorsinrealspace. Right: Reduced A.F. unit cell in momentum space. 000000 ((((((000000,,,,,,000000)))))) ((((((ππππππ//////222222,,,,,,ππππππ//////222222)))))) ((((((ππππππ,,,,,,ππππππ)))))) QQQQQQ when U is localized in real space and the sum in Eq. eff FIG.2: Static(iΩm =0) irreducibleparticle-hole interaction (11) can be cut-off at a given shell. U for different cluster momentum vectors and dopings at eff Theeffectiveparticle-holeinteractionU isshownin inverse temperatureβt=6 and U =8t. eff Fig. 2 for a variety of dopings at inverse temperature βt = 6, U/t = 8 and on an L = 8 cluster, which cor- c responds to the so-called ”8A” Betts cluster (see Refs. theirreducibletwo-particlevertexhasbeenimplemented 18,19). The U -function displays a smooth momen- by Vilk and Tremblay (Ref. 17)): eff tumdependence. Theseobservationsfurthersupportthe 1 interpolation scheme (Eq. 11). Thus, indeed, Ueff is βL χ(q)=h(Siz)2i. (10) ratherlocalizedinrealspacewith sizablereductionfrom Xq its bare U = 8t value for larger doping and a further slight reduction at q = (π,π). The reduction is partly Ofcourse,αshouldbeascloseaspossibleequaltoα=1, duetotheself-energyeffectsinthesingle-particlepropa- whichis indeedwhatwewillfindafter implementing the gator, which reduce χ¯ from its non-interacting (U =0) 0 sum rule (see below). valueχ(0). Partly,italsoreflectsboththeKanamori(see Our implementation of the DCA for the Hubbard Ref. 14) repeated particle-particle scattering and vertex model is standard6. Here, we will only discuss our inter- corrections. polationschemeaswellastheimplementationofaSU(2)- Summarizing, the new approach to two-particle prop- spin symmetry broken algorithm. Since the DCA evalu- ertiesrelies ontwoapproximationswhichrenderthe cal- ates the irreducible quantities, Σ(K) as well as Ueff(Q) culation of the corresponding Green’s function possible. for the cluster wavevectors,aninterpolationscheme has Firstly,theeffectiveparticle-holeinteractionU (Q)de- eff to be used. To this, we adopt the following strategy: for pends only on the center-of-mass momentum and fre- a fixed Matsubara frequency iΩm and for each cluster quency, i.e. Q and iΩm. Secondly, χ(Q), is extracted vector Q, the effective interaction Ueff is rewritten as a directly from the cluster and χ¯ (Q) is obtained from the 0 series expansion: bubble of the coarse-grainedGreen’s functions. TogenerateDCA resultsfor the N´eeltemperature, we Ueff(Q,iΩm)= ei∆iQAi(iΩm), (11) haveusedanSU(2)symmetrybrokencode. Thesetupis Xi X∆i illustratedinFig. 3. Weintroduceadoublingoftheunit cell — to accommodate AF ordering — which in turn with i = 0,...,n 1, where n is the number of the clus- defines the magnetic Brillouin zone. The DCA k-space − ter momentum vectors Q. The quantity ∆ represents i patchingiscarriedoutinthemagneticBrillouinzoneand vectors, where each vector from the corresponding ∆ i the Dyson equation for the single-particle propagator is belongs to the same ”shell” around the origin (0,0) in given as a matrix equation: real space, i.e. 1 Gσ(k)= , (12) ∆ = 0 ;∆ = 1 , 0 G0−1(k) Σσ(K) 0 (cid:18)0(cid:19) 1 ±(cid:18)0(cid:19) ±(cid:18)1 (cid:19) − with Gσ (k) Gσ (k) ∆2 =(cid:18)±11(cid:19),(cid:18)∓11(cid:19);∆3 =±(cid:18)20(cid:19),±(cid:18)02 (cid:19).... Gσ(k)=(cid:18)Gσdccc(k) Gσdcdd(k)(cid:19). (13) ± ± With the SU(2) symmetry broken algorithm, one can WithagivenU ,Eq. (11)canbe invertedtouniquely compute directly the staggered magnetization, i.e. m = eff determine A . With these coefficients, one can compute 1 eiQj(n n ), and thereby determine the tran- theeffectivepiarticle-holeinteractionforeverylatticemo- sLitPionj temperja,↑tu−re.j,↓Since the DCA is a conserving ap- mentum vector q. This interpolation method works well proximation, the so determined transition temperature 4 0.6 8 0)0)0)0) PM−phase (symm. breaking) 7 βt=6 ==== 0000....5555 00..45 AF−pPhAMaFse−− pp(shhyaaimnsseedm e((t22.e −−brmrppeaaiarrnkttaiiiccbnlllgeee))) iΩ=0)m 65 (π,π),iΩ(π,π),iΩ(π,π),iΩ(π,π),iΩmmmm 0000000000000000................4321432143214321 T 0.3 (π,π), 4 qqqqχ(=χ(=χ(=χ(=0000 00000000 0000....00005555 0000....1111 0000....11115555 0000....2222 0.2 Q= 3 4sitescluster ( 0.1 Ueff 2 180ssiitteesscclluusstteerr 1 16sitescluster 0 0 0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2 δ δ FIG. 5: Static irreducible particle-hole interaction U for FIG.4: Phasediagramfortheone-bandHubbardmodelwith eff theclustermomentumvectorQ=(π,π). Theinsetshowsthe U = 8t with respect to different temperatures and fillings. The calculations are carried out on an Lc = 8 cluster. The sqta=ti(cπf,rπee).laTthtiecebasruescHeputbibbailridtyinχt0erfaocrtitohnesmtroemngetnhtuismUv=ect8otr. red and gray (blue and blank) objects indicate the antifer- romagnetic (paramagnetic) phase. Shading of AF and PM regions is a guide to theeye. Details are in thetext. corresponds precisely to the temperature scale at which vertex. the corresponding susceptibility, calculated without any approximations on the irreducible vertex ΓK′,K′′(Q), di- Theblue(red)trianglesindicate thetransitionline for verges. the para- to the antiferromagnetic solutions extracted from the divergent spin susceptibility (Eq. 9) within the paramagnetic calculation. A precise estimation of III. RESULTS the N´eel temperature requires very accurate results and boils down to finding the zeros of the denominator of A. Comparison of the AF transition temperature Eq. (9). In Fig. 5, we consider the effective irreducible with a symmetry broken DCA calculation particle-hole interaction U for the static case and for eff the cluster momentum Q = (π,π) relevant for the AF A first test of the validity of our new approach is a instability. As apparent, the irreducible particle-hole in- comparison with the SU(2) symmetry broken DCA cal- teraction becomes weaker with increasing doping. On culation on an Lc = 8 cluster at U = 8t. The idea is the other hand, the susceptibility χ0(q,iΩm = 0) grows to extract the N´eel temperature T from a divergence with increasing doping. At a first glance both quantities N in the spin susceptibility as calculated in the above de- Ueff and χ0 (see Fig. 5) vary smoothly as a function of scribed (paramagnetic) scheme — see Eq. (9) — and to doping. However, in the vicinity of the phase transition, compareittotheDCAresultasobtainedfromtheSU(2) signalizedbythevanishingofthedenominatorinEq. (9), symmetry broken algorithm. This comparison provides the precise interplay between Ueff and χ0 becomes deli- informationonthe accuracyofourapproximationto the catelyimportantandrendersanaccurateestimateofthe two-particle irreducible vertex (see Eq. 4). N´eel temperature difficult. Given the difficulty in deter- mining precisely the N´eel temperature, we obtain good Using the SU(2) symmetry breaking algorithm, the agreement between both methods at δ & 10 %. Note magneticphasediagramfortheone-bandHubbardmodel that in those calculations the values of α 0.86 0.97 as a function of doping is shown in Fig. 4. The para- ≈ − are required to satisfy the sum rule in Eq. (10). At (antiferro)magnetic phase transitionis indicated here by smallerdopings,andinparticularathalf-bandfillingthe gray (blank) circles. At half-filling T 0.4t and mag- N ≃ N´eeltemperature, as determined by the vanishing of the netism survives up to approximately 15 % hole doping. denominatorin Eq. (9), underestimates the DCA result. Itisknowthattheconvergenceofthemagnetizationdur- ing the self-consistent steps in the DCA approach is ex- Hence, in this limit, the K′ and K′′ dependence of the irreducible vertex plays an important role in the deter- tremely poor near the phase transition and, therefore, mination of T and cannot be neglected. wecannotestimatethetransitiontemperaturemorepre- N cisely than shown in Fig. 4. However, our precision is sufficient for comparison. We again stress that the so Let us emphasize, that a good agreement between the determined magnetic phase diagram corresponds to the N´eeltemperaturesatδ 10%andaboveisanon-trivial exact DCAresultwherenoapproximation—apartfrom achievement lending su≃bstantial support to the above coarsegraining—ismadeontheparticle-holeirreducible newschemeforextractingtwo-particleGreen’sfunctions. 5 a) b) a) S(q,ω) S(q,ω) 〈Sz(q)Sz(-q)〉 0.1 1 0.1 1 3 2 2 2.5 2 1.5 ω/t 1 ω/t 1 1 0.5 0 0 0 6 5 (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) 4 c) d) 0 1 2 3 2 3 qy C(q,ω) C(q,ω) qx 4 5 6 0 1 0.001 0.01 0.1 0.001 0.01 0.1 b) 16 16 〈Sz(q)Sz(-q)〉 12 12 3 ω/t 8 ω/t 8 2.5 4 4 2 1.5 0 0 1 (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) 0.5 0 FIG. 6: DCA (left) versus auxiliary field QMC (BSS)(right) 6 5 fHourbtbhaerddymnoamdeilcaaltsUpi/nt a=nd8,chδar≈ge1s4tr%uctaunrde βfatct=ors3.ofTthhee 0 1 2 3 2 3 4 qy BSS data on the 8×8 lattice is essentially exact and acts as qx 4 5 6 0 1 a benchmark for the DCA approach. The DCA calculations were carried out on an Lc = 8 cluster. Here we have used FIG.7: DCAa)versusBSSb)staticspincorrelationfunction α=0.98andα=1.01tosatisfythesumruleinthespinand at U/t=8, δ≈14% and βt=3 on an Lc =8 cluster. charge sectors respectively. reproducedinthe8 8QMC-BSSdata,andqualitatively B. Dynamical Spin and Charge structure factors × in the DCA results. Asafunctionofdecreasingtemperature,the DCA dy- To further asses the validity of our approach,we com- namical spin structure factor shows a more pronounced pare it to exact auxiliary-field Blankenbecler, Scalapino, spin-wave spectrum. This is confirmed in Fig. 8 on the Sugar (BSS) QMC results (Ref. 3). This method has a left hand side. Here, we fix the temperature to βt = 6 severesign-problemespeciallyinthevicinityofδ 10% and keep the doping at δ 14 % but vary the clus- and, hence, is restricted to high temperatures. ≃ ter size. As apparent, for≈all considered cluster sizes Suchacomparisonforthedynamicalspin,S(q,ω),and (L = 4,8,10,16) a spin wave feature is indeed observ- c charge, C(q,ω) dynamical structure factors is shown in able: a peak maximum at q = (π,π) is present and the Fig. 6atβt=3,δ 14%andU/t=8. TheBSSresults correct energy scale at q = (π,0) with 2J is recovered. ≈ correspond to simulations on an 8 8 lattice. Fig. 6b) Unfortunately, a direct comparison of both calculations × depicts the BSS-QMC data in the spin sector. Due to at lower temperature is not possible due to the severe short-range spin-spin correlations, remnants of the spin- minus-sign problem in the BSS calculation. density-wave are observable, displaying a characteristic The investigation of the dynamical charge correlation energy-scale of 2J, where J is the usual exchange cou- function for the above parameters shows that the DCA pling, i.e. J = 4t2. The two-particle DCA calculations calculations, which are depicted in Fig. 6 c), can also U show spin excitations with the dominant weight concen- reproducebasic characteristicsofthe BSS chargeexcita- trated, as expected and seen in the QMC data, around tion spectrum 6 d). Both calculations show excitations the AF wavevector (π,π). As apparent from the sum- at ω U which are set by the remnants of the Mott- ≈ rule, Hubbard gap. Similar results are obtained at lower tem- peratures (βt = 6) on the right hand side of Fig. 8 for 1 Sz(q)Sz( q) = dω S(q,ω) (14) differentclustersizes(Lc =4,8,10,16). Thecorrespond- h − i π Z (seeFig. 7a)),theDCAoverestimatestheweightatthis wave vector but does very well away from q = (π,π). TABLE I: Values of α for the spectra in Fig. 8. Thedispersioninthetwo-particledatahasagainahigher Fig. 8 a,b c,d e,f g,h energybrancharound2J,butitalsoshowsfeaturesatJ. α (spin) 0.99 0.92 0.93 0.97 Since the total spin is a conserved quantity, one expects α (charge) 0.98 1.00 1.00 1.00 a zero-energy excitation at q = (0,0). This is exactly 6 a) b) a) b) S(q,ω) C(q,ω) S(q,ω) C(q,ω) 0.1 1 0.001 0.01 0.1 0.1 1 0.001 0.01 0.1 2 16 2 16 12 12 ω/t 1 ω/t 8 ω/t 1 ω/t 8 4 4 0 0 0 0 (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) c) d) c) d) S(q,ω) C(q,ω) S(q,ω) C(q,ω) 0.1 1 0.001 0.01 0.1 0.1 1 0.001 0.01 0.1 2 16 2 16 12 12 ω/t 1 ω/t 8 ω/t 1 ω/t 8 4 4 0 0 0 0 (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) e) f) S(q,ω) C(q,ω) FIG.9: Spin-andchargestructurefunctionsfordifferentdop- ings: (a,b): 27%and(c,d): 32%Thecalculationsarecarried 0.1 1 0.001 0.01 0.1 2 16 out on an Lc =8 cluster at βt=6. 12 ω/t 1 ω/t 8 and change their energy scale from J = 4t2 to an en- 4 U ergy scale set by the non-interacting bandwidth. This 0 0 effect becomes even more visible with higher dopings at (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) δ = 32% (Fig. 9 c)). Furthermore, the spectrum of the g) h) chargeresponseshowsareductionoftheweightofstates S(q,ω) C(q,ω) athighenergies(ω/t 8). Thisbehaviorcorrespondsto ≈ 0.1 1 0.001 0.01 0.1 the loss of weight of the upper Hubbard band with in- 2 16 creasing doping. The corresponding equal time spin and 12 charge correlation functions of Fig. 9 c-d) are depicted ω/t 1 ω/t 8 in Fig. 10. As in auxiliary-field QMC simulations4, the 4 equaltime spin correlationfunction showsa set of peaks 0 0 at q = (π±ǫ,π) and q = (π,π ±ǫ). Here ǫ is propor- (0,0) (π,0) (π,π) (0,0) (0,0) (π,0) (π,π) (0,0) tional to the doping. The overall trend of the doping dependence of the spin- and charge-responses is in good FIG. 8: Dynamical spin and charge structure factors of the agreementwiththepreviousfindingsofQMCsimulations Hubbardmodelatβt=6,δ≈14%andU/t=8. fordifferent (Refs.2,3): thereitwasshownthatthespin-responsehas clustersizes: (a-b): Lc =4,(c-d): Lc =8,(e-f): Lc =10and acharacteristicenergyscaleω 2J andanSDW-likedis- (g-h): Lc =16. persion up to about δ 10 ≈15 % doping, despite the ≈ − fact that at these dopings the spin-spin correlations are very short-ranged(of order of the lattice parameter). A lot of the features of the two-particle spectra have ing valuesofαarelistedinTab. I. These valuesconfirm direct influence on the single-particle spectral function the overallcorrectnessofour approachin that the corre- and vice-versa. At optimal doping, δ = 14 % the spec- sponding sum rule for the charge response is accurately tralfunctionA(q,ω)inFig. 11a)showsthreedistinguish (exactly for α=1) fulfilled. features. AnupperHubbardband(ω/t 8t)andalower The doping dependence of the spin- and charge- Hubbard band which splits in an incohe≈rent background response is examined in Fig 9. Here, we restrict our and a quasiparticle band of width set by the magnetic calculations to the Lc = 8 cluster at βt = 6 and dop- scaleJ. InagreementwithearlierQMCdata(Refs.2,3), ings between δ = 14 % and δ = 32 %. At δ = 14 % weviewthisnarrowquasiparticlebandasafingerprintof (see Fig. 8) the dynamical spin structure factor dis- a spin polaron where the bare particle is dressed by spin plays a spin wave dispersion with energy scale J. That fluctuations. The fact that the dynamical spin structure is ESDW(π,0)=2J with J =4t2. As the system is fur- factor in Fig 8 c) shows a well defined magnon desper- U therdoped(δ =27%)thedispersionisnolongersharply sion at this temperature and doping, δ = 14 %, allows peaked around q = (π,π). The excitations broaden up us to interpret the features centered around q = (0,0) 7 a) a) 〈Sz(q)Sz(-q)〉 A(q,ω) 1.1 0.01 0.1 1 0.8 8 0.5 4 0.2 ω/t 0 6 5 -4 4 0 1 2 3 2 3 qy -8 qx 4 5 6 0 1 b) (0,0) (π,0) (π,π) (0,0) 〈N(q)N(-q)〉 b) A(q,ω) 0.4 0.01 0.1 1 0.3 0.2 8 0.1 4 6 ω/t 0 5 0 1 2qx 3 4 5 6 0 1 2 3 4 qy --84 FIG. 10: Static spin a) and charge b) correlation function at (0,0) (π,0) (π,π) (0,0) U/t=8, δ≈32% and βt=6on an Lc =8 cluster. c) A(q,ω) 0.01 0.1 1 and below the Fermi energy as backfolding or shadows of the quasiparticle band at q =(π,π). A comparisonof 8 the charge response spectrum in Fig. 8 d) with the cor- 4 responding single-particle spectra in Fig. 11 a) reveals that the response in the particle-hole channel at almost ω/t 0 zero energyis causedby particle-holeexcitations around -4 the quasi-particle spin-polaron band close to the Fermi -8 energy. The high energy excitations, mentioned above, areduetotransitionsfromthequasi-particlebandtothe (0,0) (π,0) (π,π) (0,0) upper Hubbard band. Asafunctionofdoping,notablechangesinthespectral FIG.11: Angle-resolvedspectralfunctionsA(q,ω)forvarious function which are reflected in the two-particle proper- hole dopings: a): 14 %, b): 27 % and c): 32 %. Calculations ties are apparent. On one hand, the spectral weight in are carried out on an Lc =8 cluster at βt=6. the upper Hubbard band is reduced. As mentioned pre- viously, this reduction in high energy spectral weight is apparent in the dynamical charge structure factor. On quantities have remained elusive due to numerical com- the other hand, at higher dopings the magnetic fluctu- plexity. In principle, the two-particle irreducible vertex, ations are suppressed. Consequently, the narrow band ΓK′,K′′(Q), containing three momenta and frequencies changes its bandwidth from the magnetic exchange en- has to be extracted from the cluster and inserted in the ergy J to the free bandwidth. This evolution is clearly Bethe-Salpeter equation. Calculating this quantity on apparent in Figs. 11 b) and c) and is in good agreement the cluster is in principle possible, but corresponds to a with previous BSS-QMC results3. daunting task which — to the best of our knowledge — has never been carried out. To circumvent this problem, we have proposed a simplification which relies on the IV. CONCLUSIONS assumptionthatΓK′,K′′(Q)isonlyweaklydependenton K′ and K′′. Given the validity of this assumption, one Two-particle spectral functions, such as spin- and can average over K′ and K′′ and retain its dependency charge- dynamical spin structure factors, are clearly on frequency and momenta, Q, of the excitation. We quantities which are of crucial relevance for the under- have tested this idea extensively for the two dimensional standingofcorrelatedmaterials. Within theDCA,those Hubbardmodelatstrongcouplings,U/t=8. Atdopings 8 δ 10 % and higher, we have found a good agreement ≈ between the N´eel temperature on an L = 8 lattice c as calculated within a symmetry broken DCA scheme with our new two-particle approach. This finding lends Acknowledgments support to the validity of our scheme in this doping range. Furthermore, we studied the doping and tem- perature dependence of the spin and charge dynamical Discussions with D.J. Scalapino about many ideas structure factors as well as the single-particle spectral of this approach in the early stages are gratefully function at βt = 6. Our results provide a consistent acknowledge. This work has been supported by the picture of the physics of doped Mott insulators, very Deutsche Forschungsgemeinschaft within the Forscher- reminiscent of previous findings within auxiliary field gruppeFOR538,bytheLeibnizRechenzentrumMunich QMC simulations. The strong point of the method, in forprovidingcomputerresources,inpartbytheNational contrast to auxiliary field QMC approaches, is that it ScienceFoundationunderGrantNo. PHY05-51164,and can be pushed to lower temperatures above and below as well as by KONWIHR. We thank M. Potthoff, E. Ar- the superconducting transition temperature. Work in rigoni,L.Martin,S.BrehmandM.Aichhornforconver- this direction is presently under progress. sations. 1 A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Rev.B 64, 195130 (2001). Phys. 75, 473 (2003). 10 N. Bulut, D. J. Scalapino, and S.R. White, Phys. Rev. B 2 R.Preuss, W.Hanke,and W.von derLinden,Phys.Rev. 47, 2742 (1993). Lett. 75, 1344 (1995). 11 J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 3 R.Preuss, W.Hanke,C.Gr¨ober, and H.G.Evertz,Phys. (1986). Rev.Lett. 79, 1122 (1997). 12 J. Hubbard,Proc. R.Soc. London A 276, 238 (1963). 4 M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 13 M. C. Gutzwiller, Phys. Rev.Lett. 10, 159 (1963). 68, 13 (1998). 14 J. Kanamori, Prog. Theor. Phys.(Kyoto) 30, 275 (1963). 5 M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Pr- 15 A. Sandvik,Phys. Rev.B 57, 10287 (1998). uschke,andH.R.Krishnamurthy,Phys.Rev.B58,R7475 16 K. S.D. Beach, 2004, arXiv:cond-mat/0403055. (1998). 17 Y.VilkandA.-M.Tremblay,J.dePhysique7,1309(1997). 6 T.Maier,M.Jarrell,T.Pruschke,andM.H.Hettler,Rev. 18 D. Betts, H. Lin, and J. Flynn, Can. J. Phys. 77, 353 Mod. Phys.77, 1027 (2005). (1999). 7 G. Baym, Phys.Rev. 127, 1391 (1962). 19 T.Maier,M.Jarrell,T.Schulthess,P.Kent,andJ.White, 8 G.BaymandL.P.Kadanoff,Phys.Rev.124, 287(1961). Phys. Rev.Lett. 95, 237001 (2005). 9 M. Jarrell, T.Maier, C.Huscroft, andS.Moukouri,Phys.