Energy and chemical potential in the two-dimensional Hubbard model

PHYSICA ELSEVIER

Physica C 250 (1995) 184-190

Energy and chemical potential in the two-dimensional Hubbard model F. Mancini

a,*

S. Marra

a, H.

Matsumoto b

a INFM e Dipartimento di Fisica Teorica e S.M.S.A. Universith di Salerno, 1-84140 Baronissi (SA), Italy b Institute for Material Research, Tohoku University, Sendai 980, Japan

Received 5 January 1995; revised manuscript received 6 April 1995

Abstract

The two-dimensional Hubbard model is studied by means of the composite-operator approach. In a generalized mean-field approximation we calculate the ground-state energy and the chemical potential, for various values of the particle density, repulsive Coulomb interaction and temperature. Comparison with the results obtained by other authors by means of numerical analysis shows a good agreement.

1. Introduction

The discovery of new materials with unusual properties not only in ordered phases but also in the "normal phase" requires the formulation of new techniques and new methods of approximation. The picture of the band model and the use of simple mean-field approximation do not seem appropriate to describe the properties exhibited by the new materials. From a theoretical point of view, the approach to the study of these new materials is based on the investigation of some models which simulate to some extension the physical systems. Among several models, the Hubbard model [1] plays an important role, since many of the models under investigation correspond to limiting cases (Heisenberg model, t - J model, t model) or to extensions ( p - d model,

* Corresponding author. E-mail address: [email protected]. unisa.it.

Kondo-Heisenberg) of the original model. Though many theoretical formulations have been proposed, most of the progresses have been achieved by the use of numerical methods (quantum Monte Carlo methods, exact diagonalization techniques); a considerable amount of work has been dedicated to this field and we now have a large basis of numerical results to which any theoretical formulation must refer [2]. In the aspect of comparison with numerical resuits, the mean-field theory of the slave-boson method [3-6] may be considered as one of the successful methods. Despite such success, additional analytic methods are sought such as the d~ method [7], composite-operator method [8], projection-operator method [9], since the slave-boson method has the problem of constraints which are satisfied as an average in mean-field theory, leading to ambiguity in the choice of slave-boson operators, and since careful consideration is necessary in distinguishing physical and unphysical excitations. Also by developing

0921-4534/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 9 2 1 - 4 5 3 4 ( 9 5 ) 0 0 3 7 0 - 3

F. Mancini et a l . / Physica C 250 (1995) 184-190

various methods, one can cross-check the validity of approximation in the region where numerical simulation is not possible. In some recent papers [8,10] we have proposed a method, the composite-operator method (COM), for the study of strongly correlated fermion systems and we have studied various models. This method aims to identify physical excitations by certain combinations of original operators, and is specially powerful when excitations extend cluster-wise. Particularly, in Refs. [11] and [12] we considered the two-dimensional single-band Hubbard model and calculated on-site quantities showing a good agreement with the data obtained by other authors by means of numerical methods. In this work we continue our study of the Hubbard model and we present theoretical results for the ground-state energy and the chemical potential, allowing the external parameters, temperature, doping and interaction coupling to vary in a large range. The comparison with numerical data gives a good agreement showing that the method proposed and the scheme of approximation adopted might be useful for the study of strong correlation in the new materials. In Section 2 we give a brief outline of the method and the general equations for the computation of the Green's functions. In Sections 3 and 4 we present results for the ground-state energy and the chemical potential and compare them with numerical data.

2. Equations The composite-operator method is based on the observation that interacting many-particle systems are more properly described not in terms of the original particles but in terms of new particles which originate by interaction. This line of thinking is not new but is as old as the study of condensed matter physics by methods of quantum field theory. However, most of the theoretical schemes are always based on some approximate treatment on the basis of the original particles. As effect of interaction (and this becomes more relevant the stronger the interaction is) the individuality of the original particle is lost and what becomes prominent is the choice of an appropriate basis of new particles whose properties

185

are determined by the system itself according to the interaction and boundary conditions. As an application of this scheme, let us consider the one band Hubbard model, described by the hamiltonian

H = ~_,tijct(i)c(j) ij +UY'.n,(i)n,(i)

-#Ed(i)c(i),

i

(1)

i

where {c(i), c*(i)} are annihilation and creation operators for the electron field at site i in the doublet notation: C ~

,

~---

As is well known, this model describes a lattice with one orbital per site: the electrons can hop among the sites and interact with each other only when are on the same site. t 0 denotes the transfer integral and describes hopping between different sites; the U term is the repulsive Coulomb interaction between two c electrons at the same site, with

n,~(i) = c*~(i)c~(i).

(3)

/~ is the chemical potential. In the nearest-neighbor approximation the hopping matrix is given by

1

tij = - 4 t a i j = - 4 t - ~

~,ei"~R'-R,)a( k),

(4)

k

where for a two-dimensional cubic lattice with lattice constant a 1 or(k) = -~ [cos( kxa ) + coS(kya)]. (5) The scale of the energy has been fixed in such a way that t , = 0. Following the ideas mentioned above, we attack the model described by the hamiltonian (1) not in terms of the original electron fields {c(i), c*(i)}, but in terms of new fields that we treat as fundamental fields. These new fields must be expressed in terms of the set {c(i), c+(i)}, since they derive by the original fields through interaction; therefore, we call them composite fields. As a convenient choice for the new basis we consider [11] the Heisenberg field n(i)

'

(6)

186

F, Mancini et al. / Physica C 250 (1995) 184-190

which is composed in terms of the original set (c(i), c*(i)} by ~:,r(i) = c,,(i) [1 - n_~,(i)],

(7)

r b ( i } =c,,(i)n_,,(i).

(8)

The dynamics of the new field is fixed by the hamiltonian (1) and we find the equation of motion 0 i~tt ~b(i )

= iS(t i -

tj)({~b(i), ~bt(j)})

+ (R[~j(i)~bt(j)]).

(17)

Introducing the Fourier transform

S( i, j) = i( ~ I Efn \ 2'rr l

dZk dto

xeik'(R'-R')-i°a(tl-t')S(k,

=j(i),

(9)

with

j ( i ) = ( -t*~(i)-4t[c~(i)+rr(i)] --( Ix- U)Tl( i) + at~r( i) 1 t*n~(,)c " ~ (i)+c(i 7r(i) =--tr 2

)[

)

(10)

S(k,

1

to)=

'

-

(k)

(19)

to)

where the kinetic t ( k ) and the dynamical Z(k, to) parts of the self-energy are given by (11)

6(k) = m(k) i-l(k),

to)=M(k, to)S-'(k, to),

n~(i)= c*(i)tr~c(i)

being the charge (ix = 0) and spin ( / , = 1, 2, 3) density operator. In this article we use the notation c"(i) to denote field operators at the first nearest sites to Ri: (12)

The properties of the system are conveniently expressed in terms of the thermal real-time retarded Green's function:

S(i, j) = (R[~b(i)~l,t(j)]).

(18)

where O B is the first Brillouin zone, Eq. (17) can be formally solved as

to -

c ~'(i)c(i)]"

c~'(i) - ~, aijc(j ). J

to),

(13)

By means of the Heisenberg equation (9), S(i, j) satisfies the equation

i-@t S( i, j) = iS(t i - tj)({~b( i),**(j)}) (14)

We

+ (R[j(i)**(j)]). decompose j(i) in the form

j(i)

= E ( - i V ) ~b(i) + 8j(i),

(15)

(20)

(21)

with

l(k) = ({~b(i), ~bt(i)})F.T , re(k) = ({j(i), S*(j)})F.'r., M(k, to) = (R[Sj(i)*t(j)])F.T..

(22) (23)

(24)

In the present stage of computation we approximate the calculation of the Green's function (13) by neglecting the dynamical contribution to the selfenergy; we call this approximation a generalized mean-field approximation, since S(i, j) is the propagator of the composite field (6). Then, the computation requires the knowledge of the normalization matrix given by (22) and of the m matrix, defined by Eq. (23). For a paramagnetic ground state a straightforward calculation gives: 1-*(k) =

n 2

0

0 n , -~

(25)

with the constraint

({Sjn(i),

*~(j)}>E.T. = 0

Then, Eq. (14) takes the form

[i~--~i-E(-iVi)]S(i, J)

(16)

m( k)

=

mll ( k )

m12( k ) l

m21(k )

m22(k) } '

(26)

mll(k) =-/-~[1 -- 2] -4t[A+a(k)(1-n+p)],

(27)

F. Mancini et al. /Physica C 250 (1995) 184-190

[

m12(k ) = m 2 1 ( k ) = 4 t

A-c~(k)

~-p

)]

4.0

/

, 3.0

(28) 2.0 n

(29)

m2z(k ) = [ - - / ~ + U]-~ - 4 t [ A + c~(k)p], where n = (ctc} is the site occupation, n = 211 - ( ~ t )

_ (,/7/t)],

(30)

and we have defined

A- (~(i)~+(i))

187

Ref.,6 ~6x~

¢

Ref. 16 (8x8)

A

R°f ~6c~0~0)

+

R e f 16 (~2x12)

/ _fl

U=4t

Ooo -1.0

,~ , ~ , , ~ -

-2.0

- (r/'~(i)r/+(i)},

o

,

,

0

(31)

,

,

I

,

,

,

,

0.5

L , 1

,

,

,

I

,

1.5

,

,

~ 2

n

1 p =- - ~ ( n ~ ( i ) n , ( i ) }

F i g . 2. T h e e n e r g y p e r site is r e p o r t e d a s a f u n c t i o n o f the b a n d filling for

- ( [ ~ r ( i ) ~ l s ( i ) ] c ' ~ t t i X ~,t t i~, X >,

(32)

U~ t

= 4 and zero temperature. The points are numeri-

cal d a t a o n d i f f e r e n t lattices f r o m R e f s . [15] a n d [16]. T h e solid line r e p r e s e n t s the C O M result.

The quantities /~, A and p depend on the external parameters n, T, t and U and are self-consistently determined by means of Eqs. (31,32) together with the requirement

the Green's function can be computed by means of Eq. (19).

(~(i)TI*(i)) = 0,

3. G r o u n d - s t a t e energy

(33)

which recovers the Pauli principle [11]. Thus, we have a complete set of equations which allow us to compute the internal parameters without making use of any further approximation [13]. It is worth to note that intersite charge, spin and pair fluctuations are taken into account by means of Eq. (32). Once the set of coupled self-consistent equations are solved,

The energy of the ground state can be directly computed from the hamiltonian (1)

E = (H} =N[-4t(c+(i)c~(i)}

(34)

where N is the number of sites and D is the double occupancy defined as O =

2.0

+ UD],

>.

(35)

Recalling the definition (13), we can immediately express E s, the energy per site, in terms of matrix elements of S(i, j) as

1.0

Es=8t[(¢~(i)¢+(i)) E

0.0

+ 2 ( sc'~(i)r/t(i)} + (r/'~(i)~/+(i)}] + U[1 - ( ¢ ( i ) ~ t ( i ) }

-I.0

-2.0

,

~ ,

,

i , 0.5

,

,

,

I , 1 n

,

,

,

I , 1.5

,

,

When we add to the Coulomb interaction the term

, 2

F i g . 1. T h e e n e r g y p e r site is r e p o r t e d a s a f u n c t i o n o f the b a n d f i l l i n g f o r v a r i o u s v a l u e s o f t h e potential

strength

U / t . The

t e m p e r a t u r e is f i x e d at t / k a T = 6. T h e p o i n t s are t h e M o n t e C a r l o d a t a o n a 12 × 12 lattice f r o m R e f . [14]. T h e s o l i d l i n e s r e f e r to t h e o r e t i c a l r e s u l t s o b t a i n e d b y m e a n s o f the c o m p o s i t e o p e r a t o r a p p r o a c h ( C O M ) . T h e e n e r g y is in u n i t s o f t.

- 2(r/(i)7/*(i)>]. (36)

U 2 Ec+(i)c(i) i

the particle-hole symmetry becomes explicit for the new hamiltonian, and we have the relation Es(2 - n) = E s ( n ) + U(1 - n).

(37)

F. Mancini et al. / Physica C 250 (1995) 184-190

188 20

2.0

u.a(coM, | 16

U.20(COM)

......

12

a

u=s fRef. 15}

E]

U=20 (Ref. 15)

T=0 K

~///a

.~

1°f

, C]"

p.~3.p c l " t ~ ' O ' ~

E 4.08'0 I

E 0.0

"

-1.0 -

tx

~-'~

a

o o -4.0

I

~

~

,

0

I

,

~

,

,

I

0.5

,

,

,

,

1

I

,

,

t

1.5

-2.0

,

2

,

~ ,

0.6

I , 0.8

,

,

n Fig. 3. Same as Fig. 2 but for U / t = 8 , 2 0 . numerical data from Ref. [15].

The points are

By making use of expression (36) we have computed the energy for different values of temperature and Coulomb potential intensity. In Fig. 1 we report the energy per site as a function of the particle density for t/k~T = 6 and U/t = 1, 2, 4; the results are compared with the data obtained in Ref. [14] by quantum Monte Carlo methods for a 12 × 12 two-dimensional lattice. Similar calculations have been performed in Ref. [5] by means of the slave-boson method. In Fig. 2 the energy per site is plotted as a function of the particle density for T = 0 K and U/t = 4; the results are compared with the data obtained in Ref. [15] by exact diagonalization techniques on a 4 X 4 cluster and in Ref. [16] by quan1.0 I 0.0

Ref 17

• ,

I

tt

I

-4.0 0.4

,

I , 1.2

, 1.4

tum Monte Carlo methods on 6 X 6, 8 x 8, 10 X 10 and 12 X 12 lattices. The same quantity is reported in Fig. 3 for T = 0 K and U/t = 8, 20 and compared with the numerical results of Ref. [15]. The limit of strong correlation is studied in Fig. 4, where the energy per particle Epart(n)=Es(n)/n is reported as a function of the particle density for U/t= 100 and kBT/t= 0.5. For comparison we report the results taken from ReL [17], where the infinite U Hubbard model on a 4 >< 4 lattice has been studied by means of a quantum Monte Carlo algorithm, specifically designed for the study of electron systems with constraint. From the results reported in Figs. 1 - 4 we can see that the agreement between our theoretical results and numerical data is very good except for a high correlation potential in a region around n = 1. In 4.0 .

0.2

,

Fig. 5. Same as Fig. 2 but for U~ t = 8 and 0.6 < n < 1.4.

2.5

0

,

13=2

I

COM ( u - I O ~ )

I o3 o

I l rl

n

0.6

u~(co.)

[

....

I

T=0 K

a ...... I

U=20t

/

/

E 1.0

, t , , , 0.8

Fig. 4. The energy per particle is reported as a function of the particle density for U / t = 100 and kaT/t = 0.5. The points are the Monte Carlo data for the infinite U case on a 4 × 4 lattice from Ref. [17].

-0.50

-2.0 [ , 0.6

I,, 0.8

,

I 1 n

,

,

,

I 1.2

1.4

Fig. 6. Same as Fig. 2 but for U~ t = 20 and 0.6 < n < 1.4.

F. Mancini et al. / Physica C 250 (1995) 184-190 2.0

189

1.0 (coM)

1.0

.............

o.0

~

I

)l

T---OK

U=4t

,~"

I

~

~

~("°'"~

I

u.~00(coM~

/

I

0.80

i | | = ="

.Ti

• B=2

n 0.60 ~t -1.0

/

/

-2.0 -3.0

6

~

" 0.40 r

0.20

/

-4.0

,~,

, k , , 0.2

0

, I , , , I , , , I , , , 0.4 0.6 0.8 n

0.0

,

,

,

I

-4

,

,

,

-2

,

,

0

I

~

,

,

2

4

~t

Fig. 7. The chemical potential in units of t is reported as a function of the band filling for U = 4 t and zero temperature. The points are numerical data from Refs. [15,16,19]. The solid line represents the COM result•

Fig. 8. The particle density is reported as a function of the chemical potential (in units of t) for zero temperature. The solid line represents the C O M result for U = 100t. The points are the Monte Carlo data for the infinite U case on a 4 × 4 lattice from Ref. [17].

Figs. 5 and 6 a more detailed study of the energy around the half-filling point is given. Our values for the energy are higher then the numerical data. Probably, this is a manifestation of the fact that in that region the paramagnetic phase is unstable against the formation of magnetic phases. This is supported by

the results reported in Fig. 4, in the limit of infinite U. Though our results are not far from the data of Ref. [17] it seems that there is no qualitative agreement, the numerical data being close to the ferromagnetic Nagaoka state [18].

3.0

n-0.4

T=0K

(COM)

n-0.6 (COM) n = 0 . 8 --

2.0 f

gx

/ &

/ +/

1.0

n=0.6

i

--

n=0.8

(COM)

A

n=0.8 (Refs.15,20)

•

n=0.6 (Rcfs. 15,20)

0

n=,O.4 (Refs. [ S,20)

X

n=0.385 (Raf,16)

O

n=0.609 (Ref.16)

+

n=0,807 (ReC. 16)

/ 0.0

>

/

,

o" "m O

o.

O

n=0.4

-1.0 X

-2.0

i i i i

I

5

I

I

I

I

I

10

I

I

I

U

I

I

15

I

I

I

I

I

20

,

,

,

•

25

Fig. 9. The chemical potential in units of t is reported as a function of the Coulomb potential for zero temperature and different values of the particle density. The points are numerical data from Refs. [15,16,19,20].

190

F. Mancini et al. / Physica C 250 (1995) 184-190

4. Chemical potential In this section we present some results of our theoretical calculation for the chemical potential and compare them with the results obtained by other authors. According to the present scheme of calculation, the chemical potential can be computed by solving the three coupled equations ( 3 1 - 3 3 ) for the three parameters/x, A and p. This has been done for different values of the external parameters n, T and U. To compare our theoretical results with the data given in the literature, in Fig. 7 we consider the case of zero temperature and U / t = 4; the data are taken from Refs. [15,16,19]. In Fig. 8 we consider k B T / t = 0.5 and U = 100t and compare our results with the numerical data of Ref. [17]. Also for the chemical potential our approach reproduces with good accuracy the results by numerical analysis. In our previous work (cf. Ref. [12]) we showed that the double occupancy D exhibits a singularity at a certain value UD(n) of the Coulomb potential; in the region U > UD(n) D is very small. This behavior implies that for U > UD(n) the energy necessary to increase the number of particles by one is almost independent on U, and the chemical potential should remain constant. Indeed, both our results and numerical analysis seem to confirm this aspect. In Fig. 9 the chemical potential is reported as a function of the interaction potential for fixed values of the particle density. Data taken from Refs. [15,16,19,20] are also reported. The agreement with the numerical results of Refs. [15] and [20] is not good for n = 0.4, but globally we can say that there is a good agreement, by considering that a precise estimation for /x by numerical analysis is not simple owing to the fact that a finite size of lattice gives a step structure in the n-/., plot. The value of U where the chemical potential starts to be constant, coincides with UD(n) reported in Ref. [12].

Acknowledgement Two of the authors (FM and SM) wish to thank R. Micnas for valuable discussions.

References [1] [2] [3] [4]

J. Hubbard, Proc. Roy. Soc. London A 276 (1963) 238. E. Dagotto, Rev. Mod. Phys. 66 (1994) 763. S.E. Barnes, J. Phys. F 6 (1976) 1375; 7 (1977) 2637. G. Kotliar and A.E. Ruckenstein, Phys. Rev. Lett. 57 (1986) 1362. [5] L. Lilly, A. Muramatsu and W. Hanke, Phys. Rev. Lett. 65 (1990) 1379. [6] R. Fr6sard, M. Dzierzawa and P. WiAlfle,Europhys. Lett. 15 (1991) 325. [7] W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62 (1989) 324; A. Georges and G. Kotliar, Phys. Rev. B 45 (1992) 6479. [8] H. Matsumoto, M. Sasaki, I. Ishihara and M. Tachiki, Phys. Rev. B 46 (1992) 3009; M. Sasaki, H. Matsumoto and M. Tachiki, Phys. Rev. B 46 (1992) 3022; I. Ishihara, H. Matsumoto, S. Odashima, M. Tachiki and F. Mancini, Phys. Rev. B 49 (1994) 1350. [9] K.W. Becker, W. Brenig and P. Fulde, Z. Phys. B 81 (1990) 165; A.J. Fedro, Yu Zhou, T.C. Leung, B.N. Harmon and S.K. Sinha, Phys. Rev. B 46 (1992) 14785; P. Fulde, Electron Correlations in Molecules and Solids (Springer, Berlin, Heidelberg, 1993). [10] F. Mancini, S. Marra, A.M. Allega and H. Matsumoto, in: Superconductivity and Strongly Correlated Electronic Systems, eds. C. Noce, A. Romano and G. Scarpetta (World Scientific, Singapore, 1994) p. 271. [11] S. Marra, F. Mancini, A.M. Allega and H. Matsumoto, Physica C 235-240 (1994) 2253. [12] F. Mancini, S. Marra and H. Matsumoto, Physica C 244 (1995) 49. [13] L.M. Roth, Phys. Rev. 184 (1969) 451; W. Nolting and W. Borgiel, Phys. Rev. B 39 (1992) 6962. [14] A. Moreo, D.J. Scalapino, R.L Sugar, S.R. White and N.E. Bickers, Phys. Rev. B 41 (1990) 2313. [15] E. Dagotto, A. Moreo, F. Ortolani, D. Poilblanc and J. Riera, Phys. Rev. B 45 (1992) 10741. [16] N. Furukawa and M. lmada, J. Phys. Soc. Jpn. 61 (1992) 3331. [17] X.Y. Zhang, E. Abrahams and G. Kotliar, Phys. Rev. Lett. 66 (1991) 1236. [18] Y. Nagaoka, Phys. Rev. 147 (1966) 392. [19] N. Furukawa and M. Imada, Physica B 186-188 (1993) 931. [20] E. Dagotto, A. Moreo, F. Ortolani, J. Riera and D.J. Scalapino, Phys. Rev. Lett. 67 (1991) 1918.