The rotationally invariant approximation for the two-dimensional t-J model

18 Pages·0.25 MB·English

by  A. Sherman

The rotationally invariant approximation
for the two-dimensional t-J model
A. Sherman
2
0 Institute of Physics, University of Tartu, Riia 142, 51014 Tartu, Estonia
0
2

n
M. Schreiber
a
J Institut fu¨r Physik, Technische Universita¨t, D-09107 Chemnitz, Federal Republic of Germany

9
and School of Engineering and Science, International University Bremen, Campus Ring 1,
2
D-28759 Bremen, Federal Republic of Germany

]
(February 1, 2008)
l
e
-
r
t
s
. Abstract
t
a
m
-
d UsingthedescriptionintermsoftheHubbardoperatorsholeandspinGreen’s
n functions of the two-dimensional t-J model are calculated in an approxima-
o
tion which retains the rotation symmetry of the spin susceptibility in the
c
[ paramagnetic state and has no predefined magnetic ordering. In this approx-

1 imation, Green’s functions are represented by continued fractions which are
v
interrupted with the help of the decoupling corrected by the constraint of
6
2 zero site magnetization in the paramagnetic state. Results obtained in this
5
approach for an undoped 32×32 lattice (the Heisenberg model) and for one
1
0 hole in a 4×4 lattice are in good agreement with Monte Carlo and exact
2 diagonalization data, respectively. In the limit of heavy doping the hole spec-
0
trum described by the obtained formulas acquires features of the spectrum of
/
t
a weakly correlated excitations.
m
PACS numbers: 71.10.Fd, 71.27.+a, 74.25.Ha, 74.25.Jb
-
d
n
o
c
:
v
i
X
r
a
Typeset using REVTEX
1
I. INTRODUCTION
The two-dimensional t-J model is one of the most frequently used models for the de-
scription of CuO planes of perovskite high-T superconductors (for a review, see Ref. 1).
2 c
Together with the numerical methods – the exact diagonalization of small clusters,2,3 Monte
Carlo simulations4 and the density-matrix renormalization-group technique5 – a number of
analytical methods, such as the mean-field slave-boson6 and spin-wave approximations, was
used for the investigation of the model. The latter method which is based on the spin-
wave description of the magnetic excitations was shown to be remarkably accurate in the
case of small hole concentrations and zero temperature.7 This approach was extended to
the ranges of moderate hole concentrations and finite temperatures,8 in particular with the
use of the spin-wave approximation modified9 for short-range order.10 The positions, sym-
metry and size of the pseudogaps in the hole and magnon spectra, values of the magnetic
susceptibility and spin-lattice relaxation rates obtained in this approach are close to those
observed in photoemission, spin-lattice relaxation and neutron scattering experiments on
cuprate perovskites.10,11
The apparent shortcomings of the spin-wave approximation of the t-J model are the vi-
olation of the rotation symmetry of the spin susceptibility components in the paramagnetic
state, the predefined magnetic ordering in the N´eel state which serves as the reference state
of the approximation, and the neglect of the kinematic interaction. In this paper we try
to overcome these shortcomings by using the description in terms of Hubbard operators.
Green’s functions constructed from these operators are calculated with the use of the con-
tinued fraction representations following from the Mori projection procedure.12 To interrupt
these otherwise infinite continued fractions we use decouplings of the higher-order Green’s
functions arising in later stages of this calculation procedure. Following the idea of Ref. 13,
a correction parameter is introduced in these decouplings to fulfill the constraint of zero
site magnetization in the paramagnetic state. In this state the obtained components of the
spin Green’s functions are rotationally invariant. The self-energy equations are similar in
their form to the equations derived in the modified spin-wave approximation.10 In the case
of heavy doping the pole in the hole Green’s function corresponds to a weakly correlated
nearest-neighbor band. To check the validity of the obtained equations in the opposite case
of light doping we have performed calculations for conditions which allow comparison with
exact diagonalizationand Monte Carlo results. We foundgoodagreement of our results with
the results of Refs. 2,14 for spin correlations in an undoped 32×32 lattice and for the hole
spectral function of a 4×4 lattice with one hole. To gain a notion of the spectral function
in larger lattices it was calculated in a 20×20 cluster.
II. DESCRIPTION OF THE MODEL
The Hamiltonian of the two-dimensional t-J model reads
1
H = tnma†nσamσ + 2 Jnm sznszm +s+n1s−m1 +µ Xn, (1)
nXmσ Xnm (cid:16) (cid:17) Xn
where anσ = |nσihn0| is the hole annihilation operator, n and m label sites of the square
lattice, σ = ±1 is the spin projection, |nσi and |n0i are site states corresponding to the
2
absence and presence of a hole on the site. If Hamiltonian (1) is obtained from the extended
HubbardHamiltonian,15 thesestatesarelinearcombinationsoftheproductsoftherespective
3dx2−y2 copper and 2pσ oxygen orbitals.10 Wetake into account nearest neighbor interactions
only, tnm = t aδn,m+a and Jnm = J aδn,m+a where t and J are hopping and exchange
constants and the four vectors a connect nearest neighbor sites. The spin-1 operators can
P P 2
be written in the Dirac notations as sz = 1 σ|nσihnσ| and sσ = |nσihn,−σ|. The
n 2 σ n
chemical potential µ is included into Hamiltonian (1) to control the hole concentration.
P
Xn = |n0ihn0|. The term −J8 naXnXn+a is frequently included into Hamiltonian (1). For
problems considered below this term leads to an unessential renormalization of the chemical
P
potential and therefore it is omitted. The operators anσ, szn, sσn, and Xn are the Hubbard
operators in the space of states of the t-J model.
The states |nσi and |n0i satisfy the following completeness condition:
|nσihnσ|+|n0ihn0| = 1. (2)
σ
X
Usingthisconditionandtheaboveexpressionforsz theconstraintofzerositemagnetization,
n
which has to be fulfilled in the paramagnetic state, can be reduced to the form
1
hszi = (1−x)− s−1s+1 = 0, (3)
n n n
2
D E
where angular brackets denote averaging over the grand canonical ensemble and the hole
concentration x = hXni in the homogeneous state. It should be noticed that in accord
with the Mermin-Wagner theorem16 the long-range antiferromagnetic ordering is destroyed
for any nonzero temperature in the two-dimensional system. Therefore the fulfillment of
constraint (3) has to be ensured for the considered states.
The above operators satisfy the following commutation (anticommutation) relations:
s−n1,s+m1 = −2sznδnm, [sσn,szm] = −σsσnδnm,
h i 1
[anσ,szm] = − σanσδnm, anσ,sσm′ = −an,−σδnmδσ,−σ′,
2
h i
(4)
anσ,a†m,σ′ = 1−s−nσsσn δnmδσσ′ +sσnδnmδσ,−σ′, {anσ,amσ′} = 0,
ns−n1,Xm =o 0,(cid:16) [szn,Xm](cid:17)= 0, [anσ,Xm] = anσδnm.
h i
Notice that the hole creation and annihilation operators do not satisfy the fermion anticom-
mutation relations. This is the consequence of the exclusion of doubly occupied site states
due to the strong on-site repulsion [see Eq. (2)].
III. CONTINUED FRACTION REPRESENTATION OF GREEN’S FUNCTIONS
To investigate the energy spectrum and magnetic properties we shall calculate the hole
and spin retarded Green’s functions
G(kt) = akσ a†kσ = −iθ(t) akσ(t),a†kσ , D(kt) = szk sz−k = −iθ(t) szk(t),sz−k ,
t t
DD (cid:12) EE Dn oE DD (cid:12) EE Dh iE
(cid:12) (cid:12) (5)
(cid:12) (cid:12)
3
where akσ = N−1/2 nexp(−ikn)anσ, szk = N−1/2 nexp(−ikn)szn, N is the number of sites
and akσ(t) = exp(iHPt)akσexp(−iHt). In the consiPdered states G(kt) does not depend on σ.
To calculate the above Green’s functions we use their continued fraction representations
which can be obtained using the Mori projection operator technique.12 Let us consider the
inner product A·B† of the operators A and B which is defined in such a manner that the
following cond(cid:12)itions a(cid:12)re fulfilled: i) (aA+bB)·C† = a A·C† + b B ·C† , a and b are
(cid:12) (cid:12)
(cid:12) (cid:12)
∗
arbitrary numbers; ii) [A,H]·B† =(cid:12)(cid:12) A· H,B† ; ii(cid:12)(cid:12)i) A·(cid:12)(cid:12)B† =(cid:12)(cid:12) B ·A(cid:12)(cid:12) † . W(cid:12)(cid:12) e notice that
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
the inner products defi(cid:12)ned as A(cid:12),B†(cid:12) , h A,B†i(cid:12) , an(cid:12)d (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
Dn oE Dh iE
∞
A,B† = i dte−ηt A(t),B† , η → +0 (6)
0
(cid:16) (cid:17) Z Dh iE
satisfy the above properties. Let us divide the result of the commutation of some operator
A with the Hamiltonian into longitudinal and transversal parts with respect to A . The
0 0
transversal part A is determined as an operator the inner product of which with A is equal
1 0
to zero. Thus,
[A ,H] = E A +A , (7)
0 0 0 1
where E is determined from the condition A ·A† = 0,
0 1 0
(cid:12) (cid:12)
(cid:12) (cid:12)
(cid:12) † (cid:12) † −1
E = [A ,H]·A A ·A .
0 0 0 0 0
(cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
Given A0 and E0, the operator A1 m(cid:12)ay be found(cid:12)f(cid:12)rom Eq(cid:12). (7). The commutator of A1 with
the Hamiltonian will contain already three terms,
[A ,H] = E A +A +F A .
1 1 1 2 0 0
The coefficients E and F and the new operator A are determined with the use of the two
1 0 2
orthogonality conditions A ·A† = 0, i = 0,1,
2 i
(cid:12) (cid:12)
(cid:12) (cid:12)
E = [A ,H]·A† A (cid:12)·A† −1,(cid:12) F = [A ,H]·A† A ·A† −1 = A ·A† A ·A† −1,
1 1 1 1 1 0 1 0 0 0 1 1 0 0
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
where w(cid:12)e have used(cid:12)th(cid:12) e prop(cid:12)erties of the(cid:12)inner produ(cid:12)c(cid:12)t. This(cid:12)proce(cid:12)dure can(cid:12) (cid:12)be cont(cid:12)inued.
In each step of it the coefficients and a new operator are determined by the conditions of
the orthogonality of this operator to all operators obtained previously. Using the properties
of the inner product it can be shown that only the n-th, (n+1)-th and (n−1)-th operators
appear in the commutator of the n-th operator with the Hamiltonian,17
[A ,H] = E A +A +F A ,
n n n n+1 n−1 n−1
(8)
−1 −1
E = [A ,H]·A† A ·A† , F = A ·A† A ·A† .
n n n n n n−1 n n n−1 n−1
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)
As can be seen, alg(cid:12)orithm (8) i(cid:12)s(cid:12)the mod(cid:12)ification of the(cid:12)Lanczos(cid:12) (cid:12)orthogonaliz(cid:12)ationprocedure
which is well known in computational mathematics (see, e.g., Ref. 18 and references therein)
4
and physics.19 In Eq. (8), operators A play the role of mutually orthogonal wave functions
n
or vectors in the usual Lanczos procedure.
Following Mori12 we introduce the projection operator P which projects an arbitrary
n
operator Q on the operator A ,
n
−1
P Q = Q·A† A ·A† A ,
n n n n n
(cid:12) (cid:12) (cid:12) (cid:12)
and determine the time evolution of th(cid:12)e oper(cid:12)at(cid:12)or A b(cid:12)y the equation
(cid:12) (cid:12) (cid:12) n (cid:12)
d n−1
i A = (1−P )[A ,H], A = A , (9)
nt k nt n,t=0 n
dt
k=0
Y
Due to the projection operators in Eq. (9) this time dependence of the operator differs from
the conventional one except the dependence of A . To underline this difference we use the
0
subscript notation for the time dependence in Eq. (9). Notice also that in accord with this
equation A remains orthogonal to operators A , i < n for t > 0. Let us divide A into
nt i nt
two parts,
−1
A = R (t)A +A′ , R (t) = A ·A† A ·A† .
nt n n nt n nt n n n
(cid:12) (cid:12)(cid:12) (cid:12)
From this definition it follows that A′ = (1−P )A(cid:12) . Equa(cid:12)t(cid:12)ions (8)(cid:12)and (9) determine the
nt n (cid:12)nt (cid:12)(cid:12) (cid:12)
time evolution of this operator,
d n
i A′ = R (t)A + (1−P )[A′ ,H].
dt nt n n+1 k nt
k=0
Y
Solving this equation we find
t
A = R (t)A −i dτR (τ)A .
nt n n n n+1,t−τ
0
Z
This result allows us to obtain the following equation for the functions R (t):
n
d t
i R (t) = E R (t)−iF dτR (τ)R (t−τ).
n n n n n n+1
dt 0
Z
After the Laplace transformation R (ω) = −i ∞dt exp(iωt)R (t) this equation reads
n 0 n
R −1
R (ω) = [ω −E −F R (ω)] . (10)
n n n n+1
If the inner product is defined as the average of the commutator (anticommutator)
†
of operators, the function R (ω) = R (ω) A ·A coincides with the Fourier transform,
0 0 0 0
A A† = ∞ dtexp(iωt) A A† , o(cid:12)f the c(cid:12)ommutator (anticommutator) retarded
0 0 −∞ e 0 0 (cid:12) (cid:12)
ω t (cid:12) (cid:12)
DGDree(cid:12)(cid:12)n’sEfEunctioRns of the type DoDf Eq(cid:12)(cid:12). (E5E). If the inner product is defined by Eq. (6), the
funct(cid:12)ion R (ω) coincides with Kubo(cid:12)’s relaxation function
0
∞ ∞
A eA† = dteiωt A A† , A A† = θ(t) dt′ A (t′),A† . (11)
0 0 0 0 0 0 0 0
(cid:16)(cid:16) (cid:12) (cid:17)(cid:17)ω Z−∞ (cid:16)(cid:16) (cid:12) (cid:17)(cid:17)t (cid:16)(cid:16) (cid:12) (cid:17)(cid:17)t Zt Dh iE
(cid:12) (cid:12) (cid:12)
From(cid:12)Eq. (10) for all these fun(cid:12)ctions we obt(cid:12)ain the following continued fraction repre-
sentation:
5
†
A ·A
0 0
R (ω) = . (12)
0 (cid:12) (cid:12)F
ω −E −(cid:12) (cid:12) 0
0 (cid:12) (cid:12) F
e 1
ω −E −
1 ...
Thus the recursive procedure (8) in course of which the coefficients E and F of the con-
n n
tinued fraction (12) are determined allows us to calculate Green’s or Kubo’s relaxation
functions.
IV. THE SPIN GREEN’S FUNCTION
The direct application of Eqs. (8) and (12) to the spin Green’s function D(kω), Eq. (5),
meets with difficulties because the inner product sz,sz in the numerator of the contin-
k −k
ued fraction (12) is equal to zero. To overcome thDihs difficuiltEy we consider Kubo’s relaxation
function sz sz defined in Eq. (11). In this case the inner product (6) in the numerator
k −k
of the res(cid:16)p(cid:16)ect(cid:12)ive c(cid:17)o(cid:17)ntinued fraction is nonzero. After calculating the relaxation function the
(cid:12)
spin Green’s (cid:12)function can be obtained from the relation
ω sz sz = sz sz + sz,sz , (13)
k −k k −k k −k
(cid:16)(cid:16) (cid:12) (cid:17)(cid:17) DD (cid:12) EE (cid:16) (cid:17)
where we dropped the subscrip(cid:12)t ω in the rel(cid:12)axation and Green’s functions.
(cid:12) (cid:12)
We postpone the calculation of the numerator sz,sz of the continued fraction and
k −k
consider its other coefficients. From definition (6) w(cid:16)e find t(cid:17)hat E sz,sz = is˙z,sz =
0 k −k k −k
sz,sz = 0 and therefore A is the Fourier transform of the op(cid:16)erator (cid:17) (cid:16) (cid:17)
k −k 1
Dh iE
1 1
is˙zl = 2 Jmn(δln −δlm)s+n1s−m1 + 2 tmnσ(δlm −δln)a†nσamσ = Asl +Ahl. (14)
mn mnσ
X X
Here the dot over the operator indicates the time derivative. As can be seen, A contains
1
contributions from spin and hole components As and Ah. Using this result in calculat-
†
ing R (ω) which is the Laplace transform of the function A ,A we neglect the terms
1 1t 1
Ah,As† and As ,Ah† . This approximation is motivate(cid:16)d by va(cid:17)nishing values of these
1t 1t
c(cid:16)orrelatio(cid:17)ns obta(cid:16)ined with(cid:17) the decoupling. Therefore
A ,A† ≈ Ah(t),Ah† + As ,As† , (15)
1t 1 1t
(cid:16) (cid:17) (cid:16) (cid:17) (cid:16) (cid:17)
where we have additionally neglected the difference between Ah and Ah(t) (again due to
1t
zero values of the respective decoupling). In accord with our estimation the influence of
terms connected with holes in sz,sz and As As† on the spin Green’s function is
k −k k k
small in comparison with the qu(cid:16)antity (cid:17) Ah A(cid:16)h(cid:16)† (cid:12)even(cid:17)(cid:17)for moderate hole concentrations.
k k (cid:12)
(cid:12)
Therefore in the forthcoming discussion(cid:16)w(cid:16)e n(cid:12)eglec(cid:17)t(cid:17)these terms in sz,sz and As As†
(cid:12) k −k k k
(cid:12)
and consider the time evolution of operators in these quantities as(cid:16)determi(cid:17)ned so(cid:16)l(cid:16)ely b(cid:12)y th(cid:17)(cid:17)e
(cid:12)
Heisenberg part of Hamiltonian (1). In this approximation the numerator of the cont(cid:12)inued
fraction representing As As† reads
k k
(cid:16)(cid:16) (cid:12) (cid:17)(cid:17)
(cid:12)
(cid:12)
6
Ask,Ask† = is˙zk,−is˙z−k = is˙zk,sz−k = 4JC1(γk −1), (16)
(cid:16) (cid:17) (cid:16) (cid:17) Dh iE
where γk = 14 aexp(ika), Cp = N1 kγkpCk and Ck = nexp[ik(n−m)]hs+n1s−m1i. For
E we get E is˙z,−is˙z = i2s¨z,−is˙z = is˙z,−is˙z = 0. Thus, breaking off the
1 1 Pk −k k P−k k −kP
continued frac(cid:16)tion on this(cid:17)step(cid:16)we obtain f(cid:17)rom DEhqs. (12), (1i3E), (15) and (16)
ω Ahk Ahk† +4JC1(γk −1)
D(kω) = , (17)
(cid:16)(cid:16) ω2(cid:12)−2ω(cid:17)(cid:17)Π(kω)−ω2
(cid:12) k
(cid:12)
where the polarization operator and the excitation frequency are given by
1 −1 −1
Π(kω) = 2 Ahk Ahk† szk,sz−k , ωk2 = 4JC1(γk −1) szk,sz−k . (18)
(cid:16)(cid:16) (cid:12) (cid:17)(cid:17)(cid:16) (cid:17) (cid:16) (cid:17)
(cid:12)
To calculate sz,sz (cid:12)in the above formulas we notice that in the considered case A =
k −k 2
i2s¨z − is˙z,−is˙z(cid:16) sz,s(cid:17)z −1sz and
k k −k k −k k
(cid:16) (cid:17)(cid:16) (cid:17)
A ,A† = i2s¨z,−is˙z − 16J2C12(γk −1)2 = 0. (19)
2 2 k −k
sz,sz
(cid:16) (cid:17) Dh iE k −k
(cid:16) (cid:17)
We set the above result equal to zero in conformity with the approximation made above in
the continued fraction where we dropped all terms containing A and operators of higher
2
orders. Equation (19) can be used for calculating sz,sz if the value of i2s¨z,−is˙z is
k −k k −k
known. An analogous equation was obtained in R(cid:16)ef. 20 w(cid:17)ith another methDohd. iE
We calculate i2s¨z,−is˙z in Eq. (19) by decoupling terms in the second derivative of
k −k
sz, Dh iE
1
is¨zl = JlmJln 2szls+n1s−m1 −szns+l 1s−m1 −s+n1szms−l 1 +
2
mn
Xh (cid:16) (cid:17)
JlmJmn szns+m1s−l 1 −szms+n1s−l 1 +s+l 1szns−m1 −s+l 1szms−n1 . (20)
(cid:16) (cid:17)i
In the decoupling we approximate szls+n1s−m1 by the value αCnm(1−δnm)+ 21δnm szl where
Cnm = hs+n1s−m1i. In the last expression we took intohaccount that in accordaince with
Eq. (3) Cnn = 1 for x = 0 [let us remind that we neglect the influence of holes on the value
2
of sz,sz ]. Following Ref. 13 the parameter α is introduced to fulfill the constraint of
k −k
zer(cid:16)o site m(cid:17)agnetization (3) in the paramagnetic state. Before carrying out the decoupling
it has to be taken into account that terms of Eq. (20) in which the site index of the sz
operator coincides with the site index of s+1 or s−1 operators cancel each other. To verify
this statement it is necessary to take into consideration that in these terms the operators sz
n
can be substituted by −1, since for the spin-1 case sz = −1 +s+1s−1 and s+1s+1 = 0. To
2 2 n 2 n n n n
retain this exact cancellation it has to be taken into account before the decoupling. As the
result we find
i2s¨zl = α JlmJln Cmnszl −Clmszn +JlnJmn Clnszm −Clmszn
mn
Xh (cid:16) (cid:17) (cid:16) (cid:17)i
+ Jl2n (1−α)Cnn szl −szn +αCln szn −szl ,
n
X h (cid:16) (cid:17) (cid:16) (cid:17)i
7
and after the Fourier transformation
i2s¨z = ω2sz,
k k k
where
i2s¨z,−is˙z
2 k −k 2
ωk = D4hJC1(γk −1i)E = 16J α|C1|(1−γk)(∆+1+γk),
(21)
C 1−α 3
2
∆ = + − .
|C | 8α|C | 4
1 1
Combining Eqs. (16), (18) and (21) we find
(szk,sz−k)−1 = 4Jα(∆+1+γk). (22)
In the absence of holes Eqs. (17) and (21) are close to the equations for the spin Green’s
function and the excitation frequency obtained for the two-dimensional Heisenberg antifer-
romagnet in Ref. 21,22 with the use of the equations of motion for Green’s functions and
Tserkovnikov’s formalism,20 respectively. In these works, somewhat more complicated de-
couplings were used. These decouplings contain several decoupling parameters of the type
of α which depend on the site indices in the decoupled average. These additional parameters
allow one to obtain somewhat better agreement with numeric simulations. However, to fix
the additional parameters exterior data from numerical simulations or the spin-wave theory
have to be engaged and the theory ceases to be closed.
As can be shown by the analogous calculation of the transversal spin Green’s function
s−1 s+1 , in the paramagnetic state
k k
DD (cid:12) EE
(cid:12)
(cid:12) s−1 s+1 = 2 sz sz . (23)
k k k −k
DD (cid:12) EE DD (cid:12) EE
Thus, the rotation symmetry of the(cid:12)components of(cid:12)the magnetic susceptibility is retained in
(cid:12) (cid:12)
this approach. This fact can be used for the calculation of parameters C , C and α in the
1 2
above formulas. From Eq. (17) simplified for the absence of holes, Eq. (23) and the relation
∞
hsz(t)sz i = dωe−iωteβωn (ω)B(kω), (24)
k −k B
−∞
Z
we find
1
Ck = 4J|C1|(1−γk)ωk−1coth βωk , (25)
2
(cid:18) (cid:19)
where B(kE) = −π−1ImD(kE) is the spin spectral function, n (E) = [exp(βE)−1]−1 and
B
β = T−1 is the inverse temperature. Substituting this equation in the definitions of C , C
1 2
and in constraint (3) we obtain three equations for the three unknown parameters C , C
1 2
and α. This problem can be reduced to the optimization problem and solved by the steepest
descent method.
Tocheckthevalidityoftheapproximationsmadeaboveweusedtheobtainedformulasfor
calculating spin correlations in an undoped antiferromagnet. In Fig. 1 our results obtained
8
ina32×32latticeforthreetemperatures arecomparedwithdataofMonte Carlosimulations
performed for the same lattice in Ref. 14. As can be seen, the agreement is good. However,
it should be noted that at elevated temperatures in our approximation the spin correlations
are systematically overestimated in comparison with the Monte Carlo results.
As follows from Eq. (21), for low temperatures and large crystals the spectrum of ele-
mentary spin excitations is close to the spectrum of spin waves.23 For an infinite crystal and
T = 0 we found α = 1.70494 and C = −C = 0.206734. In this case in Eq. (21) the param-
2 1
eter ∆ = 0 and the excitation frequency vanishes in the two points of the Brillouin zone,
k = (0,0) and (π,π) (here and below the intersite distance is taken as the unit of length).
For any nonzero temperature ∆ becomes finite which generates a gap at the (π,π) point.
It can be shown9,21 that the gap leads to the exponential decay of spin correlations with
distance and the respective correlation length is defined by the magnitude of the gap. Thus,
in agreement with the Mermin-Wagner theorem16 for a nonzero temperature the long-range
antiferromagnetic order is destroyed in the considered two-dimensional system.
Now let us calculate the polarization operator Π(kω), Eq. (18). Using the decoupling
which is equivalent to the Born approximation8 and the relations
∞ ∞
hakσ(t)a†kσi = dωe−iωteβωnF(ω)A(kω), ha†kσakσ(t)i = dωe−iωtnF(ω)A(kω), (26)
−∞ −∞
Z Z
we find
π ∞
ImΠ(kω) = f2 dω′[n (ω′)−n (ω′ −ω)]A(k′ −k,ω′−ω)A(k′ω′),
k′k F F
ω
k′ Z−∞
X
(27)
∞ dω′ImΠ(kω′)
ReΠ(kω) = P ,
π ω′ −ω
−∞
Z
where A(kω) = −π−1ImG(kω) is the hole spectral function, n (ω) = [exp(βω) + 1]−1,
F
fk′k = 2tN−1/2(γk′ − γk′−k)(szk,sz−k)−1/2 and P indicates Cauchy’s principal value of the
integral. We notice that Eq. (27) is close in its formto the polarizationoperator obtained for
thet-J modelinthespin-waveapproximation.8,10 Wecannotdirectlycomparetheinteraction
constants, because the definition of the spin Green’s function in this paper differs from the
magnon Green’s functions in Refs. 8,10. However, we notice that the spin-wave interaction
constant and the respective quantity f2 ω−1 in Eq. (27) are of the same order of magnitude
k′k k
andtendtozerolinearlywith|k|when|k| → 0. Thespin-wave constant behavesanalogously
near the(π,π) point, while the quantity f2 ω−1 does so only inthe case of aninfinite crystal
k′k k
and zero temperature.
V. THE HOLE GREEN’S FUNCTION
Now let us consider the hole Green’s function. To use the continued fraction repre-
sentation (12) for the anticommutator Green’s function G(kt), Eq. (5), the average of the
anticommutator of operators has to be taken as the definition of the inner product in the
recursive procedure (8). From the commutation relations (4) we find for the numerator of
the continued fraction akσ,a†kσ = 21(1+x) = φ and for the time derivative
Dn oE
9
1
ia˙lσ = tlm 1−s−l σsσl sσm +sσl am,−σ − Jlm(σszmsσl +sσm)al,−σ +µalσ. (28)
2
m m
X h(cid:16) (cid:17) i X
With these results we get
−1
E0 = ia˙kσ,a†kσ akσ,a†kσ = εk +µ′,
Dn oEDn oE (29)
εk = (4tφ+6tC1φ−1 −3JF1φ−1)γk, µ′ = µ+(4tF1 −3JC1)φ−1,
where F1 = N−1 kγkFk and Fk = nexp[ik(n−m)] a†nam .
The estimation of t and J based on the parametersDof thEe extended Hubbard model24
P P
gives J/t lying in the range 0.2−0.3. For low hole concentrations we can approximate the
parameter C by its value in an undoped lattice. ForT = 0.02t in a 4×4 lattice C = 0.2119,
1 1
while in a 20×20 lattice C = 0.2068. With these parameters the unrenormalized hole
1
dispersion can be estimated as εk ≈ −0.27tγk in the former case and −0.47tγk in the latter
case. Thus thefirst approximationoftherecursive proceduredescribes abandwhich ismuch
narrowerthanthetwo-dimensionalnearest-neighborbandintheabsenceofcorrelations4tγk.
The reason for this is the antiferromagnetic alignment of spins when the hole movement is
accompanied by the spin flipping. With increasing the hole concentration C → 0 and the
1
unrenormalized dispersion tends to its uncorrelated value.
The hole Green’s function reads
φ
G(kω) = ,
ω −εk −µ′ −Σ(kω)
(30)
Σ(kω) = φ−1 A1 A†1 , A1 = ia˙kσ −(εk +µ′)akσ,
DD (cid:12) EE
(cid:12)
where the difference between A1t and(cid:12) A1(t) was neglected. Due to the mentioned small-
ness of εk for low hole concentrations and of J in comparison with t, only the term
N−1/2 lmexp(−ikl)tlm 1−s−l σsσl sσm +sσl am,−σ may be retained in A1 in the calcula-
tion ofP A1 A†1 . The tehr(cid:16)ms in A1 w(cid:17)hich are ilinear in spin operators produce the following
contribuDDtion(cid:12)toEtEhe self-energy:
(cid:12)
(cid:12)
32t2 ∞ n (−ω )+n (ω )
dω1dω2 F 1 B 2 (γk +γk−k′)2A(k−k′,ω1)B(k′,ω2).
Nφ ω −ω −ω +iη
k′ ZZ−∞ 1 2
X
Up to the prefactor this expression coincides with the respective term in the hole self-energy
calculated in the spin-wave approximation.10
The term with three spin operators inA produces terms in the self-energy which contain
1
two- and three-spin Green’s functions. To calculate these functions one would have to solve
the respective self-energy equations which could be derived in the same way as the equations
in the previous section. However, such program would essentially complicate the calculation
procedure. One of the possible ways to overcome this difficulty is to use the decoupling in
the same manner as we applied it in the previous section, this time in the term with three
spin operators in A . However, the comparison with the exact diagonalization data shows
1
that this approximation does not give satisfactory results. Another way of simplification is
10