Numerical simulations of the two-dimensional Hubbard model: Dynamic properties

PHYSICA

Physica C 185-189 (1991) 104-113 North-Holland

NUMERICAL SIMULATIONS OF THE TWO-DIMENSIONAL HUBBARD MODEL: DYNAMIC PROPERTIES Douglas J. SCALAPIN0

Department of Physics, University of California, Santa Barbara, California 93106, U.S.A. Numerical techniques have provided new insight into the physical properties of the two-dimensional Hubbard mode}. Here we review what has been learned from Monte Carlo and Lanczos calculations ~.bout the dynamic properties of the repulsive and attractive Hubbard models. 1.INTRODUCTION Numerical methods can provide insight into the physical properties of interacting manyelectron systems. This is a review of some recent progress in the study of the repulsive (positive-U) and attractive (negative-U) Hubbard models based on Lanczos 1 and Monte Carlo techniques. 2 A previous review 3 described the Monte Carlo approach and gave results for a number of equal time correlation functions and zero-frequency susceptibilities. Here we will discuss what numerical results can tell us about dynamic properties such as the density of states and the optical conductivity. In addition, the dynamic response contained in the Matsubara frequency dependence of various quantities will be examined.

,e"A.F

0

#n>

I

Z

FIGURE 1 Schematic ground state phase diagram of the repulsive-U Hubbard model.

l'c T

\ 2. PHASE DIAGRAMS

o.~

\ ,% ! ¢

Figures 1 and 2 show schematic phase diagrams for the zero temperature repulsive and the finite-temperature attractive Hubbard models. Here IUl/t characterizes the strength of the interaction, and the site occupation (n} = (nit + ni!) is a measure of the band filling. Monte Carlo simulations and finite size scaling4,2 show that at half-filling, (n} = 1, the ground state of the repulsive-U Hubbard model has long-

o.,

//q

~ _

,.. -.,4,,'

;

FIGURE 2 Schematic fipite-temperature phase diagram of the attractive-U Hubbard model.

0921-4534/91/$03.50 © 1991 - Elsevier SciencePublishers B.V. All rights reserved.

D.J. Scalapino / Numericalsimulations of the two-dimensionalHubbardmodel range anti-ferromagnetic order. It is reasonable to believe that this order is present for any finite U, onsetting as e x p ( - 2 1 r v / ~ ) for small U, but simulations have only been done for U _> 2L Away from half-filling, simulations at temperatures greater than approximately 1/40 of the b a n d w i d t h show no evidence for the onset of long-range spin or pair-field order.5,G,L8 In this work U has been typically less than the bandwidth. One of the outstanding technical problems remains the "sign" problem, 0 which limits the t e m p e r a t u r e , lattice size, and U/t regime t h a t can be explored with Monte Carlo for the positive-U Hubbard model when it is doped away from half-filling. However, away from half-filling (e.g., (n) = 0.875), the pair-field and spin correlations appear to saturate at attainable temperatures and show no measurable increase with system size. 6,7 I m a d a 8 has recently used an appropriate choice of the initial trial state to obtain canonical ground state simulations on 12 x 12 lattices away from half-filling and on 16 x 16 lattices at half-filling. All of these simulations on the doped two-dimensional repulsive H u b b a r d model find only short-range spin and pair-field correlations and no evidence for phase separation.7,10 Simulations of the attractive-U Hubbard model 11 do not have the "sign" problem and are consistent with the phase diagram of Fig. lb. At half-filling the ground state exhibits both longrange C D W ,,..~ ~.a s-wave pair-field order. W h e n the system is doped away from half-filling, the C D W correlations are suppressed and the system is believed to undergo a Kosterlitz-Thouless transition into a superconducting state with power-law pair-field correlations.

t e5

3. D E N S I T Y OF STATES AND OPTICAL CONDUCTIVITY

Using Lanczos techniques, the single-particle density of states N(+)(w)for adding an electron p,r~

and N ( - ) ( w ) for removing an electron =

2

iV,n

have been calculated on a x / ~ x x / ~ lattice for both positive and negative values of U. 12 Results for (n) -- 1 and 0.8 for U/t = 8 and -8 are shown in Figs. 3 and 4 respectively. Here the solid curves correspond to the density of states for adding electrons, N(+)(w), and the dashed curves correspond to N(-)(w) for removing electrons. At half-filling, a clear gap is evident. When the repulsive-U Hubbard model is doped away from half-filling, spectral weight is removed from both bands and transferred into the gap, producing weight at the edge of the lower band just above it. If electrons are added, the curves in Fig. 3 are simply reflected about = 0. Thus the chemical potential is found to move across the gap as one goes from electron to hole doping of the half-filled system. For the attractive-U Hubbard model, the single-particle gap remains upon doping, and the chemical pont half-filling, remaining near the middle of the gap. Monte Carlo simulations provide data on imaginary time Green's functions

G(p, T) = -{Tcps(7-)c;,(O)},

(3)

D.]. $calapino / Numerical simulations of the two.dimensional Hubbard model

106

0.60

1'l,,iil,,,,-,-,i,l,,,,

N(~)

-

;

/~

I,

0.60

L'I''"I''" , , , , l l , , r l T (a)-

N(w)

(a)-

0.40

0.40

m II II II II II #I

II fl

I I II

0.2O

II II

--

P

--

II II I I

l

II Ill { I

I fl

III II'ill" i I "

'I

I

1

II

I !

I "

0.0 -10

-5

_,I,,,,I

0 i

5

-10

10

0.60

,,11-,[,I,,,,i,

B

0.40

_'1''"1"" I I

--

I II

I

~ lI |

,

-

o o

::

,,

0.20

•

Ill I

.g,

I'L

0

5

k

10

with cps(~') = eHrcvse-Hr. Numerical techniques have been developed for inverting the integral relationl3,14

/_

~o

oo

dwA(p,w)

t ,ii:" ,:,

~'I~l'I ll-I t./ 0.0 _I [ |'~I/-I' -10 -5 0

5

10

FIGURE 4

Single-particle (single-spin) density of states N(w) = N(+)(w) on a ~ x ~ / ~ cluster for U/t -- 8. The solid line corresponds to N(+)(w) and the dashed line denotes N ( - ) ( w ) (a) (n) -- 1.0, (b) (n) --.8 (after'Re'[. 13).

>0)=-

I! II m

--

FIGURE 3

G(p,r

(b)-

It I I ii I lJl I I

,-

-5

/~

-

i

-

,

-10

I0

I

--

--

E tlt,l,,r,i,

--

0.40

-

0

-5

B

(b)-

m

0.20

II II II I

'fill

II

I

0.0 0.60

I I II II i ~

0.20

e-TW

1 + e-# ~

- (4)

to obtain the spectral v-eight A ( p , w ) This requires Monte Carlo d a t a with small statistical errors and works best if moments or other sum:'de constraints are known. Surnming .'l(p,t') o v e r p g i v e s N(w). In Fig. 5, N(w) for a 4 ~" 4

Same as Fig. 3 with

U/t --- - 8

(after Ref. 13).

lattice with U/t = 8,/3 --- 4, and (n) -- 1.0 and 0.875 are plotted as the solid and dashed curves, respectively. The chemical potential for a filling of 0.875 is #z -- -2.3. Just as in Fig. 3, one sees that the effect of doping is to remove spectral weight from both the upper and lower bands that occur for (n) = 1 and transfer it into the gap, creating weight near the edge of the lower band just above #. Results for the negative-U Hubbard model 15 with U -- - 4 t at a filling of 0.87 are shown for an 8 x 8 lattice for various temperatures in FiG. 6. At high temperatures, N ( ~ ) looks ~~~~mMli~.~-'o-

dimension al den-

D.J. Scalapino / Numerical simulations of the two
0.5

',,l',,,,]",,,,l,,,,

,

r

z 0.4

107

0s~-

+,

4x4 Hubbard. Uffi8. p=4. =l.O. 0.875

-oBT. ~ - - O . t $

0.3

A3 vZ

~|0 '

-5

0

5

I! to

0.2 l

/

05

8x8. ~-4. U=-4 <~> =0.87. /a1-0.15

"':

I

0 -10

-5

5

0

t

10

0 ; -

-lO

-5

0

FIGURE 5 Density of states N(w) for a 4 × 4 lattice with U/t ---- 8 and T "- 0.25t. T h e solid curve is for (n) -- 1 and the d a s h e d curve is for (n} = 0.875, with # = - 2 . 3 t . (After Ref. 15.) sity of states. As t h e t e m p e r a t u r e is lowered, and the s u p e r c o n d u c t i n g pair-field correlations develop, a gap is seen to open in N(w). This gap opens arou.nd/z = - 0 . 1 1 , and the a s y m m e try in N(w) reflects t h e Van Hove peak in the non-interacting density of states. T h e real part of t h e optical conductivity in t h e ground state ]¢0) is al(w)

I(¢"IJ=10°)12

= D6(w) + - ~ E he0

En

-

E0

6(w-(En-Eo)).

(5)

Jx is the p a r a m a g n e t i c current operator Jz "- it ~i~(c!+~jc~s- c!sci+xs) a n d D is the IIere

D r u d e weight o b t a i n e d from 16 2~-e2 = ~]V ( ¢ 0 1 - Tl~b0) - ~ E rt#0

{(¢=1J~l¢°))2 V :T0 ' (6)

~a 0.5

z

$

I0

8x8. B=8, U = - 4

]

=0.87. ~=-0.15

l

J o -I0

-5

0

5

tO

FIGURE 6 T h e density of states N ( ~ ) versus w for U/t -- 4 a n d (n) -- 0.87 at different temperatures, from Monte Carlo d a t a on an 8 x 8 lattice. Here w is m e a s u r e d relative to t h e chemical potential which is - 0 . 1 1 at low t e m p e r a t u r e s . (After Ref. 10.) with T t h e kinetic energy. Lanczos results 17 for cq(w) on a 4 x 4 lattice with U/t = =t=10and several values of the doping (n.) = 1 - x are shown in Figs. 7 and 8. The Drude weight D versus (n) for different values of U/t and lattice size are shown in Fig. 9 for positive U and Fig. 10 for negtive U. For tb~' repulsive-U H u b b a r d model, the Drude weight at half-filling is small and negative. We expect t h a t just as in the case of a hail-filled

D3. Scalapino / Numerical sim.,_#afionsof the two-dimensional Hubbard model

108

0.8

5O

n

0.6

20 0.4

I m

m

10

m

0.2 0.0

0 0

5

10

15

20

-5

0.20

0

5

10

15

20

D

30 (b)

-

(b)

0.15

m

20

-

x=0.125

-

0.10

_

I

m

0.05

10

0.0 -5

0 0

5

10

15

20

10

5

15

20

60 FIGURE 8

FIGURE 7 al(•) for a 4 x 4 cluster with U/t = 10 doping x = 1 - (n), (a) x = 0.0, and (b) 0.125. A small shift e = 0.01 from the real was used to plot the 6-functions (after Ref.

0

and x = axis 17).

one-dimensional Hubbard ring, 18,19 D will go exponentially to zero with the size of the system. As the repulsive-U Hubbard model is doped, spectral weight is rapidly transferred from the Hubbard o_,_~nto th,~ .... Drude peak and the midgap region, Fig. 7. The spectral weight transferred below the gap shows the same qualitative dependence on doping t h a t is seen experimentally. 2° T h e negative-U Hubbard model has a finite Drude weight at half-filling, corresponding to

Same as Fig. 7 but for U/t = - 1 0 .

the spin stiffness of the half-filled positive-U Hubbard model, which initially increases with doping as the superconducting correlations dominate. As IUl/t becei,ics l~rge, the effective mass !U!/t 2 of the pairs increases and D decreases. In addition to the excitation of quasiparticles across the gap, a low-lying collective excitation exists in the gap because of the absence of a long-range Coulomb interaction.

D.J. Scalapino / Numerical simulations of the two.dimensional Hubbard model

!09

m

0.40

0.40

Drl

0.50

0.30

0.20

0.20

0.10

0.10

0.0

1.0

0.8

0.6

0.4

0.2

0.0

0.0 1.0



0.8

0.6

0.4

0.2

0.0



FIGURE 9

F I G U R E 10

Drude weight for the repulsive Hubbard model on a 4 x 4 lattice normalized by 2ze 2, Dn -" D/2~e 2, versus x for V / t = 4 ( A ), 8 ( m ), and 2 0 ( • ) o n a 4 × 4 1 a t t i c e . Q ,~} , a n d O indicate results for a v / ~ x ~ site cluster at U/t = 8, 20, and 100, respectively. The solid line without points is for U/t = 0 in the bulk limit. (After Ref. 17.)

Drude weight fcr the attractive 4 x 4 Hubbard model with U/t = - 4 ( ~ , ) , - 8 ( A , ) , - 1 0 ( l-I ), and - 2 0 ( X ). g indicates results for U/t = - 1 9 0 on a v / ~ x ~ lattice. (After Ref. 17.)

4. MATSUBARA F R E Q U E N C Y DEPENDENCE The susceptibility x(q, iwm = 0) for the repulsive-U Hubbard model, with U/t = 4, (n) = 0.87, and T = 0.25t, obtained from Monte Carlo simulations on a 12 x 12 lattice, is plotted (solid points) versus q iti Fig. 11. Here one clearly sees the anti-ferromagnetic peak. Although a

12 × 12 lattice is limited in its q resolution, the relative peak heights suggest that the maximum is shifted from (~, n ) and moves do~n the edges to (w, u - 5) and (Tr - 6, 7r). The open squares, Fig. 11 are RPA results 21 for (n) = 0.86 and UReA = 2,

xo(q, ix-) XaeA(q, iwm) = 1 - UaP,,,xo(q, iwm)' with

(7)

1)3. Scalapino / Numerical fimulations of the two-dimensional Hubbard model

110

I

'

I

'

I

'

I

'

12x12 Lottice, T-0.25t, ==0.87 •

QMC, 8x8, U==4t, n---0.87, T=0.2Ot

Z.5

QMC: U,=4t

o RPA: U=2t

q=(w,w) A

~r •-

1.5

i 0.5

(o.o)

[

I

(,~.o)

(,.,)

(o.o)

T h e solid points are static magnetic susceptibility x ( q , 0 ) versus q for U/t = 4, {n) = 0.87, a n d T = 0.25t on a 12 x 12 lattice. The open points with t h e line to guide t h e eye are RPA results, Eq. (7) with UspA = 2. (After Ref. 21.)

=

p ia,m +

- Sp+q)"

(8)

Here e = - 2 t ( c o s p x + cospy) - #. In

Fig.

Monte Carlo d a t a for x(q = (lr, lr),iwrn) is p l o t t e d versus the M a t subara frequency win. T h e open squares show the R P A form, Eq. (7). T h e characteristic spin fluctuation energy WSF required to account for the decay shown in Fig. 12 is of order 0.3t at this t e m p e r a t u r e . T u r n i n g next to the single-particle Green's function,21 -

-

1

(9)

one can determine E(p, iwn) from Monte Carlo d a t a on G(p,7). Writing E(p, iwn) in terms of odd a n d even iwn terms, E(p, iwr~) =

I

2

,

il

B

4

II

6

, el

m,

8

10

F I G U R E 12 The solid squares are X((~r, ~r), iwm) on an 8 x 8 lattice for U/t = 4, (n) = 0.87, and T = 0.2t versus iwm. T h e open squares are RPA results. (After Ref. 21.) with Z a n d ~ even functions of iWn, we can extract the wave function renormalization factor z-l(p,i~n). In Fig. 13 we have plotted

Z(p, iwn) for p on the "fermi surface" defined

12,

G(p, iwn) =

l

~.It

F I G U R E 11

x0(q,i m)

0

o

q

(1-Z(p, iwn))iwn+f~(p, iwr~), (10)

by p such t h a t (np) is closest to 0.5. Figure 14 shows t h e "fermi surface" d e t e r m i n e d in this manner. 22 In Fig. 13, for comparison, we also shown Z(p, iwn) for the half-filled case, with p = (lr,0) a n d T = OAt. Here, with a gap, we expect t h a t Z(p, iwn) ~- A2/(Wn) 2. For (n) = 0.87 and T = 0.2t, Z(PF,iwn) a p p e a r s to be well behaved, consistent with a fermi liquid. Clearly it is necessary in this case to examine Z(p, iwn) at lower t e m p e r a t u r e s on larger lattices, and a finite-size scaling analysis is needed to determine w h e t h e r the system remains a fermi liquid as T goes to zero. We have recently been able to calculate t h e m o m e n t u m and M a t s u b a r a frequence dependence of t h e four-point scattering amplitude. 23

D.I. Sca.lq_pino / Numerical simulations of the two.dimensional Hubbard model

'

I

'

I

'

l|I

71

I

QMC. 8x8, U=4t, n=0.87, T=O.20t 1.5

p=(3=/4.0) o m

m

m

®

m

I

{B m

m

~

m

m

§

es N

0.5

I

,

0

0

I

I

,5

I

,

10

,

15

20

e,,/t I0

I

'

I

I

-71

I

kx OMC. 8x8. U,=4t. n=1.0. T=O.10t

8-

F I G U R E 14 T h e Brillouin zone for a 16 x 16 lattice, with the solid squares marking k-points for which (nit) > 0.5 for a band filling of (n) = 0.87 with U/t -- 4 and ~ = 6/t. Extrapolating the Monte Carlo data for (nk) on the 16-site !attice to the points in the Brillouin zone where {r,k) would equal 0.5 gives the crosses ( ;~ ). The solid line is the non-interacting (U = 0) fermi surface for (n) = 0.87. (From Ref. 23.)

p=(.,o) 6

a o.-~

Cb)

N

4

2 0

0

,

0

@

I 2

e

@

,

O

I

@

@

,

4

O

I 6

n"

O

O

~

@

I 8

'

O

S

i I0

¢odt F I G U R E 13

Z(pF, iwn) versus iwn on an 8 x 8 lattice with V l t = 4. (a) (n) ---- 0.87, T = O . 2 t , (b) (n) = 1.0, T = 0.1t. (After Ref. 21.)

M a t s u b a r a energy transfer is 27rT, which is well above t h e characteristic spin fluctuation frequency, a n d the scattering is suppressed. T h e interaction shown in Fig. 15 favors d-wave pairing b u t is well away from forming a bound state at t h e t e m p e r a t u r e s and lattice sizes that have been simulated.

Figure 15 shows t h e particle-particle vertex F(p 1, iWn,lp, iwn) in t h e singlet channel with zero center of mass m o m e n t u m and energy (see the d i a g r a m in the inset). In these plots, p is set to (zr, 0) and pt moves a r o u n d the non-interacting fermi surface. In Fig. 15a, the M a t s u b a r a energy transfer is zero, a n d the strong spin-fluctuation repulsive interaction at a m o m e n t u m transfer of (.x, zr) is clearly seen. In Fig. 15b, the

5. C O N C L U S I O N It is possible to obtain information about dynamic quantities for the Hubbard model using Lanczos techniques on small clusters and, if tight error control is achieved, with Monte Carlo m e t h ods on larger lattices. !n addition, the dependence on Matsubara frequencies can provide in-

D2. Scalapino / Numerical sim,_d~_~onsof the two-dimensional Hubbard model

112

I°1

'

'

'

I

Q.

o v

a.

r.*

--2

QMC, 4x4, U=4t, =0.87, T==O.50t

-4 p=(.,O), ~.=='l,rT, ~.'=wT, Qcu=O

-6 -8

-I0

(.,o)

I

I

I

(o,=)

(-.,o)

(o.-.)

(-,o)

p' I0

I

I

I

(b) 4 3 o.

o "

The work reviewed here represents the contributions of a number of faculty, postdoctoral researchers, and graduate students. I would particularly like to acknowledge N.E. Bickers, N. Bulut, E.Dagotto, R. Fye, A. Moreo, R. Noack, R. Scalettar, R. Sugar, and S.R. White. I would also like to thank L. Wilson for all of her effort and skill in putting this manuscript together. This work was partially supported by the Department of Energy under grant DE-FG0385ER45197, the National Science Foundation under grant DMR90-02492, and the Electric Power Research Institute. The numerical calculations reported in this paper were performed at the San Diego Supercomputer Center, NMFECC, and the National Center for Supercomputing Applications, Urbana, Illinois.

2

° v

ACKNOWLEDGMENTS

--2 r.."

QMC, 4x4, U=4t, =0.87, T=O.50t

-4

-6-

-10

(.,o)

p:(fr,O), ¢~.=TrT, ~.'=3vrT, Qo,=O

I

I

(o,.)

(-.,o)

(o,-.)

(.,o)

p°

FIGURE 15 Singlet particle-particle scattering vertex for U / t = 4, (n) = 0.87, and T = 0.5t, ico.) versus p' with p = (rr,0). (a)ico n, = i W n = 7rT, (b) icon, = 3~rT, icon = rrT (after Ref. 22).

REFERENCES 1. For a recent review see E. Dagotto, Int. J. Mod. Phys. B5 (1991) 77. 2. The Monte Carlo algorithm we use is described in S.R. White, D.J. Scalapino, R.L. Sugar, E.Y. Loh, Jr., J.E. Gubernatis, and R.T. Scalettar, Phys. Rev. B40 (1989) 506. 3. D.J. Scalapino, Results from numerical simulations of the 2D Hubbard Model, in: High Temperature Superconductivity Proceedings, eds. K.S. Bedell, D. Coffey, D.E. Meltzer, D. Pines, and J.R. Schrieffer (Addison Wesley, 1990) pp. 314-372. 4. J.E. Hirsch and S. Tang, Phys. Rev. Lett. 62 (1989) 591. 5. J.E. Hirsch and H.Q. Lin, Phys. Rev. B37 k..,,~,..,j

sight into characteristic response frequencies which enter susceptibilities, self-energies, and vertex functions.

uv,

v.

6. M. Imada and Y. Hatsugai, J. Phys. Soc. Jap. 58 (i989) 3752. 7. A. Moreo and D.J. Scalapino, Phys. Rev. B 43 (1991) 8211. 8. M. Imada, J. Phys. Soc. Jap. 60, No. 8 (1991). 9. E.Y. Loh, Jr., et aI., Phys. Rev. B 41 (1990) 9301.

D..?. Scalapino / Numerical simulations of the two.dimensional Hubbard model

10. H.-Q. Lin, Phys. Rev. B, to be published. 11. R.T. Scalettar et al., Phys. Rev. Lett. 62 (1989) 1407; A. Moreo and D.J. Scalapino, Phys. Rev. Lett. 66 (1991) 946. 12. E. Dagotto et al., preprint UCSBTH-91-22. 13. R.N. Silver, J.E. Gubernatis, D.S. Sivia, and M. Jarrell, Phys. Rev. Lett. 65 (1990) 496. 14. S.R. White, to appear in Computer Simulations in Condensed Matter Physics III, eds. D.P. Landau, K.K. Mon, and H.-B. Schfittler (Springer-Verlag, Heidelberg, 1990). 15. A. Moreo, D. Scalapino, S. White, preprint. 16. B. Shastry and B. Sutherland, Phys. Rev. Lett. 65 (1990) 243. 17. E. Dagotto et al., preprint NSF-ITP-91-58. 18. W. Kohn, Phys. Rev. 133 (1964) A171. 19. R.M. Fye et al., Phys. Rev. B (1991). 20. S.L. Cooper et al., Phys. Rev. B 41 (1990) 11605. 21. N. Bulut, D.J. Scalapino, and S.R. White, Proc. IWEPM 1991. 22. A. Moreo, et al., Phys. Rev. B 41 (1990) 2313. 23. S.R. White, N.E. Bickers, N. Bulut, and D.J. Scalapino, preprint.

i !3