Journal of the Physical Society of Japan

DRAFT

Staggered Flux State in Two-Dimensional Hubbard Models
Hisatoshi Yokoyama ∗ ,1 Shun Tamura2 , and Masao Ogata3
1

Department of Physics, Tohoku University, Sendai 980-8578, Japan
Department of Applied Physics, Nagoya University, Nagoya 464-8603, Japan3
Department of Physics, University of Tokyo, Bunkyo-ku, Tokyo 113-0033. Japan2

The stability and other properties of a staggered flux (SF) state or a correlated d-density wave state are studied
for the Hubbard (t-t′ -U) model on extended square lattices, as a low-lying state that competes with the d x2 −y2 -wave
superconductivity (d-SC) and possibly causes the pseudogap phenomena in underdoped high-T c cuprates and organic
κ-BEDT-TTF salts. In calculations, a variational Monte Carlo method is used. In the trial wave function, a configurationdependent phase factor, which is vital to treat a current-carrying state for a large U/t, is introduced in addition to ordinary
correlation factors. Varying U/t, t′ /t, and the doping rate (δ) systematically, we show that the SF state becomes morestable than the normal state (projected Fermi sea) for a strongly correlated (U/t & 5) and underdoped (δ . 0.16)
area. The decrease in energy is sizable, particularly in the area where Mott physics prevails and the circular current(order parameter) is strongly suppressed. These features are consistent with those for the t-J model. The eﬀect of thefrustration t′ /t plays a crucial role in preserving charge homogeneity and appropriately describing the behavior of holeand electron-doped cuprates and κ-BEDT-TTF salts. We argue that the SF state does not coexist with d-SC and is not a‘normal state’ from which d-SC arises. We also show that a spin current (flux or nematic) state is never stabilized in thesame regime.

staggered flux (SF) state (sometimes called a d-density wave
state)—as a possible pseudogap state for the Hubbard model.
We should understand such a state in the context of a doped
Mott insulator.3, 4) To respect the strong correlation, we use avariational Monte Carlo (VMC) method,31) which deals withthe local correlation factors exactly and has yielded consistent results for many aspects of cuprates.32–37) If the pseudogap phenomena are generated by a symmetry-breaking state,
it should be more stable than the (symmetry-preserved) ordinary normal state. Also, when a predominant antiferromagnetic (AF) or d-SC state is suppressed for some reason,
features of the symmetry-breaking state will manifest themselves. Note that a recent VMC calculation with a bandrenormalization eﬀect showed that an AF state is considerably stabilized compared with the d-SC state in a wide regionof the Hubbard model.38)
Since the early years of research on cuprate SCs, the SF
state has been studied by many groups from both weak- andstrong-correlation sides. In the early studies,39–45) the main
aim was to check whether the SF state becomes the groundstate, but it was shown mainly using the t-J model that the
SF state yields to other ordered states (AF and d-SC) for any
relevant parameters. Later, the SF state was mainly studied
as a candidate for a normal state that causes the pseudogap
phenomena and underlies d-SC in underdoped cuprates.46–50)
At half filling, for the Heisenberg model, owing to theSU(2) symmetry, the SF state is equivalent to the d-wave BCS
state,51, 52) which has a very low energy31, 32) comparable tothat of the AF ground state.53–55) The t-J model with finite
doping was studied using U(1) and SU(2) slave-boson meanfield theories56–58) and a perturbation theory of Hubbard X
operators,59) which revealed that the SF phase exists in phase
diagrams but is restricted to very small doping regions.58) As amore reliable treatment, VMC calculations32, 43, 44, 60) showedthat the SF state has lower energy than the projected Fermisea, although the d-SC state has even lower energy. It was

1. Introduction
Superconductivity (SC) in underdoped high-T c cuprates
should be understood through the relationship to the pseudogap phase observed for T c < T < T ∗ , where T c [T ∗ ] is thesuperconducting (SC) transition [pseudogap] temperature.1, 2)
Because the pseudogap phase appears in the proximity ofhalf filling, it is probably related to Mott insulators3, 4) (precisely, charge-transfer insulators5) ). Experimentally, the pseudogap phase presents various features distinct from an ordinary Fermi liquid.3, 4) (1) A large gap diﬀerent from the d x2 −y2 wave SC (d-SC) gap opens in the spin degree of freedom near
the momenta of (π, 0) and (0, π). (2) However, the material isconductive and does not have a charge gap. (3) Fragmentary
Fermi surfaces, i.e., Fermi arcs6, 7) or hole pockets,8, 9) appear
in the zone-diagonal direction near (π/2, π/2).
The origin of the pseudogap has often been studied as a
linkage to d-SC, although it will not be related to SC fluctuation.10–15) On the other hand, recent experimental studies argued that the pseudogap phase is accompanied bysome symmetry-breaking phase transitions at T ∗ .2) (1) Timereversal symmetry breaking16–19) is claimed from polarized
neutron scattering signals at the momentum (0, 0) as well asfrom the appearance of the Kerr eﬀect. (2) Rotational symmetry breaking (or nematic order) similar to the stripe phase
is observed, and the oxygen sites between copper atoms are
involved.2, 20) (3) Charge orders or charge density waves areobserved in resonant X-ray scattering experiments.21–23) (4)(π, π)-folded (shadow) bands appear in ARPES spectra, andso forth.24–27) Note, however, that the thermodynamic properties such as specific heat and spin susceptibility have not provided any evidence of the phase transition. It is also important
to study whether a pseudogap and other orders coexist or are
mutually exclusive.15, 28–30)
In this context, we study a symmetry-breaking state—a∗ E-mail, yoko@cmpt.phys.tohoku.ac.jp

1

J. Phys. Soc. Jpn.

DRAFT

pointed out that the δ dependence of the SC condensation energy using the SF state as a normal state becomes domelike60)
but that the SF state tends to be unstable toward phase separation.61) These VMC results claim that the strongly correlated
Hubbard model should have the same features.
For the Hubbard model, SF states have been studied using
a phenomenological theory,62) mean-field theories,63–65) andmore refined renormalization group methods66, 67) from theweak-correlation side. These studies obtained various knowledge of the SF state, but it is still unclear whether or not the
SF state is stabilized in the weakly as well as strongly correlated regions. On the other hand, a Gutzwiller approximation
study65) claimed that the SF state is not realized in the Hubbard model. A study using a Hubbard operator approach68)
showed the absence of SF order for a large U/t (= 8) unless an attractive intersite interaction is introduced. A study
using a dynamical cluster approximation for a 2 × 2 cluster69) argued that the circular-current susceptibility increases
in the pseudogap-temperature regime but does not diverge,
and there is no qualitative change as U/t and t′ /t are varied.
A study using a variational cluster approach70) concluded that
the SF phase is not stabilized with respect to the ordinary normal state for a strongly correlated region (U/t & 4). An extended dynamical-mean-field approximation showed that although the SF susceptibility is enhanced, it is dominated byd-SC and an inhomogeneous phase for t′ /t = 0.71) Thus, it
is still unclear whether the results in the Hubbard model areconsistent with those in the t-J model.
The purpose of this paper is to show that the SF state becomes considerably stable with respect to the projected Fermi
sea (an ordinary normal state) in the underdoped regime forlarge values of U/t and t′ /t ∼ −0.3 in the Hubbard (t-t′ -U)model, and to clarify various properties of this state on thebasis of systematic VMC calculations. It is essential to introduce a configuration-dependent phase factor to treat a currentcarrying state such as the SF state in the regime of Mott physics.72) Without it, the SF state is never stabilized in models
permitting double occupation such as the Hubbard model. Wechange the model parameters U/t, t′ /t, and the doping rateδ (= 1 − N/Ns ) in a wide range, with N and Ns being thenumbers of electrons and sites, respectively. Additionally, we
study the spin-current flux phase (sometimes called the spinnematic phase) using the same method.
Besides cuprates, we consider a model for layered organic
conductors, κ-(BEDT-TTF)2 X, [henceforth, abbreviated as κ(ET)2 X] with X being a univalent anion.73–75) In these compounds, SC arises for T c . 12 K, and a pseudogap behavior
similar to that of cuprates has been observed. Therefore, weneed to check whether its origin is identical to that of cuprates.
Various low-energy properties of κ-(ET)2 X are considered to
be described by the Hubbard model76) on an anisotropic twodimensional triangular lattice. The value of U/t can be controlled by applying pressure. U is estimated as U ∼ W–2W
with W being the band width.73) The degree of frustrationt′ /t can be varied by substituting X or applying uniaxial pressure. t′ /t is estimated by ab initio calculations as 0.4–0.7 forweakly frustrated compounds and ∼ 0.8 for the highly frustrated compound κ-(ET)2 Cu2 (CN)3 .73, 77) Among the former
compounds, deuterated κ-(ET)2 Cu[N(CN)2 ]Br (t′ /t ∼ 0.4)
under applied pressure has been shown to exhibit pseudogap
behavior such as a steep decrease in the NMR spin-lattice re-laxation time (1/T 1 T ) in the metallic phase (T > T c ). On theother hand, κ-(ET)2 Cu2 (CN)3 , which has a spin liquid state in
the insulating phase under ambient pressure, exhibits the Korringa relation (1/T 1 T = const.) in the metallic phase underpressure, namely, pseudogap behavior is absent.78) Furthermore, similar pseudogap behavior was observed in a holedoped κ-ET salt [κ-(ET)4 Hg2.89 Br8 ],79) in which the doping
rate is 0.11 and t′ /t ∼ 0.8. With these experimental results in
mind, we study the SF state on an anisotropic triangular lattice in the framework applied to the frustrated square latticefor cuprates.
This paper is organized as follows. In Sect. 2, we introduce
the model and method used in this paper. In Sects. 3 and 4, we
discuss the results mainly for the simple square lattice (t′ = 0)
at half filling and in doped cases, respectively, to grasp thecommon properties of the SF state. Section 5 is assigned to theeﬀect of the diagonal hopping term t′ for the frustrated squarelattice and anisotropic triangular lattice. In Sect. 6, we discuss
the results. In Sect. 7, we recapitulate this work. In Appendix
A, we summarize the fundamental features of the noninteracting SF state. In Appendix B, we briefly review the stability of
the SF phase for t-J-type models with new accurate data. In
Appendix C, we show that the spin current (flux) state is unstable toward the projected Fermi sea in t-J-type models for
any J (> 0) and δ. Preliminary results on the eﬀect of t′ terms
have been reported in two preceding publications.80, 81)
2.

Model and Wave Functions

In Sects. 2.1 and 2.2, we explain the model and variational
wave functions used in this paper, respectively. In Sect. 2.3,
we introduce a phase factor essential for treating a currentcarrying state in a strongly correlated regime. In Sect. 2.4, we
describe the numerical settings of our VMC calculations.
2.1 Hubbard model
As models of cuprates and κ-ET organic conductors, weconsider the following Hubbard model (U ≥ 0) on extended
square lattices (Fig. 1):
H==

Hkin + HU∑ (∑)−
ti j c†iσ c jσ + H.c. + U
n j↑ n j↓ ,(i, j),σ(1)jwhere n jσ = c†jσ c jσ and (i, j) indicates the sum of pairs onsites i and j. In this work, the hopping integral ti j is t for nearest neighbors (≥ 0), t′ for diagonal neighbors, and 0 otherwise(Hkin = Ht + Ht′ ) for the two lattices shown in Fig. 1. The
bare energy dispersions are

(
)
′


 −2t (cos k x + cos ky ) − 4t cos k( x cos ky), (a)(2)ϵk = 
 −2t cos k x + cos ky − 2t′ cos k x + ky . (b)In the following, we use t and the lattice spacing as the unitsof energy and length, respectively.
We refer to the former (latter) lattice as a frustrated square(anisotropic triangular) lattice for convenience. The eﬀective
values of t′ /t are considered to be −0.4–−0.1 (∼ 0.3) in holedoped (electron-doped) cuprates.3, 82) For the organic compounds, t′ /t is 0.4–0.8. Hubbard models have been extensively studied, and we have shown that a first-order Mott transition occurs at U = Uc ∼ W at half filling for nonmagnetic
2

J. Phys. Soc. Jpn.

DRAFT

a

(b)

( )

t

with
t'

t

t'
t

Γθ,k =

tSθ,k =

cases and that a doped Mott insulator is realized for U & W
near half filling.33–35) In this paper, we show that similar Mott
physics appears in the SF states. It has been shown that, in awide range of the parameter space of concern, a d-SC state
becomes stable compared with the projected Fermi sea (ordinary normal state).35)

εkk∈kF (t1 /t), σ

A

4θ

eθ
-i

B

A

eθ
i

eθ
i

A

eθ

(8)

(10)

τ

τ

(11)
where d j = n j↑ n j↓ , h j = (1 − n j↑ )(1 − n j↓ ), and τ runs over thenearest-neighbor sites of site j. g, ζd , and ζh are variationalparameters. As shown before,33, 34) the doublon-holon (D-H)
binding eﬀect is crucial for appropriately treating Mott physics. At half filling, ζd and ζh become identical because of theD-H symmetry. In addition to PG and PQ , it is vital for the SF
state to introduce a phase-adjusting factor Pϕ , which we will
explain in the next subsection.

(b)A−2t(cos k x + cos ky ),Here, PG is the fundamental onsite (Gutzwiller) projection∏
PG = j [1 − (1 − g)n j↑ n j↓ ]84) and PQ is an asymmetric projection between a nearest-neighbor doubly occupied site (doublon) and an empty site (holon),35, 85, 86)


∏(
∏(
∏ 
)
)


1 − h j+τ − ζh h j
1 − d j+τ  ,PQ =1 − ζd d jjB=

P = Pϕ PQ PG .

where t1 /t is a band-adjusting variational parameter independent of t′ /t in H and kF (t1 /t) denotes a Fermi surface obtainedby replacing ϵk (t′ ) in Eq. (2) with ϵk (t1 ). It was shown forΨN 38) that the band-renormalization eﬀect through t1 /t owing
to the electron correlation (P) is sizable for U & Uc , a finitet′ /t, and δ ∼ 0.

A

(6)

∆k = ∆θ (cos k x − cos ky ),(9)√and ∆θ = 2t (1 − cos 2θ)/(1 + cos 2θ). Equation (7) is similar to the quasiparticle dispersion of the d-wave BCS wave
function. Note that some important features of the bare ΦSF(summarized in Appendix A) survive in ΨSF = PΦSF . We
do not consider Band-renormalization eﬀects on ΨSF because
those due to diagonal currents or hopping are known to raisethe variational energy for typical cases.83)
The correlation factor P is defined as

2.2 Trial wave functions
We follow many-body variation theory using Jastrow-type
trial wave functions: Ψ = PΦ, where P indicates a product ofmany-body projection (Jastrow) factors discussed later and Φ
is a mean-field-type one-body wave function.
As a normal (paramagnetic) and reference state, we use aprojected Fermi sea, ΨN = PΦN , with∏ΦN (t1 /t) =c†kσ |0⟩,
(3)

-4θ

√
cos2 k x + 2 cos 2θ cos k x cos ky + cos2 ky .

(5)

In H SF , a Peierls phase θ is added to the hopping integrals
so that circular current flows in alternate directions in eachplaquette as shown in Fig. 2(a). In the present variational theory, θ is a variational parameter to be optimized together withthe other parameters. Because H SF breaks time-reversal andlattice-translational symmetries, ΦSF does not have these symmetries. The lower-band energy dispersion of H SF is given as√√
1 + cos 2θ
SF
E− (k) = −2tSθ,k = −ε2k + ∆2k ,
(7)
2
with

Fig. 1. (Color online) Lattice connectivity or hopping paths in extended
square lattice addressed in this study. (a) Frustrated square lattice for cuprate
SCs. (b) Anisotropic triangular lattice for organic κ-ET salt SCs. At lattice
points (solid circles), onsite repulsion U acts.

(a)eiθ cos k x + e−iθ cos ky
,
Sθ,k

B

-i

Fig. 2. (Color online) (a) Schematic figure of local current in staggered
flux state. Arrows denote the directions of currents. (b) Peierls phase factors
attached to hopping terms in H SF (a case of sublattice A); the sign of thephase depends on the relative directions of the hopping and current.

2.3 Configuration-dependent phase factor
A current-carrying state is essentially complex, because thecurrent is proportional to |Ψ|2 ∇Θ, if we represent the state as
Ψ(r) = |Ψ(r)|eiΘ(r) . It is natural that when electron correlation
is introduced, the phase part Θ(r) varies accordingly. However, the conventional correlation factors, PG and PQ , are realand do not modify the phase in ΦSF . Therefore, we need to
introduce an appropriate phase-adjusting factor into the trial
wave function. Such a phase factor was recently introduced
for calculating the Drude and SC weights in strongly correlated regimes;72) thereby, a long-standing problem proposedby Millis and Coppersmith87) —D-H binding wave functions
yield finite (namely incorrect) Drude weights even in the Mott

As a candidate for the pseudogap state, we study a correlated SF state, ΨSF = PΦSF . Here, ΦSF is the one-body
SF state, namely, the ground state of the noninteracting SF
Hamiltonian H SF [shown in Eq. (A·1)], given as∑∏ 1 ∑
ik·ri †
ik·ri † e cBiσ  |0⟩, (4)ΦSF =√  Γθ,k e cAiσ +Ns i∈Ai∈Bk∈kF ,σ3

J. Phys. Soc. Jpn.

DRAFT

parameters, but ΨSF with Pϕ has lower energy than ΨN , aswe will see below.91) Incidentally, this type of phase factor
was also recently shown to be crucial for SF states in a Bose
Hubbard model92) and a d-p model.93) In the regime of Mott
physics, the D-H binding aﬀects not only the real part but also
the phase in the wave function.

insulating regime— was solved.

Pϕ

(a)

B

A
D

B
H

Pϕ

(b)

B

e-iϕ

A

B

(a’)

e-iϕ
e

e

iϕ

iϕ

e0

e-iϕ

A
H

B
D

2.4 Variational Monte Carlo calculations
To estimate variational expectation values, we adopt a plain
VMC method.94–97) In this study, we repeat linear optimization of each variational parameter with the other ones being
fixed, typically for four rounds of iteration. The linear optimization is convenient for obtaining an energy that is discontinuous in some parameters (θ in this case). After convergence, we continue the same processes for more than 16
rounds and estimate the optimized energy by averaging thedata measured in these rounds, excluding excessively scattered data (beyond twice the standard deviation). In each optimization, 2.5×105 samples are collected, so that substantiallyabout 4 × 106 measurements are averaged. Only for ΨSF with
L = 16 and δ = 0, the sample number is reduced to 2.5×104 tosave CPU time. Typical statistical errors are 10−4 t in the totalenergy and 10−4 –2 × 10−3 in the parameters, except near theMott transition points. We use systems of Ns = L×L sites with
L = 10–18 under periodic-antiperiodic boundary conditions.

e-iϕΦSF(θ)

eiθ

B

eiθ

y

x

Fig. 3. (Color online) Illustration for assigning configuration-dependent
phase ϕ in Pϕ . Here, we assume that an electron hops in the x direction.
For the hopping in the y direction, the signs of θ and ϕ have to be reversed.
The (Peierls) phase factor assigned by ΦSF in hopping [Fig. 2(b)] is shown inthe blue dashed box. The values in the red boxes are the phase factors in Pϕcorresponding to the three-site parts shown. The ratio (e±ϕ ) indicated by red
arrows is produced by Pϕ in hopping.

We show that this type of phase factor also plays a vitalrole in the correlated SF state. In ΦSF , a phase θ or −θ isadded when an electron hops to a nearest-neighbor site depending on the direction and position (sublattice), as shownin Fig. 2(b). In the noninteracting case, such hopping occurs
equally in all directions. On the other hand, in the strongly
correlated regime, the probability of hopping depends on thesurrounding configuration (see Fig. 3). For example, whena D-H pair is created [configurations (a) and (a′ )], the next
hopping occurs probably in the direction in which the singly
occupied configuration is recovered [configuration (b)]. Thishopping process does not contribute to a global current in the> Uc ).35) According to a previous study,72) toMott regime (U ∼reduce the energy, it is important to cancel the phase attachedin this type of hopping (±θ) by introducing a phase parameter.
This hopping process can be specified by its local configurations and, correspondingly, we can attach a phase-adjusting
variational factor to the trial wave function. To be more specific, Pϕ gives e−iϕ as shown by the solid boxes (red) in Fig. 3,
with ϕ being a variational parameter. This phase assignment
can be written as[ ∑2∑
Pϕ = exp iϕ (−1)λ+1
dλ, j
λ=1

3. Staggered Flux State at Half FillingFirst, we study the unfrustrated cases (t′ = 0) in order tograsp the global features of the SF state because most of them
do not change even if t′ is introduced. In this section, we focus
on the half-filled case.
3.1 Variational energy
Figure 4(a) shows the variational energy per site of ΨSFmeasured from that of ΨN ,Ẽ = E SF (θ) − E N ,as a function of θ for five values of U/t. Here, the variationalparameters other than θ are optimized for both ΨSF and ΨN .
The size dependence in the case of U/t = 12 is also shown to
see the convergence of the values. For U/t = 6, Ẽ monotonically increases as a function of θ. This behavior is the samefor U/t = 0 shown in Fig. A·2(b) in Appendix A. Hence, ΨSF
is not stabilized for small values of U/t. The situation changes
for U/t > 6; Ẽ/t becomes considerably negative for finite θand has a minimum at θ/π ∼ 0.2 for large values of U/t (= 8–16). This behavior is qualitatively consistent with that of thet-J model, the results of which are summarized in Appendix B
for comparison. In Fig. 4(b), we plot the optimized values of
the configuration-dependent phase factor ϕ as a function of θ.
At the optimized points indicated by arrows, ϕ is very close toθ, especially for large values of U/t. As discussed in Sect. 2.3,
the Peierls phase θ in the hopping process is mostly canceledby ϕ. Although θ is canceled out, the state ΨSF preserves thenature of the original flux state, as shown shortly in Sects. 3.3
and 3.4. That is, a local staggered current flows and the momentum distribution function has a typical k-dependence.
Next, we discuss the U/t dependence of the energy gain of
the fully optimized ΨSF with respect to the reference state ΨN ,j()]× hλ, j+x + hλ, j−x − hλ, j+y − hλ, j−y ,

(13)

(12)

where x and y indicate the lattice vectors in the x and y directions, respectively, λ = 1 (λ = 2) indicates sublattice A(B), and j runs over all the lattice points in sublattice λ. By
Pϕ , a phase factor e±iϕ is assigned to a D-H creation or annihilation process, in which e∓iθ is yielded by ΦSF as shown inthe dashed box (blue) in Fig. 3. Therefore, when the relation
ϕ = θ holds, the total phase shift in a D-H process vanishes.88)
This phase cancelation is acceptable since a phase shift does
not appear in an exchange process in the t-J model. On theother hand, the phase is not canceled in the hopping processes
unrelated to doublons (or of isolated holons).89)
The configuration-dependent phase factor Pϕ is conceptually distinct from position-dependent phase factors used invarious contexts.36, 43, 90) Note that, without Pϕ , the energy of
the SF state is never reduced from that of ΨN for any model∆E(SF) = E(N) − E(SF),4(14)

J. Phys. Soc. Jpn.

DRAFT

0.04

L= 10 U/t

(a)
0.02

t/

ing study (Ref. 35)98, 99) but we adopt Eq. (11) for PQ . For
δ = 0 and U > Uc , Ψd is not SC but Mott insulating. In
Fig. 5, each state exhibits a maximum at U ∼ W (= 8t). The
system-size dependence of ∆E for each state is large near the
maximum but, as shown in the inset of Fig. 5, ∆E remains
finite and the order of the variational energy will not change
as L → ∞. ∆E(AF) is largest, i.e., the AF state has the lowest
energy for any U/t.33) For the SF state, ∆E(SF) ∼ 0 for smallvalues of U/t (. 5). Although ∆E(SF) is always smaller than∆E(d−SC), it is close to ∆E(d−SC). At U/t ∼ 5, ∆E(SF)starts to increase abruptly. The range of U/t where ΨSF is stabilized (U/t & 5) is similar to that of Ψd . In addition, thebehavior of physical quantities such as the momentum distribution function is similar between ΨSF and Ψd as shown
shortly. As mentioned, in the Heisenberg model, ΨSF (g = 0)and Ψd (g = 0) are equivalent due to the SU(2) symmetry, but
in the Hubbard model, the two states are not equivalent, probably due to the diﬀerence in the distribution of doublons andholons.

6
6.8
8
12
16

0

~E
-0.02
-0.04
-0.06

0.4

U/t=12 L

10
12
14

δ=0

t' / t = 0(b)

δ = 0.0
t' / t

φ =θ

=0

0.3
π
/
φ 0.2

=10 U/t
6
6.8
8
12
16

L

0.1
Staggered flux
0

0.1

0.2

θ/π

0.3

0.4

3.2 SF transition and Mott transition
Figure 6 shows the optimized θ and ϕ in ΨSF as a functionof U/t. We find two transition points: USF /t at ∼ 4 – 5 andUc /t at ∼ 7. The former corresponds to the SF transition atwhich ΨSF starts to have finite θ and ϕ and its variational energy becomes lower than that of ΨN . The latter corresponds
to a Mott transition at which the system starts to have a gap in
the charge degree of freedom. The symmetry does not change
at Uc /t.

0.5

Fig. 4. (Color online) (a) Variational energy per site of the staggered flux
state ΨSF (including Pϕ ) measured from that of ΨN [E(θ = 0)] as a functionof θ for several values of U/t (L = 10) in the Hubbard model at half filling.
Data for L = 12 and 14 are added by dashed lines for U/t = 12. (b) Optimized
phase parameter ϕ for the same values of U/t as in (a). The line of ϕ = θ isadded for comparison. The size dependence in (b) is small. In both panels,
the arrows indicate the optimal values of θ when Ẽ/t is minimum.

0.2

-0.3

δ=0

E/t

0.15

0.2
L=18

N
14

-0.4

0.1

U / t = 8.0∆E / t

-0.5
0.005

0.05

1 / L2

0.01

12
14

0

10

20

Uc / t30

U /t

θ /π φ /π L

USF / t

10
12
14

AF
d-waveStag. flux
10

π
φ
, 0.1
π
/
θ

δ=0

t' / t = 0

θ /π φ /π L

0.1

/

d
AF

0

0

SF 10

0

t' / t = 0
20

30

16

5

6

U/t

7

8

U/t
Fig. 6. (Color online) Optimized phase parameters (θ and ϕ) in ΨSF at
half filling. The arrows indicate the Mott transition (Uc /t) and SF transition(USF /t) for four system sizes. The inset shows the same quantities for a wider
range of U/t, with the arrows denoting Uc /t and USF /t for L = 16.

Fig. 5. (Color online) Energy gain of AF, d-SC, and SF states with respect
to projected Fermi sea (ΨN ) at half filling as functions of U/t. Data for L =14, 12, and 10 for each state are plotted as solid lines with symbols, dashdotted lines, and dashed lines, respectively. A guide curve proportional tot/U is drawn for ∆E(SF) with L = 14 (dash-dotted line). In the inset, thesystem-size dependence is shown for U/t = 8.0 and fitted by second-orderpolynomials.

At USF /t, θ and ϕ exhibit first-order-transition-like discontinuities, for example, at USF /t = 6.28 for L = 10. However,as L increases, USF /t shifts to lower values and the discontinuities become small and unclear, suggesting that the SF
transition is continuous and occurs at a small U/t. Because
an appropriate scaling function is not known, we simply perform a polynomial fit of USF /t up to the square of 1/L2 as arough estimate. This yields USF /t = 2.93 for L = ∞ with asmall error. In any case, since θ and ϕ are tiny for U/t . 5,
we consider that ΨSF is substantially not stable in a weakly
correlated regime.

where E(N) and E(SF) are the optimized (including θ) energies per site of ΨN and ΨSF , respectively. If ∆E(SF) is positive, the SF state is stabilized with respect to ΨN . In Fig. 5,
we show ∆E(SF) compared with other ordered states, i.e., theAF state, ΨAF = PΦAF , and the d-SC (projected BCS) state,Ψd = PΦd . We use the same ΦAF and Φd as in the preced5

J. Phys. Soc. Jpn.

DRAFT

1
0.8

Staggered flux
δ = 0.0
t' / t = 0

d

0.6

ζd

ζh

L
10
12

ζ

d

14
16

0.4

Uc/t

0.2
0

On the basis of the single-mode approximation,103–105) excitations in the charge sector are gapless when N(q) ∝ |q|
for |q| → 0, whereas a gap opens in the charge sector whenN(q) ∝ |q|2 . For S (q), a similar relation holds for the spin
sector, although excitations cannot be sharply divided into thecharge and spin sectors except for in one-dimensional systems. In Figs. 8(b) and 8(c), N(q) and S (q) are respectively
shown for various values of U/t. For U < USF , both N(q)and S (q) behave linearly for |q| → 0 as expected for a Fermi
liquid. For U > Uc , the behaviors of both N(q) and S (q) appear to be quadratic and consistent with a Mott insulator. For
USF < U < Uc , N(q) is linear near (0, 0), whereas S (q) is
quadratic-like,106) indicating that the state is a spin-gap metal.
Namely, the charge and spin sectors show diﬀerent tendencies
in excitation. This feature is distinct from that of the noninteracting SF state ΦSF , which has a gap common to both sectors,
and thus N(q) = S (q) holds (see Appendix A). Therefore, themetallic SF state stable for USF < U < Uc is not perturbatively connected to ΦSF . We will show in the next section that
this state is connected to the metallic SF state in the doped
case.

0.2

0.1

d
5

6

7

U/t

8

9

0
10

Fig. 7. (Color online) U/t dependences of optimized D-H binding parameter ζ (= ζd = ζh ) in ΨSF and density of doublons d shown at half filling forsome system sizes. The arrow indicates the Mott transition point.

Next, we examine the Mott transition at Uc /t, where thebehaviors of θ and ϕ change as shown in Fig. 6. Uc /t = 7.1
for L = 16.100–102) (Note that in the t-J model, the Mott transition cannot be discussed.) To confirm that Uc /t is a Mott
transition, we plot the U/t-dependences of the optimized DH binding parameter ζ (= ζd = ζh ) and the doublon density
d = EU /U in Fig. 7. These quantities are sensitive indicatorsof Mott transitions. In Fig. 7, we find abrupt changes in both ζand d at Uc /t, similarly to those in the Mott transitions in Ψdand ΨN .34) In ΨSF , discontinuities in ζ and d at Uc /t are notfound even for the largest system we treat (L = 16). However,
because the behavior of both ζ and d becomes more singularas L increases, we consider that this transition is first-order,
similarly to those in Ψd and ΨN .34)

3.4 Circular current
Now, we turn to the local circular current in a plaquettedefined as
1 ∑ ∑JC /t =(−1)ℓ̃+τ̃ Im⟨c†ℓ+τ,σ cℓ,σ − c†ℓ,σ cℓ+τ,σ ⟩, (18)
Ns ℓ∈A,σ τ
where ℓ runs over all the A sublattice sites, ℓ̃ = ℓ x + ℓy , τ
indicates the nearest-neighbor directions, and τ̃ = 1 [−1] forτ = (±1, 0) [(0, ±1)]. JC is regarded as the order parameterof the SF phase. In the main panel of Fig. 9, we show |JC |/t
at half filling as a function of U/t. In the metallic SF phase(USF < U < Uc ), a relatively large current flows. In the insulating SF phase (U > Uc ), the local current is reduced but still
finite. At U/t = 12, however, |JC |/t is 1/20 of that in ΦSF withthe same θ. The U/t-dependence of JC /t in this regime is fitted
by a curve proportional to (t/U)2 and the system-size dependence is small, as shown in Fig. 9. This suggests that JC in
this range of U/t has a localized nature. More specifically, |Jc |will be related to the local four-site ring exchange interaction,
which appears in the fourth-order perturbation with respect tot/U in the large-U expansion of the Hubbard model.

3.3 Spin-gap metal
In the intermediate regime USF < U < Uc , the present SF
state is expected to be metallic. In order to clarify the natureof ΨSF , we calculate the momentum distribution function
1∑ †n(k) =⟨c ckσ ⟩(15)
2 σ kσfor the optimized ΨSF . Figure 8(a) shows n(k) along the path(0, 0)-(π, 0)-(π, π)-(0, 0) in the original Brillouin zone for various values of U/t. In the region of U < USF (i.e., θ = 0),
we find two discontinuities (crossings of the Fermi surface)at k ∼ (π, 0) and (π/2, π/2), indicating a typical Fermi liquid. For USF < U < Uc (half-solid symbols), the discontinuity at k ∼ (π, 0) disappears, while the discontinuity atk ∼ (π/2, π/2) remains. This is qualitatively identical to thatof the noninteracting SF state ΦSF , in which there is a Diracpoint at (π/2, π/2) and a certain gap opens near the antinodal
points. On the other hand, for U > Uc (solid symbols), both
discontinuities disappear, indicating that a gap opens in thewhole Brillouin zone. This is consistent with a Mott insulator.
We can reveal the characters of the gaps to some extent byanalyzing the charge density and spin structure factors,⟩
1 ∑ iq·(Ri −R j ) ⟨eni n j − n2 ,(16)N(q) =
Ns i, jS (q)=1 ∑ iq·(Ri −R j ) ⟨ z z ⟩e
SiS j .
Ns i j4. Staggered Flux State at Finite Doping
4.1 Energy gain and optimized phase parameters
First, we show the energy gain of ΨSF with respect to thereference state ΨN [Eq. (14)] in Fig. 10(a) for four values ofthe doping rate δ. Similarly to the half-filled case (Fig. 5),∆E is zero for the weakly correlated regime (U < USF ); thevalue of the SF transition, USF /t, increases as δ increases. The
sharp peak of ∆E for δ = 0 changes to a broader peak with amaximum at U/t ∼ 12 – 16, and finally vanishes at δ ≡ δSF ∼0.1.
In Fig. 10(b), optimized values of the phase parameters areplotted. Both parameters decrease as δ increases. Although θand ϕ have a discontinuity at USF /t at this system size, thisbehavior is owing to a finite-size eﬀect.107) The SF transition
for L → ∞ will be continuous, similarly to the half-filledcase. When we compare with the results at δ = 0, we see that

(17)

6

J. Phys. Soc. Jpn.

DRAFT

1

0.6

0.8
0.6

)
(

k
n

U/t

(a)

U/t

3.0
5.0
6.0
6.5
6.9
7.2
12
16

0.5

5.0
6.0

1

(b)

δ =0
t' / t = 0
L = 14

4.0

0.8

6.5
6.9

0.4

/

Uc t

12

)
(

q



/

USF t

7.2

0.6

16

q

6.0

7.2
12

)
( 0.3

=0
L = 14t' / t

5.0

6.9

(c)

δ =0

4.0

6.5

0.4

U/t

3.0

～～ 5.05
7.1

16

S

0.4

0.4

L

U/t= 4.0

10
12

0.3

0.2

14
16

L

0.2

0.2

δ = 0.0

0.1

=0
= 14

t' / t
L

0
Γ
(0,0)

0.2

X
(π,0)

M
(π,π )

k

Γ

(0,0)

10
12

0.1

14

U/t= 6.5

0
(0,0)

(π,0)

q

(π,π )

0
(0,0)

(0,0)

16

q

0
(0,0)

(π,0)

π

( /2,0)(π,π )

q(0,0)

Fig. 8. (Color online) Correlation functions in momentum space of ΨSF at half filling for various values of U/t. (a) Momentum distribution function, (b)charge density structure factor, and (c) spin structure factor. In the inset in (c), the system-size dependence of S (q) (L = 10–16) for typical cases of a Fermi
liquid (U/t = 4) and spin-gap metal (U/t = 6.5) is shown for small |q| in the direction of (0, 0)–(π, 0). See also Ref. 106. Open (black) symbols are for U < USF(i.e., θ = 0), half-solid (brown) symbols for USF < U < Uc , and solid (red) symbols for U > Uc . For this system size (L = 14), USF /t ∼ 5.05 and Uc /t ∼ 7.1.

0.5

L
10

0.4

12

0.06

δ

L = 12

0.4

0.0

0.3

0.0556

/
|

0.1

0

0.1

0.02
10

U /t

USF /t20

30∝(t/U)

2

Uc/t

0

10

0.0
0.02780.05560.0833

t/
E∆

0.2

0.2

L=12 δ

0.04

0.0833

δ=0
t 0.3
JC

10
12
14

0.0278

14
16

|

δ=0 L

(a)

U/t

20

0

0.3

30

θ/π φ/π δ

L = 120.0
0.02780.0556

Fig. 9. (Color online) Absolute values of local circular current at half filling as a function of U/t for some system sizes. The SF and Mott transition
points are shown by arrows. A curve proportional to (t/U)2 is shown by agray dash-dotted line. The inset shows the same quantity for some doping
rates for L = 12 (discussed in Sect. 4).

π
/
φ
,π
/
θ

0.0833

0.2

0.1
(b)

ΦSF is realized in the strongly correlated region (U > W), and
it is smoothly connected to the Mott insulating state at halffilling. It is also interesting that ϕ becomes larger than θ as δ
increases, while they are close to each other when δ = 0. Thissuggests that ϕ overscreens the phase θ in the D-H processes
owing to the increasing number of free-holon processes.
Figure 11 shows the δ-dependence of ∆E/t for the casewith U/t = 16. Except for the case with δ = 0, ∆E/t monotonically decreases as a function of δ. Because the L dependence is appreciable, δSF should be somewhat larger in the
L → ∞ limit. The behavior of ∆E is consistent with that for
the t-J model shown in Appendix B.108) In the inset of Fig. 11,
the δ-dependences of the optimized θ and ϕ are plotted. Theirsystem-size dependences are very small.

0

USF/t

Staggered flux
10

U/t

20

30

Fig. 10. (Color online) (a) Energy gain of SF state with respect to ΨN as
a function of U/t for four values of doping rate δ (L = 12). For δ = 0, data
for three system sizes are shown. (b) Optimized phase parameters θ and ϕ. In
both panels, the SF transition is indicated by arrows.

a discontinuity near (π/2, π/2), indicating that ΨSF is always
metallic and there is no Mott transition. Furthermore, N(q)is linear in |q| for |q| → 0, indicating that the charge degreeof freedom is gapless. On the other hand, S (q) appears to beapproximately quadratic at small |q| for U > USF , suggesting that the SF state in the doped region has a gap in the spin
sector.
(ii) Segmented Fermi surface: The bare SF state, ΦSF , has asegmented (or small) Fermi surface around k = (π/2, π/2) asshown in Appendix A (see Fig. A·4). Here we show that thisfeature is preserved for strongly correlated cases. Shown in
Fig. 13 are contour maps of n(k) for U/t = 12 and four values

4.2 Various properties(i) Spin-gap metal: In Fig. 12, we show the behavior ofcorrelation functions in the momentum space n(k), N(q), andS (q) for δ = 0.0556 (L = 12). The U/t-dependences ofthese quantities are basically similar to those at half filling
discussed in Fig. 8. In the region of U > USF , n(k) preserves
7

J. Phys. Soc. Jpn.

DRAFT

1

4.0
6.0
7.0
7.75
8.0
10
12
16
USF / t = 7.92

0.8
0.6
)
(

k
n

1

0.6

U/t

U/t

δ = 0.0556
L = 12t' / t = 0

4.0
6.0
7
7.75
8
10
12
16

0.5
0.4
)
(

q



(c)

(b)
0.8

0.6
)
(

q

0.3

S

0.4

U/t

4.0
6.0
7
7.75
8
10
12
16

0.4

0.2

δ = 0.0556

L = 12t' / t = 0

0.2

δ = 0.0556

L = 12t' / t = 0

0.2

0.1

(a)

USF / t = 7.92

USF / t = 7.92

0
Γ
(0,0)

X
(π,0)

M
(π,π)

k

Γ

(0,0)

0
(0,0)

(π,0)

(π,π)

q

(0,0)

0
(0,0)

(π,0)

q

(π,π)

(0,0)

Fig. 12. (Color online) Behavior of (a) momentum distribution function, (b) charge density structure factor, and (c) spin structure factor of ΨSF for a finite
doping (δ = 0.0556) and various values of U/t. Open (black) symbols are for U < USF and half-solid (brown) symbols are for U > USF . For this system size
(L = 12) and doping, the SF transition is at USF /t = 7.92.

0.05

U / t = 160.04

t/ 0.03E∆(a) į = 0.0

0.05

䢲䢰䢵
䢲䢰䢶

0

δ

δ

(b) į ~ 0.04

䢲䢰䢷

φ /π

θ

0.02 ∆E/t L
10
12
14
0.01
16
Staggered flux
0

U / t = 16

0.2
π
φ
,
π
/
θ0.1
/

䢲䢰䢵

0.3

䢲䢰䢶

䢲䢰䢷

/π
0.05

0.1
θ /π

0.1

φ /π

0.7

䢲䢰䢹

0.5

L

10
12
14
16

䢲䢰䢶
䢲䢰䢵

0.15

䢲䢰䢷

䢲䢰䢸

䢲䢰䢶
䢲䢰䢵
䢲䢰䢷

(c) į ~ 0.08

(ʌ, ʌ)(d) į ~ 0.12
0.25

Fig. 11. (Color online) Energy gain of SF state with respect to ΨN as afunction of doping rate at U/t = 16. Data for four system sizes are shown. In
the inset, the optimized phase parameters θ and ϕ are shown.

0.3

0.25

(ʌ, 0)(0, 0)of δ. Here, we show the data for t′ /t = −0.3 because ΨSF is
stabilized in a wide doping range (see Fig. 15 later) and the
behavior is similar to that for t′ /t = 0. At half filling, there is
no Fermi surface, as shown in panel (a). Upon doping, however, pocket Fermi surfaces appear around (π/2, π/2) and, as δ
increases, they extend to the antinodes along the AF Brillouinzone edge. These Fermi surfaces are shown by blue dashedovals in panels (b)-(d). A gap remains open near (π, 0).
(iii) Circular currents: The local circular currents JC defined in Eq. (18) for the doped cases have already been shownin the inset of Fig. 9, where the evolution of JC with increasing δ is shown as a function of U/t. We find that JC increases
as δ increases, although the optimized phase parameters θ andϕ decrease [see Fig. 10(b)]. This is probably because the number of mobile carriers increases as δ increases in the strongly
correlated regime, whose feature is typical of a doped Mott
insulator. In contrast, as shown in Appendix A, JC decreases
as δ increases in the noninteracting ΦSF . At the phase transition point δSF , where E(SF) becomes equal to E(N), the order parameter |JC |/t drops suddenly from 0.25–0.3 (almost the
maximum value) to zero. This indicates that this transition is
first-order, in contrast to the corresponding AF and d-SC transitions, as a function of δ.

0.7 0.65

0.75

0.7

Q = (ʌ, ʌ)

Fig. 13. (Color online) Contour maps of momentum distribution functionn(k) of the optimized ΨSF in the original Brillouin zone for t′ /t = −0.3,
U/t = 12, L = 10 – 14, and four values of δ in (a)-(d). The wiggles of lines
are simply due to the small number of k points and are not important. The AF
Brillouin zone boundary is indicated by pink dotted lines, and zone-diagonallines are shown with gray dotted lines in (c). Fermi surfaces are indicatedwith blue ovals in the third quadrants. In (d), the scattering vector of q = Q
connecting the antinodes, discussed in Sect. 6.4, is shown with a blue arrow.

5.

Eﬀect of Diagonal Hopping t ′In this section, we study the eﬀect of diagonal hopping t′ inthe two cases shown in Fig. 1.
5.1 Frustrated square lattice
Figure 14 summarizes the total energies of ΨSF , ΨN , and Ψd
as functions of t′ /t. Note that the energy for ΨN without band
renormalization exhibits complicated behaviors as a functionof t′ /t. This is because the occupied k-points in the Fermi surface change discontinuously in ΦN . However, if we consider
the band-renormalization eﬀect38) for ΨN and use the optimized t1 /t, the lowest energy for ΨN becomes the black solid
8

J. Phys. Soc. Jpn.

DRAFT

(a)
-0.2 L=12 Normal
L=10
NormalStag. flux
d-waveδ = 0.0
U / t = 12

-0.3

0

NormalNormal (opt.)
Staggered flux
d-waveE/t

E/t

U / t = 12
L = 12-0.3

Normal (opt.)
Stag. flux
d-waved-wave (opt.)-0.5

δ = 0.0278

(b)

-0.35

0.5

-0.5

0

t' / t0.5

t' / t
-0.55
(d)

(c)
δ = 0.0833

U / t = 8.0
L = 12

t'SC
-0.56

t'c1

-0.5

-0.55

-0.6

U / t = 12
L = 12t'c2

δ = 0.0833

E/t

E/ t

-0.45

-0.65

NormalNormal (opt.)
Staggered flux
d-waved-wave (opt.)

-0.57

t'c1

)UXVWUDWHGVTXDUH
-0.75
-0.5

t'c2

-0.7 -0.58

0

-0.59
-0.5

-0.4

-0.3

-0.2

-0.5

0.5

-0.1

0

0

0.5

t' / tt' / tFig. 14. (Color online) Comparison of total energies among ΨN , ΨSF , and Ψd as functions of t′ /t in frustrated square lattice of L = 12. Panels (a)-(c) display
the cases of diﬀerent doping rates for U/t = 12 and (d) shows the case of U/t = 8 and δ = 0.0833. In (a), data for L = 10 are also shown. Symbols are common
to all panels. The black and red bold lines indicate the values of ΨN and Ψd , respectively, when the band-renormalization eﬀect is considered (for details, see
Ref. 38). The arrows indicate the positions of the energy crossings. The inset in (d) shows a magnification of the area of energy crossings.

14

line in Fig. 14. We use the solid lines as energies for ΨN . Although we expect some size eﬀects in ΨN , we can see generaltrends of the energy diﬀerences between ΨSF , ΨN , and Ψd .
At half filling [Fig. 14(a)], E is symmetric with respect tot′ /t = 0 owing to the electron-hole symmetry.35) E(SF) is always lower than E(N) and does not depend on t′ /t because⟨Ht′ ⟩ = 0 for any t′ /t and U/t. E(d−SC) tends to increaseas |t′ /t| increases. (When band renormalization is taken into
account, E(d−SC) also becomes constant.38, 81) )In a slightly doped case [Fig. 14(b)], E for every state becomes a decreasing function of t′ /t, but the order of the energies does not change. E(SF) slightly depends on t′ and remains a linear function of t′ . However, for large δ, the situation changes as shown in Fig. 14(c). The range of E(SF) <′′E(N) is restricted to tc2< t′ < tc1, as indicated by arrows. This
range becomes smaller when U/t decreases [Fig. 14(d)]. We
also find that this stable range of ΨSF becomes smaller as δ increases and finally vanishes at δSF ∼ 0.16 (0.12) for U/t = 12
and 16 (8).
Obtaining similar data for diﬀerent values of δ, t′ /t, and L,
we construct a phase diagram, as shown in Fig. 15(a). The
stable area of the SF state expands above the optimum dopingof cuprates (δ ∼ 0.15) for −0.4 . t′ /t . −0.1, which corresponds to the hole-doped cuprates. As L increases, the area ofΨSF tends to expand slightly. For t′ /t > 0, on the other hand,the area of ΨSF shrinks to a very close vicinity of half filling,
especially for U/t = 8.
It is useful to draw a phase diagram in the δ-U plane. Figure

12

L
10

Staggered
flux

12

14

t' / t
-0.3
0.0
0.2

t/10
U
8

Ordinary normal

6
Frustrated square

4
0

0.05

0.1

0.15

δ
Fig. 16. (Color online) Phase diagram between SF and normal states in δU plane for three values of t′ /t. The area of the SF phase may be somewhat
smaller for δ ∼ 0 because the band renormalization in ΨN is not consideredhere.

16 shows the region in which the SF state is stabilized forthe cases with t′ /t = −0.3, 0, and 0.2. Irrespective of t′ /t,
the boundary value, USF /t, increases as δ increases.109, 110) Fort′ /t = −0.3 (corresponding to the hole-doped case), the SF isstable in the whole underdoped regime (δ . 0.16) for U/t &10. In contrast, for t′ /t = 0.2 (electron-doped case), USF /t
rapidly increases as δ increases.
This stability of ΨSF in a wide range of δ for −0.4 . t′ /t .
−0.1 originates primarily from the large t′ /t dependence of9

J. Phys. Soc. Jpn.
0.5

t

0

'
t

-0.5

1

de
po
d
no
rt
ce
le
de
po
d
leo
h

/

DRAFT

Frustrated square

Ordinary
normalκt/'t

Staggered flux
Ordinary
normal

(a)

0

0.05

δ

0.1

Aniso. triangular

tsl
sa
TE
-

L=10 12 14 U / t
8
12
16

0

-1
0

0.15

Ordinary
normalStaggeredflux

L=10 12 U / t
8
12

(b)
0.05

0.1

δ

Fig. 15. (Color online) Phase diagrams in the δ-t′ plane with doping rate (δ) and frustration strength (t′ /t) for (a) the frustrated square lattice and (b) the
anisotropic triangular lattice for U/t = 8 and 12 [and 16 in (a)]. “Ordinary normal” indicates the projected Fermi sea. The symbols with vertical bars placed att′ /t ∼ −0.25 indicate that the SF state is not stabilized at these values of δ. The scale of the abscissa is identical in the two panels.

ΨN and the very small t′ /t dependence of ΨSF . This meansthat the nature of ΨSF for t′ = 0 quantitatively remains that fort′ /t , 0. For example, we show in Fig. 17 the t′ /t dependencesof the optimized phase parameters and local circular current,
JC , which do not strongly depend on t′ /t. Furthermore, we
confirm that the momentum distribution function n(k) is almost the same for t′ /t , 0. Note that this is in sharp contrast to
the d-SC state, in which n(k) in the antinodal region markedly
changes with t′ /t (see Fig. 29 in Ref. 35). The reason for thisdiﬀerence between ΨSF and Ψd will be as follows. Since ΨSFis very appropriately defined for the simple square lattice, t′change the wave function of ΨSF only slightly. On the other
hand, Ψd has a gap opening at the Fermi surface near (π, 0),
which is markedly aﬀected by t′ . In this context, it is natural
to expect that extra current, such as diagonal currents in chiralspin states,111) will not be favored.83)

0.4

χc-1

10

|JC|/ tθ /π φ /π

0.3

Stag. flux

0

/
|

J
|

0.1

θ /π

0.1

π
/
φ
,π
/
θ

metastable

χ−1
c =

L = 12

-0.6

-0.4

-0.2

0

0.2

0.2

case. To this end, we consider the charge compressibility κ orequivalently the charge susceptibility χc (= n2 κ), the inverseof which is given as

U / t = 1200.1

Fig. 18. (Color online) Inverse charge susceptibility of the SF state shown
as a function of doping rate δ for some parameter sets for the frustrated squarelattice [Fig. 1(a)]. We add guide lines (thick dashed) for t′ /t = 0 and −0.3
(U/t = 12, L = 14). The same quantity of the d-SC state is also shown withhalf-solid green symbols for comparison. Zigzags of the data for ΨSF are dueto the discontinuous change in the occupied k points and other finite-size
eﬀects.

0.08330.1111

φ /π

C0.2

U/t, t'/t, L
12 0.3 12
16 0 14
12 0 12
12 0 14
12 -0.3 10
12 -0.3 12
12 -0.3 14

δδ

|JC| / t

t

0

-10

0.2

Frustrated squared-waveU/t, t'/t, L
12 0 12
12 -0.3 12

0∂2 E(n)E(N + 4) + E(N − 4) − 2E(N)
∼ Ns2
,
2
∂n
42

(19)

with n = N/Ns . If χ−1
c < 0, the system is unstable toward
phase separation. In Fig. 18, we show the δ dependence of′−1
χ−1
c for three values of t /t and L. We find that χc is basically′negative for t /t = 0 and 0.3, indicating that ΨSF is unstabletoward phase separation (data points of δ ∼ 0 should be disregarded because they are aﬀected by the Mott singularity at
δ → 0). This result is consistent with the previous one for thet-J model for t′ /t = 0.61) In contrast, for t′ /t = −0.3, χ−1c becomes positive for δ & 0.05 and comparable to that of d-SC
(green symbols). Therefore, a homogeneous SF state is possible for the parameters of hole-doped cuprates.

t' / tFig. 17. (Color online) t′ /t dependences of the optimized phase parametersθ and ϕ and the local circular current JC in ΨSF for two values of δ. The
unstable regions of ΨSF are indicated by gray dashed lines and “metastable”.

In many studies on the t-J and Hubbard models, the instability toward phase separation near half filling has been discussed. Recently, states with AF long-range orders have beenshown to be unstable toward phase separation for t′ /t ∼ 0 for
the Hubbard model using the VMC method.35, 38, 112–114) For
the t-J model, an SF state has also been shown to be unstabletoward phase separation in a wide range of δ for t′ /t = 0.61)
Therefore, we need to check this instability in the present
10J. Phys. Soc. Jpn.

DRAFT

USF /t between ΨSF and ΨN increases as t′ /t increases.
As shown before, the Mott transition occurs at Uc ∼ 7.1t
for t′ = 0. Since the properties of ΨSF are similar to the caseof t′ = 0, the SF state is metallic for U . Uc and insulating
for U > Uc . This phase boundary is also shown in Fig. 20.

5.2 Anisotropic triangular lattice
For the anisotropic triangular lattice, Fig. 19 summarizes
the total energies of ΨSF , ΨN , and Ψd as a function of t′ /t.
Again, E(N), the lowest energy of ΨN considering band renormalization,38) is shown by the black solid lines. First, let
us consider the half-filled case [Fig. 19(a)]. E is symmetricwith respect to t′ /t = 0 owing to the electron-hole symmetry.35) Similarly to the case on the frustrated square lattice,
E(FS) and E(d−SC) (if the band renormalization is considered115, 116) ) are constant. For the ground state, we find from′
Fig. 19(a) that ΨN has a lower energy than Ψd for t′ > tSC
with′tSC /t = 0.807 for L = 12. However, the (π, π)-AF state or incommensurate AF states including the case of a 120◦ structure
has a lower energy at half filling.115–117)

12$QLVRWULDQJXODU

δ = 0.0

t1 = 0

10

ӆ0

U/t

Staggered flux
8

t1
Uc/t

6

t'c0

Ordinary normal
4(a)δ = 0.00

U / t = 8.0

E/t

L=10 12t'c2d (renorm.)$QLVRWURSLFWULDQJXODU
-0.45
0

0.2

For a doped case, we show in Fig. 19(b) the t′ /t dependenceof the total energy for the three states for typical parameters.
It is noteworthy that E/t for ΨN and Ψd decreases rapidly forlarge values of |t′ /t| (∼ 1). Obtaining similar data for various
values of t′ /t and δ, we construct a phase diagram in the δ-t′
space [Fig. 15(b)]. Compared with the case of the frustrated
square lattice [Fig. 15(a)], the area of ΨSF is restricted to thesmall doping region.

t'SC

0.4

0.6

0.8

t' / t-0.36
NormalStag. flux
d-wave

(b)

E/t

-0.38

6. Discussion

-0.4

t'c2

-0.42
-1

1.5

Fig. 20. (Color online) Phase diagram between the SF and normal states int′ -U plane at half filling on anisotropic triangular lattice. Boundaries (USF ,solid lines) are determined for three system sizes. The green arrow near thevertical axis represents the range of the metallic SF state. The dashed lines
show the boundary in ΨN regarding whether the nesting condition is restored(t1 = 0) or not (t1 , 0) in the renormalized band (see Sect. 6.2 later).

t'c0NormalStag. flux
d-wave-0.4

1t' / t-0.35
Normal(opt.)

0.5

t'SF L
1012
14

6.1 Phase cancelation mechanism
First, let us consider intuitively why ΨSF has a low energy
in the strongly correlated region of the Hubbard model. As
discussed in a previous study,35) the processes corresponding
to the J term in the t-J model are those in which a D-H pairis created or annihilated as shown in Fig. 21(a). Generally
speaking, the phase yielded in this process causes a loss ofkinetic energy. In order to reduce this kinetic energy loss, thephase θ should be eliminated by the phase ϕ by applying Pθ
with ϕ ∼ θ in the same manner as introduced in Ref. 72.
We expect a similar phenomenon in the Heisenberg interaction in the t-J model. In the J term, an ↑ spin at site i hops tosite j (= i+τ) and simultaneously a ↓ spin at site j hops to site
i. As shown in Fig. 21(b), if the former hopping yields a phaseθ, the latter yields −θ in ΦSF ; the total phase in the exchange
process precisely cancels out (shown in the square brackets in
Fig. 21). Since the two processes occur simultaneously, it isunnecessary to introduce ϕ in the t-J model to stabilize the SF
state. On the other hand, in the Hubbard model [Fig. 21(a)],a hopping resulting in D-H-pair annihilation does not necessarily occur immediately after a D-H pair is created; thesetwo processes are mutually independent events. Therefore, itis necessary to introduce the phase ϕ to eliminate ±θ in eachprocess in order to stabilize the SF state.

t'c1

δ = 0.0556
U / t = 12
L = 12
-0.5

0

0.5

t' / tFig. 19. (Color online) Comparison of total energy among ΨN , ΨSF , andΨd as functions of t′ /t. Panel (a) displays the half-filled case for U/t = 8(only the data for t′ /t ≥ 0 are shown) and L = 10 and 12, and (b) a doped casefor U/t = 12 (δ = 0.0556) and L = 12. We show the optimized values by theband renormalization for ΨN with (black) solid lines. The band-renormalized
data for Ψd are also shown in (a) with a (red) solid line.

We compare the energies between ΨN and ΨSF , assumingthat AF states are not stabilized. In Fig. 19(a), ΨSF is more′′stable than ΨN for t′ < tc2
with tc2/t = 0.749 (0.763) for
U/t = 8 and L = 12 (10). We employ the value U/t = 8 simply because it is frequently used as a plausible value for κ-ET
salts.73) Actually, ΨSF is Mott insulating at U/t = 8; a value of
U/t . 7.1 is necessary for a metallic state. However, the point
here does not change qualitatively irrespective of whether ΨSF
is insulating or metallic. On the basis of similar calculations
for various values of U/t and t′ /t, we construct a phase diagram within the SF and normal state at half filling (Fig. 20),
which is relevant for organic conductors. The boundary value
11

J. Phys. Soc. Jpn.

DRAFT(a) ΨSF (Hubbard)
D-H process~(b) ΨSF (t - J)
exhange~(c) ΨSSC (t - J)
exchangeeiθeiθei(θ-ϕ)[~ ei0]ij[= e 0]

0.1
0

t/
-0.1
E
∆

= e2 ]

i

(a)

[

t

iθ

-0.2

e-i(θ-ϕ)

e-iθ

δ

=0

δ

=0

eiθ
-0.3

Fig. 21. (Color online) Illustration of phase factors added in (a) creationor annihilation process of doublon-holon (D-H) pair in ΨSF for large-U/t
Hubbard model, (b) spin exchange process in Ψ̃SF for t-J model, and (c) spin
exchange process in Ψ̃SSC for t-J model. For details, see text.

0.3

(b)

0.2

Et (N) − Et (SF),EU (N) − EU (SF),==

Et (SF) − Et (d−SC),
EU (SF) − EU (d−SC).

0.0
0.02780.05560.0833

L=12

δ

0.0
0.02780.05560.0833

-0.1
0

10

U/t

20

30

Fig. 22. (Color online) The two components of the energy diﬀerence between the SF state and the projected Fermi sea are shown for some doping
rates (t′ /t = 0, U/t = 12): (a) kinetic energy and (b) interaction energy parts.
For δ = 0, we add data for L = 10 and 14. The arrows indicate USF /t for
L = 12.

kinetic energy gain for a suﬃciently large U/t.119)
6.3 Comparison with experiments(i) High-T c cuprates: Here, we discuss the lattice translational symmetry, which is broken in the present SF state.
The peaks arising from local loop currents in the polarized
neutron scattering spectra are found at k = (0, 0),17–19) suggesting that the lattice translational symmetry is preservedin the pseudogap phase. Some authors have argued that the
SF state breaks this symmetry, but physical quantities calculated with SF states display a (0, 0) peak in addition to a (π, π)
peak.120) The above neutron experiments appear to be consistent with more complicated circular-current states that do notbreak this symmetry.90, 121, 122) Recently, however, one of theauthors showed that this type of circular-current state is notstabilized with respect to the normal state in a wide range ofthe model parameters on the basis of systematic VMC calculations with refined wave functions for d-p-type models.
Instead, SF states are stabilized in some parameter ranges.123)
On the other hand, the shadow bands observed in the ARPES
spectra,24–27) which also characterize the pseudogap phase,
seem to require the scattering of (π, π) and a folded Brillouinzone. Therefore, the issue of translational symmetry breaking
is still controversial.
(ii) Organic conductors: In Sect. 5.2, we studied theanisotropic triangular lattice. Let us here discuss the relevance of the present results to experiments. As discussed in
Sect. 1, deuterated κ-(ET)2 Cu[N(CN)2 ]Br has t′ /t ∼ 0.4.
The present results show that the SF state is stabilized forthe case of t′ /t ∼ 0.4. Therefore, the pseudogap behavior
for T > T c observed in deuterated κ-(ET)2 Cu[N(CN)2 ]Bris probably caused by the SF state. On the other hand, κ-(20)

for four values of δ. In previous papers, we perfomed the same
analysis for Ψd and ΨN ,33, 34) and showed that the SC transition is driven by the kinetic energy gain for U & Uco with
Uco /t being the crossover value from weakly to strongly correlated regimes. In Fig. 22, we find that a similar phenomenon
emerges between ΨSF and ΨN : Kinetic energy gain occurs inthe strongly correlated region. The physical reason for thiswill be as follows. In the strongly correlated regime, the kinetic energy is dominated by the D-H pair creation or annihilation processes (not shown). Since the phases arising in theseprocesses are canceled out by ϕ, this kinetic energy gain corresponds to that in the J-term in the t-J model.
In Fig. 23. we show a similar comparison between the d-SC
and optimal SF states, i.e.,∆Et∆EU

L=12 δ

0

6.2 Kinetic energy gain
We discuss another physical reason for the stabilization of
the SF state. In Fig. 22. we show the diﬀerence in the kinetic
energy ∆Et and interaction energy ∆EU between the optimal
SF state and the projected Fermi sea,
=
=

L

10
12
14

t/
E∆U 0.1

Finally, let us apply the present mechanism to the spincurrent-carrying state. Staggered spin current (SSC) states (orsometimes called spin-nematic states63) ) have been considered to be candidates for hidden orders in various systems.118)
In these states, counter-rotating currents of ↑ and ↓ spins alternately flow in each plaquette [Fig. A·1(b)]. We have carried out similar VMC calculations for the projected SSC stateΨ̃SSC = PG (0)ΦSSC . The results are summarized in Appendix
C. We conclude that Ψ̃SSC is not stabilized for any J/t and underdoped δ. We can easily see the reason for this by considering the phase cancelation. As we can see from Fig. 21(c),the total phase added in an exchange process in Ψ̃SSC remains 2θ. We found that this phase is diﬃcult to eliminate
by configuration-dependent phase factors such as Pϕ . Therefore, we conclude that the SSC state or spin-nematic state will
never be stabilized.

∆Et
∆EU

L

10
12
14

(21)

In particular, in the regime of U > Uc (SF) at half filling andU > USF for δ > 0, the energy gain occurs exclusively inthe kinetic part (∆Et > 0 and ∆EU < 0). Thus, the causeof stabilization both in ΨN → ΨSF and in ΨSF → Ψd is the12

J. Phys. Soc. Jpn.

0

DRAFT

It revealed that the AF state is considerably stable and occupies a wide range of the ground-state phase diagram. In doped
metallic cases for t′ /t . −0.05, an AF state called type-(ii) AF
state is stabilized, while for t′ /t & −0.05, a type-(i) AF stateis stabilized. In a type-(ii) AF state, a pocket Fermi surface
arises around (π/2, π/2) and a gap opens in the antinode [near(π, 0)]. As δ increases, the Fermi surface around (π/2, π/2)
extends toward the antinodes along the AF Brillouin zoneedge. Such behavior resembles the pseudogap phenomena,as the SF state treated in this paper does. Thus, if such features are preserved when the AF long-range order is broken
into a short-range order for some reason, as actually observedin cuprates,127) a (disordered) type-(ii) AF state becomes another candidate for a pseudogap state, although the symmetry
breaking is rather diﬀerent from that in the SF state.
Let us discuss the coexistence with d-SC. The same
study38) as discussed above showed that, although type-(ii)AF states do not coexist with d-SC, metallic AF states fort′ /t & −0.05 [called type-(i) AF] coexist with d-SC; thesetype-(i) AF states have pocket Fermi surfaces in the antinodes. This corroborates the fact that the electron scatteringof q = (π, π) that connects two antinodes is crucial for theappearance of d-SC. From this result, we expect that the SF
state is unlikely to coexist with d-SC because gaps open in
the antinodes in the SF state, as shown in Fig. 13(d). As an
exception, coexistence may be possible for δ ∼ δSF , where
the Fermi surfaces extend to the antinodes, as discussed inRef. 38. Thus, the SF order probably competes with the d-SC
order rather than underlies it.128) We need to directly confirm
this by examining a mixed state of the SF and d-SC orders.
Finally, we consider the possibility of the coexistence of AF
and SF orders. Recently, a Hubbard model with an SF phase,namely, H = H SF + HU [see Eqs. (1) and (A·1)], was studied
using a VMC method with a mixed state of SF and AF orders,ΨSF+AF .129) For θ = 0 [Eq. (1) with t′ = 0], the optimized
ΨSF+AF is reduced to ΨAF , which belongs to the type-(i) AF
phase. Namely, the SF order is excluded by the type-(i) AF
order. This is probably because the AF order is energeticallydominant over the SF order, and the loci of Fermi surfaces
compete with each other.

(a)
0.03

t/
-0.1
E
∆

0.02

t

t/
E∆

∆E = ESF - Ed

USF/t
Uc/t

δ=0

(SF)

L

10
12
14

0.01

-0.2

Uc/t (d)
0

0.2

(b)

δ=0

t/
E∆U 0.1

U/t

10

20

L

L=12 δ

10
12
14

30

0
0.02780.05560.0833

0
0

10

U/t

20

30

Fig. 23. (Color online) The two components of the energy diﬀerence between the d-SC and SF states are shown for some doping rates (t′ /t = 0,
L = 12): (a) kinetic energy and (b) interaction energy parts. The symbols arecommon to all panels. For δ = 0, data for L = 10 and 14 are added. The
arrows indicate USF /t in the SF state for δ > 0. The inset in (a) shows thediﬀerence in total energy (∆E = ∆Et + ∆EU ) for four values of δ.

(ET)2 Cu2 (CN)3 with t′ /t ∼ 0.8 shows Fermi-liquid-like behavior above T c . Since the present result shows that the SF
state is not stabilized for the case of t′ /t ∼ 0.8, the normal
state of κ-(ET)2 Cu2 (CN)3 is naturally understood on the basisof ΨN . Although our results are consistent with experiments,
quantitative discussions will be necessary to determine the effective value of U/t as well as t′ /t more accurately for eachcompound.124, 125)
For the organic conductors with finite doping, we find thatΨSF is not stabilized at δ = 0.11 for both U/t = 8 and 12, regardless of the value of t′ /t. Thus, concerning the pseudogap
phenomena found in a doped κ-ET salt,79) we cannot conclude
that the SF state is a candidate for the pseudogap phase. Other
factors may be necessary to understand this pseudogap.

7. Conclusions
In this paper, we studied the stability and other propertiesof the staggered flux (SF) state in the two-dimensional Hubbard model at and near half filling. We carried out systematiccomputations for U/t, t′ /t, and δ, using a variational Monte
Carlo method, which is useful for treating correlated systems.
In the trial SF state, a configuration-dependent phase factor
was introduced, which is vital to treat a current-carrying statein the regime of Mott physics. In this SF state, we found agood possibility of explaining the pseudogap phenomena inhigh-T c cuprates and κ-ET salts. The main results are summarized as follows:(1) The SF state is not stabilized in a weakly correlated
regime (U/t . 5), but becomes considerably stable in astrongly correlated regime [Figs. 5 and 10(a)]. The physicalproperties in the latter regime are consistent with those of thet-J model.32, 43, 44, 60) The transition from ΨN to ΨSF at USF /t
is probably continuous.
(2) At half filling (δ = 0), the SF state becomes Mott insulating for U > Uc ∼ 7t. A metallic SF state is realized for

6.4 Related studies and coexistence with other orders
A decade ago, Yang, Rice, and Zhang introduced a phenomenological Green’s function that can represent various anomalous properties of the pseudogap phase.126) Their
Green’s function contains a self-energy that reproduces thed x2 −y2 -wave RVB state at half filling. For finite doping, theGreen’s function is assumed to have the same self-energy but
without the features of SC. On the other hand, the SF stateused in the present paper is also connected to the d x2 −y2 -wave
RVB state due to the SU(2) symmetry at half filling. For finite
doping, however, the SF state does not show SC. Therefore,
we expect a close relationship between the present SF state
and the phenomenological Green’s function, although the explicit correspondence is not known.
As mentioned in Sect. 1, an AF state was recently studied
by applying a VMC method with a band-renormalization effect to the Hubbard model on the frustrated square lattice.38)
13

J. Phys. Soc. Jpn.

DRAFT

USF < U < Uc , which is gapless in the charge degree of freedom but gapped in the spin sector. This gap behavior of themetallic SF state at δ = 0 continues to the doped cases of
U > USF . However, it is distinct from the case of the noninteracting SF state ΦSF in the sense that spin-charge separation
occurs. In doped cases, ΨSF has a segmentary Fermi surface
near the nodal point (π/2, π/2) but is gapped near the antinodal (π, 0) (Fig. 13). By analyzing the kinetic energies (Et ),we found that Et (Ψd ) < Et (ΨSF ) < Et (ΨN ) for a large U/t,meaning that a kinetic-energy-driven SC takes place even if
we assume that the SF state is realized above T c .
(3) Although ΨSF is unstable toward phase separation fort′ /t ∼ 0 in accordance with the feature in the t-J model,61)
ΨSF restores stability against inhomogeneity for t′ /t ∼ −0.3.
This aspect is similar to that of AF states.35, 38, 113)
(4) For the simple square lattice (t′ /t = 0), the stable SF
area is δ . 0.1. In the anisotropic triangular lattice (|t′ /t| > 0),
this area does not expand. In the frustrated square lattice, however, the t′ term makes this area expand to δ . 0.16 for−0.4 . t′ /t . −0.1 (hole-doped cases) but shrink to a very
close vicinity of half filling for t′ /t > 0 (electron-doped cases)(Fig. 15). This change is mostly caused by the sensitivity ofΨN to t′ , while ΨSF is insensitive to t′ because it is defined
suitably for the square-lattice plaquettes. This result may berelated to the fact that pseudogap behavior is not clearly observed in electron-doped cuprates.130–135)(5) On the basis of this study and another study,38) we argue
that the SF state does not coexist with d-SC as a homogeneousstate and is not an underlying normal state from which d-SC
arises. This is because the SF state has no Fermi surface inthe antinodes necessary for generating d-SC. However, when
the optimized θ becomes small (for δ ∼ δSF ), coexistence ispossible. The coexistence of SF and AF orders also does notoccur for t′ /t & −0.05;129) further study is needed to clarify
the cases of t′ /t . −0.05.
(6) The local circular current in a plaquette, which is anorder parameter of the SF phase, is strongly suppressed in thelarge-U/t region, but it does not vanish even in the insulatingphase. A so-called chiral Mott insulator is realized.
(7) We showed that the spin current state (or spin-nematicstate) is not stabilized for the t-J and Hubbard models.
Because these results are mostly consistent with the behaviors in the pseudogap phase of cuprates, the SF state should bereconsidered as a candidate for the anomalous ‘normal state’competing with d-SC in the underdoped regime. Note that the
AF state is considerably stabilized in a wide region of the
Hubbard model.38) Therefore, possible disordered AF [type(ii)] is another candidate for the pseudogap phase, although
the symmetry breaking is diﬀerent. Besides this claim, there
are relevant subjects left for future studies. (i) What will be
the phase transition between the SF and d-SC states if the
SF state is the state above T c ? (ii) In this study, we treatedthe SF and d-SC states independently. However, it is important to check directly whether the two orders coexist,15, 28, 29)
and how the coexistent state behaves, if it exists.44, 136) (iii) In
this study, we introduced a phase factor for the doublon-holon
processes. It will be worthwhile to search a useful phase factor
that controls isolated (doped) holons for δ > 0. (iv) It will be
intriguing to search for a low-lying circular-current state otherthan the SF state and the state proposed by Varma.90, 121, 122)
We thank Yuta Toga, Ryo Sato, Tsutomu Watanabe, Hiroki

Tsuchiura, and Yukio Tanaka for useful discussions and information. This work was supported in part by Grants-in-Aid
from the Ministry of Education, Culture, Sports, Science andTechnology.
Appendix A:

Details of Bare Staggered Flux State

y

(a)

(b)

B j-x+yB j-2xAj

Spin current

B
Bj

x

B

B j-x-y

A

B
B

B

A

B

B

B

↑- spin↓-spinFig. A·1. (Color online) (a) Coordinates in the extended unit cell used in
HSF . (b) Currents in staggered spin current state for up and down spins discussed in Appendix C. The arrows denote the directions of easy flow.

In this Appendix, we give a definition of the one-body SF
state ΦSF used in this study (see Sect. 2.2) and summarize its
characteristic properties. ΦSF is the ground state of a noninteracting SF Hamiltonian H SF (θ) written as∑[ ()
H SF = −teiθ c†A j,σ cB j,σ + c†A j,σ cB j−2x,σ
j∈A,σ

(
)
]
+e−iθ c†A j,σ cB j−x+y,σ + c†A j,σ cB j−x−y,σ + H.c.

(A·1)

in the sublattice (A,B) representation [see Fig. A·1(a)]. Here,
we abbreviate r j (the position of site j) as j, and x (y) is theunit vector in the x (y) direction. For θ = 0, H SF is reduced to
Ht in Eq. (1). H SF is diagonalized as∑[]
H SF =
E−SF (k) α†kσ αkσ + E+SF (k) β†kσ βkσ ,
(A·2)
k,σ

with the band dispersions given as
E±SF (k) = ±2t Sθ,k ,by applying the unitary transformation
(
)
(
)(
)
1
cAkσ
γθ,k γθ,k
αkσ
= √
,
cBkσ
1
−1
βkσ
2
where
(
)
e−ikx eiθ cos k x + e−iθ cos ky
γθ,k =
Sθ,k

(A·3)

(A·4)

(A·5)

with Sθ,k given in Eq. (6). The lower band dispersion can betransformed to the form of Eq. (7). The one-body SF state forn ≤ 1 is given by filling the lower band as
∏
∏
)
1 (
†
†
ΦSF =
α†kσ |0⟩ =√ γθ,k cAkσ + cBkσ |0⟩,2
k∈kF (θ),σ
k∈kF (θ),σ(A·6)
which leads to Eq. (4) by applying the Fourier transformation,√ ∑2
eik·r j cΛkσ .
(Λ = A, B)(A·7)
cΛ jσ =
Ns k
Because ΦSF is a current-carrying state, ΦSF is essentiallycomplex except when 4θ = 0 and π.
14

J. Phys. Soc. Jpn.

DRAFT

0
-0.2

t/-0.4
~E

δ

0
0.1
0.2
0.3
0.4
0.5

-0.6
-0.8

bard model with U = 0 (Ht ), Ẽ/t is positive because the exact
ground state of Ht is ΦN [ΦSF (θ = 0)]. Ẽ/t monotonically increases as θ increases, as shown in Fig. A·2(b). For θ ∼ 0, ESF
in Eq. (A·11) increases quadratically as()2∑ cos k x − cos ky+ · · · , (A·13)
ESF = EN + θ2 t
k∈kF (θ),σ cos k x + cos ky(a)

U / t= 0

SF modelat least for δ = 0. Hence, the SF state is unlikely to be stabilized even if U/t is added as a perturbation; this feature is inagreement with that for U < USF discussed in Sects. 3 and 4.

1.5 (b) Analytic
δ=0
t/
~E

1

0.5

VMC
δ=0

L
10
20
δ = 0.08
20

0

t/

0.1

0.2

θ/π

0.3

0.4

-1

)

k -2

U / t= 0
Hubbard model0

Stag. flux

E－

FS (

θ/π

0
0.05
0.1
0.15
0.25
0.375
0.5

-3

0.5

-4

Fig. A·2. (Color online) Total energies per site of ΦSF measured from thatof ΦN [Eq. (A·8)] are drawn as functions of θ for (a) the SF Hamiltonian[Eq. (A·1)] and (b) the Hubbard Hamiltonian [Eq. (1)]. In (a), Eqs. (A·10)and (A·12) are used with a common θ in ΦSF and H SF as the ground state.
In (b), Eqs. (A·9) and (A·11) are used with θ being varied in ΦSF ; we plot
VMC data for some cases in addition to the analytic result at half filling. We
confirmed that P has almost no eﬀect for U = 0. t′ is fixed at 0.

(0,0)

k

π π

( /2, /2)

(0,0)

Fig. A·3. (Color online) Energy dispersion of the bare SF state ΦSF forseveral values of θ along (0, 0)-(π, 0)-(π/2, π/2)-(0, 0).

The (lower) band structure of the bare SF state, E−SF (k)[Eq. (A·3)], is shown in Fig. A·3 for several values of θ. In
the ordinary Fermi sea (θ = 0), the band top is degenerate
along the AF Brillouin zone edge (π, 0)-(0, π)-(−π, 0)-(0, −π),namely, the nesting condition is completely satisfied at halffilling. By introducing θ, this degeneracy is lifted and the band
top becomes located at (π/2, π/2) and the three other equivalent points. In particular, for the π-flux state (θ = π/4), theband top forms an isotropic Dirac cone centered at (π/2, π/2).
This cone becomes elongated in the (π, 0)-(0, π) direction as θ
decreases from π/4.
This peculiar band structure brings about anomalous properties in ΦSF . At half filling, the state for θ > 0 is not a conventional metal. Although it is not explicitly shown here, n(k)is a smooth continuous function except for a discontinuity at
k = (π/2, π/2), and N(q) = S (q) becomes a quadratic function of |q| for |q| → 0. In a doped case, a Fermi surface appears
that is made of a cross section of the elongated Dirac conenear k = (π/2, π/2), which is shown in Fig. A·4 for somevalues of θ and δ and is reminiscent of a Fermi arc or hole
pocket observed in cuprates by ARPES and so forth. This is
in contrast to the d-SC state Φd , in which the Fermi surface isa point on the nodal line irrespective of the value of δ. As δ increases or θ decreases, this segmentary Fermi surface of ΦSFbecomes longer, and the gap region shrinks to the vicinity of
the antinodal points. However, the behavior of N(q) [= S (q)]for |q| → 0 in doped ΦSF is basically unchanged from thatat half filling. This gap behavior diﬀers from the case of themetallic SF state ΨSF for U > USF , as discussed in Sects. 3.3
and 4.2.
The local circular current defined by Eq. (18) is calculated

The total energy per site of ΦSF measured from that of thebare Fermi sea ΦN is written asẼ = ESF − EN .

π

( ,0)

(A·8)

Here, EN for ΨN is obtained for the Hubbard model [Eq. (1)with U = 0] through∑1 ⟨ΦN |Ht |ΦN ⟩
1
EN =
=
εk ,
(A·9)
Ns ⟨ΦN |ΦN ⟩Ns k∈k (θ=0),σF

where εk [Eq. (8)] is the bare dispersion of an ordinary Fermisea, and EN is obtained for the SF model [Eq. (A·1)] through∑1 ⟨ΦN |H SF |ΦN ⟩1
EN ==E SF (k). (A·10)
Ns ⟨ΦN |ΦN ⟩Ns k∈k (θ=0),σ −FThe energy of ΦSF for the Hubbard model is given by()2
2t cos θ ∑ cos k x + cos ky⟨ΦSF |Ht |ΦSF ⟩
=−
,
ESF =
Ns ⟨ΦSF |ΦSF ⟩Ns k∈k (θ),σ
Sθ,kF(A·11)and that for the SF Hamiltonian is given by1 ⟨ΦSF |H SF |ΦSF ⟩
1 ∑
ESF =
=
E SF (k). (A·12)
Ns ⟨ΦSF |ΦSF ⟩Ns k∈k (θ),σ −F

In Fig. A·2(a), we show Ẽ/t for H as a function of θ. Because ΦSF is the exact ground state of H SF , Ẽ/t ≤ 0 holds; athalf filling, ESF and EN are identical because the Fermi surfaces of ΦN and ΦSF are identical, but the energy of ΦSF issizably reduced as δ or θ increases. In contrast, for the HubSF

15

J. Phys. Soc. Jpn.
1

DRAFT
1

0.2

1

0.4

(a)

θ / π = 0.0625

0.1

0.8

0.8

0.3

0.2

0.4

π
/

0.6

y

k

π
/

δ=0

0.6

(b)

(c)

θ / π = 0.125

θ / π = 0.25

0.8

0.5

δ=0

π
/

y

k

0.4

0.2

0.2

0.4

0.6

0.8

0.6

0.1

δ=0

y

k

0.4

0.2

0

0.4
0.3
0.2

0.5

0.1

1

0.4

0.2

0

0.2

kx / π

0.4

0.6

0.8

1

0

0.2

0.4

kx / π

0.6

0.8

1

kx / πFig. A·4. (Color online) Fermi surfaces of the bare SF state ΦSF shown for (a) θ = 0.0625, (b) θ = 0.125, and (c) θ = 0.25 (π flux) for some doping rates. kpoints near kF (k ∈ kF ) for L = 1000 are plotted. The (π, π)-folded Brillouin zone edge is shown by a gray dash-dotted line.

1.5

δ

0.0
0.05
0.1
0.2

1

t
/
|

U=0

(L = 200)

ΦSF

U/t=12J|C

0.5

ΨSF (VMC)

0

0.1

0.2

θ/π

0.3

L = 12
δ

0.0
0.02780.05560.0833

0.4

=

J

)
∑(
1
Si · S j − ñi ñ j ,4⟨i, j⟩

H3site=−J ∑ ( †c̃ c̃ j,−σ c̃†j+τ,σ c̃ j+τ′ ,σ
4 j,τ,τ′ ,σ j,−σwhere c̃ jσ†
1 ∑α,β c jα σαβ c jβ with σ being the Pauli matrix for S = 1/2
2
spins. We call Ht- J the t-J model and H3 the three-site model.
In doped cases, the behavior of the Hubbard model should be
more similar to that of H3 . Here, we disregard diagonal hopping terms for simplicity.
To these models, we apply a VMC scheme similar to thatfor the Hubbard model. As a many-body factor, we use only
the complete Gutzwiller projector, P = PG with g = 0,
as in previous studies.32, 43, 44, 60) Thus, Ψ̃SF = PG (0)ΦSF (θ)[Ψ̃N = PG (0)ΦN ] has one [no] variational parameter. Here,we concentrate on the decrease in energy of Ψ̃SF from that of
Ψ̃N ,

0.5

for ΦSF as(A·14)

Ẽ = E SF (θ) − E N ,

F

L
10

-0.02

with
Ht

=

−t

∑ (
)
c̃†iσ c̃ jσ + H.c. ,14-0.095

J
EJ

20

/

26

~

-0.04

-0.1

J/
J
~ E -0.06

In this Appendix, we summarize the stability of the SF statein t-J-type models with calculations of reliable accuracy for acomparison with the Hubbard model treated in the main text.
For this purpose, we include the following three-site (or pairhopping) term H3site , which is the same order as H J [t2 /U(= J/4)] in the strong-coupling expansion:
H3 = Ht- J + H3site ,SF

0

Staggered Flux State for t-J Model

Ht- J = Ht + H J ,(B·5)

where E = E (θ = 0).
N

and is shown in Fig. A·5; the optimal θ is always 0 for
the Hubbard model [Fig. A·2(a)]. For comparison, data forstrongly correlated cases are also plotted. Here, we only pointout two notable features. (i) As δ increases, |JC | decreases forthe noninteracting ΦSF , while |JC | increases for the strongly
correlated ΨSF for the Hubbard model with U/t = 12. (ii) As
the interaction increases, |JC | is markedly reduced.
Appendix B:

(B·3)

)
+c̃†j+τ,−σ c̃ j,−σ c̃†j,σ c̃ j+τ′ ,σ , (B·4)∑ †= c jσ (1 − n j−σ ), ñ j =σ c̃ jσ c̃ jσ , and S j =Fig. A·5. (Color online) Local circular currents of one-body SF state (ΦSF )are shown as functions of Peierls phase for four doping rates. The data arenumerically computed using Eq. (A·14) for L = 200. For comparison, VMC
data of ΨSF for U/t = 12 and L = 12 are shown for four values of δ.

()22t sin θ ∑ cos k x − cos ky
JC = −,
Ns k∈k (θ),σ
Sθ,k

HJ

0.06

-0.08

0.08

0

0.1

0.12

δ = 0.0

Staggered flux

-0.1

(B·1)

θ/π

Square
0.1

0.2

0.3

0.4

0.5

θ/π
Fig. B·1. (Color online) Exchange energy per site Ẽ J /J of Ψ̃SF measured
from that of Ψ̃N at half filling as a function of θ. The graph is symmetric withrespect to θ/π = 1/4 (π-flux). The inset shows the magnification near theminimum (θopt ) of the smaller θ. Data for four values of L are compared.

(B·2)⟨i, j⟩,σ16

J. Phys. Soc. Jpn.

DRAFT

0.04

First, we consider the half-filled case (δ = 0), in which Et(= ⟨Ht ⟩/Ns ) and E3site (= ⟨H3site ⟩/Ns ) vanish; the total energy is given by the exchange term E = E J = ⟨H J ⟩/Ns . In
Fig. B·1, we plot Ẽ J /J (= Ẽ/J) as a function of θ. Ẽ J has
a minimum at θ/π ∼ 0.08. Because Ψ̃SF is equivalent to thed-SC state Ψ̃d = PG Φd owing to the SU(2) symmetry,39, 52)
the minimum energy of Ψ̃SF [e.g., E SF
J /J = −1.1396(4) for
L = 10] coincides with that of Ψ̃d [E dJ /J = −1.1398].32) Thisvalue is very low and broadly comparable to the minimum
energy of the AF state on the same footing, Ψ̃AF = PG ΨAF[E AF /J = −1.1412].53) Thus, the SF state is very stable athalf filling irrespective of the value of J/t.
For doped cases (δ > 0), Et and E3site also make contributions. Figure B·2 shows the θ dependence of the three energy
components of Ψ̃SF measured from those of Ψ̃N . By introducing θ, the exchange energy is lowered (Ẽ J /J < 0), similarly
to the case of half filling, in a wide range of δ [Fig. B·2(a)];|Ẽ J |/t is large, especially near half filling. In contrast, Ẽt /t andẼ3site /J monotonically increase with θ, namely, they destabilize the SF state, and become more marked as δ increases[Figs. B·2(b) and B·2(c)]. For fixed values of δ and J/t, thetotal energy Ẽ/t is the sum of these competing components.
For example, in Fig. B·3, we plot Ẽ/t for the t-J and three-site
models for typical values of δ and J/t of underdoped cuprates.
We find that ΨSF is stable with respect to ΨN in a wide rangeof J/t. Because Ẽ3site /J is disadvantageous to the SF state[Fig. B·2(c)], the decrease in Ẽ3 /t is somewhat smaller thanthat in Ẽt−J /t.

(a)

0.03

∆E

t-J

t/
E0.02
∆
0.01

∆E3 L

20
14
10

t - J

J / t = 0.3

3-site

0
(b)

0.1

t-J
t - J

π
/

t
p

θo0.05

3-site3site L
20
14
10

J / t = 0.3

Staggered flux
0

0.1

0.2

δ
Fig. B·4. (Color online) (a) Doping rate dependence of energy gain of SF
state measured from the energy of ΨN for the t-J (∆E2 ) and the three-site(∆E3 ) models. Data for three values of L are plotted for each model. (b)Optimized values of θ in SF state as functions of doping rate.

compared with Ψ̃N , defined as∆E(SF) = E(N) − E(SF),
0

∆E

t-J

for J/t = 0.3. Note that ∆E has the inverse sign to Ẽ. In
Fig. B·4(b), we show the optimized θ (θopt ). The behavior of
θopt is similar to that of Ẽ, but θopt vanishes abruptly at theboundary δSF owing to finite-size eﬀects. From the systemsize dependence, the range of the SF state seems to expand tosome extent in the thermodynamic limit.

/t

-0.01

t/
-0.02
~E
J
t

J/t
0.2
0.3
0.4
0.5

-0.03
-0.04

(a)

0

t - J

model

L = 2003-site model

(b)
0.1

θ/π

i∈A,σ

(
)
]
+e−iθs(σ) c†Ai,σ cBi−x+y,σ + c†Ai,σ cBi−x−y,σ + H.c. , (C·1)( )0.2 (20)
0.3 (13.3)
0.4 (10)0.5 ( 8)J/t

Staggered Spin Current StateIn this Appendix, we study the staggered spin current (SSC) state ΨSSC = PΦSSC ,40, 63, 64) as illustrated inFig. A·1(b). The one-body state, ΦSSC , is obtained as theground state of the noninteracting SSC model written as∑[()
H SSC = −teiθs(σ) c†Ai,σ cBi,σ + c†Ai,σ cBi−2x,σδ = 0.08

-0.02
-0.03

Appendix C:∆E3 / tt/ -0.01
~ E3

(B·6)

U/t

where s(σ) = 1 or −1 according to whether σ =↑ or ↓. H SSC
is diagonalized in the same way as H SF . The energy dispersion is identical to E±SF (k) [Eq. (A·3)]. Consequently, we have
∏ 1 [
]
†
†
(C·2)
ΦSSC =√ γkσ (θ) cAkσ + cBkσ |0⟩,2
k∈kF ,σ0.2

Fig. B·3. (Color online) Total energies of the SF state (ΨSF ) measured fromthat of ΨN plotted as a function of θ for (a) the t-J model and (b) the threesite model. The doping rate is fixed at 0.08 and the values of J/t are chosen
appropriately for cuprates. The values of U/t in (b) are converted using J =4t2 /U. Arrows in both panels indicate the energy gains ∆E/t by the SF state
for J/t = 0.3.

γkσ (θ) =()
e−ikx eiθs(σ) cos k x + e−iθs(σ) cos ky
Sθ,k

.

(C·3)

ΦSSC has a doubled unit cell but, in contrast to ΨSF , it has nomagnetic flux and preserves the time-reversal symmetry. The
SU(2) symmetry is broken in ΦSSC even at half filling.

Finally, we look at the δ dependence of the stability of Ψ̃SF .
In Fig. B·4(a), we plot the energy gain or diﬀerence of Ψ̃SF as
17

J. Phys. Soc. Jpn.

DRAFT

0

0.15

(a) exchange term

0.15(b) hopping term

Et / t~~

EJ / J

0.05

δ
0.0
0.04
0.08
0.12
0.16
0.20

0.1

~

0.1

-0.05

δ
0.0
0.04
0.08
0.12
0.16
0.20

E3site / J

L =20(c) 3-site term

L = 200.05

L = 20

-0.1

0

0.1

0.2

0.3

0

0.1

0.2

θ/π

0.3

0

0.1

θ/π

0.2

0.3

θ/π

Fig. B·2. (Color online) Energy components of SF state measured from those of ΨN as functions of θ: (a) exchange term ⟨H J ⟩, (b) hopping term ⟨Ht ⟩, and(c) three-site term ⟨H3site ⟩. In each panel, data for several doping rates are plotted. The cusplike foldings appearing for δ > 0 are caused by the discontinuouschange in the occupied k-points in ΦSF as θ varies owing to a finite L.
0.3(a) exchange term
0.6

L=14 δ
0.0

10
14
20

0.041
0.061

(c) 3-site term

hopping termt-J-type modelt - J-type model

δ=0 L

0.020

J/ 0.4
~ EJ

0.2
(b)

L= 14

δ

0.2

t/
~E

0.082
0.102

0.0

0.020

0.020
0.041

0.041

J/
0.1
~ E3

0.061
0.082

t

0.122

L= 14

δ

0.0

0.102
0.122

0.061
0.082
0.102
0.122

0.1
0.2

Stag. spin current
t-J-type model
0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

θ/π

0.3

θ/π

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

θ/π

Fig. C·1. (Color online) Energy components of staggered spin current state measured from those of ΨN = ΨSSC (θ = 0), similarly to in Fig. B·2. In (a), data
for three values of L are compared at half filling to show the small finite-size eﬀect.

For the noninteracting Hubbard model, it is trivial that Ẽ
[= E SSC (θ) − E N ] increases as θ increases because ΨN is theexact ground state and E SSC (θ) is equivalent to E SF (θ) [seeEq. (A·13)]. To consider strongly correlated cases, we applyΨ̃SSC = PG (0)ΦSSC to the t-J-type model ⟨H3 ⟩ in Eq. (B·1) inthe same manner as Ψ̃SF in Appendix B. Figure C·1 shows theθ dependence of the three energy components of H3 . The behavior of Ẽt /t and Ẽ3 /J is similar to that of Ψ̃SF [Figs. B·2(b)and B·2(c)]. In contrast to ΨSF , however, Ẽ J /J also monotonically increases with θ. The cause of this diﬀerence is discussed in Sect. 6.1. Because every component of the energy
increases as θ increases, we conclude that ΨSSC never has alower energy than ΨN for δ ∼ 0 and positive J/t.
Although we have argued that ΨSSC is not stabilized in
the square-lattice t-J and Hubbard models, a recent VMC
study125) showed that an SSC state has a lower energy than theparamagnetic state for the Heisenberg model on anisotropic
triangular lattices with J ′ ∼ J in magnetic fields.118)(2012).
3) M. Ogata and H. Fukuyama, Rep. Prog. Phys. 71, 036501 (2008).
4) P. W. Anderson, P. A. Lee, M. Randeria, T. M. Rice, N. Trivedi, andF. C. Zhang, J. Phys.: Condens. Matter 16, R755 (2004).
5) J. Zaanen, G. A. Sawatzky, and J. W. Allen, Phys. Rev. Lett. 55, 418(1985).
6) M. R. Norman, H. Ding, M. Randeria, J. C. Campuzano, T. Yokoya,
T. Takeuchi, T. Takahashi, T. Mochiku, K. Kadowaki, P. Guptasarma,and D. G. Hinks, Nature 392, 157 (1998).
7) T. Yoshida, X. J. Zhou, T. Sasagawa, W. L. Yang, P. V. Bogdanov,
A. Lanzara, Z. Hussain, T. Mizokawa, A. Fujimori, H. Eisaki, Z.X. Shen, T. Kakeshita, and S. Uchida, Phys. Rev. Lett. 91, 027001(2003); T. Yoshida, X. J. Zhou, K. Tanaka, W. L. Yang, Z. Hussain, Z.X. Shen, A. Fujimori, S. Sahrakorpi, M. Lindroos, R. S. Markiewicz,
A. Bansil, S. Komiya, Y. Ando, H. Eisaki, T. Kakeshita, and S. Uchida,Phys. Rev. B 74, 224510 (2007).
8) 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).
9) N. Baris̆ić, S. Badoux, M. K. Chan, C. Dorow, W. Tabis, B. Vignolle,
G. Yu, J. Béard, X. Zhao, C. Proust, and M. Greven, Nat. Phys. 9, 761(2013).
10) Studies in this line are broadly classified into two groups: (i) SC fluctuation or precursors to SC and (ii) resonating valence bond (RVB) theory. It was found that the SC fluctuation appears below T pair with T c <T pair < T ∗ by measurements of the Nernst eﬀect,11, 12) ARPES,13)
and optical conductivity.14, 15) Therefore, the SC fluctuation will not
be related to the pseudogap below T ∗ . In contrast, the RVB mean-field
theory successfully explains the pseudogap, particularly, in the spin

1) For a recent review, see T. Yoshida, M. Hashimoto, I. M. Vishik, Z.X. Shen, and A. Fujimori, J. Phys. Soc. Jpn. 81, 011006 (2012).
2) K. Fujita, A. R. Schmidt, E.-A. Kim, M. J. Lawler, D. H. Lee,J. C. Davis, H. Eisaki, and S. Uchida, J. Phys. Soc. Jpn. 81, 011005

18

J. Phys. Soc. Jpn.

DRAFT

sector, because the spin degrees of freedom form RVBs, which represent a singlet pairing.3, 4) In this theory, T ∗ is not a phase-transition
temperature but a crossover temperature.
11) L. Li, Y. Wang, S. Komiya, S. Ono, Y. Ando, G. D. Gu, and N. P. Ong,
Phys. Rev. B 81, 054510 (2010).
12) F. Rullier-Albenque, H. Alloul, and G. Rikken, Phys. Rev. B 84,
014522 (2011).
13) T. Kondo, A. D. Palczewski, Y. Hayama, T. Takeuchi, J. S. Wen,
Z. J. Xu, G. Gu, and A. Kaminski, Phys. Rev. Lett. 111, 157003(2013).
14) A. Dubroka, M. Rössle, K. W. Kim, V. K. Malik, D. Munzar,
D. N. Basov, A. A. Schafgans, S. J. Moon, C. T. Lin, D. Haug,
V. Hinkov, B. Keimer, Th. Wolf, J. G. Storey, J. L. Tallon, and C. Bernhard, Phys. Rev. Lett. 106, 047006 (2011).
15) E. Uykur, K. Tanaka, T. Masui, S. Miyasaka, and S. Tajima, Phys. Rev.
Lett. 112, 127003 (2014).
16) A. Kaminski, S. Rosenkranz, H. M. Fretwell, J. C. Campuzano, Z. Li,
H. Raﬀy, W. G. Cullen, H. You, C. G. Olsonk, C. M. Varma, andH. Höchst, Nature 416, 610 (2002); A. Shekhter, B. J. Ramshaw,
R. Liang, W. N. Hardy, D. A. Bonn, F. F. Balakirev, R. D. McDonald, J. B. Betts, S. C. Riggs, and A. Migliori, Nature 498, 75 (2013).
17) B. Fauqué, Y. Sidis, V. Hinkov, S. Pailhès, C. T. Lin, X. Chaud, and
P. Bourges, Phys. Rev. Lett. 96, 197001 (2006).
18) V. Balédent, B. Fauqué, Y. Sidis, N. B. Christensen, S. Pailhès, K. Conder, E. Pomjakushina, J. Mesot, and P. Bourges, Phys. Rev. Lett. 105,
027004 (2010).
19) Y. Li, V. Balédent, N. Baris̆ić, Y. Cho, B. Fauqué, Y. Sidis, G. Yu,
X. Zhao, P. Bourges, and M. Greven, Nature 455, 372 (2008); Y. Li,
V. Balédent, G. Yu, N. Baris̆ić, K. Hradil, R. A. Mole, Y. Sidis, P. Steffens, X. Zhao, P. Bourges, and M. Greven, Nature 468, 283 (2010);
Y. Li, G. Yu, M. K. Chan, V. Balédent, Y. L. N. Baris̆ić, X. Zhao,
K. Hradil, R. A. Mole, Y. Sidis, P. Steﬀens, P. Bourges, and M. Greven,Nat. Phys. 8, 404 (2012).
20) M. Fujita, H. Hiraka, M. Matsuda, M. Matsuura, J. M. Tranquada,
S. Wakimoto, G. Xu, and K. Yamada, J. Phys. Soc. Jpn. 81, 011007(2012).
21) G. Ghiringhelli, M. Le Tacon, M. Minola, S. Blanco-Canosa, C. Mazzoli, N. B. Brookes, G. M. De Luca, A. Frano, D. G. Hawthorn, F. He,
T. Loew, M. M. Sala, D. C. Peets, M. Salluzzo, E. Schierle, R. Sutarto,
G. A. Sawatzky, E. Weschke, B. Keimer, and L. Braicovich, Science
337, 821 (2012).
22) J. Chang, E. Blackburn, A. T. Holmes, N. B. Christensen, J. Larsen,
J. Mesot, R. Liang, D. A. Bonn, W. N. Hardy, A. Watenphul,M. v. Zimmermann, E. M. Forgan, and S. M. Hayden, Nat. Phys. 8,871 (2012).
23) A. J. Achkar, R. Sutarto, X. Mao, F. He, A. Frano, S. BlancoCanosa, M. Le Tacon, G. Ghiringhelli, L. Braicovich, M. Minola,
M. M. Sala, C. Mazzoli, R. Liang, D. A. Bonn, W. N. Hardy,
B. Keimer, G. A. Sawatzky, and D. G. Hawthorn, Phys. Rev. Lett.
109, 167001 (2012).
24) K. Nakayama, T. Sato, T. Dobashi, K. Terashima, S. Souma, H. Matsui, T. Takahashi, J. C. Campuzano, K. Kudo, T. Sasaki, N. Kobayashi,
T. Kondo, T. Takeuchi, K. Kadowaki, M. Kofu, and K. Hirota, Phys.
Rev. B 74, 054505 (2006).
25) S. E. Sebastian, N. Harrison, E. Palm, T. P. Murphy, C. H. Mielke,
R. Liang, D. A. Bonn, W. N. Hardy, and G. G. Lonzarich, Nature 454,
200 (2008).
26) E. Razzoli, Y. Sassa, G. Drachuck, M. Mansson, A. Keren, M. Shay,
M. H. Berntsen, O. Tjernberg, M. Radovic, J. Chang, S. Pailhès,
N. Momono, M. Oda, M. Ido, O. J. Lipscombe, S. M. Hayden,
L. Patthey, J. Mesot, and M. Shi, New J. Phys. 12, 125003 (2010).
27) R.-H. He, X. J. Zhou, M. Hashimoto, T. Yoshida, K. Tanaka, S.K. Mo, T. Sasagawa, N. Mannella, W. Meevasana, H. Yao, M. Fujita, T. Adachi, S. Komiya, S. Uchida, Y. Ando, F. Zhou, Z. X. Zhao,
A. Fujimori, Y. Koike, K. Yamada, Z. Hussain, and Z.-X. Shen, New
J. Phys. 13, 013031 (2011).
28) For instance, Y. H. Liu, Y. Toda, K. Shimatake, N. Momono, M. Oda,
and M. Ido, Phys. Rev. Lett. 101, 137003 (2008).
29) J. L. Tallon, F. Barber, J. G. Strey, and J. W. Loram, Phys. Rev. B 87,
140508 (2013).
30) T. Yoshida, W. Malaeb, S. Ideta, D. H. Lu, R. G. Moor, Z.-X. Shen,
M. Okawa, T. Kiss, K. Ishizaka, S. Shin, S. Komiya, Y. Ando,
H. Eisaki, S. Uchida, and A. Fujimori, Phys. Rev. B 93, 014513(2016).
31) H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 57, 2482 (1988);
C. Gros, Ann. Phys. (New York) 189, 53 (1989).
32) H. Yokoyama and M. Ogata, J. Phys. Soc. Jpn. 65, 3615 (1996).
33) H. Yokoyama, Y. Tanaka, M. Ogata, and H. Tsuchiura, J. Phys. Soc.
Jpn. 73, 1119 (2004).
34) H. Yokoyama, M. Ogata, and Y. Tanaka, J. Phys. Soc. Jpn. 75, 114706(2006).
35) H. Yokoyama, M. Ogata, Y. Tanaka, K. Kobayashi, and H. Tsuchiura,
J. Phys. Soc. Jpn. 82, 014707 (2013).
36) A. Paramekanti, M. Randeria, and N. Trivedi, Phys. Rev. B 70, 054504(2004).
37) D. Tahara and M. Imada, J. Phys. Soc. Jpn. 77, 114701 (2008).
38) R. Sato and H. Yokoyama, J. Phys. Soc. Jpn. 85, 074701 (2016).
39) I. Aﬄeck and J. B. Marston, Phys. Rev. B 37, 3774 (1988).
40) H. J. Schulz, Phys. Rev. B 39, 2940 (1989).
41) P. Lederer, D. Poilblanc, and T. M. Rice, Phys. Rev. Lett. 63, 1519
(1989); D. Poilblanc and Y. Hasegawa, Phys. Rev. B 41, 6989 (1990).
42) F. C. Zhang, Phys. Rev. Lett. 64, 974 (1990).
43) S. Liang and N. Trivedi, Phys. Rev. Lett. 64, 232 (1990).
44) T. K. Lee and L. N. Chang, Phys. Rev. B 42, 8720 (1990).
45) M. Ogata, B. Douçot, and T. M. Rice, Phys. Rev. B 43, 5582 (1991).
46) P. A. Lee and X.-G. Wen, Phys. Rev. B 63, 224517 (2001).
47) J. Kishine, P. A. Lee, and X.-G. Wen, Phys. Rev. Lett. 86, 5365 (2001).
48) Q.-H. Wang, J. H. Han, and D.-H. Lee, Phys. Rev. Lett. 87, 167004(2001).
49) H. Tsuchiura, M. Ogata, Y. Tanaka, and S. Kashiwaya, Phys. Rev. B
68, 012509 (2003).
50) T. Kuribayashi, H. Tsuchiura, Y. Tanaka, J. Inoue, M. Ogata, and
S. Kashiwaya, Physica C 392-396, 419 (2003).
51) I. Aﬄeck, Z. Zou, T. Hsu, and P. W. Anderson, Phys. Rev. B 38, 745(1988).
52) F. C. Zhang, C. Gros, T. M. Rice, and H. Shiba, Supercond. Sci. Technol. 1, 36 (1988).
53) H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 56, 3570 (1987).
54) N. Trivedi and D. Ceperley, Phys. Rev. B 40, 2737 (1989).
55) K. J. Runge, Phys. Rev. B 45, 12292 (1992).
56) M. U. Ubbens and P. A. Lee, Phys. Rev. B 46, 8434 (1992).
57) X.-G. Wen and P. A. Lee, Phys. Rev. Lett. 76, 503 (1996); P. A. Lee,N. Nagaosa, T.-K. Ng, and X.-G. Wen, Phys. Rev. B 57, 6003 (1998).
58) K. Hamada and D. Yoshioka, Phys. Rev. B 67, 184503 (2003).
59) E. Cappelluti and R. Zeyher, Phys. Rev. B 59, 6475 (1999).
60) D. A. Ivanov and P. A. Lee, Phys. Rev. B 68, 132501 (2003).
61) D. A. Ivanov, Phys. Rev. B 70, 104503 (2004).
62) S. Chakravarty, R. B. Laughlin, D. K. Morr, and C. Nayak, Phys. Rev.
B 63, 094503 (2001).
63) A. A. Nersesyan, G. I. Japaridze, and I. G. Kimeridze, J. Phys.: Condens. Matter 3, 3353 (1991).
64) M. Ozaki, Int. J. Quantum Chem. 42, 55 (1992).
65) B. Normand and A. M. Oleś, Phys. Rev. B 70, 134407 (2004).
66) B. Binz, D. Baeriswyl, and B. Douçot, Eur. Phys. J. B 25, 69 (2002).
67) C. Honerkamp, M. Salmhofer, and T. M. Rice, Eur. Phys. J. B 27, 127(2002). The tendency of these RG results seems consistent with thepresent results.
68) T. D. Stanescu and P. Phillips, Phys. Rev. B 64, 220509 (2001).
69) A. Macridin, M. Jarrell, and T. Maier, Phys. Rev. B 70, 113105 (2004).
70) X. Lu, L. Chioncel, and E. Arrigoni, Phys. Rev. B 85, 125117 (2012).
71) J. Otsuki, J. Hafermann, and A. I. Lichtenstein, Phys. Rev. B 90,
235132 (2014).
72) S. Tamura and H. Yokoyama, J. Phys. Soc. Jpn. 84, 064707 (2015).
73) A. Ardavan, S. Brown, S. Kagoshima, K. Kanoda, K. Kuroki, H. Mori,
M. Ogata, S. Uji, and J. Wosnitza, J. Phys. Soc. Jpn. 81, 011004(2012).
74) H. Mayaﬀre, P. Wzietek, C. Lenoir, D. Jérome, and P. Batail, Europhys. Lett. 28, 205 (1994).
75) B. J. Powell and R. H. McKenzie, Rep. Prog. Phys. 74, 056501 (2011).
76) H. Kino and H. Fukuyama, J. Phys. Soc. Jpn. 65, 2158 (1996).
77) Y. Shimizu, M. Maesato, and G. Saito, J. Phys. Soc. Jpn. 80, 074702(2011).
78) Y. Shimizu, H. Kasahara, T. Furuta, K. Miyagawa, K. Kanoda,
M. Maesato, and G. Saito, Phys. Rev. B 81, 224508 (2010).
79) Y. Eto, M. Itaya, and A. Kawamoto, Phys. Rev. B 81, 212503 (2010).
80) H. Yokoyama, S. Tamura, and M. Ogata, JPS Conf. Proc. 3, 012029
19

J. Phys. Soc. Jpn.

DRAFTzone) do not alter as U/t varies; the SF transition tends toward a continuous transition as L increases. In a doped case with finite L, however, the occupied k points suddenly change as θ changes at several
values of θ (= θa , θb , · · · ). Consequently, ∆E becomes a discontinuous function of θ with multiple local minima. Thus, the optimized θswitches, for instance, from θ1 (0 ≤ θ1 < θa ) to θ2 (θa ≤ θ2 < θb ) at
certain values of U/t. The present SF transition for finite L corresponds
to a first-order transition from θ1 (= 0) to θ2 (> 0). However, because
the change in the Fermi surface becomes continuous as L → ∞, the SF
transition is expected to become a continuous transition in this limit.
108) In the t-J model, not only δSF but also δd (the transition point from Ψdto ΨN ) tends to be large as compared with in the Hubbard model.
109) In Fig. 16, we find that the system-size dependence becomes large
when δ approaches δSF . This is because the optimal θ discontinuously
switches from the lower edge of the second continuous segment (θain Ref. 107) to zero (minimum in the first segment) at U = USF . The
value of θa depends on L; θa → 0 as L → ∞.
110) Here, USF /t is estimated without considering the band renormalizationin ΨN for ease of calculation; therefore, the area of the SF state may
shrink to some extent near half filling.
111) X. G. Wen, F. Wilczek, and A. Zee, Phys. Rev. B 39, 11413 (1989).
112) 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).
113) T. Misawa and M. Imada, Phys. Rev. B 90, 115137 (2014).
114) H. Yokoyama, R. Sato, S. Tamura, and M. Ogata, Phys. Proc. 65, 29(2015).
115) T. Watanabe, H. Yokoyama, Y. Tanaka, and J. Inoue, J. Phys. Soc. Jpn.
75, 074707 (2006); Phys. Rev. B 77, 214505 (2008).
116) L. F. Tocchio, H. Feldner, F. Becca, R. Valentı́, and C. Gros, Phys.
Rev. B 87, 035143 (2013).
117) L. F. Tocchio, A. Parola, C. Gros, and F. Becca, Phys. Rev. B 80,
064419 (2009).
118) For instance, S. Fujimoto, Phys. Rev. Lett. 106, 196407 (2011);
A. V. Chubkov and O. A. Starykh, Phys. Rev. Lett. 110, 217210(2013).
119) Generally, the band renormalization also occurs so as to reduce thekinetic energy.38) In Fig. 20, we show the boundary where the eﬀective
band in ΨN is renormalized to t1 = 0 and then the nesting condition,
which is lost for t′ , 0, is restored. Note that this boundary for t′ /t &0.3 behaves similarly to the phase boundary between ΨN and ΨSF ofconcern in the text.
120) J. D. Sau, I. Mandal, S. Tewari, and S. Chakravarty, Phys. Rev. B 87,
224503 (2013).
121) C. M. Varma, Phys. Rev. B 55, 14554 (1997); Phys. Rev. B 73, 155113(2006).
122) C. Weber, T. Giamarchi, and C. Varma, Phys. Rev. Lett. 101, 017001(2014).
123) S. Tamura, Dr. Thesis, Faculty of Science, Tohoku University, Sendai(2016) [in Japanese].
124) In this connection, we mention the renormalization flow in ΨN . As in′ /t at U/t = 8 decreases as L
Fig. 19(a), t1 < t′ holds and the value of tc2
increases for t′ /t . 0.8, indicating that the flow is directed to t′ /t → 0
in this regime. On the other hand, we find that t1 > t′ for t′ /t & 1.3 and′ /t increases with L at U/t = 12, meaning that the flow is directed
tc2
to t′ /t → ∞. Similar analysis has been performed for a paramagnetic
state in the Heisenberg model (U/t → ∞) on an anisotropic triangularlattice, and it was concluded that the boundary is J ′ /J ∼ 0.95.125)125) S. Tamura and H. Yokoyama, Phys. Proc. 58, 10 (2014).
126) K.-Y. Yang, T. M. Rice, and F.-C. Zhang, Phys. Rev. B 73, 174501(2006).
127) 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, P. Dai, S. M. Hayden, A. Hiess,
J. W. Lynn, S.-H. Lee, and F. Doğan, 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).
128) We should rectify the argument in Refs. 80 and 81 that the SF state
may be an underlying normal state from which d-SC arises.
129) Y. Toga and H. Yokoyama, Phys. Proc. 81, 13 (2016), and unpublished.
130) In electron-doped cuprates, pseudogap behavior has been found at theso-called hot spots—the intersections of a quasi-Fermi surface with
the AF-Brillouin-zone boundary—, which are detached from the anti-(2014).
81) H. Yokoyama, S. Tamura, T. Watanabe, K. Kobayashi, and M. Ogata,Phys. Proc. 58, 14 (2014).
82) T. Tohyama and S. Maekawa, Supercond. Sci. Technol. 13, R17(2000).
83) H. Shiba and M. Ogata, J. Phys. Soc. Jpn. 59, 2971 (1990).
84) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
85) T. A. Kaplan, P. Horsch, and P. Fulde, Phys. Rev. Lett. 49, 889 (1982).
86) H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 59, 3669 (1990).
87) A. J. Millis and S. N. Coppersmith, Phys. Rev. B 43, 13770 (1991).
88) Precisely speaking, θ in the wave function slightly deviates from theactually averaged value of the phase variation ⟨θ⟩ in the hopping processes as in Fig. 3. At half filling, ⟨θ⟩ is slightly larger than θ butθ ∼ ⟨θ⟩ holds as a whole. For a suﬃciently large U/t, ϕ = ⟨θ⟩ holds
rather than ϕ = θ.
89) In this process, a phase assignment similar to Eq. (12), such as







) (
)
∏ (


 ∑


j 

Pξ = exp (−1) h j iξ
1 − d j+τ  n j+x + n j−x − n j+y − n j−y 
,






τ
j
can be assumed, where n j = n j↑ + n j↓ and ξ is a variational parameter.
However, we can easily show that Pξ = 1, namely, the phase variation
θ in the hopping of isolated holons is not controlled by this type ofphase factor.
90) C. Weber, A. Läuchli, F. Milla, and T. Giamarchi, Phys. Rev. Lett. 102,
017005 (2009).
91) If we apply Pϕ to a real wave function such as ΨN , ΨAF , or Ψd , thephase factor in Pϕ is optimized as ϕ = 0, namely, Pϕ = 1. Therefore,
expectation values associated with the optimized real state do not alterregardless of whether Pϕ is introduced or not.
92) Y. Toga and H. Yokoyama, Phys. Proc. 65, 29 (2015).
93) S. Tamura and H. Yokoyama, Phys. Proc. 81, 5 (2016).
94) W. L. McMillan, Phys. Rev. 138, A442 (1965).
95) D. Ceperley, G. V. Chester, and D. H. Kalos, Phys. Rev. B 16, 3081
(1977).
96) H. Yokoyama and H. Shiba, J. Phys. Soc. Jpn. 56, 1490 (1987).
97) C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. Rev. Lett. 60,
1719 (1988).
98) Because the band-renormalization eﬀect99) is small in ΨAF as well as
in Ψd even for δ ∼ 0 (See Ref. 38) in the case of t′ = 0, we do not treat
it here. It does not directly aﬀect the present result.
99) A. Himeda and M. Ogata, Phys. Rev. Lett. 85, 4345 (2000).
100) To estimate Uc /t more accurately, the long-range part of the D-H (Jastrow) factor is important.101, 102) Therefore, here we do not pursue thequantitative accuracy of Uc /t.
101) M. Capello, F. Becca, M. Fabrizio, S. Sorella, and E. Tosatti, Phys.
Rev. Lett. 94, 026406 (2005).
102) T. Miyagawa and H. Yokoyama, J. Phys. Soc. Jpn. 80, 084705 (2011).
103) R. P. Feynman, Statistical Mechanics (Benjamin/Cummings, Reading,MA., 1972) Chap. 11; A. Auerbach, Interacting Electrons and Quantum Magnetism (Springer, New York, 1974) Chap. 9.
104) For a procedure of applying SMA to finite-size systems, see, for instance, the treatment in constructing Fig. 5 in Ref. 105, in which a conductive spin-gap state was similarly confirmed in the Bose-Einsteincondensation regime for the attractive Hubbard model.
105) S. Tamura and H. Yokoyama, J. Phys. Soc. Jpn. 81, 064718 (2012).
106) To examine the behavior of S (q) for |q| → 0, we consider the secondorder finite diﬀerence for L calculated using the smallest three |q|points in the x direction:
C(L) =∆2 S (q) S (0, 0) + S (4π/L, 0) − 2S (2π/L, 0)
=
.
∆|q|2
(2π/L)2

If C(L) is positive for L → ∞, S (q) is a quadratic (or higher-order)
function, whereas if |C(∞)| vanishes, S (q) becomes a linear function.
For a typical value of the spin-gap metal phase [U/t = 6.5; brown symbols in the inset in Fig. 8(c)], C(L) becomes 0.154, 0.203, 0.232, and0.257 for L = 10, 12, 14, and 16, respectively. C(L) is an increasing
function of L and likely to remain positive for L → ∞. This tendency
also holds for the Mott insulating regime (U > Uc ). In contrast, for atypical Fermi liquid [U/t = 4; black symbols in the inset in Fig. 8(c)],|C(L)| with C(L) < 0 tends to decrease as L increases, suggesting that
S (q) is linear.
107) At half filling, the occupied k points (within the magnetic Brillouin
20

J. Phys. Soc. Jpn.

DRAFT

nodal points.131) This pseudogap is related to AF correlations132) or
AF orders.133) The d-SC gap also has a maximum at the hot spot.134)
However, recent experiments on high-quality samples showed thatsuch an AF pseudogap is suppressed.135)
131) N. P. Armitage, D. H. Lu, C. Kim, A. Damascelli, K. M. Shen, F. Ronning, D. L. Feng, P. Bogdanov, Z.-X. Shen, Y. Onose, Y. Taguchi,
Y. Tokura, P. K. Mang, N. Kaneko, and M. Greven, Phys. Rev. Lett.
87, 147003 (2001).
132) Y. Onose, Y. Taguchi, K. Ishizaka, and Y. Tokura, Phys. Rev. Lett. 87,
217001 (2001).
133) H. Matsui, K. Terashima, T. Sato, T. Takahashi, S.-C. Wang, H.-

B. Yang, H. Ding, T. Uefuji, and K. Yamada, Phys. Rev. Lett. 94,
047005 (2005).
134) H. Matsui, K. Terashima, T. Sato, T. Takahashi, M. Fujita, and K. Yamada, Phys. Rev. Lett. 95, 017003 (2005).
135) 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. Arita, H. Namatame, M. Taniguchi, D. Ootsuki, K. Sawada,
M. Takahashi, T. Mizokawa, Y. Koike, and A. Fujimori, Nat. Commun. 7, 10567 (2016).
136) R. B. Laughlin, Phys. Rev. B 89, 035134 (2014).

21