# Perturbative description of high-$T_{c}$ superconductivity in the Hubbard
model via unitary transformation and classical spins

Kirill Alpin  Max Planck Institute for Solid State Research, Heisenbergstrasse
1, D-70569 Stuttgart, Germany

(September 1, 2024)

###### Abstract

A unitary transformation is applied to the Hubbard model, which maps the
Hubbard interaction to a single particle term. The resulting Hamiltonian
consists of unconstrained fermions, which is then mapped to a Hamiltonian of
spinless fermions coupled to pseudospins. The fermions are integrated out
using second order perturbation theory in $1/U$, resulting in an effective
spin Hamiltonian. An order parameter is identified, stabilizing d-wave
superconductivity. The groundstate energy of classical spin configurations is
minimized at a finite value of this order parameter after a critical chemical
potential, resulting in d-wave superconductivity at non-zero doping. This
suggests, that the onset of high-$T_{c}$ superconductivity is governed by the
groundstate of a classical spin system.

##  I Introduction

The Hubbard model in the strong coupling regime is believed to describe the
emergence of high-$T_{c}$ superconductivity in cuprates [1]. The pairing in
these compounds is further believed to be d-wave [2]. There have been numerous
previous theoretical studies of the Hubbard model [3, 4, 5]. This paper shows
a perturbative approach to observing this superconductivity in the strongly
repulsive Hubbard model at finite doping. The first step, a unitary
transformation, is similar to [6], where auxiliary spins are coupled to the
original Hamiltonian via a unitary Schrieffer-Wolff transformation. Unlike
this approach, here the unitary transformation is applied exactly. This
results in a Hamiltonian with unconstrained fermions, unlike the t-J model
[7]. The interaction terms now scale with $t,\mu$ instead of $U$, so
perturbative approaches in $1/U$ are viable. This is also not given in the t-J
model, where interactions are of the same order as the free part of the
Hamiltonian, due to the correlated hopping terms. Therefore, in the t-J model,
one must resort to uncontrolled slave-boson approaches [7]. The Hamiltonian is
then mapped to a system of fermions coupled to pseudospins. In this form,
unlike [6], the Hilbert space is not artificially expanded. It is solved using
variational second order perturbation theory and spins are replaced with
classical spins. Beyond [8], this perturbative expansion can be done at finite
doping. An order parameter is identified, which stabilizes d-wave
superconductivity. The groundstate energy as a function of this order
parameter shows a second order phase transition at a critical chemical
potential, resulting in the development of d-wave superconductivity. The flow
of the approach described in this paper is shown schematically in Fig. 1,
which we also apply to the Hubbard model with $t^{\prime}=-0.3t$, showing
characteristics of cuprate superconductivity.

Figure 1: Flow diagram of the $1/U$ perturbative approach described in this
paper.

Notable similar approaches are ones based on strong coupling expansions [9,
10] around the atomic limit of the Hubbard model. Methods like these, in which
class for example Dynamical Mean Field Theory (DMFT) [11, 12, 13], Self-Energy
Functional Theory [14] and Cluster Perturbation Theory [14, 15] belong to,
capture the local correlations of a given system, while it is a challenge to
incorporate non-local self-energies [16, 17]. The unitary transformation
defined in this manuscript is applied exactly, so no correlations are lost
there. The perturbative expansion following can in principle capture any non-
locality, since no diagrams are discarded and this series can be carried out
to higher order than the one in this paper, as it has been done in studies
using DiagMC [18]. Although there, the expansion is carried out in the weak
coupling regime, while here we offer to do the same at strong coupling (using
for example ordinary diagrammatics in the case of the Hamiltonian with
unconstrained fermions). Other weak coupling expansions include [19, 20].

##  II Local unitary transformation

The square lattice repulsive Hubbard model with $U>0$ and $t<0$ is given by

$\displaystyle H$ $\displaystyle=\sum_{{\bf
r},{\bm{\delta}},\sigma}t(c^{\dagger}_{{\bf r},\sigma}c_{{\bf
r}+{\bm{\delta}},\sigma}+c^{\dagger}_{{\bf r}+{\bm{\delta}},\sigma}c_{{\bf
r},\sigma})$ $\displaystyle+\sum_{{\bf r}}U\left(n_{{{\bf
r}},\uparrow}-\frac{1}{2}\right)\left(n_{{{\bf
r}},\downarrow}-\frac{1}{2}\right)+\sum_{{\bf r},\sigma}\mu n_{{{\bf
r}},\sigma}.$ (1)

In this convention half filling is at $\mu=0$ and the system is hole doped
with $\mu>0$. ${\bm{\delta}}$ runs over the two vectors $(1,0)$, $(0,1)$. We
apply a unitary transformation

$\displaystyle\mathcal{U}=\prod_{{\bf r}}U_{\bf r}$ (2)

to the Hamiltonian (1), which decomposes into a series of local unitaries
$U_{\bf r}$. We consider the following unitary

$\displaystyle U_{\bf r}=1-n_{{\bf r},2}-n_{{\bf r},2}\bar{c}_{{\bf
r},1}\tilde{c}_{{\bf r},3}$ (3)

where we relabeled $c_{{\bf r},1}\equiv c_{{\bf r},\uparrow}$, $c_{{\bf
r},2}\equiv c_{{\bf r},\downarrow}$. We have introduced a third fermion
species $c_{{\bf r},3}$ at every lattice site, which has no contribution in
the original Hamiltonian (1). It does not contribute to any observables. This
auxiliary fermion is necessary to construct a parity conserving unitary
operation. Parity conservation is needed to form a local $U_{{\bf r}}$ (a non-
parity conserving version of this transformation, dropping the third fermion
species, can still be used in 1D).

We used the following definitions in Eq. (3)

$\displaystyle\bar{c}_{{\bf r},n}$ $\displaystyle\equiv c^{\dagger}_{{\bf
r},n}-c_{{\bf r},n},$ (4) $\displaystyle\tilde{c}_{{\bf r},n}$
$\displaystyle\equiv c^{\dagger}_{{\bf r},n}+c_{{\bf r},n}.$ (5)

The unitary transformation (3) performs the following mapping from the
original basis $\ket{\text{fermions},n_{3}}$ to the $\ket{n_{1},n_{2},n_{3}}$
basis

$\displaystyle\ket{0,0}$
$\displaystyle\shortrightarrow\ket{0,0,0},\quad\ket{\uparrow,0}\shortrightarrow\ket{1,0,0},$
(6) $\displaystyle\ket{\downarrow,0}$
$\displaystyle\shortrightarrow\ket{1,1,1},\quad\ket{\updownarrows,0}\shortrightarrow\ket{0,1,1},$
(7) $\displaystyle\ket{0,1}$
$\displaystyle\shortrightarrow\ket{0,0,1},\quad\ket{\uparrow,1}\shortrightarrow\ket{1,0,1},$
(8) $\displaystyle\ket{\downarrow,1}$
$\displaystyle\shortrightarrow\ket{1,1,0},\quad\ket{\updownarrows,1}\shortrightarrow\ket{0,1,0}.$
(9)

Notice, that all, at large $U$ and half filling, prohibited states, so double
and $0$ occupancy states, are mapped to a vacuum state of fermion species 1,
while all low energy states are mapped to ones, where fermion 1 is fully
occupied. This allows us, to effectively integrate out double and $0$
occupancies by perturbation theory. With

$\displaystyle\bar{c}^{2}_{{\bf r},n}=-1,\quad\tilde{c}^{2}_{{\bf r},n}=1$
(10)

we can show that $U^{\dagger}_{\bf r}=U_{\bf r}$ and that $U^{\dagger}_{\bf
r}U_{\bf r}=1=U^{2}_{\bf r}$. So $\mathcal{U}$ is unitary
$\mathcal{U}^{\dagger}\mathcal{U}=1$. Because $U_{\bf r}$ is parity
conserving, it commutes with itself at any other site (and on the same site)
$[U_{\bf r},U_{\bf r^{\prime}}]=0$.

Applying $U_{\bf r}$ to the Hubbard interaction maps it to a single particle
term

$\displaystyle U_{\bf
r}\left(n_{\uparrow}-\frac{1}{2}\right)\left(n_{\downarrow}-\frac{1}{2}\right)U_{\bf
r}=\frac{1}{4}-\frac{1}{2}n_{{\bf r},1}$ (11)

where we used $U^{\dagger}_{\bf r}=U_{\bf r}$. The chemical potential term is
mapped to

$\displaystyle U_{\bf r}\sum_{\sigma}n_{{\bf r},\sigma}U_{\bf r}=n_{{\bf
r},1}+2n_{{\bf r},2}(1-n_{{\bf r},1}).$ (12)

Finally, the hopping term is mapped to

$\displaystyle U_{{\bf r}+{\bm{\delta}}}U_{\bf
r}\sum_{\sigma}(c^{\dagger}_{{\bf r},\sigma}c_{{\bf
r}+{\bm{\delta}},\sigma}+c^{\dagger}_{{\bf r}+{\bm{\delta}},\sigma}c_{{\bf
r},\sigma})U_{\bf r}U_{{\bf r}+{\bm{\delta}}}$ (13) $\displaystyle=$
$\displaystyle K_{{\bf r},{\bm{\delta}}}c^{\dagger}_{{\bf
r},1}c^{\dagger}_{{\bf r}+{\bm{\delta}},1}+(1+A_{{\bf
r},{\bm{\delta}}})c^{\dagger}_{{\bf r},1}c_{{\bf
r}+{\bm{\delta}},1}+\text{h.c.}$ (14)

with

$\displaystyle K_{{\bf r},{\bm{\delta}}}=$ $\displaystyle-n_{{\bf
r},2}+n_{{\bf r}+{\bm{\delta}},2}$ $\displaystyle-(c^{\dagger}_{{\bf
r},2}c_{{\bf r}+{\bm{\delta}},2}+c^{\dagger}_{{\bf r}+{\bm{\delta}},2}c_{{\bf
r},2})\tilde{c}_{{\bf r},3}\tilde{c}_{{\bf r}+{\bm{\delta}},3},$ (15)
$\displaystyle A_{{\bf r},{\bm{\delta}}}=$ $\displaystyle-n_{{\bf
r},2}-n_{{\bf r}+{\bm{\delta}},2}$ $\displaystyle+(c^{\dagger}_{{\bf
r},2}c_{{\bf r}+{\bm{\delta}},2}+c^{\dagger}_{{\bf r}+{\bm{\delta}},2}c_{{\bf
r},2})\tilde{c}_{{\bf r},3}\tilde{c}_{{\bf r}+{\bm{\delta}},3}.$ (16)

Let’s apply $\mathcal{U}$ to the whole Hubbard model (1)

$\displaystyle\mathcal{U}^{\dagger}H\mathcal{U}$ $\displaystyle=\prod_{{\bf
r}^{\prime}}U_{{\bf r}^{\prime}}H\prod_{{\bf r}}U_{\bf r}$ (17)
$\displaystyle=t\sum_{{\bf r},{\bm{\delta}},\sigma}U_{{\bf
r}+{\bm{\delta}}}U_{\bf r}(c^{\dagger}_{{\bf r},\sigma}c_{{\bf
r}+{\bm{\delta}},\sigma}+c^{\dagger}_{{\bf r}+{\bm{\delta}},\sigma}c_{{\bf
r},\sigma})U_{\bf r}U_{{\bf r}+{\bm{\delta}}}$ $\displaystyle+U\sum_{{\bf
r}}U_{\bf r}\left(n_{{{\bf r}},\uparrow}-\frac{1}{2}\right)\left(n_{{{\bf
r}},\downarrow}-\frac{1}{2}\right)U_{\bf r}$ $\displaystyle+\mu\sum_{{\bf
r},\sigma}U_{\bf r}n_{{{\bf r}},\sigma}U_{\bf r}.$ (18)

The order of the application of $U_{{\bf r}}$ in the two products above is not
important, since $U_{{\bf r}}$ commute. We further used that for any bosonic
or fermionic operators $X_{\bf a},Y_{{\bf b}}$ we have

$\displaystyle\prod_{{\bf r}^{\prime}}U_{{\bf r}^{\prime}}X_{\bf a}Y_{{\bf
b}}\prod_{{\bf r}}U_{{\bf r}}$ (19) $\displaystyle=U_{{\bf b}}U_{\bf a}X_{\bf
a}Y_{{\bf b}}U_{\bf a}U_{{\bf b}}$ (20)

because $[X_{\bf a}Y_{{\bf b}},U_{\bf r}]=0$ as long as ${\bf a}\neq{\bf r}$
and ${\bf b}\neq{\bf r}$.

So the transformed Hamiltonian becomes

$\displaystyle H^{\prime}$ $\displaystyle=\mathcal{U}^{\dagger}H\mathcal{U}$
$\displaystyle=t\sum_{{\bf r},{\bm{\delta}}}K_{{\bf
r},{\bm{\delta}}}c^{\dagger}_{{\bf r},1}c^{\dagger}_{{\bf
r}+{\bm{\delta}},1}+(1+A_{{\bf r},{\bm{\delta}}})c^{\dagger}_{{\bf
r},1}c_{{\bf r}+{\bm{\delta}},1}+\text{h.c.}$ $\displaystyle+\sum_{{\bf
r}}-\frac{U}{2}n_{{\bf r},1}+\mu(n_{{\bf r},1}+2n_{{\bf r},2}(1-n_{{\bf
r},1})).$ (21)

Since the Hubbard interaction is mapped to just a chemical potential of the
fermion species 1, the interactions in $K_{{\bf r},{\bm{\delta}}}$ and
$A_{{\bf r},{\bm{\delta}}}$ now scale with $t$ instead of $U$. Also notice,
that, unlike the t-J model, Hamiltonian (21) describes unconstrained fermions.
Therefore, in this basis, the large interaction limit $U\gg t$ can be studied
perturbatively at $\mu\neq 0$. Further, the mapping we have done is exact and
works at any $U,t,\mu$.

Note that $H^{\prime}$ commutes with any bond operator $i\bar{c}_{{\bf
r},3}\bar{c}_{{\bf r}^{\prime},3}$, so one could block diagonalize it in the
eigenbasis of all $i\bar{c}_{{\bf r},3}\bar{c}_{{\bf r}^{\prime},3}$.

Ultimately, we would like to compute observables in this Hamiltonian. To do
so, we also need to transform those too. The transformed magnetization is

$\displaystyle U_{\bf r}(n_{{\bf r},\uparrow}-n_{{\bf r},\downarrow})U_{\bf
r}=n_{{\bf r},1}(1-2n_{{\bf r},2}).$ (22)

The transformed particle density is

$\displaystyle U_{\bf r}(n_{{\bf r},\uparrow}+n_{{\bf r},\downarrow})U_{\bf
r}=2n_{{\bf r},2}(1-n_{{\bf r},1})+n_{{\bf r},1}.$ (23)

A two site singlet pair operator transformed into the new basis is

$\displaystyle\Delta=U_{{\bf r}+{\bm{\delta}}}U_{\bf
r}\sum_{\sigma}(-1)^{\sigma}(c_{{\bf r}+{\bm{\delta}},\bar{\sigma}}c_{{\bf
r},\sigma})U_{\bf r}U_{{\bf r}+{\bm{\delta}}}=$ (24) $\displaystyle c_{{{\bf
r}},1}c_{{{\bf r}+{\bm{\delta}}},1}(-c_{{{\bf r}},2}\tilde{c}_{{{\bf
r}},3}+c_{{{\bf r}+{\bm{\delta}}},2}\tilde{c}_{{{\bf r}+{\bm{\delta}}},3}$
$\displaystyle-n_{{{\bf r}},2}c_{{{\bf r}+{\bm{\delta}}},2}\tilde{c}_{{{\bf
r}+{\bm{\delta}}},3}+n_{{{\bf r}+{\bm{\delta}}},2}c_{{{\bf
r}},2}\tilde{c}_{{{\bf r}},3})+$ $\displaystyle c^{\dagger}_{{{\bf
r}+{\bm{\delta}}},1}c_{{{\bf r}},1}(c_{{{\bf r}},2}\tilde{c}_{{{\bf
r}},3}+n_{{{\bf r}},2}c_{{{\bf r}+{\bm{\delta}}},2}\tilde{c}_{{{\bf
r}+{\bm{\delta}}},3}-n_{{{\bf r}+{\bm{\delta}}},2}c_{{{\bf
r}},2}\tilde{c}_{{{\bf r}},3})+$ $\displaystyle c^{\dagger}_{{{\bf
r}},1}c_{{{\bf r}+{\bm{\delta}}},1}(c_{{{\bf
r}+{\bm{\delta}}},2}\tilde{c}_{{{\bf r}+{\bm{\delta}}},3}-n_{{{\bf
r}},2}c_{{{\bf r}+{\bm{\delta}}},2}\tilde{c}_{{{\bf
r}+{\bm{\delta}}},3}+n_{{{\bf r}+{\bm{\delta}}},2}c_{{{\bf
r}},2}\tilde{c}_{{{\bf r}},3})+$ $\displaystyle c^{\dagger}_{{{\bf
r}+{\bm{\delta}}},1}c^{\dagger}_{{{\bf r}},1}(n_{{{\bf r}},2}c_{{{\bf
r}+{\bm{\delta}}},2}\tilde{c}_{{{\bf r}+{\bm{\delta}}},3}-n_{{{\bf
r}+{\bm{\delta}}},2}c_{{{\bf r}},2}\tilde{c}_{{{\bf r}},3})$ (25)

with $\bar{\sigma}=-\sigma$, $\uparrow=0$ and $\downarrow=1$.

##  III Perturbation theory

In the groundstate and at large $U$ and $\mu=0$, the fermion species 1 is, due
to the $-\frac{U}{2}$ chemical potential, fully occupied. In this limit, we
can therefore perform perturbation theory around $\ket{\Psi}=\prod_{{\bf
r}}c^{\dagger}_{{{\bf r}},1}\ket{0}$.

The first order is zero. Second order is

$\displaystyle H_{2}=-\frac{1}{U}\bra{\Psi}H^{\prime 2}\ket{\Psi}.$ (26)

The factor $-\frac{1}{U}$ comes from the energy needed to excite two $c_{1}$
fermions. $H_{2}$ becomes

$\displaystyle H_{2}=$ $\displaystyle-\frac{2t^{2}}{U}\sum_{{\bf
r},{\bm{\delta}}}n_{{{\bf r}},2}+n_{{{\bf r}+{\bm{\delta}}},2}-2n_{{{\bf
r}},2}n_{{{\bf r}+{\bm{\delta}}},2}$ $\displaystyle-(c^{\dagger}_{{{\bf
r}},2}c_{{{\bf r}+{\bm{\delta}}},2}-c^{\dagger}_{{{\bf
r}+{\bm{\delta}}},2}c_{{{\bf r}},2})\tilde{c}_{{{\bf r}},3}\tilde{c}_{{{\bf
r}+{\bm{\delta}}},3}.$ (27)

Consider the following composite particles

$\displaystyle b^{\dagger}_{{\bf r}}=\tilde{c}_{{{\bf r}},3}c^{\dagger}_{{{\bf
r}},2}\quad b_{{{\bf r}}}=c_{{{\bf r}},2}\tilde{c}_{{{\bf r}},3}.$ (28)

We find

$\displaystyle[b_{{\bf r}},b_{{\bf r}^{\prime}}]$ $\displaystyle=0$ (29)
$\displaystyle[b^{\dagger}_{{\bf r}},b^{\dagger}_{{\bf r}^{\prime}}]$
$\displaystyle=0$ (30) $\displaystyle[b_{{\bf r}},b^{\dagger}_{{\bf
r}^{\prime}}]$ $\displaystyle=(1-2b^{\dagger}_{{\bf r}}b_{{{\bf
r}}})\delta_{{{\bf r}},{\bf r}^{\prime}}.$ (31)

These are the commutation relations of hardcore bosons. Eq. (28) is equivalent
to the drone-fermion representation of spins [21]. $n_{{{\bf r}},2}$ can be
written in $b$ like so

$\displaystyle n_{{{\bf r}},2}=$ $\displaystyle c^{\dagger}_{{{\bf
r}},2}c_{{{\bf r}},2}=c^{\dagger}_{{{\bf r}},2}c_{{{\bf
r}},2}\tilde{c}^{2}_{{{\bf r}},3}$ (32) $\displaystyle=$
$\displaystyle\tilde{c}_{{{\bf r}},3}c^{\dagger}_{{{\bf r}},2}c_{{{\bf
r}},2}\tilde{c}_{{{\bf r}},3}=b^{\dagger}_{{\bf r}}b_{{\bf r}}=n_{{{\bf
r}},b}$ (33)

where we used $\tilde{c}^{2}_{{{\bf r}},3}=1$. $H_{2}$ in hardcore bosons is

$\displaystyle H_{2}=$ $\displaystyle-\frac{2t^{2}}{U}\sum_{{\bf
r},{\bm{\delta}}}n_{{{\bf r}},b}+n_{{{\bf r}+{\bm{\delta}}},b}-2n_{{{\bf
r}},b}n_{{{\bf r}+{\bm{\delta}}},b}$ $\displaystyle-b^{\dagger}_{{{\bf
r}},2}b_{{{\bf r}+{\bm{\delta}}},2}-b^{\dagger}_{{{\bf
r}+{\bm{\delta}}},2}b_{{{\bf r}},2}.$ (34)

Further, we can write hardcore bosons as spins

$\displaystyle S^{+}_{{\bf r}}=$ $\displaystyle b_{{\bf r}}\quad S^{-}_{{\bf
r}}=b^{\dagger}_{{\bf r}}$ (35) $\displaystyle S^{z}_{{\bf r}}=$
$\displaystyle\frac{1}{2}-b^{\dagger}_{{\bf r}}b_{{\bf r}}$ (36)
$\displaystyle S^{x}_{{\bf r}}=$ $\displaystyle\frac{1}{2}(b_{{\bf
r}}+b^{\dagger}_{{\bf r}})$ (37) $\displaystyle S^{y}_{{\bf r}}=$
$\displaystyle\frac{1}{2i}(b_{{\bf r}}-b^{\dagger}_{{\bf r}})$ (38)

This reproduces the well known $U\gg 1$ limit of the Hubbard model at half
filling, that is the Heisenberg model

$\displaystyle H_{2}=$ $\displaystyle\frac{4t^{2}}{U}\sum_{{\bf
r},{\bm{\delta}}}{\bf S}_{{\bf r}}\cdot{\bf S}_{{\bf r}+{\bm{\delta}}}.$ (39)

Notice that the reverse transformation of a general spin Hamiltonian, that is
first to hardcore bosons and then to the composite fermions containing
$\tilde{c}_{{{\bf r}},3}$ can always be done and results in a Hamiltonian with
unconstrained fermions [21].

##  IV Variational perturbation theory at finite doping

In this section, we perform perturbation theory at finite $\mu$ with the goal
of calculating the expectation value of $\Delta$. In terms of hardcore bosons,
the $\Delta$ order parameter in Eq. (25) transforms to

$\displaystyle\Delta=c_{{{\bf r}}}c_{{{\bf r}+{\bm{\delta}}}}(-b_{{\bf
r}}+b_{{\bf r}+{\bm{\delta}}}$ $\displaystyle-n_{{{\bf r}},b}b_{{\bf
r}+{\bm{\delta}}}+n_{{{\bf r}+{\bm{\delta}}},b}b_{{\bf r}})+$ $\displaystyle
c^{\dagger}_{{{\bf r}+{\bm{\delta}}}}c_{{{\bf r}}}(b_{{\bf r}}+n_{{{\bf
r}},b}b_{{\bf r}+{\bm{\delta}}}-n_{{{\bf r}+{\bm{\delta}}},b}b_{{\bf r}})+$
$\displaystyle c^{\dagger}_{{{\bf r}}}c_{{{\bf r}+{\bm{\delta}}}}(b_{{\bf
r}+{\bm{\delta}}}-n_{{{\bf r}},b}b_{{\bf r}+{\bm{\delta}}}+n_{{{\bf
r}+{\bm{\delta}}},b}b_{{\bf r}})+$ $\displaystyle c^{\dagger}_{{{\bf
r}+{\bm{\delta}}}}c^{\dagger}_{{{\bf r}}}(n_{{{\bf r}},b}b_{{\bf
r}+{\bm{\delta}}}-n_{{{\bf r}+{\bm{\delta}}},b}b_{{\bf r}}).$ (40)

Hardcore bosons can be expressed as spin operators, which in turn reduces the
singlet order parameter $\Delta$ to a spin wave order parameter. Note, that an
uneven number of hardcore boson operators are present in $\Delta$. When the
expectation value of $\Delta$ is computed, even at first order in
perturbation, an uneven number of hardcore boson operators remain. Translated
into the spin picture, this means that to have a finite expectation value of
$\Delta$, that is superconductivity in any way, pseudospin U(1) symmetry must
be broken. This U(1) symmetry breaking can also be seen when $\Delta$ is
multiplied with a global U(1) gauge of the original fermionic operators
$c_{{{\bf r}},\sigma}\rightarrow\exp(i\theta/2)c_{{{\bf r}},\sigma}$. The only
place where this phase can consistently be absorbed to are the $b_{{\bf r}}$
operators. A phase in $\exp(i\theta)b_{{\bf r}}$ translates to rotations
around the z-axis of the pseudospins $\cos(\theta)S_{x}+\sin(\theta)S_{y}$,
see Eq. (35). So charge $U(1)$ symmetry is transformed using the local unitary
transformation (3) to a pseudospin $U(1)$ symmetry, which we expect to be
broken in a superconductor. We can transform the Hamiltonian (21) into one
containing hardcore bosons

$\displaystyle H^{\prime}$ $\displaystyle=t\sum_{{\bf r},{\bm{\delta}}}K_{{\bf
r},{\bm{\delta}}}c^{\dagger}_{{\bf r}}c^{\dagger}_{{\bf
r}+{\bm{\delta}}}+(1+A_{{\bf r},{\bm{\delta}}})c^{\dagger}_{{\bf r}}c_{{\bf
r}+{\bm{\delta}}}+\text{h.c.}$ $\displaystyle+\sum_{{\bf
r}}-\frac{U}{2}n_{{\bf r}}+\mu(n_{{\bf r}}+2n_{{\bf r},b}(1-n_{{\bf r}})).$
(41)

We dropped the 1 index for the $c_{1}$ fermions. The $K$ and $A$ terms
expressed in hardcore bosons are

$\displaystyle K_{{{\bf r}},{\bm{\delta}}}$ $\displaystyle=-n_{{{\bf
r}},b}+n_{{{\bf r}+{\bm{\delta}}},b}-b^{\dagger}_{{\bf r}}b_{{\bf
r}+{\bm{\delta}}}+b^{\dagger}_{{\bf r}+{\bm{\delta}}}b_{{\bf r}}$ (42)
$\displaystyle A_{{{\bf r}},{\bm{\delta}}}$ $\displaystyle=-n_{{{\bf
r}},b}-n_{{{\bf r}+{\bm{\delta}}},b}+b^{\dagger}_{{\bf r}}b_{{\bf
r}+{\bm{\delta}}}-b^{\dagger}_{{\bf r}+{\bm{\delta}}}b_{{\bf r}}$ (43)

Fourier transforming everything, we end up with the following interaction
vertices

$\displaystyle V_{K}$ $\displaystyle=\sum_{{\bm{\delta}}}tK_{{\bf
r},{\bm{\delta}}}c^{\dagger}_{{\bf r}}c^{\dagger}_{{\bf
r}+{\bm{\delta}}}+\text{h.c.}$
$\displaystyle=\sum_{{\bm{\delta}}}tb^{\dagger}_{{{\bf k}}_{1}}b_{{{\bf
k}}_{2}}{\bf c}^{\dagger}_{{{\bf k}}_{3}}\cdot(-k({\bf 0},{\bf
0},{\bm{\delta}})+k({\bm{\delta}},{\bm{\delta}},{\bm{\delta}})$
$\displaystyle-k({\bf 0},{\bm{\delta}},{\bm{\delta}})+k({\bm{\delta}},{\bf
0},{\bm{\delta}}))\cdot{\bf c}_{{{\bf k}}_{4}}$ (44) $\displaystyle V_{A+\mu}$
$\displaystyle=-2\mu n_{{{\bf r}},b}n_{{{\bf
r}}}+\sum_{{\bm{\delta}}}t(A_{{\bf r},{\bm{\delta}}}c^{\dagger}_{{\bf
r}}c_{{\bf r}+{\bm{\delta}}}+A^{\dagger}_{{\bf
r},{\bm{\delta}}}c^{\dagger}_{{\bf r}+{\bm{\delta}}}c_{{\bf r}})$
$\displaystyle=b^{\dagger}_{{{\bf k}}_{1}}b_{{{\bf k}}_{2}}{\bf
c}^{\dagger}_{{{\bf k}}_{3}}\cdot(-2\mu a({\bf 0},{\bf 0},{\bf 0})$
$\displaystyle+t\sum_{{\bm{\delta}}}a({\bf 0},{\bf
0},{\bm{\delta}})-a({\bm{\delta}},{\bm{\delta}},{\bm{\delta}})$
$\displaystyle+a({\bf 0},{\bm{\delta}},{\bm{\delta}})-a({\bm{\delta}},{\bf
0},{\bm{\delta}}))\cdot{\bf c}_{{{\bf k}}_{4}}$ (45)

with

$\displaystyle k({\bm{\delta}}_{1},{\bm{\delta}}_{2},{\bm{\delta}})$
$\displaystyle=\begin{pmatrix}0&e^{i{{\bf k}}_{4}{\bm{\delta}}}e^{-i{{\bf
k}}_{1}{\bm{\delta}}_{1}+i{{\bf k}}_{2}{\bm{\delta}}_{2}}\\\ e^{-i{{\bf
k}}_{3}{\bm{\delta}}}e^{-i{{\bf k}}_{1}{\bm{\delta}}_{2}+i{{\bf
k}}_{2}{\bm{\delta}}_{1}}&0\end{pmatrix}$ (46) $\displaystyle
a({\bm{\delta}}_{1},{\bm{\delta}}_{2},{\bm{\delta}})$
$\displaystyle=\begin{pmatrix}e^{i{{\bf k}}_{4}{\bm{\delta}}-i{{\bf
k}}_{1}{\bm{\delta}}_{1}+i{{\bf k}}_{2}{\bm{\delta}}_{2}}+e^{-i{{\bf
k}}_{3}{\bm{\delta}}-i{{\bf k}}_{1}{\bm{\delta}}_{2}+i{{\bf
k}}_{2}{\bm{\delta}}_{1}}&0\\\ 0&0\end{pmatrix}$ (47) $\displaystyle{\bf
c}_{{{\bf k}}}$ $\displaystyle=\begin{pmatrix}c_{{{\bf k}}}\\\
c^{\dagger}_{-{{\bf k}}}\end{pmatrix}$ (48)

We dropped the sum over ${{\bf r}}$ and ${{\bf k}}_{n}$ and the momentum
conserving $\delta({{\bf k}}_{a}+{{\bf k}}_{b}+...)$ in all equations above.
We see, that the chemical potential term acts as a magnetic field for the
pseudospins, making it possible to break $U(1)$ symmetry of the pseudospins,
as discussed before.

We continue to perform second order perturbation theory on this Hamiltonian at
finite $\mu$, to derive an effective hardcore boson and therefore pseudospin
Hamiltonian. We note here, that there is in general an arbitrariness in how
one can perform perturbation theory. Usually, perturbations around $H_{0}$
only contain the interaction terms $V$ of the Hamiltonian $H=H_{0}+V$. In
principle it is possible to add and subract a non-interacting term
$V_{\lambda}$ like so $H=H_{0}-V_{\lambda}+V+V_{\lambda}$. One can define a
new $H_{0}^{\prime}=H_{0}-V_{\lambda}$ and perturbation around it
$V^{\prime}=V+V_{\lambda}$. This way it is possible to describe symmetry
broken phases using a perturbative expansion, as long as $V_{\lambda}$ is much
smaller than $H_{0}$. The right $V_{\lambda}$ is chosen variationally. A
similar perturbative approach has been applied directly to the Hubbard model
at weak coupling in [20], where it was possible to optimize the introduced
fields via a self-consistency equation and where a finite d-wave
superconductive field has been found at non-zero doping.

One of the terms we vary between $H_{0}$ and $V_{\lambda}$ is the hopping term
$tc^{\dagger}_{{\bf r}}c_{{\bf r}+{\bm{\delta}}}+\text{h.c.}$ Without
variational perturbation theory (VPT), it is unclear if one has to put this
term in $H_{0}$ or $V^{\prime}$, as it is of order $t$ while still being a
free particle term.

Before going on to the results of VPT, we describe in the following the way
VPT is implemented. We first define $V_{\lambda=0}=0$, and
$H_{0}^{\prime}=H_{0}(\lambda)$ such that $H_{0}(\lambda=0)=H_{0}$. The zero
and first order terms are

$\displaystyle
H_{f}=\bra{\Psi(\lambda)}H_{0}(\lambda=0)+V_{K}+V_{A+\mu}\ket{\Psi(\lambda)}$
(49)

where $\ket{\Psi(\lambda)}$ is the groundstate of $H_{0}(\lambda)$. In second
order we must compute

$\displaystyle\bra{\Psi(\lambda)}V\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}V\ket{\Psi(\lambda)}$
(50)

where $V=V_{K}+V_{A+\mu}+V_{\lambda}$ and where we only consider connected
diagrams, such that we can drop the projection operators to excited states.
Due to the structure of all of the present vertices, the computation of the
second order terms reduce to the evaluation of

$\displaystyle\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}({\bf A}\cdot{\bf
c}^{\dagger}_{{{\bf k}}})({\bf B}\cdot{\bf c}_{{{\bf
k}}^{\prime}})\ket{\Psi(\lambda)}$ $\displaystyle=$
$\displaystyle\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}({\bf
A}\sigma_{x}\cdot{\bf c}_{-{{\bf k}}})({\bf B}\cdot{\bf c}_{{{\bf
k}}^{\prime}})\ket{\Psi(\lambda)}$ (51)

with some vectors ${\bf A}$ and ${\bf B}$ and where we used that

$\displaystyle{\bf c}^{\dagger}_{{{\bf k}}}$ $\displaystyle=\sigma_{x}{\bf
c}_{-{{\bf k}}}$ (52)

is the particle-hole transformation in the Nambu formalism. With the
commutation relation

$\displaystyle[H_{0},({\bf A}\sigma_{x}\cdot{\bf c}_{-{{\bf k}}})({\bf
B}\cdot{\bf c}_{{{\bf k}}^{\prime}})]$ $\displaystyle=(\sigma_{x}H_{{{\bf
k}}}\sigma_{x}{\bf A}\sigma_{x}\cdot{\bf c}_{-{{\bf k}}})({\bf B}\cdot{\bf
c}_{{{\bf k}}^{\prime}})$ $\displaystyle+({\bf A}\sigma_{x}\cdot{\bf
c}_{-{{\bf k}}})(\sigma_{x}H_{-{{\bf k}}^{\prime}}\sigma_{x}{\bf B}\cdot{\bf
c}_{{{\bf k}}^{\prime}})$ (53)

where $H_{{{\bf k}}}$ is defined by $H_{0}$

$\displaystyle H_{0}(\lambda)=\frac{1}{2}\sum_{{{\bf k}}}{\bf
c}^{\dagger}_{{{\bf k}}}H_{{{\bf k}}}(\lambda){\bf c}_{{{\bf k}}},$ (54)

we can show

$\displaystyle(1-E_{0}+H_{0})({\bf A}\sigma_{x}\cdot{\bf c}_{-{{\bf k}}})({\bf
B}\cdot{\bf c}_{{{\bf k}}^{\prime}})\ket{\Psi(\lambda)}=$
$\displaystyle(\sigma_{x}H_{{{\bf k}}}\sigma_{x}{\bf A}\sigma_{x}\cdot{\bf
c}_{-{{\bf k}}})({\bf B}\cdot{\bf c}_{{{\bf k}}^{\prime}})$
$\displaystyle+({\bf A}\sigma_{x}\cdot{\bf c}_{-{{\bf
k}}})(\sigma_{x}H_{-{{\bf k}}^{\prime}}\sigma_{x}{\bf B}\cdot{\bf c}_{{{\bf
k}}^{\prime}})$ $\displaystyle+({\bf A}\sigma_{x}\cdot{\bf c}_{-{{\bf
k}}})({\bf B}\cdot{\bf c}_{{{\bf k}}^{\prime}})\ket{\Psi(\lambda)}.$ (55)

We can define a two particle space

$\displaystyle({\bf A}\sigma_{x}\otimes{\bf B})({\bf c}_{-{{\bf
k}}}\otimes{\bf c}_{{{\bf k}}^{\prime}})\ket{\Psi_{1}}.$ (56)

Acting with $1-E_{0}+H_{0}$ on this state gives

$\displaystyle(1+\sigma_{x}H_{{{\bf
k}}}\sigma_{x}\otimes\mathbb{1}+\mathbb{1}\otimes\sigma_{x}H_{-{{\bf
k}}^{\prime}}\sigma_{x})$ $\displaystyle{\bf A}\sigma_{x}\otimes{\bf
B}\cdot{\bf c}_{-{{\bf k}}}\otimes{\bf c}_{{{\bf
k}}^{\prime}}\ket{\Psi(\lambda))}.$ (57)

So

$\displaystyle\frac{1}{E_{0}-H_{0}}{\bf A}\sigma_{x}\otimes{\bf B}\cdot{\bf
c}_{-{{\bf k}}}\otimes{\bf c}_{{{\bf k}}^{\prime}}\ket{\Psi(\lambda)}$
$\displaystyle=$ $\displaystyle\sum_{n=0}^{\infty}(1-E_{0}+H_{0})^{n}{\bf
A}\sigma_{x}\otimes{\bf B}\cdot{\bf c}_{-{{\bf k}}}\otimes{\bf c}_{{{\bf
k}}^{\prime}}\ket{\Psi(\lambda)}$ $\displaystyle=$
$\displaystyle\frac{1}{-\sigma_{x}H_{{{\bf
k}}}\sigma_{x}\otimes\mathbb{1}-\mathbb{1}\otimes\sigma_{x}H_{-{{\bf
k}}^{\prime}}\sigma_{x}}{\bf A}\sigma_{x}\otimes{\bf B}\cdot{\bf c}_{-{{\bf
k}}}\otimes{\bf c}_{{{\bf k}}^{\prime}}\ket{\Psi(\lambda)}.$ (58)

$\sigma_{x}H_{{{\bf
k}}}\sigma_{x}\otimes\mathbb{1}+\mathbb{1}\otimes\sigma_{x}H_{-{{\bf
k}}^{\prime}}\sigma_{x}$ can be seen as the free two particle Hamiltonian
propagating two holes. In general, vertices have a ${{\bf k}},{{\bf
k}}^{\prime}$ dependence

$\displaystyle{\bf A}\sigma_{x}\otimes{\bf B}\cdot{\bf c}_{-{{\bf
k}}}\otimes{\bf c}_{{{\bf k}}^{\prime}}={\bf c}_{-{{\bf k}}}(\sigma_{x}V({{\bf
k}},{{\bf k}}^{\prime})){\bf c}_{{{\bf k}}^{\prime}}.$ (59)

So in practice, to apply $\frac{1}{E_{0}-H_{0}}$ to this vertex,
$-\sigma_{x}H_{{{\bf
k}}}\sigma_{x}\otimes\mathbb{1}-\mathbb{1}\otimes\sigma_{x}H_{-{{\bf
k}}^{\prime}}\sigma_{x}$ is first diagonalized and $\sigma_{x}V({{\bf
k}},{{\bf k}}^{\prime})$ is transformed into a basis where the two particle
propagator is diagonal. Then it is straightforward to apply
$\frac{1}{E_{0}-H_{0}}$ and the transformation is reverted. Afterwards we
arrive at

$\displaystyle\frac{1}{E_{0}-H_{0}}{\bf c}_{-{{\bf k}}}(\sigma_{x}V({{\bf
k}},{{\bf k}}^{\prime})){\bf c}_{{{\bf k}}^{\prime}}\ket{\Psi(\lambda)}$
$\displaystyle=$ $\displaystyle{\bf c}_{-{{\bf k}}}(\sigma_{x}\tilde{V}({{\bf
k}},{{\bf k}}^{\prime})){\bf c}_{{{\bf k}}^{\prime}}\ket{\Psi(\lambda)}.$ (60)

Then we only need to contract expectation values of the following form

$\displaystyle\bra{\Psi(\lambda)}{\bf c}^{\dagger}_{{{\bf k}}_{1}}(V({{\bf
k}}_{1},{{\bf k}}_{2})\sigma_{x}){\bf c}^{\dagger}_{-{{\bf k}}_{2}}{\bf
c}_{-{{\bf k}}_{3}}(\sigma_{x}\tilde{V}({{\bf k}}_{3},{{\bf k}}_{4})){\bf
c}_{{{\bf k}}_{4}}\ket{\Psi(\lambda)}$ (61)

which can be done using Wick’s theorem. Leaving only connected diagrams, this
results with Einstein summation in

$\displaystyle(V({{\bf k}}_{4},{{\bf
k}}_{3})\sigma_{x})_{ab}(\sigma_{x}\tilde{V}({{\bf k}}_{3},{{\bf
k}}_{4}))_{cd}$ $\displaystyle\rho({{\bf
k}}_{4})^{ad}\delta_{k_{1},k_{4}}\rho(-{{\bf
k}}_{3})^{bc}\delta_{k_{2},k_{3}}$ $\displaystyle-$ $\displaystyle(V({{\bf
k}}_{3},{{\bf k}}_{4})\sigma_{x})_{ab}(\sigma_{x}\tilde{V}({{\bf k}}_{3},{{\bf
k}}_{4}))_{cd}$ $\displaystyle\rho(-{{\bf
k}}_{3})^{ac}\delta_{k_{1},-k_{3}}\rho({{\bf
k}}_{4})^{bd}\delta_{-k_{2},k_{4}}$ (62)

where $\rho({{\bf k}})=\langle{\bf c}^{\dagger}_{{{\bf k}}}{\bf c}_{{{\bf
k}}}\rangle$ are the one particle density matrices of $\frac{1}{2}H_{{\bf
k}}$. In the following, $\rho$ is computed at a temperature of $T=10^{-4}$.
Note that for interacting vertices $V$ contains $b$ operators, resulting in an
effective hardcore boson Hamiltonian, with terms of the form

$\displaystyle\sum_{{\bf k}}\epsilon_{f}({{\bf k}})b^{\dagger}_{{\bf
k}}b_{{\bf k}}$ (63)

and interactions

$\displaystyle\sum_{{\bf k}}v_{b}({{\bf k}},{{\bf k}}^{\prime},{\bf
q})b^{\dagger}_{{\bf k}}b_{{{\bf k}}+{\bf q}}b^{\dagger}_{{{\bf
k}}^{\prime}}b_{{{\bf k}}^{\prime}-{\bf q}}.$ (64)

Both terms are Fourier transformed back to real space, after which hardcore
bosonic commutation relations are applied to result in the simplest form of a
given expression. For example

$\displaystyle b^{\dagger}_{{{\bf r}}_{1}}b_{{{\bf r}}_{2}}b^{\dagger}_{{{\bf
r}}_{1}}b_{{{\bf r}}_{3}}=0$ (65)

and

$\displaystyle b^{\dagger}_{{{\bf r}}_{1}}b_{{{\bf r}}_{2}}b^{\dagger}_{{{\bf
r}}_{2}}b_{{{\bf r}}_{2}}=b^{\dagger}_{{{\bf r}}_{1}}b_{{{\bf r}}_{2}}.$ (66)

The resulting Hamiltonian is then transformed to a pseudospin Hamiltonian
using Eq. (35). To compute the groundstate energy and expectation values,
spins are replaced by classical spins of length $1/2$. The spins are chosen to
minimize the classical groundstate energy. This approximation is the first
step of spin wave theory and is valid in an ordered phase with low densities
of excitations around the classical groundstate. All of the above steps are
done algorithmically. Integrations over the BZ are done on a regular grid and
the minimization of the classical spins energy is implemented using the
Pytorch library [22] and an Adam optimizer [23].

###  IV.1 Pure Hubbard model at $U=8$

We apply the above described algorithms and the local unitary transformation
to study the pure Hubbard model without $t^{\prime}$ and investigate first the
variation of $t$ in $H_{0}$, so

$\displaystyle H_{0}(\lambda)$ $\displaystyle=t(1-\lambda_{1})\sum_{{\bf
r},{\bm{\delta}}}c^{\dagger}_{{\bf r}}c_{{\bf
r}+{\bm{\delta}}}+\text{h.c.}+\sum_{{\bf
r}}\left(-\frac{U}{2}+\mu\right)n_{{\bf r}},$ (67)

such that when $\lambda_{1}=1$, $H_{0}$ contains no hopping terms and at
$\lambda_{1}=0$ it does. Therefore at $\lambda_{1}=1$, the vertex
$V_{\lambda}=H_{0}(\lambda=0)-H_{0}(\lambda)$ contains the hopping whereas at
$\lambda_{1}=0$ it does not, where $V_{\lambda}=0$. In the following we set
$t=-1$. We consider the $U=8$ strong interaction regime relevant to cuprates.
Computing the free energy by minimizing the classical pseudospin groundstate
energy we find at $\mu=0$ an optimum at $\lambda_{1}=1$, see Fig. 2. This
optimum exactly replicates the Heisenberg model Eq. (39). Even though this is
a maximum, this groundstate is chosen by the system, since there the
sensitivity (first derivative or second derivative of the energy) to changes
in $H_{0}$ is minimized, see the principle of minimal sensitivity [24]. At
$\mu=1.5$ (hole doping regime due to the sign convention) this maximum
persists, see Fig. 3. So in the following, we fix $\lambda_{1}=1$. Doing so,
we find that at any $0<\mu<U/2$ and $T=0$, the system is unable to fill in
holes, see Eq. (23), since $\ket{\Psi(\lambda)}$ remains a state fully filled
with $c$ fermions. We introduce a variational field into $H_{0}$ and
$V_{\lambda}$, which at the same time breaks $C_{4}$ symmetry and is able to
introduce holes into the system

Figure 2: Groundstate energy over $\lambda_{1}$ for $\mu=0$ and $U=8$, which
shows an optimum at $\lambda_{1}=1$. Integrations for interaction vertices
were done on a 16x16 grid BZ. For all other integrations a 64x64 grid was
used. Figure 3: Groundstate energy over $\lambda_{1}$ for $\mu=1.5$ and $U=8$.
The optimum adiabatically connected to the AFM groundstate at $\mu=0$
persists. At $\lambda_{1}=0.4$ a phase transition to a long range pseudospin
stripe groundstate appears, where no optimum can be found. The resulting
minimum at the phase transition is discontinuous and therefore not stable.
Integrations for interaction vertices were done on a 16x16 grid BZ. For all
other integrations a 64x64 grid was used.

$\displaystyle V_{\lambda_{2}}$ $\displaystyle=\sum_{{{\bf
r}},{\bm{\delta}}}\lambda_{2}e^{i({\bf Q}{{\bf
r}}+\pi\delta_{y})}c^{\dagger}_{{\bf r}}c^{\dagger}_{{\bf
r}+{\bm{\delta}}}+\text{h.c.}$ (68)

with ${\bf Q}=(\pi,\pi)$. We find that this field stabilizes d-wave
superconductivity. We note that in the pure Hubbard model, we do not only
expect to see superconductivity but also for example magnetic stripe orders
[25], which are preferred by the system by $0.01t$ over superconductivity at a
doping of $12.5\%$ [26]. However, even if a different state is lower in
energy, we still expect the system to prefer d-wave superconductivity over
pure AFM and therefore to see an optimum appearing at a finite $\lambda_{2}$,
which is exactly what we find, see Fig. 4. A mexican hat potential develops at
a critical $\mu\approx 0.7$ and remains at any finite doping. This critical
chemical potential is comparable to the one found in [27] for the same system
parameters. We suspect that at higher dopings, the approximation of classical
spins no longer holds and quantum fluctuations lead to the disappearance of
this phase. Nevertheless, the onset of superconductivity is well described by
the classical groundstate of a spin Hamiltonian. At $\mu=1.4$, the optimal
$\lambda_{2}=0.37$ corresponds to a doping of $2.4\%$ (the expectation value
was computed up to first order in perturbation) and d-wave superconductivity
is preferred over pure AFM by an energy of $0.011t$, so $163.4\,\text{K}$ at
$t=-1\,\text{eV}$. We find as expected, that a finite chemical potential acts
as a magnetic field for the pseudospins, resulting in a $U(1)$ symmetry
breaking. The expectation value $|\langle\Delta_{{{\bf r}},{{\bf
r}+{\bm{\delta}}}_{x}}\rangle|$ is $1.9\cdot 10^{-3}$ in zeroth order on a
bond in the $x$ direction, while $\langle\Delta_{{{\bf r}},{{\bf
r}+{\bm{\delta}}}_{y}}\rangle=-\langle\Delta_{{{\bf r}},{{\bf
r}+{\bm{\delta}}}_{x}}\rangle$, realizing uniform d-wave superconductivity. We
have not yet found an order parameter implementing s-wave or stripe order. We
note, that the presented method severely undershoots groundstate energies at
higher dopings [26], which we attribute to the missing energy from the
groundstate quantum fluctutations of the pseudospins.

To test the convergence, the same calculation, which was done on a 16x16 BZ
grid for the integration of interaction vertices, was carried out for
$\mu=1.4$ on a 32x32 BZ grid. The average error between the 16x16 and the
32x32 calculation is $5\cdot 10^{-6}$, meaning that the results are well
converged, even at a temperature of $T=10^{-4}$ for the density matrix
computations. This signifies, that in all integrations carried out here, no
significant singularities have been integrated over. This is the result of
modes being integrated out, that are with $U/2$ strongly gapped. This
suggests, that the perturbative series described in this paper is controlled
and convergent for large $U$ and that higher orders contribute less than the
ones we accounted for here.

$\mu=0$ |  $\mu=0.5$
---|---
|

$\mu=1.0$ |  $\mu=1.4$
---|---
|

$\mu=2.0$

Figure 4: Groundstate energy over $\lambda_{2}$ for different $\mu$. We see
that at a critical $\mu$ a new minimum at a finite $\lambda_{2}$ is developed,
leading to nonzero hole densities and d-wave superconductivity. Figure 5: Hole
densities $1-n$ and expectation values of $\Delta$ (singlet x/y) over $\mu$ at
optimized $\lambda_{2}$ for $U=8$ in the pure Hubbard model. Since the $U(1)$
phase of both singlets is arbitrary, singlet x has been rotated to be purely
real and positive. The rotation has also been applied to singlet y. The
imaginary part of singlet y is numerically small. We see a phase transition at
$\mu\approx 0.7$ to a d-wave superconducting phase. At smaller $\mu$, the hole
densities remain at $0$. For every $\mu$, $\lambda_{2}$ has been optimized to
minimize the free energy up to an error of $10^{-4}$.

In Fig. 5, hole densities and singlet Cooper pair expectation values on
nearest neighbor bonds in x and y direction are shown over chemical potential
$\mu$, where at every $\mu$ a $\lambda_{2}$ is found to minimize the free
energy up to an error of $10^{-4}$. Hole densities were computed up to first
and the singlet expectation values up to zero order in perturbation theory. We
find a critical $\mu\approx 0.7$ for $U=8$, where at the same time holes start
to fill into the system while d-wave superconductivity develops. At smaller
$\mu$, the system remains at half-filling.

Figure 6: Comparison of groundstate energy over $\lambda_{2}$ for the full and
simplified model at $\mu=1.4$. The terms left out still contribute
significantly. Nevertheless, the included terms are enough to observe the
development of an energy minimum at a non-zero $\lambda_{2}$. See Appendix A
for a description of the simple model.

We have identified the most crucial terms to see the minimum in $\lambda_{2}$
at a non-zero value and therefore the development of superconductivity. This
simplified model is described in Appendix A. A comparison of groundstate
energies between the full algorithmically build model and this simple model is
shown in Fig. 6. We see, that more terms contribute significantly, so this
simple model is just a rough approximation of the full one.

###  IV.2 Hubbard model at $U=8$ with $t^{\prime}=-0.3t$

Figure 7: Groundstate energy over order parameter $\lambda_{2}$ for TBCO/YBCO
hopping strengths $t=-1$, $t^{\prime}=0.3$ and at $U=8$ [28]. The minimum at a
finite $\lambda_{2}$ remains, signifying the development of superconductivity.
Figure 8: Hole densities $1-n$ and expectation values of $\Delta$ (singlet
x/y) over $\mu$ at optimized $\lambda_{2}$ for $U=8$ in the Hubbard model with
$t^{\prime}=0.3$, at TBCO/YBCO parameters [28]. The development of d-wave
superconductivity is observed. Figure 9: Energy difference between the
groundstate energy at optimal $\lambda_{2}$ and $\lambda_{2}=0$ in Kelvin (for
$t=0.43\,\text{eV}$ [28]) and singlet pairing expectation values over hole
doping for TBCO/YBCO parameters.

We apply the same methods described in the previous sections to TBCO/YBCO
hopping parameters with $t^{\prime}=-0.3t$ being the next nearest neighbor
hopping amplitude [28]. $U$ is set to $8$. The $\lambda_{2}$ groundstate
energy potential is shown in Fig. 7. The minimum in the free energy remains at
a finite value of $\lambda_{2}=0.33$. In Fig. 8, hole densities and singlet
pairing expectation values are shown. The latter correspond to the development
of d-wave superconductivity. Fig. 9 shows the singlet expectation value with
$t^{\prime}=-0.3t$ and the groundstate energy difference between optimal
$\lambda_{2}$ and $\lambda_{2}=0$ in Kelvin at $t=0.43\,\text{eV}$ [28]. We
see that this newly developed energy scale is on the order of $T_{c}$ in these
compounds and that the energy difference scales proportionally with the
superconducting order parameter. Energy differences are likely overestimated
as we expect that quantum fluctuations in the pseudospins weaken the states
found at $\lambda_{2}>0$.

##  V Discussion

In this paper, we propose a new theoretical description of high temperature
superconductivity, based on the application of a unitary transformation to the
Hubbard model, which maps the Hubbard interaction to a single particle term.
This unitary transformation can be applied exactly. The resulting Hamiltonian
contains, unlike the t-J model, unconstrained fermions and can be transformed
into one with a single fermion species coupled to hardcore bosons/pseudospins,
which suggests that this transformation implements a spin charge separation in
any dimension, where the fermions carry charge and the pseudospins spin
information. The Hamiltonian is solved by variational second order
perturbation theory (VPT), resulting in an effective pseudospin Hamiltonian.
The problem of fermions coupled to spins generating an effective spin-spin
interaction is reminiscent of the RKKY mechanism [29]. The breakdown of the
phases found in this paper might be equivalent to a Kondo-RKKY transition. To
compute the groundstate energy, the spins are replaced by classical spins and
their classical energy is minimized. We find that in zero order, the
introduction of a chemical potential $\mu$ acts as a magnetic field for the
pseudospins, breaking $U(1)$ symmetry. Mapping back to the original system,
this $U(1)$ symmetry breaking corresponds to the $U(1)$ symmetry breaking
found in a superconducting state, assuming the superconducting state is the
one minimizing the free energy. An order parameter is proposed which at the
same time fills in holes into the system and breaks $C_{4}$ symmetry. The free
energy is minimized at a finite value of this order parameter, resulting in a
finite expectation value of the d-wave superconductor pairing operator. This
perturbative approach is possible, because the fermionic modes that are
integrated out are strongly gapped with a gap of $U/2$. We therefore encourage
further studies of Hamiltonians (21) or (41) with non-perturbative
calculations (for example DMFT [11, 12, 13] or DMRG [30]), which should even
increase the accuracy of the presented results. Future work also involves
exploring similar unitary transformations like (3), for example ones that map
a whole cluster of sites to single particle terms. It can also be interesting
to vary different and more variational parameters than just one in VPT, like
the chemical potential. It is also possible to connect auxiliary bath sites
with VPT. Exploring the influence of quantum fluctuations around the classical
spin groundstate found in this paper should also give more insights. This is
especially important for studying magnetism in the models found in this paper,
since classical spins can not describe paramagnetic states at $T=0$, which is
why we did not explore magnetism here. In light of numerous potential paths
for future work, this paper already suggests that the emergence of
high-$T_{c}$ superconductivity is governed by the ground state of a classical
spin system, showcasing promising directions for further research.

## Acknowledgements

I am thankful for Robin Scholle, Lukas Debbeler, Steffen Bollmann, Silvia
Neri, Pietro Maria Bonetti, Elio König, Laura Classen, Nikolaos Parthenios,
Sida Tian, Raffaele Mazzilli, Henrik Müller-Groeling, Paulo Forni, Janika
Reichstetter, Andreas Schnyder and Walter Metzner for the many enlightening
discussions. I would also like to thank Thomas Schäfer, Andreas Schnyder and
Walter Metzner for reviewing this paper. Kirill Alpin is funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – TRR 360 –
492547816.

##  APPENDIX A Simplest effective pseudospin model

We have identified the most crucial terms of the effective pseudospin
Hamiltonian in second order perturbation theory to observe the development of
d-wave superconductivity. The spin Hamiltonian takes on this form

$\displaystyle H_{S}$
$\displaystyle=\bra{\Psi(\lambda)}H_{0}(\lambda=0)\ket{\Psi(\lambda)}$ (69)
$\displaystyle-\bra{\Psi(\lambda=0)}H_{0}(\lambda=0)\ket{\Psi(\lambda=0)}$
(70) $\displaystyle+2\mu\sum_{{\bf r}}\left(\frac{1}{2}-S^{z}_{{\bf
r}}\right)(1-\langle n_{{\bf r}}\rangle)$ (71)
$\displaystyle+X_{1}+X_{2}t^{2}\sum_{{\bf r},{\bm{\delta}}}(1-4{\bf S}_{{\bf
r}}\cdot{\bf S}_{{\bf r}+{\bm{\delta}}})$ (72)
$\displaystyle+(-2X_{3}\mu+4X_{4}\mu^{2})\sum_{{\bf
r}}\left(\frac{1}{2}-S^{z}_{{\bf r}}\right)$ (73)
$\displaystyle+4X_{5}\mu^{2}\sum_{{\bf
r},{\bm{\delta}}}\left(\frac{1}{2}-S^{z}_{{\bf
r}}\right)\left(\frac{1}{2}-S^{z}_{{\bf r}+{\bm{\delta}}}\right)$ (74)
$\displaystyle+X_{6}t^{2}\sum_{{\bf r}_{1},{\bf r}_{2},{\bf r}_{3}}W_{{\bf
r}_{1},{\bf r}_{2},{\bf r}_{3}}$ (75)

with

$\displaystyle X_{1}$
$\displaystyle=\bra{\Psi(\lambda)}V_{l}\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}V_{l}\ket{\Psi(\lambda)}$
(76) $\displaystyle X_{2}$ $\displaystyle=\bra{\Psi(\lambda)}c^{\dagger}_{{\bf
r}}c^{\dagger}_{{\bf
r}+{\bm{\delta}}}\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}c_{{\bf
r}+{\bm{\delta}}}c_{{\bf r}}\ket{\Psi(\lambda)}$ (77) $\displaystyle X_{3}$
$\displaystyle=\bra{\Psi(\lambda)}V_{l}\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}n_{{\bf
r}}\ket{\Psi(\lambda)}+\text{h.c.}$ (78) $\displaystyle X_{4}$
$\displaystyle=\bra{\Psi(\lambda)}n_{{\bf
r}}\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}n_{{\bf r}}\ket{\Psi(\lambda)}$ (79)
$\displaystyle X_{5}$ $\displaystyle=\bra{\Psi(\lambda)}n_{{\bf
r}}\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}n_{{{\bf
r}+{\bm{\delta}}}}\ket{\Psi(\lambda)}$ (80) $\displaystyle X_{6}$
$\displaystyle=\bra{\Psi(\lambda)}c^{\dagger}_{{\bf r}_{1}}c^{\dagger}_{{\bf
r}_{2}}\frac{1}{E_{0}(\lambda)-H_{0}(\lambda)}c_{{\bf r}_{2}}c_{{\bf
r}_{3}}\ket{\Psi(\lambda)}$ (81)

and

$\displaystyle W_{{\bf r}_{1},{\bf r}_{2},{\bf r}_{3}}$ $\displaystyle=K_{{\bf
r}_{1},{\bf r}_{2}}K^{\dagger}_{{\bf r}_{3},{\bf r}_{2}}$ (82)
$\displaystyle=-\frac{1}{4}+{\bf S}_{{\bf r}_{1}}\cdot{\bf S}_{{\bf
r}_{2}}+{\bf S}_{{\bf r}_{2}}\cdot{\bf S}_{{\bf r}_{3}}$ $\displaystyle-{\bf
S}_{{\bf r}_{1}}\cdot{\bf S}_{{\bf
r}_{3}}-2i\sum_{nmp}\epsilon_{nmp}S^{n}_{{\bf r}_{3}}S^{m}_{{\bf
r}_{2}}S^{p}_{{\bf r}_{1}}.$ (83)

$\epsilon_{nmp}$ is the Levi-Civita tensor and the sum $\sum_{{\bf r}_{1},{\bf
r}_{2},{\bf r}_{3}}$ runs over all nearest-neighbor L shaped 3 spin clusters,
where ${\bf r}_{1}$ and ${\bf r}_{3}$ are the furthest apart. All $X_{n}$ are
isotropic in ${\bf r}$ and ${\bm{\delta}}$, so it does not matter where they
are evaluated. The term (70) is not necessary and is applied to measure
relative energies. You arrive at the Heisenberg term (72) by $K_{{{\bf
r}},{\bm{\delta}}}K^{\dagger}_{{{\bf r}},{\bm{\delta}}}=1-4{\bf S}_{{\bf
r}}\cdot{\bf S}_{{\bf r}+{\bm{\delta}}}$.

## References

  * Anderson [1987] P. W. Anderson, Science 235, 1196 (1987), https://www.science.org/doi/pdf/10.1126/science.235.4793.1196 .
  * Tsuei  _et al._ [1994] C. C. Tsuei, J. R. Kirtley, C. C. Chi, L. S. Yu-Jahnes, A. Gupta, T. Shaw, J. Z. Sun, and M. B. Ketchen, Phys. Rev. Lett. 73, 593 (1994).
  * Qin  _et al._ [2022] M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, Annual Review of Condensed Matter Physics 13, 275 (2022).
  * Arovas  _et al._ [2022] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, Annual Review of Condensed Matter Physics 13, 239 (2022), https://doi.org/10.1146/annurev-conmatphys-031620-102024 .
  * Scalapino [2012] D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Chowdhury  _et al._ [2022] D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, Reviews of Modern Physics 94, 035004 (2022).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * MacDonald  _et al._ [1988] A. H. MacDonald, S. M. Girvin, and D. Yoshioka, Phys. Rev. B 37, 9753 (1988).
  * Pairault  _et al._ [1998] S. Pairault, D. Sénéchal, and A.-M. Tremblay, Physical review letters 80, 5389 (1998).
  * Pairault  _et al._ [2000] S. Pairault, D. Senechal, and A.-M. Tremblay, The European Physical Journal B-Condensed Matter and Complex Systems 16, 85 (2000).
  * Metzner and Vollhardt [1989] W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * Georges and Kotliar [1992] A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Potthoff and Balzer [2007] M. Potthoff and M. Balzer, Physical Review B 75, 125112 (2007).
  * Sénéchal  _et al._ [2002] D. Sénéchal, D. Perez, and D. Plouffe, Physical Review B 66, 075129 (2002).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev. Mod. Phys. 77, 1027 (2005).
  * Rohringer  _et al._ [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Rev. Mod. Phys. 90, 025003 (2018).
  * Šimkovic  _et al._ [2021] F. Šimkovic, Y. Deng, and E. Kozik, Phys. Rev. B 104, L020507 (2021).
  * Šimkovic  _et al._ [2016] F. Šimkovic, X.-W. Liu, Y. Deng, and E. Kozik, Phys. Rev. B 94, 085106 (2016).
  * Neumayr and Metzner [2003] A. Neumayr and W. Metzner, Phys. Rev. B 67, 035112 (2003).
  * Khaliullin [1994] G. Khaliullin, Journal of Experimental and Theoretical Physics 79, 473 (1994).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in  _Advances in Neural Information Processing Systems 32_ (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980 (2014).
  * Stevenson [1984] P. Stevenson, Physical Review D 30, 1712 (1984).
  * Qin  _et al._ [2020] M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 10, 031016 (2020).
  * Zheng  _et al._ [2017] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science 358, 1155 (2017), https://www.science.org/doi/pdf/10.1126/science.aam7127 .
  * Akiyama  _et al._ [2022] S. Akiyama, Y. Kuramashi, and T. Yamashita, Progress of Theoretical and Experimental Physics 2022, 023I01 (2022), https://academic.oup.com/ptep/article-pdf/2022/2/023I01/42931059/ptac014.pdf .
  * Pavarini  _et al._ [2001] E. Pavarini, I. Dasgupta, T. Saha-Dasgupta, O. Jepsen, and O. Andersen, Physical review letters 87, 047003 (2001).
  * Ruderman and Kittel [1954] M. A. Ruderman and C. Kittel, Phys. Rev. 96, 99 (1954).
  * White [1992] S. R. White, Phys. Rev. Lett. 69, 2863 (1992).