PhysRevX 5 041041

PHYSICAL REVIEW X 5, 041041 (2015)
Solutions of the Two-Dimensional Hubbard Model: Benchmarks and Results from

a Wide Range of Numerical Algorithms
J. P. F. LeBlanc,1 Andrey E. Antipov,1 Federico Becca,2 Ireneusz W. Bulik,3 Garnet Kin-Lic Chan,4 Chia-Min Chung,5
Youjin Deng,6 Michel Ferrero,7 Thomas M. Henderson,3,8 Carlos A. Jiménez-Hoyos,3 E. Kozik,9 Xuan-Wen Liu,6
Andrew J. Millis,10 N. V. Prokof’ev,11,12 Mingpu Qin,13 Gustavo E. Scuseria,3,8 Hao Shi,13 B. V. Svistunov,11,12
Luca F. Tocchio,2 I. S. Tupitsyn,11 Steven R. White,5 Shiwei Zhang,13 Bo-Xiao Zheng,4
Zhenyue Zhu,5 and Emanuel Gull1,*
(Simons Collaboration on the Many-Electron Problem)

1
Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA
2
CNR-IOM-Democritos National Simulation Centre and International School for Advanced Studies
(SISSA), Via Bonomea 265, I-34136 Trieste, Italy
3
Department of Chemistry, Rice University, Houston, Texas 77005, USA
4
Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA
5
Department of Physics and Astronomy, University of California—Irvine, Irvine, California 92617, USA
6
Hefei National Laboratory for Physical Sciences at Microscale and Department of Modern Physics,
University of Science and Technology of China, Hefei, Anhui 230026, China
7
Centre de Physique Theorique, Ecole Polytechnique, CNRS, 91128 Palaiseau Cedex, France
8
Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA
9
Department of Physics, King’s College London, Strand, London WC2R 2LS, United Kingdom
10
Department of Physics, Columbia University, New York, New York 10027, USA
11
Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003, USA
12
Russian Research Center “Kurchatov Institute,” 123182 Moscow, Russia
13
Department of Physics, College of William and Mary, Williamsburg, Virginia 23187, USA
(Received 9 May 2015; revised manuscript received 30 October 2015; published 14 December 2015)
Numerical results for ground-state and excited-state properties (energies, double occupancies, and
Matsubara-axis self-energies) of the single-orbital Hubbard model on a two-dimensional square lattice are
presented, in order to provide an assessment of our ability to compute accurate results in the
thermodynamic limit. Many methods are employed, including auxiliary-field quantum Monte Carlo, bare
and bold-line diagrammatic Monte Carlo, method of dual fermions, density matrix embedding theory,
density matrix renormalization group, dynamical cluster approximation, diffusion Monte Carlo within a
fixed-node approximation, unrestricted coupled cluster theory, and multireference projected Hartree-Fock
methods. Comparison of results obtained by different methods allows for the identification of uncertainties
and systematic errors. The importance of extrapolation to converged thermodynamic-limit values is
emphasized. Cases where agreement between different methods is obtained establish benchmark results
that may be useful in the validation of new approaches and the improvement of existing methods.
DOI: 10.1103/PhysRevX.5.041041 Subject Areas: Computational Physics,

Condensed Matter Physics,
Strongly Correlated Materials
I. INTRODUCTION systems of large numbers of interacting electrons is one

of the grand scientific challenges of the present day.
The “many-electron problem” of providing a useful and
Improved solutions are needed both for the practical
sufficiently accurate calculation of the properties of
problems of materials science and chemistry and for the
basic science questions of determining the qualitative
behaviors of interacting quantum systems.
*
[email protected] While many problems of implementation arise, including
calculation of the multiplicity of orbitals and interaction
Published by the American Physical Society under the terms of
the Creative Commons Attribution 3.0 License. Further distri- matrix elements needed to characterize real materials, the
bution of this work must maintain attribution to the author(s) and fundamental difficulties are that the dimension of the
the published article’s title, journal citation, and DOI. Hilbert space needed to describe an interacting electron
2160-3308=15=5(4)=041041(28) 041041-1 Published by the American Physical Society

J. P. F. LEBLANC et al. PHYS. REV. X 5, 041041 (2015)
system grows exponentially in the system size so that direct This improved understanding of the Hubbard model will
diagonalization is not practical except for small systems serve as a tool to analyze methods for solving the general
and that the minus sign associated with the Fermi statistics many-electron problem.
of electrons leads to exponentially slow convergence of We take the view that the only meaningful points of
straightforward Monte Carlo calculations. It is generally comparison between methods are results for the actual
accepted that a complete solution to the many-electron thermodynamic limit of the Hubbard model, with the
problem cannot be obtained in polynomial time. uncertainties arising from the needed extrapolations (to
The difficulties associated with obtaining a complete infinite system size, to all diagrams, to infinite statistical
solution have motivated the development over the years precision in Monte Carlo calculations, etc.) quantified.
of approximate methods, and comparison of the different However, examination of results obtained at different
approximations remains a crucial open question. In this stages of the extrapolation sequence for a given method
paper, we address this issue in the context of the provide considerable information. Therefore, we present,
repulsive-interaction Hubbard model [1–3] defined on when needed, both converged values and the intermediate
a two-dimensional (2D) square lattice. The Hubbard results from which these were obtained.
model is one of the simplest models of interacting The methods considered are auxiliary-field quantum
fermions, but despite its simplicity, it exhibits a wide Monte Carlo (AFQMC) [21–23], bare and bold-line dia-
range of correlated electron behavior including interac- grammatic Monte Carlo (DiagMC) [24–26], the dual
tion-driven metal-insulator transitions, superconductiv- fermion method (DF) [27], density matrix embedding
ity, and magnetism. The precise behavior depends theory (DMET) [28,29], density matrix renormalization
delicately on parameters, creating an interesting chal- group theory (DMRG) [30,31], cluster dynamical mean-
lenge for theory and computation. field theory in the dynamical cluster approximation (DCA)
Exact solutions are available for one-dimensional [4] and [32], diffusion Monte Carlo based on a fixed-node approxi-
infinite-dimensional cases [5,6]. High-temperature series mation (FN) [33–39], unrestricted coupled cluster theory
expansions provide numerically exact results but only for including singles and doubles (UCCSD), and in certain
temperatures too high to be relevant for physically inter- cases, higher excitations [40–42], and multireference pro-
esting situations [7]. In general dimensions at relevant jected Hartree-Fock (MRPHF) [43,44].
temperatures, only approximate solutions are available. In We examine energies and double occupancies, which are
some cases, these provide rigorous bounds to the ground- single numbers and can be obtained by essentially all
state energies or other thermodynamic properties [8,9]. methods, enabling straightforward comparison. We also
consider properties related to the electron Green’s function,
Analytical perturbative methods can provide information
which, at this stage, are only available from a few methods.
about the behavior at very small interaction strengths
The results obtained from different techniques enable us to
[10–16] and at very large interaction strengths (for the
identify regimes of phase space that are well understood, in
special case of nearly one electron per site) [17,18], but
the sense that several different methods provide results that
outside of these limits, obtaining results requires numerical are converged and agree within (reasonable) errors, and
methods [19]. Other techniques such as diagrammatic regimes that are not well understood, in the sense that there
resummation are expected to work well in the weak is as yet no agreement between different methods. We show
coupling regime [20]. The known numerical methods are excellent agreement and small uncertainties between
based on approximation schemes. Among the approxima- numerically exact techniques at half filling (all coupling
tions employed are the study of finite systems (either strengths), weak coupling (all carrier concentrations),
directly or via embedding constructions), use of variational and for carrier concentrations far from half filling (most
wave functions, and evaluation of subsets of all possible interaction strengths). For carrier concentrations near half
Feynman diagrams. Controlling these approximations and filling and for intermediate interaction strengths, results can
assessing the remaining uncertainties is a challenging but be obtained, but the resulting uncertainties are much larger,
essential task, requiring analysis of results obtained from in general, and more difficult to eliminate. We surmise that
different methods. The past decade has seen the develop- at least part of this uncertainty has a physical origin related
ment of interesting new methods and substantial improve- to the presence of several competing phases, leading to
ments in capabilities of previously developed approaches, sensitive dependence on parameters.
suggesting that the time is ripe for a careful comparison. The rest of this paper is organized as follows. In Sec. II,
In this paper, we undertake this needed task. Our aim we define the model, delineate parameter regimes, and
is to assess the state of our knowledge of the Hubbard define and discuss the observables of interest in this paper.
model, identifying parameter regimes where reliable In Sec. III, we summarize the methods, giving brief
results have been established and regimes where further descriptions and focusing on issues most relevant to this
work is needed. In regimes where reliable results have paper while referring the reader to the literature for detailed
been established, our results will serve as benchmarks to descriptions. Section IVA demonstrates the importance of
aid in the development and validation of new methods. the extrapolation of results to the thermodynamic limit and
041041-2
SOLUTIONS OF THE TWO-DIMENSIONAL HUBBARD … PHYS. REV. X 5, 041041 (2015)
TABLE I. Parameters studied for the Hubbard model. Here, t III. METHODS
denotes the nearest-neighbor hopping, t0 the next-nearest-neigh-
bor hopping, U the interaction strength, n the density, and T the A. Overview
temperature. Many numerical methods provide solutions to the
Hubbard Hamiltonian on a finite-size lattice. In this work,
Parameters Parameters studied
we restrict attention to techniques that can access systems
t0 =t −0.2 0 0.2 large enough that an extrapolation to the infinite system can
U=t 2 4 6 8 12 be performed, and our aim is to obtain results for the
n 1.0 0.875 0.8 0.6 0.3 thermodynamic limit. Also important to our analysis is the
T=t 0 1=8 1=4 1=2 assessment of uncertainties, either by providing an
unbiased error bar or an error bar that contains all errors
discusses the issues involved in the extrapolations while except for a systematic contribution that may be assessed
Sec. IV B summarizes sources of uncertainty. In Secs. V, VI, by comparison to other methods or reference systems.
and VII we present static observables in the strong coupling, We consider three broad classes of methods: ground-state
intermediate coupling, and weak coupling regimes, respec- wave function methods, embedding methods, and Green’s
tively. Section VIII presents momentum and frequency function methods. The distinction between these methods is
dependence at finite T. A conclusion summarizes the work, not sharp; several of the algorithms fit into multiple catego-
outlines the important areas where our present-day knowl- ries, but the categorization provides a useful way to organize
edge is inadequate, and indicates directions for future a discussion of the different kinds of uncertainties.
research. Supplementary material presents the thermody- Wave-function methods construct an approximation to
namic-limit values obtained here. A database of results is the ground-state wave function of a system. Expectation
made available electronically [45,46]. values of observables (energies and correlation functions)
are then evaluated by applying operators to this ground-
II. THE HUBBARD MODEL state wave function. The issues are the accuracy of the wave
function for a given system size and the accuracy of the
The Hubbard model is defined by the Hamiltonian
extrapolation to the thermodynamic limit. AFQMC (with
X X and without constrained path CP), UCCSD, FN, MRPHF,
H ¼ − tij ðc†iσ cjσ þ H:c:Þ þ U ni↑ ni↓ ; ð1Þ and DMRG are wave-function methods.
i;j;σ i
Embedding methods approximate properties of the full
where c†iσ (ciσ ) creates (annihilates) an electron with spin system (for example, the self-energy or density matrix) by the
σ ¼ ↑, ↓ on site i, niσ ¼ c†iσ ciσ is the number operator, and solution of a finite cluster self-consistently embedded in an
tij denotes the hopping term. In this work, we restrict appropriately designed infinite lattice. The full solution of the
our discussion to the repulsive Hubbard model (U > 0) original problem is recovered as the cluster size is taken to
defined on a two-dimensional square lattice; we further infinity. Errors in embedding methods arise from three
assume that the hopping contains only nearest-neighbor sources: the solution of the cluster problem, the convergence
(tij ¼ t) and second-nearest-neighbor (tij ¼ t0 ) terms. of the self-consistency loop that performs the embedding,
Here, t is used to set the scale of all energies presented and finite cluster size, with maximum cluster sizes depending
in this work. on the method and ranging from 16 to 100. The DMET,
We consider interaction strengths ranging from U=t ¼ 2 DCA, and DF are embedding methods.
to U=t ¼ 12 and focus on temperatures where high- Green’s function methods are defined here as methods
temperature expansion methods fail [7]. We give represen- based on stochastic evaluations of many-body perturbation
tative results for the ground state and temperatures series. They provide many-body self-energies Σðk; iωn Þ and
T=t ¼ 0.125, 0.25, and 0.5. Table I contains a complete Green’s functions Gðk; iωn Þ, as functions of momenta k and
list of the parameters studied. fermionic Matsubara frequency iωn, from which other
At zero temperature, we systematically compute the observables (energies and densities) can be calculated.
energy per site and the double occupancy, and we also Techniques that produce real-frequency (rather than
present some data on the nature of the order and the order Matsubara frequency) information [47–50] are either
parameter, where an ordered phase is found. At nonzero restricted to small systems, molecules or impurity models,
temperature, we also present dynamical information, in or work best at weak coupling. DiagMC is a Green’s function
particular, the Matsubara self-energy. We focus on values at method formulated directly in the infinite lattice; the main
a given density, rather than at a given chemical potential. issues for this method are the accuracy of the stochastic
This implies that methods based on a grand-canonical approximation to the full diagrammatic expansion and the
formulation need to adjust the chemical potential to find the systematic truncation and extrapolation of the series. The DF
right density, leading to additional uncertainty in computed and DCA techniques use Green’s function techniques to
quantities. evaluate the impurity problem and the expansion around it;
041041-3
they therefore are subject both to embedding uncertainties reported). The systematic error from the constrained-path
and to the uncertainties arising from the evaluation of the method is not variational but depends on jΨT i in the sense
diagrammatic expansion. that it vanishes if jΨT i is exact. Its magnitude will be
quantified below by comparison to other techniques.
B. Auxiliary-field quantum Monte Carlo A Trotter decomposition is used in the imaginary-time
evolution. The Trotter error from the finite time step,
The AFQMC method is a ground-state wave-function Δτ ¼ β=n, must be extrapolated to zero. This extrapolation
method based on the idea that in the limit β → ∞ the can be controlled and does not make a significant con-
operator e−βH applied to an initial wave function jψ ðβ¼0Þ i tribution to the error budget. Most calculations reported
projects out the ground state of the Hamiltonian H. The here use a time step fixed at Δτ ¼ 0.01 in units of t.
projection is formulated as an imaginary-time path integral Results obtained for finite L are averaged over twist
that is stochastically evaluated with the help of auxiliary
angle Θ to remove one-body finite-size effects. For small
fields introduced by a Hubbard-Stratonovich transforma-
systems, a large number of Θ values are needed (about 200
tion. The method is applied to finite-size lattices, and an
for 4 × 4 or 6 × 6) [59,61,62], but for larger systems, far
extrapolation to the infinite lattice case is required. If the
fewer Θ are required to reach the same level of accuracy
calculation is converged, the exact ground-state energy and
(for an L × L system with L ¼ 20, averaging over five
wave function for the lattice are obtained. The issues are the
convergence of the stochastic evaluation of the projector twist angles is sufficient). These results are then extrapo-
and, when particle-hole symmetry is broken, the presence lated to L → ∞. The extrapolation requires care because
of a sign problem. The sign problem is managed using a the ground state depends on the system geometry. For
constrained-path approximation, which introduces a poten- n ¼ 0.875, we used rectangular supercells (mostly 8 × 32,
tial systematic error that must be quantified by comparison checked with sizes 8 × 16, 16 × 16, and 8 × 48 for con-
to other methods. For an introduction to the basics of sistency) to accommodate spin- and charge-density wave
AFQMC methods, see, e.g., Ref. [21]. orders. The extrapolation also requires careful attention to
In this manuscript, we present results obtained from two the functional form of the leading finite volume correction
ground-state AFQMC methods. The first [22] is based on [62–66]. Our final results at the thermodynamic limit
the ground-state path integral form of AFQMC [51–53] but include all statistical errors and a conservative estimate
introduces several new algorithmic ingredients including an of the uncertainty resulting from the extrapolation of
acceleration technique [54] (with force bias [21,55]) in the L → ∞ in order to remove the two-body finite-size effects.
Metropolis sampling and control of Monte Carlo variance (The fit includes 1=L3 and 1=L4 terms [61–63] for the
divergence [22]. This approach is applied to systems at half energies and 1=L for magnetization m2 .) To provide a
filling with t0 ¼ 0, where the sign problem [56] is absent concrete example, at n ¼ 1, t0 ¼ 0, T ¼ 0, U=t ¼ 4, the
because of particle-hole symmetry [57]. The algorithmic energy per site for a 20 × 20 system after Θ averaging is
improvements allow us to reach longer imaginary time in E=t ¼ −0.86038 3 × 10−5 . After a weighted least-
the calculations, achieve a higher acceptance ratio, and square fitting with L ¼ 4; 6; …; 18; 20, the final result in
greatly reduce the Monte Carlo variance [22]. The second the thermodynamic limit is E=t ¼ −0.8603 2 × 10−4 .
approach we employ, to treat cases where the sign problem For n ≠ 1 or t0 ≠ 0, a sign problem appears. The sign
does occur, is referred to as the constrained-path problem makes it impossible to converge the ground-state
Monte Carlo method [58,59]. This approach controls the projection method for the system sizes and propagation
sign problem with a constraint (implemented via a choice lengths β needed, and an alternative method, the CP
of trial wave function [21] jΨT i) on the paths in auxiliary- approximation, is used. The results reported in this paper
field space, which allows stable calculations for arbitrarily follow Ref. [21], using a trial wave function jΨT i to apply a
long imaginary time and system size. boundary or gauge condition on the paths that are included
Both methods obtain the ground state of the Hamiltonian in the path integral in auxiliary field space. All results
for a supercell of size L × L under twist-averaged boundary reported here used single-determinant jΨT i with no release.
conditions [60,61]. The ground state is obtained by use of In these calculations, jΨT i is taken to be a mean-field
Monte Carlo methods to estimate jψ ðβÞ i ¼ e−βH ψ ðβ¼0Þ i. solution for the Hubbard model for given U, L, and Θ with
The total projection length is typically β ¼ 64 in the ground- a U value U eff ¼ minfU; 4tg. The order parameter in the
state projection method, although test calculations with mean-field solution is chosen to be orthogonal to the spin
imaginary-time lengths several times larger were performed. quantization axis. This choice is found to help preserve
The convergence error from finite values of β is negligible. In symmetry in jΨT i, improving the CP results [23,67]. The
the constrained-path method, the runs are open ended and accuracy of the constrained-path approach has been exten-
tend to correspond to much larger values of β. In both sively benchmarked [23,61,62,67]. We have carried out
the ground-state projection and constrained-path methods, additional comparisons with exact diagonalization on 4 ×
the statistical error from the Monte Carlo calculation can 4 systems. At U ¼ 4t, the relative error on the energy,
be reliably estimated (1-standard-deviation error bars are averaged over 60 randomly chosen Θ values, is þ0.018% for
041041-4
n ¼ 0.25 and −0.15% for n ¼ 0.625. At U=t ¼ 8 and chain using the Metropolis algorithm, where walkers are
n ¼ 0.875, the relative error is −0.51% averaged over 20 defined by many-body configurations having electrons on
random Θ values. We have also verified these estimates in a lattice sites with given spin along the z axis. After performing
few systems of larger L, using multideterminant trial wave this optimization, a further improvement can be obtained by
functions and constraint release [23,67]. applying the Green’s function Monte Carlo projection
technique [38] to the optimal trial state within a fixed-node
C. Fixed-node diffusion Monte Carlo with nodes from approximation [39]. This procedure allows accurate calcu-
variational Monte Carlo lations of the energy and diagonal correlation functions, such
The variational Monte Carlo method constructs a trial as double occupancies or density-density correlation func-
wave function that approximates the exact ground state tions. The Ansatz on the nodal structure given by the
of a correlated Hamiltonian [33–37]. The wave function variational wave function induces a systematic error, which
depends on parameters that are optimized by minimizing cannot be determined a priori but can be estimated from the
the expectation value of the Hamiltonian, which requires a change in energy as the trial wave function is improved. We
Monte Carlo sampling whenever the trial state is correlated point out that there is a difference between continuum and
(i.e., it is not a simple Slater determinant). We remark that lattice fixed-node approaches. In the continuum, only the
the variational Monte Carlo energy gives an upper bound to signs of the trial function are important: If the nodes are
the exact value and that, with this approach, it is possible to correctly placed, the exact energy is obtained. By contrast, on
access quite large clusters, with all relevant spatial sym- the lattice, both the signs and the relative magnitudes of the
metries (translations, rotations, and reflections) preserved. trial function in configurations that are connected by a sign
However, it is difficult to quantify the systematic errors flip must be correct in order to have the exact energy [39]. The
introduced by the choice of the trial state. Moreover, while error bars for finite systems are given as the statistical errors of
spatial correlations may be correctly captured, dynamical the Green’s function Monte Carlo technique and do not
properties are missed. include any estimates of the systematic errors coming from
We generated simple variational wave functions by the fixed-node approximation.
applying a density Jastrow factor on top of uncorrelated The finite-size results are then extrapolated to the infinite
states that are built from local (mean-field) Hamiltonians, system size by using a scaling that depends on the carrier
containing only a few parameters, where the physical concentration. At half filling, we use the 1=N 3=2 scaling
properties are reflected in a transparent way as different (where N is the system size) that is appropriate for two-
terms inside the variational state [68,69]. At finite doping, the dimensional ordered antiferromagnets [64,65]. In this
uncorrelated states have been obtained from the Bardeen- regime, the error bar for the infinite system size is given
Copper-Schrieffer (BCS) Hamiltonian, including supercon- by a fitting error of the linear regression. At finite dopings, the
ducting pairing on top of electron hopping; in addition, size scaling may suffer from shell effects: A smooth behavior
collinear antiferromagnetism with Néel order parallel to the z can be obtained only when a sequence of closed-shell
spin quantization axis is also included. At half filling, where configurations are taken (i.e., electron fillings for which
the system exhibits long-range magnetic order, the uncorre- the noninteracting case corresponds to a unique ground
lated state contains only magnetism in the x-y plane; in this state). In the generic case, size effects may be dominated by
case, an additional spin Jastrow factor involving the z the ones present at U ¼ 0. This is the case for large dopings
component of the spin operator is also taken. This term (e.g., n ¼ 0.8) and all interactions, and intermediate dopings
couples spins along a direction orthogonal to the magnetic (e.g., n ¼ 0.875) and small interactions (U ≲ 4). Here, for
ordering plane, reproducing the spin-wave fluctuations every available size, the ratio between the energy at finite U
above the mean-field state [70]. All these variational wave and the one at U ¼ 0 is roughly constant and the thermo-
functions with Jastrow factors generalize the so-called
dynamic limit can be assessed by fitting this ratio; namely, the
Gutzwiller state [2], allowing a description of metals, super-
infinite-size energy is obtained by multiplying the afore-
conductors, and also Mott insulators [68,69]. Nevertheless,
they do not give an accurate approximation to the exact mentioned ratio by the thermodynamic value at U ¼ 0. The
ground state in two spatial dimensions, especially close to extrapolated value is assumed to be normally distributed with
half filling. We obtain a substantial improvement by includ- an error bar taken as the difference between the estimated
ing the backflow correlations inside the uncorrelated state. thermodynamic limit and the largest available size. For
On the lattice, this corresponds to a redefinition of the single- intermediate dopings (n ¼ 0.875), the size scaling starts to
particle orbitals, and it is particularly efficient at strong deviate from the U ¼ 0 case around U=t ¼ 4, and we
coupling [37,71]. decided to take the point at the largest size as the infinite
To determine the variational parameters (i.e., the ones size limit, with an uncertainty of twice the difference to the
related to the Jastrow factors, the ones included in the mean- next-lowest system size. We remark that, in this case, a linear
field Hamiltonian, and also those of the backflow correla- regression with the 1=N scaling gives an estimate of the
tions), we minimize the expectation value of the Hamiltonian. thermodynamic limit that is compatible (within one error bar)
This minimization is performed by constructing a Markov with the point at the largest size.
041041-5
D. Multireference projected Hartree-Fock method us to keep the number of determinants roughly constant, so
The multireference projected Hartree-Fock method the quality of the solution decreases as the system size is
[43,44,72] is a ground-state wave-function approach based increased, and, consequently, the energy increases, pre-
on a trial wave function jΨi characterized by the quantum cluding a thermodynamic-limit extrapolation. We have
numbers Θ, K that is constructed out of a set of broken- used expansions with 4, 24, and 32 determinants for
symmetry Hartree-Fock wave functions (fjDi ig) via pro- half-filled 4 × 4, 6 × 6, and 8 × 8 lattices, respectively.
jection operators. The idea is that a broken-symmetry In the lightly doped regime, larger expansions have been
determinant includes the dominant correlation physics, used: 48 determinants in a 10 × 4 (hni ¼ 0.8) lattice and 80
while the projection restores the physical symmetries. determinants in a 16 × 4 (hni ¼ 0.875) one. A linear
The wave function takes the form extrapolation in the reciprocal of the number of determi-
nants, i.e., in 1=n, has been performed to the infinite
X
n X configuration limit for the ground-state energies.
Θ
jΨΘ
Ki ¼ K 0 P̂KK 0 jD i;
f i;Θ i
ð2Þ The calculations presented could, in principle, be
i¼1 K0 improved in a number of ways: Additional symmetries
could be broken and restored (such as Ŝz or particle
and the parameters f i;Θ
K 0 are determined by minimizing the number) in the reference configurations, more configura-
energy. The projector P̂ΘKK 0 restores the symmetries (char- tions could be included, and/or a full optimization of all
acterized by the quantum numbers Θ, K) in jΨi and can be determinants could be performed, in the spirit of the
formally written as [73,74] resonating Hartree-Fock approach [75]. The accuracy of
any one result can therefore be increased, as shown by
h X Θ Rodríguez-Guzmán et al. [76] or by Mizusaki and Imada
P̂Θ
KK 0 ¼ Γ 0 ðmÞR̂ðmÞ ð3Þ
L m KK [77] in the closely related path integral renormalization
group (PIRG) approach. Recently, Tahara and Imada [78]
in terms of the rotation operators R̂ðmÞ and the irreducible have combined the symmetry-projected determinant expan-
representation matrices ΓΘ ðmÞ associated with the elements sion with short- and long-range Jastrow factors within a
m of the symmetry group of the problem. Here, h is the variational Monte Carlo framework, which may be used to
dimension of the irreducible representation, while L is the further increase the accuracy. These techniques, and others,
volume of the group. The character of the broken- have been used to explore, for example, spin and charge
symmetry determinant is optimized in the presence of stripe phases, which we do not explore in this work
the projection operator (i.e., a variation-after-projection [79–81].
approach), which results in broken-symmetry determinants
with well-defined defects [44]. E. Unrestricted coupled cluster—singles and doubles
Our expansion employs Slater determinants that break Coupled cluster (CC) theory [40–42] is a ground-state
the space group and spin (Ŝ2 ) symmetries of the lattice but wave-function technique. It is widely used in quantum
preserve Sz symmetry. All the broken symmetries are chemistry and often considered the best source of precise
restored using the appropriate projection operators. The data for molecules that are neither too large nor too strongly
series of i determinants in Eq. (2) is constructed through a correlated. Its application to the Hubbard model has been
chain of variational calculations, using the few-determinant more limited, where it has been used in two different
(FED) approach [43,44,74]. In this procedure, after a wave formulations. In the first form, which directly exploits the
function with n − 1 intrinsic determinants is available, a translational invariance of the lattice to work in the
wave function with n determinants is variationally opti- thermodynamic limit, the theory is formulated in the site
mized by adjusting the Thouless coefficients determining basis, starting from an infinite Néel-ordered reference, from
the last-added determinant. The full set of linear coeffi- which clusters of excitations that change occupancy and
cients f iΘ
K 0 is readjusted. The MRPHF approach becomes flip spins are created [82,83]. In the second form, the theory
exact as the number of determinants retained tends to the starts from a single-determinant reference state on a finite
size of the Hilbert space, i.e., as Eq. (2) becomes a lattice, from which clusters of particle-hole excitations are
coherent-state representation of the exact ground-state created [84]. This is similar to how the theory is used in
wave function. quantum chemistry and is the formulation discussed further
If the number of determinants is fixed at a finite, not-too- here. As this second form does not work in the thermo-
large value, these calculations can be performed for large dynamic limit, the energies must be extrapolated.
lattices with a polynomial cost [OðNÞ4 or so, where N is the The CC wave function is written as jΨi ¼ expðTÞj0i,
number of sites in the lattice]. However, the number of wherePj0i is a single-determinant
P reference state, and
determinants required to obtain results with a given T ¼ n T n , where T n ¼ qn tqn A†qn is the cluster operator.
accuracy increases exponentially with increasing system The operator A†qn creates an excited determinant jΦqn i
size. In this work, implementation aspects have compelled which contains n particle-hole pairs relative to the reference
041041-6
state. In its standard and simplest version (used here), the thermodynamic limit even for restricted excitations nmax
energy and coefficients tqn are obtained by solving the [40–42]. Thus, as the lattice size increases, the energy per
Schrödinger equation projectively: site approaches the thermodynamic limit for the given nmax
and the exact thermodynamic limit as nmax is increased. For
E ¼ h0je−T HeT j0i; ð4aÞ smaller values of U where convergence to the thermody-
namic limit is slower, we have converged second-order
0 ¼ hΦqn je−T HeT j0i ∀ qn : ð4bÞ perturbation theory out to the thermodynamic limit and
added a correction for the difference between CC theory
CC theory thus diagonalizes a similarity-transformed and perturbation theory, which, for small U, converges
Hamiltonian H̄ ¼ expð−TÞH expðTÞ in a subspace of quickly with respect to the system size. Double occupan-
states defined by j0i and jΦqn i. Note that because T is a cies have been computed by numerical differentiation of the
pure excitation operator, the commutator expansion used to CC energy with respect to U. We also provide UCCSDTQ
evaluate H̄ truncates after four commutators. estimates of the ground-state energy (labeled as
If the sum over n in defining the cluster operator is UCCSDTQ*) for n ¼ 1 systems, where the triples (T)
carried out to all orders, the exact ground-state wave correction is obtained as UCCSDT-UCCSD energies for a
function is reproduced. In practice, the operator T is 6 × 6 system, and the quadruples (Q) correction is obtained
typically restricted to terms involving a small number of from UCCSDTQ-UCCSDT energies for a 4 × 4 lattice. For
particle-hole pair excitations above the reference state n ¼ 0.8 (n ¼ 0.875), our UCCSDT estimates (labeled as
(n ≤ nmax ), and jΨi ¼ expðTÞj0i is projected onto the UCCSDT*) are obtained from 10 × 4 (16 × 4) lattices.
space of states with up to nmax particle-hole pairs; the The UCCSD calculations reported here can be com-
accuracy then depends upon nmax and the choice of pleted in a few hours using standard quantum chemistry
reference state j0i. In a lattice model such as the packages [89]. Even at half filling where there are not many
Hubbard model with N sites, the computational cost is, Hartree-Fock solutions to be concerned with, we find large
roughly speaking, proportional to N 2ðnmax þ1Þ . effects due to single excitations, which suggests that the
The calculations reported in this manuscript primarily coupled cluster calculations could be substantially
limit T to n ≤ 2, i.e., to the creation of only singly and improved by optimizing the identity of the reference
doubly excited determinants, giving what we refer to as CC determinant. Similarly, we generally see significant cor-
with single and double excitations (CCSD) [85,86]. For rections due to triple and higher excitations; these can also
select examples, we have included corrections for triple and be computed with standard packages [90]. However, while
occasionally quadruple excitation effects. The accuracy of optimizing the reference determinant and including higher
CC theory, and the need for higher excitations, depends on excitation effects will increase the accuracy of the coupled
how well the reference j0i captures the qualitative physics. cluster calculations, they also increase the cost.
When the reference is accurate, then single and double
excitations may be sufficient; however, when the reference F. Density matrix renormalization group
determinant bears little resemblance to the exact wave The density matrix renormalization group [30] is a
function (as may happen in strongly correlated systems), a variational ground-state wave-function technique. It con-
much higher degree of excitation is required to recover the structs the ground state of a system by diagonalizing the
correct physics. For this reason, the calculations reported Hamiltonian in a finite subspace spanned by an iteratively
here use a symmetry-broken unrestricted Hartree-Fock constructed basis that is optimized via a Schmidt decom-
(UHF) reference determinant because UHF can provide position that minimizes the spatial extent of the quantum-
a better mean-field description, particularly near half filling mechanical entanglement between basis states.
where antiferromagnetic correlations dominate; this defines DMRG is generally believed to be the optimal method
the UCCSD method used here. Note, however, that for finding ground states of one-dimensional lattice models.
particularly for doped systems with large U, there are a In the application to one-dimensional systems, two sources
plethora of nearly degenerate UHF states and finding the of error must, in principle, be controlled. Results for a fixed
best reference for UCCSD is not straightforward. We have system length of L must be converged with respect to basis
prepared the UHF solution following the prescription of size m, and then the converged results must be extrapolated
Ref. [87]. In principle, the deficiencies of the reference to L → ∞. However, for most of the one-dimensional
determinant can be corrected in what is known as Hamiltonians of current interest, the convergence is very
Brueckner CC [88], where the reference determinant is rapid; in practice, large enough m and L are accessible
adjusted to eliminate single excitation effects. These numerically, so extrapolation is not required.
calculations are more computationally demanding, and Application of DMRG to a finite-size 2D system
we have not pursued them here. proceeds by defining an effective one-dimensional prob-
An important virtue of the exponential parametrization lem, to which the standard one-dimensional DMRG is
of the wave function is that the CC energy has a nontrivial applied. Two-dimensional tensor network generalizations
041041-7
of the DMRG ideas have attracted tremendous interest were able to sort out stripe fillings, finding the low-energy
recently, but these methods have not yet produced results states and avoiding metastability. For the 45° rotated
for the two-dimensional Hubbard model that can be lattices, the patterns seem more delicate, and we have
included in the present comparison [91,92]. not yet sorted out the lowest-energy configurations. We
Most current implementations of DMRG require open thus present only the results for standard orientation for
boundary conditions. The canonical method for creating an doped systems.
effective one-dimensional system from a finite-size two- Convergence issues pose more severe problems than in
dimensional one is to impose periodic boundary conditions the standard one-dimensional cases, both because of the
in one direction and open boundary conditions in the intrinsic limits on system size discussed above and because,
other, thereby defining a cylinder of finite length and in some cases, the presence of several metastable states can
finite circumference. One then defines an effective one- cause trouble for extrapolation and can lead to the appear-
dimensional problem by indexing the sites along a one- ance of states that are not the ground state and may be
dimensional path that covers all of the sites on the cylinder important only for finite systems. In the systems studied
and taking the matrix elements of the Hamiltonian in this here, metastability was traced to the presence of physically
basis. The price is that two sites separated by a small different “striped” states for different hole doping in the
distance along the cylinder axis in the physical system are DMRG cylinder.
separated by a distance of order the cylinder circumference We obtain converged results as follows. For each
in the effective one-dimensional model. The effective one- cylinder of a specific length L, we extrapolate the energy
dimensional problem thus has long-ranged terms, which in the number of basis states, m. To make this extrapolation
imply longer-ranged entanglement and require that more reliable, m is slowly increased, but each m is used for two
states are kept in the optimal basis. The number of states consecutive sweeps. The truncation error and the energy are
needed grows exponentially with the circumference of the measured on the second sweep for each m. Then, a linear
extrapolation of energy versus truncation error is used to
cylinder (width of the original finite lattice), meaning that
obtain the ground-state energy with error bars. The deter-
there is a sharp cutoff in the accessible system widths,
ministic nature of DMRG can result in uncertainties due to
typically around width 6 in the Hubbard systems; however,
fitting which do not appropriately represent the uncertainty
systems of very large cylinder length can be studied. The
in the choice of extrapolation procedure. We therefore
extrapolation to the thermodynamic limit must thus be
assume a normally distributed error of 1=5 the difference
handled with care. between the last point and the extrapolated value, which we
The DMRG calculations reported here were performed justify by comparison to the accuracy of previous DMRG
with the standard DMRG finite-system algorithm [30,31]. data [93,94]. Metastability is signaled by lack of linearity in
Two types of cylinders were considered: one with the axial this extrapolation. When this is found, we determine which
and circumferential directions aligned parallel to the bonds
state had the lowest energy. The system is then rerun with
of the square lattice and one rotated by 45°. When one cuts
an initial state favoring the lower-energy state producing
a cylinder in two, the 45° rotated system has fewer sites on
results in the lower-energy configuration such that the
the boundary per unit length, and thus one expects a smaller
extrapolations are linear in the truncation error.
growth of the entanglement with the length of the cut
For a fixed width, we then extrapolate the energies
(governed by the area law). For the undoped antiferro-
linearly in 1=L to get an energy per site for an infinite
magnetic system, the rotated system is also unfrustrated,
cylinder. Errors are estimated statistically using the error
both for odd and even circumferences, reducing shell
bars on each point. To reduce the finite-size effects from
effects in the finite-size results. For half-filled systems,
different widths, we employ a simple version of the phase-
both types of orientations are considered, and they show
averaging trick [60] by taking an average over periodic and
good agreement; error bars are estimated to incorporate
antiperiodic boundary conditions around the cylinder. The
both the error bars on the data points for specific widths and
phase averaging eliminates an oscillation in the energy as a
the differences between the two orientations. However, we
function of system width in the 45° rotated systems and, in
note that for the half-filled systems, better results were
general, provides a nicely accelerated convergence. We
obtained with the rotated system [93,94]. With doped
then analyzed the results as a function of cylinder width.
systems, we see striping behavior at stronger coupling.
For half filling, we found that extrapolating the energy to
A stripe is a line of holes that act as a domain wall in the
the thermodynamic limit by 1=width3 works well. This is
antiferromagnetism on either side. These stripes have been
the finite-size behavior in the Heisenberg model, where it is
seen in the t-J model with DMRG starting in the late 1990s
well understood [95].
[66]. The stripes typically wrap around the cylinder, with a
specific even number of holes in the ring stripe. With
doping, the optimal number of holes can change. Striped G. Density matrix embedding theory
configurations with the wrong number of holes in a stripe The DMET [28,29] is a ground-state embedding method
are typically metastable. For the ordinary orientation, we formulated in terms of wave-function entanglement. Given
041041-8
an impurity cluster of size N c , DMET maps an L × L bulk renormalized states m in DMRG. This error is only nonzero
many-body problem (L is chosen to be very large) to an in clusters larger than 2 × 2. The energy and observables
impurity plus bath problem with N c correlated (impurity) are extrapolated to m ¼ ∞ using the standard linear
sites and the same number of bath orbitals. The mapping is relation between energy and DMRG truncation error
constructed from the Schmidt decomposition [96] of a bulk [93,102,105]. For 4 × 4 impurity clusters, the truncation
wave function jΨi. The formulation is exact if jΨi is the error is large enough to also contribute to the converged
exact bulk ground state, or in the limit of impurity cluster DMET self-consistent correlation potential uðmÞ. To take
size N c → ∞. Since the exact bulk state is not known this into account, we extrapolate using (a) self-consistent
a priori, an approximate bulk state is used for the impurity DMET results converged at different m, and (b) non-self-
mapping. Recently, DMET has been applied to both consistent DMET results using different m in the DMRG
ground-state and linear response spectral properties of impurity solver at a fixed correlation potential uðmmax Þ,
the Hubbard model [28,97–99]. In this work, we use a where mmax is the maximum m used in the DMET self-
general BCS bulk state, the ground state of the DMET consistency. The difference between the two extrapolations
lattice Hamiltonian given by the hopping part of the is then added to the total DMET uncertainty. For the cluster
Hubbard Hamiltonian augmented with the DMET corre- sizes in this study, the errors due to (i) and (ii) are small and
lation potential u, which is applied in each cluster supercell easily controlled. Therefore, the finite-size impurity error is
of N c sites in the bulk lattice. This bulk state is allowed to the main source of uncertainty at almost all points in the
spontaneously break spin and particle-number symmetry phase diagram. It is estimated as the standard deviation of
through the self-consistency cycle that determines u (which the finite-size extrapolation. The quality of the approximate
contains both particle-number conserving and nonconserv- DMET impurity mapping depends on the approximate
ing terms). In this cycle, the bulk state jΨi is updated from lattice wave function jΨi and decreases as the coupling
the interacting impurity and bath solution jΦi, by minimiz- strength U increases, especially for carrier concentrations
ing the difference between jΨi and jΦi [as measured by near half filling. This slows down the cluster-size con-
their (generalized) one-particle density matrices] with vergence of the results for large U, increasing the uncer-
respect to the potential u. tainty. In the strong coupling, weakly doped region, we find
For the bulk lattice, we use L ¼ 72. From calculations on competing homogeneous and inhomogeneous orders that
larger L, we find that the finite lattice error associated with become very sensitive to the cluster size and shapes (similar
this choice is negligible, on the scale of the significant digits to the stripes observed in DMRG). It is difficult to reliably
reported. The BCS bulk state is obtained by solving the extrapolate these results to the thermodynamic limit. As a
lattice spin-unrestricted Bogoliubov-de Gennes equation result, the total DMET uncertainties range from about 10−4 t
[100,101] with the correlation potential u. The impurity at half filling and low densities to a maximum of about
and bath problem is solved in the BCS quasiparticle basis, 10−2 t in the strong coupling, underdoped region. A detailed
with general one-body and two-body interactions that do not description of the methodology and extrapolation proce-
conserve particle number or locality. We use the DMRG [30] dures for the calculations is contained in Ref. [106].
as an impurity solver (adapted from the quantum chemistry
DMRG code BLOCK [102–104]), with a maximum number H. Dynamical cluster approximation
of renormalized states m ¼ 2000 (DMET self-consistency is The DCA [32,107,108] is an embedding technique in
performed up to m ¼ 1200). The DMET lattice and impurity which an approximation to the electron self-energy is
Hamiltonians are augmented with a chemical potential μ, to obtained from the solution of a quantum impurity model
ensure that the relative error in particle number is less than consisting of some number N c of interacting sites coupled
0.05%. We use impurity clusters of dimensions 2 × 2, 4 × 2, to an infinite bath of noninteracting electrons. In applica-
4 × 4, and 8 × 2 at each point in the phase diagram. The tions to the Hubbard model, the interactions in the impurity
energies and observables are then extrapolated to the thermo- model are the interactions of the original problem, while the
dynamic limit using a linear relationship with N −1=2 c , as one-body terms are determined from a self-consistency
appropriate for a nontranslationally invariant cluster embed- condition relating the Green’s functions of the impurity
ding theory. model to those of the lattice. The DCA is a particular
The total DMET uncertainty is estimated by combining generalization to N c > 1 of the “single-site” dynamical
the errors from three sources: (i) convergence of DMET mean-field method [109,110]. Other generalizations have
self-consistency, (ii) solution of the impurity many-body also been introduced [111,112], but results from these
problem using DMRG, and (iii) extrapolation to the limit methods are not considered here. The single-site method
of infinite impurity size. We estimate the self-consistency was motivated by the observation that in an appropriately
error (i) using the difference of the last two DMET self- defined infinite coordination number limit, an exact sol-
consistent iterations. The average self-consistency error is ution of the Hubbard model can be found [5] and can be
below 5 × 10−4 t in the energy for all cluster sizes. The recast in terms of the solution of a single-site impurity
impurity solver error (ii) is from using a finite number of model [109]. It was later understood that generalization to
041041-9
N c > 1 impurity-model representations allows a treatment uncertainty in the extrapolation procedure. We include
of the finite coordination number model that becomes exact both a statistical uncertainty in the extrapolation in 1=N c
as the number of impurity sites N c → ∞. and an additional uncertainty, which we take as half the
The DCA formulation partitions the Brillouin zone into difference of the extrapolated value from the largest N c
N c equal area tiles and approximates the self-energy Σ as a value explored. This gives a robust representation of an
piecewise constant function of momentum taking a differ- extrapolation error, which is larger when finite-size effects
ent value in each tile: are large but that also vanishes as N c increases. As such, the
X error bars for extrapolation of our DCA results to the
Σðk; ωÞ ¼ ϕa ðkÞΣa ðωÞ; ð5Þ thermodynamic limit contain both stochastic and finite-size
a¼1.::N c uncertainties, and values for finite system sizes with
stochastic error bars are provided.
with ϕa ðkÞ ¼ 1 for k ∈ a and ϕa ðkÞ ¼ 0 for k∉a. The
tiles a map directly onto impurity-model single-particle I. Dual fermion ladder approximation
levels, the impurity-model Green’s function and self-energy The dual fermion approach [27,122,123] is a diagram-
are diagonal in the impurity model a basis, and the DCA matic extension of the single-site dynamical mean-field
self-consistency equation theory (DMFT). The DF technique is motivated by the idea
Z that nonlocal corrections to DMFT can be captured by a
d2 k 1 perturbative expansion around a solution of the dynamical
Gimp
a ðωÞ ¼ 2 ω − ε − Σ ðωÞ
ð6Þ
a ð2πÞ k a mean-field equations. In formal terms, the expansion
requires summing a series of diagrams for two and higher
determines the remaining parameters of the impurity particle correlations, with vertices defined in terms of the
model. The self-consistency loop is solved by iteration; fully interacting but reducible vertices of the impurity
an initial guess for the impurity-model parameters produces model. In this regard, the DF technique is similar in spirit to
a set of Σa , which are then used to update the impurity- methods such as the dynamical vertex approximation
model parameters. The loop typically converges well, and (DΓA) and dynamical mean-field extensions of fRG
errors associated with the self-consistency are smaller than (DMF2 RG), which approximate interactions on the level
errors in the solution of the impurity model. of two-particle vertex functions [124–126]. In practical
We obtain results for different N c in the paramagnetic implementations to date, the dual fermion expansion is
phase and extrapolate to the thermodynamic limit by truncated at the two-particle level (higher than two-particle
exploiting the known [32] N −d=2c scaling for momentum- interactions are omitted), and the series of two-particle
averaged quantities in d dimensions and systematically corrections is approximated by a few low-order terms or a
increasing N c [113,114]. ladder resummation. One of its strengths lies in the ability
To solve the N c -site impurity problem, we use a CT- to describe phase transitions of the system by employing
AUX algorithm [115,116] with submatrix updates [117]. resummations of the relevant diagrams [127,128].
Our codes are based on the ALPS libraries [118,119]. The DF results presented here are obtained using the
Selected points have been compared to a CT-INT [120] open-source opendf code [129], starting from a single-
implementation based on the TRIQS libraries [121]. site (N c ¼ 1) dynamical mean-field solution obtained with
In this work, we provide extrapolated DCA results from a continuous-time auxiliary-field (CT-AUX) [115–117]
clusters of sizes N c ¼ 16, 20, 32, 34, 50, 64, 72, and 98, impurity solver. The method is limited by the accuracy
depending on temperature and densities, in order to give a to which reducible impurity vertex functions can be
reliable estimate of the properties of the 2D Hubbard model obtained, which is a polynomial (cubic) complexity prob-
in the thermodynamic limit. lem. Within the approximation of neglecting higher-order
The CT-AUX method is a type of diagrammatic vertices and only considering ladder diagrams, a systematic
Monte Carlo method. In the T > 0, impurity-model con- estimation of deviation from the true interacting system is
text, the diagrammatic series is absolutely convergent, and not possible, and we omit error bars altogether.
the issues discussed below for Diag-MC in the infinite-
lattice context are not important. However, the CT-AUX J. Diagrammatic Monte Carlo method
method is restricted at low T by the existence of a The diagrammatic Monte Carlo method (DiagMC)
Monte Carlo sign problem in the auxiliary field solver. begins from the observation that within standard many-
The sign problem worsens rapidly as U is increased, as T is body perturbation theory, any quantity Q that depends on
decreased, or as N c is increased, and this is particularly some set of arguments y (which may include both con-
evident in the density range n ¼ ½0.8; 1.0Þ. Furthermore, as tinuous components such as frequency and momentum and
temperature is decreased, the length scale of correlations in discrete components such as spin) may be expressed as an
the system increases, resulting in larger finite-size effects. infinite series of terms, each of which involves multidi-
We take a conservative approach to determining the mensional integrals and sums over many variables:
041041-10
∞ XZ
X Z
order, the contributions of diagrams with plus and minus
QðyÞ ¼ … dx1 …dxα Dðα; ξ; x1 ; …; xα ; yÞ: signs tend to cancel. This cancellation is, in fact, respon-
α¼0 ξ
sible for the convergence of the many-body perturbation
ð7Þ theory [26,131]. To manage this issue, it is useful to
Here, the D are known functions defined by the Feynman consider a Monte Carlo process for an approximate series
rules of diagrammatic perturbation theory. The series order in which the maximum diagram order is limited to some
α controls the number of internal integration or summation finite value via a hard or soft cutoff, and then the cutoff is
variables fx1 ; …; xα g, and ξ labels different terms of the increased until convergence is reached.
same order in the series. The quantity Q may be the electron Diagrammatic Monte Carlo techniques can also take
Green’s function G, the self-energy Σ, the screened inter- advantage of known field-theoretical techniques to run the
action W, the polarization operator Π, the pair propagator calculation in a self-consistent mode in which certain
(for contact interaction) Γ and its self-energy ΣΓ, or the infinite series of diagrams are summed and then automati-
current-current or other correlation functions. Basic cally absorbed into the renormalized propagators and
thermodynamic properties (particle, kinetic, and potential interaction lines using Dyson-type equations. One example
energy densities) in the grand canonical ensemble are is the skeleton expansion [131]; another is a “bold”
readily obtained from G and Σ (see Ref. [130]). expansion in perturbative corrections to an analytic partial
The most widely used formulation of perturbation theory resummation [132,133]. This flexibility allows for an
is in terms of Feynman diagrams. In this case, standard additional control over systematic errors coming from
rules relate the graphical representation of a given term in series extrapolation as well as convergence issues—
the series to the corresponding mathematical expressions different schemes should produce consistent final results.
for the D, which are typically given (up to a sign or phase In this work, we employ four complementary techniques:
factor) by aQproduct of functions associated with graph (i) A ½Gð0Þ 2 U scheme based on a Taylor series expan-
lines, D ¼ lines Fline ðxline Þ. In perturbation theory for sion of Σ in powers of U with fixed shifted chemical
particles interacting via pairwise forces, these lines are potential μ~ ¼ μ − Un=2 (see P Refs. [26,134]). The
associated with the interaction potential U and bare single- total electron density n ¼ σ nσ and the chemical
particle propagators Gð0Þ . If we denote the collection of all potential are computed a posteriori (after results are
external and internal variables that allow for a complete extrapolated to the αmax → ∞ limit).
characterization of the diagram (diagram topology and (ii) A G2 W scheme based on skeleton series for Σ and Π
internal and external variables) as ν ≡ ðα; ξ; x1 ; …; xα ; yÞ, in which all lines in the diagram are understood as
then Eq. (7) can be viewed as the weighted fully dressed Green’s functions and screened inter-
P average P over the
configuration space fνg: i.e., Q ¼ ν Dν ≡ ν eiφν jDν j, actions. Self-consistency is implemented by feed-
where the modulus of Dν defines the configuration back loops when G and W are obtained by solving
“weight,” and φν ¼ arg Dν . The basic idea of the algebraic Dyson equations, G−1 ¼ ½Gð0Þ −1 − Σ and
DiagMC method is to use stochastic sampling of the W −1 ¼ U−1 − Π, in momentum-frequency repre-
configuration space to compute Q by treating ∝ jDν j as sentation (see Refs. [135–137] for more details).
the probability density for the contribution of point ν to the (iii) A G2 Γ scheme based on the skeleton series for Σ and
sum. In DiagMC, the diagram order, its topology, and all ΣΓ when all lines in the graph are understood as fully
internal and external variables are treated on equal footing, dressed single-particle (Green’s functions) and two-
and each diagram represents a point, not an integral, in fνg. particle propagators. This compact formulation is
The MC process of generating diagrams with probabilities possible only for a contact interaction potential when
proportional to their weight is based on the conventional the sum of ladder-type diagrams for spin-up and spin-
Markov-chain updating scheme [24–26] implemented down particles has the same functional structure as the
directly in the space of continuous variables. single-particle propagator (see Refs. [132,138,139]).
Since connected Feynman diagrams are formulated Again, self-consistency is implemented by feedback
directly in the asymptotic limit, there is no infinite- loops using Dyson equations, G−1 ¼ ½Gð0Þ −1 − Σ
system-size limit to take. The main numerical issue con- and Γ−1 ¼ ½Γð0Þ −1 − ΣΓ , where Γð0Þ is the sum of
cerns the convergence of the Monte Carlo process, which is bare ladder diagrams.
complicated by the exponential proliferation of the number (iv) A ½Gð0Þ 2 Γð0Þ scheme based on diagrams expressed
of diagrams with perturbation order α. This leads to in terms of bare single and pair propagators with
exponential computational complexity since final results shifted chemical potential μ~ ¼ μ − Un=2; this is
are recovered only after extrapolation to the infinite similar in spirit to the ½Gð0Þ 2 U expansion but with
diagram-order limit. The fermion sign enters the calculation one extra geometrical series (bare ladder diagrams)
in an interesting way: Different diagrams have different being accounted for analytically.
signs (arising from the different orderings of fermion To establish the parameter region where DiagMC
creation and annihilation operators), and at large diagram works, we performed calculations using all four schemes.
041041-11
_
The results were compared to each other and to those 1/√N c
0.5 0.4 0.3 0.2 0.1 0
obtained by DCA. Additional insight was also gained by -0.522
DMET
doing calculations in the atomic limit [140]. We find that 2x2
-0.515
for bare coupling U=t < 4 and temperature T=t > 0.1, all 4x2
8x2
schemes produce consistent results within statistical errors. FN
E/t
-0.523
4x4 -0.52 242
At half filling, n ¼ 1, and U=t ≤ 6, the ðGð0Þ Þ2 U and 162
U/t=8, n=1
ðGð0Þ Þ2 Γð0Þ schemes still produce results consistent with 98
-0.525
those obtained by DCA and the determinant Monte Carlo _
-0.524
E/t
method. In the dilute region (small filling factors), the G2 Γ _ 6√2
5√2
and ðGð0Þ Þ2 Γð0Þ schemes can be applied for larger values of U. DMRG 10
14 16
For dilute systems, our benchmarks include points (U=t ¼ 6, _

4√2
6 12 18 -0.525
n ¼ 0.6), (U=t ¼ 6, n ¼ 0.3), and (U=t ¼ 8, n ¼ 0.3). 8
-0.526
IV. EXTRAPOLATIONS AND UNCERTAINTIES U/t=8, n=1 AFQMC
0.005 0.0025 0
A. Extrapolations
1/L3
All of the numerical methods we have considered rely
on the extrapolation of results to a thermodynamic or FIG. 1. Extrapolations of the ground-state energy at U=t ¼ 8,
asymptotic limit. For DiagMC, which is formulated n ¼ 1. Main panel: AFQMC and FN extrapolated as a function of
directly on an infinite system, the extrapolation is in the inverse cube of the system’s linear dimension L along with
diagrammatic order. All other methods are extrapolated DMRG extrapolated in the cube of the inverse cylinder circum-
from a finite embedded system, finite cluster, or cylinder ference (also pdenoted
ffiffiffi L). DMRG data are presented both for
rotated (with 2) and unrotated wrapping of cylinders. Inset:
with finite width to the infinite-system-size limit. In many
DMET data for clusters of size and geometry indicated, plotted
cases, a considerable contribution to our errors comes against the reciprocal of the square root of the total number of
from this extrapolation procedure, which differs from sites in the cluster N c.
method to method. In some cases, determining stochastic
uncertainties in extrapolation is not possible, in which
case we produce estimates of uncertainties by choosing a For DMRG, results from both unrotated cylinders (filled
reference system for a given technique. We then assume a symbol; other smaller-L data are not shown but lie on the
normal distribution of uncertainty with respect to the same scaling curve) and rotated cylinders (crosses) are
reference. Specifics of the extrapolation procedure for consistent, both scaling as 1=L3 , although with different
each system (and of the associated procedure for estimat- slopes, allowing for a clean extrapolation to the thermo-
ing extrapolation uncertainties) are described in Sec. III. dynamic limit. On the other hand, the AFQMC data
All methods have therefore defined procedures to estimate indicate a change in scaling for system sizes larger than
error in TL quantities as accurately as possible through the 10 × 10 geometry under twist-averaged boundary condi-
use of known reference systems. Additional uncertainties tions. This could mean either that unidentified complica-
due to extrapolation, curve fitting, and truncation in tions occur in the large-system AFQMC calculations or that
excitation order are addressed on a per-technique basis. deviations from the 1=L3 size dependence might occur in
These added uncertainties are assumed to be normally the DMRG data at larger cylinder sizes (i.e., that the
distributed and defined such that they can, in principle, be DMRG error bar is underestimated). In this regard, it is
made arbitrarily small by adding additional data. This important to note that the ground-state energy of the largest
section illustrates these extrapolations and presents some system examined in DMRG, rotated 6 × ∞, is within the
of the challenges encountered in performing them. uncertainty of the extrapolated AFQMC data.
We start our discussion with ground-state properties. In The FN data also demonstrate a systematic dependence
Fig. 1, DMRG, FN, and AFQMC results are presented in of the energy on system size, allowing a precise thermo-
the main panel. The system-size dependence in all three dynamic-limit extrapolation. The deviation of the FN
techniques is clearly visible, and the difference between the results from the AFQMC and DMRG results is caused
estimated thermodynamic-limit value and the largest size by a systematic fixed-node error, which by comparison to
considered is, in some cases, outside of the error bars of the other methods, seems to be no more than 2 × 10−3 t.
thermodynamic-limit value; in other words, extrapolation is Shown in the inset
pffiffiffiffiffiffiof Fig. 1 are the extrapolations in DMET
essential for obtaining the thermodynamic-limit value. For which scale as 1= N c . Because of the restricted small system
this reason, in the results sections, we typically present both sizes in DMET and large U, the resulting uncertainty is
the thermodynamic-limit value and the sequence of finite- dominated by the extrapolation. The value, also shown in
size results which led to it, so that the reader can see how the main panel, is in good agreement with DMRG and FN,
large an extrapolation is required. and only slightly outside error bars of the AFQMC result.
041041-12
2
0.07 0.06 0.05 0.04
1/α
0.03 0.02 0.01 0
(i) AFQMC: At n ¼ 1, error bars include all sources of
-1.06 -1.06
uncertainty—stochastic errors and extrapolation to
TL. For n ≠ 1, uncertainty from the constrained-path
-1.062 -1.062
approximation is not estimated by the error bar.
32
34 98 (ii) FN: Error bars account for stochastic Monte Carlo
-1.064 -1.064
20 64 errors and for extrapolation to the TL. Uncertainties
72 due to the fixed-node approximation are not in-
E/t
E/t
-1.066 -1.066
cluded in the error bar.
U/t=4, T/t=0.25, n=0.8 (iii) MRPHF: Results are not extrapolated to the TL, and
-1.068 DCA -1.068
DCA - TL on each finite system, an estimate of the uncertainty
DiagMC - G Γ
2 due to truncation in the number of Slater determi-
-1.07 -1.07
(0)
DiagMC - [G ] Γ
(0) 2
nants is not included.
0.04 0.02 0
(iv) UCCSD: Error bars do not include an estimate of
1/N c uncertainty for truncation of excitation order to
doubles.
FIG. 2. Energies obtained at U=t ¼ 4, T=t ¼ 0.25 and n ¼ 0.8 (v) DMRG: Error bars include all sources of
from G2 Γ (magenta) and ½Gð0Þ 2 Γð0Þ (turquoise) as a function of uncertainty—the extrapolation in the number of
the inverse square of the diagram order parameter α (upper axis basis states and extrapolation to TL.
label) along with DCA results obtained from finite clusters
(vi) DMET: Error bars include all sources of
plotted against the inverse of the cluster size N c (lower axis
labels).
uncertainty—uncertainty due to extrapolation in
the number of basis states of the impurity solver
and extrapolation to TL, as well as estimates of
In Fig. 2, we show data for two DiagMC methods at DMET self-consistency convergence.
U=t ¼ 4, T=t ¼ 0.25, and n ¼ 0.8, along with DCA data. (vii) DCA: Error bars include all sources of uncertainty—
DiagMC is done directly in the thermodynamic limit, and stochastic Monte Carlo uncertainties and an addi-
the results become successively more precise as more and tional estimate of uncertainty due to extrapolation to
more expansion orders are added to the series. The results the TL.
from the two diagrammatic series we show agree, within (viii) DF: Values are presented without error bars; the
error bars, with the G2 Γ series converging more smoothly effect of neglecting nonladder and higher-order
than ½Gð0Þ 2 Γð0Þ . The convergence with expansion order in diagrams is not quantified.
the regimes we present is very rapid, so the value at order (ix) DiagMC: Error bars include all sources of
α ¼ 6 or 7 can be taken as representative of the infinite uncertainty—the stochastic Monte Carlo uncertainty
order series, with error bars estimated by statistics and by at each expansion order and an estimate of the
comparison to the results at the second largest order; in uncertainty in convergence of expansion order.
other words, extrapolation to α → ∞ is not needed.
DCA for the 2D Hubbard model approaches the thermo- V. RESULTS AT INTERMEDIATE TO STRONG
dynamic limit ∼1=N c . However, in the parameter regime INTERACTION STRENGTHS
considered here, the many-body physics is converged with We begin our discussion of results with an analysis of
respect to N c and deviations from the thermodynamic limit an intermediate to strongly coupled case, namely, U=t ¼ 8.
are dominated by single-particle shell effects. In other Throughout all the figures, we use common legends,
regimes, especially at larger U, extrapolation in 1=N c is distinguishing techniques by symbol and color. We present
required (see, for example, Refs. [113,114,117]). both results for the thermodynamic-limit and the finite-
The key result of this section is that, in many cases, system-sizedata fromwhichthethermodynamic-limit results
extrapolation to the infinite-system-size limit is needed to were obtained. This information is useful in assessing both
obtain accurate results, with the value obtained by extrapo- the importance of the extrapolation and other aspects of the
lation significantly different from the value obtained by the performance of the method. We use open symbols to denote
largest size studied. For this reason, we typically display values in the thermodynamic limit and filled symbols for
below both the extrapolated thermodynamic-limit results finite-size values from which the extrapolations are obtained.
and the finite-size results that produced the extrapolation.
A. Half-filled, particle-hole symmetric case
B. Sources of uncertainty (U=t ¼ 8, n ¼ 1, t0 =t ¼ 0)
For clarity, we repeat the main sources of uncertainties We begin our discussion with an analysis of the energy
and the meaning of the error bars shown in the graphs for per site for a widely accessible parameter choice, half
each technique; further details can be found in the sections filling, showing in Fig. 3 the temperature dependence of the
on each method. energy and in Fig. 4 an expanded view of the T ¼ 0 energy.
041041-13
order of 10−3 ) can be systematically reduced with addi-

tional computation. Thermodynamic-limit data are
20 obtained from the 1=N c extrapolation. Computational
-0.35 32 34
50 scaling towards low temperatures results in an increase
72
of uncertainty for fixed computational time, and this is
reflected in the uncertainty in the extrapolated values. At
-0.4
T=t ¼ 0.5, our results agree within error bars with high-
20
temperature series and lattice Monte Carlo data (see
E/t
32 34 Ref. [114]).
-0.45 50 72
TL The results of a DF calculation are shown at T=t ¼ 0.5
(lower T data are not available). The DF technique neglects
vertex functions of higher order than two-particle vertices.
-0.5
Furthermore, at the two-particle level, we sum only a ladder
n=1, U=8 series in the spin and charge channels. Despite these
0 0.1 0.2 0.3 0.4 0.5 approximations, we see that the DF technique provides
T/t an energy which falls on top of that of the extrapolated
thermodynamic-limit DCA result.
FIG. 3. Temperature dependence of the energy for n ¼ 1 for Results from a variety of algorithms are available at zero
U=t ¼ 8 obtained by DCA (black circles) and DF (red cross) and temperature. Figure 4 presents an expanded view of the T ¼
compared to zero-temperature results compiled from various
0 results, with the energy on the vertical axis and data for
techniques. Solid symbols represent finite systems; open symbols
represent extrapolations to the TL. each method offset in the x axis. Note that in some cases, the
thermodynamic-limit results are further offset for clarity.
We start our discussion of zero-T results with a
Monte Carlo technique, AFQMC, which is extrapolated
to the thermodynamic limit. In this case, finite-size results
are averaged over twisted-boundary conditions, which
-0.51 allows a smooth and rapid convergence to the thermody-
namic limit. These results, obtained at half filling from
2x2
Monte Carlo, are unbiased and therefore expected to be
-0.515 4x2 exact within a quoted uncertainty of 0.0002t.
8x8
8x2 DMRG results on cylinders of infinite length but finite
6x6
4x4 8x8 width of 3, 4, 5, and 6 for 45° rotated systems and width of
E/t
-0.52 10x10 4 and 6 for nonrotated systems are shown. All the finite-size
6x6 SDT(Q) data are after phase averaging, showing only very weak
_ _ 98, 162, 242
5√2×∞ 6√2×∞ finite-size effects so that an extrapolation to the thermo-
_
-0.525 6x6 4√2×∞
6x∞ dynamic limit is feasible. In this case, the estimation of
_ UCCSDTQ*
3√2×∞
n=1, U=8 uncertainty (as discussed in Sec. III F) contains the uncer-
4x4 4x∞ T=0 tainty of each extrapolation and the difference between the
-0.53 two orientations (rotated and nonrotated), both of which are
on the order of 10−4 t. The resulting energy is close to, but
FIG. 4. Thermodynamic-limit ground-state energy for n ¼ 1 slightly outside of, the AFQMC results. This extrapolation
for U=t ¼ 8 as obtained by various algorithms (open symbols). issue was discussed in more detail in Sec. IVA.
Also shown are the finite-size systems (filled symbols with
For DMET, we show results obtained for finite clusters of
adjacent labels) from which the TL ground-state energy was
obtained. Data are from AFQMC (red crosses), DMET (blue size 2 × 2, 4 × 2, 8 × 2, and 4 × 4. The thermodynamic
triangles), UCCSD (maroon diamonds), MRPHF (purple trian- limit is obtained by extrapolating the 2 × 2, 4 × 2, and 4 × 4
pffiffiffiffiffiffi
gles), DMRG (orange squares), and FN (green triangles). clusters in 1= N c . Errors from the solution of the finite
Horizontal thin dotted lines show the best estimates for the impurity are on the order of 10−4 t. DMET cluster-size
ground-state energy. convergence is slower at large U; thus, U=t ¼ 8 corresponds
to the largest half-filling DMET error bar discussed here.
The total thermodynamic-limit uncertainty is estimated to be
First, we discuss the DCA results in Fig. 3; in this 0.001t and comes entirely from the thermodynamic-limit
particle-hole symmetric parameter regime, the impurity extrapolation. The lower end of the DMET error bar lies at
solvers have no sign problem, and have cubic complexity, the average of the DMRG and AFQMC estimates.
meaning that reliable results can be obtained on relatively For the FN technique, a diffusion Monte Carlo calcu-
large clusters and that the Monte Carlo errors (here, on the lation based on the nodal structure of a trial wave function
041041-14
obtained with variational Monte Carlo, we show finite-size

results for a sequence of 45-degree rotated clusters with
size 98, 162, and 242, which have the property of being 0.055
closed shells at U ¼ 0. The results show only weak size
dependence, so the thermodynamic-limit value shown is
close to the finite-size results. However, the results are
systematically above the values obtained by AFQMC and
0.0545
DMRG, while they are consistent with DMET. This is a 8x8
consequence of the fixed-node approximation, which, in
D
0.05
this particular case, resulted in a fixed-node error of 5xinf
6x6
about 0.0015t. 6xinf 0.054
The results of UCCSD are shown for systems of size 4x2
4x4
6 × 6, 8 × 8, and 10 × 10, and exhibit weak finite-size 2x2
effects at this U value. We see that the result is accurate xxxx

4x4
2x2 n=1, U=8
to roughly the 1% level. The deviation is caused by 0.045
0.0535
correlations that are not captured by singles and doubles. 0 0.1 0.2 0.3 0.4 0.5
T/t
Higher-order excitations (triples, quadruples, etc.) will even-
tually recover the remaining energy. To support this claim, we FIG. 5. Double occupancy data for U=t ¼ 8 and n ¼ 1. Main
show a single case in Fig. 4 labeled as UCCSDT(Q), which panel: Temperature dependence of double occupancy, obtained
includes all triples and a subset of quadruples. In this higher- from DCA (finite T, black circles) and DF (finite T, red plus
order approximation, the deviation from other techniques is sign), and the T ¼ 0 techniques AFQMC (red crosses), DMET
reduced by a factor of 2. The higher-order corrections are (blue triangles), UCCSD (maroon diamonds), MRPHF (purple
triangles), DMRG (orange squares), and FN (green triangles).
more important for these coupled cluster calculations than
Solid symbols represent finite systems; open symbols represent
extrapolations in cluster geometry size; however, the extrapolations to the thermodynamic limit. Inset: Data at T ¼ 0
improvement with increased excitation order converges reproduced with an arbitrary x-axis offset, from MRPHF,
slowly. Also shown are approximate results including qua- UCCSD, FN, DMET, DMRG, and AFQMC.
druples from small system sizes (4 × 4), which we label
UCCSDTQ*. While not exact, this approximation scheme DMRG results for the double occupancy. DMET obtains a
produces results that deviate from AFQMC by only 0.15%. value (after thermodynamic-limit extrapolation) compa-
MRPHF calculations have been performed for several rable to AFQMC and, overall, shows a weaker system-
finite systems (4 × 4, 6 × 6, and 8 × 8). As summarized in size dependence than for the energy.
the methods description, reaching a constant level of accuracy UCCSD and FN produce a double occupancy that is
would require a successively larger MR expansion. Results underestimated as compared to AFQMC and DMET.
for larger systems are therefore solved less precisely; in Finite-size effects of FN are on the order of 0.0001.
particular, the energy of the 8 × 8 lattice in Fig. 4 is too high. Finally, for MRPHF, we quote two values for 4 × 4 and
More sophisticated implementations and additional optimi- 6 × 6 systems which show a system-size dependence on the
zations may make it possible to reach the accuracy needed to order of 0.001. This makes a thermodynamic-limit extrapo-
perform extrapolations to the thermodynamic limit. lation impractical.
We now discuss the results for the double occupancy in
Fig. 5 at U=t ¼ 8 and n ¼ 1. Open symbols denote results B. Doped strongly correlated regime
in the thermodynamic limit, filled symbols results on finite (U=t ¼ 8, n ¼ 0.875, t0 =t ¼ 0)
systems. The finite-T DCA results show that the double The half-filled particle-hole symmetric case of Sec. VA
occupancy contribution rises as the temperature is lowered. is, in many ways, ideally suited for numerical algorithms: A
The finite-T results are consistent with the T ¼ 0 values large charge gap allows methods like the DMRG to quickly
obtained by AFQMC, DMRG, DMET, FN, and MRPHF. converge, and particle-hole symmetry makes Monte Carlo
At T=t ¼ 0.5, the double occupancy obtained from the DF simulations without a sign problem possible. We now turn
technique is also shown. Unlike the total energy, the DF to a case that is particularly difficult to simulate, where we
double occupancy shows deviations from the DCA result, expect results to be substantially less accurate than for the
suggesting a cancellation of errors in the kinetic and half-filled case. This parameter regime shows metallic
potential energy terms. As for all other points we have behavior, strong particle-hole asymmetry, and interesting
examined, the DF method produces results that lie between inhomogeneous phases in the ground state. In Fig. 6, we
single-site DMFT values (not shown) and the extrapolated plot the total energy per site and in Fig. 7 the double
DCA results. occupancy per site at U=t ¼ 8 and n ¼ 0.875.
The inset shows the various T ¼ 0 values. Within error The main panel presents data as a function of temper-
bars, there is agreement between AFQMC, DMET, and ature. The DCA results remain consistent with T ¼ 0, but
041041-15
-0.74 0.05
TL
2
TL 20 4x4
-0.64 16
2
2 -0.75
12 4x2
8x2
2x2 0.04
-0.76
-0.68 6
T=0
E/t
D
-0.77 0.045
2x2
-0.72 0.03
8x2
6 4x2 0.0425
2 2
4 12 - 20
>
4x4
n=0.875, U=8 0.04

-0.76
n=0.875, U=8
0 0.1 0.2 0.3 0.4 0.5 0.02
T/t 0 0.1 0.2 0.3 0.4 0.5
T/t
FIG. 6. Data for n ¼ 0.875 for U=t ¼ 8. Main panel: Temper-
FIG. 7. Data for n ¼ 0.875 for U=t ¼ 8. Main panel: Temper-
ature dependence of E=t compiled from various techniques. Solid
ature dependence of double occupancy D, compiled from various
symbols represent finite systems; open symbols represent extrap-
techniques. Solid symbols represent finite systems; open symbols
olations to the thermodynamic limit. Finite-T results are shown
represent extrapolations to the thermodynamic limit. Finite-T
for DCA (black circles), and zero-T data from AFQMC (red
results are shown for DCA (black circles), and zero-T data from
crosses), DMET (blue triangles), UCCSD (maroon diamonds),
DMET (blue triangles), MRPHF (purple triangles), UCCSD
DMRG (orange squares), and FN (green triangles). Top left inset:
(maroon diamonds), DMRG (orange squares), and FN (green
Zoom-in of the zero-T data from DMET, AFQMC, FN, and
triangles). Inset: Zoom-in of the zero-T data from FN, AFQMC,
DMRG including finite-system-size data (as labeled).
DMET, and DMRG.
results are not available at the lowest benchmark temper-

ature, T=t ¼ 0.125, because of a large sign problem in the which does not allow for inhomogeneous order, yields a
Monte Carlo impurity solver. The finite-size effects, at value of E=t ¼ −0.737ð5Þ. Since the energy changes
the system sizes accessible in DCA, are smaller than in the nonmonotonically—the 8 × 2 energy lies above the
n ¼ 1 case. 2 × 2 energy but below the 4 × 2 energy—the uncertainty
The inset to Fig. 6 presents data at T ¼ 0 with an in the thermodynamic-limit extrapolation is very large and
arbitrary x-axis offset added for clarity. The AFQMC does not provide any more information than the results
simulation, using a (nonvariational) constrained-path obtained from the largest clusters.
approximation in the absence of particle-hole symmetry, The FN method shows a clear finite-system-size depend-
yields a result for the total energy that is lower than the one ence. The infinite system value is estimated from the 16 ×
obtained from DMET, FN, and DMRG. The total energy 16 and 20 × 20 values, and finite-size errors are on the
difference is about 1% when compared to finite-size order of 0.001t, much larger than the stochastic errors of
DMRG, and about 1.4% (2.1%) when compared to the 0.00001t. Here, the FN results are consistent with DMET
DMET 8 × 2 (thermodynamic-limit) cluster, respectively, extrapolation, which omits the 4 × 4 cluster, and these are
and about 2.1% in comparison to FN. considerably higher than AFQMC. This is suggestive of a
DMRG shows the results for cylinders of infinite length fixed-node error of ≈0.015t, indicating that a uniform
and finite widths of 4 and 6 lattice sites after using phase variational wave function may not fully account for the
averaging. The energy is higher for the wider cylinder, and nature of the ground state. Indeed, the variational
for the width-6 cylinder, the energy is above the energy Monte Carlo (VMC) error is of the order of 0.022t, much
from AFQMC. Given that the extrapolation is performed larger than the one obtained at half filling, which was
with only two widths, we consider the extrapolated DMRG ≈0.004t. In both cases, the FN projection improves the
value to be unreliable in this case and omit it entirely. VMC results by the same order of magnitude.
DMET shows a large-system-size dependence and a UCCSD and MRPHF results are much higher in energy
dependence of the thermodynamic-limit value on the [E=t ¼ −0.7094ð5Þ for MRPHF (16 × 4 system) and
cluster sequence chosen for the extrapolation. We show E=t ¼ −0.7122 for UCCSD, barely visible on the main
an extrapolation based on 2 × 2, 4 × 2, and 8 × 2 clusters. panel], an indication that correlated metallic states are
The use of the 8 × 2 cluster allows inhomogeneous order to difficult to capture with these methods. Although not
develop, giving an extrapolated value of E=t ¼ −0.749ð7Þ. shown, data are available for UCCSD(T) (perturbative
The extrapolation using the 4 × 4 rather than 8 × 2 cluster, inclusion of triples) in Ref. [46], which improves upon
041041-16
the value from UCCSD and gives E=t ¼ −0.7272

(−0.7281) for a 16 × 4 (16 × 8) cluster. Full inclusion of
triples (UCCSDT) lowers the 16 × 4 estimate to -0.3
2x2
E=t ¼ −0.7427. The MRPHF results indicate the need T=0
8x8
4x2
for a much larger MR expansion than that afforded in 2
2
6 8 102
-0.35 4x4 -0.52
this work.
162 242
In this parameter regime, ordered “stripe” phases might 98 6xinf
6x6
exist. However, the precise form of these stripes is strongly -0.4 4xinf 4x4
E/t
-0.53
influenced by choice of finite-size systems (e.g., width and
orientation of the cylinder in DMRG and shape of the cluster -0.45 T=0
in DMET) that are used for the thermodynamic extrapolation -0.524
and the approximations used to solve that finite system. The
finite-temperature algorithms have not reached the onset of -0.5
-0.526
inhomogeneous states at the lowest temperature accessible. n=1, U=8, t’=-0.2
The precise nature of the inhomogeneities in the ground state -0.55
0 0.1 0.2 0.3 0.4 0.5
in this parameter regime is still open. T/t
Finally, we briefly mention the results for double occu-
pancy in Fig. 7 for U=t ¼ 8 and n ¼ 0.875. As was the case FIG. 8. Data for n ¼ 1 for U=t ¼ 8 with t0 =t ¼ −0.2. Main
for the energies, the finite-T results smoothly connect to the panel: Temperature dependence of E=t compiled from various
zero-T values. MRPHF and UCCSD overestimate the double techniques. Solid symbols represent finite systems; open symbols
occupancy by close to 15%. The remaining ground-state represent extrapolations to the thermodynamic limit. Finite-T
methods (DMRG, AFQMC, DMET, and FN) present con- results are shown for DCA (black circles), and zero-T data from
sistent values in the range from 0.04 to 0.043. Both FN and AFQMC (red crosses), DMET (blue triangles), UCCSD (maroon
diamonds), MRPHF (purple triangles), DMRG (orange squares),
AFQMC values contain additional (fixed-node and con-
and FN (green triangles). Top left inset: Zoom-in of the zero-T
strained-path) errors that are not estimated by the error bar. data from MRPHF, UCCSD, DMET, FN, DMRG, and AFQMC,
including finite-system-size data (as labeled) for MRPHF, FN,
C. Half-filled, non-particle-hole symmetric case DMET, and DMRG. Bottom right panel: Enlarged region of
(U=t ¼ 8, n ¼ 1, t0 =t ¼ −0.2) the top left inset showing DMET, DMRG, FN, and AFQMC
data at T ¼ 0, including error bars, extrapolated to the infinite
We now turn our attention to a case of half filling without system size.
particle-hole symmetry, by adding a second nearest-
neighbor hopping t0 . An overview of the energies from
several algorithms for U=t ¼ 8, n ¼ 1, and t0 ¼ −0.2 is finite-system-size dependence of the fixed-node results is
shown in Fig. 8. small on this scale.
The main panel shows the temperature dependence of the UCCSD results show only small finite-size effects and
data. The DCA results available at finite T show almost no an overall energy ≈1% higher than other techniques. The
sign problem for T=t ¼ 0.5 and T=t ¼ 0.25 but are MRPHF results obtained on finite systems show an energy
hampered by a severe sign problem at T=t ¼ 0.125. The that rises rapidly as the system size is increased. As in the
results are consistent within error bars with the zero- case of t0 ¼ 0, a systematic extrapolation to the thermo-
temperature results. dynamic limit is not possible, and we only present results
As at U=t ¼ 8, n ¼ 0.875, t0 ¼ 0, the AFQMC is on finite systems.
approximate because of a constrained-path approximation
due to the lack of particle-hole symmetry. Despite this, the VI. RESULTS IN THE INTERMEDIATE
results are in agreement with both the DMET and DMRG COUPLING REGIME
results.
In this section, we repeat the previous discussion for an
The DMET results are obtained on clusters of size 2 × 2,
interaction strength of half the size, U=t ¼ 4. As before, we
2 × 4, and 4 × 4. Errors of the individual finite-size systems
start our discussion at half filling. We then discuss a
are substantially smaller than the system-size dependence.
correlated metallic case with 20% doping.
The DMET thermodynamic limit is consistent with the
thermodynamic estimates obtained from DMRG (from
cylinders of width 4 and 6) and from AFQMC. This is A. Half-filled, particle-hole symmetric case
even more evident in the bottom right inset, which displays (U=t ¼ 4, n ¼ 1, t0 =t ¼ 0)
the thermodynamic-limit estimates on a smaller scale. In Fig. 9, we report the energy as a function of temper-
FN results are higher in energy than AFQMC and DMET ature. At finite T and U=t ¼ 4, both DCA and the
(well within 2 joint standard deviations) and are higher diagrammatic Monte Carlo method for the ½Gð0Þ 2 U series
than DMRG by 0.0013t. As seen in previous plots, the provide results in the thermodynamic limit. DCA results in
041041-17
the thermodynamic limit are extrapolated from finite

clusters; DiagMC results are extrapolated in the expansion
order. The results are consistent within the error bars of the
respective methods. The large error bars of the extrapola-
tion in DiagMC-½Gð0Þ 2 U mainly come from a conservative 0.14
20
32
estimate of the diagram-order extrapolation error. DCA 34
shows surprisingly large finite-size effects which persist

above N c ¼ 72, unlike at U=t ¼ 8. While each individual 50
72
N c result has uncertainties in the energy on the order of 0.13 98
10−4 t, the spread in values results in a large uncertainty
D
when extrapolated to the thermodynamic limit. 2x2 T=0
4x2
We now move to the zero-temperature methods, which 0.12
0.127
4x4
12x12 162
are shown in the inset of Fig. 9. AFQMC provides 0.126 6

10x10
8x8
242
numerically exact ground-state energies for this system. 4 8x8
The value quoted is E=t ¼ −0.8603ð2Þ, which is in agree- 0.125

n=1, U=4
ment with the DMET value of E=t ¼ −0.8604ð3Þ and the 0.11
0 0.1 0.2 0.3 0.4 0.5
DMRG value of E=t ¼ −0.8605ð5Þ. DMET values are T/t
obtained from an extrapolation of 2 × 2, 2 × 4, and 4 × 4
clusters. DMRG values are obtained from an extrapolation FIG. 10. Data for n ¼ 1 for U=t ¼ 4. Main panel: Temperature
of widths 3, 4, and 5 for 45-degree rotated cylinders and of dependence of double occupancy D, compiled from various
widths 4 and 6 for nonrotated cylinders. techniques. Solid symbols represent finite systems; open symbols
represent extrapolations to the thermodynamic limit. Finite-T
The results obtained by AFQMC, DMRG, and DMET
results are shown for DCA and DiagMC (turquoise stars), and
are in excellent agreement with recent calculations zero-T data from DMET (blue triangles), UCCSD (maroon
obtained from linearized auxiliary-field Monte Carlo diamonds), MRPHF (purple triangles), DMRG (orange squares),
(LAQMC) available in the literature [141], which gives and FN (green triangles). Inset: Zoom-in of the zero-T data from
E=t ¼ −0.85996ð5Þ. FN results are higher in energy DMRG, FN, UCCSD (only finite-system data), DMET, and
AFQMC.
[E=t ¼ −0.8575ð3Þ], and unlike in previous cases for

stronger interaction, a clear dependence on the finite
4x4 system studied is visible. The FN projection technique
-0.7
6x6 leads to an energy gain of ≈0.002t with respect to the
8x8
8x8
12x12 -0.855 VMC result of E=t ¼ −0.8558ð5Þ. This number can be
6x6 SDTQ* compared with a previous estimation of the thermody-
-0.75 242
162 -0.86 namic limit in VMC obtained with a slightly less accurate
4x4
4x2 6 98 variational state (see Ref. [78]).
2x2
E/t
The UCCSD thermodynamic limit overestimates the

-0.8 energy by about 0.7%. MRPHF values, shown as purple
triangles in the main panel and inset, show large finite-size
effects and are higher than the values obtained with other
-0.85 methods. We see that as the system size is increased, the
n=1.0, U=4 energy increases rapidly.
In Fig. 10, we report the double occupancy vs T. At finite
0 0.1 0.2 0.3 0.4 0.5
T/t T, the DCA results show a clear rise in D as T decreases
from 0.5t → 0.25t. However, this trend reverses as T
FIG. 9. Data for n ¼ 1 for U=t ¼ 4. Main panel: Temperature decreases further. A similar behavior is also obtained by
dependence of E=t compiled from various techniques. Solid using DiagMC, demonstrating that this is a genuine effect
symbols represent finite systems; open symbols represent extrap- present in the Hubbard model.
olations to the thermodynamic limit. Finite-T results are shown These trends are consistent with the T ¼ 0 data, which
for DCA (black circles) and DiagMC (turquoise stars), and zero-
lie below all of the DCA data points at finite T. For clarity
T data from AFQMC (red crosses), DMET (blue triangles),
UCCSD (maroon diamonds), MRPHF (purple triangles), DMRG of presentation, we again display this data in the inset and
(orange squares), and FN (green triangles). Top left inset: Zoom- add an arbitrary x-axis offset. We see that finite-size
in of the zero-T data from MRPHF, DMET, FN, DMRG, and effects in DMET are very small and that the extrapolation
AFQMC, including finite-system-size data (as labeled) for of DMET agrees perfectly with AFQMC. Finite-size FN
MRPHF, FN, DMRG, and DMET. results produce values comparable to DMRG, DMET,
041041-18
-1.07
T=0.25 -1.06 -1.08 10x8

α =3
50 10x12
-1 10x16 10x4 SD
64 α =4
72 -1.065
α >4 -1.09 10x4
SD(T)
E/t
-1.07 SDT(Q) SDT
-1.04 iω n
-1.1 SDTQ* 400
E/t
0 4 8
0
6×∞
4×∞
T=0.25 3×∞ 100
-0.1 -1.11
200
Im Σ (iωn)
-1.08 -0.2 4x4
α =5 8x2
α =3 4x2
-0.3 -1.12 n=0.8, U=4 T=0
α =4 2x2
-0.4
n=0.8, U=4
α =7
-1.12 -0.5
0 0.1 0.2 0.3 0.4 0.5 FIG. 12. Thermodynamic-limit ground-state energy for n ¼ 0.8
T/t for U=t ¼ 4 as obtained by various algorithms (open symbols).
Also shown are the finite-size systems (filled symbols with
FIG. 11. Data for n ¼ 0.8 for U=t ¼ 4. Main panel: Temper- adjacent labels) from which the thermodynamic-limit ground-
ature dependence of E=t compiled from various techniques. Solid state energy was obtained. Data are from AFQMC (red crosses),
symbols represent finite systems; open symbols represent extrap- DMET (blue triangles), UCCSD (maroon diamonds), UCC on a
olations to the thermodynamic limit. Finite-T results are shown 10 × 4 lattice (shaded diamonds), MRPHF (purple triangles),
for DCA (black circles) and DiagMC (pink and turquoise DMRG (orange squares), and FN (green triangles). Horizontal
asterisks), and zero-T data from AFQMC (red crosses), DMRG thin dotted lines show the best estimates for the ground-state
(orange squares), FN (green triangles), and DMET (blue trian- energy.
gles). Top left inset: Zoom-in of the T=t ¼ 0.25 data from DCA
and two types of DiagMC for different orders α ¼ 3; 4; ….
Bottom right inset: Plot of the imaginary part of the local self- resummations, α. The values obtained from the three
energy ImΣðiωn Þ from DiagMC for different expansion orders α techniques are within error bars.
and from DCA (black circles covered by magenta stars).
The T > 0 values smoothly connect to T ¼ 0 (although
a precise comparison is not possible because we lack a
and AFQMC. However, extrapolation in FN results in an quantitative functional form to extrapolate the T > 0
underestimate of the double occupancy, although within values to T ¼ 0), which we display separately in
uncertainties. In the DMRG simulations, phase averaging Fig. 12, where data from DMET, AFQMC (constrained
has greatly reduced finite-size effects, and the DMRG path), UCCSD, MRPHF, FN, and DMRG are shown.
error bars are determined by the truncation errors. Within In this case, MRPHF and UCCSD are systematically
those error bars, DMRG results are consistent with higher in energy from the other techniques. In the case of
AFQMC and DMET. UCCSD, we see larger finite-size effects than at U=t ¼ 8.
Inclusion of higher orders of excitation [perturbative
triples (T), triples T, and perturbative quadruples (Q)]
B. Doped case (U=t ¼ 4, n ¼ 0.8, t0 =t ¼ 0) suggests that the dominant error is associated with the
Away from half filling (with t0 ¼ 0), we can perform truncation of the excitation order and not with finite size
further comparisons at finite T between DCA and effects. The FN result is in agreement with the value from
DiagMC at U=t ¼ 4. We begin the discussion with the DMRG. At slightly lower energy, AFQMC (constrained
inset of Fig. 11, which shows the convergence of path) and DMET are in close agreement. Overall, the
the imaginary part of the local Matsubara self-energy spread in energies is similar to that shown away from half
of the G2 Γ and G2 W DiagMC series as a function of filling at U=t ¼ 8 (see Fig. 6) but smaller in magnitude,
evaluation order. The values are compared to DCA results. perhaps due to better convergence for weak coupling in
We see that the first six orders of the series are precise some techniques.
enough to get good agreement of the Matsubara self-
energy in the thermodynamic limit, and convergence is
VII. RESULTS IN THE WEAK COUPLING
rapid. While deviations are visible in the energy, these are
REGIME
attributed to differences of the chemical potential, i.e., the
real part of the self-energy. In the weak coupling limit, we restrict the presentation of
The top left inset of Fig. 11 shows the convergence of the data to the half-filled case since the correlated metallic phase
energy in DCA and two different series of DiagMC, is not qualitatively distinct from U=t ¼ 4. Data sets for
½Gð0Þ 2 U and G2 Γ, for increasing order of the diagrammatic doped, weakly correlated systems are available in Ref. [46].
041041-19
A. Half-filled, particle-hole symmetric case agree precisely, while FN is slightly higher in energy but
(U=t ¼ 2, n ¼ 1, t0 =t ¼ 0) compatible within two error bars. In the upper left inset,
In Figs. 13 and 14, we present results for U=t ¼ 2 and we explore the finite-size effects of the methods. In the case
n ¼ 1, the half-filled weak coupling regime. This regime is of DMET, these finite-size effects are small and can be
particularly easy for methods based on a weak coupling extrapolated with small error bars. FN shows much larger
expansion around a noninteracting system, and many of finite-size effects, approaching the thermodynamic-limit
the algorithms show uncertainties that are much smaller energy from below (only the largest system size is visible
than in the intermediate or strong interaction limit. on the scale of the plot). DMRG results with phase
At nonzero T, the data from two types of DiagMC and averaging are precise even at U=t ¼ 2, though much larger
from DCA in the thermodynamic limit agree within uncertainties than at U=t ¼ 8 are present.
uncertainty. The values smoothly connect to the T ¼ 0 The results of MRPHF, outside of the scale shown by the
values, except for MRPHF energies, that are higher than the inset, show a gradual decrease of the energy with increasing
ground-state energies obtained by the other methods and system size: 4 × 4 yields E=t ¼ −1.1260, 6 × 6 yields
higher than the energies obtained for the lowest-temper- E=t ¼ −1.1515, and 8 × 8 yields E=t ¼ −1.1629.
ature point obtained from both finite-T methods. UCCSD results, shown here for cluster sizes of 12 × 12,
At T=t ¼ 0.5, we show a result from DF. As was the case 10 × 10, 8 × 8, and 6 × 6, show large finite-size effects. With
in the strong coupling regime, the DF procedure produces the aid of an extrapolation, the value in the thermodynamic
an energy estimate consistent with DCA results. In this limit is estimated to be higher than other techniques. The
case, with only weak finite-size dependence, the underlying deviation of the thermodynamic-limit value from AFQMC
DMFT approximation differs from DCA by only 0.4%. and DMET is on the order of 2 × 10−3 t, suggesting that
The DF value improves on the DMFT and differs from the excitations beyond the single and double levels are important
extrapolated DCA results by only 0.07%. even in this relatively weak coupling regime.
In the lower right inset of Fig. 13, we present T ¼ 0 Finally, Fig. 14 shows the double occupancy for these
extrapolations. As in the case of larger interactions, DMET parameters. We see an increase in finite-size effects in DCA
and AFQMC (which is numerically exact in this situation) as we progress to lower temperatures. Reasonable agree-
ment with DiagMC is achieved in the double occupancy.
At temperatures lower than our lowest temperature, the
double occupancy will need to dip, as was the case at
T=0
6xinf -1.175
SDTQ*
-1.05 4x4
4x2
4xinf
2x2
-1.18
n=1, U=2
242
0.2
E/t
-1.1
-1.174 T=0
8x8
98
162
-1.175 0.2
-1.15 6x6 242
0.18
D
4x4
-1.176 0.195 2x2
n=1, U=2 4x2
4x4
0.19
0 0.1 0.2 0.3 0.4 0.5 12x12
T/t 0.16 _
10x10
0.185 4√2_×∞
3√2×∞
FIG. 13. Data for n ¼ 1 for U=t ¼ 2. Main panel: Temperature 8x8 8x8 T=0
dependence of E=t compiled from various techniques. Solid 0 0.1 0.2 0.3 0.4 0.5
symbols represent finite systems; open symbols represent extrap- T/t
olations to the thermodynamic limit. Finite-T results are shown
for DCA (black circles) and DiagMC (blue asterisks), and zero-T FIG. 14. Data for n ¼ 1 for U=t ¼ 2. Main panel: Temperature
data from AFQMC (red crosses), DMET (blue triangles), dependence of double occupancy D, compiled from various
UCCSD (maroon diamonds), MRPHF (purple triangles), DMRG techniques. Solid symbols represent finite systems; open symbols
(orange squares), and FN (green triangles). Top left inset: Zoom- represent extrapolations to the thermodynamic limit. Finite-T
in of the zero-T data from UCCSD, DMET, FN, DMRG, results are shown for DCA (black circles) and DiagMC (turquoise
AFQMC, and DMET. Bottom right panel: Enlarged region of asterisks), and zero-T data from DMET (blue triangles), UCCSD
the top left inset showing thermodynamic-limit data for DMRG, (maroon diamonds), MRPHF (purple triangles), DMRG (orange
DMET, UCCSD, FN, and AFQMC at T ¼ 0, including error squares), and FN (green triangles). Inset: Zoom-in of the zero-T
bars, extrapolated to the infinite system size. data from FN, UCCSD, MRPHF, DRMG, DMET, and AFQMC.
041041-20
U=t ¼ 4, in order for the finite-T data to be consistent with

T ¼ 0. Similar to the strong coupling case, DF provides
only a slight shift to the double occupancy, a minimal
improvement over DMFT alone.
The ground-state double occupancies are very precisely
determined by AFQMC and DMET, which are in agree-
ment. The FN value is somewhat overestimated. Results
from DMRG fall below AFQMC and DMET, and the error
bar underestimates the uncertainty. The larger error appears
to be consistent with the difficulty in treating the small U
limit in the DMRG calculations. The results from MRPHF
show an improvement in D as the system size is increased,
consistent with the behavior for the energy. In the case of
UCCSD, since it is an expansion in the coupling strength, at
weak coupling the procedure is more reliable, and while FIG. 16. Comparison of the real part of the self-energy,
there are substantial finite-size effects, the extrapolation ReΣðk; iω0 Þ, at the lowest Matsubara frequency iω0, obtained
produces a result within error bars of AFQMC and, in from DF (red) compared to 72-site DCA calculations (black)
general, agreement with DMET. plotted as a function of momentum k, throughout the Brillouin
zone for n ¼ 1.0, U=t ¼ 8, T=t ¼ 0.5. The dual fermion and
DCA self-energies are plotted as step functions. Also included are
VIII. FREQUENCY AND MOMENTUM interpolated results obtained by diagrammatic determinantal
DEPENDENCE Monte Carlo (DDMC) [120,140,142,143] (dashed black) with
a gray shading to indicate the level of uncertainty.
Next, we discuss single-particle finite-temperature prop-
erties. All finite-temperature algorithms discussed in this
In this metallic regime, self-energies are small. Black circles
work are based on approximations of the single-particle
denote the imaginary part of the local self-energy from an
self-energy. We show three characteristic plots for this
N c ¼ 20 DCA calculation, which for these parameters
quantity: Figure 15 shows the imaginary part of the local
shows essentially no finite-size effects. The data agree
self-energy as a function of Matsubara frequency, Fig. 16
perfectly with DiagMC-G2 W data shown as red dashed
shows the dependence of the real part of the lowest
lines, and convergence of the DiagMC-G2 Γ method to the
Matsubara frequency on k space, and Fig. 17 shows the
result of the other two methods (stars, magenta dotted line)
frequency dependence of the imaginary part of the self-
is observed as a function of expansion order α. This
energy for a specific momentum. Any discrepancy in the
agreement implies that the local physics is captured well
energy or double occupancy is the result of discrepancies in
by all three algorithms.
the single-particle self-energy.
The data shown in Fig. 15 are obtained for weak 0
interaction strength U=t ¼ 2 and for a density n ¼ 0.3. DCA
DF
DΓA
-1
-0.01
ImΣ(k=(π ,0) ,iωn)
U=2, T=0.5, n=0.3, t’=0

-2
-0.02
Im Σloc(iωn)
2
DiagMC G W 64x64 -3
DCA
2
α =1 DiagMC-G Γ
-4
-0.03
α =6
-5
0 5 10 15 20
iω n
-0.04
0 10 20 30 40 FIG. 17. Comparison of the frequency dependence of the
iω n imaginary part of the self-energy, ImΣðk; iωn Þ, at fixed k ¼
ðπ; 0Þ obtained from DF (red) compared to 72-site DCA
FIG. 15. Imaginary part of the local self-energy, ImΣloc ðiωn Þ at calculations (black) plotted for n ¼ 1.0, U=t ¼ 8, T=t ¼ 0.5.
U=t ¼ 2, T=t ¼ 0.5 and n ¼ 0.3 from DCA and DiagMC. In the Also shown are results from the dynamical vertex approximation
case of DiagMC-G2 Γ, we label α, the series order from Eq. (7). DΓA (blue) [124,144].
041041-21
In Fig. 16, we examine momentum-dependent data. limit, obtained from AFQMC, DMRG, DMET, and FN.
We show a path of (kx , ky ) through the Brillouin zone The MRPHF and UCCSD values presented show the value
and plot the real part of the self-energy at the lowest for the largest finite-size system studied. The uncertainties
Matsubara frequency at U=t ¼ 8, β ¼ 2, and n ¼ 1. presented are the best uncertainties available within each
DiagMC data are not available in this regime, but for algorithm and do not contain an assessment of systematic
comparison, we plot large DCA cluster results (N c ¼ 72) errors (e.g., fixed node or truncation in expansion order).
and results from continuous-time lattice Monte Carlo We see that for much of phase space, errors are a few times
simulations (see Refs. [120,140,142,143]). The DCA 10−4 t, and values between the techniques are remarkably
approximation (blue lines) produces a step- consistent.
discretized self-energy, which is in approximate agreement Table III shows the double occupancy for the same
with the momentum dependence from other techniques. values. Relative errors for the double occupancy are of
Discrepancies between the approximate DF method and the the same magnitude as for the total energy, and values in the
(essentially converged) DDMC method are visible but thermodynamic limit are consistent within error bars.
within the uncertainty of the DDMC comparison data. Tables IV and V present ground-state energies for the
Any discrepancies are expected to rapidly disappear at densities n ¼ 0.8 and 0.875 for values of U=t considered in
higher Matsubara frequencies, which can be seen in a this work. The full table of values for the data presented in
comparison of the imaginary part of the self-energy this paper is available online (see Ref. [46]).
(Fig. 17) of DF and DCA at fixed momentum k ¼ ðπ; 0Þ. Not all quantities are as consistent as the energies. This is
For comparison or verification purposes, we include results especially true for order parameters and correlation func-
from the dynamical vertex approximation (DΓA, see tions, where discrepancies outside of error bars are present.
Refs. [124,144,145]), which, in a spirit similar to DF, solves Presumably, many competing phases exist in a narrow
the model in an expansion of two-particle vertex functions. energy window near the ground state, and the most
We find that the results from DF and DΓA are consistent. favorable state found in each method will depend on
details of the finite-size system and the approximation.
Table VI shows the comparison between the magnetization
IX. TABLES FOR GROUND-STATE PROPERTIES
that DMET observes and the magnetization found in
We conclude the discussion of our results with a list of AFQMC for the full range of U=t at half filling. At weak
the thermodynamic-limit estimates for the half-filled t0 ¼ 0 interaction strength, DMET finds a larger polarization than
case. Table II shows a list of energies in the thermodynamic AFQMC even though the energies agree to all significant
TABLE II. Zero-temperature energy and uncertainty for n ¼ 1, T ¼ 0, for a range of interaction strengths U, obtained from AFQMC,
DMET, DMRG, FN, MRPHF, and UCCSD. Where extrapolations to the TL are not available, finite-size geometries are listed in lieu of
uncertainties. UCCSDTQ* data estimate higher-order corrections by including triples from a ½6 × 6 and quadruples from a ½4 × 4.
UCCSD data for U=t > 4 provide nearly converged energy estimates with respect to system size.
U 2 4 6 8 12
AFQMC −1.1763 0.0002 −0.8603 0.0002 −0.6568 0.0003 −0.5247 0.0002 −0.3693 0.0002
DMET −1.1764 0.0003 −0.8604 0.0003 −0.6562 0.0005 −0.5234 0.0010 −0.3685 0.0010
DMRG −1.176 0.001 −0.8605 0.0005 −0.6565 0.0001 −0.5241 0.0001 −0.3689 0.0001
FN −1.175 0.001 −0.8575 0.0003 −0.6551 0.0001 −0.52315 0.00005 −0.36835 0.00005
MRPHF −1.1628 ½8 × 8 −0.8554 ½8 × 8 −0.6512 ½8 × 8 −0.5169 ½8 × 8 −0.3626 ½8 × 8
UCCSD −1.1735 0.0004 −0.8546 ½14 × 14 −0.6510 ½10 × 10 −0.5191 ½10 × 10 −0.3647 ½10 × 10
UCCSDTQ* −1.1749 −0.8610 −0.6582 −0.5255 −0.3696
TABLE III. Zero-temperature double occupancy and uncertainty for n ¼ 1, for a range of interaction strengths U, obtained from
AFQMC, DMET, DMRG, FN, MRPHF, and UCCSD. Where extrapolations to the TL are not available, finite-size geometries are listed
in lieu of uncertainties.
U 2 4 6 8 12
AFQMC 0.1923 0.0003 0.1262 0.0002 0.0810 0.0001 0.0540 0.0001 0.0278 0.0001
DMET 0.1913 0.0004 0.1261 0.0001 0.08095 0.00004 0.05398 0.00007 0.02780 0.00003
DMRG 0.188 0.001 0.126 0.001 0.0809 0.0003 0.0539 0.0001 0.0278 0.0001
FN 0.198 0.001 0.125 0.001 0.0803 0.0002 0.0535 0.0001 0.0278 0.0002
MRPHF 0.1824 ½8 × 8 0.1262 ½8 × 8 0.0818 ½8 × 8 0.0544 ½8 × 8 0.0275 ½8 × 8
UCCSD 0.194 0.002 0.1268 ½12 × 12 0.0807 ½10 × 10 0.0537 ½10 × 10 0.0276 ½10 × 10
041041-22
TABLE IV. Zero-temperature energy and uncertainty for n ¼ 0.8, T ¼ 0, for a range of interaction strengths U,
obtained from AFQMC (constrained path), DMET, DMRG, FN, MRPHF, and UCCSD. Where extrapolations to the
TL are not available, finite-size geometries are listed in lieu of uncertainties.
n ¼ 0.8
U 2 4 6 8
AFQMC −1.306 0.002 −1.110 0.003
DMET −1.3062 0.0004 −1.108 0.002 −0.977 0.004 −0.88 0.03
DMRG −1.104 0.0014
FN −1.3044 0.0007 −1.1032 0.0007 −0.967 0.001 −0.877 0.001
MRPHF −1.2931 ½10 × 4 −1.0910 ½10 × 4 −0.9454 ½10 × 4 −0.8415 ½10 × 4
UCCSD −1.3065 ½10 × 16 −1.0868 ½10 × 16 −0.9300 ½10 × 16 −0.8233 ½10 × 16
UCCSDT* −1.3078 ½10 × 4 −1.0981 ½10 × 4 −0.9607 ½10 × 4 −0.8641 ½10 × 4
TABLE V. Zero-temperature energy and uncertainty for n ¼ 0.875, T ¼ 0, for a range of interaction strengths U,
obtained from AFQMC (constrained path), DMET, DMRG, FN, MRPHF, and UCCSD. Where extrapolations to the
TL are not available, finite-size geometries are listed in lieu of uncertainties.
n ¼ 0.875
U 2 4 6 8
AFQMC −1.026 0.001 −0.766 0.001
DMET −1.2721 0.0006 −1.031 0.003 −0.863 0.013 −0.749a 0.007
DMRG −1.028 ½6 × ∞ −0.759 0.004
FN −1.270 0.002 −1.0225 0.0015 −0.854 0.002 −0.749 0.002
MRPHF −1.2855 ½16 × 4 −1.0195 ½16 × 4 −0.8318 ½16 × 4 −0.7094 ½16 × 4
UCCSD −1.2667 ½16 × 12 −1.0093 ½16 × 12 −0.8298 ½16 × 12 −0.7123 ½16 × 12
UCCSDT* −1.2681 ½16 × 4 −1.0253 ½16 × 4 −0.8570 ½16 × 4 −0.7434 ½16 × 4
a
Due to strong spatial inhomogeneity observed at this filling, the TL extrapolated number in the table is not
meaningful, as different cluster shapes show different orders. Out of the clusters used here, the 8 × 2 impurity cluster
is likely the best estimate with E=t ¼ −0.755ð0.007Þ. An exhaustive study of the DMET cluster size dependence
at this filling will be carried out in the future.
TABLE VI. Comparison of magnetization data from DMET and AFQMC at n ¼ 1.
U 2 4 6 8
Technique m δm m δm m δm m δm
AFQMC 0.094 0.004 0.236 0.001 0.280 0.005 0.26 0.03
DMET 0.133 0.005 0.252 0.009 0.299 0.012 0.318 0.013
digits, as a result of DMET finite-size scaling from small various carrier concentrations, obtained using state-of-the-
clusters. At large interaction strength, AFQMC gives a art numerical methods. We believe the results are useful for
polarization with large statistical fluctuation despite very two reasons. First, the two-dimensional Hubbard model is
accurate energies. Similar behavior (not shown here) is one of the paradigm models of quantum condensed-matter
apparent for other variables, e.g., the d-wave order param- theory, and it is therefore important to determine, as reliably
eter or the stripe geometry observed at U=t ¼ 8 and as possible, the state of our knowledge about it. Second,
n ¼ 0.875. solving the grand-challenge problem of determining the
physics of interacting many-electron systems will require
numerics, and as no one technique is likely to provide
X. CONCLUSIONS
solutions in all regimes or for all quantities of physical
In this paper, we have presented a detailed examination interest, it is important to develop tools for assessing the
of results for static and dynamic properties of the two- strengths and weaknesses of different approaches.
dimensional Hubbard model at correlation strengths rang- We argue that the only quantities that can meaningfully
ing from weak to intermediate to strong coupling, and at be compared between different approaches are estimates,
041041-23
with error bars, for the thermodynamic-limit values of ACKNOWLEDGMENTS

observables, including local operators such as the energy,
We acknowledge the Simons Foundation for funding.
double occupancy, density (or chemical potential), and
The work at Rice University was supported by Grant
magnetization, as well as correlation functions such as the
No. NSF-CHE-1462434. G. K. C. acknowledges funding
electron self-energy. We restrict our attention to methods
from the U.S. Department of Energy (DOE) for the develop-
and regimes for which large enough systems can be studied
ment of the DMET method through Grant No. DE-
such that reasonable extrapolations to the thermodynamic
SC0010530, and for its application to the Hubbard model
limit can be performed. Care is required in performing the
and superconductivity through Grant No. DE-SC0008624.
extrapolations, and we have found it useful to present both
the extrapolated results and (in most cases) the finite-size F. B. and L. F. T. acknowledge support from PRIN
data that led to the extrapolation. 2010_2010LLKJBX. S. Z. and H. S. acknowledge support
Comparison of results obtained from different methods from the National Science Foundation Grant No. DMR-
shows that the ground-state properties of a substantial part 1409510 for AFQMC method development. M. Q. was also
supported by DOE Grant No. DE-SC0008627, and A. E. A.
of the Hubbard model phase space are now under numerical
by DOE Grant No. ER 46932. Y. D. and X. W. L. acknowl-
control (see, e.g., Figs. 4 and 12). Moreover, where there is
edge support from NSFC Grant No. 11275185. The
agreement on the ground-state properties, the nonzero
AFQMC calculations were carried out at the Oak Ridge
temperature methods appear to connect smoothly to the
Leadership Computing Facility at the Oak Ridge National
ground state as the temperature is decreased, although a
Laboratory, and at the computational facilities at the College
quantitative extrapolation to T ¼ 0 is not yet available. The of William and Mary.
most substantial uncertainties exist at intermediate corre-
lations (e.g., U=t ≈ 4 → 8) and at dopings near but not
equal to the half-filling value n ¼ 1. In this intermediate
coupling–near half-filling regime, several physically differ-
ent states seem to have very similar energies, and small [1] J. Hubbard, Electron Correlations in Narrow Energy
effects can favor the choice of one state over the other, Bands, Proc. R. Soc. A 276, 238 (1963).
leading to substantial uncertainties in physical quantities. [2] M. C. Gutzwiller, Effect of Correlation on the Ferromag-
Also, the issues associated with fermion sign problems netism of Transition Metals, Phys. Rev. Lett. 10, 159
(1963).
seem to be most severe. Interestingly, it is this regime that is [3] J. Kanamori, Electron Correlation and Ferromagnetism of
of most physical interest in connection with high-T c Transition Metals, Prog. Theor. Phys. 30, 275 (1963).
superconductivity in the copper-oxide materials. [4] E. H. Lieb and F. Y. Wu, Absence of Mott Transition in an
Where two or more methods produce results that agree Exact Solution of the Short-Range, One-Band Model in
within reasonable error bars, we take the result to be One Dimension, Phys. Rev. Lett. 20, 1445 (1968).
established and appropriate for use as a benchmark. Tables [5] W. Metzner and D. Vollhardt, Correlated Lattice Fermions
of our benchmark results are contained in Ref. [46] and in d ¼ ∞ Dimensions, Phys. Rev. Lett. 62, 324 (1989).
made available online. We expect that these results will be [6] E. Müller-Hartmann, Correlated Fermions on a Lattice in
useful in validating new methods or new implementations High Dimensions, Z. Phys. B 74, 507 (1989).
of existing methods. [7] J. Oitmaa, C. Hamer, and W. Zheng, Series Expansion
Methods for Strongly Interacting Lattice Models
Turning now to prospects and open questions, we first
(Cambridge University Press, Cambridge, England, 2006).
observe that all of the methods we have considered to date [8] R. Valent, J. Stolze, and P. J. Hirschfeld, Lower Bounds for
have difficulty in the physically interesting intermediate the Ground-State Energies of the Two-Dimensional Hub-
coupling, near the half-filling regime. Development of new bard and t-J Models, Phys. Rev. B 43, 13743 (1991).
methods, or improvement of existing methods to deal with [9] R. Valent, C. Gros, P. Hirschfeld, and W. Stephan,
this regime, is urgently needed. Furthermore, we remark Rigorous Bounds for Ground-State Properties of Corre-
that our understanding of dynamical correlation functions, lated Fermi Systems, Phys. Rev. B 44, 13203 (1991).
even ones as simple as the electron Green’s function, is [10] M. Y. Kagan and A. V. Chubukov, The Enhancement of a
much less advanced than our understanding of ground-state Superfluid Transition Temperature in a Polarised Fermi-
properties and simple expectation values. Finally, we Gas with Repulsion, JETP Lett. 50, 517 (1989).
observe that the process of producing this paper, in [11] M. A. Baranov and M. Y. Kagan, D-Wave Pairing in the
Two-Dimensional Hubbard Model with Low Filling, Z.
particular, the confrontation of each method with the body
Phys. B 86, 237 (1992).
of related work produced by other methods, led in several [12] A. V. Chubukov and J. P. Lu, Pairing Instabilities in the
cases to substantial improvements in algorithm and error Two-Dimensional Hubbard Model, Phys. Rev. B 46,
analysis. We suggest that as quantitative numerics in the 11163 (1992).
many-electron field continues to evolve, intercomparison of [13] A. V. Chubukov, Kohn-Luttinger Effect and the Instability
different methods, leading to benchmarking on important of a Two-Dimensional Repulsive Fermi Liquid at T ¼ 0,
model problems, will significantly advance the field. Phys. Rev. B 48, 1097 (1993).
041041-24
[14] R. Hlubina, Phase Diagram of the Weak-Coupling Two- [35] H. Yokoyama and H. Shiba, Variational Monte-Carlo
Dimensional t-t0 Hubbard Model at Low and Intermediate Studies of Hubbard Model. I, J. Phys. Soc. Jpn. 56, 1490
Electron Density, Phys. Rev. B 59, 9600 (1999). (1987).
[15] D. Zanchi and H. J. Schulz, Superconducting Instabilities [36] H. Yokoyama and H. Shiba, Variational Monte-Carlo
of the Non-half-filled Hubbard Model in Two Dimensions, Studies of Hubbard Model. II, J. Phys. Soc. Jpn. 56, 3582
Phys. Rev. B 54, 9509 (1996). (1987).
[16] C. J. Halboth and W. Metzner, d-wave Superconductivity [37] L. F. Tocchio, F. Becca, A. Parola, and S. Sorella, Role of
and Pomeranchuk Instability in the Two-Dimensional Backflow Correlations for the Nonmagnetic Phase of the
Hubbard Model, Phys. Rev. Lett. 85, 5162 (2000). t − t0 Hubbard Model, Phys. Rev. B 78, 041101(R) (2008).
[17] P. W. Anderson, New Approach to the Theory of Super- [38] N. Trivedi and D. M. Ceperley, Ground-State Correlations
exchange Interactions, Phys. Rev. 115, 2 (1959). of Quantum Antiferromagnets: A Green-Function
[18] Y. Nagaoka, Ferromagnetism in a Narrow, Almost Half- Monte Carlo Study, Phys. Rev. B 41, 4552 (1990).
Filled s Band, Phys. Rev. 147, 392 (1966). [39] D. F. B. ten Haaf, H. J. M. van Bemmel, J. M. J. van
[19] D Scalapino, Numerical Studies of the 2D Hubbard Model, Leeuwen, W. van Saarloos, and D. M. Ceperley, Proof
in Handbook of High-Temperature Superconductivity, for an Upper Bound in Fixed-Node Monte Carlo for
edited by J. Schrieffer and J. Brooks (Springer, New York, Lattice Fermions, Phys. Rev. B 51, 13039 (1995).
2007), pp. 495–526. [40] J. Paldus and X. Li, A Critical Assessment of Coupled
[20] J. Gukelberger, L. Huang, and P. Werner, On the Dangers Cluster Method in Quantum Chemistry, Adv. Chem. Phys.
of Partial Diagrammatic Summations: Benchmarks for the 110, 1 (1999).
Two-Dimensional Hubbard Model in the Weak-Coupling [41] R. J. Bartlett and M. Musiał, Coupled-Cluster Theory in
Regime, Phys. Rev. B 91, 235114 (2015). Quantum Chemistry, Rev. Mod. Phys. 79, 291 (2007).
[21] S. Zhang, Auxiliary-Field Quantum Monte Carlo for [42] I. Shavitt and R. J. Bartlett, Many-Body Methods in
Correlated Electron Systems, in Emergent Phenomena Chemistry and Physics (Cambridge University Press,
New York, 2009).
in Correlated Matter Modeling and Simulation, edited by
[43] R. Rodríguez-Guzmán, K. W. Schmid, C. A. Jiménez-
E. Pavarini, E. Koch, and U. Schollwöck (Verlag des
Hoyos, and G. E. Scuseria, Symmetry-Projected Varia-
Forschungszentrums Jülich, 2013), Vol. 3.
tional Approaches for Ground and Excited States of the
[22] H. Shi and S. Zhang, Infinite Variance in Fermion
One-Dimensional Hubbard Model, Phys. Rev. B 85,
Quantum Monte Carlo Calculations, arXiv:1511.04084.
245130 (2012).
[23] H. Shi and S. Zhang, Symmetry in Auxiliary-Field Quan-
[44] R. Rodríguez-Guzmán, C. A. Jiménez-Hoyos, R. Schutski,
tum Monte Carlo Calculations, Phys. Rev. B 88, 125132
and G. E. Scuseria, Multireference Symmetry-Projected
(2013).
Variational Approaches for Ground and Excited States of
[24] N. V. Prokof’ev and B. V. Svistunov, Polaron Problem by
the One-Dimensional Hubbard Model, Phys. Rev. B 87,
Diagrammatic Quantum Monte Carlo, Phys. Rev. Lett. 81, 235129 (2013).
2514 (1998). [45] J. P. F. LeBlanc et al., Github Repository, doi:10.5281/
[25] A. Mishchenko, N. Prokofev, A. Sakamoto, and B. zenodo.34433 (2015).
Svistunov, Diagrammatic Quantum Monte Carlo Study [46] See Supplemental Material at http://link.aps.org/
of the Fröhlich Polaron, Phys. Rev. B 62, 6317 (2000). supplemental/10.1103/PhysRevX.5.041041 for complete
[26] K. Van Houcke, E. Kozik, N. Prokof’ev, and B. Svistunov, data sets.
Diagrammatic Monte Carlo, Phys. Procedia 6, 95 (2010). [47] R. Bulla, T. A. Costi, and T. Pruschke, Numerical Re-
[27] A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, normalization Group Method for Quantum Impurity
Dual Fermion Approach to Nonlocal Correlations in the Systems, Rev. Mod. Phys. 80, 395 (2008).
Hubbard Model, Phys. Rev. B 77, 033101 (2008). [48] E. Dagotto, Correlated Electrons in High-Temperature
[28] G. Knizia and G. Kin Lic Chan, Density Matrix Embed- Superconductors, Rev. Mod. Phys. 66, 763 (1994).
ding: A Simple Alternative to Dynamical Mean-Field [49] M van Schilfgaarde, T. Kotani, and S. Faleev, Quasipar-
Theory, Phys. Rev. Lett. 109, 186404 (2012). ticle Self-Consistent GW Theory, Phys. Rev. Lett. 96,
[29] G. Knizia and G. Kin Lic Chan, Density Matrix Embed- 226402 (2006).
ding: A Strong-Coupling Quantum Embedding Theory, J. [50] M. S. Hybertsen and S. G. Louie, Electron Correlation in
Chem. Theory Comput. 9, 1428 (2013). Semiconductors and Insulators: Band Gaps and Quasi-
[30] S. R. White, Density Matrix Formulation for Quantum particle Energies, Phys. Rev. B 34, 5390 (1986).
Renormalization Groups, Phys. Rev. Lett. 69, 2863 (1992). [51] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar,
[31] S. R. White, Density-Matrix Algorithms for Quantum Monte Carlo Calculations of Coupled Boson-Fermion
Renormalization Groups, Phys. Rev. B 48, 10345 (1993). Systems. I, Phys. Rev. D 24, 2278 (1981).
[32] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, [52] G. Sugiyama and S. E. Koonin, Auxiliary Field Monte-
Quantum Cluster Theories, Rev. Mod. Phys. 77, 1027 Carlo for Quantum Many-Body Ground States, Ann. Phys.
(2005). (N.Y.) 168, 1 (1986).
[33] C. Gros, Superconductivity in Correlated Wave Functions, [53] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E.
Phys. Rev. B 38, 931(R) (1988). Gubernatis, and R. T. Scalettar, Numerical Study of the
[34] C. Gros, Physics of Projected Wavefunctions, Ann. Phys. Two-Dimensional Hubbard Model, Phys. Rev. B 40, 506
(N.Y.) 189, 53 (1989). (1989).
041041-25
[54] H. Shi, S. Chiesa, and S. Zhang, Ground-State Properties [73] P. Ring and P. Schuck, The Nuclear Many-Body Problem
of Strongly Interacting Fermi Gases in Two Dimensions, (Springer-Verlag, Berlin, 1980).
Phys. Rev. A 92, 033603 (2015). [74] K. W. Schmid, On the Use of General Symmetry-Projected
[55] S. Zhang and H. Krakauer, Quantum Monte Carlo Method Hartree–Fock–Bogoliubov Configurations in Variational
Using Phase-Free Random Walks with Slater Determi- Approaches to the Nuclear Many-Body Problem, Prog.
nants, Phys. Rev. Lett. 90, 136401 (2003). Part. Nucl. Phys. 52, 565 (2004).
[56] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, [75] N. Tomita, Many-Body Wave Functions Approximated by
D. J. Scalapino, and R. L. Sugar, Sign Problem in the the Superposition of Spin-Projected Nonorthogonal Slater
Numerical Simulation of Many-Electron Systems, Phys. Determinants in the Resonating Hartree–Fock Method,
Rev. B 41, 9301 (1990). Phys. Rev. B 69, 045110 (2004).
[57] J. E. Hirsch, Two-Dimensional Hubbard Model: Numeri- [76] R. Rodríguez-Guzmán, C. A. Jiménez-Hoyos, and G. E.
cal Simulation Study, Phys. Rev. B 31, 4403 (1985). Scuseria, Variational Description of the Ground State of
[58] S. Zhang, J. Carlson, and J. E. Gubernatis, Constrained the Repulsive Two-Dimensional Hubbard Model in Terms
Path Monte Carlo Method for Fermion Ground States, of Nonorthogonal Symmetry-Projected Slater Determi-
Phys. Rev. B 55, 7464 (1997). nants, Phys. Rev. B 90, 195110 (2014).
[59] H. Nguyen, H. Shi, J. Xu, and S. Zhang, CPMC-Lab: A [77] T. Mizusaki and M. Imada, Quantum-Number Projection
Matlab Package for Constrained Path Monte Carlo in the Path-Integral Renormalization Group Method,
Calculations, Comput. Phys. Commun. 185, 3344 (2014). Phys. Rev. B 69, 125110 (2004).
[60] C. Lin, F. H. Zong, and D. M. Ceperley, Twist-Averaged [78] D. Tahara and M. Imada, Variational Monte Carlo Method
Boundary Conditions in Continuum Quantum Monte Carlo Combined with Quantum-Number Projection and Multi-
Algorithms, Phys. Rev. E 64, 016702 (2001). Variable Optimization, J. Phys. Soc. Jpn. 77, 114701 (2008).
[61] C.-C. Chang and S. Zhang, Spatially Inhomogeneous [79] C. A. Jiménez-Hoyos and G. E. Scuseria, Cluster-Based
Phase in the Two-Dimensional Repulsive Hubbard Model, Mean-Field and Perturbative Description of Strongly
Phys. Rev. B 78, 165101 (2008). Correlated Fermion Systems: Application to the One-
[62] C.-C. Chang and S. Zhang, Spin and Charge Order in the and Two-Dimensional Hubbard Model, Phys. Rev. B
Doped Hubbard Model: Long-Wavelength Collective 92, 085101 (2015).
Modes, Phys. Rev. Lett. 104, 116402 (2010). [80] P. Corboz, T. M. Rice, and M. Troyer, Competing States in
[63] D. A. Huse, Ground-State Staggered Magnetization of the t-j Model: Uniform d-wave State versus Stripe State,
Two-Dimensional Quantum Heisenberg Antiferromag- Phys. Rev. Lett. 113, 046402 (2014).
nets, Phys. Rev. B 37, 2380 (1988). [81] P. Corboz, S. R. White, G. Vidal, and M. Troyer, Stripes in
[64] H. Neuberger and T. Ziman, Finite-Size Effects in the Two-Dimensional t-J Model with Infinite Projected
Heisenberg Antiferromagnets, Phys. Rev. B 39, 2608 Entangled-Pair States, Phys. Rev. B 84, 041108 (2011).
(1989). [82] F. Petit and M. Roger, Coupled-Cluster Method Applied to
[65] D. S. Fisher, Universality, Low-Temperature Properties, the Motion of a Single Hole in a Hubbard Antiferromag-
and Finite-Size Scaling in Quantum Antiferromagnets, net, Phys. Rev. B 49, 3453 (1994).
Phys. Rev. B 39, 11783 (1989). [83] R. F. Bishop, Y. Xian, and C. Zeng, A Microscopic
[66] S. R. White and D. J. Scalapino, Density Matrix Renorm- Coupled-Cluster Treatment of Electronic Correlations in
alization Group Study of the Striped Phase in the 2D t-J Hubbard Models, Int. J. Quantum Chem. 55, 181 (1995).
Model, Phys. Rev. Lett. 80, 1272 (1998). [84] Y. Asai and H. Katagiri, Coupled-Cluster Approach to
[67] H. Shi, C. A. Jiménez-Hoyos, R. Rodríguez-Guzmán, G. Electron Correlations in the Two-Dimensional Hubbard
E. Scuseria, and S. Zhang, Symmetry-Projected Wave Model, Phys. Rev. B 60, R13946 (1999).
Functions in Quantum Monte Carlo Calculations, Phys. [85] G. D. Purvis and R. J. Bartlett, A Full Coupled-Cluster
Rev. B 89, 125129 (2014). Singles and Doubles Model—The Inclusion of Discon-
[68] M. Capello, F. Becca, M. Fabrizio, S. Sorella, and E. nected Triples, J. Chem. Phys. 76, 1910 (1982).
Tosatti, Variational Description of Mott Insulators, Phys. [86] G. E. Scuseria, A. Sc. Scheiner, T. J. Lee, J. E. Rice, and
Rev. Lett. 94, 026406 (2005). H. F. Schaeffer III, The Closed-Shell Coupled Cluster
[69] M. Capello, F. Becca, S. Yunoki, and S. Sorella, Uncon- Single and Double (CCSD) Excitation Model for the
ventional Metal-Insulator Transition in Two Dimensions, Description of Electron Correlation: A Comparison with
Phys. Rev. B 73, 245116 (2006). Configuration Interaction (CISD) Results, J. Chem. Phys.
[70] F. Becca, M. Capone, and S. Sorella, Spatially Homo- 86, 2881 (1987).
geneous Ground State of the Two-Dimensional Hubbard [87] J. Xue, C.-C. Cheng, E. J. Walter, and S. Zhang, Spin- and
Model, Phys. Rev. B 62, 12700 (2000). Charge-Density Waves in the Hartree-Fock Ground State
[71] L. F. Tocchio, F. Becca, and C. Gros, Backflow Correla- of the Two-Dimensional Hubbard Model, J. Phys. Con-
tions in the Hubbard Model: An Efficient Tool for the dens. Matter 23, 505601 (2011).
Study of the Metal-Insulator Transition and the Large-U [88] N. C. Handy, J. A. Pople, M. Head-Gordon, K.
Limit, Phys. Rev. B 83, 195138 (2011). Raghavachari, and G. W. Trucks, Size-Consistent Brueck-
[72] A. Leprévost, O. Juillet, and R. Frésard, Intertwined ner Theory Limited to Double Substitutions, Chem. Phys.
Orders from Symmetry Projected Wavefunctions of Lett. 164, 185 (1989).
Repulsively Interacting Fermi Gases in Optical Lattices, [89] M. J. Frisch et al., Gaussian Development Version,
New J. Phys. 17, 103023 (2015). Revision H.21 (Gaussian Inc., Wallingford, CT, 2010).
041041-26
[90] M. Kállay, Z. Rolik, J. Csontos, I. Ladjánszki, L. Szegedy, [108] M. H. Hettler, M. Mukherjee, M. Jarrell, and H. R.
B. Ladóczki, and G. Samu, www.mrcc.hu; See also, Z. Krishnamurthy, Dynamical Cluster Approximation: Non-
Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. local Dynamics of Correlated Electron Systems, Phys.
Kállay, An Efficient Linear-Scaling CCSD(T) Method Based Rev. B 61, 12739 (2000).
on Local Natural Orbitals, J. Chem. Phys. 139, 094105 [109] A. Georges and G. Kotliar, Hubbard Model in Infinite
(2013). Dimensions, Phys. Rev. B 45, 6479 (1992).
[91] F. Verstraete and J. I. Cirac, Renormalization Algorithms [110] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,
for Quantum-Many Body Systems in Two and Higher Dynamical Mean-Field Theory of Strongly Correlated
Dimensions, arXiv:cond-mat/0407066. Fermion Systems and the Limit of Infinite Dimensions,
[92] G. Vidal, Entanglement Renormalization, Phys. Rev. Lett. Rev. Mod. Phys. 68, 13 (1996).
99, 220405 (2007). [111] A. I. Lichtenstein and M. I. Katsnelson, Antiferromagnet-
[93] S. R. White and A. L. Chernyshev, Neel Order in Square ism and d-wave Superconductivity in Cuprates: A Cluster
and Triangular Lattice Heisenberg Models, Phys. Rev. Dynamical Mean-Field Theory, Phys. Rev. B 62, R9283
Lett. 99, 127004 (2007). (2000).
[94] E. M. Stoudenmire and S. R. White, Studying Two- [112] G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli,
Dimensional Systems with the Density Matrix Renormal- Cellular Dynamical Mean Field Approach to Strongly
ization Group, Annu. Rev. Condens. Matter Phys. 3, 111 Correlated Systems, Phys. Rev. Lett. 87, 186401 (2001).
(2012). [113] S. Fuchs, E. Gull, L. Pollet, E. Burovski, E. Kozik, T.
[95] A. W. Sandvik, Finite-Size Scaling of the Ground-State Pruschke, and M. Troyer, Thermodynamics of the 3D
Parameters of the Two-Dimensional Heisenberg Model, Hubbard Model on Approaching the Néel Transition,
Phys. Rev. B 56, 11678 (1997). Phys. Rev. Lett. 106, 030401 (2011).
[96] I. Peschel, Special Review: Entanglement in Solvable [114] J. P. F. LeBlanc and E. Gull, Equation of State of the
Many-Particle Models, Braz. J. Phys. 42, 267 (2012). Fermionic Two-Dimensional Hubbard Model, Phys. Rev.
[97] G. H. Booth and G. Kin-Lic Chan, Spectral Functions of B 88, 155108 (2013).
Strongly Correlated Extended Systems via an Exact [115] E. Gull, P. Werner, O. Parcollet, and M. Troyer, Continuous-
Quantum Embedding, Phys. Rev. B 91, 155107 (2015). Time Auxiliary-Field Monte Carlo for Quantum Impurity
[98] I. W. Bulik, G. E. Scuseria, and J. Dukelsky, Density Models, Europhys. Lett. 82, 57003 (2008).
Matrix Embedding from Broken Symmetry Lattice Mean [116] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M.
Fields, Phys. Rev. B 89, 035140 (2014). Troyer, and P. Werner, Continuous-Time Monte Carlo
[99] Q. Chen, G. H. Booth, S. Sharma, G. Knizia, and G. Kin Methods for Quantum Impurity Models, Rev. Mod. Phys.
Lic Chan, Intermediate and Spin-Liquid Phase of the Half- 83, 349 (2011).
Filled Honeycomb Hubbard model, Phys. Rev. B 89, [117] E. Gull, P. Staar, S. Fuchs, P. Nukala, M. S. Summers, T.
165134 (2014). Pruschke, T. C. Schulthess, and T. Maier, Submatrix
[100] P.-G. de Gennes, Superconductivity of Metals and Alloys Updates for the Continuous-Time Auxiliary-Field Algo-
(Benjamin, New York, 1966). rithm, Phys. Rev. B 83, 075122 (2011).
[101] D. Yamaki, T. Ohsaku, H. Nagao, and K. Yamaguchi, [118] B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire,
Formulation of Unrestricted and Restricted Hartree- S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler
Fock-Bogoliubov Equations, Int. J. Quantum Chem. 96, et al., The ALPS Project Release 2.0: Open Source
10 (2004). Software for Strongly Correlated Systems, J. Stat. Mech.
[102] G. Kin-Lic Chan and M. Head-Gordon, Highly Correlated Theor. Exp. (2011) P05001.
Calculations with a Polynomial Cost Algorithm: A Study [119] E. Gull, P. Werner, S. Fuchs, B. Surer, T. Pruschke, and M.
of the Density Matrix Renormalization Group, J. Chem. Troyer, Continuous-Time Quantum Monte Carlo Impurity
Phys. 116, 4462 (2002). Solvers, Comput. Phys. Commun. 182, 1078 (2011).
[103] G. K. L. Chan, An Algorithm for Large Scale Density [120] A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein,
Matrix Renormalization Group Calculations, J. Chem. Continuous-Time Quantum Monte Carlo Method for
Phys. 120, 3172 (2004). Fermions, Phys. Rev. B 72, 035122 (2005).
[104] G. Kin-Lic Chan and S. Sharma, The Density Matrix [121] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I.
Renormalization Group in Quantum Chemistry, Annu. Krivenko, L. Messio, and P. Seth, TRIQS: A Toolbox
Rev. Phys. Chem. 62, 465 (2011). for Research on Interacting Quantum Systems, Comput.
[105] Ö. Legeza and G. Fáth, Accuracy of the Density-Matrix Phys. Commun. 196, 398 (2015).
Renormalization-Group Method, Phys. Rev. B 53, 14349 [122] A. N. Rubtsov, M. I. Katsnelson, A. I. Lichtenstein, and A.
(1996). Georges, Dual Fermion Approach to the Two-Dimensional
[106] B.-X. Zheng and G. Kin-Lic Chan, Ground-State Phase Hubbard Model: Antiferromagnetic Fluctuations and
Diagram of the Square Lattice Hubbard Model from Fermi Arcs, Phys. Rev. B 79, 045133 (2009).
Density Matrix Embedding Theory, arXiv:1504.01784. [123] H. Hafermann, F. Lechermann, A. N. Rubtsov, M. I.
[107] M. H. Hettler, A. N. Tahvildar-Zadeh, M. Jarrell, T. Katsnelson, A. Georges, and A. I. Lichtenstein, Modern
Pruschke, and H. R. Krishnamurthy, Nonlocal Dynamical Theories of Many-Particle Systems in Condensed Matter
Correlations of Strongly Interacting Electron Systems, Physics, in Modern Theories of Many-Particle Systems in
Phys. Rev. B 58, R7475 (1998). Condensed Matter Physics, Lecture Notes in Physics,
041041-27
edited by D. C. Cabra, A. Honecker, and P. Pujol (Springer, [135] S. A. Kulagin, N. Prokofev, O. A. Starykh, B. Svistunov,
Berlin-Heidelberg, 2012), Vol. 843, pp. 145–214. and C. N. Varney, Bold Diagrammatic Monte Carlo
[124] A. Toschi, A. A. Katanin, and K. Held, Dynamical Vertex Method Applied to Fermionized Frustrated Spins, Phys.
Approximation: A Step Beyond Dynamical Mean-Field Rev. Lett. 110, 070601 (2013).
Theory, Phys. Rev. B 75, 045118 (2007). [136] S. Kulagin, N. Prokof’ev, O. A. Starykh, B. V. Svistunov,
[125] C. Taranto, S. Andergassen, J. Bauer, K. Held, A. Katanin, and C. N. Varney, Bold Diagrammatic Monte Carlo
W. Metzner, G. Rohringer, and A. Toschi, From Infinite to Technique for Frustrated Spin Systems, Phys. Rev. B
Two Dimensions through the Functional Renormalization 87, 024407 (2013).
Group, Phys. Rev. Lett. 112, 196402 (2014). [137] A. S. Mishchenko, N. Nagaosa, and N. Prokofev, Dia-
[126] N. Wentzell, C. Taranto, A. Katanin, A. Toschi, and S. grammatic Monte Carlo Method for Many-Polaron Prob-
Andergassen, Correlated Starting Points for the Functional lems, Phys. Rev. Lett. 113, 166402 (2014).
Renormalization Group, Phys. Rev. B 91, 045120 (2015). [138] K. Van Houcke, F. Werner, N. Prokof’ev, and B. Svistunov,
[127] A. E. Antipov, E. Gull, and S. Kirchner, Critical Ex- Bold Diagrammatic Monte Carlo for the Resonant Fermi
ponents of Strongly Correlated Fermion Systems from Gas, arXiv:1305.3901.
Diagrammatic Multiscale Methods, Phys. Rev. Lett. 112, [139] Y. Deng, E. Kozik, N. V. Prokof’ev, and B. V. Svistunov,
226401 (2014). Emergent BCS Regime of the Two-Dimensional Fermionic
[128] J. Otsuki, H. Hafermann, and A. I. Lichtenstein, Super- Hubbard Model: Ground-State Phase Diagram, Euro-
conductivity, Antiferromagnetism, and Phase Separation phys. Lett. 110, 57001 (2015).
in the Two-Dimensional Hubbard Model: A Dual-Fermion [140] E. Kozik, M. Ferrero, and A. Georges, Nonexistence of the
Approach, Phys. Rev. B 90, 235132 (2014). Luttinger-Ward Functional and Misleading Convergence
[129] A. E. Antipov, J. P. F. LeBlanc, and E. Gull, Opendf-An of Skeleton Diagrammatic Series for Hubbard-Like
Implementation of the Dual Fermion Method for Strongly Models, Phys. Rev. Lett. 114, 156402 (2015).
Correlated Systems, Phys. Procedia 68, 43 (2015). [141] S. Sorella, Linearized Auxiliary Fields Monte Carlo Tech-
[130] A. L. Fetter and J. D. Walecka, Quantum Theory of Many- nique: Efficient Sampling of the Fermion Sign, Phys. Rev.
Particle Systems (Dover Publications Inc., New York, 2003). B 84, 241110 (2011).
[131] N. Prokofev and B. Svistunov, Bold Diagrammatic [142] E. Burovski, N. Prokof’ev, B. Svistunov, and M. Troyer,
Monte Carlo Technique: When the Sign Problem Is The Fermi-Hubbard Model at Unitarity, New J. Phys. 8,
Welcome, Phys. Rev. Lett. 99, 250201 (2007). 153 (2006).
[132] K. Van Houcke, F. Werner, E. Kozik, N. Prokofev, B. [143] E. Kozik, E. Burovski, V. W. Scarola, and M. Troyer, Néel
Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Temperature and Thermodynamics of the Half-Filled Three-
Schirotzek, and M. W. Zwierlein, Feynman Diagrams Dimensional Hubbard Model by Diagrammatic Determi-
versus Fermi-Gas Feynman Emulator, Nat. Phys. 8, 366 nant Monte Carlo, Phys. Rev. B 87, 205102 (2013).
(2012). [144] T. Schäfer, F. Geles, D. Rost, G. Rohringer, E. Arrigoni, K.
[133] E. Gull, D. R. Reichman, and A. J. Millis, Bold-Line Held, N. Blümer, M. Aichhorn, and A. Toschi, Fate of the
Diagrammatic Monte Carlo Method: General Formu- False Mott-Hubbard Transition in Two Dimensions, Phys.
lation and Application to Expansion Around the Non- Rev. B 91, 125109 (2015).
crossing Approximation, Phys. Rev. B 82, 075109 (2010). [145] T. Schaefer, A. Toschi, and K. Held, Dynamical Vertex
[134] E. Kozik, K. Van Houcke, E. Gull, L. Pollet, N. Prokof’ev, Approximation for the Two-Dimensional Hubbard Model,
B. Svistunov, and M. Troyer, Diagrammatic Monte Carlo J. Magn. Magn. Mater., doi: 10.1016/j.jmmm.2015.07.103
for Correlated Fermions, Europhys. Lett. 90, 10004 (2010). (2015).
041041-28

PhysRevX 5 041041

Uploaded by

Copyright:

Available Formats

PhysRevX 5 041041

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

PhysRevX 5 041041

Uploaded by

Copyright:

Available Formats

PHYSICAL REVIEW X 5, 041041 (2015)

Solutions of the Two-Dimensional Hubbard Model: Benchmarks and Results from

(Simons Collaboration on the Many-Electron Problem)

DOI: 10.1103/PhysRevX.5.041041 Subject Areas: Computational Physics,

I. INTRODUCTION systems of large numbers of interacting electrons is one

2160-3308=15=5(4)=041041(28) 041041-1 Published by the American Physical Society

schemes produce consistent results within statistical errors. FN

For dilute systems, our benchmarks include points (U=t ¼ 6, _

n ¼ 0.6), (U=t ¼ 6, n ¼ 0.3), and (U=t ¼ 8, n ¼ 0.3). 8

order of 10−3 ) can be systematically reduced with addi-

obtained with variational Monte Carlo, we show finite-size

6 × 6, 8 × 8, and 10 × 10, and exhibit weak finite-size 2x2

effects at this U value. We see that the result is accurate xxxx

n=0.875, U=8 0.04

results are not available at the lowest benchmark temper-

the value from UCCSD and gives E=t ¼ −0.7272

the thermodynamic limit are extrapolated from finite

shows surprisingly large finite-size effects which persist

10−4 t, the spread in values results in a large uncertainty

are shown in the inset of Fig. 9. AFQMC provides 0.126 6

numerically exact ground-state energies for this system. 4 8x8

The value quoted is E=t ¼ −0.8603ð2Þ, which is in agree- 0.125

[E=t ¼ −0.8575ð3Þ], and unlike in previous cases for

The UCCSD thermodynamic limit overestimates the

T=0.25 -1.06 -1.08 10x8

U=t ¼ 4, in order for the finite-T data to be consistent with

U=2, T=0.5, n=0.3, t’=0

TABLE VI. Comparison of magnetization data from DMET and AFQMC at n ¼ 1.

with error bars, for the thermodynamic-limit values of ACKNOWLEDGMENTS

You might also like