# Phase diagram for the extended Hubbard model on a square lattice

Sebastião dos Anjos Sousa-Júnior  Instituto de Física, Universidade Federal do
Rio de Janeiro Cx.P. 68.528, 21941-972 Rio de Janeiro RJ, Brazil  Natanael C.
Costa  Instituto de Física, Universidade Federal do Rio de Janeiro Cx.P.
68.528, 21941-972 Rio de Janeiro RJ, Brazil  Raimundo R. dos Santos  Instituto
de Física, Universidade Federal do Rio de Janeiro Cx.P. 68.528, 21941-972 Rio
de Janeiro RJ, Brazil

###### Abstract

The Extended Hubbard model describes electron-electron interactions and is
commonly used to understand various phenomena in strongly correlated systems.
In this paper, we focus on investigating the phase diagram of this model,
which considers a nearest neighbor interaction, $V$, in addition to the on-
site term $U$. To perform our study, we utilize the determinant quantum Monte
Carlo method simulations, where the average sign of the determinant is used to
extract information about the quantum critical points. When both $U$ and $V$
are positive, the on-site term favors antiferromagnetism (AFM), while the
extended interaction favors the formation of a charge density wave (CDW). In
cases where at least one of the interactions is negative, superconductivity
competes with the above-mentioned phases, and the system displays phase
separation for intense attractive extended interaction. Our findings indicate
that only $s$-wave and $d$-wave pairing symmetries exist in such scenarios.

Soon after the discovery of high-temperature cuprate superconductors Bednorz
and Müller (1986); Varma (2020); Zhou  _et al._ (2021), a widespread consensus
was formed around the idea that the basic physical mechanism leading to
superconductivity was contained in the two-dimensional Hubbard model: pairing
would emerge from strong antiferromagnetic (AFM) fluctuations arising from the
competition between itinerancy (hopping) and localisation, the latter driven
by an on-site repulsion of strength $U$ Anderson (1987); Emery (1987). Despite
its simplicity, the repulsive Hubbard model has eluded an unambiguous
characterisation of superconductivity. Subsequently, it was proposed Scalapino
(2012) that antiferromagnetic fluctuations could actually lead to an
additional effective nearest-neighbour attractive interaction, $V$. The ground
state phase diagram for the extended Hubbard model (EHM) in one dimension and
at half filling has been known for some time; see, e.g. Ref. Lin  _et al._
(1995). For different values and signs of $U$ and $V$, one finds phases such
as charge-density wave (CDW), spin-density wave (SDW), and $s$-wave
superconductivity. A phase separated (PS) state appears for sufficiently large
$-V$, e.g., half of the lattice with doubly occupied sites, and the other half
empty, in the strong $-V$ regime. Less conventional phases such as bond-
ordered wave Lin  _et al._ (1995) and $p$-wave superconducting Xiao  _et al._
(2022) have also been proposed to fit into the diagram.

Figure 1: Ground state phase diagram for the extended Hubbard model at half
filling, obtained through our DQMC simulations; $U$ and $V$ are on-site and NN
couplings, respectively. Filled symbols are estimates for critical points
obtained through the joint analyses of structure factors and
$\langle\text{s}\rangle$; the latter played no part in determining the
critical points represented by empty symbols. CDW stands for charge density
wave, AFM for antiferromagnetic, $s$-wave and $d$-wave are superconducting
phases, and PS denotes phase separation. Lines through data points are guides
to the eye.

In two dimensions, while one expects the ground state phase diagram at half
filling to share similarities with that for the one-dimensional model, the
picture is far from settled. Indeed, different methods agree with the
existence of a CDW-AFM transition near $U=4V$ Zhang and Callaway (1989); Huang
_et al._ (2013); Vandelli  _et al._ (2020); Yao  _et al._ (2022). Nonetheless,
the square lattice topology in principle allows for a wider range of pairing
symmetries to be stabilised in the $V<0$ region, which has proved very hard to
probe theoretically. Weak-coupling methods Huang  _et al._ (2013); Wolf  _et
al._ (2018) may not fully capture the competition between the tendencies of
phase separating and forming pairs of different symmetries, such as $s$ and
$d_{x^{2}-y^{2}}$. Numerically exact diagonalisation of the Hamiltonian is
currently restricted to $4\times 4$ systems, so that the prediction of a
$p$-wave paired state Chen  _et al._ (2022) may not stand for larger systems.
Quantum Monte Carlo (QMC) simulations Blankenbecler  _et al._ (1981); Hirsch
(1985), on the other hand, suffer from the infamous ‘minus-sign problem’
Hirsch (1985); Scalettar  _et al._ (1989); Kawashima (2002); dos Santos
(2003); Becca and Sorella (2017): when the effective Boltzmann factor (given
by a product of fermionic determinants) becomes negative, averages are taken
with it in absolute value, at the expense of dividing by its own average,
$\langle\text{s}\rangle$, thus introducing excessive noise when
$\langle\text{s}\rangle\ll 1$. Recent studies of the EHM through QMC
simulations were restricted to parameter regimes where the minus-sign problem
was not too severe, so that the focus was mostly on the CDW-AFM transition
Sushchyev and Wessel (2022). The inescapable conclusion is that so far the
overall knowledge of the EHM phase diagram is, at best, semi-qualitative.

In light of the recent suggestion Mondaini  _et al._ (2022a, b) that a
severely degraded $\langle\text{s}\rangle$ may actually be used to pinpoint
phase transition points and boundaries, a thorough investigation of the EHM
through determinant QMC (DQMC) simulations is certainly in order. With this in
mind, here we report on the results for the ground state phase diagram at
half-filling in the whole $V$-$U$ parameter space with an unprecedented
accuracy; our main findings are highlighted in Fig. 1.

Figure 2: a) Contour plot of the average sign for $U=4$ on a $4\times 4$
lattice at half filling. b) Contour plot of the average sign for $T=0.1$ in
the $V$-$U$ plane for a $10\times 10$ lattice.

The Hamiltonian for the EHM reads,

$\displaystyle\mathcal{H}=$
$\displaystyle-t\sum_{\langle\textbf{i,j}\rangle}(c_{\textbf{i}\sigma}^{\dagger}c^{\phantom{\dagger}}_{\textbf{j}\sigma}+\text{H.c.})-\mu\sum_{\textbf{i},\sigma}n_{\textbf{i},\sigma}$
$\displaystyle+U\sum_{\textbf{i}}n_{\textbf{i}\uparrow}n_{\textbf{i}\downarrow}+V\sum_{\langle\textbf{i,j}\rangle}n_{\textbf{i}}n_{\textbf{j}},$
(1)

where $\mathbf{i}$ and $\mathbf{j}$ denote sites of a square lattice, with
$\langle\mathbf{i},\mathbf{j}\rangle$ restricting the sums to nearest
neighbour (NN) sites. In standard second-quantised notation, the first term
describes fermionic hopping (energy scale $t$), the second controls the band
filling through the chemical potential, $\mu$, while the third and fourth
terms describe the on-site and NN interactions, with strengths $U$ and $V$,
respectively. Hereafter, the chemical potential is set to $\mu=U/2+4V$ to
yield a half-filled band (by virtue of particle-hole symmetry), and energies
is expressed in units of $t$.

In the DQMC method Hirsch (1985); Scalettar  _et al._ (1989); Kawashima
(2002); dos Santos (2003); Becca and Sorella (2017), discrete Hubbard-
Stratonovich (HS) transformations Hirsch (1983) are used to express the
quartic interactions in quadratic forms. This leads to 5 auxiliary fields on
any given site $\mathbf{i}$: 1 from the on-site term, and 4 from the bonds
emanating from it Zhang and Callaway (1989); see, e.g., the Supplemental
Materials. The non-commutation between the one-body and the two-body terms of
the Hamiltonian is taken care of through a Suzuki-Trotter decomposition, which
adds an imaginary-time dimension, $L_{\tau}=\beta/\Delta\tau$, with $\beta$
being the inverse temperature, and $\Delta\tau$ the discrete time step. The
trace over the fermionic degrees of freedom can then be performed, leading to
the product of determinats alluded to before, which weigh the configurations
of the HS Ising fields sampled by importance, as in usual Monte Carlo methods
Hirsch (1985); Scalettar  _et al._ (1989); Kawashima (2002); dos Santos
(2003); Becca and Sorella (2017).

Figure 3: (a) Spin structure factor (filled symbols, left vertical scale) and
scaled charge structure factor (empty symbols, right vertical scale), and (b)
average sign as functions of $V$. All data are for fixed $U=2$ and $\beta=4$,
and different linear system sizes, $L$. The grey bar highlights the region
where $\langle\text{s}\rangle$ dips.

The type of ordering is characterised by quantities such as the (staggered)
charge structure factor,

$\displaystyle
S_{\text{cdw}}=\frac{1}{N}\sum_{\textbf{i},\textbf{j}}(-1)^{|\mathbf{i}-\mathbf{j}|}\langle(n_{\textbf{i},\uparrow}+n_{\textbf{i},\downarrow})(n_{\textbf{j},\uparrow}+n_{\textbf{j},\downarrow})\rangle,$
(2)

the antiferromagnetic structure factor,

$\displaystyle
S_{\text{afm}}=\frac{1}{N}\sum_{\textbf{i},\textbf{j}}(-1)^{|\mathbf{i}-\mathbf{j}|}\langle
S^{z}_{\mathbf{i}}S^{z}_{\mathbf{j}}\rangle,$ (3)

and the pairing structure factor,

$\displaystyle
P_{\text{sc}}(\alpha)=\frac{1}{N}\sum_{\textbf{i,j}}\langle\Delta_{\alpha}(\textbf{i})\Delta^{\dagger}_{\alpha}(\textbf{j})\rangle,$
(4)

with

$\Delta_{\alpha}(\textbf{i})=\sum_{\textbf{a}}f_{\alpha}(\textbf{a})c_{\textbf{i}\downarrow}c_{\textbf{i+a}\uparrow}.$
(5)

In the above equations, $N=L\times L$ is the number of sites for a linear size
$L$, $(-1)^{|\mathbf{i}-\mathbf{j}|}$ is $\pm 1$ if $\mathbf{i}$ and
$\mathbf{j}$ are on the same or opposite sublattices,
$S^{z}_{\mathbf{i}}=(n_{\mathbf{i}\uparrow}-n_{\mathbf{i}\downarrow})$ is the
$z$-component of the spin operator, and $f_{\alpha}(\textbf{a})$ is the form
factor for a given pair-wave symmetry, $\alpha={s,d,p}$ White  _et al._
(1989). In some circumstances, it is more appropriate to calculate the pair
susceptibility

$\displaystyle\chi^{\alpha}_{\rm
sc}(\beta)=\frac{1}{N}\sum_{\textbf{i,j}}\int^{\beta}_{0}\mathrm{d}\tau\,\langle\Delta_{\mathbf{i},\alpha}^{\phantom{\dagger}}(\tau)\Delta^{\dagger}_{\mathbf{j},\alpha}(0)\rangle~{},$
(6)

with
$\Delta_{\alpha}(\mathbf{i},\tau)=\sum_{\mathbf{a}}f_{\alpha}(\mathbf{a})c^{\phantom{\dagger}}_{\mathbf{i}\downarrow}(\tau)c^{\phantom{\dagger}}_{\mathbf{i}+\mathbf{a}\uparrow}(\tau)$,
and
$c_{\mathbf{i}\sigma}(\tau)=e^{\tau\cal{H}}c_{\mathbf{i}\sigma}e^{-\tau\cal{H}}$,
which provides a stronger signal of pairing properties. As mentioned before,
$\langle\text{s}\rangle$, which is automatically calculated in the
simulations, will also play a crucial role in our analyses.

When $V\neq 0$ one cannot guarantee that $\langle\text{s}\rangle=1$ even at
half filling. As shown in Fig. 2 (a), the average sign decreases as the
temperature is lowered, and it is seen to be particularly worse near $V=1$.
Recalling that the phase boundary for the SDW-CDW transition is $V_{c}\approx
U/4$, the minimum of $\langle s\rangle$ at $V=1$ can hardly be regarded as
fortuitous. A mapping of $\langle\text{s}\rangle$ in the $V$-$U$ plane at a
fixed low temperature is shown in Fig. 2 (b). Apart from the second quadrant
($U<0$,$V>0$), one sees that there are regions where
$\langle\text{s}\rangle\ll 1$, surrounded by less severe ones. In what follows
we combine analyses of $\langle\text{s}\rangle$ with the quantities defined by
Eqs. (2)-(4) to obtain the phase boundaries in Fig. 1.

We start with the transition between the AFM (Mott) phase and the CDW phase.
As indicated in Fig. 3 (a), increasing $V$ causes a sharp decrease of
$S_{\text{afm}}$ and a sharp increase of $S_{\text{cdw}}$ thus indicating a
phase transition at $V_{c}=0.6\pm 0.05$, where the error reflects the $V$
increments in Fig. 3 (a). This is accompanied by a dip in
$\langle\text{s}\rangle$ at the same value of $V$; see Fig. 3 (b). The
behaviour for $U=1$ is similar, leading to $V_{c}=0.2\pm 0.05$, so that,
similarly to the 1D case Hirsch (1984), the AFM-CDW critical line lies
slightly above the strong-coupling prediction, $V_{c}=U/4$.

Figure 4:  (a) $s$-pairing structure factor (full lines, left vertical axis)
and charge structure factor (dashed line, right vertical axis), and (c)
average sign as functions of $V$ at fixed $U$, for different temperatures; (b)
charge structure factor [notice the different scale from (a)] and (d) average
sign as functions of $V$, at fixed inverse temperature, for different values
of the on-site attraction.

Let us now focus on the second ($U<0$, $V>0$) and third ($U<0$, $V<0$)
quadrants of the parameter space. Figure 4 (a) shows that the charge structure
factor decreases steadily as one crosses the $V=0$ line towards $V<0$, at
fixed $U=-2$; by contrast, the $s$-wave pairing structure factor increases
steadily with increasing $\beta$ on the $V<0$ side. This behaviour is
consistent with the fact that exactly at $V=0$ one reaches the attractive
Hubbard model, which definitely displays a CDW state coexisting with an
$s$-wave superconducting state Micnas  _et al._ (1990); Fontenele  _et al._
(2022): one deals effectively with a three-component order parameter, one for
CDW and two for the superconductivity. Further decrease in $V$ causes a
suppression of $s$-wave pairing correlations, as one enters the PS region, at
$V=-0.45\pm 0.05$ for $U=-2$. Interestingly, Figure 4 (c) shows that as the
temperature decreases, the dip in $\langle\text{s}\rangle$ deepens and widens
throughout the superconducting phase; in the PS region,
$\langle\text{s}\rangle$ returns to 1. Thus the sharp drops in
$\langle\text{s}\rangle$ indicate the boundaries of the superconducting phase
with the CDW and PS regions. Figure 4 (b) provides an interesting insight into
the behaviour of the charge structure factor, by keeping the temperature fixed
at $\beta=6$, and examining how the plots evolve with $U$. Two regimes can be
distinguished: one in which charge correlations are enhanced as $|U|$
increases, and another in which they decrease as $|U|$ increases. The analysis
of Fig. 4 (a) hence suggests that this change indicates the entrance into the
PS region. This is again in accordance with the behaviour of
$\langle\text{s}\rangle$ at fixed $\beta$ and for different values of $U$,
depicted in Fig. 4 (d): the dip occurs at smaller values of $V$ as $U$
increases. These findings are summarised in Fig. 1.

Figure 5: Density probability distribution for several values of $V$ with
$U=-0.5$ and $\beta=6.0$. As the intensity of the attractive extended
interaction increases the shape of the histograms undergoes from a single
peaked to a double peaked distribution.

The transition to the PS state deserves a complementary look, by examining the
density distribution shown in Fig. 5, generated by collecting the values of
$n_{\mathbf{i}}$ over the DQMC runs. For $V=0$, the distributions are
represented by singly peaked histograms (not shown), centred at $n=1$. As $V$
decreases, the distributions first broaden, still with a peak at $n=1$, but a
change to doubly-peaked at $n=0$ and $n=2$ takes place, interpreted as a
signature of a phase separated state. In Fig. 5, this occurs at
$V\approx-0.6$, which also marks the dip in $\langle\text{s}\rangle$ for
$U=-0.5$; see Fig. 4 (d).

Returning to the SC state, we must settle the issue of pairing symmetries. To
this end, we start with the difference between the effective pairing
susceptibilities with $s$ and $d_{x^{2}-y^{2}}$ symmetries, which presumably
dominate the third quadrant. At this point, we investigate the effective pair
(vertex) susceptibility, i.e., $\chi^{\rm eff}_{\rm sc}(\alpha)=\chi_{\rm
sc}(\alpha)-\bar{\chi}_{\rm sc}(\alpha)$, with $\bar{\chi}_{\alpha{\rm(sc)}}$
being the noninteracting susceptibility White  _et al._ (1989). A positive
(negative) response of $\chi^{\rm eff}_{\rm sc}(\alpha)$ signals an
enhancement (weakening) of pair correlations for the corresponding
$\alpha$-wave symmetry. From Fig. 6 we see that the region above the line
$V=U/2$ is dominated by $s$-wave pairing, while the region below the line is
dominated by $d_{x^{2}-y^{2}}$-wave pairing. Interestingly,
$\langle\text{s}\rangle$ shows no dip at this transition; we may attribute
this to the lack of change in the number of components of the order parameter
on either side of the transition.

Figure 6: a) Contour plot of the difference between $s$-wave and $d$-wave
effective pairing susceptibilities for $\beta=5.0$ and $L=8$.

Moving on to the fourth quadrant, $V<0$ and $U>0$, Fig. 7(a) follows how the
spin structure factor changes as $V$ is increased. Below $V\approx-0.7$,
$S_{\text{afm}}$ is quite insensitive to the temperature. For $-0.7\lesssim
V\lesssim-0.25$, $S_{\text{afm}}$ actually decreases as the temperature
increases, while the superconducting structure factors are enhanced in this
interval [see Fig. 7(b) and (d)], with $P_{\text{sc}}(d_{x^{2}-y^{2}})$
tending to dominate over $P_{\text{sc}}(p_{x})$. The dip in
$\langle\text{s}\rangle$ at $V\approx-0.7$ provides additional support to the
interpretation of a PS–$d$-wave transition at this point in the diagram.
Beyond $V\approx-0.25$, $S_{\text{afm}}$ increases as the temperature is
lowered, which is also accompanied by a dip in $\langle\text{s}\rangle$: this
signals a $d$-wave–AFM transition. By repeating these analyses for other
values of $U>0$, we obtain the critical curve in Fig. 1, in which the error
stems from the resolution of the crossings in Fig. 7(a).

In view of the recent suggestion that a $p$-wave SC state could be stabilised,
we examined pairing structure factors and susceptibilities. Typical data are
shown in Fig. 8: while for $L=4$, $d$-wave and $p$-wave structure factors are
degenerate, and dominate over $s$-wave, for larger systems $d$-wave pairing
become dominant as $\beta\to\infty$. This tendency is confirmed by the
behaviour of the pairing susceptibilities (both bare and effective), for
$L=10$; see Fig. 8(d) As discussed in Ref. Huang  _et al._ , 2013, $d$-wave
pairing should dominate in the $V<0$ region, due to the nesting of the Fermi
surface, while the $p$-wave would be favoured in absence of this feature.
Interestingly, this is somehow observed in our results of Fig. 8: the smaller
the system size is, the weaker nesting effects are, leading to a spurious
$p$-wave enhancement.

Figure 7: a) Spin structure factor, b) $d$-wave pairing correlation function
c) Average sign and d) $p$-wave pairing correlation function as a function of
$V$ for $U=0.8$ with $L=6$.  Figure 8:  (a)-(c) Pairing structure factors and
(d) pairing susceptibilities as functions of the inverse temperature for
different symmetries, for the sizes shown. All data are for $U=0.5$, and
$V=-0.5$.

In summary, we have have resolved long standing issues relative to the ground
state phase diagram of the half-filled extended Hubbard model on a square
lattice, such as the pairing symmetries of the superconducting phase, and
locations of phase boundaries in all four quadrants of the $U-V$ plane.
Indeed, through our DQMC simulations, accurate boundaries involving
antiferromagnetic, charge-density wave, $s$-wave and $d$-wave superconducting,
and phase-separated phases were determined. The presence of a $d$-wave
superconducting phase over a reasonably wide region in the parameter space
adds credence to the use of the EHM as a minimal single-band model to describe
the high-$T_{c}$ cuprates Scalapino (2012). Also worth noticing is that the
‘minus-sign’ problem of QMC simulations was overcome by performing extensive
simulations, and by using the recently proposed connection between critical
points and strong dips in $\langle\text{s}\rangle$. Indeed, the multitude of
phases in the diagram allowed us to verify that dips in
$\langle\text{s}\rangle$ only occur at transitions involving different
universality classes, such as AFM-CDW, AFM-SC, SC-PS, and CDW-SC, but not
between $s$\- and $d$-waves.

## ACKNOWLEDGMENTS

The authors are grateful to the Brazilian Agencies Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq), Coordenação de
Aperfeiçoamento de Pessoal de Ensino Superior (CAPES), and Instituto Nacional
de Ciência e Tecnologia de Informação Quântica (INCT-IQ) for funding this
project. N.C.C. acknowledges support from FAPERJ Grant No. E-26/200.258/2023 -
SEI-260003/000623/2023, and CNPq Grant No. 313065/2021-7.

## References

  * Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, “Possible high-Tc superconductivity in the Ba-La-Cu-O system,” Zeitschrift für Physik B Condensed Matter 64, 189–193 (1986).
  * Varma (2020) Chandra M. Varma, “Colloquium: Linear in temperature resistivity and associated mysteries including high temperature superconductivity,” Rev. Mod. Phys. 92, 031001 (2020).
  * Zhou  _et al._ (2021) Xingjiang Zhou, Wei-Sheng Lee, Masatoshi Imada, Nandini Trivedi, Philip Phillips, Hae-Young Kee, Päivi Törmä, and Mikhail Eremets, “High-temperature superconductivity,” Nature Reviews Physics 3, 462–465 (2021).
  * Anderson (1987) P. W. Anderson, “The resonating valence bond state in La2CuO4 and superconductivity,” Science 235, 1196–1198 (1987).
  * Emery (1987) V. J. Emery, “Theory of high-${\mathrm{T}}_{\mathrm{c}}$ superconductivity in oxides,” Phys. Rev. Lett. 58, 2794–2797 (1987).
  * Scalapino (2012) D. J. Scalapino, “A common thread: The pairing interaction for unconventional superconductors,” Rev. Mod. Phys. 84, 1383–1417 (2012).
  * Lin  _et al._ (1995) H. Q. Lin, E. R. Gagliano, D. K. Campbell, E. H. Fradkin, and J. E. Gubernatis, “The phase diagram of the one-dimensional extended Hubbard model,” in  _The Hubbard Model: Its Physics and Mathematical Physics_, edited by Dionys Baeriswyl, David K. Campbell, Jose M. P. Carmelo, Francisco Guinea, and Enrique Louis (Springer US, Boston, MA, 1995) pp. 315–326.
  * Xiao  _et al._ (2022) Bo Xiao, Javier Robledo Moreno, Matthew Fishman, Dries Sels, Ehsan Khatami, and Richard Scalettar, “Extracting Off-Diagonal Order from Diagonal Basis Measurements,” (2022), arXiv:2209.10565 [cond-mat.str-el] .
  * Zhang and Callaway (1989) Y. Zhang and J. Callaway, “Extended Hubbard model in two dimensions,” Phys. Rev. B 39, 9397–9404 (1989).
  * Huang  _et al._ (2013) Wen-Min Huang, Chen-Yen Lai, Chuntai Shi, and Shan-Wen Tsai, “Unconventional superconducting phases for the two-dimensional extended Hubbard model on a square lattice,” Phys. Rev. B 88, 054504 (2013).
  * Vandelli  _et al._ (2020) M. Vandelli, V. Harkov, E. A. Stepanov, J. Gukelberger, E. Kozik, A. Rubio, and A. I. Lichtenstein, “Dual boson diagrammatic Monte Carlo approach applied to the extended Hubbard model,” Phys. Rev. B 102, 195109 (2020).
  * Yao  _et al._ (2022) Meng Yao, Da Wang, and Qiang-Hua Wang, “Determinant quantum Monte Carlo for the half-filled Hubbard model with nonlocal density-density interactions,” Phys. Rev. B 106, 195121 (2022).
  * Wolf  _et al._ (2018) Sebastian Wolf, Thomas L. Schmidt, and Stephan Rachel, “Unconventional superconductivity in the extended Hubbard model: Weak-coupling renormalization group,” Phys. Rev. B 98, 174515 (2018).
  * Chen  _et al._ (2022) Wei-Chih Chen, Yao Wang, and Cheng-Chien Chen, “Superconducting phases of the square-lattice extended hubbard model,” (2022), arXiv:2206.01119 [cond-mat.supr-con] .
  * Blankenbecler  _et al._ (1981) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, “Monte Carlo calculations of coupled boson-fermion systems. I,” Phys. Rev. D 24, 2278–2286 (1981).
  * Hirsch (1985) J. E. Hirsch, “Two-dimensional Hubbard model: Numerical simulation study,” Phys. Rev. B 31, 4403–4419 (1985).
  * Scalettar  _et al._ (1989) R. T. Scalettar, N. E. Bickers, and D. J. Scalapino, “Competition of pairing and Peierls–charge-density-wave correlations in a two-dimensional electron-phonon model,” Phys. Rev. B 40, 197–200 (1989).
  * Kawashima (2002) Naoki Kawashima, “Quantum Monte Carlo Methods,” Progress of Theoretical Physics Supplement 145, 138–149 (2002).
  * dos Santos (2003) Raimundo R. dos Santos, “Introduction to quantum Monte Carlo simulations for fermionic systems,” Braz. J. Phys 33, 63–54 (2003).
  * Becca and Sorella (2017) Federico Becca and Sandro Sorella,  _Quantum Monte Carlo Approaches for Correlated Systems_ (Cambridge University Press, 2017).
  * Sushchyev and Wessel (2022) Alexander Sushchyev and Stefan Wessel, “Thermodynamics of the metal-insulator transition in the extended Hubbard model from determinantal quantum Monte Carlo,” Phys. Rev. B 106, 155121 (2022).
  * Mondaini  _et al._ (2022a) R. Mondaini, S. Tarat, and R. T. Scalettar, “Quantum critical points and the sign problem,” Science 375, 418–424 (2022a).
  * Mondaini  _et al._ (2022b) Rubem Mondaini, Sabyasachi Tarat, and Richard T. Scalettar, “Universality and Critical Exponents of the Fermion Sign Problem,” (2022b), arXiv:2207.09026 [cond-mat.str-el] .
  * Hirsch (1983) J. E. Hirsch, “Discrete Hubbard-Stratonovich transformation for fermion lattice models,” Phys. Rev. B 28, 4059–4061 (1983).
  * White  _et al._ (1989) S. R. White, D. J. Scalapino, R. L. Sugar, N. E. Bickers, and R. T. Scalettar, “Attractive and repulsive pairing interaction vertices for the two-dimensional Hubbard model,” Phys. Rev. B 39, 839–842 (1989).
  * Hirsch (1984) J. E. Hirsch, “Charge-Density-Wave to Spin-Density-Wave Transition in the Extended Hubbard Model,” Phys. Rev. Lett. 53, 2327–2330 (1984).
  * Micnas  _et al._ (1990) R. Micnas, J. Ranninger, and S. Robaszkiewicz, “Superconductivity in narrow-band systems with local nonretarded attractive interactions,” Rev. Mod. Phys. 62, 113–171 (1990).
  * Fontenele  _et al._ (2022) Rodrigo A. Fontenele, Natanael C. Costa, Raimundo R. dos Santos, and Thereza Paiva, “Two-dimensional attractive Hubbard model and the BCS-BEC crossover,” Phys. Rev. B 105, 184502 (2022).
  * Rademaker  _et al._ (2013) Louk Rademaker, Steve Johnston, Jan Zaanen, and Jeroen van den Brink, “Determinant quantum Monte Carlo study of exciton condensation in the bilayer Hubbard model,” Phys. Rev. B 88, 235115 (2013).
  * Gubernatis  _et al._ (2016) James Gubernatis, Naoki Kawashima, and Philipp Werner, “Determinant method,” in  _Quantum Monte Carlo Methods: Algorithms for Lattice Models_ (Cambridge University Press, 2016) p. 180–213.

Supplemental Material for: Phase diagram for the extended Hubbard model on a
square lattice

The DQMC algorithm: To take into account the interaction terms of (1) one
first performs a Suzuki-Trotter decomposition,

$e^{-\beta\mathcal{H}}\approx\left(e^{-\Delta\tau\mathcal{H}_{k}}e^{-\Delta\tau\mathcal{H}_{U}}e^{-\Delta\tau\mathcal{H}_{V}}\right)^{M}+\mathcal{O}(\Delta\tau^{2})$
(S1)

where the 2D problem is mapped into a 3D system with the new dimension being
the imaginary time axis, and the inverse of the temperature $\beta$ is cut
into $M$ discrete intervals with length $\Delta\tau=\beta/M$. The quartic
terms of $\mathcal{H}_{U}$ and $\mathcal{H}_{V}$ are expressed in a quadratic
form using the Hubbard-Stratonovich transformation Hirsch (1983); Zhang and
Callaway (1989),

$e^{-\Delta\tau
Wn_{\mathbf{i}\sigma}n_{\mathbf{j}\sigma^{\prime}}}=\frac{1}{2}\sum_{x_{\nu}=\pm
1}e^{\left[\alpha
x_{\nu}(n_{\mathbf{i}\sigma}-n_{\mathbf{j}\sigma^{\prime}})-\frac{\Delta\tau
W}{2}(n_{\mathbf{i}\sigma}+n_{\mathbf{j}\sigma^{\prime}})\right]}$ (S2)

with $\cosh{\alpha}=\exp{(\Delta\tau W/2)}$, and $W=\\{U,V\\}$ depending if it
is the on-site or the extended Hubbard term. The second term in the argument
of the exponential of Eq. (S2) is cancelled by setting the chemical potential
$\mu=U/2+4V$ at half-filling. One auxiliary field $x_{0}(\mathbf{i},l)$ is
defined for each lattice site $\mathbf{i}$ and time slice $l$, for the on-site
interaction of $\mathcal{H}_{U}$. Since $\mathcal{H}_{V}$ is expanded as

$Vn_{\mathbf{i}}n_{\mathbf{j}}=V(n_{\mathbf{i}\uparrow}n_{\mathbf{j}\uparrow}+n_{\mathbf{i}\uparrow}n_{\mathbf{j}\downarrow}+n_{\mathbf{i}\downarrow}n_{\mathbf{j}\uparrow}+n_{\mathbf{i}\downarrow}n_{\mathbf{j}\downarrow}),$
(S3)

four auxiliary fields for each bond between nearest neighbors sites
$\mathbf{i}$ and $\mathbf{j}$ are defined, namely
$x_{1}(\mathbf{i},\mathbf{j},l)$, $x_{2}(\mathbf{i},\mathbf{j},l)$,
$x_{3}(\mathbf{i},\mathbf{j},l)$, and $x_{4}(\mathbf{i},\mathbf{j},l)$.
Working on a square lattice with periodic boundary conditions with $N$ sites,
and consequently $2\times N$ bonds, results in a total of $8\times N\times M$
auxiliary fields for the extended term in addition to the $N\times M$
auxiliary fields for the on-site term. The partition function of the system
then becomes,

$\mathcal{Z}=\dfrac{1}{2^{9NM}}\sum_{x_{\nu}=\pm
1}\prod_{\sigma}\det{A^{\sigma}(x_{0},x_{1},x_{2},x_{3},x_{4})},$ (S4)

with,

$A^{\sigma}(x_{0},x_{1},x_{2},x_{3},x_{4})=\mathbb{1}+B^{\sigma}(M)B^{\sigma}(M-1)\cdots
B^{\sigma}(1)$ (S5)

and

$B^{\sigma}(l)=e^{-\Delta\tau\hat{H}_{k}}e^{-\Delta\tau\hat{H}_{U}}e^{-\Delta\tau\hat{H}_{V}}.$
(S6)

The weight $P(x_{\nu})$ of each configuration of auxiliary fields is given by

$P(x_{\nu})=\left|\det\left[A^{\uparrow}A^{\downarrow}\right]\right|.$ (S7)

One can also compute Green’s functions Rademaker  _et al._ (2013),

$G_{\mathbf{i},\mathbf{j}}^{\sigma}=\dfrac{1}{2^{9NM}}\sum_{x_{\nu}=\pm
1}\left[A^{\sigma}\right]^{-1}_{\mathbf{i},\mathbf{j}}\prod_{\sigma}\det{A^{\sigma}(x_{0},x_{1},x_{2},x_{3},x_{4})},$
(S8)

and calculate the relevant physical quantities. Since the number of Hubbard-
Stratonovich fields (HSF) grows as $2^{9NM}$, Ising importance sampling Monte
Carlo over $\\{x_{\nu}\\}$ is used to analyze the Extended Hubbard model
properties. Details of the sampling techniques are discussed in the next
topic.

Sampling the Hubbard-Stratonovich fields: The Ising sampling of the HSF
consists in sweeps over the sites and bonds trying to make the flip
$x_{\nu}\rightarrow-x_{\nu}$ in each imaginary time slice $l$. From time to
time the physical quantities are measured. For instance, a flip of the $x_{0}$
field at a single site and time slice is proposed: if the change is accepted,
the new Green’s function is computed by $O(N^{2})$ operations using the
Sherman-Morrison update Gubernatis  _et al._ (2016), since the evaluation of
the Green’s function from Eq. (S8) requires $O(N^{3})$ operations. The same
approach can be used to update the Green’s functions by flipping the $x_{2}$
and $x_{3}$ fields, which couples the terms with $\sigma\neq\sigma^{\prime}$
in Eq. (S2), and consequently allows us to use Sherman-Morrison update for
$G^{\uparrow}$ and $G^{\downarrow}$ separately. The challenge is to update the
Green’s functions by flipping the $x_{1}$ and $x_{4}$ fields since they change
the Green’s functions related to different sites for $\sigma=\sigma^{\prime}$.
Following the procedure described in Rademaker  _et al._ (2013) we use the
Woodbury matrix identity, which is a generalized Sherman-Morrison update. The
steps to implement this update are detailed in this section.

Starting with a known Green’s function at a time slice $l$ given the fields
$x_{1(4)}$, we propose a flip $x_{1(4)}\rightarrow
x_{1(4)}^{\prime}=-x_{1(4)}$ in a bond between sites $\mathbf{i}$ and
$\mathbf{j}$ for a time slice $l$ individually. Under this change the
$A^{\sigma}$ matrix becomes,

$A^{\sigma}\rightarrow\mathbb{1}-B^{\sigma}(M)\cdots\left(\mathbb{1}+\Delta^{\sigma}\right)B^{\sigma}(l)\cdots
B^{\sigma}(1)=A^{\prime\sigma}.$ (S9)

The label $\sigma$ corresponds to spin up and down for $x_{1}$ and $x_{4}$
respectively. The matrix $\Delta^{\sigma}$ is diagonal and has only two
nonzero elements, namely

$\begin{split}\Delta_{i,i}^{\sigma}&=\exp\left[-2x_{\nu}(\mathbf{i},\mathbf{j},l)\right]-1\\\
\Delta_{j,j}^{\sigma}&=\exp\left[2x_{\nu}(\mathbf{i},\mathbf{j},l)\right]-1.\end{split}$
(S10)

One can decide whether or not to accept a change using the ratio of
determinants Gubernatis  _et al._ (2016),

$R^{\sigma}=\dfrac{\det\left[A^{\prime\sigma}\right]}{\det\left[A^{\sigma}\right]}=\det\left[\mathbb{1}-\Delta^{\sigma}(\mathbb{1}-G^{\sigma})\right].$
(S11)

If the flipping is accepted, by using matrix identities it is possible to
write a simple expression for the updated $G^{\sigma\prime}$ after the flip
Gubernatis  _et al._ (2016),

$G^{\prime\sigma}\rightarrow
G^{\sigma}-G^{\sigma}\dfrac{1}{\mathcal{R}}\Delta^{\sigma}(\mathbb{1}-G^{\sigma}).$
(S12)

Due to the sparseness of the matrices $\mathcal{R}$ and $\Delta^{\sigma}$, Eq.
(S12) can be written in terms of matrix elements Rademaker  _et al._ (2013),

$G^{\prime\sigma}_{r,s}\rightarrow
G^{\sigma}_{r,s}-\sum_{i,j}G^{\sigma}_{r,i}\mathcal{D}_{i,j}(\delta_{j,s}-G^{\sigma}_{j,s}),$
(S13)

where the matrix $\mathcal{D}$ may be cast in a $2\times 2$ form,

$\mathcal{D}=\dfrac{\Delta^{\sigma}_{i,i}\Delta^{\sigma}_{j,j}}{R^{\sigma}}\begin{pmatrix}1-G_{j,j}^{\sigma}+\dfrac{1}{\Delta^{\sigma}_{j,j}}&G^{\sigma}_{i,j}\\\
&\\\
G^{\sigma}_{j,i}&1-G_{i,i}^{\sigma}+\dfrac{1}{\Delta^{\sigma}_{i,i}}\end{pmatrix}$
(S14)

with $i$ and $j$ being the indices of the sites related with the flip of
$x_{\nu}(\mathbf{i},\mathbf{j},l)$, and,

$\begin{split}R^{\sigma}=\det\mathcal{R}=&(\Delta^{\sigma}_{i,i}(1-G_{i,i}^{\sigma})+1)(\Delta^{\sigma}_{j,j}(1-G_{j,j}^{\sigma})+1)\\\
&-\Delta^{\sigma}_{i,i}G_{i,j}^{\sigma}\Delta^{\sigma}_{j,j}G_{j,i}^{\sigma}.\end{split}$
(S15)

In order to decide whether a proposed change is accepted we used a combination
of the Metropolis and the heat-bath algorithm, similar to the one used in Ref.
Rademaker  _et al._ (2013),

$P_{T}=\left\\{\begin{array}[]{ll}\dfrac{R^{\sigma}}{1+\gamma
R^{\sigma}}&\text{if }R^{\sigma}\leq 1\\\ &\\\
\dfrac{R^{\sigma}}{\gamma+R^{\sigma}}&\text{if
}R^{\sigma}>1\end{array}\right.$ (S16)

The parameter $\gamma$ is tuned self-consistently to achieve an acceptance
ratio of approximated $50\%$.

The order of the flipping attempts is also an important detail in the DQMC
method. To achieve better results, four different strategies were tested and
the one which provided the best outcome is explained below. A sweep is done
over all sites $\mathbf{i}$ trying to flip $x_{0}(\mathbf{i},l)$ for
$l=0,1,2,\cdots,M-1$, after this, a sweep is done over all bonds trying to
flip $x_{1}(\mathbf{i},\mathbf{j},l)$ for $l=0,1,2,\cdots,M-1$ and the same is
done for $x_{2}$, $x_{3}$ and $x_{4}$ to complete one Monte Carlo step. Trying
to flip all five types of auxiliary fields at the same time, while speeds up
the sweep, leads to spurious results for our simulation.

Using the Woodbury update we have performed a few tests in order to verify the
accuracy of the results, where the main source of data for comparison is from
Ref. Zhang and Callaway (1989). All results were obtained on a $4\times 4$
square lattice with periodic boundary conditions. The parameter of the Trotter
decomposition was chosen to be $\Delta\tau=0.05$. A comparison with exact
diagonalization was also done.