# Multiplicity, localization, and domains in the Hartree-Fock ground state of
the two-dimensional Hubbard model

Kazue Matsuyama<EMAIL_ADDRESS>Jeff Greensite<EMAIL_ADDRESS>Physics and
Astronomy Department, San Francisco State University,
1600 Holloway Ave, San Francisco, CA 94132 USA

###### Abstract

We explore certain properties of the Hartree-Fock approximation to the ground
state of the two-dimensional Hubbard model, emphasizing the fact that in the
Hartree approach there is an enormous multiplicity of self-consistent
solutions which are nearly degenerate in energy, reminiscent of a spin glass,
but which may differ substantially in other bulk properties. It is argued that
this multiplicity is physically relevant at low temperatures. We study the
localization properties of the one-particle wavefunctions comprising the
Hartree-Fock states, and find that these are unlocalized at small and moderate
values of $U/t$, in particular in the stripe region, but become highly
localized at values corresponding to strong repulsion. We also find
rectangular domains as well as stripes in the stripe region of the phase
diagram, and study pair correlations in the neighborhood of half-filling.

###### keywords:

Hubbard model, Hartree-Fock, localization, strongly correlated systems

††journal: Annals of Physics

##  1  Introduction

It is well known that the wave function of a particle in a stochastic
potential may be localized; this is called Anderson localization. One of the
questions we would like to address here is whether, in a many-body quantum
system with a translation-invariant Hamiltonian, such a phenomenon could also
occur at zero or very low temperature. The rough idea is that the potential
seen by one particle, due to all particles in the system, appears to be
disordered, perhaps sufficiently so as to induce localization. It seems
obvious that if a translation invariant Hamiltonian in a finite volume has a
unique ground state, then spatial localization in that state is impossible;
the ground state would have to share the translation invariance of the
Hamiltonian, and the expectation value of any observable must likewise be
translation invariant. But a different, spin-glass scenario is possible, in
which the low temperature behavior is not necessarily dominated by the lowest
energy state. It may be that the landscape consists of very many states which,
as in a spin glass, are nearly degenerate in energy. As a many-body system is
cooled it may become stuck in one of these states, which need not be
translation invariant, and in this scenario some analog of Anderson
localization could be possible. This would answer in the affirmative the
question raised in ref. [1], of whether localization could occur in a system
described by a translation-invariant many-body Hamiltonian, with a stochastic
element introduced only in the initialization.

The picture of one particle moving in a random potential due to all other
particles of the system immediately suggests a Hartree-Fock approach, and the
two-dimensional Hubbard model, which has been extensively studied by the
Hartree-Fock method, see e.g. [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18] seems like a good framework in which to address the questions we
are interested in. It is known [10, 11, 16] that there is no unique solution
of the self-consistent Hartree-Fock equations; rather there is great
multiplicity of self-consistent solutions generated by an iterative procedure.
Some authors (e.g. [16]) look for states of lower energy out this multiplicity
of states via some simulated annealing procedure. We believe this is not the
right approach if the system at very low temperature exhibits localization or
translation-symmetry breaking of any sort, and a multiplicity of states is
physically relevant. The situation is then better described by analogy to spin
glasses, where the search for the global minimum of the energy seems futile,
and we accept that the system, cooled from some finite temperature, will get
stuck in one or another local minimum. It is possible that none of the
numerous self-consistent solutions in the Hartree-Fock approach represent the
true, translation-invariant ground state, but rather correspond to very long-
lived metastable states in the exact theory. The example of the double well
potential with a very high barrier between the two minima comes to mind. The
true ground state has the symmetry of the Hamiltonian, but there exist states
which remain concentrated in one of the two wells for a long period which
depends on the tunneling amplitude. In case of 2D many body systems with a
translation invariant Hamiltonian, if spatial dependence in some observables
is seen experimentally at very low temperatures, then it must be that a
multiplicity of either stable, or long-lived metastable states, such as the
states identified by the Hartree-Fock procedure, are the states which are
physically relevant.

Localization of one particle wavefunctions in the Slater determinant has not,
to our knowledge, been studied in the context of the Hartree-Fock
approximation to the standard Hubbard model (i.e. the model with no additional
random potential). In fact the standard 2D Hubbard model is known to exhibit
periodic patterns such as stripes and a checkerboard in certain regions of the
phase diagram, and it is natural to attribute those patterns to localized spin
up/down electrons at each site. As we will see, this is not the case; electron
localization in the Hubbard model, as quantified by the inverse participation
ratio (IPR), appears only at rather large values of $U/t$.

Apart from localization there are of course many other observables which can
be computed in the framework of the Hartree-Fock approach, and in this article
we will display a selection which we feel are either relevant to the
multiplicity of solutions, or which may not have been emphasized in previous
studies. In particular we point out the existence of rectangular domains in
spin density distributions, in addition to the stripe patterns noted long ago.

We motivate the Hartree-Fock approximation to the Hubbard Hamiltonian in a
slightly unconventional manner, drawn from our picture of a single electron
moving in a stochastic potential generated by all other electrons. Let us
suppose that we have $M$ electrons on a two dimensional ${L\times L}$ square
lattice with periodic boundary conditions, interacting according to the usual
Hubbard model Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle-t\sum_{<xy>}\sum_{s=\uparrow,\downarrow}c^{\dagger}(x,s)c(y,s)$
(1)
$\displaystyle+U\sum_{x}c^{\dagger}(x,\uparrow)c(x,\uparrow)c^{\dagger}(x,\downarrow)c(x,\downarrow)\
,$

with the first term a sum over nearest neighors, and define the Hartree-Fock
state $|\Phi\rangle=\Phi|0\rangle$ where

$\Phi=\prod_{i=1}^{M}\sum_{\mathbf{x}_{i}}\phi_{i}(\mathbf{x}_{i})c^{\dagger}(\mathbf{x}_{i})\
,$ (2)

and where the single-particle states $\\{\phi_{i}\\}$ are all orthogonal. Our
notation is that ${\bf x}=(x,s)$ where $x$ labels the lattice site, and
${s=\\{\uparrow,\downarrow\\}}$ is the spin. Now we introduce an operator for
an auxiliary (and fictitious) electron interacting with the other $M$
electrons, and define the reduced one-particle Hamiltonian $H[\Phi]$ with
matrix elements

$H_{\mathbf{w}_{1},\mathbf{w}_{2}}[\Phi]=\langle
0|\Phi^{\dagger}c(\mathbf{w}_{1})Hc^{\dagger}(\mathbf{w}_{2})\Phi|0\rangle_{conn}\
,$ (3)

where the subscript $conn$ means that we keep only contributions to the
expectation value which contain the anticommutators of $c(\mathbf{w}_{1})$ and
$c^{\dagger}(\mathbf{w}_{2})$ with creation/destruction operators in $H$. This
is our candidate for the Hamiltonian describing the propagation of a single
electron in the average field of all other electrons. The approximation to the
Hubbard model ground state is arrived at by an iterative procedure. Given
$\Phi^{(n)}$ at the $n$-th iteration, we solve numerically the eigenvalue
problem

$H[\Phi^{(n)}]\phi_{i}({\bm{x}})=\epsilon_{i}\phi_{i}({\bm{x}})$ (4)

for the $M$ lowest energy states, and insert those eigenstates
$\\{\phi_{i},~{}i=1,2,...M\\}$ into (2) to obtain the next approximation
$|\Phi^{n+1}\rangle$ to the ground state at the ${(n+1)}$-th iteration. The
procedure terminates when a certain convergence criterion, described below, is
satisfied.

##  2 Procedure

It is straightforward, given $\Phi$, to find the non-zero matrix elements:

$w_{1},w_{2}$ nearest neighbors

$\displaystyle H[\Phi]_{\mathbf{w}_{1},\mathbf{w}_{2}}=-t\delta_{s_{1},s_{2}}$
(5)

and we note again that periodic boundary conditions are imposed.

$w_{1}=w_{2}=w$ same site

$\displaystyle H[\Phi]_{w\uparrow,w\uparrow}$ $\displaystyle=$ $\displaystyle
U\rho(w,\downarrow,\downarrow)$ $\displaystyle
H[\Phi]_{w\downarrow,w\downarrow}$ $\displaystyle=$ $\displaystyle
U\rho(w,\uparrow,\uparrow)$ $\displaystyle H[\Phi]_{w\uparrow,w\downarrow}$
$\displaystyle=$ $\displaystyle-U\rho(w,\downarrow,\uparrow)$ $\displaystyle
H[\Phi]_{w\downarrow,w\uparrow}$ $\displaystyle=$
$\displaystyle-U\rho(w,\uparrow,\downarrow)\ .$ (6)

where

$\displaystyle\rho(w,\downarrow,\downarrow)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}\phi_{i}^{*}(w,\downarrow)\phi_{i}(w,\downarrow)$
$\displaystyle\rho(w,\uparrow,\uparrow)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}\phi_{i}^{*}(w,\uparrow)\phi_{i}(w,\uparrow)$
$\displaystyle\rho(w,\downarrow,\uparrow)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}\phi_{i}^{*}(w,\downarrow)\phi_{i}(w,\uparrow)$
$\displaystyle\rho(w,\uparrow,\downarrow)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}\phi_{i}^{*}(w,\uparrow)\phi_{i}(w,\downarrow)\ .$
(7)

The problem is then to find the eigenvalues and corresponding eigenvectors of
a $2L^{2}\times 2L^{2}$ sparse matrix, which can be handled by standard
numerical software.111We use the Matlab eigs function. It is important to
point out here that the iterative procedure does not converge to a unique
ground state, but depends instead on the (random) initialization. Energy
densities are only very weakly dependent on the initialization; magnetization
and spatial distributions have a far stronger dependence; we will return to
this point below. We initialize the system with a small, random choice of
$\rho$, i.e. at each site $w$ we generate three uniformly distributed random
numbers $r_{1},r_{2},r_{3}$ in the range $[0,1]$, and let

$\displaystyle\rho(w,\uparrow,\uparrow)$ $\displaystyle=$ $\displaystyle
ar_{1}$ $\displaystyle\rho(w,\downarrow,\downarrow)$ $\displaystyle=$
$\displaystyle ar_{2}$ $\displaystyle\rho(w,\uparrow,\downarrow)$
$\displaystyle=$ $\displaystyle\rho(w,\downarrow,\uparrow)=b(r_{3}-0.5)\ ,$
(8)

with $a,b$ taken to be small constants, e.g. $a=0.01,b=0.001$. The choice, and
even the order of magnitude of these constants is not critical. However, while
each set of iterations converges to a solution with very nearly degenerate
energy densities, typically differing by fractional deviations of order
$O(10^{-4})$, solutions obtained with different stochastic starting points
vary widely in the spatial distribution of spin up and spin down electron
densities. If this Hartree-Fock result reflects a true property of the 2D
Hubbard model, then the enormous multiplicity of nearly degenerate ground
states is reminiscent of a spin glass.

After the first iteration, one usually finds equal numbers of electrons of
either spin

$\sum_{w}\rho(w,\uparrow,\uparrow)=\sum_{w}\rho(w,\downarrow,\downarrow)=\frac{1}{2}M\
,$ (9)

but any small deviation from this sum rule is corrected by adding or
subtracting a small constant of order $1/L^{2}$ to the electron densities,
prior to the next iteration. Thus we are selecting for self-consistent
solutions with equal numbers of up and down spins. The convergence criterion
is that after $n$ iterations the mean square deviations of electron spin
densities satisfy

$\displaystyle{1\over
L^{2}}\sum_{w}(\rho^{(n)}(w,\uparrow,\uparrow)-\rho^{(n-9)}(w,\uparrow,\uparrow))^{2}$
$\displaystyle<$ $\displaystyle\delta^{2}$ $\displaystyle{1\over
L^{2}}\sum_{w}(\rho^{(n)}(w,\downarrow,\downarrow)-\rho^{(n-9)}(w,\downarrow,\downarrow))^{2}$
$\displaystyle<$ $\displaystyle\delta^{2}\ ,$ (10)

where the $\rho^{(n)}$ denote densities obtained after the $n$-th iteration.
We choose ${\delta=0.001}$. We have investigated, at some points in the
$U-$density plane, the effect of decreasing $\delta$ by one or two orders of
magnitude. This will be discussed below, but in brief the only effect of
decreasing $\delta$ by an order of magnitude is to increase the number of
iterations required by convergence modestly, at $U/t=3$, or by a factor of 4
or 5 at $U/t=30$. But strengthening the convergence criterion in this way
makes very little difference quantitatively, e.g. on the energy density, and
no difference whatever qualitatively.

In all our numerical computations we have taken $t=1$, so $U/t=U$ in the
results reported below.

###  2.1 Relation to the mean field decomposition

With the subtraction of an (irrelevant) constant, the operator

$H_{mf}=\sum_{\mathbf{w}_{1},\mathbf{w}_{2}}c^{\dagger}(\mathbf{w}_{1})H[\Phi]_{\mathbf{w}_{1},\mathbf{w}_{2}}c(\mathbf{w}_{2})-E_{0}$
(11)

is the standard mean field approximation to the Hubbard model Hamiltonian,
where

$E_{0}=U\sum_{x}\Bigl{(}\rho(x,\uparrow,\uparrow)\rho(x,\downarrow,\downarrow)-\rho(x,\uparrow,\downarrow)\rho(x,\downarrow,\uparrow)\Bigr{)}\
,$ (12)

cf. Fazekas [15], and Lechermann in [19]. Some authors, however, adopt a
simpler mean-field decomposition. Writing

$n(x,\uparrow)n(x,\downarrow)=c^{\dagger}(x,\uparrow)c(x,\uparrow)c^{\dagger}(x,\downarrow)c(x,\downarrow)\
,$ (13)

where

$n(x,s)=c^{\dagger}(x,s)c(x,s)\ ,$ (14)

the simpler decomposition is

$\displaystyle
c^{\dagger}(x,\uparrow)c(x,\uparrow)c^{\dagger}(x,\downarrow)c(y,\downarrow)$
$\displaystyle\qquad\rightarrow n(x,\uparrow)\langle
n(x,\downarrow)\rangle+\langle n(x,\uparrow)\rangle n(x,\downarrow)-\langle
n(x,\uparrow)\rangle\langle n(x,\downarrow)\rangle\ ,$

and in our notation this amounts to dropping the spin-fllip terms
$H[\Phi]_{x\uparrow,x\downarrow}$ and $H[\Phi]_{x\downarrow,x\uparrow}$ in
$H_{mf}$. This is what is done in the Hartree-Fock calculations of ref. [16],
whose approach is in other respects quite similar to our own. Dropping the
spin-flip terms is, however, a severe truncation of the standard mean field
expansion, and is analogous to dropping the exchange term in atomic physics

We also make no attempt to select the “best” self-consistent solution via some
annealing procedure, which aims to select solutions which are slightly lower
in energy than solutions obtained without this procedure. This reflects our
opinion, already discussed above, that a multiplicity of physically relevant
states, analogous to a spin glass, is a requirement if translation non-
invariance is seen at low temperatures, and therefore a state which is very
slightly lower in energy than most Hartree-Fock states (which are all nearly
degenerate in energy) is of no special significance.

##  3 Observables

Once a self-consistent solution is obtained, we compute:

  * 1.

The energy density

$\mathcal{E}={1\over L^{2}}\sum_{i=1}^{M}\epsilon_{i}\ .$ (16)

  * 2.

The energy gap (at zero temperature) between the highest energy occupied and
next unoccupied states

$\Delta\epsilon=\epsilon_{M+1}-\epsilon_{M}\ .$ (17)

  * 3.

Charge and spin densities

$\displaystyle C(x)$ $\displaystyle=$
$\displaystyle\rho(x,\uparrow,\uparrow)+\rho(x,\downarrow,\downarrow)$
$\displaystyle D(x)$ $\displaystyle=$
$\displaystyle\rho(x,\uparrow,\uparrow)-\rho(x,\downarrow,\downarrow)\ .$ (18)

  * 4.

Local magnetization

$\displaystyle m$ $\displaystyle=$ $\displaystyle{1\over
2L^{2}}\sum_{x}\sum_{\mu=1}^{2}D(x)D(x+\hat{\mu})\ .$ (19)

To be clear, $m$ is not the global magnetization; it is simply the corellator
of nearest neighbor spin densities. Antiferromagnetic regions of the lattice
will have negative $m$, ferromagnetic are positive $m$, but it is understood
that non-zero $m$ does not necessarily imply long-range order.

  * 5.

Long range order. In principle any regular arrangement of spins in the 2D
Hubbard model, at any non-zero temperature, is a violation of the Mermin-
Wagner theorem. Nevertheless, such arrangements have been observed on finite
lattices in quantum Monte Carlo simulations at half-filling [20]; this must be
attributed to the very low (zero) temperature and finite volume. A useful
observable to probe periodicity, at least on the scale of the finite lattice,
is

$\displaystyle S(k)$ $\displaystyle=$ $\displaystyle{1\over
L^{2}}\sum_{x,y}D(x)D(y)e^{ik\cdot(x-y)}$ (20) $\displaystyle=$
$\displaystyle{1\over L^{2}}\widetilde{D}(k)\widetilde{D}^{*}(k)\ ,$

where $\widetilde{D}(k)$ is the Fourier transform of the spatial spin
distribution $D(x)$. We have long range order (with the caveat just
mentioned), when $S(k)$ is concentrated at only a few values of $k$, typically
resulting in either a checkerboard, stripe, or rectangular domain pattern.

  * 6.

Localization. This is our original motivation. We can judge whether individual
energy eigenstates $\phi_{i}(x,s)$ are localized by computing the inverse
participation ratio (IPR)

$\displaystyle IPR_{i}=\sum_{x}\sum_{s}(\phi^{*}_{i}(x,s)\phi_{i}(x,s))^{2}\
,$ (21)

where $IPR_{i}=1$ means that the one electron state is localized to a single
site and definite spin, while $IPR_{i}\sim 1/L^{2}$ indicates that the
positional probability density is spread evenly over most of the lattice area.

  * 7.

Momentum distribution. We compute momentum-space occupation numbers

$\displaystyle n(k)$ $\displaystyle=$ $\displaystyle{1\over
2L^{2}}\sum_{x,y,s}e^{ik\cdot(x-y)}\langle\Phi|c^{\dagger}(x,s)c(y,s)|\Phi\rangle$
(22) $\displaystyle=$ $\displaystyle{1\over
2L^{2}}\sum_{i,s}\phi_{i}^{*}(k,s)\phi_{i}(k,s)\ ,$

where

$\phi_{i}(k,s)=\sum_{x}e^{ik\cdot x}\phi_{i}(x,s)\ ,$ (23)

and we also compute the magnitude of the discretized gradient

$\displaystyle\nabla_{x}n(k)$ $\displaystyle=$ $\displaystyle{L\over
4\pi}(n(k_{x}+1,k_{y})-n(k_{x}-1,k_{y}))$ $\displaystyle\nabla_{y}n(k)$
$\displaystyle=$ $\displaystyle{L\over
4\pi}(n(k_{x},k_{y}+1)-n(k_{x},k_{y}-1))$ $\displaystyle|\nabla n(k)|$
$\displaystyle=$ $\displaystyle\sqrt{(\nabla_{x}n)^{2}+(\nabla_{y}n)^{2}}\ .$
(24)

  * 8.

Pairing correlations. We search for d-wave arrangements in the correlator

$\Delta(k^{\prime},k)=\langle
c^{\dagger}(k^{\prime},\uparrow)c^{\dagger}(-k^{\prime},\downarrow)c(k,\uparrow)c(-k,\downarrow)\rangle\
,$ (25)

where $k^{\prime}\neq k$, and choosing the $x,y$ components of $k$ to be the
$y,x$ components of $k^{\prime}$.

##  4 Results

###  4.1 Momentum space distributions

A good starting point is to compare our results for momentum space
distributions $n(k)$ and gradients $|\nabla n(k)|$ with the results of Monte
Carlo simulations, particularly where those Monte Carlo results exist away
from half-filling. Because of the sign problem, Monte Carlo simulations away
from half-filling must rely on a reweighting procedure of some kind [21], and
this apparently does not work down to zero temperature. Thus our comparison
with the Monte Carlo simulations of Varney et. al [20], away from half-
filling, is necessarily a comparison of our zero temperature results with
Monte Carlo simulations at finite temperature. The comparison is nonetheless
interesting, even if only at the qualitative level.

Figure 1: 2D color plots of the momentum distribution $n(k_{x},k_{y})$ and the
magnitude $|\nabla n|$ of its gradient. This plot should be compared with the
Monte Carlo results displayed in Figure 2 of Varney et al. [20]. Densities,
increasing horizontally from left to right, are $f=0.23,0.41,0.61,0.79,1.0$,
where $f$ is the density. (a) $n(k)$ at $U/t=2$; (b) $|\nabla n|$ at $U/t=2$;
(c) $n(k)$ at $U/t=4$; (d) $|\nabla n|$ at $U/t=4$.

In Fig. 1 we display our results (on a $24^{2}$ lattice) for momentum
distributions $n(k)$ and gradients $|\nabla n(k)|$ at the same densities
(electrons/site)
${f=0.23,0.41,0.61,0.79,1.0}$ used in ref. [20], and Fig. 1 should be compared
with Fig. 2 in that reference. Despite the finite temperatures used in [20]
the figures are very similar, even at the quantitative level. Note in
particular the absence of a sharp boundary in momentum space, separating
occupied and unoccupied states at the larger filling fractions, as $U$ is
increased from $U=2$ to $U=4$.

###  4.2 Multiplicity of self-consistent states

One significant aspect of the Hartree-Fock approach to the 2D Hubbard model is
that there is no unique self-consistent solution for the ground state, as
noted many times in the literature. Instead there are very many such states,
obtained from different random starts (8), except, of course, at $U=0$. In
Fig. 2 we show our results at density $f=0.8$ on $24^{2}$ lattices for the
energy density $\epsilon$, the local magnetization $m$, and the energy gap
$\Delta\epsilon$, together with their standard deviations. The standard
deviations are obtained at each $U$ from 20 separate self-consistent states,
obtained as described in section II. These are plotted in Fig. 2 as error
bars, but we stress that in these figures the “error bar” represents the
standard deviation, rather than standard deviations of the mean.

The standard deviation of energy density, $\sigma_{\epsilon}$, is so small
that it is not discernible in Fig. 2(a). As in a spin glass, with a landscape
of local minima of near-degenerate energies, the different self-consistent
solutions are nearly degenerate in energy. The difference between these self-
consistent states is quite apparent, however, in the local magnetization $m$
and the energy gap $\Delta\epsilon$, seen in Figs. 2(b) and 2(c) respectively.
In these cases, at $f=0.8$ and $U\geq 3$ (where local antiferromagnetic order
becomes apparent), the standard deviations are quite substantial, indicating
that despite their near-degeneracy in energy, these different self-consistent
solutions are physically distinct.

(a)

(b)

(c)

Figure 2: Values and standard deviations vs. $U$ at density $f=0.8$ of (a)
energy density $\mathcal{E}$; (b) local magnetization $-m$; and (c) energy gap
$\Delta\epsilon$. Standard deviations, displayed as error bars, were taken
from a set of twenty independent, self-consistent Hartree-Fock solutions. Note
that no deviation is visible among these different solutions in the energy
density, but there are significant deviations in the magnetization and energy
gap for $U\geq 3$.

(a)

(b)

Figure 3: A histogram of energy vs. state number, with the bar for the last
occupied state $\epsilon_{M}$ shown in red. $U=4,f=0.9$ on a $24^{2}$ lattice.
(a) A closeup of energies in the region of the last occupied state; here the
gap to the first unoccupied state is quite clear. (b) Overview of energies
$\epsilon_{i}$ up to, and a little beyond, the last occupied state.

###  4.3 Local magnetization and the energy gap

In general there is a correlation, at zero temperature, between local
magnetization $m$ and the gap $\Delta\epsilon$ between the energies of the
last occupied and first unoccupied states. These gaps increase with $|m|$,
which in turn is largest at or near half-filling, and large $U$. In Fig. 3(a)
we display a histogram of energies $\epsilon_{n}$ vs. $n$ on a $24\times 24$
lattice at $U=4$ and density $f=0.9$, taken from a typical configuration which
has converged from a random initialization, as described above. This is a
closeup view of energies in the neighborhood of $n=M$, with the energy of the
last occupied state ($n=M=518$ for these parameters) shown in a different
color. Note the significant jump in energy between the last occupied state and
the first unoccupied state ($n=M+1$). Figure 3(b) is a more global display of
the energies, beginning at the lowest energy.

Figure 4: Energy density and magnetization vs. density $f$ at ${U=4}$ on a
$24^{2}$ lattice. In this plot we have multiplied the magnetization by $-4$,
for ease of comparison.

(a) $(-1)\times$ local magnetization

(b) energy gap

(c) energy density

Figure 5: Surface plots of (a) local magnetization $\times(-1)$; (b) energy
gap $\Delta\epsilon$; (c) energy density $\mathcal{E}$; vs. coupling $U$ and
density $f$. All data taken on $30^{2}$ lattices. Note, in (a,b), that the
region of significant local antiferromagnetic order is also the region where
the energy gap is significant.

At moderate values of $U$, in a region centered around half-filling ($f=1$),
we find that the local magnetization defined in (19) is significantly non-zero
and negative, indicating local antiferromagnetism. It is interesting that the
local magnetization is closely correlated with the existence of an energy gap
$\Delta\epsilon>0$. An example of this correlation, at $U=4$ on a $24^{2}$
lattice, is shown in Fig. 4. In order that both the magnetization and energy
gap are clearly visible on the same figure, we have multiplied the
magnetization by a factor of $-4$.

Figure 6: Color plot of local magnetization $m$ in the $f-U$ plane up to
$U=12$. The bright region is a region of local antiferromagnetism, surrounded
by a region of negligible magnetization.

Figure 7: Local magnetization $m$ vs. $U$ at fixed filling density $f=0.7$. We
note that the local magnetization switches from antiferromagnetic to
ferromagnetic around $U=11$. The “error bars” shown actually represent the
standard deviation taken from 20 independent, self consistent Hartree-Fock
solutions.

A more complete picture of the local magnetization $m$ and energy gap
$\Delta\epsilon$ in the $U-f$ parameter plane up to $U=6$ is shown in Figs.
5(a) and 5(b). Again we see that local antiferromagnetic order in some region
of the phase diagram is correlated with the existence of an energy gap
$\Delta\epsilon$. The energy density in the $U-f$ plane is shown in Fig. 5(c).
There is no obvious indication in this figure of any non-analyticity
indicating a thermodynamic phase transition, although other types of phase
transitions (see, e.g., [22]) may exist. We note in passing that the energy
density at $U=0$ is symmetric around half-filling, as can be verified
analytically.

Fig. 6 is a 2D color plot of the local magnetization in the $f-U$ plane up to
$U=12$. Data is taken from the average of results obtained from 20 independent
Hartree-Fock ground state solutions on a $30^{2}$ lattice. The brightly
colored oval region is a region of local antiferromagnetism; the local
magnetization is negligible outside this region. At the upper border of the
diagram, at $U=12$, the local magnetization turns positive away from half-
filling, e.g. at $f=0.7,1.3$ which marks the beginning of ferromagnetic
regions in the phase diagram. In Fig. 7 we plot magnetization vs. $U$ at fixed
$f=0.7$, also averaging results taken from 20 solutions on a $30^{2}$ lattice
(again with “error bars” representing standard deviations) and it is clear
that at this density, from $U=11$ onwards, the local magnetization is
ferromagnetic rather than antiferromagnetic. All in all, the picture presented
in Figs. 5 and 6 is consistent, at least qualitatively, with much earlier
Hartree-Fock explorations of the phase structure of the 2D Hubbard model, e.g.
[3]

####  4.3.1 Large $U$ limit

At $U=0$ the energies of one-particle states can be computed analytically,

$\displaystyle\epsilon(m,n)$ $\displaystyle=$ $\displaystyle-2(\cos(2\pi
n/L)+\cos(2\pi m/L))$ (26) $\displaystyle\text{where}~{}~{}m,n=0,1,...,L-1\ ,$

and, after sorting these values from lowest to highest, agreement with the
numerical calculation at $U=0$ is simply a modest check of our code. A more
interesting limit is the computation of energies and energy gap at half-
filling and large $U$. If we ignore the hopping term by comparison with the
potential then the energy of the ground state at the classical level is zero
(each site occupied by a single electron), and the energy of the first excited
state is simply $U$, corresponding to having one doubly-occupied site. The
actual values for the $\epsilon_{i}$ at $f=1$ and $U=30$ are shown in Fig. 8.
The last occupied state on this $24^{2}$ lattice is at state number
$M=24^{2}=576$, and the energy gap in this case is found to be
$\Delta\epsilon=24.6$, which is not so far from $U=30$. The energy density of
the ground state is very nearly zero, due to a near-exact cancellation of the
positive and negative energies of the occupied states. The convergence to
$\mathcal{E}=0$ at half-filling with increasing $U$, and the near linear
increase of $\Delta\epsilon$ with $U$ for $U>4$ is seen in Fig. 9.

Figure 8: One particle energies at $U=30,~{}f=1$. The energy gap between the
last occupied and first unoccupied states is $\Delta\epsilon=24.6$.

Figure 9: Convergence to $\mathcal{E}=0$, and linear growth of the energy gap
$\Delta\epsilon$ with increasing $U$, at half-filling.

Subsequent figures will display the results (spin distribution, IPR values,
etc.) taken from “typical” self-consistent solutions of the Hartree-Fock
equations although, as seen above and also in the next section, “typical”
solutions can vary a great deal in certain observables, spatial ordering in
particular.

###  4.4 Spatial patterns in the spin distribution

Examination of the spatial distribution of charge and spin densities, $C(x)$
and $D(x)$ respectively in (18), reveal some interesting geometric patterns in
the neighborhood of half-filling. The existence of stripes in Hartree-Fock
treatments of the 2D Hubbard model goes back to [4, 5, 6, 7], and there is
also experimental evidence of stripe order in strongly correlated materials
[23, 24]. We will focus here on the spin densities $D(x)$, which are displayed
in Fig. 10 on a $24^{2}$ lattice at $U/t=3$ at zero temperature near half
filling.

At exactly half-filling one observes an antiferromagnetic “checkerboard”
pattern, as seen in Fig. 10(d). Despite this apparently antiferromagnetic
order, one should refrain from the interpretation that each electron is
localized at one lattice site, in a pattern of alternating up/down spins. In
fact this is far from the case, as is best seen from an examination of the IPR
of each $\phi_{i}$ at $U/t=3,~{}f=1$, which is shown in Fig. 11. The great
majority of these IPR values are $\approx 0.0025$, and we recall that IPR=1
means that a particle is localized at a single lattice site, while
IPR$=1/L^{2}$, which is 0.0017 on a $24^{2}$ lattice, is complete
delocalization. It is evident that almost all single electron states extend
over the entire lattice, and the picture of electrons localized at sites in an
up/down alternating pattern, as strongly suggested by the checkerboard
pattern, is completely untenable. It is remarkable that such highly
unlocalized electrons nonetheless contrive to produce such a regular geometric
pattern.

Figure 10: Display of $D(x)$ at densities $f=0.7$ to $f=1.3$ and $U/t=3$. The
checkerboard pattern at $f=1$ (d) indicates an antiferromagnetic order across
the $24^{2}$ lattice area. Within the antiferromagnetic region, e.g. at
$f=0.8$ and $f=1.2$, both “wavelike” in (b) and “domains” shown in (c) and (e)
are found, depending on the initialization. Outside the antiferromagnetic
region there is no obvious pattern, and also the magnitude of $D(x)$ (note the
colorbar scale at (a) and (f)) is reduced by two orders of magnitude.

Figure 11: Localization (IPR values) of one particle energy eigenstates
$\phi_{i}$ vs. $i$ at $U/t=3,f=1$, $24^{2}$ lattice. Despite the checkerboard
pattern seen at these values in Fig. 10, the one particle wavefunctions are
either not at all, or only very weakly localized. Localization is a feature of
strong repulsion (see Fig. 15 below), and is not associated with the stripe,
domain, or checkerboard spatial patterns seen at moderate repulsion.

Stripe order in the spin density observable is seen in Fig. 10(b) at $f=0.8$.
This type of pattern is observed in many of the self-consistent Hartree-Fock
solutions generated by our iterative procedure with random initial conditions,
at moderate couplings in the neighborhood of half-filling. But it is not the
only pattern found. An equally common pattern is the periodic arrangement of
quasi-rectangular domains seen in Fig. 10(c) again at $f=0.8$. We emphasize
that these different patterns are obtained at the same coupling parameters and
at the same zero temperature. This is simply another feature of the
multiplicity of self-consistent Hartree-Fock solutions, as already discussed,
and is again indicative of just how different these different solutions can
be. The possibility of domain structures of this kind was suggested previously
in [25].

We have examined, at $f=0.8$ and $U/t=3$, the sensitivity of our results to
the convergence criterion. Of course the number of iterations required to
reach convergence varies significantly from one initialization to the next,
but the following numbers are typical: 120 iterations to convergence on a
$24^{2}$ lattice at $\delta=10^{-3}$, 170 iterations at $\delta=10^{-4}$, 280
iterations at $\delta=10^{-4}$. What is crucial is that the same geometric
patterns, with the same energy density and very nearly the same amplitudes,
are found in all three cases. There is no qualitative and very little
quantitative difference.

(a) stripe

(b) domain

(c) checkerboard

Figure 12: Varieties of long-range order in momentum space as displayed by the
momentum space spin correlator $S(k)$, defined in eq. (20), on a $24^{2}$
lattice at $U=3$ . (a) stripe pattern, $f=0.8$. (b) domain pattern, also at
$f=0.8$. (c) checkerboard pattern, $f=1$.

Figure 13: A plot of $S_{4}/S_{tot}$ (purple squares, left hand $y$-axis), and
$S_{4}$ alone (filled circles, right hand $y$-axis) vs. $U$ at density
$f=0.8$. Here $S_{tot}$ is the sum of all $S(k)$ values, and $S_{4}$ is the
sum of the four largest $k$ values. Long range order seems most evident near
$U=3$, where the sum is almost saturated by the top four values.

A more quantitative measure of long range order in magnetization is the
momentum space correlator $S(k_{x},k_{y})$ defined in (20). In Fig. 12 we
display this quantity for the stripe, rectangular domain, and checkerboard
patterns, corresponding to Figs. 10(b), 10(c), 10(d), respectively, again at
$U=3$ on a $24^{2}$ lattice. Concentration of $S(k_{x},k_{y})$ at just a few
values (even just one $k_{x}=k_{y}=\pi$ for the checkerboard pattern) is
obviously evidence of long range order, and it is of interest to study the
$U$-dependence. So let us choose $f=0.8$, where we see both stripes and
rectangular domains, and let $S_{4}$ be the sum of the largest four values of
$S(k_{x},k_{y})$, with $S_{tot}$ the sum of all values. In Fig. 13 we plot the
ratio $S_{4}/S_{tot}$ (scale on the left-hand y-axis), and $S_{4}$ alone
(scale on the right-hand y-axis). The long range order at $U=3$ decays away on
either side. In addition, for $U<2$, $S_{4}$ is hardly distinguishable from
zero on the scale of the plot.

###  4.5 Pair correlation and the checkerboard

It is obvious that in an eigenstate of particle number, ${\langle
c(k,\uparrow)c(-k,\downarrow)\rangle=0}$. On the other hand, the correlator of
a pair creation operator in the vicinity of point $x$, and a pair destruction
operator in the vicinity of point $y$, where $R=|x-y|\gg 1$, could be non-
zero. Transforming to momentum space, as in eq. (25), it is of special
interest to see where in $k$-space this operator is non-zero, and whether
there is any indication of d-wave symmetry. In the ground state we find

$\displaystyle\Delta(k^{\prime},k)$ $\displaystyle=$
$\displaystyle\langle\Phi|c^{\dagger}(k^{\prime},\uparrow)c^{\dagger}(-k^{\prime},\downarrow)c(k,\uparrow)c(-k,\downarrow)|\Phi\rangle$
(27) $\displaystyle=$
$\displaystyle\sum_{j=2}^{M}\sum_{i=1}^{j-1}\\{\phi_{i}(-k^{\prime},\uparrow)\phi_{j}(k^{\prime},\downarrow)-\phi_{i}(k^{\prime},\downarrow)\phi_{j}(-k^{\prime},\uparrow)\\}$
$\displaystyle\times\\{\phi_{i}(k,\uparrow)\phi_{j}(-k,\downarrow)-\phi_{i}(-k,\downarrow)\phi_{j}(k,\uparrow)\\}\
.$

Define

$\displaystyle\omega_{ij}(k)$ $\displaystyle=$
$\displaystyle\phi_{i}(k,\uparrow)\phi_{j}(-k,\downarrow)-\phi_{i}(-k,\downarrow)\phi_{j}(k,\uparrow)$
(28) $\displaystyle=$
$\displaystyle\phi_{i}(k,\uparrow)\phi^{*}_{j}(k,\downarrow)-\phi^{*}_{i}(k,\downarrow)\phi_{j}(k,\uparrow)\
.$

Then

$\Delta(k^{\prime},k)=\sum_{j=2}^{M}\sum_{i=1}^{j-1}\omega^{*}_{ij}(k^{\prime})\omega_{ij}(k)\
.$ (29)

The space of all $k,k^{\prime}$ is four-dimensional, and we choose a two
dimensional slice by taking $k^{\prime}$ to be the wavevector $k$ with $x,y$
components interchanged, i.e $k=(k_{x},k_{y}),k^{\prime}=(k_{y},k_{x})$. If
pairing follows a D-wave pattern, then in this two-dimensional slice we would
expect something like

$\Delta(k^{\prime},k)\sim-(\cos(k_{x})-\cos(k_{y}))^{2}\ ,$ (30)

which is negative everywhere, and most negative at $k_{x}=0,k_{y}=\pm\pi$ and
$k_{x}=\pm\pi,k_{y}=0$.

These features are seen, roughly, but only for checkerboard patterns with the
momentum space distribution shown in Fig. 14(a), in the immediate neighborhood
of $f=1$. The correlation is also only significant at the edge of the $n(k)$
occupation zone. An example at $f=1,U/t=1$ is shown in Fig. 14(a), but the
very same correlation appears at half-filling over a range of two orders of
magnitude, from $U/t=0.1$ to $U/t=10$, and the corresponding plots are very
similar to Fig. 14(a).222 We are only interested in the correlator for
$k^{\prime}\neq k$, and have (arbitrarily) set $\Delta(k,k)=0$ in the figures
shown.

The pair correlation at half-filling only starts to disintegrate below
$U/t=0.03$. The situation is quite different away from half-filling, e.g. the
checkerboard order and similar pair correlation at $f=0.95$ is only seen at a
rather strong repulsion of $U/t=10$ (Fig. 14(b)). The picture seen in Figs.
14(a) and 14(b) is not exactly what we have in eq. (30); for one thing the
correlation is sharply concentrated at the boundary of the $n(k)$ occupation
zone, and the numerical values do not really agree with (30). But the fact
that the correlator is everywhere negative, and most negative at $(0,\pm\pi)$
and $(\pm\pi,0)$ is definitely reminiscent of d-wave correlation. We
emphasize, however, that the checkerboard pattern and this type of pair
correlation are only found in the immediate neighborhood of half-filling, and
even then (e.g. at $f=0.95$) may only be seen at comparatively large values of
$U/t$, as in Fig. 14(b). When the checkerboard pattern is absent the
correlation function $\Delta(k,k^{\prime})$ is still generally concentrated at
the edge of the $n(k)$ distribution, but may be everywhere positive, or a
mixture of positive and negative values. Or the pair correlation may be
negligible, as seen in Fig. 14(c) (note the scale) at $f=0.8,U=3$, where there
is a stripe order.333We note that a weak-coupling analysis of the 2D Hubbard
model in ref. [26] also found d-wave pairing strongest near half-filling.

(a) f=1, U=1

(b) f=0.95, U=10

(c) f=0.8, U=3

Figure 14: Display of the pair correlation function $\Delta(k^{\prime},k)$ in
the $k_{x},k_{y}$ plane, where $k^{\prime}=(k_{x},k_{y}),k=(k_{y},k_{x})$. (a)
The pattern shown, at half-filling and $U=1$, has some of the features of
d-wave correlation. This pattern is virtually unchanged at half-filling across
two orders of magnitude, from $U=0.1$ to $10$. (b) Pair correlation at
$f=0.95$ and $U=10$. Even slightly away from half-filling, the pattern seen
here at $U=10$ is not seen at moderate and small values of $U$. (c) The
correlator is negligible (note the scale), compared to the previous
subfigures, at $f=0.8,U=3$. The pattern seen in (a) and (b) is always
associated with a checkerboard pattern, and in case (c) there is instead a
stripe pattern.

###  4.6 Localization and particle-hole symmetry

We return, finally, to the question which motivated this work. The stripe,
domain and checkerboard patterns do not imply that one electron wave functions
are localized, and in fact the opposite is true, as we have already noted. The
Hartree-Fock state may break translation invariance, as is obvious from the
geometric patterns, but according to our investigation of the IPR values of
the constituent one-particle wave functions in $|\Phi\rangle$, these wave
functions spread over most of the lattice. The situation is different at
strong $U/t$, as can be seen in Fig. 15, at densities $f=0.9,1.0,1.1$ and
$U/t=30,100$. In this figure we display the IPR values of all $2L^{2}$
particle states ($L=24$), both filled and unfilled. An interesting feature is
that it is the hole (unfilled) states which are localized at $f<1$, and the
electron (filled) states which are localized at $f>1$, with about an equal
degree of localization among hole and particle states at $f=1$. The dramatic
shift from localized holes at $f<1$ to localized electrons at $f>1$, no doubt
a consequence of particle-hole symmetry, is most obvious at an extremely
strong coupling $U/t=100$, but it also quite evident at, e.g., $U/t=30$. As
$U/t$ is further reduced to moderate values, the asymmetry in particle/hole
localization on either side of $f=1$, and the localization itself, gradually
disappears.

It is a general feature of the iterative procedure that the number of
iterations required for convergence increases with $U/t$, and it is of
interest to repeat the calculation with more stringent convergence criteria.
The computation of IPR values shown in Fig. 15(d), with convergence parameter
$\delta=10^{-3}$ was repeated at $\delta=10^{-4}$ and $\delta=10^{-5}$. At
$U/t=30,f=0.9$, the number of iterations required for convergence increases
very significantly, with 2030 iterations required at $\delta=10^{-3}$, 10810
iterations at $\delta=10^{-4}$, and 40510 at $\delta=10^{-5}$. Nevertheless,
the qualitative picture of localization is unchanged by making the convergence
criterion more stringent. In Fig. 16 we compare IPR values at $U=30,f=0.9$ for
different values of the convergence parameter $\delta$, where each subfigure
is derived from a different self-consistent solution of the Hartree-Fock
equations. Here it should be understood that there is always some variation in
the IPR values of one-particle eigenstates $\phi_{i}$ from one self-consistent
solution to the next, and the variation seen by decreasing $\delta$ is not
much different from the variation among different self-consistent solutions at
fixed $\delta$.

Figure 15: Localization of particle and hole states at very strong repulsion
$U=30,100$ as quantified by their IPR values. The $x$-axis is the one-particle
state number, i.e. the $i$ of $\phi_{i}(x,s)$. Note that it is hole states
which are localized for $f<1$, and particle states which are localized for
$f>1$, with a symmetric distribution at $f=1$. Subfigures (a,b,c): $U=100$ and
$f=0.9,1.0,1.1$ respectively. Subfigures (d,e,f): $U=30$ and $f=0.9,1.0,1.1$
respectively.

(a) $\delta=10^{-3}$

(b) $\delta=10^{-4}$

(c) $\delta=10^{-5}$

Figure 16: A comparison of IPR values at $U=30,f=0.9$ (same as Fig. 15(d)),
but with different convergence criteria parametrized by $\delta$. Each figure
is taken from a different self-consistent solution. While there is some
variation in the IPR values shown in these figures, it is consistent with the
variation at fixed $\delta=10^{-3}$ among different self-consistent solutions.

Anderson localization is a well known phenomenon for particle propagation in a
random potential, and we feel it is significant that here, by contrast, there
is localization in a translationally invariant system, where the degree of
localization among the different states depends on the random initialization.
This seems to answer in the affirmative the question raised in [1], which asks
whether this phenomenon is possible.

##  5 Conclusions

The Hartree-Fock approach yields a multiplicity of self-consistent solutions
of surprising complexity; more complexity than one might expect given that any
Hartree-Fock state is just a single Slater determinant, where entanglement is
limited to what is required by Fermi statistics. We see emergent spin patterns
in the form of stripes, checkerboards, and rectangular domains, band gaps
correlated with antiferromagnetic order, localization of particles or holes
for strong repulsion at $f>1$, $f<1$ respectively, and pair correlation
functions near half-filling which have properties similar to d-wave
distributions. Perhaps most important is the multiplicity itself; i.e. the
enormous number self-consistent solutions of the Hartree-Fock equations,
derived from an iterative procedure with only slightly different random
initializations, that are very nearly degenerate in energy.

One can take different views of this near degeneracy. One view is that the
true ground state, so far as it can be determined in the Hartree-Fock
approach, is the state with the lowest energy. This is in a landscape of
states which have very nearly the same energy, but are distinguished from one
another by other observables, such as local magnetization and energy gap,
which can differ widely among different solutions. If it were true that the
low temperature behavior is dominated by a unique lowest energy ground state,
then in our opinion this would rule out translation non-invariance in the
expectation value of any observable. The other possible view, if we take the
multiplicity of near degenerate states seriously, is that as the system is
cooled from finite temperature it finds itself in one of those near-degenerate
states, rarely if ever falling into the true quantum ground state. If that is
so, then the 2D Hubbard model may have other features analogous to spin
glasses, e.g. some degree of non-ergodicity, which deserve further study.

We have also found localization of one-particle hole and electron states just
below and just above half-filling, at strong repulsion, with IPR values which
again depend on the particular solution obtained from the particular
initialization. Unlike in ordinary Anderson localization there is no random
potential intrinsic to the many body Hamiltonian, but a random effective
potential, as seen by a single particle, may arise from the random
initialization. There has been speculation [1] whether small random
differences in intialization might drive a system to become localized, and our
results would seem to support that possibility.

The Hartree-Fock approach to the Hubbard model must obviously be viewed with
caution; as in any mean-field theory it simply ignores the correlation of the
average field acting on one particle with the position/spin of that particle,
and this neglect is well known to lead to errors, particularly in the
neighborhood of a transition. In its defense, we have seen that apart from
very strong repulsion, the one-electron wave functions are unlocalized. That
means that a single electron is effectively interacting with all other
electrons in the system, and not just with a handful of nearest neighbors. The
fact that one degree of freedom is interacting with many is the usual mean-
field (and in particular Hartree-Fock) justification for replacement of the
“many” by their average. On the other hand, since entanglement in a Hartree-
Fock state is limited to what is required by Fermi statistics, any phenomenon
which depends on entanglement beyond that requirement will simply be invisible
in the framework of the Hartree-Fock approximation. Among the important
observables which may depend crucially on entanglement we would include
electron pairing.

All of the work presented in this article concerns the zero temperature
Hartree-Fock approximation to the 2D Hubbard model. We leave an investigation
of localization by more sophisticated methods, perhaps including the effects
of finite temperature, to later investigation.

Acknowledgements
This work is supported by the U.S. Department of Energy under Grant No. DE-
SC0013682.

## References

  * [1] R. Nandkishore, D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6 (1) (2015) 15–38.  doi:10.1146/annurev-conmatphys-031214-014726.
  * [2] D. R. Penn, Stability theory of the magnetic phases for a simple model of the transition metals, Phys. Rev. 142 (1966) 350–365.  doi:10.1103/PhysRev.142.350.
URL https://link.aps.org/doi/10.1103/PhysRev.142.350

  * [3] J. E. Hirsch, Two-dimensional hubbard model: Numerical simulation study, Phys. Rev. B 31 (1985) 4403–4419.  doi:10.1103/PhysRevB.31.4403.
URL https://link.aps.org/doi/10.1103/PhysRevB.31.4403

  * [4] D. Poilblanc, T. M. Rice, Charged solitons in the hartree-fock approximation to the large-u hubbard model, Phys. Rev. B 39 (1989) 9749–9752.  doi:10.1103/PhysRevB.39.9749.
URL https://link.aps.org/doi/10.1103/PhysRevB.39.9749

  * [5] J. Zaanen, O. Gunnarsson, Charged magnetic domain lines and the magnetism of high-${T}_{c}$ oxides, Phys. Rev. B 40 (1989) 7391–7394.  doi:10.1103/PhysRevB.40.7391.
URL https://link.aps.org/doi/10.1103/PhysRevB.40.7391

  * [6] K. Machida, Magnetism in la2cuo4 based compounds, Physica C: Superconductivity 158 (1-2) (1989) 192–196.
  * [7] H. Schulz, Domain walls in a doped antiferromagnet, Journal de Physique 50 (18) (1989) 2833–2849.
  * [8] H. J. Schulz, Incommensurate antiferromagnetism in the two-dimensional hubbard model, Phys. Rev. Lett. 64 (1990) 1445–1448.  doi:10.1103/PhysRevLett.64.1445.
URL https://link.aps.org/doi/10.1103/PhysRevLett.64.1445

  * [9] M. Ichimura, M. Fujita, K. Nakao, Spin density wave states in two dimensional hubbard model, Journal of the Physical Society of Japan 61 (6) (1992) 2027–2039.  arXiv:https://doi.org/10.1143/JPSJ.61.2027, doi:10.1143/JPSJ.61.2027.
URL https://doi.org/10.1143/JPSJ.61.2027

  * [10] J. Vergés, E. Louis, P. Lomdahl, F. Guinea, A. Bishop, Holes and magnetic textures in the two-dimensional hubbard model, Physical Review B 43 (7) (1991) 6099.
  * [11] M. Inui, P. Littlewood, Hartree-fock study of the magnetism in the single-band hubbard model, Physical Review B 44 (9) (1991) 4415.
  * [12] C. Dasgupta, J. Halley, Phase diagram of the two-dimensional disordered hubbard model in the hartree-fock approximation, Physical Review B 47 (2) (1993) 1126\.
  * [13] V. Bach, E. H. Lieb, J. P. Solovej, Generalized hartree-fock theory and the hubbard model, Journal of statistical physics 76 (1) (1994) 3–89.
  * [14] M. Imada, A. Fujimori, Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70 (1998) 1039–1263.  doi:10.1103/RevModPhys.70.1039.
URL https://link.aps.org/doi/10.1103/RevModPhys.70.1039

  * [15] P. Fazekas, Lecture Notes on Electron Correlation and Magnetism, WORLD SCIENTIFIC, 1999.  arXiv:https://www.worldscientific.com/doi/pdf/10.1142/2945, doi:10.1142/2945.
URL https://www.worldscientific.com/doi/abs/10.1142/2945

  * [16] J. Xu, C.-C. Chang, E. J. Walter, S. Zhang, Spin-and charge-density waves in the hartree–fock ground state of the two-dimensional hubbard model, Journal of Physics: Condensed Matter 23 (50) (2011) 505601.
  * [17] R. Rodríguez-Guzmán, C. A. Jiménez-Hoyos, R. Schutski, G. E. Scuseria, Multireference symmetry-projected variational approaches for ground and excited states of the one-dimensional hubbard model, Phys. Rev. B 87 (2013) 235129.  doi:10.1103/PhysRevB.87.235129.
URL https://link.aps.org/doi/10.1103/PhysRevB.87.235129

  * [18] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, E. Gull, Solutions of the two-dimensional hubbard model: Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5 (2015) 041041.  doi:10.1103/PhysRevX.5.041041.
URL https://link.aps.org/doi/10.1103/PhysRevX.5.041041

  * [19] E. Pavarini, E. Koch, A. Lichtenstein, D. Vollhardt (Eds.), The LDA+DMFT approach to strongly correlated materials, Vol. 1 of Schriften des Forschungszentrums Julich. Reihe modeling and simulation, Forschungszenrum Julich GmbH Zentralbibliothek, Verlag, Julich, 2011.
URL https://juser.fz-juelich.de/record/17645

  * [20] C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell, R. T. Scalettar, Quantum monte carlo study of the two-dimensional fermion hubbard model, Phys. Rev. B 80 (2009) 075116.  doi:10.1103/PhysRevB.80.075116.
URL https://link.aps.org/doi/10.1103/PhysRevB.80.075116

  * [21] R. Blankenbecler, D. J. Scalapino, R. L. Sugar, Monte carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24 (1981) 2278–2286.  doi:10.1103/PhysRevD.24.2278.
URL https://link.aps.org/doi/10.1103/PhysRevD.24.2278

  * [22] J. Kertesz, Existence of weak singularities when going around the liquid-gas critical point, Physica A161 (1989) 58–62.
  * [23] S. A. Kivelson, I. P. Bindloss, E. Fradkin, V. Oganesyan, J. M. Tranquada, A. Kapitulnik, C. Howald, How to detect fluctuating stripes in the high-temperature superconductors, Rev. Mod. Phys. 75 (2003) 1201–1241.  doi:10.1103/RevModPhys.75.1201.
URL https://link.aps.org/doi/10.1103/RevModPhys.75.1201

  * [24] J. M. Tranquada, Topological doping and superconductivity in cuprates: An experimental perspective, Symmetry 13 (12) (2021) 2365.
  * [25] B. V. Fine, Hypothesis of two-dimensional stripe arrangement and its implications for the superconductivity in high-${T}_{c}$ cuprates, Phys. Rev. B 70 (2004) 224508.  doi:10.1103/PhysRevB.70.224508.
URL https://link.aps.org/doi/10.1103/PhysRevB.70.224508

  * [26] S. Raghu, S. A. Kivelson, D. J. Scalapino, Superconductivity in the repulsive hubbard model: An asymptotically exact weak-coupling solution, Phys. Rev. B 81 (2010) 224505.  doi:10.1103/PhysRevB.81.224505.
URL https://link.aps.org/doi/10.1103/PhysRevB.81.224505