# Bands renormalization and superconductivity in the strongly correlated
Hubbard model using composite operators method

L. Haurie  Institut de Physique Théorique, Université Paris Saclay, CEA CNRS,
Orme des Merisiers, 91190 Gif-sur-Yvette Cedex, France  M. Grandadam
Department of Physics and Physical Oceanography, Memorial University of
Newfoundland, St. John’s, Newfoundland & Labrador, Canada A1B 3X7  E. Pangburn
Institut de Physique Théorique, Université Paris Saclay, CEA CNRS, Orme des
Merisiers, 91190 Gif-sur-Yvette Cedex, France  A. Banerjee  Department of
Physics, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel  S.
Burdin  Université de Bordeaux, CNRS, LOMA, UMR 5798, F-33400 Talence, France
C. Pépin  Institut de Physique Théorique, Université Paris Saclay, CEA CNRS,
Orme des Merisiers, 91190 Gif-sur-Yvette Cedex, France

(August 27, 2024)

###### Abstract

We use the composite operator method (COM) to analyze the strongly correlated
repulsive Hubbard model, investigating the effect of nearest-neighbor hoppings
up to fourth order on a square lattice. We consider two sets of self-
consistent equations, one enforcing the Pauli principle and the other imposing
charge-charge, spin-spin, and pair-pair correlations using a decoupling scheme
developed by L. Roth [Phys. Rev. 184, 451–459, (1969)]. We extract three
distinct solutions from these equations: COM1 and COM2 by imposing the Pauli
principle and one from Roth decoupling. An overview of the method studying the
validity of particle-hole symmetry and the Luttinger theorem for each solution
is presented. Additionally, we extend the initial basis to study
superconductivity, concluding that it is induced by the Van Hove singularity.
Finally, we include higher-order hoppings using realistic estimates for tight
binding parameters and compare our results with ARPES measurements on
cuprates.

##  I  Introduction

Exotic behavior and multiple phases exhibited by strongly correlated materials
have attracted significant interest in the quantum condensed matter
community[1, 2]. One prominent example is cuprates, where superconductivity
survives up to very high critical temperatures compared to conventional
superconductors described by the BCS theory [3]. Since their discovery in 1987
by G. Bednorz and K.A. Müller [4], extensive research has focused on modeling
and understanding the pairing mechanism in these materials to obtain a
microscopic theory of cuprate superconductors. The complexity of their phase
diagram and the numerous unexplained quantum phases [5] have resulted in
various theoretical models [6].

The repulsive Hubbard model remains the paradigmatic tool for studying
strongly correlated quantum matter even after sixty years of its proposal. It
offers a minimal description of such systems with just two ingrediants -
hopping of electrons between lattice sites $t$ and on-site electron repulsion
$U$ creating an energy cost for double occupancy on the same lattice site. It
is given by

$H=-\sum\limits_{i,j\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum\limits_{i}n_{i\uparrow}n_{i\downarrow}+\mu\sum\limits_{i\sigma}n_{i\sigma},$
(1)

where the chemical potential $\mu$ fixes the total number of electrons. In Eq.
(1) $t_{ij}$ is such that $t_{ij}=t$ if $i$ and $j$ are nearest neighbouring
sites, and zero otherwise. One can extend the model to longer ranged hoppings.
$c_{i\sigma}/c_{i\sigma}^{\dagger}$ destroys/creates an electron on site $i$
with spin $\sigma\in\\{\uparrow,\downarrow\\}$ on the lattice and
$n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ is the number operator for site
$i$ and spin $\sigma$. However, this seemingly straightforward model has
eluded unbiased theoretical solutions.

The Hubbard Hamiltonian is known to be solvable only in one dimension (1d) and
infinite dimensions [7]. In the intermediate dimensions, the equations of
motion of the electronic Green’s function involve higher order Green’s
function and cannot be solved exactly [8]. In the past, approximations such as
Hubbard I and II (Ref. [8]) have been made to truncate the higher-order terms
and close the equations of motion.

We revisit a perturbative scheme called ”Composite Operators Method” (COM) [9]
[10], which was developed to close the equations of motion of the electronic
Green’s function of the Hubbard model to derive sets of self-consistent
equations. We consider two sets of such self-consistency equations. The first
enforces charge-charge, spin-spin, and pair-pair correlations using L. Roth’s
decoupling scheme [11], and the second imposes the Pauli principle instead of
correlations. We perform a numerical minimization of these two sets of
equations. The self-consistency equations enforcing the Pauli principle yield
two distinct solutions named COM1 and COM2. The other sets of equations using
Roth decoupling only exhibit one solution, which is referred to as the Roth
solution. Despite violating the Pauli principle [9], the Roth solution shows a
band structure in agreement with quantum Monte-Carlo simulations [10]. In
contrast, COM1 and COM2, which have been studied more widely, exhibit a Fermi
surface consistent with cuprate ARPES experiments only under some
approximations [12].

Although the self-consistency equations enforcing the Pauli principle have
been previously minimized with next-nearest neighbor hoppings [13], this is
not the case for the Roth decoupling. The latter has been studied with next-
nearest neighbor hoppings for a three-band Hubbard model in cuprates [14], but
never for the one-orbital Hubbard model.

Superconductivity has also been studied using the COM[10]. Past studies
indicate that the proximity to the Van Hove singularity enhances SC within
Roth minimization scheme [15]. The enhancement of the density of the state
associated with this singularity allows electrons to form more pairs readily.
However, a detailed study comparing and contrasting different self-consistent
solutions and including realistic hopping parameters to model the cuprates is
still lacking.

The objective of this study is to analyze the strongly correlated regime of
the Hubbard model where $|U|\gg t$ using the composite operators method. We
aim to compare and contrast the solutions obtained in previous studies using
the COM. We evaluate their physical consistency with experiments in cuprates
by testing the particle-hole symmetry and the validity Luttinger theorem.
Finally, we include the superconductivity with longer-ranged hoppings. The
plan for the rest of the paper is as follows: Section II details the formalism
and self-consistency we use in the manuscript. In section III, we restrict to
the nearest neighbor hoppings and compare the solutions obtained by the
different minimizations [10] [16] while discussing their physical
implications, particularly with regard to particle-hole symmetry and the
Luttinger theorem. Finally, we allow for the superconductivity in the model,
and in the last section IV, we study the impact of the longer-ranged hopping
orders.

##  II Composite operator formalism

Composite operators are introduced in order to solve exactly the equations of
motion of composite Green’s functions (Green’s functions composed of two
composite operators) in some limits. Since we are interested in a strongly
correlated regime of Hubbard model, we introduce composite operators solving
exactly the Hubbard model in the atomic limit (no hopping term). This is a
good starting point to treat the hopping term $t$ with a perturbation
expansion.

###  II.1 Presentation of the method

We start by introducing the following composite operators

$\begin{cases}\xi_{i\sigma}&=c_{i\sigma}-c_{i\sigma}n_{i\bar{\sigma}}\\\
\eta_{i\sigma}&=c_{i\sigma}n_{i\bar{\sigma}}\end{cases}$ (2)

Where $\bar{\sigma}$ means we take a spin $\uparrow$ if $\sigma=\downarrow$
and a spin $\uparrow$ if $\sigma=\downarrow$. We introduce the following
spinors

$\displaystyle\psi_{i\sigma}=$ $\displaystyle\begin{pmatrix}\xi_{i\sigma}\\\
\eta_{i\sigma}\end{pmatrix}$ (3)

Hereafter, $\psi^{1}_{i\sigma}=\xi_{i\sigma}$ and
$\psi^{2}_{i\sigma}=\eta_{i\sigma}$ denotes the first and second component of
the spinor at site i and with spin $\sigma$.

####  II.1.1 Atomic limit

To illustrate the interest of introducing these operators, we first consider
the atomic limit of Eq. (1) : we put the hopping term $t$ to zero. The
Hamiltonian is now local (ie electrons are not hopping anymore and each site
is independent) and is given by

$H_{loc}=U\sum\limits_{i}n_{i\uparrow}n_{i\downarrow}+\mu\sum\limits_{i\sigma}n_{i\sigma}$
(4)

We introduce the $2\times 2$ composite Green’s function matrix at sites $i$
and $j$, with spins $\sigma$ and $\sigma^{\prime}$ and defined for an
imaginary time $\tau$ by

$\textbf{S}_{ij\sigma\sigma^{\prime}}^{loc}(\tau,\tau^{\prime})=\langle\langle\psi_{i\sigma}(\tau);\psi_{j\sigma^{\prime}}^{\dagger}(\tau^{\prime})\rangle\rangle_{loc}$
(5)

Where, for two operators X and Y,

$\langle\langle
X(\tau);Y(\tau^{\prime})\rangle\rangle_{loc}=\theta_{H}(\tau-\tau^{\prime})\langle\\{X(\tau);Y(\tau^{\prime})\\}\rangle_{loc}$
(6)

Where $\theta_{H}(\tau-\tau^{\prime})$ is one if $\tau>\tau^{\prime}$ and zero
otherwise (Heaviside function). $\langle...\rangle_{loc}$ denotes the thermal
average taken with the Hamiltonian $H_{loc}$ and
$\\{X(\tau);Y(\tau^{\prime})\\}$ is the anticommutator of X and Y. Since we
are at equilibrium we have
$\textbf{S}^{loc}_{ij\sigma\sigma^{\prime}}(\tau,\tau^{\prime})=\textbf{S}^{loc}_{ij\sigma\sigma^{\prime}}(\tau-\tau^{\prime})$.
By differentiating with respect to time, we get the following equations of
motion for the composite Green’s function matrix

$\displaystyle\frac{d}{d\tau}\textbf{S}^{loc}_{ij\sigma\sigma^{\prime}}(\tau)=$
$\displaystyle\delta(\tau)\delta_{\sigma\sigma^{\prime}}\langle\\{\psi_{i\sigma}(\tau);\psi_{j\sigma}^{\dagger}(0)\\}\rangle_{loc}$
(7)
$\displaystyle+\delta_{\sigma\sigma^{\prime}}\langle\langle[\psi_{i\sigma}(\tau);H_{loc}];\psi_{j\sigma}^{\dagger}(0)\rangle\rangle_{loc}$

Where $[A;B]$ is the usual commutator between two operators A and B. We
enforce a paramagnetic solution by adding $\delta_{\sigma\sigma^{\prime}}$
prefactor. The currents in the atomic limit are given by

$J^{loc}_{i\sigma}(\tau)=\frac{d}{d\tau}\psi_{i\sigma}(\tau)=[\psi_{i\sigma}(\tau);H_{loc}]=\textbf{A}\psi_{i\sigma}(\tau)$
(8)

With

$\textbf{A}=\begin{pmatrix}\mu&0\\\ 0&U-\mu\end{pmatrix}$ (9)

The equations of motion becomes

$\displaystyle\frac{d}{d\tau}$
$\displaystyle\textbf{S}^{loc}_{ij\sigma\sigma^{\prime}}(\tau)=$ (10)
$\displaystyle\delta(\tau)\delta_{\sigma\sigma^{\prime}}\langle\\{\psi_{i\sigma}(\tau);\psi_{j\sigma}^{\dagger}(0)\\}\rangle_{loc}-\textbf{A}\
\textbf{S}^{loc}_{ij\sigma}(\tau)$

By time fourier transform we get

$\textbf{S}^{loc}_{ij\sigma\sigma^{\prime}}(\omega)=\delta_{\sigma\sigma^{\prime}}(\omega-\textbf{A}+i0^{+})^{-1}\textbf{I}_{i\sigma}^{loc}\delta_{ij}$
(11)

With
$\textbf{I}_{i\sigma}^{loc}=\langle\\{\psi_{i\sigma};\psi_{i\sigma}^{\dagger}\\}\rangle_{loc}$
the normalization matrix and $0^{+}$ a small positive parameter used for
analytic continuation since we are working with Matsubara time $\tau$. This
matrix can be explicitly computed. A bit of algebra leads to

$\textbf{I}_{i\sigma}^{loc}=\begin{pmatrix}1-\langle
n_{i\sigma}\rangle_{loc}&0\\\ 0&\langle n_{i\sigma}\rangle_{loc}\end{pmatrix}$
(12)

We finally obtain

$\textbf{S}_{ij\sigma\sigma}^{loc}(\omega)=\delta_{ij}\delta_{\sigma\sigma^{\prime}}\begin{pmatrix}1-\frac{\langle
n_{i\sigma}\rangle_{loc}}{\omega-\mu+i0^{+}}&0\\\ 0&\frac{\langle
n_{i\sigma}\rangle_{loc}}{\omega-U+\mu+i0^{+}}\end{pmatrix}$ (13)

In the atomic limit the equations of motion can therefore be closed. The
solution is given in Eq. (11), and by using the relation between composite and
electronic operators $\xi_{i\sigma}+\eta_{i\sigma}=c_{i\sigma}$, we can deduce
the electronic Green’s function

$\displaystyle G_{ij\sigma\sigma^{\prime}}^{loc}(\tau)=$
$\displaystyle\delta_{\sigma\sigma^{\prime}}\langle\langle
c_{i\sigma}(\tau);c_{j\sigma^{\prime}}^{\dagger}\rangle\rangle_{loc}$ (14)
$\displaystyle=$ $\displaystyle\delta_{\sigma\sigma^{\prime}}(S^{11\
loc}_{ij\sigma}(\tau)+S^{12\ loc}_{ij\sigma}(\tau)$ $\displaystyle+S^{21\
loc}_{ij\sigma}(\tau)+S^{22\ loc}_{ij\sigma}(\tau))$

Where $S^{nm\
loc}_{ij\sigma}=\langle\langle\psi^{n}_{i\sigma}(\tau);\psi_{j\sigma}^{m}(0)\rangle\rangle_{loc}$.
Therefore we have showed that the composite operators we introduced solve
exactly the atomic limit. The electronic Green’s function can be directly
recovered, allowing to extract all of the information we want from the system.

####  II.1.2 General case

Let us consider the full Hamiltonian Eq. (1) that includes both the local term
$H_{loc}$ and the hopping term t. We consider the limit $U\gg t$ and we build
an approximation from the atomic limit. From now on, $\langle...\rangle$ are
the thermal averages taken with the full Hamiltonian. We then introduce

$\delta J_{i\sigma}=[\psi_{i\sigma},H-H_{loc}]=J_{i\sigma}-J_{i\sigma}^{loc}$
(15)

Where the current operator in the atomic limit is given by Eq. (8) and
$J_{i\sigma}=[\psi_{i\sigma},H]$ is the current operator taken with the full
Hubbard Hamiltonian Eq. (1). Because of the tight-binding term, higher order
Green’s functions will appear in the equations of motion of the composite
Green’s function. We will not be able to solve the problem exactly as in the
atomic limit. We thus need to do an approximation to be able to obtain the
composite Green’s function by truncating the equations of motion. Indeed, let
us consider the composite Green’s functions $2\times 2$ matrix
$\textbf{S}_{ij\sigma\sigma^{\prime}}(\tau)=\delta_{\sigma\sigma^{\prime}}\langle\langle\psi_{i\sigma}(\tau);\psi_{j\sigma^{\prime}}(0)\rangle\rangle$
with thermal average on the full Hamiltonian. Its equations of motion can be
written as

$\displaystyle\frac{d}{d\tau}\textbf{S}_{ij\sigma\sigma^{\prime}}(\tau)=$
$\displaystyle\delta_{\sigma\sigma^{\prime}}(\delta(\tau)\delta_{ij}\textbf{I}_{i\sigma}+\theta_{H}(\tau)\textbf{M}_{ij\sigma}(\tau))$
(16)

Where we introduced the normalization matrix I and the overlap matrix M
respectively as

$\displaystyle\textbf{I}_{i\sigma}=$
$\displaystyle\langle\\{\psi_{i\sigma},\psi_{i\sigma}^{\dagger}\\}\rangle=\delta_{ij}\begin{pmatrix}1-\langle
n_{i\sigma}\rangle&0\\\ 0&\langle n_{i\sigma}\rangle\end{pmatrix}$ (17)
$\displaystyle\textbf{M}_{ij\sigma}=$
$\displaystyle\langle\\{J_{i\sigma},\psi_{j\sigma}^{\dagger}\\}\rangle=\begin{pmatrix}m_{ij}^{11}&m_{ij}^{12}\\\
m_{ij}^{12}&m_{ij}^{22}\end{pmatrix}$

To solve the equations of motion for the composite Green’s function matrix, we
need to compute the I and M matrix. We directly computed the I matrix.
Hereafter we focus on the determination of M. The current of the total
Hamiltonian is given by

$\displaystyle J_{i\sigma}=$
$\displaystyle\sum\limits_{l}\textbf{E}_{il\sigma}\psi_{l\sigma}+\delta\phi_{i\sigma}$
(18) $\displaystyle\textbf{E}_{il\sigma}=$
$\displaystyle\textbf{A}\delta_{il}+\textbf{P}_{il\sigma}$

The E contains all the terms proportional to $\psi$, and
$\delta\phi_{i\sigma}$ contains all terms which are not. The A matrix
appearing in E is the contribution of the atomic limit terms. It is given by

$\textbf{A}=\begin{pmatrix}\mu&0\\\ 0&U-\mu\end{pmatrix}$ (19)

The P matrix appearing in E is defined by

$\textbf{P}_{ij\sigma}=\langle\\{\delta
J_{i\sigma};\psi^{\dagger}_{j\sigma}\\}\rangle\textbf{I}^{-1}_{j\sigma}$ (20)

P is the contribution of the terms proportional to $\psi$ in $\delta
J_{i\sigma}$.

With this rewriting, the M matrix is now given by the following expression

$\displaystyle\textbf{M}_{il\sigma}(\tau)=$
$\displaystyle\sum\limits_{j}\textbf{E}_{ij\sigma}\langle\\{\psi_{j\sigma}(\tau);\psi^{\dagger}_{l\sigma}\\}\rangle$
(21)
$\displaystyle+\langle\\{\delta\phi_{i\sigma}(\tau);\psi_{l\sigma}^{\dagger}\\}\rangle$

The first term is proportional to $\textbf{S}_{jl\sigma\sigma^{\prime}}$.
However the second term is not and will introduce higher-order Green’s
function in the equations of motion

$\displaystyle\frac{d}{d\tau}\textbf{S}_{ij\sigma\sigma^{\prime}}(\tau)$
$\displaystyle=\delta_{\sigma\sigma^{\prime}}(\delta(\tau)\delta_{ij}\textbf{I}_{i\sigma}+\sum\limits_{l}\textbf{E}_{il\sigma}\textbf{S}_{lj\sigma}(\tau)$
(22)
$\displaystyle+\theta_{H}(\tau)\langle\\{\delta\phi_{i\sigma}(\tau);\psi_{l\sigma}^{\dagger}\\}\rangle$

An approximation is needed: we will assume that $\delta\phi_{i\sigma}$ is
negligible. Therefore, we have

$\displaystyle\textbf{S}_{ij\sigma\sigma^{\prime}}(\omega)$
$\displaystyle\approx\delta_{\sigma\sigma^{\prime}}(\omega-\textbf{E}_{ij\sigma}+i0^{+})^{-1}\textbf{I}_{j\sigma}$
(23a) $\displaystyle J_{i\sigma}(\tau)$
$\displaystyle\approx\sum\limits_{l}\textbf{E}_{il\sigma}\psi_{l\sigma}(\tau)$
(23b) $\displaystyle\textbf{E}_{ij\sigma}$
$\displaystyle\approx\textbf{M}_{ij\sigma}(0)\textbf{I}_{jj\sigma}^{-1}$ (23c)

Where these three equations are direct consequences from the approximation.
The first one is obtained by applying a time Fourier transform on the equation
of motion. The second equation is the current from Eq. (18)and the last
equation is Eq. (21) at equal imaginary time (notice that
$\textbf{S}_{ij\sigma}(0)=\delta_{ij}\textbf{I}_{i\sigma}$).

In appendix A, this approximation is studied in depth and we detail the
physical consequences of neglecting $\delta\phi$. E acts as an effective
energy matrix. Note that Eq. (23b) is similar to a Schrödinger equation for
the composite operators.We can also include higher order terms in the basis to
go further in the approximation. This has been done in Ref. [16].

The equations of motion of the composite Green’s function only depends on the
E and I matrices under our approximation. In order to perform a self-
consistent scheme, we introduce the $2\times 2$ correlation function matrix

$\textbf{C}_{ij\sigma}=\langle\psi_{i\sigma}\psi_{j\sigma}^{\dagger}\rangle$
(24)

We want to find an expression of $\textbf{C}_{ij\sigma\sigma^{\prime}}$ as a
function of the eigenvalues of the E matrix. Using the spectral representation
we get

$\textbf{C}_{ij\sigma}=\int
d\omega(1-f_{D}(\omega))\left(-\frac{1}{\pi}\right)Im(\textbf{S}_{ij\sigma}(\omega))$
(25)

With $f_{D}=\frac{1}{Exp(\beta\omega)+1}$ the Fermi Dirac distribution. From
there, we can use Eq. (23a) to apply the residue theorem on S and express it
as a function of the eigenvalues of the E matrix (proof in appendix B)

$\textbf{S}_{ij\sigma}(\omega)=\sum\limits_{a=1}^{2}\frac{\bm{\kappa}^{a}_{ij\sigma}}{(\omega-\epsilon_{ij}^{a}+i0^{+})}$
(26)

With

$\bm{\kappa}^{a}_{ij\sigma}=\frac{(-1)^{a+1}Com(\epsilon_{ij\sigma}^{a}-\textbf{E}_{ij\sigma})^{T}\textbf{I}_{j\sigma}}{(\epsilon_{ij\sigma}^{1}-\epsilon_{ij\sigma}^{2})}$
(27)

Where $Com(\textbf{A}_{ij})=(-1)^{i+j}det(\textbf{A}^{\prime}_{ij})$ is the
cofactor matrix of A, and $\textbf{A}^{\prime}$ is the matrix that results
from deleting row i and column j of A. In these equations, $a\in\\{1,2\\}$ and
$\epsilon_{ij}^{1}$ and $\epsilon_{ij}^{2}$ are the two eigenvalues of the
energy matrix $\textbf{E}_{ij\sigma}$. Finally, combining Eq. (25) and (26)
leads to an expression of the correlation function as a function of
$\epsilon^{a}$

$\displaystyle\textbf{C}_{ij\sigma}=$
$\displaystyle\frac{1}{2}\left(1+\tanh{\frac{\beta\epsilon^{1}_{ij\sigma}}{2}}\right)\bm{\kappa}^{1}_{ij\sigma}$
(28)
$\displaystyle\frac{1}{2}\left(1+\tanh{\frac{\beta\epsilon^{2}_{ij\sigma}}{2}}\right)\bm{\kappa}^{2}_{ij\sigma}$

Where $\beta=\frac{1}{k_{B}T}$ comes from the Fermi Dirac distribution. To
summarize, from the M and I matrices we can obtain the energy matrix E. We
showed that its eigenvalues are directly related to the composite correlation
function $\textbf{C}_{ij\sigma\sigma}$. To close the system, in the next
section we explicitly compute the M and I matrices and express them as a
function of the correlation functions.

###  II.2 Self-consistent scheme

####  II.2.1 Self-consistent equations of the one-body parameters

We can compute the algebra associated to our composite operators defined in
Eq. (2). We assume spin rotational symmetric solutions,, allowing us to have
$\delta_{\sigma\sigma^{\prime}}$ in our Green’s and correlation functions.
This implies

$\langle n_{i\uparrow}\rangle=\langle n_{i\downarrow}\rangle=\frac{n_{i}}{2}$
(29)

We compute explicitly the currents

$\displaystyle j_{i\sigma}^{1}=$
$\displaystyle-\mu\xi_{i\sigma}-\sum\limits_{l}t_{il}\left(c_{l\sigma}-n_{i\bar{\sigma}}c_{l\sigma}+S_{i}^{sign(\bar{\sigma})}c_{l\bar{\sigma}}+sign(\bar{\sigma})\Delta_{i}c^{\dagger}_{l\bar{\sigma}}\right)$
(30) $\displaystyle j_{i\sigma}^{2}=$
$\displaystyle-(\mu-U)\eta_{i\sigma}+\sum\limits_{l}t_{il}\left(-n_{i\bar{\sigma}}c_{l\sigma}+S_{i}^{sign(\bar{\sigma})}c_{l\bar{\sigma}}+sign(\bar{\sigma})\Delta_{i}c^{\dagger}_{l\bar{\sigma}}\right)$

With $S_{i}^{-}=c^{\dagger}_{i\downarrow}c_{i\uparrow}$,
$S_{i}^{+}=c^{\dagger}_{i\uparrow}c_{i\downarrow}$ and
$\Delta_{i}=c_{i\uparrow}c_{i\downarrow}$. We take the following convention:
$sign(\uparrow)=1$ and $sign(\downarrow)=-1$.

The coefficients of the M matrix are (some details of the algebra can be found
in appendix C)

$\displaystyle m_{ij}^{11}=$
$\displaystyle-\mu\left(1-\frac{n_{i}}{2}\right)\delta_{ij}-\delta_{ij}\sum\limits_{l}t\alpha_{il}^{1}e_{il}-t\alpha_{ij}^{1}(1-n_{i}+p_{ij})$
(31) $\displaystyle m_{ij}^{12}=$
$\displaystyle\delta_{ij}\sum\limits_{l}t\alpha_{il}^{1}e_{il}-t\alpha_{ij}^{1}\left(\frac{n_{i}}{2}-p_{ij}\right)$
$\displaystyle m_{ij}^{22}=$
$\displaystyle-(\mu-U)\frac{n_{i}}{2}\delta_{ij}-\delta_{ij}\sum\limits_{l}t\alpha_{il}^{1}e_{il}-t\alpha_{ij}^{1}p_{ij}$

With

$\displaystyle e_{ij}$
$\displaystyle=\langle\xi_{j\sigma}\xi^{\dagger}_{i\sigma}\rangle-\langle\eta_{j\sigma}\eta_{i\sigma}^{\dagger}\rangle$
(32) $\displaystyle p_{ij}$ $\displaystyle=\langle
n_{i\sigma}n_{j\sigma}\rangle+\langle
S_{i}^{-}S_{j}^{+}\rangle-\langle\Delta_{i}\Delta_{j}^{*}\rangle$

$\alpha_{il}^{1}=\frac{t_{il}}{t}$ is equal to 1 if $i$ and $l$ are nearest
neighbors, and 0 otherwise.

The parameter $e$ contains correlations between neighboring composite
operators and will mainly re-normalize the chemical potential because it
always appear in M in front of a $\delta_{ij}$. The $p$ parameter contains
charge-charge, spin-spin and pair-pair correlations and will affect the
bandwidth.

We now use the translational invariance of the system: each parameters with a
dependence in two sites $i$ and $j$ will only depend on the distance between
these two sites $i-j$. Hence we get

$\displaystyle e_{ij}=$ $\displaystyle e_{i-j}$ (33) $\displaystyle p_{ij}=$
$\displaystyle p_{i-j}$ $\displaystyle\textbf{M}_{ij}=$
$\displaystyle\textbf{M}_{i-j}$

We assume lattice inversion and rotation ($C_{4}$) symmetries and we take n, p
and e as constants

$\displaystyle n_{i}=n$ (34) $\displaystyle p_{i\pm\delta_{x}}=$
$\displaystyle p_{i\pm\delta_{y}}=p$ $\displaystyle e_{i\pm\delta_{x}}=$
$\displaystyle e_{i\pm\delta_{y}}=e$

Under these assumptions the coefficients of the M matrix become

$\displaystyle m_{ij}^{11}=$
$\displaystyle-\mu\left(1-\frac{n}{2}\right)\delta_{ij}-\delta_{ij}4te-t\alpha_{ij}^{1}(1-n+p)$
(35) $\displaystyle m_{ij}^{12}=$
$\displaystyle\delta_{ij}4te-t\alpha_{ij}^{1}\left(\frac{n}{2}-p\right)$
$\displaystyle m_{ij}^{22}=$
$\displaystyle-(\mu-U)\frac{n}{2}\delta_{ij}-\delta_{ij}4te-t\alpha_{ij}^{1}p$

And the E matrix defined by Eq. (23c) is given by

$\textbf{E}_{ij}=\begin{pmatrix}\frac{2}{2-n}m_{ij}^{11}&\frac{2}{n}m_{ij}^{12}\\\
\frac{2}{2-n}m_{ij}^{12}&\frac{2}{n}m_{ij}^{22}\end{pmatrix}$ (36)

This matrix is diagonalizable and Eq. (28) allows to express the correlation
function
$\textbf{C}_{ij\sigma}=\langle\psi_{i\sigma}\psi_{j\sigma}^{\dagger}\rangle$
as a function of its eigenvalues. Thus in order to close the system and be
able to solve it self-consistently we need to express the parameters in the M
matrix as a function of the correlation functions. A diagram of the self-
consistent loop is given in Fig. 1 for clarity.

Since $e$ and $n$ are one-body parameters, therefore they can directly be
expressed as

$\displaystyle n=$ $\displaystyle 2(1-C^{11}_{0}-2C^{12}_{0}-C^{22}_{0})$ (37)
$\displaystyle e=$ $\displaystyle C^{11}_{ij}-C^{22}_{ij}$

In this equation $\textbf{C}_{0}=\textbf{C}_{ii}$ are constants by
translational invariance. and $\textbf{C}=\textbf{C}_{ij}=\textbf{C}_{i-j}$ is
treated as e and p in Eq. (34) using lattice inversion and rotation ($C_{4}$)
symmetry, so it is also a constant but different from $\textbf{C}_{0}$.

Expressing $p$ as a function of the correlations functions is not so direct
since $p$ is composed of two-bodies operators while $e$, $n$ and $C$ are one-
body operators.

####  II.2.2 Self-consistent equations of two-bodies operators : Roth
decoupling

Following the step of L. Roth [11], we express p as a function of correlation
functions by mean of equations of motion. We give here the main idea of the
decoupling, taking the pair-pair term $\langle\Delta_{i}\Delta_{j}^{*}\rangle$
as an example. The full details of the computation can be found in Appendix D.
First, notice we can write

$\langle\Delta_{i}\Delta_{j}^{*}\rangle=\langle\xi_{i\uparrow}c_{i\downarrow}\Delta_{j}^{*}\rangle+\langle\eta_{i\uparrow}c_{i\downarrow}\Delta_{j}^{*}\rangle$
(38)

Note we illustrate this decoupling with $\xi_{i\uparrow}$ and
$\eta_{i\uparrow}$ but the idea is exactly the same with $\psi_{i\downarrow}$.
We introduce the following Green’s functions

$\begin{cases}F_{ijl}(\tau)&=\langle\langle\xi_{i\uparrow}(\tau);c_{j\downarrow}\Delta_{l}^{*}\rangle\rangle\\\
G_{ijl}(\tau)&=\langle\langle\eta_{i\uparrow}(\tau);c_{j\downarrow}\Delta_{l}^{*}\rangle\rangle\end{cases}$
(39)

Figure 1: The self consistency loop. Our hypothesis on the currents allows us
to get a relation between the M and I matrices (Eq. (23c)). The energy matrix
E is diagonalizable. We can express correlation functions as a function of its
eigenvalues with Eq. (28). Then, we can rewrite the M and I matrices in term
of these correlations functions and do a self-consistency (Eq. (45)).

We then consider the equations of motion for these Green’s functions

$\displaystyle\partial_{\tau}\begin{pmatrix}F_{ijl}(\tau)\\\
G_{ijl}(\tau)\end{pmatrix}$
$\displaystyle=\theta_{H}(\tau)\begin{pmatrix}\langle\\{\partial_{t}\xi_{i\uparrow}(\tau);c_{j\downarrow}\Delta_{l}^{*}\\}\rangle\\\
\langle\\{\partial_{\tau}\eta_{i\uparrow}(\tau);c_{j\downarrow}\Delta_{l}^{*}\\}\rangle\end{pmatrix}$
(40) $\displaystyle+\delta(\tau)\begin{pmatrix}f^{1}_{ijl}\\\
f^{2}_{ijl}\end{pmatrix}$

where
$f^{n}_{ijl}=\langle\\{\psi_{i}^{n};c_{j\downarrow}\Delta_{k}^{*}\\}\rangle$.
We now use Eq. (23b) to obtain

$\partial_{\tau}\begin{pmatrix}F_{ijl}(\tau)\\\
G_{ijl}(\tau)\end{pmatrix}=\sum\limits_{k}E_{ik}\begin{pmatrix}F_{kjl}(\tau)\\\
G_{kjl}(\tau)\end{pmatrix}+\delta(\tau)\begin{pmatrix}f^{1}_{ijl}\\\
f^{2}_{ijl}\end{pmatrix}$ (41)

We then time Fourier transform and use Eq. (23a) to obtain

$\begin{pmatrix}F_{ijl}(\omega)\\\
G_{ijl}(\omega)\end{pmatrix}=\sum\limits_{k}S_{ik}(\omega)I^{-1}\begin{pmatrix}f^{1}_{kjl}\\\
f^{2}_{kjl}\end{pmatrix}$ (42)

Finally, we can extract $\langle\Delta_{i}\Delta_{j}^{*}\rangle$ by summing F
and G, integrating over $\omega$ and taking the imaginary part to use Eq. (25)
in order to replace the composite Green’s functions by correlation functions.
We get

$\displaystyle\langle c_{i\sigma}c_{j\bar{\sigma}}\Delta_{l}^{*}\rangle$
$\displaystyle=\frac{2}{2-n}\sum\limits_{k}(C^{11}_{ik}+C^{12}_{ik})f^{1}_{kjl}$
(43)
$\displaystyle+\frac{2}{n}\sum\limits_{k}(C^{12}_{ik}+C^{22}_{ik})f^{2}_{kjl}$

The computation of
$f^{n}_{kjl}=\langle\\{\psi_{k}^{n};c_{j\downarrow}\Delta_{l}^{*}\\}\rangle$
is done in appendix D. Setting $i=j$, we finally get

$\langle\Delta_{i}\Delta_{l}^{*}\rangle=\frac{4}{n(2-n)}\frac{(C^{11}_{il}+C^{12}_{il})(C^{22}_{il}+C^{21}_{il})}{1-\frac{2}{2-n}(C^{11}_{0}+C^{12}_{0})+\frac{2}{n}(C^{21}_{0}+C^{22}_{0})}$
(44)

The same method can be applied for $\langle n_{i\sigma}n_{l\sigma}\rangle$ and
$\langle S_{i}^{+}S_{l}^{-}\rangle$, the details are in appendix D. Note that
we always obtain numerically $C^{21}_{0}$ = $C^{12}_{0}$.
Finally, the expressions for the parameters appearing in M and I as a function
of the $C^{nm}_{ij}=\langle\psi^{n}_{i}(\psi_{j}^{m})^{\dagger}\rangle$ are
the following

$\begin{cases}e&=C^{11}-C^{22}\\\ n&=2(1-C^{11}_{0}-2C_{0}^{12}-C_{0}^{22})\\\
p&=\frac{n^{2}}{4}-\frac{\rho_{1}}{1-\phi^{2}}-\frac{\rho_{1}}{1-\phi}-\frac{\rho_{3}}{1+\phi}\end{cases}$
(45)

With

$\begin{cases}\phi&=\frac{2}{2-n}(C^{11}_{0}+C^{12}_{0})+\frac{2}{n}(C^{12}_{0}+C^{22}_{0})\\\
\rho_{1}&=\frac{2}{2-n}(C^{11}+C^{12})^{2}+\frac{2}{n}(C^{22}+C^{12})^{2}\\\
\rho_{3}&=\frac{4}{n(2-n)}(C^{11}+C^{12})(C^{22}+C^{12})\end{cases}$ (46)

With the self consistent equations Eq. (45) on the parameters n, e and p we
just closed the system. In Fig. 1, we represented the self-consistency
pattern. Starting from initial guess for e, n and p, we compute the M and I
matrices. We can then obtain E and diagonalize it. Then, using Eq. (28), we
can express the correlation functions from its eigenvalues. Finally, using the
self consistent Eq. (45), we compute again e, p and n. We stop when
$f(x)-x<\delta$ where $x=(e,p,\mu)$ and f are given in Eq. (45). We chose
$\delta=10^{-8}$. Once the system converges, we use the parameters $(e,p,\mu)$
to compute the electronic Green’s function

$\displaystyle G_{ij}(\omega)=$ $\displaystyle
S_{ij}^{11}(\omega)+S_{ij}^{12}(\omega)+S_{ij}^{21}(\omega)+S_{ij}^{22}(\omega)$
(47) $\displaystyle=$
$\displaystyle\sum\limits_{l=1}^{2}\frac{(\kappa^{l})^{11}_{ij}+(\kappa^{l})^{12}_{ij}+(\kappa^{l})^{21}_{ij}+(\kappa^{l})^{22}_{ij}}{\omega-\epsilon^{l}_{ij}+i0^{+}}$

Where the $\bm{\kappa}$ act as spectral weights and are defined in Eq. (27)
and $\epsilon^{1/2}$ are the eigenvalues of the E matrix, and can be built
from the parameters $(e,p,\mu)$.

##  III  Nearest-neighbour study

This self-consistent scheme has a drawback. With the numerical solution of the
self consistent equations we obtain a non vanishing $C^{12}_{0}$. However,
analytically this correlation function is zero because of the Pauli principle.

$\displaystyle C_{0}^{12}$
$\displaystyle=\langle\xi_{i\sigma}\eta^{\dagger}_{i\sigma}\rangle$ (48)
$\displaystyle=\langle(c_{i\sigma}-c_{i\sigma}n_{i\bar{\sigma}})c^{\dagger}_{i\sigma}n_{i\bar{\sigma}}\rangle$
$\displaystyle=\langle
c_{i\sigma}c^{\dagger}_{i\sigma}n_{i\bar{\sigma}}+c_{i\sigma}c^{\dagger}_{i\bar{\sigma}}c^{\dagger}_{i\sigma}c_{i\bar{\sigma}}n_{i\bar{\sigma}}\rangle$
$\displaystyle=\langle-
c_{i\sigma}(c^{\dagger}_{i\bar{\sigma}})^{2}c^{\dagger}_{i\sigma}(c_{i\bar{\sigma}})^{2}\rangle$
$\displaystyle=0$

((a)) First Pauli Solution (COM1)

((b)) Second Pauli solution (COM2)

((c)) Roth Minimization.

Figure 3: Bands along high symmetry points with $t=1,U=8t,T=0,n=0.8$. The non
interacting (NI) band (red line) is splitted in 2 Hubbard bands obtained for
every solutions with the composite operator method (blue lines). Therefore at
half filling the chemical potential lies in between the two Hubbard bands and
we get a Mott insulator. (a) One of the two solutions obtained with Pauli
minimization (using Eq. (48) instead of the p parameter in Eq. (45)). It has
no renormalization from the interaction since the shape of the Hubbard bands
is similar to the non-interacting band. (b) The second solution obtained with
the Pauli minimization. It has a minimum at M=($\pi,\pi$) and exhibits two
hole pockets. (c) The solution obtained with the Roth minimization. We observe
a flattening of the bands around X=$(\pi,0)$.

This numerical violation of the Pauli principle gets smaller as U becomes
larger than t [9].

Figure 4: Fermi Surfaces obtained from the spectral function (imaginary part
of Eq. (47)) at $\omega=0$, associated to the bands depicted in (3). The
parameters are $t=1,n=0.8,U=8t$ and $T=0$. top left: COM1, top right: COM2,
bottom left: Roth, bottom right: Non interacting (tight binding)

It is possible to solve a different set of self-consistent equations by
imposing $C^{12}_{0}=0$ [16] as the third self consistent equation instead of
$p$ in Eq. (45). From now on, we will call the minimization with $p$ the “Roth
minimization” since it uses the decoupling formalized by L. Roth and the
minimization with $C^{12}_{0}$ the “Pauli minimization”.

Figure 5: Density of state in energy with parameters $t=1,U=8t$ and $T=0$,
using Roth minimization with nearest neighbors. From top to bottom the filling
is respectively: n=0.75, 0.88, 0.93, 0.98, 1 (half filling), 1.07 and 1.25.
The blue line on each panel corresponds to the chemical potential.

In the following section we consider both minimizations and discuss the bands
and Fermi surfaces obtained for each solution we found. We compare to the non
interacting tight-binding model, and study the particle-hole symmetry, the
Luttinger theorem as well as superconductivity.

###  III.1 Comparison of Pauli and Roth minimization

By varying the initial conditions and using a minimizer we isolate two
distinct solutions with the Pauli minimization, and one unique solution with
the Roth minimization. Following notations from Ref [9], we call these three
solutions respectively COM1, COM2 and Roth solutions. In Fig. 3 we plot the
bands along high symmetry points for these solutions.

These bands correspond to the eigenvalues of the E matrix and act as the poles
of the electronic Green’s function (cf Eq. (47)). All of these solutions have
two bands associated to the two eigenvalues of the energy matrix, which are
splitted by U. COM2 and Roth solution lead to a Mott insulator physics at
half-filling because the chemical potential lies in between the two bands.
This is not the case for COM1 which is not in an insulating phase at half-
filling for t=1, U=8t. The bands obtained by COM1 in figure 3a look exactly
like the non interacting (NI) tight-binding band for a square lattice
$\epsilon_{tb}=2t(cos(k_{x})+cos(k_{y}))$ (the red line), but splitted by U.
Therefore this solution is not really renormalized by the interaction. COM2 is
changed a lot compared to the non interacting case (Fig. 3b), but always
presents two holes pockets, leading to two Fermi Surfaces in Fig. 4. This is
unexpected, as this has never been observed by ARPES experiments for strongly
correlated materials such as cuprates where this treatment of the Hubbard
model is relevant. However, in Ref. [16] the basis has been extended to take
into account dynamical corrections of the self-energy (cf appendix A for
details). The lifetime of the second hole pocket is then computed and happens
to be small, which can explain why it is not observed experimentally. Finally,
the Roth solution (Fig. 3c) has a very different shape from the non
interactive solution. It also has the advantages of presenting only one hole
pocket and a maximum at $(\pi,\pi)$.

Note we plotted the bands we obtained at half filling for the Roth solution in
appendix E, on Fig. 20. Close to half filling (around 3% hole doping), the
Roth solution exhibits a second small hole pocket at $(\pi,\pi)$ (In appendix
E we plot the bands and Fermi Surface of Roth decoupling at 2% hole doping in
Fig. 19). This second hole pocket around $(\pi,\pi)$ may be the consequence of
our paramagnetic assumption $\langle n_{i\uparrow}\rangle=\langle
n_{i\downarrow}\rangle=\frac{n}{2}$. It appears around half filling where we
know the antiferromagnetic phase dominates [17]. The wave-vector $(\pi,\pi)$
is also associated with antiferromagnetism, so this second hole pocket might
be an instability of the system because we neglected it.

For comparison, in the Hubbard I approximation (the first and simplest
approximation developped by Hubbard in Ref. [8]) a simple factorization
procedure of the 2-bodies Green functions is used. According to Ref. [9],
Hubbard I approximation in the composite operator framework would be
equivalent to set

$\displaystyle e$ $\displaystyle=0$ (49) $\displaystyle p$
$\displaystyle=\frac{n^{2}}{4}$

It is known that the Hubbard I approximation does not exhibit any Mott-
insulator transition as long as $U\neq 0$ (cf Ref. [8]). Therefore the
composite operator method is a more refined approximation as Hubbard I.

Figure 6: Fermi Surface obtained with Roth solution, at n=0.7 and n=0.9. The
Van Hove singularity is at a chemical potential associated to electron density
n=0.8. Below the singularity the Fermi surface is centered around $(0,0)$, and
above the singularity it is centered around $(\pi,\pi)$

It is also instructive to consider the density of states. In Fig. 5, we plot
it as a function of energy for several doping with the Roth solution. Only
Roth solution is considered because it is the one who has the most reasonable
Fermi Surface, but similar behaviors are expected for COM1 and COM2. The
density of state has been computed from the spectral function using the
following formula

$\mathcal{D}(\omega)=\frac{1}{N_{k}^{2}}\sum\limits_{k_{x},k_{y}=1}^{N_{k}}\left(-\frac{1}{2\pi}\right)Im(G(\textbf{k},\omega))$
(50)

Where $N_{k}$ denotes the number of considered points for sampling $k_{x}$ and
$k_{y}$. At half filling we do not have any state at the Fermi energy, since
the model leads to an insulator for this doping. Around the gap two peaks can
be distinguished. For every doping except half filling, a third peak is
observable, and corresponds to the Van Hove singularity. It is associated to a
flattening of the band, meaning there are a lot of states associated to this
energy. In term of Fermi surface for a square lattice the Van Hove singularity
corresponds to the doping below which the Fermi Surface is centered on $(0,0)$
and above which it is centered on $(\pi,\pi)$, as shown in Fig. 6. On Fig. 3,
this flatness of the bands can be found for every solutions around $(\pi,0)$.
The Fermi energy is exactly at the Van Hove singularity for the Roth solution
at n=0.8 or n=1.2. On Fig. 3, the Roth solution is plotted exactly at this
doping and we can check the flat band around $(\pi,0)$ lies exactly at the
Fermi energy (corresponding to the dotted black line).

Regarding the Mott transition, we see at half filling no quasiparticle peak is
observed around the Fermi energy. Instead, the density of state is closer to
what was observed with determinantal quantum Monte Carlo simulation (DQMC) in
Ref. [18]. DQMC is a stochastic algorithm which allows under some limitations
to perform direct studies of complex condensed matter problems. As predicted
by DQMC in Ref [18], at high doping we have only one peak corresponding to a
Van Hove singularity (Fermi liquid behavior), and when approaching the Mott
transition a transfer of spectral weight occurs, changing the density of
states, without creating a quasiparticle peak at the Fermi energy at half
filling. Therefore in this regime where $U\gg t$ the density of states of the
lower and upper Hubbard bands are the only contribution.

((a)) Chemical potential

((b)) Parameter e

((c)) Parameter p

((d)) Luttinger Violation

Figure 8:  (a), (b) and (c): parameters as a function of doping for each
solution. The dash lines in the electron dope ($n>$1) region are the particle-
hole symmetric (54) of their equivalent in the hole doped region. Bottom
right: area of Fermi surface as a function of doping. It is computed with the
use of Eq. (57). We observe a violation of Luttinger theorem for all
solutions. All of these figures are obtained at $T=0$K ($\beta\approx
10^{5}t$) and $U=8t$.

Figure 9: Luttinger violation as a function of electron density. (left):
Figure from Ref. [18], obtained from Determinantal Monte Carlo, with
$\beta=2t$ and U=10t. (right): Figure we obtain with the method at the same
parameters for the three solutions. $n_{c}$ corresponds to the electron
density under which the Luttinger theorem is violated. Above $n_{c}$ the Fermi
surface volume is proportional to the electron density (the curves and the
line “FS=n” are parallel). Under $n_{c}$ this isn’t true anymore and we enter
in a “Luttinger breaking phase”.

###  III.2 Particle-hole symmetry and Luttinger theorem

####  III.2.1 Particle-hole symmetry

It is well known that the Hubbard model with nearest neighbour hoppings only
is particle-hole symmetric. This symmetry exchanges particles and holes with
the transformation

$c_{i\sigma}\rightarrow(-1)^{i}c_{i\sigma}^{\dagger}\ \ \ \ \ \ \ \ \ \ \ \ \
\ \ c_{i\sigma}^{\dagger}\rightarrow(-1)^{i}c_{i\sigma}$ (51)

We could also have taken another convention for this transformation without
any $(-1)^{i}$, as long as we change t to -t to keep the Hamiltonian
invariant. Applied to the composite operators one can show that it becomes

$\eta_{i\sigma}\rightarrow(-1)^{i}\xi_{i\sigma}^{\dagger}\ \ \ \
\xi_{i\sigma}\rightarrow(-1)^{i}\eta_{i\sigma}^{\dagger}$ (52)

Let us study the behaviour of the parameters $\mu$, p and e under the
particle-hole symmetry. With the paramagnetic assumption, the particle-hole
transformation is rewritten

$\displaystyle\langle c_{i\sigma}^{\dagger}c_{i\sigma}\rangle$
$\displaystyle\rightarrow\langle c_{i\sigma}c_{i\sigma}^{\dagger}\rangle$ (53)
$\displaystyle\Leftrightarrow\frac{n}{2}$ $\displaystyle\rightarrow
1-\frac{n}{2}$

Hence applying the particle-hole transformation on electronic filling gives
$n\rightarrow 2-n$. The transformation changes p and e (Eq. (32)) as follow
(details are given in appendix G)

$\displaystyle e(2-n)=$ $\displaystyle-e(n)$ (54) $\displaystyle p(2-n)=$
$\displaystyle p(n)+(1-n)$ $\displaystyle\mu(2-n)=$ $\displaystyle U-\mu(n)$

The relation of the chemical potential can be obtained by using the fact that
the Hubbard Hamiltonian must stay invariant under this symmetry for nearest-
neighbour hoppings. Finally, applying this transformation on the composite
bands $\epsilon^{1}$ and $\epsilon^{2}$ leads to

$\displaystyle\epsilon^{1}(2-n,\textbf{k})=$
$\displaystyle-\epsilon^{2}(n,\textbf{k}+(\pi,\pi))$ (55)
$\displaystyle\epsilon^{2}(2-n,\textbf{k})=$
$\displaystyle-\epsilon^{1}(n,\textbf{k}+(\pi,\pi))$

In Fig. 8a, 8b and 8c we plot e, p and $\mu$ as a function of doping for the
three solutions we studied (COM1, COM2 and Roth). The dashed-lines on the
electron-doped region is the value the parameter must have to satisfy the
particle-hole relations Eq. (54). We see that particle-hole symmetry is
respected for every parameters for the three solutions. The chemical potential
presents a jump of the order of U at half filling for the COM2 and Roth
solutions. This is because the lower Hubbard band is filled and the upper
Hubbard band starts to be occupied at half filling. COM1 does not exhibit this
feature.

We obtain a different result from Ref. [9]. It is possible to have a particle-
hole symmetric solution which violates Pauli principle. This is indeed the
case of the Roth solution which has a non vanishing $C^{12}_{0}$ despite the
fact it is zero analytically. The solution is particle-hole symmetric as long
as $C^{12}_{0}$ is not put to zero by hand in the self-consistent equations it
must appear both in the equation of n and in the $\phi$ term in the equation
of p.

####  III.2.2 Luttinger theorem

We now turn our attention to the Luttinger theorem. This theorem states that
the volume enclosed by the Fermi surface is proportional to the electron
density [19]. The Luttinger theorem is only demonstrated for systems where the
interaction can be treated perturbatively. For this reason it is common that
it can be violated in strongly correlated systems. To compute the volume
enclosed by the Fermi surface, we need to remember the relation between the
composite and electronic Green’s function Eq. (47) The $\bm{\kappa}$ act as
spectral weights of the electronic Green’s function. The Fermi surface is
given by the imaginary part of this electronic Green’s function at the Fermi
energy

$\displaystyle A_{ij}(\omega=0)=$
$\displaystyle-\frac{1}{\pi}Im(G_{ij}(\omega=0))$ (56) $\displaystyle=$
$\displaystyle\frac{1}{\pi}\sum\limits_{l}\frac{0^{+}((\kappa^{l})^{11}_{ij}+(\kappa^{l})^{12}_{ij}+(\kappa^{l})^{21}_{ij}+(\kappa^{l})^{22}_{ij})}{(\epsilon_{ij}^{l})^{2}+(0^{+})^{2}}$

((a)) Bands for the Roth solution

((b)) Bands for the COM2 solution

((c)) Fermi Surface for the Roth solution

((d)) Fermi Surface for the COM2 solution

Figure 10: Bands and Fermi Surface with superconductivity at $t=1$, $U=8t$ and
$T=0$. The plots on the left correspond to the Roth solution at $n=0.8$,
whereas the plot on the right correspond to COM2 at $n=0.9$. These are the
respective doping at which superconductivity is maximum. We observe a doubling
of the bands associated to the doubling of the basis. The insets on the top
plots correspond to a zoom around zero energy of the bands: we see a gap
opening at k=$(\pi,0)$. There is no gap opening between k=$(0,0)$ and
k=$(\pi,\pi)$ because of d-wave symmetry. Both Fermi surfaces on the bottom
plots have a loss of spectral weight around k=$(\pi,0)$.

Thus, the Fermi contour is the set of points associated to a vanishing
$\epsilon_{ij}^{l}$. Therefore a way to compute the enclosed area $A_{FS}$ is
simply by considering the following equation

$\displaystyle A_{FS}^{l}=$
$\displaystyle\frac{1}{N^{2}}\sum\limits_{k_{x},k_{y}=1}^{N}\theta_{H}(\epsilon^{l}(\textbf{k}))$
(57)

Note in this last equation we do not have the sum over the two eigenvalues.
Indeed, if hole doped, the upper Hubbard band $\epsilon^{2}(\textbf{k})$ is
never negative so the heaviside function returns zero. In the electron doped
regime, the lower Hubbard band is always filled, so it does not affect the
Fermi surface. Therefore, if hole doped we take l=1, and if electron doped we
take l=2. Fig. 8d reveals that the Luttinger theorem is violated. This
violation is analogous to what was observed by Ref. [18] using determinantal
quantum Monte-Carlo simulations. In Fig. 9, we plot side by side Luttinger
violation taken from this reference and the Luttinger violation we observed
with the same parameters (U=10t and $\beta$=2). While none of the solutions we
considered are exactly similar to what was observed with determinantal quantum
Monte-Carlo, we have the same behavior and conclusion : above a critical
electron density $n_{c}$, Luttinger theorem is verified and all the solutions
are proportionnal to the electron density (the curves are parallel to the
black dashed line). Under $n_{c}$ this is not the case anymore, we are in a
“Luttinger breaking phase”. This is explained in Ref. [18] as a topological
order which is a consequence of the proximity to the Mott transition.

###  III.3 Superconductivity and Van Hove singularity

####  III.3.1 Method

((a)) $\theta$ and $\Delta^{d}$ parameters as a function of electron density
$n$

((b)) Density of states at the Fermi Energy

Figure 11: (a): Anomalous superconducting mean field parameter $\theta$ as a
function of doping for Roth minimization and Pauli minimization (COM2
solution). $\theta$ satisfy particle-hole symmetry in both cases. The dashed
lines corresponds to the d-wave superconducting order parameter $\Delta^{d}$.
(b): Density of states at the Fermi energy with no superconductivity. We see a
clear correlation between enhancement of the density of states and
superconductivity. $\theta$ and $\Delta^{d}$ are maximum at the Van Hove
singularity. It lies at n=0.8 for Roth and around n=0.9 for COM2, as shown in
Fig. 10.

Superconductivity can be studied by extending the initial basis. The new
spinor $\psi$ to consider is

$\psi=\begin{pmatrix}\xi_{i\sigma}\\\ \eta_{i\sigma}\\\
\xi_{i\bar{\sigma}}^{\dagger}\\\ \eta_{i\bar{\sigma}}^{\dagger}\end{pmatrix}$
(58)

The method described before remains the same except for the I and M matrix
which are now $4\times 4$ matrices. Since we are interested in cuprates, we
consider only the case of d-wave superconductivity. Therefore, $\langle
c_{i\uparrow}c_{i\downarrow}\rangle=0$. In this framework, the M and I
matrices are given by

$\textbf{I}_{i}=\begin{pmatrix}1-\frac{n_{i}}{2}&0&0&0\\\
0&\frac{n_{i}}{2}&0&0\\\ 0&0&1-\frac{n_{i}}{2}&0\\\
0&0&0&\frac{n_{i}}{2}\end{pmatrix}\\\
\textbf{M}_{ij}=\begin{pmatrix}m_{ij}^{11}&m_{ij}^{12}&m_{ij}^{13}&-m_{ij}^{13}\\\
m_{ij}^{12}&m_{ij}^{22}&-m_{ij}^{13}&m_{ij}^{13}\\\
m_{ij}^{13}&-m_{ij}^{13}&-m_{ij}^{11}&-m_{ij}^{12}\\\
-m_{ij}^{13}&m_{ij}^{13}&-m_{ij}^{12}&-m_{ij}^{22}\end{pmatrix}$ (59)

The expression of $m_{ij}^{11}$, $m_{ij}^{12}$ and $m_{ij}^{22}$ are the same
as before (cf Eq. (31)). The off diagonal coefficient is given by

$m_{ij}^{13}=-t\gamma_{ij}\theta_{ij}$ (60)

Where $\gamma_{ij}$ is a d-wave coefficient such that
$\gamma_{i,i\pm\delta_{y}}=-\gamma_{i\pm\delta_{x},i}=1$ and
$\delta_{x}$/$\delta_{y}$ is the lattice constant along x/y axis. The
parameter $\theta$ is given by

$\theta_{ij}=\langle c_{i\sigma}c_{i\bar{\sigma}}n_{j\sigma}\rangle$ (61)

$\theta_{ij}$ can be expressed in several symmetry channels. For the sake of
giving an intuition of this, we apply Wick theorem on $\theta_{ij}$ (this
cannot be done since Wick theorem is only valid for weak correlations but it
will give an insight of the physics)

$\displaystyle\theta_{ij}=$ $\displaystyle\langle
c_{i\sigma}c_{i\bar{\sigma}}\rangle n_{j\sigma}-\langle
c_{i\sigma}c^{\dagger}_{j\sigma}\rangle c_{i\bar{\sigma}}c_{j\sigma}$ (62)
$\displaystyle+\langle c_{i\sigma}c_{j\sigma}\rangle
c_{i\bar{\sigma}}c^{\dagger}_{j\sigma}-\langle
c_{i\bar{\sigma}}c_{j\sigma}\rangle c_{i\sigma}c^{\dagger}_{j\sigma}$

Since we impose d-wave symmetry by mean of the d-wave $\gamma_{ij}$ factor in
Eq.(60), all the symmetry channels except the one associated to the d-wave
symmetry $\langle c_{i\bar{\sigma}}c_{j\sigma}\rangle$ will be assumed to be
negligible. $\theta_{ij}$ can therefore been seen as an anomalous d-wave
superconductivity mean field parameter. Since we consider singlet pairing, we
have

$\langle c_{i\sigma}c_{i\bar{\sigma}}n_{j\sigma}\rangle=\langle
c^{\dagger}_{i\bar{\sigma}}c^{\dagger}_{i\sigma}n_{j\sigma}\rangle$ (63)

We can still apply translationnal invariance to treat p, n and e as a
constant. We can do the same for $\theta$, but because of d-wave symmetry
$\theta$ is not the same along the x and y axis ($\theta_{x}=-\theta_{y}$). We
will therefore only average on one axis in order to have a non zero $\theta$.
The self-consistent equations remain the same for n and e since they are only
one body correlations. However extending the basis changes the self-consistent
equations of $p_{ij}$ and $\theta_{ij}$. We obtain (cf appendix F for details)

$\begin{cases}p&=\frac{n^{2}}{4}-\frac{\rho_{1}+\phi\rho_{2}}{1-\phi^{2}}-\frac{\rho_{1}-\rho_{2}}{1-\phi}-\frac{\rho_{3}}{1+\phi}\\\
\theta&=\frac{\zeta}{1+\phi}\end{cases}$ (64)

With

$\begin{cases}\phi&=\frac{2}{2-n}(C^{11}_{0}+C^{12}_{0})+\frac{2}{n}(C^{12}_{0}+C^{22}_{0})\\\
\rho_{1}&=\frac{2}{2-n}(C^{11}+C^{12})^{2}+\frac{2}{n}(C^{22}+C^{12})^{2}\\\
\rho_{2}&=\frac{2}{2-n}(C^{13}+C^{14})^{2}+\frac{2}{n}(C^{23}+C^{24})^{2}\\\
\rho_{3}&=\frac{4}{n(2-n)}(C^{11}+C^{12})(C^{22}+C^{12})\\\
\zeta&=\frac{2}{2-n}(C^{11}+C^{12})(C^{13}+C^{14})\\\
&+\frac{2}{n}(C^{12}+C^{22})(C^{23}+C^{24})\end{cases}$ (65)

Let us note that this decoupling is not unique. Several choices can be made.
These choices give similar results but tend to overestimate or underestimate
some quantities, depending on which regime we consider [4]. Except for the
larger M and I matrices and these changes in the self-consistent equation,
everything else remains the same. The expression of $\bm{\kappa}$ from Eq.
(27) remains unchanged, but will involve $4\times 4$ matrices.

####  III.3.2 Results

On Fig. 10 we plot the bands for the Roth and COM2 solutions. There is a
doubling of the bands due to the particle-hole symmetry of the basis: we have
four distincts eigenvalues $\epsilon^{l}$ of the E matrix verifying the
property $\epsilon^{1}=-\epsilon^{3}$ and $\epsilon^{2}=-\epsilon^{4}$. Beside
this doubling, the bands are almost unmodified compared to what we have
without superconductivity. Only one difference can be seen: a gap opening at
$(\pi,0)$. We performed a zoom around zero energy in order to see the gap
better on the insets. The presence of the gap also appears on the Fermi
surface: there is less weight near the ($\pi,0$) compared to what we had in
Fig. 3 without superconductivity.

On Fig. 11a, we plot the parameter $\theta$ as a function of the electron
density n for the COM2 and the Roth solutions (COM1 solution has a vanishing
$\theta$ for every electron density n). The dashed-line corresponds to the
usual d-wave superconducting order parameter $\Delta^{d}_{ij}=\langle
c_{i\sigma}c_{j\bar{\sigma}}\rangle$. We can recover it directly from the
correlation functions involving nearest-neighbors
$C_{ij}^{nm}=\langle\psi^{n}_{i}(\psi^{m}_{j})^{\dagger}\rangle$ using the
following equation

$\Delta_{ij}^{d}=C_{ij}^{13}+C_{ij}^{14}+C_{ij}^{23}+C_{ij}^{24}$ (66)

The maximum of $\theta$ and $\Delta^{d}$ are at the same electron density. For
the Roth solution this corresponds to n=0.8, while for the COM2 solution it is
around n=0.9. We already showed that n=0.8 corresponds to the Van Hove
singularity for the Roth solution in the discussion of Fig. 5 and 6. This is
in agreement with other studies [15] [10]. We claim the same phenomenon occurs
for the COM2 solution. In Fig. 10 we plotted the bands and Fermi surfaces for
the COM2 solution at n=0.9 where superconductivity is maximum. Beside the gap
opening, the band for COM2 exhibits some flatness at $(\pi,0)$ (it is at least
flatter than the Roth solution) and its Fermi Surface is almost diamond-like.
We justified this claim by plotting the density of state at the Fermi energy
as a function of electron density on Fig. 11b. Let us note we did not
considered superconductivity to compute this density of state (in order to see
the Van Hove singularity): this is why we do not have any superconducting gap.
The density of states was computed using Eq. (50) at $\omega=0$. The maximum
of the density of states can be seen at n=0.9 for the COM2 solution. For the
Roth solution a maximum can be found near half-filling but this maximum is
associated with the extremum of the band and does not improve
superconductivity because it is too close to half-filling. Another maximum is
found at n=0.8 and corresponds to the flatness at $(\pi,0)$ we already
associated this with the Van Hove singularity. This is in agreement with our
physical intuition: a high density of states means there are a lot of
available electrons available to form Cooper pairs.

We checked that the $\theta$ parameter also verifies the following particle-
hole symmetry

$\theta\rightarrow\theta^{*}$ (67)

In order to satisfy particle-hole symmetry, there must be another maximum of
$\theta$, therefore another Van Hove singularity in the electron doped regime.
On Fig. 11b, we indeed see another peak both for superconductivity and the
density of states in the $n>1$ area. They correspond to the particle-hole
symmetric of the peaks in the hole doped region.

In this method, the gap opening observed on the bands in Fig. 10 is of the
order of $\Delta^{d}$, as it is expected. The value of $\theta$ affects both
superconductivity and the density n, since $\theta$ involves both quantities.

Figure 12: Square lattice with the considered hopping. We will consider long
ranged hoppings up to $t_{5}$. The associated values of these hoppings can be
found in Table. 1

((a)) tb1: ARPES fit Bi2212

((b)) tb2: fit of bonding surface of Bi2212

((c)) tb 3: Modified tb2 to get better incommensurability

((d)) tb 4: underdoped LSCO

Figure 13: Bands renormalization using Roth minimization with hoppings up to
order 5. The tight bindings parameters are taken from Ref [20] and are given
in table I. The red line is the non-interactive, tight binding dispersion. We
are at 20% hole doping.

##  IV Further hoppings study

We now want to consider the effect of higher hopping terms in order to get
closer to more realistic materials. We want to see if the results we had with
nearest-neighbours hoppings on the bands, the validity of the Luttinger
theorem and superconductivity are modified by further hoppings. Including at
least next nearest neighbors in the model is enough to break the particle-hole
symmetry. We will only consider the Roth solution in this section, since COM2
and COM1 solutions have been studied with next-nearest neighbors in Ref. [13],
and Roth presents Fermi surface closer to what is observed by ARPES for
cuprates.

###  IV.1 Effect of further nearest neighbours and bands analysis

From now on we will consider four different sets of tight-binding parameters
all corresponding to a square lattice as in Fig. 12.

These tight-binding parameters correspond to cuprates Bi2212 and LSCO, which
are strongly correlated. Their values, taken from Ref. [20],are given in Table
1. They are such that the energy for a tight binding model of a square lattice
is given by

$\displaystyle\epsilon_{tb}(k)$ $\displaystyle=2t(cos(k_{x})+cos(k_{y}))$ (68)
$\displaystyle+4t_{2}cos(k_{x})cos(k_{y})$
$\displaystyle+2t_{3}(cos(2k_{x})+cos(2k_{y}))$
$\displaystyle+4t_{4}(cos(k_{x})cos(2k_{x})+cos(2k_{x})cos(k_{y}))$
$\displaystyle+4t_{5}cos(2k_{x})cos(2k_{y})$

In the following we will normalize every plots so we have t=1 (we will divide
every tight binding parameter by t in absolute value).

| t | $t_{2}$ | $t_{3}$ | $t_{4}$ | $t_{5}$
---|---|---|---|---|---
tb1 | -0.2956 | 0.0818 | -0.0260 | -0.0280 | 0.0255
tb2 | -0.3399 | 0.1184 | -0.0397 | 0.0086 | 0.0006
tb3 | -0.2941 | 0.0731 | 0.0048 | -0.0325 | 0.0035
tb4 | -0.3912 | 0.0370 | -0.0294 | -0.0350 | -0.0087

Table 1: Values of the 4 tight bindings we are going to consider. tb1 is based
on an ARPES fit of Bi2212. tb2 corresponds to the bonding surface of Bi2212,
tb3 is a modified version of tb2 to get a flatter band and tb4 corresponds to
underdoped LSCO.

In this section we include hoppings up to $t_{5}$ and study the bands and
Fermi surface behavior. With additional hopping terms, only the expression of
the matrix M changes. In addition to the $\alpha_{il}^{1}$ parameter appearing
in Eq. (31) higher hopping terms will appear in the M matrix. It becomes

$\displaystyle m_{ij}^{11}=$
$\displaystyle-\mu(1-\frac{n_{i}}{2})\delta_{ij}-\sum\limits_{k=1}^{N}t\left[\alpha_{ij}^{k}(1-n_{i}+p_{ij})+\delta_{ij}\sum\limits_{l}\alpha^{k}_{il}e_{il}\right]$
(69) $\displaystyle m_{ij}^{12}=$
$\displaystyle\sum\limits_{k=1}^{N}t\left[\alpha_{ij}^{k}(\frac{n_{i}}{2}-p_{ij})-\delta_{ij}\sum\limits_{l}\alpha_{il}^{k}e_{il}\right]$
$\displaystyle m_{ij}^{22}=$
$\displaystyle-(\mu-U)\frac{n_{i}}{2}\delta_{ij}+\sum\limits_{k=1}^{N}t\left[\alpha_{ij}^{k}p_{ij}-\delta_{ij}\sum\limits_{l}\alpha_{il}^{k}e_{il}\right]$

where

$\begin{cases}\alpha^{N}_{il}&=1\ \ \ \ \text{if i and l are
}\overbrace{\text{next-...}}^{N-1\text{ times}}\text{-nearest neighbour}\\\
\alpha^{N}_{il}&=0\ \ \ \ \text{Otherwise}\end{cases}$ (70)

Each new hopping considered adds a term in the tight-binding Hamiltonian which
is then added in the M matrix. The parameters p and e depend on i-j, so we
should make a distinction between $e^{1}_{ij}$ with i and j nearest neighbors
(NN), $e^{2}_{ij}$ with i and j next nearest neighbors (NNN)… $e^{n}_{ij}$ and
$p^{n}_{ij}$ will be associated with their corresponding hopping as in Fig.
12. Translationnal invariance (cf Eq. (33)) still allows us to treat $e^{1}$,
$e^{2}$,…,$e^{5}$, $p^{1}$, $p^{2}$,…,$p^{5}$ as constants. Correlation
functions $\textbf{C}_{ij}=\langle\psi_{i};\psi^{\dagger}_{j}\rangle$ are at
different sites too, so we will also have to make a distinction in the self-
consistent equations between $\textbf{C}^{1}$ for NN, $\textbf{C}^{2}$ for NNN
and so on.

Figure 14: Fermi Surfaces renormalized by the composite operator methods using
Roth minimization with the parameters of Ref [20]. We are at 20% hole doping.
Top left: tb1, top right: tb2, bottom left: tb3, bottom right: tb4.

Figure 15: Non interacting (NI) Fermi Surfaces for the parameters of Ref [20]
at 20% hole doping. Top left: tb1, top right: tb2, bottom left: tb3, bottom
right: tb4.

Fig. 13 presents the bands we obtain for Roth solutions for the four sets of
tight-binding parameters in Table. (1). In Fig. 14 and 15 we plotted
respectively the Fermi surfaces obtained with the method and the Fermi
surfaces of the non interactive tight-binding dispersions (corresponding to
Eq. (68)).

((a)) Chemical potential

((b)) Parameters e

((c)) Parameters p

((d)) Luttinger Violation

Figure 17: (a), (b) and (c) :Parameters as a function of doping for each tight
binding with next nearest neighbors. The circles corresponds to p1 and e1 (NN)
while the squares corresponds to p2 and e2 (NNN). (d) : Luttinger violation
for each tight binding.

In the hole doped regime, the Fermi surface we obtained from the Roth solution
has the same general shape as the non interacting Fermi Surface. The composite
operators method produces Fermi surfaces that appear to be smaller/larger than
the tight-binding ones. This is in agreement with the violation of the
Luttinger theorem observed with nearest-neighbour hopping, and indicates that
it is still violated with further hoppings. We checked the opposite situation
is accordingly observed in electron doped area: the Fermi Surface obtained
with the method is at a lower doping than the non interacting one. In the
following we study the particle-hole symmetry and the Luttinger theorem
violation with next-nearest neighbours.

((a)) $\theta$ and $\Delta^{d}$

((b)) Density of states at the Fermi energy

((c)) Bands for tb3 at n=0.6 with no superconductivity

((d)) Bands for tb3 at n=0.6 with superconductivity

Figure 18: (a): Anomalous superconducting order parameter $\theta$ as a
function of doping for the four sets of tight bindings parameters with Roth
minimization. The dashed line corresponds to the usual d-wave superconducting
order parameter we rebuilt from the correlation function. We are at U=8t. (b):
Density of state (with no superconductivity in order to see the peaks) at the
Fermi energy. The 2 peaks corresponds to the 2 Van Hove singularity in hole
and electron doping respectively. To illustrate this the bottom plots (c) and
(d) are the bands at the Van Hove singularity with tight binding parameters
tb3 at n=0.6 respectively with: (c) no superconductivity, where we clearly see
the flat band associated to Van Hove singularity (d) with superconductivity,
where we see on the zoom in the inset the gap of the order of $\Delta^{d}$.

###  IV.2 Particle-hole symmetry and Luttinger theorem with further neighbours

For the sake of simplicity we restrict ourselves to next-nearest neighbors in
this part. In order to keep the full Hamiltonian invariant under particle-hole
symmetry the following relation is needed

$t_{2}\rightarrow-t_{2}$ (71)

This is because the hopping term of the Hamiltonian
$\sum\limits_{ij}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}$ transforms into
$\sum\limits_{ij}t_{ij}(-1)^{i+j}c_{i\sigma}c^{\dagger}_{j\sigma}$. For i and
j next-nearest neighbour we have $(-1)^{i+j}=1$. We obtain
$-\sum\limits_{\langle\langle
i,j\rangle\rangle}t_{2}c^{\dagger}_{j\sigma}c_{i\sigma}$, which has an
additional minus sign, breaking the symmetry unless $t_{2}$ changes its sign
with the transformation.

Hence, we know that at least next-nearest neighbors will break particle-hole
symmetry (because $t_{2}$ is an external parameter that will always keep the
same sign, therefore we do not verify Eq. (71)). The particle-hole
transformation works the same for n and $\mu$

$\displaystyle n$ $\displaystyle\rightarrow 2-n$ (72) $\displaystyle\mu$
$\displaystyle\rightarrow U-\mu$

The main difference with the nearest-neighbors case comes from the necessity
to differentiate the e and p parameters depending on the considered hopping.
We already know the transformation for $p^{1}$ and $e^{1}$. The only
difference for $p^{2}$ and $e^{2}$ is $(-1)^{i+j}=1$. $p^{2}$ will not be
affected by the $(-1)^{i+j}$ coming from the transformation because it
involves only two bodies operators. Thus only $e^{2}$ has an additional minus
sign under the transformation, and we obtain the following relations

$\displaystyle p^{1}$ $\displaystyle\rightarrow 1-n+p^{1}$ (73) $\displaystyle
p^{2}$ $\displaystyle\rightarrow 1-n+p^{2}$ $\displaystyle e^{1}$
$\displaystyle\rightarrow-e^{1}$ $\displaystyle e^{2}$
$\displaystyle\rightarrow e^{2}$

In Figs. 17a, 17b and 17c we plot the parameters $\mu$, e and p as a function
of doping for the four considered tight-binding parameters of Table. (1), only
considering t and $t_{2}$. Parameters $e^{2}$ and $p^{2}$ indeed break the
particle-hole symmetry from Eq. (73). The chemical potential, as well as
$e^{1}$ and $p^{1}$ behave as in the nearest neighbor case. We also checked
that Eq. (73) is respected if we impose $t_{2}\rightarrow-t_{2}$ when
$n\rightarrow 2-n$.

Finally we can again study the Luttinger theorem. In Fig. 17d, the area of the
Fermi surface is plotted as a function of electron density. Interestingly we
observe an analogous behavior as in the nearest-neighbor case. The Luttinger
violation does not seem to be modified by next-nearest neighbors hoppings and
is barely affected when we modify tight-binding parameters. As already
mentioned before, the Luttinger theorem is strongly violated around half
filling and is recovered further away. This confirms why in the previous
section the Fermi surfaces seemed to be at a higher doping than the electron
density we considered when we are in the hole doped regime (and conversely, at
a lower doping in the electron doped regime).

###  IV.3 Superconductivity, Van Hove singularity and density of state

Superconductivity can also be included in the model with more hoppings. In
this section we include hoppings up until $t_{5}$. We assume again a d-wave
symmetry and only consider nearest-neighbor pairing. As before, the M and I
matrices become 4x4 and have the same symmetries as the nearest neighbours
case. The I matrix is independent of t and t’ and is thus the same as before.
The main difference is that the coefficients $m_{ij}^{11}$, $m_{ij}^{12}$ and
$m_{ij}^{22}$ are now given by Eq. (69). $m_{ij}^{13}$ stays identical to its
expression in the nearest neighbor case since we only consider
superconductivity for nearest neighbors. Therefore $m_{ij}^{13}$ is
proportionnal to the $\theta$ parameter, defined the same way as before
($\theta_{ij}=\langle c_{i\uparrow}c_{i\downarrow}n_{j\sigma}\rangle$).

In Fig. 18a, the $\theta$ and $\Delta^{d}$ parameters are plotted as a
function of electron density for the four tight-binding parameters with the
Roth solution. Since we have included further neighbors, $\theta$ is not
particle-hole symmetric anymore. In Fig. 18b, we plotted the density of states
at the Fermi energy without superconductivity using Eq. (50) at $\omega=0$ (in
order to see the peaks with no superconducting gap). In Fig. 18c and 18d, we
plotted without and with superconductivity the bands at electron densities
corresponding to the maximum of the density of states for tb3 when the system
is hole doped ($n\approx 0.6$). The bands with no superconductivity on 18c are
flats at the Fermi energy. This proves that the maximum of Fig. 18b correspond
to the Van Hove singularity. On 18d we see in the inset the gap at $(\pi,0)$
is again of the order of $2\Delta^{d}$. The maximum of the $\theta$ and
$\Delta^{d}$ parameters are thus at the same electron density as the Van Hove
singularity. Hence, the situation is the same as in the nearest-neighbor case.
In Table 2 we give the electron densities associated to the Van Hove
singularities for the four sets of tight-binding parameters. Let us note there
is no proportionality between the peak in the density of state and
$\Delta^{d}$: it is different for every set of tight-binding parameters. This
is seemingly a consequence of the electronic correlations which are treated
differently for each tight-binding parameters as a consequence of the main
approximation of the method Eq. (23b).

| tb1 | tb2 | tb3 | tb4
---|---|---|---|---
Van Hove (electron) | 1.04 | 1.05 | 1.10 | 1.11
Van Hove (hole) | 0.57 | 0.45 | 0.6 | 0.58

Table 2: Electron density at which there is a Van Hove singularity at the
chemical potential. At this values a flat band lies at the Fermi energy and
the associated Fermi surface is diamond like. These singularities occur both
in electron and in hole.

The most striking feature is that superconductivity in the method seems to be
induced by the proximity of the Van Hove singularity. This was already the
case with nearest neighbors but this property seems unaffected by further
hopping terms. For cuprates there exist a consensus that antiferromagnetism is
the interaction necessary to explain the pairing mechanism [21]. This has also
the advantage to explain why experimentally superconductivity is observed
around 15% hole doping. Since this method predicts superconductivity only
close to the Van Hove singularity, it is non zero at some doping values which
does not correspond to what is expected. This flaw is maybe a consequence of
the main hypothesis of the method (detailed in Appendix A): it is maybe
necessary to consider dynamical corrections to the self-energy in order to
observe a different behavior for superconductivity.

##  V Conclusion

Considering only nearest neighbors first, three solutions have been isolated
and studied: COM1, COM2 and Roth solutions. The Roth solution is the unique
solution obtained if the self-consistency is performed on a parameter we
called p and which enforced charge-charge, spin-spin and pair-pair
correlations. This solution violates the Pauli principle because the
correlation function
$C^{12}_{0}=\langle\xi_{i\uparrow}\eta_{i\uparrow}^{\dagger}\rangle$ obtained
after convergence is numerically non zero, even though it should be
analytically because of the Pauli principle. COM1 and COM2 are two independant
solutions obtained by performing a self-consistency imposing the correlation
function $C^{12}_{0}$ instead of the parameter p to enforce the Pauli
principle.

We have performed a systematic comparison of the bands and Fermi surfaces for
these three solutions. The three solutions present two Hubbard bands
associated to the two eigenvalues of the method. COM1 is not renormalized by
the interactions. COM2 has two hole pockets whatever the doping. Roth only has
one hole pocket but a second small one appears close to half filling at
$(\pi,\pi)$, possibly because antiferromagnetism has been neglected. Both Roth
and COM2 solutions are insulators at half filling because the chemical
potential lies in between the two Hubbard bands. The density of states for the
Roth solution does not present the usual quasiparticle peak at $\omega=0$ at
half-filling and is in agreement with Determinantal quantum Monte-Carlo
simulations.

The particle-hole symmetry is respected with nearest-neighbour hoppings: we
checked that all solutions verify particle-hole relations, including the Roth
solution. This result is different from what was predicted by ref. [9]: a
solution of this method can violate the Pauli principle and still be particle-
hole symmetric. All solutions break the Luttinger close to half filling. This
results is also in agreement with the determinantal quantum Monte-Carlo
simulation [18]. The Luttinger theorem is only proven for weakly interacting
systems. We think the Roth solution is the most physical. Although there is a
violation of the Pauli principle, this solution exhibits Fermi surfaces
typical of strongly correlated materials such as cuprates, contrary to COM1
and COM2. The COM1 solution consists of two tight-binding dispersions splitted
by U and is not renormalized by interactions, whereas COM2 presents a second
hole pocket which is not observed in ARPES experiments for Bi2212 and LSCO.

We extended the method to study d-wave superconductivity. Four composite
operators are necessary to perform this study. The doubling of the basis leads
to four eigenvalues, which are particle-hole symmetric because of the symmetry
of the basis. We observed a gap opening at $(\pi,0)$ for the COM2 and Roth
solutions, of the order of the usual d-wave superconducting order parameter
$\Delta_{ij}^{d}=\langle c_{i\uparrow}c_{j\downarrow}\rangle$. We only observe
non zero superconductivity close to the Van Hove singularity. While this has
already been predicted for the Roth solution [15], we showed that the COM2
solution behaves the same. We observed two Van Hove singularities for the COM2
and Roth solutions: one in the hole-doped and one in the electron-doped
regime. Both singularities are a consequence of a flat band at $(\pi,0)$.
Lastly we noticed $\Delta^{d}=\langle c_{i\uparrow}c_{j\downarrow}\rangle$ is
three or four times bigger than $\theta$. $\theta$ has no proportionnality to
$\Delta^{d}$ because it includes both a superconducting and a charge channel
when we try to decouple it.

We studied the effect of longer ranged tight-binding parameters fitted from
ARPES experiments on cuprates Bi2212 and LSCO from ref.[20], with the Roth
minimization. If hole doped, the Fermi surfaces obtained with the method are
similar to the non interacting Fermi surface, but at a higher doping. The
opposite situation happens if the system is doped in electron: the Fermi
surfaces appear similar but at a lower doping than the non interacting one.
This is in agreement with the violation of the Luttinger theorem we still
observe around the Mott transition. Adding further hoppings does not seem to
modify this property. We checked that adding further hoppings breaks particle-
hole symmetry. Finally, superconductivity also behaves mostly the same with
further neighbors compared to the nearest neighbor case. The main difference
is that Van Hove singularities as well as superconductivity peaks are moved at
different electron density compared to the nearest-neighbors case. While this
result is not at all in agreement with what was observed experimentally [5],
it may be an indication that the main approximation done in this scheme needs
to be refined. The dynamical corrections to the self-energy need to be
considered in order to maybe have a behavior for superconductivity which is
more in agreement with experimental observations.

## References

  * [1] E. Morosan, D. Natelson, A. H. Nevidomskyy, and Q. Si, “Strongly Correlated Materials,” Advanced Materials, vol. 24, pp. 4896–4923, Sept. 2012.  arXiv:1309.4473 [cond-mat].
  * [2] S. Sachdev, “Order and quantum phase transitions in the cuprate superconductors,” Reviews of Modern Physics, vol. 75, pp. 913–932, July 2003.  arXiv:cond-mat/0211005.
  * [3] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of Superconductivity,” Physical Review, vol. 108, pp. 1175–1204, Dec. 1957\.
  * [4] J. G. Bednorz and K. A. Müller, “Perovskite-type oxides—the new approach to high-${T}_{c}$ superconductivity,” Rev. Mod. Phys., vol. 60, pp. 585–600, Jul 1988.
  * [5] M. R. Norman and C. Pepin, “The Electronic Nature of High Temperature Cuprate Superconductors,” Reports on Progress in Physics, vol. 66, pp. 1547–1610, Oct. 2003.  arXiv:cond-mat/0302347.
  * [6] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Reviews of Modern Physics, vol. 68, pp. 13–125, Jan. 1996\.
  * [7] E. H. Lieb and F. Y. Wu, “The one-dimensional Hubbard model: A reminiscence,” Physica A: Statistical Mechanics and its Applications, vol. 321, pp. 1–27, Apr. 2003.  arXiv:cond-mat/0207529.
  * [8] J. Hubbard, “Electron correlations in narrow energy bands,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, vol. 276, pp. 238–257, Nov. 1963.
  * [9] A. Avella, F. Mancini, D. Villani, L. Siurakshina, and V. Y. Yushankhai, “The Hubbard Model in the Two-Pole Approximation,” International Journal of Modern Physics B, vol. 12, pp. 81–97, Jan. 1998.
  * [10] J. Beenen and D. M. Edwards, “Superconductivity in the two-dimensional Hubbard model,” Physical Review B, vol. 52, pp. 13636–13651, Nov. 1995\.
  * [11] L. M. Roth, “Electron Correlation in Narrow Energy Bands. I. The Two-Pole Approximation in a Narrow S Band,” Physical Review, vol. 184, pp. 451–459, Aug. 1969.
  * [12] A. Avella and F. Mancini, “Anomalous Self-Energy Features in the 2D Hubbard Model,” Acta Physica Polonica A, vol. 113, pp. 395–398, Jan. 2008.
  * [13] A. Avella, F. Mancini, D. Villani, and H. Matsumoto, “The two-dimensional t-t’-U model as a minimal model for cuprate materials,” The European Physical Journal B, vol. 20, pp. 303–311, Apr. 2001.
  * [14] E. J. Calegari, A. C. Lausmann, S. G. Magalhaes, C. M. Chaves, and A. Troper, “Pseudogap and the specific heat of high T$_c$ superconductors: a Hubbard model in a n-pole approximation,” Journal of Physics: Conference Series, vol. 592, p. 012075, Mar. 2015.  arXiv:1412.4692 [cond-mat].
  * [15] E. J. Calegari, S. G. Magalhaes, and A. A. Gomes, “Superconductivity in a two dimensional extended Hubbard model,” The European Physical Journal B, vol. 45, pp. 485–496, June 2005.  arXiv:cond-mat/0505427.
  * [16] A. Avella and F. Mancini, “The Composite Operator Method (COM),” vol. 171, pp. 103–141, 2012.  arXiv:1112.1103 [cond-mat].
  * [17] D. Bergeron, D. Chowdhury, M. Punk, S. Sachdev, and A.-M. S. Tremblay, “Breakdown of Fermi liquid behavior at the (\pi,\pi)=2k_f spin-density wave quantum-critical point: the case of electron-doped cuprates,” Physical Review B, vol. 86, p. 155123, Oct. 2012.  arXiv:1207.1106 [cond-mat].
  * [18] I. Osborne, T. Paiva, and N. Trivedi, “Fermi-surface Reconstruction in the Repulsive Fermi-Hubbard Model,” June 2020.  arXiv:2001.07197 [cond-mat].
  * [19] J. M. Luttinger, “Fermi Surface and Some Simple Equilibrium Properties of a System of Interacting Fermions,” Physical Review, vol. 119, pp. 1153–1163, Aug. 1960.
  * [20] M. R. Norman, “Linear response theory and the universal nature of the magnetic excitation spectrum of the cuprates,” Physical Review B, vol. 75, p. 184514, May 2007.
  * [21] L. Taillefer, “Scattering and Pairing in Cuprate Superconductors,” Annual Review of Condensed Matter Physics, vol. 1, pp. 51–70, Aug. 2010.
  * [22] T. D. Stanescu, I. Martin, and P. Phillips, “d x 2 – y 2 pairing of composite excitations in the two-dimensional Hubbard model,” Physical Review B, vol. 62, pp. 4300–4308, Aug. 2000.

## Appendix A : physical meaning of the approximation on the current

In order to be able to compute the M matrix with the full Hubbard Hamiltonian
in the strongly coupled regime, we neglected all terms in the current (Eq.
(8)) that are not along $\psi$. Though afterwards we computed explicitly the
currents, the terms orthogonal to $\psi$ are neglected through the self
consistency because we use the relation $\textbf{M}=\textbf{E}\textbf{I}$,
which is a consequence of our hypothesis.

To better understand the consequences of such approximation we need to
consider the self energy. To do so, we follow the step of Ref. [12]. Let’s
therefore consider the full expression of the current

$J_{i}=\sum\limits_{l}\textbf{E}_{il}\psi_{l}+\delta\phi_{i}$ (74)

Our approximation, $\langle\\{\delta\phi_{i},\psi_{j}^{\dagger}\\}\rangle=0$
allows us to write
$\textbf{E}_{ij}=\sum\limits_{l}\textbf{M}_{il}\cdot\textbf{I}^{-1}_{lj}$.
From this, we know the 0-th order Green’s function is defined for
$\delta\phi=0$ by $\textbf{S}^{0}=\frac{\textbf{I}}{\omega-\textbf{E}}$.

Without this approximation, the full equations of motion for the composite
Green’s function is

$\frac{d}{d\tau}\textbf{S}_{ij}(\tau)=\delta(\tau)\langle\\{\psi_{i}(\tau);\psi^{\dagger}_{j}\\}\rangle+\theta_{H}(\tau)\langle\\{J_{i}(\tau);\psi_{j}^{\dagger}\\}\rangle$
(75)

In Fourier space, and with Eq. (74), this becomes

$\textbf{S}_{ij}(\omega)=\textbf{S}^{0}_{ij}(\omega)+\textbf{I}^{-1}_{ij}\textbf{S}^{0}_{ij}(\omega)\langle\\{\delta\phi_{i}(\omega);\psi^{\dagger}_{j}(\omega)\\}\rangle$
(76)

Now we know that

$\frac{d}{d\tau}\psi_{i}^{\dagger}=\sum\limits_{l}\textbf{E}_{il}\psi_{l}^{\dagger}+\delta\phi_{i}$
(77)

Therefore in Fourier space

$\psi^{\dagger}_{i}(\omega)=(\omega-\textbf{E}_{ij}+i0^{+})^{-1}\delta\phi_{i}=\textbf{S}^{0}_{ij}(\omega)\textbf{I}^{-1}_{ij}\delta\phi_{i}$
(78)

Hence we can inject this in the equations of motion on $\textbf{S}_{ij}$ to
get

$\textbf{S}_{ij}(\omega)=\textbf{S}^{0}_{ij}(\omega)+\textbf{S}^{0}_{ij}(\omega)\textbf{I}_{ij}^{-1}\langle\\{\delta\phi_{i};\delta\phi_{j}^{\dagger}\\}\rangle\textbf{I}_{ij}^{-1}\textbf{S}^{0}_{ij}(\omega)$
(79)

We can introduce the scattering matrix T

$\textbf{T}=\textbf{I}^{-1}\langle\\{\delta\phi;\delta\phi^{\dagger}\\}\rangle\textbf{I}^{-1}$
(80)

We therefore obtain the familiar form

$\textbf{S}=\textbf{S}^{0}+\textbf{S}^{0}\textbf{T}\textbf{S}^{0}$ (81)

We can introduce the self-energy through the relation
$\textbf{T}\textbf{S}^{0}=\textbf{I}^{-1}\bm{\Sigma}\textbf{S}$, and obtain in
reciprocal space the Dyson equation

$\textbf{S}(k,\omega)=\frac{\textbf{I}}{\omega-\textbf{E}(k)-\bm{\Sigma}(k,\omega)+i0^{+}}$
(82)

From this equation we clearly see the consequences of our approximation.
Neglecting $\delta\phi$ in the current, which are all the contributions
orthogonal to $\psi$, is equivalent to neglecting $\bm{\Sigma}(k,\omega)$,
therefore working with a static self energy. This approximation neglects all
dynamical dependencies of the self energy.

## Appendix B : Spectral representation and residue theorem of the composite
Green’s function

We start from the following expression for the composite Green’s function

$\textbf{S}_{ij}(\omega)=\frac{\textbf{I}_{ij}}{\omega-\textbf{E}_{ij}+i0^{+}}$
(83)

E is diagonalizable, its eigenvalues are called $\epsilon^{1}$ and
$\epsilon^{2}$, and we can rewrite this as

$\displaystyle\textbf{S}_{ij}(\omega)=$ $\displaystyle\frac{\textbf{I}_{ij}\
Com(\omega-\textbf{E}_{ij}+i0^{+})^{T}}{det(\omega-\textbf{E}_{ij}+i0^{+})}$
(84) $\displaystyle=$ $\displaystyle\frac{\textbf{I}_{ij}\
Com(\omega-\textbf{E}_{ij}+i0^{+})^{T}}{(\omega-\epsilon^{1}_{ij}+i0^{+})(\omega-\epsilon^{2}_{ij}+i0^{+})}$

Where Com(A) is the cofactors matrix of A. We now want to apply residue
theorem.

We consider the function $\frac{\textbf{S}_{ij}(z)}{z-\omega+i\eta}$ with
$\eta$ small. Since we know any Green’s functions, including S, is analytical
on the upper half-circle of the complex plane we integrate over this contour
$C_{u}$. The poles of S are real and we have another pole at $z=\omega-i\eta$,
which is on the lower half-circle of the complex plane. Therefore no poles
lies in $C_{u}$ and we have

$\oint\limits_{C_{u}}\frac{\textbf{S}_{ij}(z)}{z-\omega-i\eta}dz=0$ (85)

We now use the previous expression of S to get

$\oint\limits_{C_{u}}\frac{\textbf{I}_{ij}\
Com(z-\textbf{E}_{ij})^{T}}{det(z-\textbf{E})(z-\omega-i\eta)}dz=0$ (86)

Replacing the integrals by the sum of all the residues, we obtain

$\sum\limits_{z_{0}\in\mathcal{P}}Res\left(\frac{\textbf{I}_{ij}\
Com(z-\textbf{E}_{ij})^{T}}{(z-\epsilon^{1}_{ij})(z-\epsilon^{2}_{ij})(z-\omega-i\eta)},z\rightarrow
z_{0}\right)=0$ (87)

In this equation, $\mathcal{P}$ denotes the poles of
$\frac{\textbf{S}_{ij}(z)}{z-\omega+i\eta}$. It has three poles

$\displaystyle z_{1}$ $\displaystyle=\omega+i\eta$ (88) $\displaystyle z_{2}$
$\displaystyle=\epsilon^{1}_{ij}$ $\displaystyle z_{3}$
$\displaystyle=\epsilon^{2}_{ij}$

These poles are all non degerated, meaning the residue can easily be computed
using

$Res(f(z),z\rightarrow z_{0})=\lim\limits_{z\rightarrow z_{0}}(z-z_{0})f(z)$
(89)

Doing so leads to

$\displaystyle 0=$ $\displaystyle\frac{\textbf{I}\
Com(\omega-\textbf{E}_{ij}+i\eta)^{T}}{(\omega-\epsilon^{1}_{ij}+i\eta)(\omega-\epsilon^{2}_{ij}+i\eta)}$
(90) $\displaystyle+$ $\displaystyle\frac{\textbf{I}_{ij}\
Com(\epsilon^{1}_{ij}-\textbf{E}_{ij})^{T}}{(\epsilon^{1}_{ij}-\epsilon^{2}_{ij})(\epsilon^{1}_{ij}-\omega-i\eta)}$
$\displaystyle+$ $\displaystyle\frac{\textbf{I}_{ij}\
Com(\epsilon^{2}_{ij}-\textbf{E}_{ij})^{T}}{(\epsilon^{2}_{ij}-\epsilon^{1}_{ij})(\epsilon^{2}_{ij}-\omega-i\eta)}$

We recognize the first term can be rewritten as

$\frac{\textbf{I}_{ij}\
Com(\omega-\textbf{E}_{ij}+i\eta)^{T}}{det(\omega-\textbf{E}_{ij}+i\eta)}=\textbf{I}_{ij}(\omega-\textbf{E}_{ij}+i\eta)^{-1}=\textbf{S}_{ij}(\omega+i\eta)$
(91)

Finally we rearrange the equation (and evaluate S in $\omega-i\eta$ to get
retarded composite Green’s functions) to recover the desired decomposition
(replacing $\eta$ by $0^{+}$)

$\displaystyle\textbf{S}_{ij}(\omega)=$ $\displaystyle\frac{\textbf{I}_{ij}\
Com(\epsilon^{2}_{ij}-\textbf{E}_{ij})^{T}}{(\epsilon^{1}_{ij}-\epsilon^{2}_{ij})(\omega-\epsilon^{2}_{ij}+i0^{+})}$
(92) $\displaystyle-$ $\displaystyle\frac{\textbf{I}_{ij}\
Com(\epsilon^{1}_{ij}-\textbf{E}_{ij})^{T}}{(\epsilon^{1}_{ij}-\epsilon^{2}_{ij})(\omega-\epsilon^{1}_{ij}+i0^{+})}$

We pose

$\bm{\kappa}^{a}_{ij}=(-1)^{a+1}\frac{\textbf{I}_{ij}\
Com(\epsilon^{a}_{ij}-\textbf{E}_{ij})^{T}}{(\epsilon^{1}_{ij}-\epsilon^{2}_{ij})}$
(93)

With $a\in\\{1,2\\}$ to obtain the form given in Eq. (26).

## Appendix C : Computations of M and I matrices

In this section we first derive from the composite operator algebra the
currents, then we use them to obtain the M and I matrices. Let us start by
writing few of the most useful commutation relation that we used for the
computations

$\displaystyle\\{\eta_{i\sigma};\eta_{j\sigma^{\prime}}^{\dagger}\\}=$
$\displaystyle\delta_{ij}(\delta_{\sigma\sigma^{\prime}}n_{i\bar{\sigma}}-\delta_{\sigma\bar{\sigma^{\prime}}}c^{\dagger}_{i\sigma}c_{i\bar{\sigma}})$
(94) $\displaystyle\\{\xi_{i\sigma};\xi_{j\sigma^{\prime}}^{\dagger}\\}=$
$\displaystyle\delta_{ij}(\delta_{\sigma\sigma^{\prime}}(1-n_{i\bar{\sigma}})+\delta_{\sigma\bar{\sigma^{\prime}}}c^{\dagger}_{i\sigma}c_{i\bar{\sigma}})$
$\displaystyle\\{\xi_{i\sigma};\eta_{j\sigma^{\prime}}^{\dagger}\\}=$
$\displaystyle 0$
$\displaystyle\\{c_{i\sigma};\xi_{j\sigma^{\prime}}^{\dagger}\\}=$
$\displaystyle\delta_{ij}(\delta_{\sigma\sigma^{\prime}}(1-n_{j\sigma})+\delta_{\sigma\bar{\sigma^{\prime}}}c^{\dagger}_{i\sigma^{\prime}}c_{i\bar{\sigma^{\prime}}})$
$\displaystyle\\{c_{i\sigma};\eta_{j\sigma^{\prime}}^{\dagger}\\}=$
$\displaystyle\delta_{ij}(\delta_{\sigma\sigma^{\prime}}n_{i\bar{\sigma}}-\delta_{\sigma\bar{\sigma^{\prime}}}c^{\dagger}_{i\sigma^{\prime}}c_{i\bar{\sigma^{\prime}}})$
$\displaystyle\\{c_{i\sigma};\xi_{j\sigma^{\prime}}\\}=$
$\displaystyle\delta_{ij}\delta_{\sigma\bar{\sigma^{\prime}}}c_{i\sigma^{\prime}}c_{i\bar{\sigma^{\prime}}}$
$\displaystyle\\{c_{i\sigma};\eta_{j\sigma^{\prime}}\\}=$
$\displaystyle-\delta_{ij}\delta_{\sigma\bar{\sigma^{\prime}}}c_{i\sigma^{\prime}}c_{i\bar{\sigma^{\prime}}}$

From this we can explicitly compute the commutators of the composite operators
$\psi$ with the hamiltonian to get the following currents

$\displaystyle j_{i}^{1}=$
$\displaystyle-\mu\xi_{i\sigma}-\sum\limits_{l}t_{il}\left(c_{l\sigma}-n_{i\bar{\sigma}}c_{l\sigma}+S_{i}^{-}c_{l\bar{\sigma}}-\Delta_{i}c^{\dagger}_{l\bar{\sigma}}\right)$
(95) $\displaystyle j_{i}^{2}=$
$\displaystyle-(\mu-U)\eta_{i\sigma}+\sum\limits_{l}t_{il}\left(-n_{i\bar{\sigma}}c_{l\sigma}+S_{i}^{-}c_{l\bar{\sigma}}-\Delta_{i}c^{\dagger}_{l\bar{\sigma}}\right)$

The $\textbf{I}=\langle\\{\psi;\psi^{\dagger}\\}\rangle$ and
$\textbf{M}=\langle\\{j;\psi^{\dagger}\\}\rangle$ matrices can be explicitly
computed from these expressions. Note that

$\displaystyle\textbf{I}_{ij}^{12}=$ $\displaystyle\textbf{I}_{ij}^{21}$ (96)
$\displaystyle\textbf{M}_{ij}^{12}=$ $\displaystyle\textbf{M}_{ij}^{12}$

## Appendix D : Roth decoupling and computation of p

In this appendix we derive the self-consistent equation of

$p(i-j)=\langle n_{i\uparrow}n_{j\uparrow}\rangle+\langle
S_{i}^{+}S_{j}^{-}\rangle-\langle\Delta_{i}\Delta_{j}^{*}\rangle$ (97)

We will not detail the main decoupling scheme as it is already done in the
main text but rather focus on the algebraic computations and give the form of
the equations for the charge-charge and spin-spin term.

###  V.1 Pair-pair term

Since the general method has already been detailed in the text for this term
we start directly from eq. (43) and compute
$f_{kjl}^{n}=\langle\\{\psi^{n}_{k};c_{j\downarrow}\Delta_{l}^{*}\\}\rangle$.
Computing the anticommutators using the algebraic relations given in appendix
4 leads to

$\displaystyle f^{1}_{kjl}=$
$\displaystyle\delta_{kj}\langle\Delta_{k}\Delta_{l}^{*}\rangle+\delta_{kl}(C^{21}_{kj}+C^{22}_{kj})$
(98) $\displaystyle f^{2}_{kjl}=$
$\displaystyle-\delta_{kj}\langle\Delta_{k}\Delta_{l}^{*}\rangle+\delta_{kl}(C^{11}_{kj}+C^{12}_{kj})$

Eq. (43) then becomes

$\displaystyle\langle c_{i\uparrow}c_{j\downarrow}\Delta_{l}^{*}\rangle$
$\displaystyle=\frac{2}{2-n}\sum\limits_{k}(C^{11}_{ik}+C^{12}_{ik})$ (99)
$\displaystyle[\delta_{kj}\langle\Delta_{k}\Delta_{l}^{*}\rangle+\delta_{kl}(C^{21}_{kj}+C^{22}_{kj})]$
$\displaystyle+\frac{2}{n}\sum\limits_{k}(C^{12}_{ik}+C^{22}_{ik})$
$\displaystyle[\delta_{kj}\langle\Delta_{k}\Delta_{l}^{*}\rangle+\delta_{kl}(C^{21}_{kj}+C^{22}_{kj})]$

Simplifying this equation and setting i=j leads to

$\langle\Delta_{i}\Delta_{l}^{*}\rangle=\frac{4}{n(2-n)}\frac{(C^{11}_{il}+C^{12}_{il})(C^{22}_{il}+C^{21}_{il})}{1-\frac{2}{2-n}(C^{11}_{0}+C^{12}_{0})+\frac{2}{n}(C^{21}_{0}+C^{22}_{0})}$
(100)

Which is the form in the main text. We pose

$\phi=-\frac{2}{2-n}(C^{11}_{0}+C^{12}_{0})+\frac{2}{n}(C^{21}_{0}+C^{22}_{0})$
(101)

Replacing $C^{11}_{0}$ and $C^{22}_{0}$ by their definitions allows us to
express these correlations function explicitly as a function of n. We don’t
explicit $C^{12}_{0}$ however, else it will be zero while it is not
numerically: $C^{12}_{0}$ should stay in the numerical minimization process to
obtain our results. Doing so leads to the following expression for $\phi$

$\phi=\frac{n^{2}-4(\frac{n}{2}-\langle
n_{i\uparrow}n_{i\downarrow}\rangle-C^{12}_{0})}{n(2-n)}$ (102)

With our notations the pair-pair term becomes

$\langle\Delta_{i}\Delta_{l}^{*}\rangle=\frac{\rho_{3}}{1+\phi}$ (103)

###  V.2 Spin-Spin term

The spin-spin term is defined as $\langle S_{i}^{+}S_{l}^{-}\rangle=\langle
c^{\dagger}_{i\uparrow}c_{i\downarrow}c^{\dagger}_{l\downarrow}c_{l\uparrow}\rangle$.
In order to have our basis element as the first term, we rewrite it as

$\langle S_{i}^{+}S_{l}^{-}\rangle=-\langle
c^{\dagger}_{i\uparrow}c_{i\downarrow}c_{l\uparrow}c^{\dagger}_{l\downarrow}\rangle=-\langle
c_{l\uparrow}c^{\dagger}_{l\downarrow}c^{\dagger}_{i\uparrow}c_{i\downarrow}\rangle$
(104)

We therefore introduce the following Green’s functions (setting
$\tau^{\prime}=0$)

$\begin{cases}F_{ijl}(\tau)=&\langle\langle\xi_{i\uparrow}(\tau);c^{\dagger}_{j\downarrow}S_{l}^{+}\rangle\rangle\\\
G_{ijl}(\tau)=&\langle\langle\eta_{i\uparrow}(\tau);c^{\dagger}_{j\downarrow}S_{l}^{+}\rangle\rangle\end{cases}$
(105)

The next steps are the same as with the pair-pair term. The only difference
lies in the definition of $f^{n}$ in the resulting equations of motion. For
the spin-spin term it is defined as
$f^{n}_{kjl}=\langle\\{\psi^{n}_{k};c^{\dagger}_{j\downarrow}S_{l}^{+}\\}\rangle$.
Hence we arrive at the following equation

$\displaystyle\langle c_{i\uparrow}c^{\dagger}_{j\downarrow}S_{l}^{+}\rangle$
$\displaystyle=\frac{2}{2-n}\sum\limits_{k}(C^{11}_{ik}+C^{12}_{ik})f_{kjl}^{1}$
(106)
$\displaystyle+\frac{2}{n}\sum\limits_{k}(C^{12}_{ik}+C^{22}_{ik})f^{2}_{kjl}$

A computation of the $f^{n}$ leads to

$\displaystyle f^{1}_{kjl}$ $\displaystyle=\delta_{kj}\langle
S_{k}^{-}S_{l}^{+}\rangle-\delta_{kl}(C_{kj}^{11}+C_{kj}^{12})$ (107)
$\displaystyle f^{2}_{kjl}$ $\displaystyle=-\delta_{kj}\langle
S_{k}^{-}S_{l}^{+}\rangle-\delta_{kl}(C_{kj}^{12}+C_{kj}^{22})$

Therefore by replacing f by their explicit expression we obtain

$\displaystyle\langle c_{i\uparrow}c^{\dagger}_{j\downarrow}S_{l}^{+}\rangle$
$\displaystyle=\frac{2}{2-n}\sum\limits_{k}(C_{ik}^{11}+C_{ik}^{12})$ (108)
$\displaystyle[\delta_{kj}\langle
S_{k}^{-}S_{l}^{+}\rangle-\delta_{kl}(C_{kj}^{11}+C_{kj}^{12})]$
$\displaystyle+\frac{2}{n}\sum\limits_{k}(C_{ik}^{12}+C_{ik}^{22})$
$\displaystyle[-\delta_{kj}\langle
S_{k}^{-}S_{l}^{+}\rangle-\delta_{kl}(C_{kj}^{12}+C^{kj}_{22})]$

Which we can again simplify and set i=j to obtain

$\langle
S_{i}^{-}S_{l}^{+}\rangle=-\frac{\frac{2}{2-n}(C^{11}_{il}+C^{12}_{il})^{2}+\frac{2}{n}(C^{12}_{il}+C_{22}^{il})^{2}}{1-\frac{2}{2-n}(C^{11}_{0}+C^{12}_{0})+\frac{2}{n}(C^{12}_{0}+C^{22}_{0})}$
(109)

Which become with our notations

$\langle S_{i}^{-}S_{l}^{+}\rangle=\frac{\rho_{1}}{1-\phi}$ (110)

###  V.3 Charge-Charge term

As we did for the $\langle S_{i}^{-}S_{l}^{+}\rangle$ term we need to commute
the charge charge term so the first element can be decomposed using our
spinor. We then rewrite

$\displaystyle\langle
c^{\dagger}_{i\uparrow}c_{i\uparrow}n_{l\uparrow}\rangle=\frac{n}{2}-\langle
c_{i\uparrow}c^{\dagger}_{i\uparrow}n_{l\uparrow}\rangle$ (111)

We introduce the following Green’s functions

$\begin{cases}F_{ijl}(\tau)=&\langle\langle\xi_{i\uparrow}(\tau);c^{\dagger}_{j\uparrow}n_{l\uparrow}\rangle\rangle\\\
G_{ijl}(\tau)=&\langle\langle\eta_{i\uparrow}(\tau);c^{\dagger}_{j\uparrow}n_{l\uparrow}\rangle\rangle\end{cases}$
(112)

Once again the general form of the equation for $\langle
c_{i\sigma}c^{\dagger}_{i\sigma}n_{l\sigma}\rangle$ will be the same as for
the other 2 terms. However the definition of the involved $f^{n}_{kjl}$ will
not be the same. We compute
$f^{n}_{kjl}=\langle\\{\psi^{n}_{k};c^{\dagger}_{j\uparrow}n_{l\uparrow}\\}\rangle$

$\displaystyle f^{1}_{kjl}$ $\displaystyle=\delta_{kj}(\frac{n}{2}-\langle
n_{k\downarrow}n_{l\uparrow}\rangle)+\delta_{kl}(C^{11}_{kj}+C^{12}_{kj})$
(113) $\displaystyle f^{2}_{kjl}$ $\displaystyle=\delta_{kj}\langle
n_{k\downarrow}n_{l\uparrow}\rangle+\delta_{kl}(C^{12}_{kj}+C^{22}_{kj})$

We hence obtain

$\displaystyle\langle n_{i\uparrow}n_{l\uparrow}\rangle=$
$\displaystyle\frac{n}{2}-\rho_{1}-\phi\langle
n_{i\uparrow}n_{l\downarrow}\rangle+\frac{n}{2-n}(C^{11}_{0}+C^{12}_{0})$
(114)

We don’t know how to express $\langle n_{i\downarrow}n_{l\uparrow}\rangle$ as
a function of the correlations functions. So we need to redo this decoupling
on this term. This time we use

$\langle n_{i\downarrow}n_{l\uparrow}\rangle=\frac{n}{2}-\langle
c_{l\uparrow}c^{\dagger}_{l\uparrow}n_{i\downarrow}\rangle$ (115)

We therefore introduce

$\begin{cases}F_{ijl}(\tau)=&\langle\langle\xi_{i\uparrow}(\tau);c^{\dagger}_{j\uparrow}n_{l\downarrow}\rangle\rangle\\\
G_{ijl}(\tau)=&\langle\langle\eta_{i\uparrow}(\tau);c^{\dagger}_{j,\uparrow}n_{l\downarrow}\rangle\rangle\end{cases}$
(116)

The
$f^{n}_{kjl}=\langle\psi^{n}_{k},c^{\dagger}_{j\uparrow}n_{l\downarrow}\rangle$
are given by

$\displaystyle f^{1}_{kjl}$ $\displaystyle=\delta_{kj}(\frac{n}{2}-\langle
n_{i\downarrow}n_{l\downarrow}\rangle)$ (117) $\displaystyle f^{2}_{kjl}$
$\displaystyle=\delta_{kj}\langle n_{i\downarrow}n_{l\downarrow}\rangle$

Using the paramagnetic assumption we have $\langle
n_{i\uparrow}n_{l\uparrow}\rangle=\langle
n_{i\downarrow}n_{l\downarrow}\rangle$, leading to

$\displaystyle\langle n_{i\uparrow}n_{l\downarrow}\rangle=$
$\displaystyle\frac{n}{2}-\phi\langle
n_{i\uparrow}n_{l\uparrow}\rangle+\frac{n}{2-n}(C^{11}_{0}+C^{12}_{0})$ (118)

Figure 19: Fermi Surfaces and bands obtained by Pauli minimization in the
neighboorhood of half filling. Top row: 2% hole doped. Bottom row: 2% electron
doped.

If we inject this in the equation of $\langle
n_{i\uparrow}n_{l\uparrow}\rangle$, we obtain

$\displaystyle\langle n_{i\uparrow}n_{l\uparrow}\rangle$
$\displaystyle=-\rho_{1}+\phi^{2}\langle
n_{i\uparrow}n_{l\uparrow}\rangle+\frac{n}{2-n}(C^{11}_{0}+C^{12}_{0})(1-\phi)$
(119)

The last term can be simplified. An explicit computations of the $C_{0}$ leads
to $C^{11}_{0}+C_{0}^{12}=1-n+\langle n_{i\uparrow}n_{i\downarrow}\rangle$,
allowing us to show the last term is in fact just equal to $\frac{n^{2}}{4}$.
We therefore obtain

$\langle
n_{i\uparrow}n_{l\downarrow}\rangle=\frac{n^{2}}{4}-\frac{\rho_{1}}{1-\phi^{2}}$
(120)

Which is the self consistent equation we have. Combining the three terms,
since $p=\langle n_{i\uparrow}n_{l\uparrow}\rangle+\langle
S_{i}^{+}S_{l}^{-}\rangle-\langle\Delta_{i}\Delta_{l}^{*}\rangle$, we obtain
the following self consistent equation

$p=\frac{n^{2}}{4}-\frac{\rho_{1}}{1-\phi^{2}}-\frac{\rho_{1}}{1-\phi}-\frac{\rho_{3}}{1+\phi}$
(121)

## Appendix E: Bands at and around half-filling

Here we plot the bands for nearest-neighboors at half filling for the 3
minimization. We see the Fermi energy in between the two bands therefore we
have an insulator. Note we get the same result with next-nearest neighboors.
We plot the bands for Roth solution at 2% hole doping with its Fermi Surface:
we see the formation of a small hole pocket at $(\pi,\pi)$ when approaching
half filling. This can be interpreted as a magnetic instability [17], our
assumption of paramagnetism becomes invalid.

Figure 20: Band at half filling with t=1, U=8t. Note that all three solutions
(COM2 and Roth) are alike.

## Appendix F: Roth decoupling with superconductivity

In this appendix we show how Roth decoupling changes when we include
superconductivity with 4x4 basis. We will redo one of the terms of p as an
example, and perform the decoupling for the $\theta$ parameter. Let us start
by considering $\langle S_{i}^{-}S_{l}^{+}\rangle$ in p for example. The
superconducting basis is of size 4 and define by

$\psi=\begin{pmatrix}\xi_{i\uparrow}\\\ \eta_{i\uparrow}\\\
\xi^{\dagger}_{i\downarrow}\\\ \eta^{\dagger}_{i\downarrow}\end{pmatrix}$
(122)

To include the full basis we need to introduce 4 Green’s functions now

$\displaystyle A_{ijl}=$
$\displaystyle\langle\langle\xi_{i\uparrow};c_{j\downarrow}^{\dagger}S_{l}^{+}\rangle\rangle$
(123) $\displaystyle B_{ijl}=$
$\displaystyle\langle\langle\eta_{i\uparrow};c_{j\downarrow}^{\dagger}S_{l}^{+}\rangle\rangle$
$\displaystyle F_{ijl}=$
$\displaystyle\langle\langle\xi_{i\downarrow}^{\dagger};c_{j\downarrow}^{\dagger}S_{l}^{+}\rangle\rangle$
$\displaystyle G_{ijl}=$
$\displaystyle\langle\langle\eta_{i\downarrow}^{\dagger};c_{j\downarrow}^{\dagger}S_{l}^{+}\rangle\rangle$

Now the equations of motion has to be defined for the 4 Green’s functions.
Hence

$\displaystyle\partial_{\tau}\begin{pmatrix}A_{ijl}\\\ B_{ijl}\\\ F_{ijl}\\\
G_{ijl}\end{pmatrix}(\tau)$
$\displaystyle=\sum\limits_{k}E_{ik}\begin{pmatrix}A_{ijl}\\\ B_{ijl}\\\
F_{ijl}\\\
G_{ijl}\end{pmatrix}(\tau)+\delta(\tau)\begin{pmatrix}f^{1}_{ijl}\\\
f^{2}_{ijl}\\\ f^{3}_{ijl}\\\ f^{4}_{ijl}\end{pmatrix}$ (124)

With
$f^{n}_{ijl}=\langle\\{\psi^{n}_{i};c_{j\downarrow}^{\dagger}S_{l}^{+}\\}\rangle$.
As before we can use $(\omega-E_{ik})=S_{ik}(\omega)I^{-1}$ and integrate over
$\omega$ after we took the imaginary part to replace the composite Green’s
function matrix by a correlation function matrix. Thus

$\displaystyle-\langle
c_{j\downarrow}^{\dagger}c_{i\uparrow}S_{l}^{+}\rangle=$
$\displaystyle\frac{2}{2-n}\sum\limits_{k}[(C_{ik}^{11}+C_{ik}^{12})f^{1}_{kjl}+(C_{ik}^{12}+C_{ik}^{22})f^{2}_{kjl}$
(125)
$\displaystyle+(C^{13}_{ik}+C^{23}_{ik})f^{3}_{kjl}+(C^{14}_{ik}+C^{24}_{ik})f^{4}_{kjl}]$

In the last equation, the first line is the same as before, while the second
line are additionnal terms appearing with the superconducting basis. A bit of
algebra yields to

$\displaystyle f^{3}_{kjl}=$
$\displaystyle\delta_{kl}(C^{13}_{kj}+C^{14}_{kj})$ (126) $\displaystyle
f^{4}_{kjl}=$ $\displaystyle\delta_{kl}(C^{23}_{kj}+C^{24}_{kj})$

Finally the spin-spin term in p becomes

$\displaystyle\langle S_{i}^{-}S_{l}^{+}\rangle=$
$\displaystyle-\frac{\rho_{1}+\rho_{2}}{1-\phi}$ (127)

With
$\rho_{2}=\frac{2}{2-n}(C^{13}_{il}+C^{14}_{il})^{2}+\frac{2}{n}(C^{23}_{il}+C^{24}_{il})^{2}$.

Now we move on to the Roth decoupling for $\theta$. As Beenen and Edwards [10]
have mentioned, we have several ways of decoupling $\theta$ depending on
whether we consider $\langle c_{i\uparrow}c_{i\downarrow}n_{l\sigma}\rangle$
or for example $\langle
c^{\dagger}_{i\downarrow}c^{\dagger}_{i\uparrow}n_{l\sigma}\rangle$. Depending
on the decoupling scheme we will over estimate or underestimate the real value
of $\theta$ but the behaviour will remain globally the same [22]. Here we
consider a decoupling starting from $\langle\langle
c_{i\downarrow}^{\dagger},c_{i\uparrow}^{\dagger}n_{l\sigma}\rangle\rangle$.
We introduce the following Green’s functions

$\begin{pmatrix}A_{ijl}=&\langle\langle\xi_{i\uparrow};c_{j\uparrow}^{\dagger}n_{l\uparrow}\rangle\rangle\\\
B_{ijl}=&\langle\langle\eta_{i\uparrow};c_{j\uparrow}^{\dagger}n_{l\uparrow}\rangle\rangle\\\
F_{ijl}=&\langle\langle\xi_{i\downarrow}^{\dagger};c_{j\uparrow}^{\dagger}n_{l\uparrow}\rangle\rangle\\\
G_{ijl}=&\langle\langle\eta_{i\downarrow}^{\dagger};c_{j\uparrow}^{\dagger}n_{l\uparrow}\rangle\rangle\\\
\end{pmatrix}$ (128)

The decoupling is identical and by considering F+G we arrive to

$\displaystyle\langle
c^{\dagger}_{i\downarrow}c^{\dagger}_{j\uparrow}n_{l\sigma}\rangle$
$\displaystyle=\frac{2}{2-n}\sum\limits_{k}[(C_{ik}^{13}+C_{ik}^{14})f^{1}_{kjl}+(C_{ik}^{33}+C_{ik}^{44})f^{3}_{kjl}]$
(129)
$\displaystyle+\frac{2}{n}\sum\limits_{k}[(C^{23}_{ik}+C^{24}_{ik})f^{2}_{kjl}+(C^{34}_{ik}+C^{44}_{ik})f^{4}_{kjl}]$

With
$f^{n}_{kjl}=\langle\\{\psi^{n}_{k};c^{\dagger}_{j\uparrow}n_{l\sigma}\\}\rangle$.
Computing the f gives

$\displaystyle f^{1}_{kjl}=$ $\displaystyle\delta_{kj}(\frac{n}{2}-\langle
n_{k\uparrow}n_{l\uparrow}\rangle)+\delta_{kl}(C^{11}_{kj}+C^{12}_{kj})$ (130)
$\displaystyle f^{2}_{kjl}=$ $\displaystyle\delta_{kj}\langle
n_{k\uparrow}n_{l\uparrow}\rangle+\delta_{kl}(C^{12}_{kj}+C^{22}_{kj})$
$\displaystyle f^{3}_{kjl}=$ $\displaystyle-\delta_{kj}\frac{\theta_{kl}}{2}$
$\displaystyle f^{4}_{kjl}=$ $\displaystyle\delta_{kj}\frac{\theta_{kl}}{2}$

We therefore obtain

$\displaystyle\frac{\theta}{2}=$ $\displaystyle\frac{\langle
n_{i\uparrow}n_{l\uparrow}\rangle[\frac{2}{n}(C^{23}_{0}+C^{24}_{0})-\frac{2}{2-n}(C^{13}_{0}+C^{14}_{0})]+\zeta}{1+\phi}$
(131) $\displaystyle+\frac{n}{2}\frac{C^{13}_{0}+C^{14}_{0}}{1+\phi}$

With
$\zeta=\frac{2}{2-n}(C^{11}_{il}+C^{12}_{il})(C^{13}_{il}+C^{14}_{il})+\frac{2}{n}(C^{12}_{il}+C^{22}_{il})(C^{23}_{il}+C^{24}_{il})$.
Finally by noting that $C^{23}_{0}+C^{24}_{0}=0$ and
$C^{13}_{0}+C^{14}_{0}=0$, we obtain the equation we used

$\frac{\theta}{2}=\frac{\zeta}{1+\phi}$ (132)

## Appendix G: Effects of particle-hole transformation

In this appendix we give some details on how we derived the particle-hole
relations in Eq. (54). The relation for the chemical potential is obtained by
using the fact that the Hubbard hamiltonian stays invariant under this
transformation.

$\displaystyle H=$
$\displaystyle\sum\limits_{ij\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum\limits
n_{i\uparrow}n_{i\downarrow}+\mu\sum\limits_{i\sigma}n_{i\sigma}$ (133)
$\displaystyle\rightarrow\sum\limits_{ij\sigma}t_{ij}(-1)^{i+j}c_{i\sigma}c^{\dagger}_{j\sigma}+U\sum\limits_{i}(-1)^{4i}c_{i\uparrow}c^{\dagger}_{i\uparrow}c_{i\downarrow}c^{\dagger}_{i\downarrow}$
$\displaystyle-\mu\sum\limits_{i\sigma}(-1)^{2i}c_{i\sigma}c^{\dagger}_{i\sigma}$
$\displaystyle=-\sum\limits_{ij\sigma}t_{ij}c_{i\sigma}c^{\dagger}_{j\sigma}+U\sum\limits_{i}\left(c_{i\downarrow}c^{\dagger}_{i\downarrow}-n_{i\uparrow}c_{i\downarrow}c^{\dagger}_{i\downarrow}\right)+\mu\sum\limits_{i\sigma}n_{i\sigma}$
$\displaystyle=\sum\limits_{ij\sigma}t_{ij}c^{\dagger}_{j\sigma}c_{i\sigma}+U\sum\limits_{i}\left(1-n_{i\downarrow}-n_{i\uparrow}+n_{i\uparrow}n_{i\downarrow}\right)+\mu\sum\limits_{i\sigma}n_{i\sigma}$
$\displaystyle=\sum\limits_{ij\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum\limits_{i}n_{i\uparrow}n_{i\downarrow}+(\mu-U)\sum\limits_{i\sigma}n_{i\sigma}+cste$

Thus, to keep the Hamiltonian invariant and therefore have under the particle-
hole transformation $H\rightarrow H$, we need to impose

$\mu(2-n)=-(\mu(n)-U)$ (134)

Which is the relation we get in the main text. For e and p we work directly
with their definitions

$\displaystyle e=$
$\displaystyle\langle\xi_{i\sigma}\xi_{j\sigma}^{\dagger}\rangle-\langle\eta_{i\sigma}\eta^{\dagger}_{j\sigma}\rangle$
(135) $\displaystyle p=$ $\displaystyle\langle
n_{i\sigma}n_{j\sigma}\rangle+\langle
S_{i}^{+}S_{j}^{-}\rangle-\langle\Delta_{i}\Delta_{j}^{*}\rangle$

Under the particle-hole transformation we have

$\xi_{i\sigma}\rightarrow(-1)^{i}\eta_{i\sigma}^{\dagger}\ \ \ \ \ \ \ \ \ \ \
\ \eta_{i\sigma}\rightarrow(-1)^{i}\xi_{i\sigma}^{\dagger}$ (136)

Hence

$\displaystyle e(2-n)$
$\displaystyle\rightarrow(-1)^{i+j}(\langle\eta_{i\sigma}^{\dagger}\eta_{j\sigma}\rangle-\langle\xi_{i\sigma}^{\dagger}\xi_{j\sigma}\rangle)$
(137)
$\displaystyle=\langle\eta_{j\sigma}\eta^{\dagger}_{i\sigma}\rangle-\langle\xi_{j\sigma}\xi^{\dagger}_{i\sigma}\rangle$
$\displaystyle=-(\langle\xi_{i\sigma}\xi^{\dagger}_{j\sigma}\rangle-\langle\eta_{i\sigma}\eta^{\dagger}_{j\sigma}\rangle)$
$\displaystyle=-e(n)$

We didn’t kept the terms with $\delta_{ij}$ because e and p always appear with
a $t_{ij}$ prefactor and $t_{ij}=0$ if $i=j$. We used the fact that i and j
are always nearest-neighbours to get $(-1)^{i+j}=-1$. For p we have

$\displaystyle p$ $\displaystyle\rightarrow(-1)^{2i+2j}\langle
c_{i\sigma}c^{\dagger}_{i\sigma}c_{j\sigma}c^{\dagger}_{j\sigma}+c_{i\uparrow}c^{\dagger}_{i\downarrow}c_{j\downarrow}c^{\dagger}_{j\uparrow}$
(138)
$\displaystyle-c^{\dagger}_{i\uparrow}c^{\dagger}_{i\downarrow}c_{j\downarrow}c_{j\uparrow}\rangle$
$\displaystyle=\langle
c_{j\sigma}c^{\dagger}_{j\sigma}-n_{i\sigma}c_{j\sigma}c^{\dagger}_{j\sigma}+S_{i}^{-}S_{j}^{+}-\Delta_{j}^{*}\Delta_{i}\rangle$
$\displaystyle=\langle
1-n_{j\sigma}-n_{i\sigma}+n_{i\sigma}n_{j\sigma}+S_{j}^{+}S_{i}^{-}-\Delta_{i}\Delta_{j}^{*}\rangle$
$\displaystyle=1-\frac{n}{2}-\frac{n}{2}+p(n)$ $\displaystyle=(1-n)+p(n)$

Which is the other relation we had.