Geometry Dependence of the Sign Problem

0.65 MB·English

by  V. I. Iglovikov

Geometry Dependence of the Sign Problem
V. I. Iglovikov,1 E. Khatami,2 and R. T. Scalettar1
1Department of Physics, University of California, Davis, California 95616, USA
2Department of Physics and Astronomy, San Jose State University, San Jose, CA 95192, USA
The sign problem is the fundamental limitation to quantum Monte Carlo simulations of the
statistical mechanics of interacting fermions. Determinant quantum Monte Carlo (DQMC) is one
of the leading methods to study lattice models such as the Hubbard Hamiltonian, which describe
strongly correlated phenomena including magnetism, metal-insulator transitions, and (possibly)
5
exotic superconductivity. Here, we provide a comprehensive dataset on the geometry dependence
1
of the DQMC sign problem for different densities, interaction strengths, temperatures, and spatial
0
lattice sizes. We supplement these data with several observations concerning general trends in
2
the data, including the dependence on spatial volume and how this can be probed by examining

y decoupledclusters,thescalingofthesigninthevicinityofaparticle-holesymmetricpoint,andthe
a correlation between thetotal sign and thesigns for the individualspin species.
M
PACSnumbers: 71.10.Fd,02.70.Uu

7
1
I. INTRODUCTION accessingthelowtemperaturepropertiesofHamiltonians

] like the Fermi-Hubbard model with QMC.
l
e
Monte Carlo simulations of classical systems have Fortunately, there are some situations where the sign
-
r generated very precise information about phases and problem is not manifest. For example, the Boltzmann
t
s phase transitions in statistical mechanics. One dramatic weight often takes the form of the product of two
t. example of the power of the methodology is that of the determinants, one for each of the two electron spin
a
Isingmodel,wherethetransitiontemperatureonacubic species,anditcanhappenthatthesignsoftheindividual
m
lattice is now known to six decimal places1 and critical determinants perfectly match, so that the Boltzmann
- exponentshavebeenevaluatedtofourdecimalplaces.2A weight remains positive. This occurs in the complete
d
n roughlysimilarsituationholdsforunfrustratedquantum parameter range of the Hubbard model with attractive
o spin and boson systems. For example, in the Bose- interactions, enabling a study of superconductivity and
c HubbardHamiltonian,3 the criticalinteractionstrengths charge density wave physics. It also occurs in the
[ at a density of one boson per site for the superfluid Hubbard model with repulsive interactions on bipartite
3 to Mott insulator transition in the ground state, are lattices in the limit of average one particle per lattice
v available4 to an accuracy of better than one part in site (half-filling), so that the Mott transition and long-
2 103. Spatial lattice sizes are somewhat smaller than for range antiferromagnetism can be explored. Other
3 classicalproblems since writing the partitionfunction as instances of situations where the sign problem is absent
8 apathintegralintroducesanadditional‘imaginarytime’ are mentioned in Sec. IIC. However, many of the
2
dimension, but are nevertheless quite large, e.g., up to most interesting questions concerning strong correlation
0
. 104 sites for the Bose-Hubbard example cited above. physics remain inaccessible, most notably the question
1 Fermions (in more than one dimension) are more of whether the two-dimensional (2D) square lattice
0
challenging for two reasons. First, the fermionic action repulsiveHubbardHamiltonianhasalowtemperatured-
5
is non-local: the Boltzmann weight typically takes the wavesuperconductingtransition,sothatitwouldprovide
1
: form of a determinant. Thus, updating all the degrees a good description of cuprate superconductivity.9
v
of freedom has a computation time which scales as We have two main goals in this manuscript. The
i
X a nonlinear power of the system size N. The cost first is to present a set of data for the sign problem in
r of a method like determinant quantum Monte Carlo DQMC for different geometries. These include bipartite
a (DQMC)5 scales as N3, as opposed to the naive linear [one-dimensional chain, ladder, 2D square, and three-
in N scaling (ignoring such complications as critical dimensional (3D) cubic, Lieb, honeycomb and the 1/5-
slowingdown)inmanyclassicalandquantumspin/boson depleted square] lattices, and non-bipartite Kagome
applications. Second,andfarworse,thereisnoguarantee and triangular lattices. We consider DQMC because
the sign of the determinant, which is used as the it is a powerful and widely utilized approach to the
probability, is positive. Although one can formally use correlated electron problem. Our second goal is to
the absolute value as a weight, and include the sign in discern trends in the DQMC sign problem. We will
the measurement, in practice one ends up evaluating consider, for example, several new issues: how the sign
the ratio of two numbers, which becomes dominated by depends on the spatial lattice size, the scaling of the
statisticalerroratlowtemperaturesastheybothbecome sign in the vicinity of particle-hole symmetric (PHS)
very small. This situation is known as the “fermion sign points, and the ‘entanglement’ of the sign as probed
problem”,6,7 whose solution is conjectured to be NP- by the consideration of coupled cluster geometries.
hard.8Atpresent,therefore,thereisnoknownmethodof Althoughstatementsconcerningthefirstofthesepoints,
2
the spatial size dependence, have been made in the II. GENERAL CONSIDERATIONS
literature, numerical data are rather scanty, owing to CONCERNING THE SIGN PROBLEM
the computational limitations existing in the initial
investigations. Specifically, the data in Refs. 7 and 10 A. The Hubbard Hamiltonian
were restricted to 4x4, 6x6, and 8x8 lattices. As a
consequence,thescalingregimewasnotreachedformany
Ourfocusisonthe single-bandHubbardHamiltonian,
of the parameter sets. For example, the average sign
sometimes increases as the system size grows, rather
than decreasing. This situation is rectified here. The Hˆ =− tij c†iσcjσ +c†jσciσ −µ niσ
lattices we study also have various unique features, such hXijiσ (cid:0) (cid:1) Xiσ
as nesting of the Fermi surface, van-Hove singularities 1 1
ibnanthdes,dwehnossiteypoofsssitbalteeesffNec(tωs)on=thPeksiδg(nωp−roǫbkle),mawndewflailtl +UXi (cid:18)ni↑− 2(cid:19)(cid:18)ni↓− 2(cid:19) (1)
examine.
Here c† (c ) is the creation (destruction) operator for
iσ iσ
a fermion with spin σ on site i = 1,2, N and
n = c† c is the number operator. t is th·e··hopping
iσ iσ iσ ij
amplitude between nearest-neighbor sites i and j, U is
the interactionstrength,andµ is the chemicalpotential.
For geometries where there is only one type of hopping,
It is worth noting what we will not cover: There
we set t = t = 1 as the unit of energy. We will
are a number of methods which are closely related to ij
denote the first line of the Hamiltonian Hˆ in Eq. (1)
DQMC in that they involve a Hubbard-Stratonovich
by Kˆ, and the second line by Vˆ. This latter term is
decoupling of the interaction, and a Boltzmann weight
writteninPHSform. (SeeSec.IIC.)Thedifferentmodels
builtfromfermiondeterminants. Theseincludeimpurity
algorithms,11 as well as the dynamical mean field considered in this paper are distinguished solely by the
theory (DMFT)12–14 and its cluster extensions, the geometry encoded in the near-neighbor designation ij
h i
dynamical cluster approximation (DCA),15 and the in the kinetic energy term. Even with a common choice
cellular DMFT.16 A strength of some of these methods t = 1, different geometries have distinct bandwidths W,
the spread of eigenvalues of the U = 0 (single particle)
is that the sign problem is greatly mitigated relative
Hamiltonian. Although it is sometimes the case23 that
to DQMC, at least if the cluster size is not too large.
using W as the scale of kinetic energy, rather than t,
This is true even if the bath degrees of freedom are
discretized.17 Also closely related to DQMC are zero- produces better comparisons across different models, we
temperature algorithms which use e−βHˆ as a projection did not find that to be useful here. We retain the
standard convention of normalizing to t.
operatoronatrialwavefunction.10,18Constraintscanbe
There is, of course, much interest in generalizations
introduced in these ground state methods to eliminate
of the Hubbard Hamiltonian, e.g., to multiple bands,
the sign problem, at the expense of systematic errors
in the solution.19 Despite their relations to DQMC, we longer-range density-density interactions, and Hund’s
rule type interactions. Multiple bands, can in fact be
will not discuss these approaches here. Similarly, within
written in the form of Eq. (1), with the understanding
DQMC itself there are different choices of the manner in
which the Hubbard-Stratonovich field is decoupled.20–22 that the label i incorporates both spatial and band
indices. Thus, from the viewpoint of a DQMC
Here, we base our calculations on only the “density
simulation, setting up the 2D ‘periodic Anderson model’
decoupling”, described in Sec. IIB. Though we do
(PAM) which has a square lattice of spatial sites
not explicitly consider the above related approaches, we
and two orbitals per site, is formally identical to a
expectthatsomeofourresultsandgeneralanalysismay
two layer geometry, in which there is a single orbital
have applications to them as well.
on each site. Thus, with the freedom to choose
U (t ) to be site/orbital (bond) dependent, Eq. (1)
i ij
incorporates Hamiltonians like the PAM. Concerning
intersite (interorbital) and Hund’s rule interactions, the
sign problem is typically much worse than for an on-
site U between fermions of different spin species only.
The remainder of this paper is organized as follows: For example, for a model of the CuO planes of
2
In Sec. II we review the Hubbard Hamiltonian and the cuprate superconductors,24 it was found that the sign
basic formulation of DQMC, followed by some general, problemrestrictedsimulationstointeractionsU .1at
pd
well-known properties of the resulting sign problem. In temperatures where local spin correlations were seen to
Sec.IIIwerecordvaluesfortheaveragesignfordifferent begin to develop. Like DQMC, Hund’s rule interactions
lattice geometries. SectionIV examines generalpatterns also present grave sign problem difficulties in DMFT.25
inthisdata. Finally,Sec.Vcontainsconcludingremarks. We do not explicitly consider them here.
3
B. Determinant Quantum Monte Carlo to perform the trace over the Hilbert space analytically.
The sum over the HS configurations (i,l) is performed
X
The fundamental idea of DQMC26 is to take stochastically using Monte Carlo techniques. The
advantage of the fact that it is possible to compute corresponding Boltzmann weight takes the form of the
analytically the trace of a product of the exponentials productoftwodeterminants(oneforeachspinspecie)of
of quadratic forms of fermion creation and destruction matrices Mσ( ) of dimension N. As this determinant
X
operators. If we denote the vector of creation operators product may be negative for some HS configurations,
(c†,c†,c† c† )by~c†,andA are(symmetric)N N the sampling is done using the absolute values of the
1 2 3 ··· N j × determinant product, and measured expectation values
matrices of real numbers, then
are adjusted accordingly.
Tr(e~c†A1~c e~c†A2~c e~c†A3~c e~c†AL~c )
The average sign S is then defined to be the ratio
···
h i
=det(I +B B B B ) (2) of the integral of the product of up and down spin
1 2 3 L
···
determinants,to the integralofthe absolute value of the
Here B = eAi. It is important to emphasize that product. An analogous definition holds for the average
i
the trace in the left hand side of Eq. (2) is over a sign S of the individual determinants:
σ
2N dimensional Hilbert space of the fermionic operators h i
while the determinant on the right hand side is taken
over a real matrix of dimension N. detM ( )detM ( )
The interaction term in the Hubbard Hamiltonian S = X ↑ X ↓ X
h i P detM ( )detM ( )
Eq. (1) is not quadratic in the fermionic operators, X | ↑ X ↓ X |
P detM ( )
but can be made so by first discretizing the inverse S = X σ X . (6)
temperature β = Lτ and then employing the Trotter h σi PX |detMσ(X)|
approximation27–29 P
Tre−βHˆ =Tr(e−τHˆ e−τHˆ e−τHˆ ) (3) In the case we consider here, with no external magnetic
···
field, by symmetry S = S .
Tr(e−τKˆ e−τVˆ e−τKˆ e−τVˆ e−τKˆ e−τVˆ ) h ↑i h ↓i
≈ ···
As a practical matter, these quantities are obtained
Herethe exponentialofthe fullHamiltonianHˆ =Kˆ +Vˆ by generating configurations with the (non-negative)
is approximatedby the productofthe exponentialsofKˆ weight detM ( )detM ( ) and measuring the ratios
↑ ↓
and Vˆ, a well-controlled procedure which can be made detM (| )detMX( )/ deXtM| ( )detM ( ) and
↑ ↓ ↑ ↓
X X | X X |
arbitrarily accurate by taking τ 0. detM ( )/ detM ( ) for each configuration . In
σ σ
→ X | X | X
The purpose of this procedure is the isolation of Eqs. (4) and (5), we have coupled the HS variable to
the exponential of the interaction term Vˆ, which can the z component of fermionic spin, n n . It is
i↑ i↓
−
then be rewritten using a Hubbard-Stratonovich(HS) also possible to write transformations which involve the
transformation: xy components of a spin, c† c , or even local pairing
i↑ i↓
e−τU(ni↑−12)(ni↓−21) = 12e−Uτ/4 eλX(ni↑−ni↓), (4) opproebraletmor.s2,1,2c2†i↑cIn†i↓.theThatetsreacitnivegeHneurbablawrdorsHeanmtilhteonsiaignn,
XX=±1 the HS variable couples to the charge n + n on
i↑ i↓
site i. This makes the matrices A identical for up
where coshλ = eUτ/2. Because one needs to transform i
and down fermions, so that their determinants are also
the interaction term on every spatial site i = 1,2, N
and also for each of the l =1,2, L exponentials o·f·τ·Vˆ identical, and thus, there is no sign problem. If a charge
··· decoupling is used for the repulsive model, the HS
in Eq. (3), there are a total of NL HS variables (i,l).
X transformationwould involve complex numbers, and the
In Eq. (4), we have employed the discrete HS
determinants would be complex as well, leading to an
transformation introduced by Hirsch,37 but one could
even more challenging ‘phase problem’.
alsouseacontinuousvariable andaGaussianintegral,
X
We note that there are many details omitted in this
e−τU(ni↑−12)(ni↓−21) = e−Uτ/4 d e−X2+2γX(ni↑−ni↓) brief description, including methods to stabilize the
√π Z X product of the B matrices so that round-off errors
i
(5) do not accumulate, the precise Monte Carlo update
procedure(howmany variablesarealteredin eachstep),
with γ = Uτ/2. There are some differences in the more rapid procedures for obtaining the ratio of new
efficiency opf the exploration of phase space between the to old determinants after a HS variable is updated,
discrete and continuous cases.20 how to evaluate non-equal time observables, analytic
Once the HS transformation is introduced, all the continuation to obtain dynamic behavior, and so forth.
exponentials in the trace of Eq. (3) are quadratic in the The reader is referred to Refs. [10,18,30–36] for more
fermionoperators,so,theidentityinEq.(2)canbeused complete discussions.
4
C. Particle-Hole Symmetry and down determinants tend to be rather uncorrelated,
so that the signs of the individual components lend no
On a bipartite lattice, and at µ = 0, the Hamiltonian further information, and can be qualitatively inferred
Eq. (1) is PHS. That is, the Hamiltonian is invariant from the square root of the total sign.
under the transformation c† ( 1)ic . Here ( 1)i = Some ofthese geometrieshaveunique featuresintheir
iσ → − iσ − non-interacting densities of states (DOS). The square
+1( 1)ontheA(B)sublattice. Thissymmetryispresent
− lattice possesses a van-Hove DOS singularity at ρ = 1.
even if t and U vary spatially across the lattice. As one
In contrast, the honeycomb lattice has a DOS which
physical consequence, the density ρ = n +n = 1
i i↑ i↓
h i vanisheslinearlythere. TheLieblatticehasaflatenergy
(half-filling) for all values of t,U and T. Correlation
bandbetweentwodispersingones,while the flatbandin
functions at density ρ and 2 ρ have the same values,
− the Kagomelattice canbe chosento be either the lowest
or are trivially related. Particle-hole transformations
orhighestsetofenergylevels,dependingonthesignoft.
involving only one spin specie can also be used to relate
One of our goals is to examine how such features might
the attractive and repulsive Hubbard Hamiltonians in
affect the fermion sign, an issue to which we will return
this limit.
in the conclusion.
PHS has profound implications for DQMC. When
it is present, the determinants of the up and down
matrices M have the same sign. Thus, although they
σ
individually can go negative, their product is always A. Hypercubic Lattices
positive. As a consequence, low temperature physics
can be accessed at half-filling. This fact enabled DQMC In this subsection we present data for hypercubic
to establish rigorously37,38 that the single band, square lattices; linear chains, ladders, the square lattice, and
lattice Hubbard model has long-range Ne´el order at the cubic lattice. In all cases, we use periodic boundary
T = 0, as opposed to a disordered (resonating valence conditions except for the rungs of the ladder geometry.
bond) ground state. Although the sign does not cause a problem for
Thesignproblemcanbeabsentinsomeothertypesof Hubbard world line methods48,49 in one dimension,
Hamiltonians with a similar symmetry requirement, for DQMCdoeshaveasignprobleminthiscase50. InFig.1,
example in a model with an interaction which takes the the average fermion sign, S (also denoted sign ), for
form of the square of near-neighborhopping,39 for a low the chain geometry is showhnifor fixed U = 4hand βi = 8.
energy theory of the onset of antiferromagnetism,40 and Figure 1(a) shows the doping dependence of the sign,
in spin polarized Fermi systems.41 Indeed, the number whichhasanon-trivialstructure. Mostnotably,itshows
of special situations where the sign problem is absent a peak around ρ = 0.5. One may wonder if such a local
is rapidly growing, including, for example, in quasi-1D maximum arises due to remnants of the ‘shell’ effect.
condensed matter models of ferromagnetism,42 and via That is, at U = 0 the k space grid is sufficiently coarse
the ‘fermion bag’ approach, in lattice gauge theory.43,44 such that ρ(µ) shows distinct plateaus where S tends
h i
The sign problem is also absent in Hubbard models to be closer to one. This phenomenon is well known,
with a larger number of spin components,45 and, very for example, on square lattices that are not too large
interestingly, in a class of spinless fermion models46, a (e.g., 4 4). However, as will be seen in Fig. 3, this
×
unique situation where positivity is not dependent on occurs only on small lattices, and is unlikely the origin
having an even number of fermionic species. ofthe maximumatquarterfillingherewherethek space
grid is much more refined. Another interesting feature
is that S remains small at low density. This is, again,
III. SIGN PROBLEM DATASETS rather dhiffierent from what happens on a square lattice
where S 1asρ 0(Fig.3)orevenladders(Fig.2).
h i→ →
Ourgoalinthissectionistopresentaunifiedandeasily The low value of S as ρ 0 does not appear to be
h i →
comparable collection of data for the sign problem for connectedtothedivergenceofthedensityofstatesatthe
differentlattice geometries,including hypercubic lattices bottom of the band since the same divergence occurs in
in dimensions d = 1,2,3; other bipartite structures like the ladder geometrywhere S recoversto one as ρ 0.
h i →
the honeycomb, Lieb and 1/5-depleted square lattices; Figures 1(b) and 1(c) show the scaling of the average
and finally two non-bipartite lattices: triangular and sign with spatial lattice size N and inverse temperature
Kagome. For each case we will exhibit the average sign β respectively. After plateaux at small N and β, where
for a range of temperatures T, interaction strengths U, S =1, the averagesigndecreases in a manner which is
h i
lattice sizes, and densities. largelyconsistentwithexponential. ln S isnotperfectly
h i
We focus on the product of the signs of the up and linearinβ orN,but exhibitssomedownwardcurvature,
down determinants, since that is what is relevant for whichwebelieve indicatesthe scalingregimehasnotyet
extractingphysicsfromtheDQMCsimulation. However, been fully attained. We will remark on this more fully
in Sec. IVB we will present a brief analysis of the later in this section. As we shall see, this exponential
individual spin components. Among other things, we decreaseisalsothecaseinothergeometries,althoughthe
observe that, except at PHS points, the signs of the up decrease with N once one exits the plateau is in general
5
1.0 1.0
(a) (a)
0.8 0.8
ladder
U=6 β=4
(cid:1)0.6 (cid:2)0.6
n n t /t=1
g g ⟂
si chain si
(cid:0)0.4 U=4 β=8 (cid:0)0.4 2 20
×
2 32
N=48 N=128 2×40
0.2 N=64 N=160 0.2 2×60
N=80 N=200 2×80
N=100 N=300 2×100
0.0 0.0 ×
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 ρ 0.6 0.8 1.0
ρ
0 0
0 0 (b) ρ=0.625 (c)
(b) −1(c) ρ=0.875 −1
lnsign(cid:0)(cid:1)−−−312 Uβ=ρρ===8400..682755 lnsign(cid:0)(cid:1)−−−−4352 UNρρ====41000..6802755 lnsign(cid:0)(cid:1)−−12 Uβ==46 lnsign(cid:0)(cid:1)−−−432 UN==ρρ==61000..0682755
ρ=0.2 ρ=0.2 −3 −5
−4 −6 0 60 120 180 0 2 4
0 75 150 225 0 4 8 12 N β
N β
FIG. 2: (Color online) Same as Fig. 1 but for the ladder
FIG. 1: (Color online) Sign problem for chains of different
geometry with the rung hopping equal to the hopping along
lengths at U = 4 (which equals the bandwidth in this
thechainsanddifferentvaluesoftheinteractionstrengthand
geometry). (a) The dependence of hSi on the density. (b)
β, which are noted in thepanels.
lnhSi vs the system size, N, and (c) vs β at several fixed
densities ρ. The densities used in panels (b) and (c) are
indicatedbyverticaldashedlinesinpanel(a). Theerrorbars
in this figure and in the remaining figures in this paper can system sizes at fixed U,β while the bottom two panels
be inferred from the scatter in the data, and hence, may not show that the decay is consistent with exponential in N
have been shown. Axis label hsigni is referred to as hSi in andβ. (See remarksto follow onchallengesin capturing
text. fully linear behavior of ln S (β).) S has a minimum
h i h i
at filling ρ 0.8, similar to what is known to occur also
≈
for a square lattice. (See Fig. 3.) For t /t > 2 the
⊥
lessabruptthanwithβ. InFig.1(b) atρ=0.625,0.875, noninteracting system is a band insulator (BI). Because
the average sign, despite its exponential decay, remains the BI occursat the particle-hole symmetric density ρ=
quite manageable out to N & 200. Even at ρ = 0.2 the 1, S ispinnedatunity. Therefore,forthis case,wewill
spatial size must be tripled from N 70 before ln S exahmiine the signs of the determinants of the individual
∼ h i∼
4. InFig.1(c),ontheotherhand,thedecaytoln S spin matrices S = S in Sec. IVB.
−4takesplaceafteronlya50%increaseinβ (fromhβ =i∼8 h ↑i h ↓i
The square lattice is the most well-studied Hubbard
−
to β =12).
model geometry, owing to its possible relevance as a
Next, we turn to ladder geometries, which are natural simple model of cuprate superconductivity and d-wave
extensions of chains, before studying the square lattice. pairing driven by antiferromagnetic fluctuations.9 The
Ladders are of interest for several reasons. First, they total sign for the square lattice at U = 6,β = 4 is
have been extensively studied by DMRG51 as a stepping shown as a function of filling ρ for different lattice sizes
stone to 2D. Second, by changing the ratio t⊥/t of the in Fig. 3(a). As mentioned earlier, the peak in the
rung hopping to the hopping along the chains, one can 4 4 lattice occurs as five of the 16 allowed k points
access U = 0 states that are metallic or band-insulating (c×orresponding to a density ρ = 10/16 = 0.625) fill up
at half-filling. The effect of these phase changes on the prior to half-filling. This peak is even more evident10
sign problem for U =0 is one goalof the data presented at U = 4 providing further evidence that, in this case,
6
here and in Sec. IVB. S is connected to the shell structure at U =0, though
h i
Theresultsfortheaveragesignondifferentladdersare the connection appears to diminish with larger U. It is
plotted in Fig. 2. The rung hopping is set to t = 1, so possible that the shell structure would appear on larger
⊥
thattheU =0bandstructureismetallic. AswithFig.1, lattices if lower temperatures were accessible. A rather
the top panel shows the density dependence for different universal feature is the minimum in S at ρ 0.85,
h i ≈
6
1.0 0.0
(a) ρ=0.625
−0.5 ρ=0.875
0.8
square −1.0
square
U=6 β=4 −1.5
0.6 (cid:1)
(cid:1) n
n g N=100
g i−2.0
si(cid:0)0.4 N=16 ns(cid:0)−2.5 β=4
l
N=36
N=64 −3.0
0.2 N=100
N=144 −3.5
N=196
0.0 −4.0
0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 7
ρ U
0
0
(c) FIG.4: (Coloronline)TheU dependenceofthesignproblem
−1 for the squarelattice.
n(cid:1) U=6 n(cid:1)−2 U=6
sig(cid:0)−2 (b) β=4 sig(cid:0)−4 N=100
n n
l−3 ρ=0.875 l ρ=0.625 Here, we briefly generalize our ladder results by
ρ=0.625 −6 ρ=0.875 providing a more complete description of the behavior
−4 of S as the aspect ratio of the lattice goes from one to
0 60 120 180 1 2 3 4 5 6 h i
N β twodimensions. Fig.5showsdataatU =6andβ =6for
lattices with a fixed N = 96 but different aspect ratios.
As mentioned earlier, the chain geometry seems to be
FIG. 3: (Color online) Same as Fig. 1 but for the square
rather unique. For all other cases, the average sign is
geometry and with different parameters, which are noted in
closeto unityforarangeoflowdensities,andthensteps
thepanels.
downwardstosmallvaluesinaregioncenteredatρ 0.8,
∼
recovering only at the PHS point ρ=1.
The precipitous nature of the decrease in S near
which is shared by the ladder geometry. h i
ρ = 1, evident in Fig. 5, is quantified in Fig. 6. Panel
The very nature of the sign problem makes it 6(a) shows the logarithm of the sign versus the doping
challenging to provide a completely compelling linear awayfromhalf-filling. Thelinearbehaviorindicatesthat
plot of ln S as a function of β. For the data in S ea|ρ−1|. Thedecayconstantaislargeandnegative.
Fig. 3c, forhexiample, at β = 5.33 a run with 105 sweeps Ihtsiβ∼and U dependences are given in Fig. 6(b). To our
through the space-time lattice of Hubbard-Stratonovich knowledge, this behavior of S had not been studied
h i
variables (100 spatial sites and 64 time slices) takes before. However, scaling forms for physical observables
severalhoursonawork-station,andgives S =0.0095 like the compressibility κ = ∂ρ/∂µ as one exits the
h i ±
0.0012 (corresponding to ln S = 4.66. The slope Mott phase at ρ = 1 have been suggested.60,61 In these
h i −
dln S /dβ 4, so we can roughly estimate ln S theories, κ follows a power law κ (1 ρ)−η, so that it
h i ∼ − h i ∼ ∝ −
8.9 at β = 6.33, and hence S 0.00013. A diverges just before it vanishes. This also occurs in the
− h i ∼
measurementofthisvalueto10%accuracywouldrequire boson Hubbard model.3,62
anerrorbarofabout0.00001,afactorof 102 lessthan We conclude this section on hypercubic geometries by
∼
the error obtained at β =5.33. Since error bars only go showing the behavior of S for a cubic lattice in Fig. 7.
downasthesquarerootofthe numberofmeasurements, Becausethe numberoflahttiicesitesgrowssorapidlywith
sucharunwouldentail 104timesasmanysweeps,and linearsize,weconsidercaseswhereL =L =L (while
a cpu time of severalmo∼nths. keeping all linear lengths even to axvo6 id fyru6 straztion of
The U dependence of the average sign at fixed β = 4 antiferromagneticcorrelations.) Thequalitativebehavior
and N = 100 for a square lattice is shown in Fig. 4 at is almost identical to that of rectangular lattices, with
ρ = 0.625 and 0.875. The evolution of S with U is a deep minimum in S upon doping from half filling,
h i h i
rathersimilartothatwithβ: Aplateauatweakcoupling followed by a recovery at ρ . 0.6. Curves for different
where S 1 is followedby anabruptdownturn. Thus, sizesalmostcoincidefortheselargedopings,indicatinga
h i≈
in practice, once the sign begins deviating from unity very slow decay with N, as seen in Fig. 7(b). Indeed,
thereisonlyanarrowwindowofstrongercouplingswhere for some densities S even increases as N increases.
h i
data can be acquired. This is demonstrated in Fig. 4 as Presumably, this is a transient phenomenon associated
S asafunctionofU isshowntodecreaseexponentially with the rather small linear lattice lengths which are
h i
once U &4. accessiblein3D.Thedecaywithβ,Fig.7(c),is,asusual,
7
1.0 1.0
(a)
0.8 0.8
cubic
rectangular U=5 β=4.5
(cid:2)0.6 (cid:2)0.6
n U=6 n
g g
si β=6 si
(cid:0)0.4 (cid:0)0.4
1 96 4 4 4
2×48 4×4×6
0.2 4×24 0.2 4×6×6
6×16 6×6×6
8×12 6×6×8
0.0 × 0.0 × ×
0.0 0.2 0.4 ρ 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
ρ
FIG.5: (Coloronline)Theaveragesignasafunctionofρfor 0 0 (c)
therectangularlatticeswiththesamelatticesizeN =96but
different lengths of thesides. n(cid:1)−1(b) n(cid:2)−1 U=5
nsig(cid:0)−2 ρρ==00..682755 nsig(cid:0)−2 N=4×4×6
rapid. l−3 l−3 ρ=0.625
U=5 β=4.5 ρ=0.875
−40 100 200 300 −40 2 4
B. Other Bipartite Lattices N β
FIG. 7: (Color online) Same as Fig. 1 but for the cubic
Hypercubic lattices are just one instance of bipartite
geometry and with different parameters, which are noted in
geometries which are free of the sign problem at half-
thepanels.
filling owing to the PHS. Here we present data for
three additional bipartite geometries, all of which are of
interest because of their materials applications.
Weconsiderfirstthe“Lieblattice”whichconsistsofan
underlying square array of sites with additional two-
A
fold coordinated sites on each bond (see the left panel
B
of Fig. 8). This structure is a more chemically realistic
a
multiband representation of the CuO planes of the t 2
2
a
cupratesuperconductors,withtheCuatomsformingthe t’ 1
square array and bridging O atoms. The relevant filling
of such a three-band model for the cuprates consists of
oneholeperCuO2 unitcell,thatis,wellawayfromhalf- FIG. 8: (Color online) The Lieb lattice (left) and the 1/5-
filling. Indeed, a realistic model of the cuprates would depleted square lattice (right). The latter is taken from
Ref.23. Arrows(dashedsquares)showtheunitvectors(cells)
for each geometry.
0
0 (a)
(b)
−1
−40
n(cid:3)−2 incorporate a substantial energy ǫpd which represents
lnsig(cid:0)−−−435 UUU===754,,,βββ===466 sNqu=a1r0e0 a−−18200 ββββ====4568 cthoemNepavaderdreitdthioetloneasals,Cctouhsetd-hooafrlbfhi-tfioallell.esdtcoasoec,cuwpityhatnhrOee pfe-romrbioitnasl
−6 U=3,β=10 β=10 N=100 per unit cell and ǫ =0, has considerableinterest: Lieb
pd
−70.01 0.04 0.08 0 2 4 6 8 10 has given63 a rigorous demonstration of a ferrimagnetic
ρ 1 U ground state in this situation. The crucial observation
| − |
is that the numbers of sites on the and sublattices
A B
FIG. 6: (Color online) (a) By doping the system away from (NA and NB, respectively) are unequal. Lieb showed
half filling on the square lattice, the average sign decreases that for any bipartite lattice with NB > NA there is a
exponentially, i.e., hSi ∼ ea|ρ−1|. In (b) the decay constant, ‘flat band’ with N N zero energy eigenstates64. A
B A
−
a, is plotted as a function of U for different values of β. The recentDQMCstudyhasexploredtheattractiveHubbard
lattice size is N =10×10. model in this geometry65. The presence of a flat
8
band was found to have important effects on physical 1.0
propertieslikethelocalmomentandpairingcorrelations. (a) honeycomb U=8
Figure9examinesthesignproblemfortherepulsivecase. β=4
0.8
No qualitative difference is discernible from the square
lattice. Indeed, the sign shows no signal whatsoever as
it passes through ρ = 2/3, the filling which corresponds 0.6
(cid:1)
to entry into the flat band. n
g
The honeycomb lattice is another bipartite lattice we si N=32
(cid:0)0.4 N=50
study here. It has an interesting semi-metallic density
N=72
of states which vanishes linearly at ω 0. Like the N=98
→
hypercubic lattices, it has N =N . Figure 10 exhibits 0.2 N=128
A B
the average sign in the usual array of panels. S is a N=162
bit reduced in densities ρ . 0.6 compared to thheiother N=200
0.0
bipartite lattices, but otherwise behaves in a manner 0.0 0.2 0.4 0.6 0.8 1.0
ρ
rather similar to them.
Our final bipartite geometry is a 1/5-depleted square 0
(b) U=8 (c)
lattice. This is a cousinofthe Lieb lattice, inthat it can
beregardedasasquarelatticewith1/4(ratherthan1/3) n(cid:1)−1 β=4 n(cid:1)−1 U=8
othfethmeasgitneestircemVoavteodm, asnindCisatVhe40g9eo(smeeettrhyearpigphrotppraianteeltoof nsig(cid:0)−2 nsig(cid:0)−2 N=98
Fig. 8). As with the Hubbard model on a Lieb lattice, l ρ=0.625 l ρ=0.625
this model exhibits interesting magnetic orderings. In −3 ρ=0.875 −3 ρ=0.875
particular, at half-filling and U = 6, as the ratio t′/t of 0 60 120 180 0 2 4
the inter- to intra-plaquette hopping is increased, one N β
goes from a plaquette singlet phase to a phase with
antiferromagnetic long-range order at (t′/t)c1 ≈ 0.7 and FIG.10: (Coloronline)SameasFig.1butforthehoneycomb
then to a dimer singlet phase where long-range order geometry and with different parameters, which are noted in
is again absent at (t′/t) 1.3. Figure 11 shows thepanels.
c2 ≈
the doping dependence of the average sign for the 1/5-
1.0
depleted square lattice for t′ = t, which corresponds to
(a)
the ordered phase at ρ=1.
0.8 Figure 12 gives the dependence on t′/t for two fixed
chemical potentials which correspond to ρ 0.9 and
≈
Lieb 0.6.ρ.0.7. Althoughthedensityisvaryingabitwith
0.6
(cid:1) t′/t, there is a steady decrease in S as t′/t decreases.
n U=5.6 h i
g The dimer singlet phase at large t′/t has a well-behaved
i
s(cid:0)0.4 β=6 sign,whiletheplaquettesingletphasehasamuchsmaller
average sign, presumably as a consequence of the fact
that the sign problem of an isolated 2 2 plaquette is
0.2 N=12 N=75 much worse than a dimer. ×
N=27 N=108
N=48
0.0
0.0 0.2 0.4 ρ 0.6 0.8 1.0 C. Triangular and Kagome Lattices
0
(c) We conclude our surveyof lattice geometrieswith two
−1 non-bipartite cases: the triangular and Kagome lattices.
−1 U=5.6
nsign(cid:0)(cid:1)−2 (b) Uβ==65.6 nsign(cid:0)(cid:1)−2 N=108 SpwirenocbmeleutmshtearpteraeislslendnteondsPaitHtiaeSs,o,viwneercluaexdfpiunelgclthraathnlfegfireelloiwnfgifil.lllMibneogrsae,o0sviegrn,
l−3 ρ=0.625 l−3 ρ=0.625 ρ 2. The sign of the hopping is also relevant for thes≤e
ρ=0.875 ρ=0.875 str≤uctures. Ourchoicesaret=1forthetriangularlattice
−40 50 100 0 2 4 6 and t = 1 for the Kagome lattice. In terms of the
N β −
density of states, these choices mean that the DOS is
nonzeroin the ranges[-6,3]and[-2,4]for the triangular
FIG.9: (Coloronline)SameasFig.1butfortheLieblattice and Kagome geometries, respectively.
andwithdifferentparameters, whicharenotedinthepanels. Figure13showsresultsforthetriangularlattice. Shell
9
1.0 1.0
(a) (a) triangular
U=6
0.8 0.8
1/5 depleted β=4
t=t
(cid:2)0.6 ′ (cid:1)0.6
n U=6 n
g g
si β=4 si
(cid:0)0.4 (cid:0)0.4
N=16
0.2 N=64 0.2 N=9
N=144 N=36
N=256 N=81
0.0 0.0
0.0 0.2 0.4 ρ 0.6 0.8 1.0 0.0 0.5 1.0 1.5 2.0
ρ
0 0
U=6 β=4 (c) 0 0
−1 t=t′ U=6 (c)
lnsign(cid:0)(cid:2)−−12 (b)ρρ==00..682755 t=t′ lnsign(cid:0)(cid:2)−−−−4352 UNρρ====0061..68427455 lnsign(cid:0)(cid:1)−−−−4312 (b)ρ=0.875 β=4 lnsign(cid:0)(cid:1)−−−312 UN==ρρ68==100..867255
ρ=1.125 ρ=1.125
0 80 N160 240 0 2 β 4 6 −50 40 80 −40 2 4
N β
FIG. 11: (Color online) Same as Fig. 1 but for the 1/5-
FIG.13: (Color online) SameasFig. 1butforthetriangular
depleted square lattice with equal inter- and intra-plaquette
geometry and with different parameters, which are noted in
hoppingamplitudesanddifferentparametersforU,β andN,
thepanels.
which are noted in thepanels.
structure is evident for the smallest cluster (N = 9),
previous cases. The same is true of the Kagome lattice
however, as with the square lattice, is absent for larger
in Fig. 14, except that there is a persistent bump in S
lattices (N = 36,81). The most marked difference from h i
for ρ slightly less than quarter filling. (A similar feature
thebipartitecasesisthat S isnotpinnedatoneforρ=
was noted for chains.) The structure in the average sign
h i
1,butotherwisethebehaviorof S isquitesimilartothe
issomewhatmoreasymmetricabouthalffilling thanthe
h i
triangular case. One feature which does not seem to be
sharedwithothergeometriesistheexistenceofanabrupt
1.0 change in S appearing here at ρ 1.1, so that it is
h i ≈
well behaved for most densities ρ & 1.1. The scalings
0.8
ρ with the cluster size and β of the average sign for these
, 0.6 1/5 depleted non-bipartite geometries, shown in the bottom panels of
(cid:2)
gn0.4 U=6 N=144 Figs. 13 and 14, follow a similar pattern as for the other
si bipartite geometries.
(cid:0) µ= 1.5
0.2 β=4 −
µ= 2.5
−
0.0
0.0 0.2 0.4 0.6 0.8 1.0 1 0.8 0.6 0.4 0.2 0
t/t t/t
′ ′
FIG. 12: (Color online) Evolution of the average sign on the IV. FURTHER ANALYSIS
1/5-depletedsquarelatticeastheratiooftheinter-andintra-
plaquettehoppingamplitudes(t′/t)increases. Theresultsare
obtainedfortwodifferentfixedchemicalpotentials,µ=−1.5 We now discuss possible patterns which emerge from
and -2.5. The filled symbols show the average sign and the these datasets. We focus on three areas: the role
empty symbols denote the evolution of the average density of the density of states, the contribution of individual
ρ. Note that ρ varies slowly in most of the range of t′/t and spin components to the sign problem, and spatial
only increases significantly towards one for t/t′ . 0.5. The entanglement. A fourth feature of S , scaling in the
bandwidth is kept fixedat 6 for all valuesof t′/t. h i
vicinity of the PHS point, was discussed previously.
10
1.0 1.0
(a)
0.8 0.8
kagome (cid:0)
↓
t= 1 n ladder
n(cid:2)0.6 U=−6 sig 0.6 ρ=1
g (cid:1)
i (cid:0)
s β=4 ↑ U=4
(cid:0)0.4 n 0.4
g
i N=100
s
(cid:1)
0.2 N=27 0.2 β=4
N=54 β=5
N=108 β=6
0.0 0.0
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0
ρ
t /t
⟂
0 0
(b) U=6 β=4 (c)
−1 U=6 FIG. 15: The product of the expectation values of the signs
lnsign(cid:0)(cid:1)−−−432 ρρ==00..682755 lnsign(cid:0)(cid:1)−−12 N=ρρ1==0800..682755 oahTftohptehphaieslnfi-ggifinnfldosliirnvatgirhdeueasalwaladesldplf-euibnrnegchuetaipoovmneadenotdrfiynt.shp(teinhhrSead↑tobiiwoa=nnodhfmSri↓aunitnsrbugiyclaetstsoyinmiosgnms-pechhhtoraawyis.nne)
−5 ρ=1.125 ρ=1.125 t⊥/t>2.
−3
0 50 100 0 2 4
N β
FIG. 14: (Color online) Same as Fig. 1 but for the Kagome
geometry and with different parameters, which are noted in andhence,roughlyspeaking,atleasthalfthebandwidth
thepanels. W. These sorts of interaction strengths are the ones
typicallystudiedinexaminingmagnetic,pair,andcharge
correlations in Hubbard models. The conclusion of the
A. Role of the U =0 Density of States observations above seems to be that U & W/2 brings
the system far enough from the U = 0 limit that
It seems plausible that the noninteracting density of most features of the noninteracting DOS are no longer a
states could play an important role in the sign problem. controlling factors in the sign problem. Apparently, the
In this subsection we make a few observations on that space and imaginary time fluctuations of the Hubbard-
possibility. Stratonovich field (i,t), whose effect on the fermions
X
The square and honeycomb lattices have quite increases with U through the parameter λ of Eq. 4,
dramaticallydifferentU =0densitiesofstates,especially smear out effects of the U = 0 energy levels on S .
h i
near ω = 0 where the DOS diverges logarithmically ThisappearsconsistentwithacomparisonofFig.3(a)of
for the square lattices, and vanishes linearly for the this paper with Fig. 10(a) in Ref. 10 for 4 4 Hubbard
×
honeycomb lattice. Yet, if we compare the behaviors of lattices. With βU constant, the sharpness of the feature
S as a function of filling in Figs.3 and 10, we see little near ρ=0.6 is significantly reduced when U is increased
h i
qualitative difference. Both evolutions exhibit a rapid from 4 to 6.
fall-off from S = 1 at the PHS point ρ = 1, a broad
h i
minimum centered at ρ 0.8, followed by a recovery to This lack of dependence on the DOS is the case even
∼
S = 1 in the dilute limit. The differences in DOS are for a flat band, where the very large delta function in
h i
even more diverse among the other geometries studied the DOSmighthavebeenexpectedtohaveanespecially
here. However, the special features of the DOS, which discernible impact on S . However, in the case of the
h i
could trend in completely opposite directions, appear to Lieb lattice (Fig. 9), S behaves completely smoothly
h i
have little to no effect on the behavior of S here. through the edge of the flat band at ρ = 2/3. We do
h i
Indeed, we have noted already that, more generally, note that broadfeatures in the U =0 DOSdo appearto
the different geometries and their associated distinct have some correlation with the behavior of S = S .
↑ ↓
h i h i
densities of states all share a qualitatively similar Forexample,forthechain,square,andcubiclattices,we
behaviorof S withdoping. Theonly‘unique’geometry observethat S tendstobesmallerwhena“smoothed”
σ
h i h i
wastheone-dimensionalcasewhere,forexample S did DOS is larger (not shown). Behavior near half-filling for
h i
not recover to one at low densities. bipartite lattices is additionally mediated, as discussed
Most of the data we present are for interaction previously, by the fact that at half-filling S = 1 by
h i
strengths U at least four times the fermion hopping t, symmetry.