Journal of the Physical Society of Japan

FULL PAPERS

Band-Renormalization Eﬀects and Predominant Antiferromagnetic
Order in Two-Dimensional Hubbard Model
Ryo Sato ∗ and Hisatoshi Yokoyama
Department of Physics, Tohoku University, Sendai 980-8578, Japan
Band renormalization eﬀects (BRE) are comprehensively studied for a mixed state of d x2 −y2 -wave superconducting (dSC) and antiferromagnetic (AF) orders, in addition to simple d-SC, AF, and normal (paramagnetic) states, by applying avariational Monte Carlo method to a two-dimensional Hubbard (t-t′ -U) model. In a weakly correlated regime (U/t ≲ 6),
BRE are negligible for all the states studied. As previously shown, the eﬀective band of d-SC is greatly renormalized but
the modifications of physical quantities, including energy improvement, are negligible. In contrast, BRE on the AF state
considerably aﬀects various features of the system. Because the energy is markedly improved for t′ /t < 0, the AF state
occupies almost the whole underdoped regime in phase diagrams. A doped metallic AF state undergoes a type of Lifshitz
transition at t′ = tL′ ∼ −0.05t as t′ /t varies, irrespective of the values of U/t and δ (doping rate). Pocket Fermi surfaces
arise around (π, 0) [(π/2, π/2)] for t′ > tL′ [t′ < tL′ ], which corresponds to the electron-hole asymmetry observed in angleresolved photoemission spectroscopy (ARPES) spectra. The coexistent state of the two orders is possible basically fort′ > tL′ , because the existence of Fermi surfaces near (π, 0) is a requisite for the electron scattering of q = (π, π). Actually,
the coexistent state appears mainly for tL′ /t < t′ /t ≲ 0.2 in the mixed state. Nevertheless, the AF and coexisting states
become unstable toward phase separation for −0.05 ≲ t′ /t ≲ 0.2 but become stable at other values of t′ /t owing to theenergy reduction by the diagonal hopping of doped holes. We show that this instability does not directly correlate withthe strength of d-SC.

and variational Monte Carlo (VMC) methods4, 5, 13–20) are useful tools to quantitatively treat strong local correlations. One
also needs to consider the eﬀects of antiferromagnetism (AF)because it is crucial even for subject (C). The results regarding (A)-(C) of the above studies do not seem unified but are
rather scattered at first glance. Although inconsistencies exist among them, we feel that the main source of confusion
resides in insuﬃcient consideration of the diﬀerence in the
diagonal hopping term (t′ ). In most of the above studies, t′ /t(and t′′ /t) was set to specific values, say 0 and/or −0.3, but
we are apt to read the results associated with (A)-(C) withoutcare in the value of t′ /t. If we arrange the results by specifying
the value of t′ /t, they are often consistent beyond our expectation, as shown later in Table IV for some results obtained by
the VMC method. This also applies to many results of DMFT.
From this point of view, the results of recent studies with highaccuracy5–7) are consistent. In fact, a small number of studies have considered the diﬀerence in the features of (A)-(C)between the cases of t′ /t = 0 and other cases, although they
were not suﬃciently elaborate or analytic.4, 7, 16, 20)
To study (A)-(C) in an ordinary VMC framework, one has
to use a mixed state that represents the AF and SC orders simultaneously. The properties associated with (B) have beenstudied for the t-J-type13–17) and Hubbard17–20) models. In
addition, it is crucial to take account of the eﬀects of band
renormalization (BR) owing to strong correlations in the onebody part of the wave function. To date, band renormalizationeﬀects (BRE) have been introduced into d-SC states21–25) or
the d-SC part of mixed states.17, 19, 20, 26, 27) Because BRE weredisregarded in the AF part in these studies, an AF order does
not arise for t′ /t ∼ −0.3, or it vanishes rapidly upon doping for t′ /t ∼ 0. Such features are inconsistent with recent
research.5–7) Unexpectedly, BRE have not been introduced
into normal (paramagnetic) and AF states and the AF part of

1. Introduction
To clarify the physics of cuprate superconductors (SCs),1, 2)
we have to know the fundamental properties of the t-J and
Hubbard models on a square lattice with an extension in the
kinetic part (t-t′ and t-t′ -t′′ , etc.) as basic models.3) In this paper, we mainly focus on the following subjects in the Hubbard(t-t′ -U) model:(A) The primary subject is the ground-state phase diagramin the model-parameter space. Although a typical view to dateis that the antiferromagnetic (AF) order arising at half filling rapidly vanishes upon doping holes and the d x2 −y2 -wave
superconductivity (d-SC) appears,3, 4) in accordance with the
behavior of cuprates, in recent studies using advanced techniques it was argued that AF orders or inhomogeneous phases
prevail in wider ranges of δ (doping rate).5–7)(B) In phase diagrams of cuprates, the areas of superconducting (SC) and AF phases are in proximity. In the SC phase,appreciable AF correlation or short-range AF orders are observed,8) but the coexistence of two long-range orders has notbeen detected except for in multilayered systems. In theory,
it is still unclear in what parameter range the two long-rangeorders coexist and why they are coexisting or mutually exclusive.
(C) Another subject is whether or not homogeneous states
are stable against phase separation. Actually, signs of inhomogeneous electronic states or phase separation are often noticed
in cuprates such as a stripe structure of charge and spin and amosaic distribution of the gap magnitude. Theoretically, it isagain unclear as to the ranges of U/t, t′ /t, and δ and the causeof the state becoming unstable toward phase separation.
So far, these subjects have been addressed by many researchers with a variety of methods, in particular, dynamical mean field theories (DMFTs) with some extensions6, 7, 9–12)∗ satoryo@cmpt.phys.tohoku.ac.jp

1

J. Phys. Soc. Jpn.

FULL PAPERS
Table I. Elements modified by band renormalization in one-body part for
finite systems (indicated by circles). The two elements merge for L → ∞.

(a) H

(b) Ȍ

Modified elements{k}occ or Fermi surface
Direct modification of εk

t’

t1

t

t

t3
t2

t4

j

τ

{k}occ , σ

If {k}occ is determined in accordance with the bare band dispersion ε̃k in Eq. (2), ΦN is the exact ground state of H for
U = 0. When the interaction is introduced, εk will be modified
by its self-energy. In the framework of many-body variation
theory, εk should be optimized along with the other parameters so as to reduce the total energy E = ⟨H⟩/Ns (Ns : numberof sites). Note that in ΦN [Eq. (4)], εk does not explicitly appear but has the eﬀect of determining {k}occ or the Fermi surface (see Table I for comparison). Namely, the operation of
BR for ΦN is simply reduced to the choice of {k}occ . To obtainfull BRE, we need to find the {k}occ that yields the lowest E/tamong all the {k}occ , but the number of choices of {k}occ grows
exponentially [roughly as Ns /4C N/8 (N: number of electrons)]
as the system size grows. In this work, we optimize E/t withinthe {k}occ that are generated by a tight-binding form of εk withdiagonal transfer:
(
)
εN
(5)
k = −2t cos k x + cos ky − 4t1 cos k x cos ky ,2.1 Hubbard modelWith cuprate SCs in mind, we consider the Hubbard model(U ≥ 0) on a square lattice with diagonal hopping:(i, j),στ(3)where d j = n j↑ n j↓ , h j = (1 − n j↑ )(1 − n j↓ ), and τ runs over thenearest-neighbor sites of site j. As shown before,35, 36) the DH binding eﬀect included in PQ is crucial for properly treating
Mott physics. The projector P has three variational parameters, g, ζd , and ζh , which trigger BR in Φ.
We turn to the one-body part Φ, which is the main point
for BRE. We start with the normal (paramagnetic) state. Let{k}occ denote the set of k points occupied by electrons in Φwith εk ≤ εkF (or symbolically k ∈ kF ). Then, the one-bodynormal state we use (a Fermi sea) is written as∏ΦN =c†k,σ |0⟩.
(4)

After introducing the model in Sect. 2.1, in Sect. 2.2 wedescribe the setup of trial wave functions, which is the coreof variation theory. In Sect. 2.3, we comment on a methodof computing expectation values with the present wave functions.

=

Φmix
⃝
⃝

2.2 Trial wave functions
Because our interest here is to grasp the nature of BRE
rather than obtain accurate numerical values, we employ
forms of trial functions that capture the essence of physics
but are as simple as possible. As many-body trial states, weuse a Jastrow type, Ψ = PΦ, where P is a two-body correlation factor (projector) and Φ is a one-body (mean-field-type)wave function. We use a simple form of P common to all trial
states, P = PG PQ , where PG is the well-known onsite Gutz∏
willer projector PG = j [1 − (1 − g)n j↑ n j↓ ]32) and PQ is thenearest-neighbor doublon-holon (D-H) binding factor,4, 33, 34)


∏(
∏(
∏ 
)
)

1 − h j+τ − ζh h j
1 − d j+τ  ,PQ =1 − ζd d j

2. Formulation

Hkin + HU∑ (∑)−
ti j c†iσ c jσ + H.c. + U
n j↑ n j↓ ,

ΦAF
⃝
⃝

As we will see, the diagonal hopping term Ht′ plays a crucialrole in the present theme. We use t and the lattice spacing as
the units of energy and length, respectively.

mixed states,28) probably because optimization is technicallybothersome, as mentioned in Sect. 2.3 and the Appendices.
In this paper, we study ground-state properties of the Hubbard (t-t′ -U) model by applying a VMC method with BRE of
up to fifth-neighbor hopping to a mixed state Ψmix in addition
to normal (paramagnetic), pure d-SC, and pure AF states. InAFΨmix , we renormalize the energy dispersions εSCk and εk independently. This parametrization is a key to finding correctfeatures of a mixed state. The present results are quantitatively
consistent with those in recent research.5–7) As the merits ofthe present study, we stress the following points: (a) We systematically study the dependence on the model parameters, inparticular, t′ /t and δ. (b) We clarify the physics underlying theproperties of Ψmix (or the Hubbard model) by comparing various levels of wave functions. Through these merits, we will
acquire a better understanding of subjects (A)-(C).
This paper is organized as follows. In Sect. 2, we explain
the model and method used in this study. In Sect. 3, we discuss
the results of BRE on the d-SC state. In Sect. 4, the results of
BRE on the normal (or paramagnetic) state are presented. In
Sect. 5, we consider the BRE on an AF state, referring to the
Lifshitz transition arising at t′ /t ∼ −0.05. In Sect. 6, we study
BRE on a mixed state of d-SC and AF orders, and discuss
prerequisites for the appearance of d-SC. In Sect. 7, we recapitulate the main results and make additional comments. In
Appendices A and B, details of the calculations and analysesof the normal and AF states are described, respectively. The
preliminary results referred to in this paper were presented inthree preceding publications.29–31)

=

Φd
−
⃝

In this work, the hopping integral ti j is t for nearest neighbors(≥ 0), t′ for diagonal neighbors, and 0 otherwise (Hkin = Ht +Ht′ ) [Fig. 1(a)]. The bare energy dispersion becomes()ε̃k = −2t cos k x + cos ky − 4t′ cos k x cos ky .
(2)

Fig. 1. (a) Hopping processes in Hamiltonian [Eqs. (1) and (2)] and (b)those corresponding to band-adjusting parameters tη (η = 1–4) in trial wave
functions [Eqs. (13)-(17)]. In both figure, t is the unit.

H

ΦN
⃝
−

(1)

j

where c jσ annihilates an electron of spin σ at site j, n jσ =c†jσ c jσ , and (i, j) indicates the sum of pairs on sites i and j.
2

J. Phys. Soc. Jpn.

FULL PAPERSof ϕ(k) is also modified by {k}occ determined by εAF
k , as summarized in Table I.
A pure one-body AF state ΦAF is given by the ∆d → 0 limitof Φmix as∏ΦAF =a†k,σ |0⟩,(18)

where t1 is varied. This form of εk has often been used for dSC states in previous studies21–25) and also seems reasonable
as a first setting for ΦN . Details of optimizing ΨN = PΦN are
described in Appendix A. The ordered states Φd , ΦAF , andΦmix introduced below are reduced to ΦN in the limit of ∆AF
and/or ∆d → 0.
We move on to the mixed state of AF and d-SC orders of a
fixed electron number, Φmix . This state is written as a d-wave
BCS state composed of AF quasiparticles:18)

 N2
∑

†
†
Φmix =  ϕ(k) a a  |0⟩,(6){k}occ , σwhere the AF quasiparticles are given by Eqs. (9) and (10)and {k}occ is determined by εAF
k in Eq. (12). There are fivevariational parameters (tηAF , ∆AF ) in ΦAF . A pure d x2 −y2 -wave
singlet pairing (BCS) state of a fixed electron number38) is
given by the ∆AF → 0 limit of Φmix as
 N2


∑
†
†
Φd =  ϕ(k) c c  |0⟩,(19)k↑ −k↓kwith
ϕ(k) =

εSC
k −µ+

√

∆k
2
2
(εSC
k − µ) + ∆k.

k↑ −k↓(7)kwith ϕ(k) given by Eq. (7). There are six variational parameters (tηSC , ∆d , µ) in Φd .

Here, µ is a variational parameter, which is reduced to thechemical potential for U/t → 0, and a d x2 −y2 -wave gap is assumed as∆k = ∆d (cos k x − cos ky ),2.3 Variational Monte Carlo calculations
In general, it is impossible to accurately calculate variational expectation values of a many-body wave function ⟨O⟩,
with O being an operator, by analytical means. Instead, inmany cases, the expectation values can be accurately numerically estimated using VMC methods.39–42) Recently, many
parameters (up to more than 106 ) in ⟨H⟩ have been eﬃciently optimized by newly introduced algorithms.43) In thepresent cases, however, we cannot adopt ordinary optimization schemes using derivatives of energy because E({γ}) is
constant (ΨN ) or nearly constant (ΨAF and Ψmix ) as a functionof the band parameters (tηΛ ) in the parameter set {γ} and has irregularly distributed discontinuities. To address BR in ΨN , we
combine a VMC method with the extrapolation scheme described in Appendix A. For ΨAF and Ψmix , we repeat a primitive linear optimization method in this study until optimization becomes successful, although better ways are applicable.
Details are described in Appendix B. For Ψd , ordinary optimization algorithms are applicable unless ∆d approaches zero.
For ∆d ∼ 0, a diﬃculty similar to that for ΨAF manifests itself.
We calculate physical quantities using more than 2.5 × 105
samples. The accuracy of the total energy of 10−4 t is preserved, similarly to in previous studies. It is laborious to accurately converge ∆AF (or ∆d ) and the band parameters to specific values because there is redundancy among these parameters. However, this aﬀects the calculations of physical quantities only slightly in most cases.
We use systems of Ns = L × L sites with L = 10–18 underperiodic-antiperiodic boundary conditions. The closed-shellcondition is not satisfied because we allow {k}occ to be optimized automatically, although the total momentum is preserved at zero. In this paper, we often consider rough systemsize dependence for δ ∼ 0.08 using L = 10, 12, 14, 16, and18 with N = 92, 132, 180, 236, and 296 (δ = 0.08, 0.0833,
0.0816, 0.0781, and 0.0864), respectively.

(8)with ∆d being a d-wave pairing gap parameter. As the AF
quasiparticles in Eq. (6), we employ a form of an AF HartreeFock solution at half filling with t′ /t = 0:a†k,σ = αk c†k,σ + sgn(σ) βk c†k+Q,σ ,(9)a†k+Q,σ = −sgn(σ) βk c†k,σ + αk c†k+Q,σ ,(10)where Q is the AF nesting vector (π, π), sgn(σ) = 1 (−1) forσ =↑ (↓), and
v
u
t
εAF
1
αk (βk ) = √1 − (+) ( )2k
.
(11)
2
2
εAF
+
∆
k
AF
Here, ∆AF corresponds to the AF gap parameter in the senseof mean-field theory.
To introduce BRE into Φmix , we extend the band disperAF
sions εSCk in Eq. (7) and εk in Eq. (11) independently by
including tight-binding hopping terms up to three-step processes shown in Fig. 1(b),
Λ
Λ
Λ
Λ
εΛ
k = γk + ε1 (k) + ε2 (k) + ε3 (k) + ε4 (k).

(12)

with Λ = SC or AF and
γk = −2t(cos k x + cos ky ),

(13)

Λ
εΛ
1 (k) = −4t1 cos k x cos ky ,

(14)

Λ
εΛ
2 (k) = −2t2 (cos 2k x + cos 2ky ),

(15)

Λ
εΛ
3 (k) = −4t3 (cos 2k x cos ky + cos k x cos 2ky ),

(16)

Λ
εΛ
4 (k) = −2t4 (cos 3k x + cos 3ky ).

(17)Here, the eight band-adjusting parameters tηΛ /t (Λ = SC or
AF, η = 1–4) are independent of t′ /t in H and are optimized
along with the other variational parameters (g, ζd , ζh , ∆d , µ,∆AF ). Note that the k points used in Eqs. (9) and (10) belong37)
to {k}occ determined by εAF
As a result, if {k}occ ink (not γk ).
cludes k points outside the folded AF Brillouin zone, ϕ(k) forthe corresponding k in the sum in Eq. (6) is doubled, and ϕ(k)for k (< kF ) inside the AF Brillouin zone becomes null. InAFΦmix , εSCk and εk are explicitly renormalized, and the weight3.

BRE on Pure d-Wave Pairing StateIn this section, we discuss BRE on the pure d-wave pairing
state without an AF order, Ψd = PΦd . In Sect. 3.1, we confirm that there is a large BRE in Ψd , as found in previous studies.21–25) In Sect. 3.2, however, we show that the improvement
in energy is unexpectedly small. In Sect. 3.3, we also find that3

J. Phys. Soc. Jpn.

FULL PAPERS

0.1

Table II. Rough estimate of coeﬃcients in fitting function Eq. (20) for
U/t = 12 estimated from data for L = 10-14.

(a)

δ
α+
α−

0

t1 / t-0.1-0.3
d-wavet' / t = -0.3

Normal

-0.5
0.3

(b)
0.2

Uc / tt2 / t0.1
0L=10 12

-0.1
-0.2

14

-0.3
0

δ
0.0
0.04
~0.08
0.12
0.16

d-wavet' / t = -0.3
10

20

0.04
0.41
0.41

∼ 0.08
0.55
0.67

0.12
0.25
0.90

0.16
0.05
0.95

In Fig. 2(a), we show the optimized values of t1 /t as functions of U/t, fixing the model parameter at t′ /t = −0.3. For
U/t ≲ 5, t1 preserves the bare value t1 ∼ t′ irrespective of δ;no substantial BR exists. This is because ∆d is very small inthis range of U/t and the state is reduced to the normal state,
as in the cases without BRE.4, 36) ΨN also shows no substantial BR in this range of U/t as shown in Sect. 4. The relatively
large statistical fluctuation in the case of ∆d ∼ 0 stems fromthe same diﬃculty as in ΨN in optimizing the band parameters [see Appendix A]. On the other hand, at U = Uc ∼ 6.5t36)(Uc /t: Mott transition point in Ψd ), t1 /t abruptly increases, inparticular, t1 /t approaches 0 at δ = 0. As previously pointed
out,4, 21) this BR of Ψd occurs so that the quasi-Fermi surface overlaps with or approaches antinodal points [(π, 0), etc.],
where the van Hove singularity exists for |t′ /t| ≤ 0.5 and thed-wave gap becomes maximum. Furthermore, the elastic electron scattering of q = Q connects these points with opposite
signs of ∆d . Restoration of the nesting condition, which is theprincipal cause of BRE for the AF state, seems a subordinate
aspect for Ψd . As δ increases, t1 slowly approaches the valueof t′ for the same reason (see Fig. 3). In contrast to t1 /t, theoptimized t2 /t remains almost zero (the bare value) for all U/tand δ for t′ /t = −0.3, as shown in Fig. 2(b).
Next, we consider the t′ /t dependence of t1 /t and t2 /t for
U > Uc . We find that the optimized t1 /t is roughly fitted byseparate linear functions of t′ /t for the hole- and electrondoped cases:

-0.2

-0.4

0.0
0.14
0.14

30

U/t
Fig. 2. (Color online) Optimized band parameters (a) t1 and (b) t2 of dwave singlet pairing state (BR2) as functions of U/t for several doping rates.
In (a), the area where Ψd is reduced to ΨN is shown by an arrow labeled‘Normal’. In (b), the Mott transition point at half filling is indicated by a grayarrow.

t1 /t = α± (δ) × t′ /t,
U / t = 12

0.4

t' / t= 0.3 0.0

t1 / t㸫0.3 10L

0
-0.2
-0.4

d-wave (BR2)
0

0.1

′

where α+ (δ) [α− (δ)] is the coeﬃcient for t /t > 0 [t /t < 0]at a fixed δ. If tη (η ≥ 2) is ineﬀective (this is actually thecase, as reported shortly), BRE are nonexistent for α± = 1′
and, inversely, εSC
k is renormalized to the case of t = 0 for
α± = 0. The values of α± depend on U/t only slightly andare shown for U/t = 12 in Table II. Although the magnitudesof BR exhibit opposite tendencies between α+ and α− for δ ≳0.08, α± is always positive. As a result of this positiveness, theconvexity (t′ /t > 0) or concavity (t′ /t < 0) of the bare Fermi
surface near (π/2, π/2) is preserved in the renormalized quasiFermi surface of εSC
k . As a result, the locus of a hot spot —theintersection of a (quasi-) Fermi surface and the AF Brillouinzone boundary, where scattering of q = (π, π) takes place—44)
is near (π, 0) for t′ /t < 0 but approaches (π/2, π/2) to someextent for t′ /t > 0.24, 45, 46) As we will see in Sect. 6.2, the lociof hot spots can be used to give a condition that a coexistentstate arises.
In contrast to t1 /t, t2 /t is again found to be almost zerofor any t′ /t and δ. The eﬀect of t3 and t4 is considered in Ψdwith four band parameters t1 –t4 in εk [Eq. (12)] (BR4). The
behavior of t1 /t and t2 /t for BR4 is basically similar to thatfor BR2 mentioned above. We found that both the optimized
t3 and t4 have small positive values (t3 /t ≲ 0.11, t4 /t ≲ 0.095,
with the largest values at half filling) almost independent oft′ /t. These values decrease as δ increases and almost vanishfor δ ≳ 0.1. As we will see in Sect. 3.2, the eﬀects of t3 and t4

12
14

0.2

(20)
′

0.2

δ
Fig. 3. (Color online) Optimized values of t1 /t for the d-wave pairing state(BR2) plotted as functions of doping rate for three values of t′ /t in the regimeof a doped Mott insulator (U/t = 12).

the modification of relevant physical quantities is negligible.
3.1 Large BRE for (doped) Mott insulators
First, we attempted to optimize Ψd only with two band parameters t1 and t2 by putting t3 = t4 = 0 for simplicity. We
abbreviate this two-band-parameter optimization to BR2.
4

J. Phys. Soc. Jpn.

FULL PAPERS

0.004

(a)

L= 10 120.0015

δ

0.0
0.04
~0.08
0.12
0.16

0.002

∆E=Ed

㸫E (BR2) t'/t=0 -0.3 10Ldd-waveU / t = 12

12
14

0.001

∆E / t0.003∆E / t14

d-wavet' / t = -0.3

0.0005

0.001
0

0

0

0

10

(b)

∆E / t

0.006

20

U/t

Uc / t0.008∆E = Ed30

Fig. 5. (Color online) Energy improvement [Eq. (21)] owing to BRE forthe d-wave pairing state (BR2) plotted as functions of doping rate for twovalues of t′ /t in a strongly correlated regime (U/t = 12). Broad statisticalerrors are indicated by bars for some data points.

d

14

δ
0.0
0.04
~0.08
0.12
0.16

0.004

0.002

Table III. Examples of total energy per site shown for a specific case(t′ /t = −0.3, U/t = 12, L = 10) for comparison among four states with
diﬀerent BR levels and three doping rates. The brackets denote errors in thelast digits.

d-waveU / t = 12

State

0
-0.5

0

Normal0.5

t' / td-waveFig. 4. (Color online) Energy improvement [Eq. (21)] owing to BRE forthe d-wave pairing state (BR2) shown for some values of δ and L, (a) as
functions of correlation strength with t′ /t = −0.3 and (b) as functions of t′ /twith U/t = 12. The Mott transition point at half filling is indicated by a thickgray arrow in (a). In (b), plausible areas of t′ /t for hole-doped (t′ /t < 0) andelectron-doped (t′ /t > 0) cuprates are indicated with gray arrows.

AF
Mixed

3.2 Slight improvement in energy by BRE
Here and in some later sections, we consider the improvement in the total energy per site owing to BRE, representedas(Λ = d, N, or AF)

Conditionof εkδ = 0.0

E/t
0.04

0.08

no BR
BR
no BR
BR2
BR4
no BR
BR4
BR 4+4−0.1855(2)−0.2660(1)−0.3222(2)−0.3235(2)−0.3241(2)−0.1879(2)−0.35319(2)−0.3559(2)−0.3230(2)−0.3360(2)−0.3816(4)−0.3827(1)−0.3828(4)−0.3288(2)−0.4201(3)−0.4211(2)−0.4259(2)−0.4310(1)−0.4602(2)−0.4606(3)−0.4606(10)−0.4259(2)−0.4881(1)−0.4915(2)∆E/t, which again corresponds to the degree of BRE on t1 /tshown in Fig. 3.
We have shown that the energy is basically improved according to the degree of BRE on t1 /t for every model parameter. Nevertheless, what is important here is that the magnitudeof ∆E/t is unexpectedly small. The precision (statistical error) of the energy in the present VMC calculations for Ψd ison the order of 10−4 t as shown by bars in Fig. 5, while themaximum value of ∆E/t is only ∼ 10−3 t (only slightly larger
than the errors). In Table III, E/t for Ψd is compared among
the cases of without BR, BR2, and BR4 for typical model parameters. We also find that the diﬀerence between BR2 andBR4 is very small. Furthermore, the diﬀerence in Ψd is an order (two orders) of magnitude smaller than that in ΨN (ΨAF )for all the above values of δ. This diﬀerence is later illustrated
in Fig. 11.

on energy and other quantities are also slight.
To summarize, BRE on Ψd are large for U ≳ Uc , δ ∼ 0
and large |t′ /t|. If these conditions are satisfied, the eﬀective
band tends to the bare band of a square lattice (εk → γk or|t1 /t| → 0). This feature of BRE on the d-wave pairing state
has already been pointed out in previous studies.21–27)∆E = EΛ − EΛ (BR),0.2

δ㸫 E (BR)

L=10 120.1(21)

where Ed [Ed (BR)] is the energy of Ψd without [with] BRE;
∆E/t ≥ 0 holds except for statistical errors. In Fig. 4(a), theU/t dependence of ∆E/t is shown for some values of δ fort′ /t = −0.3. The regime of finite ∆E/t for U > Uc corresponds to that of the finite BR of t1 /t shown in Fig. 2(a).
As δ increases, both the magnitude of BR and ∆E/t decrease
and almost vanish in the overdoped regime (δ ≳ 0.15). Figure 4(b) shows the t′ /t dependence of ∆E/t for U/t = 12,
which mostly corresponds to the degree of BR of t1 /t given
by Eq. (20) with α± in Table II. The exception for t′ /t > 0 andlarge δ is caused by the vanishing of hot spots, which BRE
alone cannot control. Shown in Fig. 5 is the δ dependence of

3.3 Small modification of quantities by BRE
First, we consider a d-wave pairing correlation function,
⟨
⟩
1 ∑ ∑
′
Pd (r) =(−1)1−δ(τ,τ ) ∆†τ (Ri )∆τ′ (Ri + r) , (22)
Ns i τ,τ′ =x̂,ŷ
where x̂ (ŷ) denotes the lattice vector in the x (y) direction,δ(τ, τ′ ) is the Kronecker delta, and ∆†τ (Ri ) is the creation operator of a nearest-neighbor singlet pair at site Ri ,√(23)∆†τ (Ri ) = (c†i↑ c†i+τ↓ + c†i+τ↑ c†i↓ )/ 2.

5

J. Phys. Soc. Jpn.

FULL PAPERS

8

0.03

㸫

Pd∞

t' / t = 0.3

6㸫4

10 (BR2)
12 (BR2)
14 (BR2)

0.01

Uc / tδ=0
no BR BR L
10
12
14

S(Q) S0(Q)δ = 0.00.02

㸫d-wavet' / t = 0.3

L
10 (no BR PQ')
12 (no BR PQ')
14 (no BR PQ')
10 (no BR)
12 (no BR)
14 (no BR)(a)δ ~ 0.08
no BR BR L
10
12

2

10 (BR4)

Pd∞

0

4

6

U c/ t

U/t

Pd ∞

0.02

0.01

10

6

8

10

Fig. 7. (Color online) Spin structure factor at Q = (π, π) measured fromthe bare (U = 0) value S 0 (Q) (= 1) compared between BR2 and no-BR casesfor δ = 0 and ∼ 0.08 for t′ /t = −0.3 as functions of U/t. A few cases with
diﬀerent values of L are shown. The Mott transition point at half filling areindicated by a thick gray arrow. Near the horizontal axis, the areas where
the d-wave correlation function P∞
d becomes sizable are indicated by dashedarrows for δ = 0 (blue) and ∼ 0.08 (green). See Fig. 6 for δ = 0.

L
12 (no BR, PQ')
10 (no BR)
12 (no BR)δ ~ 0.08
U / t =12

10 (BR4)

-0.5

4

U/t

(b)

10 (BR2)
12 (BR2)

0

8

0

0.8
0

0.5

d-waveno BR PQ'
no BR
BR2
BR4

t' / t0.6
n(k)Fig. 6. (Color online) Behavior of the d-wave SC correlation functions inΨd compared between BR cases and no-BR cases (a) at half filling and t′ /t =−0.3 as functions of U/t and (b) for δ ∼ 0.08 and U/t = 12 as a function oft′ /t. On the horizontal axis in (a), the Mott transition point is indicated by athick gray arrow. The data for “no BR P′Q ” are adopted from Ref. 4, in whicha slightly diﬀerent D-H factor is used.

0.4

Normalno BR
BR㸫t' / t = 0.30
δ = 0.0278
U / t = 12
L = 12

0.2
Γ
(0,0)

If Pd (r) remains finite for |r| → ∞ (P∞d ), a d-wave oﬀdiagonal long-range order exists; P∞roughly
represents thedsquare of the SC gap. For Ψd , we estimate P∞
d in the same
way as discussed in Appendix C in Ref. 4. As an example, inFig. 6(a), we show P∞
d at half filling for some levels of BR(and P) for L = 10-14. As discussed in Ref. 4, P∞
d is negligible for small values of U/t. As U/t increases, P∞d abruptly
increases at U/t ∼ 5, exhibits a sharp peak near the Mott
transition point Uc /t ∼ 6.5, and vanishes in the Mott insulator regime U > Uc . Although the peak value of P∞d tends
to be slightly decreased by BRE, the behavior does not varyas a whole. For δ > 0, the area where P∞
d is sizable extends
to large values of U/t, but the modification of P∞
d by BRE
remains small (not shown). The modification of P∞
d is also
small when t′ /t is varied, as shown in Fig. 6(b). Furthermore,
the diﬀerence between BR2 and BR4 is negligible.
Next, we consider the q = Q element of the spin structure
factor1 ∑ iq·(Ri −R j ) ⟨ z z ⟩e
SiS j .
(24)
S (q) =
Ns i j

kF
X
(π,0)

k

M
(π,π)

Γ
(0,0)

Fig. 8. (Color online) Comparison of momentum distribution function ofd-wave pairing state (red and brown) among various levels of BR (and PQ )in the regime of doped Mott insulators (U/t = 12) along the path (0, 0) →(π, 0) → (π, π) → (0, 0). The Fermi surface in the nodal (Γ-M) direction isindicated by an arrow labeled kF . For comparison, we add n(k) for the normal
state (blue) with and without BRE as discussed in Sect. 4.

is necessary for an increase in P∞
d because the electron scattering of Q yields an attractive force for pairing. Anyway, themodification of S (Q) by BRE is also small and only quantitative even at half filling.
Finally, we discuss the momentum distribution function
1∑ †n(k) =⟨c ckσ ⟩.
(25)
2 σ kσIt seems that n(k) sensitively reflects the variation of the effective band εk , which is considerably renormalized depending on the case. Figure 8 depicts n(k) in such cases with red

The U/t dependence of S (Q) is shown in Fig. 7 for δ = 0 and∼ 0.08. As previous studies pointed out, an increase in S (Q)6

J. Phys. Soc. Jpn.

FULL PAPERS(without BRE) and brown (with BRE) symbols for Ψd andblue symbols for the normal state ΨN (see Sect. 4). In accordance with the above expectation, the locus of kF (discontinuity) near the X point in ΨN is shifted to a neighbor k pointby BRE. Nevertheless, in Ψd , the modification by BRE is verysmall for the gap behavior in the antinodal area (k ∼X) as well
as for the discontinuity in the nodal direction [k ∼ (π/2, π/2)],
despite the large BRE (t1 /t ∼ −0.05 for t′ /t = −0.3 in Fig. 3).
This is probably because the choice of {k}occ , which is controlled by εk in ΨN and ΨAF , is unnecessary in Ψd as shownin Table I.
In summary, the BR of εSC
k itself is large (t1 /t → 0) forU ≳ Uc , large |t′ /t|, and δ ∼ 0, as previous studies elucidated.
Notwithstanding, BRE on relevant quantities as well as on theenergy in Ψd are very small and insignificant compared withthose on the normal and AF states discussed below.

-0.41

(a)
A3

A2

E/t

A1
A0
A4

-0.42

Optim. BR
U / t = 7.0
-0.43
-0.36

(b)

E/t

A3

no BR
(t1 = t')

A2

-0.37

A1

A4

4. BRE on Normal (Paramagnetic) State

A0
-0.38

In this section, we discuss BRE on the normal or paramagnetic state (projected Fermi sea),∏ΨN = PΦN = P
c†k,σ |0⟩.
(26)

-0.4

Optim. BR
-0.3

-0.2

U / t = 7.5-0.1

0t' / t{kocc },σWe cannot apply ordinary optimization procedures to ΨN that
use the gradients of E/t with respect to band parameters because E/t for ΨN with a finite N is constant in a certain areaof the band-parameter space. Hence, we must resort to a different way of optimizing ΨN , which is described in Appendix
A. Here, we focus on the features of the optimized ΨN .
Before discussing BRE, we briefly review some aspects ofΨN without BRE (t1 = t′ ).4, 36) At half filling, a Mott transition occurs at Uc /t ∼ 8.5 for t′ /t = 0; Uc /t increases as |t′ /t|increases: Uc /t ∼ 11.2 for |t′ /t| = 0.3. Although the Mott
transition does not exist for δ > 0, the nature of ΨN markedly
changes at U ∼ Uc . For U ≳ Uc , ΨN does not behave as asimple metal but takes on a typical feature of Mott physics(D-H binding eﬀect), behaving as a doped Mott insulator.
Now, we consider BRE. Because BRE are ineﬃcient or
weak for t′ /t ∼ 0, similarly to the case of Ψd , we first consider the moderate case t′ /t = −0.3. We start with half-filledcases. Similarly to in Ψd , the energy reduction by BRE iszero or very small for U/t ≲ 6, even if the optimized t1(to be precise, the area including t1 ) somewhat shifts fromt′ (the area including t′ ). As shown in Fig. 9 for L = 14, theoptimized energy indicated by an arrow is given by {k3 }occ
(A3 = [−0.27, −0.19]) for U/t = 7.0, while it is given by{k0 }occ (A0 = [−0.125, 0.125]) for U/t = 7.5. Namely, theoptimized band parameter t1 rapidly varies from ∼ t′ (∈ A3 )to ∼ 0 (∈ A0 ) between U/t = 7.0 and 7.5 in this case, andthe nesting condition is restored. For U/t ≳ 7, the optimized{k}occ remains equal to {k0 }occ , or the optimized value of t1 /tremains approximately 0 (εNk = γk ). Also the renormalizedstate becomes identical to the normal state without BR of thesimple square lattice, whose behavior is reviewed above.36)
Owing to BRE, the Mott transition point for t′ /t = −0.3 shifts
from Uc /t ∼ 11.2 to ∼ 8.5. In fact, the optimal energy at
U/t = 8.5 for |t′ /t| ≲ 0.5 (L = 12) is given by {k0 }occ ; thus,
if BRE are introduced, the properties of the Mott transitions
and Mott insulators for |t′ /t| ≲ 0.5 are reduced to those for thesimple square-lattice case (t′ = 0) without BRE.
In Fig. 10(a), we show the energy reduction owing to BRE

Fig. 9. (Color online) Total energy of ΨN at half filling compared between
the cases of (a) U/t = 7.0 and (b) 7.5 as a function of t′ /t. Solid circles
indicate E/t without BRE, namely t1 = t′ . Aℓ (ℓ: integer) indicates the area oft′ /t corresponding to {kℓ }occ . The optimized energy owing to BRE is given by
the lowest value among all extrapolated lines. For t′ /t = −0.3, the optimized
value of E/t is indicated by an arrow in each panel. A detailed explanation of
the optimization is given in Appendix A.

[Eq. (21)] for t′ /t = −0.3 as a function of U/t. At half filling, as U/t increases, ∆E/t abruptly increases at U/t ∼ 7 owing to the reason mentioned above, approximately as ∆E/t =α exp(−βt/U) with α and β being positive constants. Then,∆E/t exhibits a peak at U/t ∼ 11, which corresponds to Uc /t
for the case without BRE, then slowly decreases (proportionally to t/U for U/t → ∞). As the doping rate increases from
δ = 0, the overall behavior of the U/t dependence is preserved but the magnitude rapidly decreases. For all the doping
rates shown, ∆E/t is negligible for a weakly correlated regime(U/t ≲ 7), meaning that appreciable BRE are also a characteristic of strong correlation for the normal state. In Fig. 10(b),the t′ /t dependence (t′ /t < 0) of ∆E/t is shown in the regimeof Mott physics (U/t = 12) for some doping rates. The BRE
are largest at |t′ /t| ∼ 0.3-0.4 and small for |t′ /t| ∼ 0.
We turn to the doping dependence of ∆E/t. Shown in
Fig. 11 is ∆E/t for U/t = 12 and t′ /t = −0.3; these values
are marked with vertical gray lines in Fig. 10. As δ increases,∆E/t rapidly decreases as∆E/t ∝ α exp(−δ/δN ),(27)

with δN ∼ 0.022 (α : positive constant) in this case, as shownwith a thick dash-dotted line in Fig. 11. Thus, BRE almost
completely vanishes for δ ≳ 0.1. However, it should be emphasized that ∆E for ΨN is much larger than that for Ψd nearhalf filling, as shown in Fig. 11, and the BRE on ΨN are never
negligible.
Finally, we analyze ∆E by dividing it into three compo7

J. Phys. Soc. Jpn.

FULL PAPERS

0.1
L=14

(a)

∆E / t0.08

12Normalt' / t = -0.3

10 δ
0.0
0.02780.04~0.08

nents: ∆E = ∆Et + ∆Et′ + ∆EU . We find that ∆Et is positiveand becomes large for a large |t′ /t|, while ∆EU is negative andits magnitude is relatively small. Namely, the eﬀective band istransformed so as to gain kinetic energy Et at the cost of theinteraction energy EU . This corresponds to a general tendency
for a state in a strongly correlated regime to undergo a transition to reduce the kinetic energy.35, 36) This feature applies
to ΨAF and Ψmix . For ∆Et′ , the magnitude is small as compared with those of the other two components, except whent′ /t ∼ −0.1.

0.06
0.04
0.02

5. BRE on Pure Antiferromagnetic State
0
0.1

10

U/t

Uc/t

20

∆E / t0.08
0.06 L=16 14

12

10

30

In this section, we consider the features of BRE on the AF
state without a SC order,

Normal (b)ΨAF = PΦAF .

U / t = 12

In Sect. 5.1, we discuss the optimized parameters. In Sect. 5.2,
a large improvement in energy due to BRE is revealed. In
Sect. 5.3, topics associated with the Lifshitz transition are
considered. Details of the optimization of ΨAF are given in
Appendix B.

δ

0.0
0.0278
0.04
~0.08

0.04

5.1 Optimized band parameters
We start by clarifying the features of the optimized band parameters in ΨAF , for which we always use tη (η = 1–4) (BR4).
In Fig. 12, the U/t dependence of the optimized value of tη isshown for t′ /t = −0.3. For a small U/t (< UAF /t ∼ 2.75–3.5
for t′ /t = −0.3), no AF order exists and ΨAF is reduced toΨN . At U = UAF (AF transition point), ∆AF and the sublatticemagnetization (AF order parameter)2 ∑ iQ·r j zm=e⟨S j ⟩
(29)
Ns j

0.02
0
-0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0

t' / tFig. 10. (Color online) Energy improvement [Eq. (21)] owing to BRE forthe normal state shown for some doping rates and values of L, (a) as functionsof correlation strength with t′ /t = −0.3, and (b) as a function of t′ /t (< 0) with
U/t = 12. In (a), the Mott transition points at half filling are indicated by thickarrows (brown for BR case, gray for no-BR case). Guide lines proportional
to t/U are added (dashed lines). In (b), data for each δ are well fitted by∆E/t = −(α/x) exp(β/x) + γx with x = t′ /t and α, β, and γ being positive
constants, as shown with gray dashed lines.

㸫∆E = EΛ EΛ(BR)0.15∆E / t( Λ= d, N, AF )

AF

0.1

AF

suddenly become finite (not shown), probably as a result of afirst-order transition, for t′ /t , 0. For U > UAF , marked BRE
appears and tη becomes almost constant as a function of U/t.
Although tη is almost invariant as a function of U/t, it varies
with δ to some extent, at least for t′ /t = −0.3, as shown in
Fig. 12. In fact, this feature depends on t′ /t, as described in
the next paragraph. Anyway, we find that εAF
k is renormalized
so as to restore the nesting condition, irrespective of δ.
Shown in Fig. 13 is the t′ /t dependence of the optimizedtη /t for U/t = 12. At half filling [(a)], the renormalized valuesof tη /t and the other variational parameters (not shown) areconstant with respect to t′ /t. The optimized AF state is independent of t′ /t; this feature is common to all values of U/t(> UAF /t). In contrast, for δ > 0 [(b) and (c)], tη /t discontinuously changes at t′ = tL′ ∼ −0.05t, and the other parameters(not shown) also exhibit singular behaviors (a cusp or discontinuity) there. Checking various cases, we find that the valueof tL′ /t slightly varies as the model parameters (δ, L) vary but
is necessarily situated in the range −0.1 < tL′ /t < 0. Thus, in
doped cases, the AF phase is divided into two subphases according to whether t′ > tL′ [type (i)] or t′ < tL′ [type (ii)]. In
each subphase, the t′ /t dependence of tη′ /t is weak. However,the δ dependence of tη is weak in the type-(i) AF, whereas tη′changes markedly as δ increases in the type-(ii) AF. Thus, theeﬀective band εAF
k is distinct for the two subphases. As wewill discuss in Sects. 5.3 and 6, this transition is regarded as a
Lifshitz transition in the AF phase.

U / t = 12
t' / t = -0.3
N d-wave L
10
12
14
16

0.05
Normal(paramag.)0d-wave

0.1

(28)

0.2

δ
Fig. 11. (Color online) Energy improvement [Eq. (21)] by BRE comparedamong the normal, AF (Sect. 5), and d-wave pairing (Sect. 3) states as functions of doping rate. Data for some system sizes are plotted. The thick pink
dash-dotted line is a curve fitted using Eq. (27) with ΨN for all values of L.

8

J. Phys. Soc. Jpn.

FULL PAPERS

0.4

(a)
0.4

0.2

(a)

0
L=10

tη / ttη / t0.3

12
t1 / tt2 / tt3 / tt4 / t-0.2 UAF/tAF
t' / t = -0.30
δ = 0.0

-0.4
0.4

δ

(b)

0.2

AF
U / t = 12

δ = 0.0

0.1
0

-0.1
0.4

㹼0.08

(b)

0.2

tη / ttη / t

0.3

0

δ

0.2

㹼 0.080.1
00.4

tη / t

(c)
-0.1
0.4

δ = 0.16
L = 10

0.2

L=10

0

10

20

tη / t

0
30

12

14
t1
t2
t3
t4

0.3
0.2

δ = 0.12

0.1

U/t
0

(c)

Fig. 12. (Color online) Optimized band parameters in pure AF state fort′ /t = −0.3 plotted as functions of correlation strength. The doping rate is
diﬀerent among the three panels. The legends displayed in (a) are common to(b) and (c). The arrow in each panel indicates the AF transition point.

-0.1

-0.5

0

0.5

t' / tFig. 13. (Color online) Optimized band parameters in pure AF state for
U/t = 12 plotted as functions of t′ /t. The doping rate is diﬀerent among thethree panels. In (a), (b), and (c), data for L = 10–14, 10 and 12, and 10 are
plotted respectively. The legends displayed in (c) are common to (a) and (b).
The arrows in (b) and (c) indicate the values of tL′ /t.

AF
Finally, let us compare the optimized εSCk and εk . t1 /t,
which is the sole eﬀective band parameter in d-SC, behaves
as a linear function of t′ /t with a positive coeﬃcient [Eq. (20)and Table II], indicating that BRE are mild, and a trace of thebare band remains [see Fig. 31(a) later]. On the other hand, for
the AF part, BRE are prominent in that tη /t is almost independent of t′ /t, and for t′ < tL′ , the sign of t1 /t becomes oppositethat of t′ /t [Figs. 13(b) and 13(c)]. Generally, the optimized
AF
forms of εSCk and εk are distinct, particularly, in the case oft′ < tL′ . This applies to the mixed state.

In Fig. 14, we show the U/t dependence of the energy reduction by BRE [Eq. (21)]. In contrast to the case of d-wave
pairing, a very large energy improvement is brought about by
BRE for U > UAF , consistent with the large BR shown in
Fig. 12 for t′ /t = −0.3. Energy improvement occurs even fort′ /t = 0 [Fig. 14(b)] because the BRE on t3 and t4 are notsmall, as shown in Fig. 13, although ∆E/t is an order of magnitude smaller than that for t′ /t = −0.3. The U/t dependenceof ∆E/t for t′ /t = 0.3 (not shown) is quantitatively similar(somewhat smaller for a large δ) to the case of t′ /t = −0.3.
We find from Fig. 14(a) that ∆E/t is a monotonically decreasing function of δ for a large |t′ /t| (actually, when t′ /t ≲ −0.2
and t′ /t ≳ 0.3). This has already been illustrated in Fig. 11 for
U/t = 12 and t′ /t = −0.3. ∆E/t decreases as δ increases, but
the area of finite ∆E/t is considerably extended up to δ ∼ 0.22
for this parameter set. We repeat that the energy reduction by
BRE in ΨAF is much larger than that in the normal and dwave states. Such an improvement occurs in wide ranges of

5.2 Large energy reduction by BRE
As shown in a previous VMC study without BRE,4) the energy of the AF state is not lowered relative to the paramagnetic state even at half filling for |t′ /t| > 0.35–0.41 (dependingon L) and U/t = 12 (see Fig. 15 later). For δ > 0 and t′ /t < 0,
this boundary value of |t′ /t| tends to decrease; for example, fort′ /t = −0.3, the optimized AF gap ∆AF substantially vanishesfor δ ≳ 0.03. However, as discussed in preceding reports,29, 30)
the AF state ΨAF with BRE [Eq. (28)] is stabilized with respect to ΨN up to δ ∼ 0.16 (0.21) for t′ /t = 0 (±0.3). First, welook at this great improvement by BRE more systematically.
9J. Phys. Soc. Jpn.

FULL PAPERSN N(BR) SF(a)t' / t = -0.3

L=12 10-0.2

EAF - EAF(BR)0.1

d(BR) AF AF(BR) L
10
12
14

δ
0.0
0.04
~0.08
0.12
0.16
0.20

E/t

∆E / t0.2

dN

δ = 0.0
U / t = 12

AF

-0.3

0
0.015

(b)

t' / t = 0.0

L=12∆E / t

0.01

-0.5

10 δ
0.0
0.04
~0.08
0.12
0.16

10

20

0.5

t' / tFig. 15. (Color online) Total energy at half filling for U/t = 12 comparedamong various states with or without BR as a function of t′ /t. Dark (pale)
symbols indicate cases with (without) BRE. Open, filled, and half-filled symbols indicate the data of L = 10, 12, and 14, respectively. Arrows denote theenergy reductions brought about by BRE for the diﬀerent states. ‘SF’ indicates a staggered flux state.

0.005

0

0

d

30

U/t

|t′ /t| increases, mainly owing to the decrease in Et′ . Consequently, the order in Eq. (31) does not change over a wide
range of δ except for is the SF state, where it rapidly becomes
unstable, especially for t′ /t > 0.47) Incidentally, by analyzing
the charge-density structure factor, we find that ΨAF becomes
metallic for δ < 0.29) At any rate, ΨAF with BRE has much
lower energy than Ψd in the whole range of t′ /t in Fig. 16.
This is not the case for ΨAF without BRE (see also Table IV).
By drawing similar figures for various values of U/t, t′ /t,and δ, we construct the phase diagram in the U/t-δ space
shown in Fig. 17. It is notable that, in contrast to previous
studies, the AF area for t′ /t = −0.3 becomes wider than thosefor t′ /t = 0 and 0.3 and covers a very wide range of model
parameters U/t, t′ /t, and δ.

Fig. 14. (Color online) Energy gain brought about by BRE in ΨAF plotted
as a function of U/t for (a) t′ /t = −0.3 and (b) t′ /t = 0. Note that the scaleon the vertical axis in (a) is 10 times larger than that in (b). Data for several
doping rates are shown.

U/t (> UAF /t) and t′ /t (≲ −0.15).
Now, we compare the total energy among various states in
the regime of Mott physics (U > Uc ). In Fig. 15, we compare the t′ /t dependence of E/t at half filling among ΨN , Ψd ,and ΨAF . For each, the values with BRE and without BRE
are plotted. Without BRE, the total energy depends on t′ /t,whereas if BRE are introduced, ΨΛ (t′ /t) (Λ =N, d, AF) is
optimized at ΨΛ (0) for any t′ /t. Consequently, E/t becomes
independent of t′ /t because the diagonal hopping energy vanishes:Et′ ≡ ⟨Ht′ ⟩ = 0(for δ = 0),

5.3 Lifshitz transition and electron-hole asymmetry
Before discussing the Lifshitz transition, we mention the
behavior of the staggered magnetization [Eq. (29)] in ΨAF .
We find that m gradually increases as U/t increases for UAF <U ≲ 12t and is almost constant for U ≳ 12t, irrespective of δand t′ /t (not shown). Shown in Fig. 18 is the t′ /t dependenceof m for some values of δ and U/t. At half filling, m is constantand ∼ 0.88 (m becomes 1 for the Néel state) because ΨAF isinvariant for t′ /t, as mentioned in Sect. 5.1. For δ > 0, ananomaly appears at t′ = tL′ , and the diﬀerence in the two areas
becomes more conspicuous as δ increases.
To confirm that the transition arising at tL′ /t is a type of Lifshitz transition, we plot in Fig. 19 the momentum distributionfunction [Eq. (25)] in ΨAF for U/t = 12 along the path in theoriginal Brillouin zone mentioned in the caption. In panel (a),n(k) at half filling is drawn for L = 10–16, which is smoothalong the whole path, indicating that the state is insulating.
The system-size dependence is very small. On the other hand,in doped cases with δ ∼ 0.08 shown in panels (b) [type (i)] and(c) [type (ii)], pocket Fermi surfaces appear and the state be-(30)

and Et and EU become constant with respect to t′ /t. As a result, the energy in ΨAF (and ΨN ) is greatly reduced for large
values of |t′ /t|. The order of the energy becomes
EAF < Ed < ESF < EN(31)

for a wide range of |t′ /t| (at least < 0.7) at a fixed U/t (> Uc /t).
Here, ‘SF’ indicates a staggered flux state, which is a candidate pseudogap state in cuprates3, 47) and will be discussed in
Sect. 6.3.
To consider doped cases (δ > 0), E/t for various states arecompared in Fig. 16 for δ = 0.0833 (L = 12). Similar resultsfor other values of δ were presented in Fig. 2 in a preceding report.31) The energy reduction in ΨAF brought about by
BRE for large |t′ /t| is still sizable, and E/t exhibits diﬀerent
linear behaviors on opposite sides of the Lifshitz transition
point tL′ /t. In ΨN and Ψd , E/t tends to decrease for t′ /t > 0
as t′ /t increases, and also decreases for t′ /t ≲ −0.4–−0.5 as
10J. Phys. Soc. Jpn.

FULL PAPERS

t1 / t
-0.5

-0.3

0

0.5

1

Normal (BR)E(t1) / t
( t' / t = 0.3 )

δ=0

0.6

δ ~ 0.08

U / t = 12

10
12 L
14

m

㸫

0.8

0.4

-0.4

L= 10

E/t

0.2

12

δ = 0.16

14 U/t
6.0
12

U / t = 12
L = 10t'L / t0

δ = 0.0833
U / t = 12
L = 12-0.5
type (ii)

AF (BR4)

0t' / t0.5
type (i)Fig. 18. (Color online) Sublattice magnetization [Eq. (29)] in the AF phaseplotted as a function of t′ /t for various values of δ and U/t. The red arrow
indicates the Lifshitz transition point.

-0.5
Normal (No BR)
Normal (BR)Normal (BR, optim.)
Staggered flux
d-wave (No BR)
d-wave (BR2)
AF (No BR)
AF (BR4)

t'L / t-0.5

0t' / t = -0.3

comes metallic.48) In each panel, we plot data for various values of t′ /t (L = 12) and for various system sizes for a typicalt′ /t (0 or −0.3) at the same time. In the type (i) [(ii)] regime, apocket Fermi surface arises around the antinodal point (π, 0)[around (π/2, π/2) in the nodal direction]. To visualize thisfeature, we constructed corresponding contour maps of n(k)as shown in Fig. 20. The location of the pocket Fermi surface
suddenly jumps from ∼ (π, 0) to ∼ (π/2, π/2) at t′ = tL′ , although the behavior of n(k) other than at the Fermi surfacechanges only slightly with tL′ /t. Note that the form of thepocket is almost preserved for a fixed δ as t′ /t is varied. It
is notable that the pocket is narrow but very deep, suggestingthat the advantages of half filling, such as the nesting condition, are well preserved by filling this narrow pocket with
doped carriers and leaving the other parts intact. Anyway, thisfirst-order transition occurs with a topological change in the
Fermi surface.
The source of this topological transition may have alreadyarisen in the bare tight-binding dispersion or at the mean-fieldlevel. In Fig. 21(a), we show the Fermi surface at half filling for t′ = 0, namely, the AF Brillouin zone boundary, on
which ε̃k = γk = 0 as shown in Fig. 21(b) in red. If we add
an infinitesimal diagonal hopping term [ε1 (k)] (blue), the degeneracy in the region (π, 0)–(π/2, π/2) is lifted and the band
maximum appears at (π, 0) or (π/2, π/2) according to whether
t′ /t > 0 or t′ /t < 0. As shown in green in Fig. 21(b), the thirdneighbor hopping term ε2 (k) has a similar eﬀect, if the sign oft′′ is opposite the sign of t′ , although we do not treat it here.
If we consider ordinary AF mean-field theory, the situation issimilar because the quasi-particle dispersion√
U(32)− γk2 + ∆2AF
EkAF =
2is degenerate in the region (π, 0)–(π/2, π/2). When we add
Ht′ as a perturbation to this framework, the leading diﬀerencein the dispersion relation is again ⟨ΦAF |Ht′ |ΦAF ⟩ ∝ ε1 (k). In
these examples, the boundary of the topological change is at

0.5

t' / tFig. 16. (Color online) Comparison of total energy among normal, AF, dwave, and staggered flux47) states with various levels of BR as a functionof t′ /t for L = 12, δ = 0.0833, and U/t = 12. The blue arrow indicatesthe Lifshitz transition point of the AF state. In addition, we illustrate theprocedure for obtaining the optimized energy with BRE for the normal state(red line) from raw data without BRE (black circles) for t′ /t = −0.3. The
green line denotes the variational energy for BRE as a function of t1 /t (upperaxis) for t′ /t = −0.3, which corresponds to the dark-green line in Fig. A·2.
For details, see Appendix A.

U/t

20

t' / t㸫0.3
0.0
0.3
L = 10, 12

d-wave

10

AF

0

0.1

Almost normal

δ

0.2

0.3

Fig. 17. (Color online) Rough phase diagram in U/t-δ space constructedwithin ΨAF and Ψd , both with BRE. Because the dashed border lines indicate
the locus of vanishing AF orders, the region of the d-wave state may somewhat extend to the AF side. The region of the d-wave is schematic, especially,on the large-δ side.

11

J. Phys. Soc. Jpn.

FULL PAPERS

1(a)t' / t = 0(b)AF (BR4)
L
16
14
12
10

n(k)

0.6

0.4

㸫0.3L

t' / t =(c)

L = 12t' / t
0.6
0.5
0.4
0.3
0.2
0.1

㸫
㸫
㸫
㸫
㸫
㸫

18
16
14
12
10

δ ~ 0.08

t' / t : arbitrary
δ = 0.0
0.2

L = 12t' / t
0.0
0.1
0.2
0.3
0.4
0.5
0.6

L
18
16
14
12
10

0.8

δ ~ 0.08

U / t =12

Insulating
0
Γ
(0,0)

X
(π,0)

Type (i)Γ(0,0)

M(π,π )

k

X
(π,0)

k

Type (ii)
Γ
(0,0)

M
(π,π)

X
(π,0)

Γ
(0,0)

M
(π,π)

k

Fig. 19. (Color online) Momentum distribution function plotted along the path Γ → X → M → Γ for U/t = 12 in three cases: (a) half-filled case, in whichΨAF becomes independent of t′ /t, (b) doped case (δ ∼ 0.08) with t′ > tL′ [type (i)], and (c) doped case (δ ∼ 0.08) with t′ < tL′ [type (ii)]. Data for L = 10-18
are plotted together. Pocket Fermi surfaces for the doped cases are indicated by thick arrows.

AF(BR) (a) t′/t : arbitrary, δ=0(b) t′/t=0, δ ≈ 0.08(c) t′/t=-0.3, δ ≈ 0.08

0.25
0.25

0.70
0.75

-3

-2

-1

0

1

0.20

0.75
0.70

2

3

-3

-2

-1

kx

0

1

kx

0.75

ky

0.70

2

3

-3

-2

-1

0

1

2

3

kx

Fig. 20. (Color online) Contour maps of momentum distribution function n(k) of the optimized pure AF state for U/t = 12 shown for (a) δ = 0 with
arbitrary t′ /t and for δ ∼ 0.08 with (b) t′ /t = 0 [type (i)] and (c) −0.3 [type (ii)]. The parameters in (a), (b), and (c) correspond to those in (a), (b), and(c) in Fig. 19, respectively. The maps are constructed using data for L = 10–18. In these contour maps (and similar ones displayed henceforth), the fourfold
rotational symmetry is somewhat smeared on account of anisotropic boundary conditions, open-shell conditions, and functions of the graphic software used.

(a)t′ /t = 0. Nevertheless, it is not trivial whether this topologicalchange is connected to those in strongly correlated cases (andeven with a large δ), and, if it is connected, why tL′ /t slightly
deviates to the negative side of t′ /t.
A topological change equivalent to the present result was
found in the spectral function A(k, ω) for the cases in which
a few carriers are doped in the t-t′ -J model and its extensions
using various methods.49, 50) In particular, Refs. 51 and 52
clearly argued, by means of a self-consistent Born approximation and a VMC method, respectively, that the location of theband maximum is diﬀerent between hole- and electron-doped
cases for typical parameters of cuprates. Actually, angleresolved photoemission spectroscopy (ARPES) experiments
revealed that the evolution of the Fermi surface with doping is diﬀerent in hole-doped53) and electron-doped54) casesin lightly doped systems, in accordance with the results of theabove theoretical study. Our result for the Hubbard model directly corresponds to these results for slightly doped t-J-typemodels.
As we will discuss in Sects. 6.2 and 6.3, this topological

(b)
(0,ʌ)

(ʌ,ʌ)

Antinodal

4

Nodal

ε1 /t'(-ʌ,0)(ʌ,0)ī

εk / t

2
0

γk /t

-2

ε2 /t''-4(0,-ʌ)(π,0)

k(π/2,π/2)

Fig. 21. (Color online) (a) Bare Fermi surface at half filling in the tightbinding model with t′ /t = 0 shown with a pink dashed line in the first
Brillouin zone. The nodal (π/2, π/2) and antinodal (π, 0) areas are marked
by shadows. (b) Elements of bare band dispersion relations along (π, 0) →(π/2, π/2): γk /t = −2(cos k x + cos ky ), ε1 /t′ = −4 cos k x cos ky , and ε2 /t′′ =−2(cos 2k x + cos 2ky ). t′′ indicates the hopping integral to the third-neighborsites (±2, 0) and (0, ±2), which is disregarded in this paper.

12

J. Phys. Soc. Jpn.

FULL PAPERS

Table IV. Relative and intrinsic stabilities of pure AF states and mixed
states of AF and d-SC orders obtained in recent studies using the Hubbard
model summarized according to the level of BR and to whether t′ /t ∼ 0 or
−0.3. In the U/t column, a typical target value is given. The first row denotes
the range of finite AF orders. The second row indicates whether the system ishomogeneous or phase-separated (P. S.). The third row for the mixed statesshows whether AF and d-SC orders are coexisting or mutually exclusive (andthe dominant order) in the main (or small-δ) area of m > 0.
Trial statesU/tt′ /t ∼ 0t′ /t ∼ −0.3

References

AF (no BR)

8, 12
12

no AF
—
δ ≲ 0.22
homogeneous

4

AF (BR)δ ≲ 0.15
P. S.
δ ≲ 0.16
P. S.
δ ≲ 0.2—
coexistingδ ≲ 0.15
P. S.
coexistingδ ≲ 0.16
P. S.
coexistingδ ≲ 0.18
P. S.
coexisting———δ ≲ 0.15
P. S.
exclusive, AFδ ≲ 0.25
homogeneous
exclusive, AFδ ≲ 0.24
homogeneous
exclusive, AF

Mixed (no BR)
10
Mixed (BR onlyin SC)

10

Mixed (BR in
AF & SC)

12

Mixed (manyparameters)

10

-0.3

2
t′/t L

-0.4
-0.5

E/t

10
12
10
12
10
12
14
10
12

29, 30 &this work

1

m

-0.7

18

1.5

m

E/t

-0.6

-0.3
-0.3
-0.1
-0.1
0.0
0.0
0.0
0.3
0.3

-0.8
0.5

20, 23-0.9

31 &this work

-1

5

0

0.05

0.1

δ

0.15

0.2

0
0.25

Fig. 22. (Color online) Total energy per site and staggered magnetization(right axis) obtained for Ψmix with U/t = 12 plotted as a function of doping
rate. Data for four values of t′ /t and L = 10-14 are shown.

diﬀerence in the Fermi surfaces in ΨAF determines whether
or not the d-wave SC order coexists with AF orders.
Table V. Second-order coeﬃcient c2 estimated by the least-squares
method for E(δ)/t [Eq. (34)] in the AF phase (U/t = 12) of Ψmix . For positive (negative) c2 , Ψmix is stable against (unstable toward) phase separation.
Digits in round brackets indicate the error in the last digit.

6. BRE on Mixed State of AF and SC Orders
In this section, we study a mixed state of AF and d-SC orders in a strongly correlated regime (U > Uc ):Ψmix = PΦmix ,

t′ /t
−0.4
−0.3
−0.1
0.0
0.3

(33)

SC
where εAFk and εk are independently optimized. In Sect. 6.1,
we study the stability against phase separation (PS) and discuss whether charge fluctuation thereby correlates with theenhancement of d-SC. In Sect. 6.2, we consider the mechanism for the coexistence or mutual exclusivity of AF andd-SC orders. In Sect. 6.3, the notion treated in Sect. 6.2 is applied to the relationship between the staggered flux and d-SC
states. In Sect. 6.4, we discuss the relationship between the
pocket Fermi surfaces in the type-(ii) AF state and the Fermi
arcs observed in the pseudogap phase of cuprates.

c2 (L = 10)
2.42(8)1.85(9)0.509(6)−0.551(5)0.830(4)c2 (L = 12)—
1.95(8)0.323(8)−0.553(7)—are shown as a function of the doping rate. First, we discuss
the range in which the finite AF order occurs. As comparedwith the pure AF state ΨAF ,29) the value of δ at which m vanishes (δAF ) is hardly changes for t′ /t = 0: δAF ∼ 0.16, whereasδAF somewhat increases for a large |t′ /t|. This small change inδAF stems from the small energy diﬀerence between Ψmix andΨAF (or Ψd ), as shown in Table III.
We next consider the stability against PS. This property
is often judged by the sign of the charge compressibility κ[= (1 − δ)2 χc ] or charge susceptibility χc [= (∂2 E/∂δ2 )−1 ].
For χc > 0 (χc < 0), the state is stable against (unstable toward) PS. Thus, we need to consider the δ dependence of E/t(Fig. 22). Similarly to for ΨAF ,29) we find for Ψmix that E(δ)/tis fitted well by the parabolic form

6.1 Stability against phase separation
Before discussing Ψmix , we refer to known aspects of the
intrinsic stability of ΨN , Ψd , and ΨAF against PS. Except for
the limit of δ → 0, at which the anomaly of the Mott transition appears, the normal state ΨN is stable against PS.29) Asfor Ψd , E/t is a linear function of δ (χc → ∞) for a small δ, as
we will discuss later, indicating that the stability against PS is
marginal. However, this is distinct from the apparent instability of ΨAF toward PS. In the second row of Table IV, we summarize the conclusions of related VMC studies on the stability
against PS of the AF and mixed states. The pure (not mixed)
AF state is known to be unstable toward PS for t′ /t ∼ 04, 29)
but stable for t′ /t ∼ ±0.3.29) A mixed state in which BRE areAF
20)
introduced into εSCk but the AF part is fixed as εk = γk′exhibits instability toward PS for both t /t = 0 and −0.3. To
summarize, states with AF orders exhibit a tendency toward
PS depending on the value of t′ /t.
We study this property for Ψmix [Eq. (33)]. In Fig. 22, thetotal energy and sublattice magnetization [Eq. (29)] in Ψmix

E(δ)/t ≃ c0 + c1 δ + c2 δ2(34)

in the whole AF range (δ < δAF ); we have a unique value χc =c−1
2 in the AF phase. The values of c2 thus estimated for somevalues of t′ /t and L are summarized in Table V. It reveals thatc2 (namely χc ) becomes negative only for a narrow range neart′ /t = 0, minutely tL′ < t′ ≲ 0.2t (see Fig. 27 later). This aspect
is basically the same as that of the pure AF state.29) Thus, theinstability toward charge inhomogeneity originates in the AForder and is not directly connected with SC, as we will discussshortly.
13

J. Phys. Soc. Jpn.

FULL PAPERS

0.04

0.02

҄

SC
χ c=

Ed /t

-0.4

-0.02

E/t

Ed /t, Eh /t

0

-0.04

҄ L

E / t Pd

0.02

10
12
14
16

Pd

-0.3 10-0.3 12-0.1 10-0.1 120.3 100.3 12d-wave

0.01

҄

-0.6

δχ
-0.8

t′/t L-0.3 10-0.3 12-0.1 10-0.1 120.3 100.3 12

-0.08

-0.1

(a) t' / t = +0.30

Eh /t

-1

0

χc=
0

0.05

0.1

δ

0.15

0.2

҄

SC
E/t

-0.4

E/t

Fig. 23. (Color online) Two elements [Ed (blue) and Eh (red)] of diagonalhopping energy Et′ (= Ed + Eh ) plotted as functions of doping rate for threevalues of t′ /t. Data for L = 10 and 12 with U/t = 12 are shown.

Now, we identify the origin of the stability against PS forlarge values of |t′ /t|. First, we analyze E/t by dividing it into
its components EU /t, Et /t, and Et′ /t. Because Et /t (EU /t) isalmost linear (somewhat convex) as a function of δ for any
value of t′ /t and U > Uc (not shown), these components do
not contribute to phase stability. On the other hand, Et′ /t isconcave for any t′ /t, but, of course, the degree of concavity
diminishes as |t′ /t| decreases and vanishes at t′ /t = 0. We further analyze Et′ by dividing it into the two components Ed andEh (Et′ = Ed + Eh ), where Ed (Eh ) is the contribution of diagonal hopping that changes (does not change) the number ofdoublons.4) In other words, Ed is generated by the creation or
annihilation of D-H pairs, while Eh is generated by the hopping of doped (isolated) holes. In Fig. 23, we show the δ dependences of Ed and Eh for three values of t′ /t. We find thatboth Ed and Eh are concave but the curvature is much sharperfor Eh . To summarize, diagonal hopping (t′ term), especiallythat of doped holes, brings about intrinsic stability against PS.
A recent VMC study5) argued that the increase in χc has a
one-to-one correspondence with the enhancement of SC order in the wave function used. We check this point for thepresent Ψd and Ψmix . First, we discuss the pure SC state, Ψd ,whose δ dependence of E/t for t′ /t = 0.3 and 0 is shownin Fig. 24. Aside from a Mott anomaly for δ → 0, E/t becomes almost linear (χc tends to diverge) for δ ≲ δχ (spinodal point), while E/t becomes concave (χc remains moderate) for δ > δχ . Note that χc does not become negative unlike the case of ΨAF . Such behavior of E/t is preserved if
t′ /t is varied, but the range of χc → ∞ shrinks as t′ /t decreases; δχ ∼ 0.17, 0.15, and 0.12 for t′ /t = 0.3, 0, and −0.3,
respectively. On the other hand, the SC correlation function
exhibits the opposite behavior. As shown in Fig. 24, P∞
d exhibits a well-known dome shape and the SC order is perceptible for 0 < δ < δSC . Because the statistical fluctuation ofP∞
d becomes large for δ ∼ δSC , we estimate δSC very roughly
by the condition that the optimized gap parameter ∆d /t be-

δχ
-0.6

Pd҄ L
10
12
14
16

0.02

d-wave

0.01

҄

Pd

-0.06

-0.8

(b) t' / t = 0.0
0

0.1

0.2

0
0.3

δ
Fig. 24. (Color online) Total energy and d-wave SC correlation function(right axis) obtained for Ψd with U/t = 12 plotted for (a) t′ /t = 0.3 and(b) t′ /t = 0 as functions of doping rate. Data for four values of L are shown.
Above each panel, the ranges of appreciable SC and where χc = ∞ are shown.
The straight dash-dotted line for E/t is a guide for determining δχ .

comes 0.03 (∆d /t < 0.03 for δ > δSC ). We confirmed aknown tendency that δSC increases as t′ /t decreases [for instance, see Fig. 25(d) in Ref. 4]; δSC ∼ 0.20, 0.27, and 0.31
for t′ /t = 0.3, 0, and −0.3, respectively. Thus, the behaviorsof δχ and δSC as functions of t′ /t are opposite; the increase inχc rather has a negative correlation with the magnitude of SCin Ψd .
Next, we consider the case of Ψmix . As mentioned above,
the range of χc < 0 is included in the regime of type-(i) AF,−0.05 ≲ t′ /t ≲ 0.2. As an example, we show in Fig. 25 the δdependence of E/t for t′ /t = 0. We repeat that the area where
E/t is convex precisely coincides with that of finite m (δ <δAF ) indicated by a green arrow. For δ > δAF , where the stateis SC, E/t is concave. Furthermore, Pd is smooth at δ = δAF
and not particularly enhanced in the area of χc < 0. Anyway,as δ increases, after the AF order (or instability toward PS)
vanishes at δAF ∼ 0.16, SC survives up to δSC ∼ 0.27 fort′ /t = 0. In contrast, for t′ /t = 0.1 [as in Fig. 26(a)], the SC
first becomes weak at δ ∼ 0.12, but the area of χc < 0 (andAF order) continues up to δ ∼ 0.18. The relative sizes of δ
where χc < 0 and Pd > 0 are reversed as t′ /t varies.
Through the above analyses, we can conclude that the instability toward PS does not directly correlate with d-SC, al14

J. Phys. Soc. Jpn.

FULL PAPERS

0.025

(a)

Pd

0.02

m
0.02
L

10
12
10
10
12

1

0.8

10

Pd 12

0.015

14
16

0.01

0.6

10

m

E/t

0.01

E/t12

0.4

14
16

0.005

Pd

0.005
0

0.0
0.0
0.1
0.3
0.3

m

Pd

Pd

0.015

t′/t L

0.2

m
0

0.05

0.1

0.15

δ

0.2

0.25

0.3

0.35

0
0.025

Fig. 25. (Color online) Total energy (arbitrary unit) and d-wave SC correlation function obtained for Ψmix with U/t = 12 and t′ /t = 0 plotted as afunction of doping rate. Instead of E/t, we plot (E − c˜0 − c˜1 δ)/t with c˜0 andc˜1 appropriately adjusted to emphasize the curvature of E/t. The area of finite
staggered magnetization is shown by a green arrow.

0

0.05

0.1

0.15

0.2

(b)

m

Pd

0.02

0.005

0

0

t′/t L
-0.4
-0.3
-0.3
-0.2
-0.1
-0.1

0.05

m

Pd

0.01

though the ranges of t′ /t where SC and PS arise are similar as shown in Fig. 27. As discussed in Refs. 4 and 55, we
consider that the AF spin correlation and the suppression of
charge fluctuation owing to the Mott physics are responsible
for the behavior of the d-wave SC. We will return to this topic
in Sect. 6.2.
Finally, we emphasize the importance of BRE again. As
shown in Table IV, a mixed state in which BRE are considered only in εSCk exhibits instability toward PS even fort′ /t = −0.3.20) In this mixed state, εAF
k is fixed at γk [Eq. (13)],′which is similar to the optimized εAF
k for t /t = 0 (t1 , t2 ∼ 0,
see Fig. 13 for instance) in the PS area. This means that theSC
BRE on εAFk (independent of the BRE on εk ) are crucial forthis property.

0.8

0.6

Pd

0.015

0
0.25
1

m

10
10
12
10
10
12

0.4

0.2

0.1

δ

0.15

0.2

0
0.25

Fig. 26. (Color online) Doping-rate dependence of d-wave SC correlation
function Pd (left axis) and staggered magnetization m (right axis) shown for
U/t = 12. (a) Cases of t′ > tL′ [type-(i) regime] and (b) those of t′ < tL′ [type(ii) regime] for both L = 10 and 12. The Lifshitz transition point of the pure
AF state is tL′ /t ∼ −0.05.

is almost zero for δ < δAF and grows after the AF order vanishes (δ > δAF ). Thus, the two orders are mutually exclusive.
More accurately, in panel (b), a narrow range of coexistence
exists near the boundary δ = δAF for small |t′ /t|, typically fort′ /t = −0.1. The boundary between coexistence and mutual
exclusivity is situated at t′ = tL ∼ −0.05t, which is consistent
with the previous results.5, 7, 16, 19, 20) In the present results, it
seems that the AF state is always more robust than the d-SC
state and that the features of the underlying AF state control
whether d-SC appears or not. We will return to these points
shortly.
On the basis of the results for Ψmix above, we constructed
the phase diagram in the δ-t′ space shown in Fig. 27. In accordance with Fig. 17 for the pure states, the AF state occupies
a wide area. Except for the range of −0.1 ≲ t′ /t ≲ 0.2, SC
does not appear for low doping rates (δ ≲ 0.2). Furthermore,as mentioned, χc becomes negative for tL′ < t′ ≲ 0.2t. The
state phase separates into an AF state at half filling and a statein the overdoped regime (δ ≳ 0.15). Therefore, homogeneous
SC does not appear in the underdoped regime for any valueof t′ /t. This result greatly modifies the results of previous
VMC studies without BRE, in which d-SC widely prevails fort′ /t < 0, but is consistent with recent results of studies aapplying many-parameter VMC methods to Hubbard-type models5) and a d-p model56) and a study employing DMET.7) Such
predominance of the long-range AF phase is inconsistent with

6.2 Coexistence or mutual exclusivity of AF and d-SC orders
Previous studies using various mixed states with
BRE5, 16, 19, 20) and a recent study using density matrix
embedding theory (DMET)7) argued that the orders of AF
and d-SC are coexisting or mutually exclusive according towhether t′ /t ∼ 0 or t′ /t ≲ −0.1. Here, we systematically study
this point for Ψmix and deduce the origin of the coexistenceof the two orders, which is closely related to the mechanismof d-SC.
In Fig. 26, we show the δ dependence of the d-SC correlation function and staggered magnetization [Eq. (29)] measured in Ψmix . For Ψmix , we represent the d-SC correlation
function by Pd ≡ Pd (R) [Eq. (22)] with R being the vector
connecting the distant points in the system used [For instance,
R = (5, 5) for a system with L = 10], because we focus ona strongly correlated regime (See Appendix C in Ref. 4). We
show the results separately for the type-(i) AF and type-(ii)AF regimes in panels (a) and (b), respectively, because thefeatures are distinct in the two regimes. In the type-(i) regime[panel (a)], the SC order (Pd ) arises or vanishes regardlessof whether the AF order (m) is present or absent. For example, for t′ /t ∼ 0, AF and SC long-range orders coexist forδ < δAF and SC remains for δAF < δ < δSC as a pure SC order. On the other hand, in the type-(ii) regime [panel (b)], Pd
15

J. Phys. Soc. Jpn.

FULL PAPERS(ii) A d-SC gap ∆k with a similar form to Eq. (8) has a large
magnitude at (π, 0).
(iii) The scattering of q = (π, π), which is induced by the
AF exchange correlation between nearest-neighbor sites, ispossible by connecting two antinodal points with opposite
signs of ∆k .
In Fig. 29(b), we plot n(k) obtained in Ψmix for the same
parameter sets as in Fig. 29(a). Corresponding contour maps
are displayed in Fig. 30. For t′ /t = ±0.3, the results for Ψmix
are almost the same as those for ΨAF shown in Figs. 19(b)and 19(c) because the SC order does not appear. The results
are also similar in Figs. 20(c) and 30(a) for t′ /t = −0.3.
However, for t′ /t = 0, where SC appears, the pocket Fermi
surfaces at the antinodes in ΨAF in Fig. 19(b) are replacedwith gap behavior (green) similar to the decreasing slope in
Ψd in Fig. 29(a). It is clearer to compare Fig. 20(b) for ΨAF
with Fig. 30(b) for Ψmix . This reveals that for the d-SC order,
Fermi surfaces in the nodal directions are not necessary but
gap formation in the antinodes is vital. Incidentally, the resultant SC in the coexisting state, if any, does not have a featureof cuprate SCs, namely, nodal Fermi surfaces [Fig. 29(b)];they are smeared out by an AF gap. To provide an overview
of this topic, we summarize in Table VI the locations of thelocal Fermi surface centers of the three states for typical values of t′ /t. From this table with the above discussion, we may
derive two requisites for d-SC in the mixed state:(I) In the underlying pure AF (or normal) state, Fermi
surfaces exist in the antinodes [around (π, 0) and equivalentpoints].
(II) The hot spots determined by εSCk (see Sect. 3.1) are situated in the Fermi surface area mentioned in (I).
On the basis of these conditions, we can explain the evolution of the states realized in Ψmix mentioned in Table VI. We
show the main point schematically in Fig. 31. For t′ < tL′ , item(I) is not satisfied for a small δ, and d-SC does not emerge
as shown in Fig. 31(d). However, as δ approaches δAF , theedge of the Fermi surface centered at (π/2, π/2) extends tothe antinodes, as will be shown in Fig. 34(c). The scattering therein possibly yields a narrow window of coexistence,for example, for t′ /t = −0.1 [δ ∼ 0.12 and ∼ 0.139 for
L = 10 and 12, respectively] in Fig. 26(b). Regarding item
(II), the hot spots stay near the antinodes in this range of t′ /t[Fig. 31(a)]. On the other hand, for t′ > tL′ , item (I) is satisfied. For a small |t′ /t|, item (II) is also satisfied [Fig. 31(c)],so that a coexisting state appears as in Fig. 26(a). However,as t′ /t increases, the hot spots shift toward the nodal area [redin Fig. 31(a)] and deviate from the Fermi surface range in theantinodes [Fig. 31(b)], which is relatively narrow as shown
later in Fig. 35. Consequently, d-SC does not appear appreciably for t′ /t = 0.3, as shown in Fig. 26(a). This behavior
contrasts with that of the pure d-SC state (Fig. 24), in which
d-SC becomes weak more slowly because the hot spots are always situated at the Fermi surface of the underlying state ΨN ,and the scattering intensity becomes weak as the hot spots
move away from the antinodes.
To summarize, because the AF state underlies the d-SC order, substantial d-SC arises only when the scattering of Q in
the antinodes is compatible with the AF behavior. The requisites for this are given by (I) and (II) above.

SC

AF

AF

type (i)type (ii)PSt’/tFig. 27. (Color online) Phase diagram in δ-t′ space for U/t = 12 constructed for mixed state Ψmix . The purple shaded area indicates the regime
unstable toward phase separation, which is limited to within the type-(i) AF
phase. The bold red dotted line indicates the Lifshitz transition boundary tL′ /t.

Table VI. Locations of the centers of local Fermi surfaces in the state realized for a small δ (leftmost state for the mixed state) and U/t = 12 summarized for the AF, d-SC, and mixed states. For the mixed state, the evolutionof the realized states as δ increases is shown for δ ≲ 0.3. ‘Co’ (‘N’) indicatesa coexisting state with AF and d-SC orders (normal state).
t′ /t

AF

d-SC

−0.3
−0.1
0.0
0.1
0.3

(π/2, π/2)
(π/2, π/2)
(π, 0)(π, 0)(π, 0)

Always(π/2, π/2)

Mixed
Evolution of stateAF(ii) → SC
AF(ii) → (Co) → SC
Co → SC → N
Co → AF(i) → N
AF(i) → N

(π/2, π/2)
(π/2, π/2)NoNo
(π, 0)

the results of experiments on hole-doped cuprates as well as
recent advanced studies on electron-doped cuprates.57–59) We
will discuss this point in Sect. 7.
Next, we consider why a d-SC order can coexist with a
type-(i) AF order but is incompatible with a type-(ii) AF order. We can deduce the reason by considering the location of
the Fermi surface in the underlying pure AF state. First, we
review relevant properties of the d-SC state. In Figs. 28(b)and 28(c), we show contour maps of n(k) for Ψd with t′ /t = 0
and t′ /t = −0.3, respectively. The steep slope of n(k), indicative of a Fermi surface, exists only near (±π/2, ±π/2), and thegentle slopes around (±π, 0) and (0, ±π) indicate gaps, in contrast with the feature of the normal state shown in Figs. 28(a),which clearly exhibits a Fermi surface in all directions. In
Fig. 29(a), n(k) in Ψd [corresponding to Figs. 28(b) and 28(c)]is shown along the same path as in Fig. 19 for three values oft′ /t. As t′ /t varies, n(k) around (π, 0) greatly varies but thenodal Fermi surface near (π/2, π/2)21) hardly changes.4) This
indicates that the electronic states near (π, 0) are closely related to SC, because properties associated with SC such as Pd
greatly change with t′ /t. Actually, antinodal Fermi surfaces
have the following advantages for d-SC on the square lattice:(i) The density of states diverges at (π, 0) owing to a van
Hove singularity for |t′ /t| ≤ 0.5.
16

J. Phys. Soc. Jpn.

FULL PAPERSδ ≈ 0.08 (a) Normal(BR), t′/t=-0.3(b) d-wave, t′/t=0(b) d-wave, t′/t=-0.3

0.25

0.20
0.25

-3

-2

-1

0

0.70

0.70

1

2

0.75

0.75

3

-3

-2

-1

0

kx

1

kx

2

3

-3

-2

-1

0

ky

0.70

0.70

1

2

3

kx

Fig. 28. (Color online) Contour maps of n(k) at U/t = 12 and δ ∼ 0.08 shown for (a)optimized normal (paramagnetic) state ΨN with t′ /t = −0.3 (L = 10–18)
and for (b) and (c) optimized pure d-wave pairing state Ψd with t′ /t = 0 (b) and −0.3 (c) (L = 10–16).

1(a) d-wave
dispersion of a staggered flux (or d-density wave) state.3, 47)
Although this state has been extensively studied as a candidate
for the pseudogap state as well as the ground state of cuprates,
here we avoid referring to various interesting aspects of thisstate and focus on its bare dispersion:
√
2
2
εSF
k = −2t cos k x + 2 cos 2θ cos k x cos ky + cos ky , (35)

L t′/t0.8

12 -0.3
16 -0.3
12 0.016 0.012 0.3
16 0.3

n(k)

0.6

0.4

where θ corresponds to a quarter of the magnetic flux penetrating each plaquette of the square lattice and is treated asa variational parameter here. For θ = 0, εSFk is reduced toγk [Eq. (13)]; for θ = π/4 (π-flux state), εSFk at half fillingyields a Dirac cone with a linear dispersion having apices at(±π/2, ±π/2). In Fig. 32, we show the Fermi surfaces generated by εSF
k in the first quadrant of the Brillouin zone for twovalues of θ and some values of δ for each θ. At half filling,
the Fermi surface is the apex of an elongated Dirac cone at(π/2, π/2). For δ > 0, a Fermi surface appears as a slice ofan elongated Dirac cone around the nodal point (π/2, π/2).
Gaps open in the antinodes around (π, 0) and (0, π). The formof the pocket Fermi surfaces and the antinodal gaps resembles the features in the pseudogap phase of cuprates. Notethat the pocket Fermi surface becomes slender and its edge
approaches the antinodes as θ decreases and/or δ increases.
Here, we study how the energy in Ψd [Eq. (19)] varies whenSC
we use εSFk instead of γk as εk . If the coexistence of staggered
flux and d-wave SC orders is favored, the energy in Ψd may be
reduced at a finite value of θ. In Fig. 33, we show the increasein energy per site ∆E as compared with that in Ψd with γk asa function of θ for t′ /t = 0 and U/t = 12. For large valuesof θ (≳ 0.05π), the energy markedly increases regardless of δ.
On the other hand, for a small θ and large δ, ∆E is small orslightly negative, as indicated by the arrow, meaning that thetwo orders possibly coexist. In these cases, the Fermi surfaces
reach the antinodes. This is consistent with the notion that the
gap in the antinodes in εk for the underlying state is unfavorable to the d-SC order. To summarize, a robust staggered flux
order and a d-wave SC order are unlikely to coexist, although,to ensure this conclusion, we should investigate an appropriate mixed state of the two orders.

0.2

U/t=12 δ ≈ 0.08
01(b) Mixed
L t′/t

n(k)

0.8

0.6

12 -0.3
16 -0.3
12 0.016 0.012 0.3
16 0.3

AF
AF+SC
AF

0.4

0.2

0

(0,0)

(π,0)

k

(π,π)

(0,0)

Fig. 29. (Color online) Momentum distribution functions for (a) d-SC stateand (b) mixed state for U/t = 12 and δ ∼ 0.08 compared among three valuesof t′ /t. In (a), a nodal Fermi surface is indicated by an arrow. Data for L = 12
and 16 are plotted.

6.3 Coexistence of d-wave SC and staggered flux orders
To highlight the importance of Fermi surfaces in the antinodes for inducing a d-wave SC order, we consider the bare17

J. Phys. Soc. Jpn.

FULL PAPERS

Mixed, δ ≈ 0.08 (a) t′/t=-0.3(b) t′/t=0(c) t′/t=0.3

0.25

0.25
0.25

0.70
0.75

-3

-2

-1

0

1

0.75
0.70

2

3

-3

-2

-1

kx

0

1

0.75
0.70

2

3

-3

-2

-1

0

kx

1

ky

2

3

kx

Fig. 30. (Color online) Contour maps of n(k) for the optimized mixed state at U/t = 12 and δ ∼ 0.08 shown for (a) t′ /t = −0.3, (b) 0, and (c) 0.3. The datain these maps include the data used in Fig. 29(b). Systems with L = 10–16 are used.

(a) Quasi-Fermi surface(b) t' / t = 0.3, AF (i)1㸫㸩t' / t
0.3
0.0
0.3

0.8

ky /π

0.6

1

Q

0.4

Gap 0.4

d-waveU / t = 12
δ = 0.08

0

0.2

0.8

(c) t' / t = 0, AF (i)

1

(a)

0.3
0.2

0.8

0.4
0.6
kx /π

0.4

Impossible

ky / π0.2(d) t' / t = ࡳ 0.3, AF (ii)θ / π = 0.125
δ = 0.1
δ = 0 Magnetic BZ

0.6

boundaryδ =0.5
0.4

Fermi
surface

0.2

Q

Q

Gap

0
1

PossibleY

Impossible

0.2

0.8

Fig. 31. (Color online) (a) Quasi-Fermi surface in Ψd obtained using εSCkin the first quadrant of the Brillouin zone for U/t = 12 and δ = 0.08 for three′
typical values of t /t. The gray dash-dotted line indicates the AF Brillouinzone boundary. (b)-(d) Possibility of scattering of q = Q in Ψmix in the full
Brillouin zone. The hot spots in the SC part are indicated by circles. Data forΨAF with U/t = 12 and δ = 0.0816 (L = 14) are used to draw the contours.

M(b)δ = 0.1

θ / π = 0.0625

ky / πUnoccupied
0.6

δ=0

0.4

Occupied
0.2
Staggered flux stateΓ

6.4 Possible relation with pseudogap
One of the anomalous features arising in the pseudogap
phase (T c < T < T ∗ ) of underdoped cuprates is the Fermi
arcs60) observed in ARPES spectra, namely, unclosed Fermi
surfaces whose centers are situated in the nodal directions
near (±π/2, ±π/2), and similarly Fermi surface pockets.61, 62)
If T is fixed, a Fermi arc becomes longer as δ increases andbecomes connected to other arcs in adjacent quadrants of theBrillouin zone to form an ordinary closed Fermi surface at thephase boundary (T = T ∗ ). The origin of the pseudogap has
not yet been elucidated. First, as a possible candidate for the
Fermi arcs, we consider the pocket Fermi surface of a doped
ΨAF .

0

0.2

0.4
0.6
kx / πX
0.8

1Fig. 32. (Color online) Fermi surfaces of a staggered flux state for two values of θ and various doping rates drawn in the first quadrant of the Brillouinzone: Γ = (0, 0), X = (π, 0), M = (π, π), and Y = (0, π). The thickness of the
Fermi lines [e.g. for δ = 0.2 in (b)] indicates that εSF
k is relatively flat.

As shown in Fig. 20(c), a pocket Fermi surface of a type(ii) AF state is formed around (π/2, π/2) and is similar to the
Fermi arc observed by ARPES. ΨAF has energy gaps around18

J. Phys. Soc. Jpn.

∆E / t0.006

L=12 10

FULL PAPERS

0.0
0.04
~0.08
0.12
0.16

0.004

following subjects: (A) Ground-state phase diagrams in thespace of U/t, t′ /t, and δ. (B) In what regime and through whatmechanisms does the coexistence of AF and d-SC arise? (C)In what regime and from what cause does instability toward
inhomogeneous phases occur? First, we itemize the main results in this work:(1) In the d-SC state, the eﬀective band εSC
k is markedly
renormalized for the model parameters of U/t ≳ 6, a large|t′ /t|, and a small δ (≲ 0.1) (Figs. 4 and 5), as known previously. We found, however, owing to BRE, not only is the
improvement in energy much smaller than those in the normal and AF states, but also quantities related to SC [Pd , S (q),n(k)] are modified only very slightly (Figs. 6-8).
(2) In the normal state, BRE apply U ≳ Uc , δ ≲ 0.05, and|t′ /t| ≳ 0.1 with Uc /t being the Mott transition point (Fig. 10).
The improvement in energy is an order of magnitude larger
than that of the d-SC state but an order of magnitude smaller
than that of the AF state (Fig. 11).
(3) In all the states studied, band renormalization takes
place to reduce the kinetic energy (Et ) at the cost of the interaction energy (EU ), which corresponds to the tendency ofa strongly correlated state to undergo a phase transition to
reduce the kinetic energy.4, 35) In the resultant renormalized
band, the nesting condition tends to be restored (t1 /t → 0).
(4) For the AF state, BRE are useful for reducing the energy, especially for t′ /t < 0 (Fig. 14); the qualitative features
are almost independent of U/t for U > UAF . As a result, theAF state occupies a wide area (δ ≲ 0.2) in the phase diagrams (Figs. 17 and 27). The AF area is considerably wider
for t′ /t = −0.3 than for t′ /t = 0, in contrast to the results without BRE. In a doped metallic AF state, as t′ /t is varied, a typeof first-order Lifshitz transition takes place at t′ = tL ∼ −0.05t
regardless of the values of U/t and δ. In the type-(i) [(ii)] AF
regime (t′ > tL ) [(t′ < tL )], local pocket Fermi surfaces arisearound (π, 0) [(π/2, π/2)] and equivalent points (Figs. 19 and20). This diﬀerence plays a critical role in inducing the d-SC
mentioned in (6) before. The Fermi surface in the type-(ii) AF
is possibly related to the Fermi arcs found in cuprates.
(5) In the mixed state, the range of instability toward phaseseparation (PS) is found to be tL′ /t < t′ /t ≲ 0.2, similarly
to in the AF states.4, 29) The AF order is responsible for thisinstability, which does not directly correlate with d-SC. Elsewhere, the state is stable against PS. This stability is mainly
due to the diagonal hopping of doped carriers.
(6)The coexistence or mutual exclusivity of AF and d-SC
orders was studied in the mixed state (Fig. 26). The AF order has preferentially exhibits this property because the AFpart greatly reduces the energy compared with the SC part(Fig. 16). By examining various cases, we found two requisites for the d-SC order to arise (Fig. 31): (i) In the underlying
pure AF (or normal) state, Fermi surfaces exist in the antinodes [near (π, 0) and equivalent area]. (ii) The hot spots determined by εSC
k are situated in the Fermi surface area. These
requisites indicate that the scattering of q = (π, π) in the antinodes is vital for d-SC. Thus, the coexistence basically occurs
in the type-(i) AF regime. The range of t′ /t in which coexistence occurs is similar to that for the instability toward PS(Fig. 27), but this similarity is accidental. These requisites
seem to apply to the coexistence of d-SC and staggered flux
orders.
The present results are quantitatively consistent with recent

δd-waveSF(ε k )

U / t =12
t' / t = 00.002∆E = E(θ) - E(0)
0
0

0.05

θ/π

0.1

0.15

Fig. 33. (Color online) Diﬀerence in total energy upon introducing a stag′
gered flux 4θ in εSFk [Eq. (35)] for U/t = 12 and t /t = 0. The thick arrow
indicates the minimum (at θ/π ∼ 0.02) for δ = 0.16.

the antinodes in the sense that n(k) is smooth with a finite|∇n(k)|. In Fig. 34, we show contour maps of n(k) for different doping rates, where the other conditions are the same(t′ /t = −0.3, U/t = 12). This figure reveals how the pocket
Fermi surface evolves as δ increases; a small pocket Fermi
surface appears around (π/2, π/2) for very light doping, thearc length becomes long along the AF Brillouin zone boundary (π, 0)–(0, π), finally forming a connected Fermi surface
centered at Γ = (0, 0) for δ = 0.2245 (not shown), where the
AF order vanishes. This behavior is consistent with that of
the Fermi arc of cuprates. For the appearance of such behavior at a finite temperature, it is also important that the type-(ii)ΨAF has a very low energy and is stable against phase separation. Furthermore, the type-(ii) ΨAF does not coexist withd-SC except for δ ∼ δAF . Although this result cannot be directly applied to the pseudogap phase of cuprates because an
AF long-range order has not been observed, it is intriguingthat short-range AF orders of 20-30 lattice constants were observed up to high temperatures.8)
In Fig. 35, we show the evolution of contour maps of n(k)as δ increases in the type-(i) AF state (t′ /t = 0). In contrast to the type-(ii) AF state, a pocket Fermi surface grows
from the antinodes in the nodal directions and finally forms
a closed Fermi surface centered at Γ = (0, 0) for δ = 0.1633(not shown). Because energy gaps open for δ < δAF in the
nodal directions, the type-(i) AF state, which corresponds toelectron-doped cuprates, is not directly related to the Fermi
arc phenomena.
7. Summary and Discussion
In this paper, we studied band renormalization eﬀects(BRE) owing to electron correlation on a mixed state of d x2 −y2 wave pairing (d-SC) and antiferromagnetic (AF) orders, aswell as normal (paramagnetic), pure d-SC, and pure AF states,by applying a variational Monte Carlo (VMC) method to the
Hubbard (t-t′ -U) model. For the mixed state, BRE were introduced into the AF and d-SC parts independently; BRE on AF
orders, previously not investigated,28) markedly change theprevious knowledge of the Hubbard model. By searching ina wide model-parameter space with wave functions at various
levels, we obtained systematic insights, particularly into the19

J. Phys. Soc. Jpn.

FULL PAPERS

AF(BR), t′/t=-0.3 (a) δ = 0.04082(b) δ = 0.10204(c) δ = 0.14286

0.25
0.20
0.25

0.70

0.70

0.75

-3

-2

-1

0

0.70

0.75

1

2

3

-3

-2

-1

kx

0

1

ky

0.75

2

3

-3

-2

-1

kx

0

1

2

3

kx

Fig. 34. (Color online) Contour maps showing the evolution of n(k) as the doping rate increases in the type-(ii) AF state for t′ /t = −0.3 (U/t = 12). Three
typical doping rates are selected. Data for L = 14 are used.

AF(BR), t′/t=0.0 (a) δ = 0.04082(b) δ = 0.10204(c) δ = 0.14286

0.25
0.20
0.25

0.70

0.70

0.70

0.75

-3

-2

-1

0

kx

1

2

3

-3

-2

-1

0

1

ky

0.75

0.75

2

3

-3

kx

-2

-1

0

1

2

3

kx

Fig. 35. (Color online) Contour maps showing the evolution of n(k) as the doping rate increases in the type-(i) AF state for t′ /t = −0.3 (U/t = 12). Data forthe same doping rates as in Fig. 34 are displayed for comparison. Systems with L = 14 are used.

studies with advanced techniques,5–7) and make it possible toreasonably interpret individual features of previous studies.
Finally, we discuss the relationship with cuprates. The
present results that the AF order is predominant for a wide
range of model parameters (U/t ≳ 6, δ ≲ 0.2, most t′ /t)and that uniform d-SC disappears in the underdoped regime
are consistent with those of recent VMC,5) DMFT,6) andDMET7) studies based on the Hubbard model. Furthermore,
recent VMC studies on the t-J 63) and d-p56) models display
the same tendency. Nevertheless, these results are inconsistent with properties common to hole-doped cuprate SCs: theAF long-range order is broken by less than 5% doping withcarriers and high-T c d-SC appears in the underdoped regime.
In addition, it was recently shown that well-annealed electrondoped samples with small doping rates (5-10%) exhibit no AF
long-range orders but metallic or SC behavior57, 58) with entirely closed Fermi surfaces.59) Assuming that the AF orderis excluded for some reason, most properties of the remaining d-SC derived by theories so far are basically consistent
with those of cuprates. Thus, it is important to clarify why AF
long-range order is robust in the theory. It seems that the ap-

proximations applied are not responsible for the predominant
AF orders, but the models are lacking in certain important factors that destabilize AF orders. These factors may be disorder
or impurities inherent in cuprate SCs. It seems that theoreticalresearch on cuprate SCs may proceed to this direction.
After the submission of this paper, we noticed that BRE on
AF states were already considered in a VMC study of Watanabe, Shirakawa and Yunoki for three-band as well as singleband Hubbard models.64) They used the optimization method
mentioned as ‘an alternative approach’ in Appendix B. Theirresults are basically consistent with ours.
We thank Kenji Kobayashi, Masao Ogata, Shun Tamura,
Junya Otsuki, Yuta Toga, Hiroshi Watanabe, Kentaro Sato,and Masaki Fujita for useful discussions and information.
This work was supported in part by Grants-in-Aid from the
Ministry of Education, Culture, Sports, Science and Technology, Japan.
Appendix A:

Details of Optimization in Normal StateIn this Appendix, we explain how to actually deal with
the BR of the normal (paramagnetic) state ΨN [Eq. (26)] for20

J. Phys. Soc. Jpn.

-0.42
BR
(optim.)FULL PAPERS{k3}occ
Normal{k0}occ㸫t' / t 0.3δL δ ~ 0.08
0.08
10
0.0833 12
0.0781 16
0.0864 18

a3

a4

-0.35

a6
-0.44

Original-0.45

a5

Original
Extrapolation
BR

E/t

E/t

-0.43

a2
a1

δ = 0.08
U / t = 12
L = 10
-0.6

-0.4
A4
A3

a0

-0.4

{k3}occ -0.2

L=10

δ=0.08 t' / t
0.0
+0.3

-0.45 U / t = 12-0.5{k0}occ0{k2}occ 0

-0.5

0.5

t1 / tt' / tFig. A·2. (Color online) Energy expectation values of the normal state ΨNfor some model-parameter sets (L, δ, t′ /t, U/t = 12) plotted as functionsof the band parameter t1 /t. The data plotted by light-green diamonds anddenoted by Aℓ correspond to the set discussed in the text (10, 0.08, −0.3, 12).
Those plotted as dark circles are identical to those shown in Fig. 16.

Fig. A·1. (Color online) Illustration of how to obtain the bandrenormalized energy in ΨN within a single variational band parameter t1 /tfor a specific model parameter set (L = 10, δ = 0.08, U/t = 12). The totalenergy per site is plotted as a function of t′ /t (< 0). For details, see text.

finite-size systems. As mentioned in Sect. 2, ΨN depends only
on the choice of {k}occ (Fermi surface) and not explicitly on
εk . In the thermodynamic limit (L = ∞), where k is a continuous variable, {k}occ continuously changes as εk or the band
parameters therein (t1 /t, etc.) gradually vary. This means that{k}occ directly depends on the band parameters. On the other
hand, in the finite systems we treat here, for which the available k are discrete, {k}occ (namely ΨN ) is invariable in a certain range of band parameters (or εk ) and discontinuously
changes at the edges of the range. Generally, this range becomes wider for a smaller L. To begin with, we illustrate
this point assuming that the eﬀective band is given by εNk in
Eq. (5). Even for this simple form of εk , we believe that full
BRE are achieved in most cases.
To avoid confusion between t′ /t in H (model parameter)
and the variational band parameter t1 /t in ΨN , we start with
the noninteracting case (U = 0). In this case, the exact groundstate is given by Eq. (4), in which {k}occ is determined by thebare band dispersion ε̃k in Eq. (2), indicating that t1 = t′ andεk = ε̃k for U = 0 according to variation theory. If we decrease the sole band parameter t′ /t in H from zero, {k}occ is
switched from one configuration to another at certain discretevalues of t′ /t. In Fig. A·1, we show such evolution of {k}occ ,for L = 10 and δ = 0.08 as an example, with alternate red andblue arrows near the lower horizontal axis; {k}occ is switchedas{k0 }occ → {k1 }occ → {k2 }occ → {k3 }occ → · · ·

-0.45

A0

A2

-0.4

rectly on t1 /t, ΨN should exhibit completely flat energy as afunction of t1 /t in Aℓ and discontinuities at the edges of Aℓ .
In Fig. A·2, we show the t1 /t dependence of the total energy
for the above parameter set with light-green diamonds, alongwith the same quantity for other parameter sets. Because theeﬀective band dispersion εNk [Eq. (5)] in ΨN is assumed tohave the same form as the bare band dispersion ε̃k [Eq. (2)],65)
the division of the areas (Aℓ ) for t′ /t discussed above directly
corresponds to the division of t1 /t, as also marked by Aℓ in
Fig. A·2. The energy minimum for the above model parameterset (t′ /t = −0.3) is obtained not in A3 (including t1 /t = −0.3)
but in A2 , meaning that BRE manifest themselves.
Owing to this locally flat behavior of E/t, ordinary optimization techniques that require information about the gradient of E/t are inapplicable to ΨN . Here, we use another way of
optimization. Below, we describe its outline with an illustration in Fig. A·1 for a model-parameter set (L = 10, δ = 0.08,
U/t = 12) as an example. (i) Calculate the total energy E/tdensely as a function of t′ /t for a fixed set of the other parameters (L, δ, U/t) without introducing BRE, namely by putting
t1 = t′ . In Fig. A·1, the E/t thus obtained are plotted withsmall solid circles joined by a thick dashed line. We find that
E/t is described by a distinct nearly straight curve for eachAℓ . (ii) Each segmented curve (say in Aℓ ) can be well extrapolated using a first- or second-order least-squares method:
(ℓ) ′
(ℓ) ′
2
Eℓ (t′ /t) = c(ℓ)
0 + c1 (t /t) + c2 (t /t) .

(A·1)

at t′ /t ∼ −0.107, −0.137, −0.213, −0.357, · · · . Let Aℓ (ℓ: integer) denote the area of t′ /t where {k}occ = {kℓ }occ as shownin Fig. A·1, for example, A2 = [−0.213, −0.137]. Note thatwithin each Aℓ , the ground-state wave function ΨN (=ΦN ) is
unchanging but the energy changes with t′ /t owing to the diagonal hopping term.
Next, we consider interacting cases (U > 0). Let the model
parameters be fixed, for example, at L = 10, δ = 0.08,
t′ /t = −0.3, and U/t = 12. For such a parameter set, weneed to optimize ΨN by adjusting the band parameter t1 /t independently of t′ /t together with the correlation parameters.
Because ΦN in Eq. (4) depends only on {kℓ }occ and not di-(A·2)

The extrapolated curves are shown with thin dashed lines inFig. A·1 and practically coincide with the values of E/t actually calculated with {kℓ }occ (BRE) outside Aℓ , whose values
are shown with open circles joined by thin dull-green curves.
Therefore, we may substitute such extrapolated values for the
results of actual BRE calculations to save labor. (iii) The optimized energy allowing for BRE for a fixed value of t′ /t isgiven by the lowest extrapolated value among the all the Aℓ .
For t′ /t = −0.3, for example, the lowest energy is given by{k2 }occ , and the improvement in energy owing to BRE (∆E/t)
is indicated by a brown arrow. We actually estimated the optimized BRE energies of ΨN for most model parameter sets
21

J. Phys. Soc. Jpn.

0.06

FULL PAPERS

vant spikes. A given configuration {k}AFocc yields a plateau orplateaus with the same E. We determined Eopt by averaging E
in the lowest plateau and verifying that the estimated value issmoothly connected to those obtained using other model parameter sets. For t′ ∼ tL′ , the statistical fluctuations become
very large because multiple {k}AFocc have a value of E comparable to Eopt . Therefore, in this regime, we carried out up to
fifty calculations for a single model-parameter set, especiallyfor Ψmix .
As an alternative approach, we may optimize ΨAF andΨmix with a fixed {k}AFocc using the stochastic reconfiguration
method. By carrying out such operations for various values of{k}AFocc , we can identify the Ψ with the lowest E/t. Because thenumber of choices of {k}AF
occ rapidly increases as L increases,
we may adopt the method of choosing {k}occ used for ΨN in
Appendix A. Anyway, the task of optimizing ΨAF and Ψmix ismuch more burdensome than that for Ψd .

㸫U / t = 12
t' / t = 0.3
δ = 0.0833
L = 12

AF (BR4)(E Eopt) / t

0.04

㸫0.02
0

-0.02
0

1000

2000

iteration
Fig. B·1. (Color online) Evolution of the energy expectation value in ΨAF
obtained by VMC process using a simple linear optimization method. The
results of eight calculations successively performed are plotted in sequence,in each of which 320 linear optimizations were carried out. The initial parameter values in each calculation were set to those that yielded the lowest
plateau energy in the previous calculations. We estimated the optimized energy, in this case, by averaging the final results indicated by the arrow. In
averaging, we exclude scattered data that are more than twice the standard
deviation from the mean.

1) J. G. Bednorz and K. A. Müller, Z. Phys. B 64, 189 (1986).
2) For the recent status of experimental research, see, for instance, articles
in the special issue of J. Phys. Soc. Jpn. 81 vol. 1 (2012) on “Recent
developments in superconductivit”.
3) For the t-J model, see P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod.
Phys. 78, 17 (2006); M. Ogata and H. Fukuyama, Rep. Prog. Phys. 71,
036501 (2008). For the Hubbard model, see, for instance, Table I in
Ref. 5.
4) H. Yokoyama, M. Ogata, Y. Tanaka, K. Kobayashi, and H. Tsuchiura,
J. Phys. Soc. Jpn. 82, 014707 (2013).
5) T. Misawa and M. Imada, Phys. Rev. B 90, 115137 (2014).
6) J. Otsuki, H. Hafermann, and A. I. Lichtenstein, Phys. Rev. B 90,
235132 (2014).
7) B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93, 035126 (2016).
8) Y. Sidis, C. Ulrich, P. Bourges, C. Bernhard, C. Niedermayer, L. P. Regnault, N. H. Andersen, and B. Keimer, Phys. Rev. Lett. 86, 4100 (2001);
H. A. Mook, Pengcheng Dai, S. M. Hayden, A. Hiess, J. W. Lynn, S.-H.
Lee, and F. Dogan, Phys. Rev. B 66, 144513 (2002); J. A. Hodges, Y.
Sidis, P. Bourges, I. Mirebeau, M. Hennion, and X. Chaud, Phys. Rev.
B 66, 020501(R) (2002).
9) M. Capone and G. Kotliar, Phys. Rev. B 74, 054513 (2006).
10) D. Sénéchal, P.-L. Lavertu, M.-A. Marois, and A.-M. S. Tremblay,
Phys. Rev. Lett. 94, 156404 (2005).
11) M. Aichhorn, E. Arrigoni, M. Potthoﬀ, and W. Hanke, Phys. Rev. B 76,
224509 (2007).
12) S. S. Kancharla, B. Kyung, D. S 辿 n 辿 chal, M. Civelli, M. Capone,
G. Kotliar, and A.-M. S. Tremblay, Phys. Rev. B 77, 184516 (2008).
13) T. K. Lee and S. Feng, Phys. Rev. B 38, 11809 (1988).
14) A. Himeda and M. Ogata, Phys. Rev. B 60, R9935 (1999).
15) D. A. Ivanov, Phys. Rev. B 70, 104503 (2004).
16) C. T. Shih, Y. C. Chen, C. P. Chou, and T. K. Lee, Phys. Rev. B 70,
22502(R) (2004); C. T. Shih, J. J. Wu, Y. C. Chen, C. Y. Mou, C. P.
Chou, R. Eder, and T. K. Lee, Low Temp. Phys. 31, 757 (2005).
17) S. Pathak, V. B. Shenoy, M. Randeria, and N. Trivedi, Phys. Rev. Lett.
102, 027002 (2009).
18) T. Giamarchi and C. Lhuillier, Phys. Rev. B 43, 12943 (1991).
19) K. Kobayashi and H. Yokoyama, J. Phys. Chem. Solids 69, 3274
(2008); Physica C 469, 974 (2009).
20) K. Kobayashi and H. Yokoyama, JPS Conf. Proc. 1, 012120 (2014); JPS
Conf. Proc. 3, 015012 (2014); Phys. Proc. 58, 22 (2014); Phys. Proc. 65,
9 (2015).
21) A. Himeda and M. Ogata, Phys. Rev. Lett. 85, 4345 (2000).
22) C. T. Shih, T. K. Lee, R. Eder, C.-Y. Mou, and Y. C. Chen, Phys. Rev.
Lett. 92, 227002 (2004).
23) K. Kobayashi and H. Yokoyama, Physica C 463-465, 141 (2007).
24) T. Watanabe, H. Yokoyama, K. Shigeta, and M. Ogata, New J. Phys. 11,075011 (2009).
25) L. F. Tocchio, F. Becca, and C. Gros, Phys. Rev. B 86, 035102 (2012).
26) J. Liu, J. Schmalian, and N. Trivedi, Phys. Rev. Lett. 94, 127003 (2005).

through this procedure. To obtain other quantities, however,calculations using the optimized parameters are necessary.
Under the upper horizontal axis in Fig. A·1, we show theareas of {kℓ }occ that yield the optimized E/t with red and blue
arrows. It reveals that these areas of {kℓ }occ with BRE often
deviate from the areas of {kℓ }occ for bare cases shown near the
lower horizontal axis. Thus, in this model parameter set, theenergy reduction owing to BRE is brought about discontinuously as a function of t′ /t [see Fig. 10(b)]. In Fig. 16, we actually illustrate the above process of optimization associatedwith BRE for ΨN with L = 12, δ = 0.0833, and U/t = 12. The
red line indicates the optimized line for ΨN . In this parameter,
BRE are ineﬀective for or small −0.573 ≲ t′ /t ≲ 0.343.
Appendix B:

Details of Optimizing AF and Mixed States

In optimizing ΨAF and Ψmix , a similar diﬃculty exists inthe case of ΨN . Namely, if we determine {k}AF
occ according to
εAF
,
as
t
is
graduallyvaried,totalenergyE/tdiscontinuouslyηkchanges at a value where {k}AFisswitched
toanother configocc
uration. In contrast to ΨN , we have to optimize εAFk in addition to {k}AF
for
Ψ
and
Ψ
,
as
shown
in
Table
I. An addiAF
mix
occ
tional problem is that E/t depends on tη only very weakly. For
this reason, the stochastic reconfiguration method and quasiNewton methods did not work eﬀectively, and we returned toa primitive linear optimization method in most cases.
We show an example of optimizing ΨAF in Fig. B·1, where
the expectation value of E/t obtained in each linear optimization of the parameters is plotted for the specified model parameter set. Typically, 2.5 × 105 samples are used for the linear optimization. The expectation value of E does not monotonically decrease to the optimized value Eopt but irregularly
fluctuates, exhibiting wide and narrow plateaus and irrele22

J. Phys. Soc. Jpn.

FULL PAPERS

49) For instance, T. K. Lee and C. T. Shih, Phys. Rev. B 55, 5983 (1997),
and references therein.
50) For instance, see Ref. 3 and T. Tohyama and S. Maekawa, Supercond. Sci. Technol. 13, R17 (2000).
51) T. Xiang and J. M. Wheatley, Phys. Rev. B 54, R12653 (1996).
52) T. K. Lee, C.-M. Ho, and N. Nagaosa, Phys. Rev. Lett. 90, 067001(2003).
53) For instance, F. Ronning, C. Kim, D. L. Feng, D. S. Marshall, A. G.
Loeser, L. L. Miller, J. N. Eckstein, I. Bozovic, and Z.-X. Shen, Science
282, 2067 (1998).
54) N. P. Armitage, F. Ronning, D. H. Lu, C. Kim, A. Damascelli, K. M.
Shen, D. L. Feng, H. Eisaki, Z.-X. Shen, P. K. Mang, N. Kaneko, M.
Greven, Y. Onose, Y. Taguchi, and Y. Tokura, Phys. Rev. Lett. 88,
257001 (2002); N. P. Armitage, P. Fournier, and R. L. Greene, Rev.
Mod. Phys. 82, 2421 (2010).
55) F. C. Zhang, C. Gros, T. M. Rice, and H. Shiba, Supercond. Sci. Technol. 1, 36 (1988).
56) S. Tamura, Dr. Thesis, Faculty of Science, Tohoku University, Sendai(2016) [in Japanese].
57) A. Tsukada, Y. Krockenberger, M. Noda, H. Yamamoto, D. Manske,
L. Alﬀ, and M. Naito, Solid State Commun. 133, 427 (2005); O. Matsumoto, A. Utsuki, A. Tsukada, H. Yamamoto, T. Manabe, and M.
Naito, Physica C 469, 924 (2009).
58) T. Adachi, Y. Mori, A. Takahashi, M. Kato, T. Nishizaki, T. Sasaki, N.
Kobayashi, and Y. Koike, J. Phys. Soc. Jpn. 82, 063713 (2013).
59) M. Horio, T. Adachi, Y. Mori, A. Takahashi, T. Yoshida, H. Suzuki, L.
C. C. Ambolode II, K. Okazaki, K. Ono, H. Kumigashira, H. Anzai, M.

27) T. Watanabe, H. Yokoyama, Y. Tanaka, and J. Inoue, J. Phys. Soc. Jpn.
75, 074707 (2006); Phys. Rev. B 77, 214505 (2008).
28) In a recent VMC study with a diﬀerent formulation (Ref. 5), BRE on
AF parts seem to be implicitly introduced.
29) H. Yokoyama, R. Sato, S. Tamura, and M. Ogata, Phys. Proc. 65, 17(2015).
30) H. Yokoyama, R. Sato, and K. Kobayashi, to be published in Phys. Proc.
31) R. Sato and H. Yokoyama, to be published in Physica C.
32) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
33) T. A. Kaplan, P. Horsch, and P. Fulde, Phys. Rev. Lett. 49, 889 (1982).
34) H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 59, 3669 (1990).
35) H. Yokoyama, Y. Tanaka, M. Ogata, and H. Tsuchiura, J. Phys. Soc.
Jpn. 73, 1119 (2004).
36) H. Yokoyama, M. Ogata, and Y. Tanaka, J. Phys. Soc. Jpn. 75, 114706(2006).
37) The usual Hartree-Fock approximation for the completely nested case is
not suited to the case of t′ /t , 0, and the quasiparticle band dispersion√
U± ε2k + ∆2AF
Ek =2
is not appropriate for determining kF . However, the form of Eq. (18) isuseful for the AF part of the variational wave function even for t′ /t , 0.
38) J. P. Bouchaud, A. Georges, and C. Lhuillier, J. Phys. (Paris) 49, 553(1988).
39) W. L. McMillan, Phys. Rev. 138, A442 (1965).
40) D. Ceperley, G. V. Chester, and M. H. Kalos, Phys. Rev. B 16, 3081
(1977).
41) H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 56, 1490 (1987).
42) C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. Rev. Lett. 60,
1719 (1988).
43) S. Sorella, Phys. Rev. B 64, 024512 (2001); E. Neuscamman, C. J. Umrigar, and G. K.-L. Chan, Phys. Rev. B 85, 045103 (2012).
44) To consider hot spots more seriously, it is necessary to optimize the dwave gap form ∆d appropriately, as carried out in Refs. 24 and 45. In
this sense, the discussion here is rather broad.
45) H. Yoshimura and D. S. Hirashima, J. Phys. Soc. Jpn. 74, 712 (2005);
T. Watanabe, T. Miyata, H. Yokoyama, Y. Tanaka, and J. Inoue, J. Phys.
Soc. Jpn. 74, 1942 (2005).
46) G. Blumberg, A. Koitzsch, A. Gozar, B. S. Dennis, C. A. Kendziora, P.
Fournier, and R. L. Greene, Phys. Rev. Lett. 88, 107002 (2002); H. Matsui, K. Terashima, T. Sato, T. Takahashi, M. Fujita, and K. Yamada,
Phys. Rev. Lett. 95, 017003 (2005).
47) For instance, H. Yokoyama, S. Tamura, T. Watanabe, K. Kobayashi, andM. Ogata, Phys. Proc. 58, 14 (2014), and references therein.
48) This change from an insulator to a metal is also confirmed by the small|q| behavior of the charge density structure factor N(q).

Arita, H. Namatame, M. Taniguchi, D. Ootsuki, K. Sawada, M. Takahashi, T. Mizokawa, Y. Koike, and A. Fujimori, Nat. Commun. 7, 10567(2016).
60) T. Yoshida, M. Hashimoto, I. M. Vishik, Z.-X. Shen, and A. Fujimori,J. Phys. Soc. Jpn. 81, 011006 (2012).
61) N. Doiron-Leyraud, C. Proust, D. LeBoeuf, J. Levallois, J.-B. Bonnemaison, R. Liang, D. A. Bonn, W. N. Hardy, and L. Taillefer, Nature
447, 565 (2007).
62) N. Barišić, S. Badoux, M. K. Chan, C. Dorow, W, Tabis, B. Vignolle, G.
Yu, J. B 辿 ard, X. Zhao, C. Proust, and M. Greven, Nat. Phys. 9, 761(2013).
63) R. Sato, unpublished.
64) H. Watanabe, T. Shirakawa, and S. Yunoki, Phys. Rev. Lett. 110,027002 (2013).
65) It is possible to extend εk to more refined forms such as εΛk in Eq. (12)in a similar manner. We leave this for a future study.

23