# Thermodynamics of an Exactly Solvable Model for Superconductivity in a Doped
Mott Insulator

Jinchao Zhao  Luke Yeo  Edwin W. Huang  Philip W. Phillips  Department of
Physics and Institute of Condensed Matter Theory, University of Illinois at
Urbana-Champaign, Urbana, 61801 Illinois, USA

###### Abstract

Computing superconducting properties starting from an exactly solvable model
for a doped Mott insulator stands as a grand challenge. We have recently shown
that this can be done starting from the Hatsugai-Kohmoto (HK) model which can
be understood generally as the minimal model that breaks the non-local
$\mathbb{Z}_{2}$ symmetry of a Fermi liquid, thereby constituting a new
quartic fixed point for Mott physics [Phillips et al., Nature Physics 16, 1175
(2020); Huang et al., Nature Physics (2022)]. In the current work, we compute
the thermodynamics, condensation energy, and electronic properties such as the
NMR relaxation rate $1/T_{1}$ and ultrasonic attenuation rate. Key differences
arise with the standard BCS analysis from a Fermi liquid: 1) the free energy
exhibits a local minimum at $T_{p}$ where the pairing gap turns on
discontinuously above a critical value of the repulsive HK interaction,
thereby indicating a first-order transition, 2) a tri-critical point emerges,
thereby demarcating the boundary between the standard second-order
superconducting transition and the novel first-order regime, 3) Mottness
changes the sign of the quartic coefficient in the Landau-Ginzburg free-energy
fuctional relative to that in BCS, 4) as this obtains in the strongly
interacting regime, it is Mott physics that underlies the generic first-order
transition, 5) the condensation energy exceeds that in BCS theory suggesting
that multiple Mott bands might be a way of enhancing superconducting, 6) the
heat-capacity jump is non-universal and increases with the Mott scale, 7)
Mottness destroys the Hebel-Slichter peak in NMR, and 8) Mottness enhances the
fall-off of the ultrasonic attenuation at the pairing temperature $T_{p}$. As
several of these properties are observed in the cuprates, our analysis here
points a way forward in computing superconducting properties of strongly
correlated electron matter.

##  I Introduction

A truly remarkable feature of superconductivity in elemental metals is that
although the superconducting order parameter, $\Psi=\Delta e^{i\theta}$, has
two components, a phase, $\theta$ and an amplitude, $\Delta$, both turn on at
the same temperature. This behaviour is captured by the mean-field treatment
of BCSBardeen  _et al._ (1957) which predicts a second-order phase transition
in which superconducting fluctuations of the pair amplitude satisfy the
Ginzburg criterion, $\delta(\Delta)^{2}/\Delta\ll 1$. Should this criterion be
met, pairing and phase coherence are synonymous. From a technical standpoint,
satisfying the Ginzburg criterion is quite surprising as the latter is
governed by a divergence in the pair susceptibility while the former stems
from a solution to a mean-field integral equation. Experimentally, no more
than a 1K difference between these temperatures in BCS materials is observed.
Hence, the Ginzburg criterion is really a statement about the accuracy of
mean-field theory.

It is well known that the cuprates violate the Ginzburg criterion or
equivalently the BCS dictum that pairing without long-range order is
impossibleCorson  _et al._ (1999); Bollinger  _et al._ (2011); Uemura  _et
al._ (1989); Emery and Kivelson (1995). While this trend was thought to vanish
in the overdoped regime, recent spectroscopic measurementsHe  _et al._ (2020)
on the bilayer cuprate (Pb,Bi)2Sr2CaCu2O8+δ (Bi-2212) have verified that even
in such samples remnants of a particle-hole symmetric superconducting gap
persist up to $T_{\rm pair}=86K$ where the superconducting transition is
$T_{c}=63K$. The focus on a particle-hole symmetric gap as the signature of
the pairing gap is designed to disentangle superconductivity-related gaps from
the wider range of phenomena associated with the pseudogapNorman  _et al._
(2005) that need notMiller  _et al._ (2017); Ghosh  _et al._ (2019); Anderson
and Haldane (2001) have this symmetry. Within the cuprate family, Bi-2212 and
YBa2Cu3O7-δ (YBCO) have the highest ratio of $T_{\rm pair}/T_{c}$ of 1.6 while
in (La,Sr)2CuO4 (LSCO) the ratio is 1.2. Although 2D disordered materials are
expected to have a wide range where thermal superconducting fluctuations
obtainJaeger  _et al._ (1989); Yazdani and Kapitulnik (1995), the efficient
cause of the disconnect between phase coherence and the gap turn-on remains
unsettled in the cuprates. For example, the discrepancy has been attributedHe
_et al._ (2020) to the existence of a flat band at $(\pi,0)$. Further, it is
unclear how to think about a pairing gap without a phase transition. The
question arises, what is the order of the transition in which pairing and
phase coherence are decoupled? It is questions of this type that we address in
this work.

In trying to understand the source of the discrepancy between pairing and
phase coherence, it is worth cataloguing other instances of deviations from
BCS superconductivity in the cuprates. Two features stand out: 1) the color
changeMolegraaf  _et al._ (2002) and 2) a violation of the Glover-Ferrell-
TinkhamFerrell and Glover (1958); Tinkham and Ferrell (1959) (GFT) sum rule.
Regarding the latter, in a standard BCS superconductor, condensation leads to
loss of spectral weight at energy scales no more than ten times the pairing
energy. Not so in the cuprates. Bontemps and colleaguesSantander-Syro  _et
al._ (2003) have directly observed that in underdoped (but not overdoped)
Bi-2212, the Glover-Ferrel-Tinkham sum rule is violated and the optical
conductivity must be integrated to 20,000$cm^{-1}$ to recover the spectral
weight lost upon condensation into the superconducting state. Similarly,
Rübhaussen, et. alRübhausen  _et al._ (2001) and othersHomes  _et al._ (2004)
have shown that changes in the optical conductivity occur at energies 3eV
(roughly 100$\Delta$ where $\Delta$ is the maximum superconducting gap) away
from the Fermi energy at $T_{c}$. Finally, van der Marel and
colleaguesMolegraaf  _et al._ (2002) have seen an acceleration in the
depletion of the high energy spectral weight accompanied with a compensating
increase in the low-energy spectral weight at and below the superconducting
transition. Specifically the integrated weight of the optical conductivity
over the lower Hubbard band scales ($<1eV$) increases below the
superconducting temperature, whereas the high-energy component ($[1eV,2eV]$)
decreases. Since the integrated weight determines the number of charge
carriers, the color change indicates that high-energy scales contribute to the
superfluid density in contrast to the standard BCS picture. Consequently, the
superfluid density in the cuprates is derived not just from low energy
physics. In essence, it involves UV-IR mixing.

What all of this seems to indicate is that the departures from BCS
superconductivity in the cuprates are tethered to the Mott state. What is
difficult then is to solve a reasonable model for a Mott insulator which
captures this range of deviations from the BCS paradigm. It is in attempting
to answer this array of questions that we have focusedPhillips  _et al._
(2020) on the Hatsugai-Kohmoto (HK)Hatsugai and Kohmoto (1992) model, an
exactly solvable model for a Mott insulator. The HK model,

$\displaystyle
H_{HK}=\sum_{{\mathbf{k}}\sigma}\xi_{{\mathbf{k}}}n_{{\mathbf{k}}\sigma}+U\sum_{{\mathbf{k}}}n_{{\mathbf{k}}\uparrow}n_{{\mathbf{k}}\downarrow},$
(1)

is essentially the Hubbard model in momentum space. As we pointed out
previouslyHuang  _et al._ (2022), this model is important because it
represents the simplest way of breaking the hidden $\mathbb{Z}_{2}$ symmetry
of a Fermi liquid. We illustrate this from the basic Hamiltonian for a Fermi
liquid

$\displaystyle H_{\rm
FL}=\sum_{{\mathbf{p}}}\psi_{{\mathbf{p}}}^{\dagger}(\epsilon_{{\mathbf{p}}}-\epsilon_{F})\tau_{3}\psi_{{\mathbf{p}}}+\cdots.$
(2)

Here
$\psi^{\dagger}_{{\mathbf{p}}}=(c^{\dagger}_{{\mathbf{p}}\uparrow},c_{-{\mathbf{p}}\downarrow})$
and $\tau_{3}$ is the standard z-Pauli matrix. For electrons at the Fermi
surface, $\epsilon_{{\mathbf{p}}}=\epsilon_{F}$, this Hamiltonian obeys the
symmetry $n_{{\mathbf{k}}\uparrow}\rightarrow-n_{{\mathbf{k}}\uparrow}$ where
only one of the spin currents changes sign. Consequently, the interaction term
of the form in the HK model maximally breaks this symmetry. Based on this, we
showed previouslyHuang  _et al._ (2022) that the HK model represents a fixed
point for Mott physics that even encompasses Hubbardology. A straightforward
Fourier transform of the Hubbard on-site interaction reveals that it contains
the HK interaction term. As we have shown previouslyHuang  _et al._ (2022), it
is this term that is most relevant in the renormalization sense and the only
one that maximally breaks the hidden $Z_{2}$ symmetry of a Fermi liquid. For
these reasons, we have focused on revealing its superconducting
propertiesPhillips  _et al._ (2020). To accomplish this, we appended to the HK
Hamiltonian, which only describes a doped Mott insulator, a pairing term, in
analogy with the Cooper programBardeen  _et al._ (1957) in which a pairing
term was added to a Fermi liquid, to investigate if a corresponding pairing
instability obtains. Indeed it didPhillips  _et al._ (2020) and hence we found
a computationally tractable starting point for superconductivity and Mottness.
While our previous work revealed $T=0$ properties of the HK model appended
with a pairing term with couplinng constant, $g$, none of the thermodynamics
were obtained. Nonetheless, several non-BCS properties were apparent: 1)
$\lim_{g\rightarrow 0}\frac{2\Delta}{T_{c}}\rightarrow\infty$ rather than the
BCS ratio of $3.52$, 2) Composite quasiparticle excitations consisting of
doublons and holons rather than the standard particle-hole excitations of BCS,
and 3) a suppression of the superfluid density relative to that of BCS.

Because the model we explore is highly amenable to computation, we can with
certainty catalog how the finite temperature properties derived from Mottness
lead to substantial deviations from BCS theory. First, we establish that the
turn-on of the gap and the divergence of the susceptibility occur at distinct
temperatures. The former appears to obtain at a first-order transition while
the latter tends to a global second-order transition of the superconducting
state. We show that this effect vanishes when $U=0$, thereby making this a
true consequence of Mottness. We make contact with earlier results on multi-
band superconductors that found a first-order transitionSellin and Babaev
(2016); Huang and Lin (2020); Padilha and Continentino (2009); Sarasua and
Continentino (2002). We trace the first-order nature of the pairing transition
to a singly-occupied holon band that arises purely from Mott physics. Second,
we compute the heat capacity jump at $T_{c}$ as well as the condensation
energy. Unlike BCS theory, we find that $E_{\rm cond}/\Delta$ is stronger than
in BCS theory. Finally, we compute the ultrasonic attenuation as well as
$T_{1}$ near the superconducting transition. We are able to show that the
Hebel-SlichterHebel and Slichter (1957) peak in BCS theory vanishes in the
strongly correlated limit as seen widely in the cupratesOhsugi  _et al._
(1994); Yoshimura  _et al._ (1989); Walstedt  _et al._ (1994). Although the
absence of this peak in the cupratesZoli (1991); Cavanagh and Powell (2021)
has been attributed to spin fluctuations, we argue here that it is just a
consequence of Mottness, the splitting of the spectral weight over two
correlated bands. Subsequent experiments are discussed. While it is possible
to use this model to address a possible BCS/BEC crossover, we do not explore
this here. Such a crossover has been explored previously from a Luttinger
surfaceSetty (2021).

##  II Superconductivity in the HK model

Superconductivity in the cuprates necessitates a solution to at least the
Cooper instability in a doped Mott insulator. We have shownPhillips  _et al._
(2020) previously that this can be done exactly by solving Cooper instability
equation that arises from the HK analogue

$\displaystyle H=H_{\mathrm{HK}}-H_{p},\qquad
H_{p}=\frac{g}{V}\sum_{{\mathbf{k}},{\mathbf{k}}^{\prime}}b_{{\mathbf{k}}}^{\dagger}b_{{\mathbf{k}}^{\prime}}$
(3)

of the pairing Hamiltonian for a doped Mott insulator. Here
$b_{{\mathbf{k}}}=c_{-{\mathbf{k}}\downarrow}c_{{\mathbf{k}}\uparrow}$ is the
$s$-wave pair creation operator at zero total momentum. As is well known for
Mott systems, the single-particle Green function for $H_{\rm HK}$,

$\displaystyle G_{{\mathbf{k}}\sigma}(i\omega_{n})$
$\displaystyle\equiv-\int_{0}^{\beta}d\tau\;\langle{c_{{\mathbf{k}}\sigma}(\tau)c_{{\mathbf{k}}\sigma}^{\dagger}(0)}\rangle
e^{i\omega_{n}\tau}$ (4) $\displaystyle G_{{\mathbf{k}}\sigma}(i\omega_{n}\to
z)$
$\displaystyle=\frac{1-\langle{n_{{\mathbf{k}}\overline{\sigma}}}\rangle}{z-\xi_{{\mathbf{k}}}}+\frac{\langle{n_{{\mathbf{k}}\overline{\sigma}}}\rangle}{z-(\xi_{{\mathbf{k}}}+U)}.$
(5)

exhibits a bifurcation of the spectral weight between lower (l) and upper (u)
bands with weights that are determined by the electron filling, $1-\langle
n_{{\mathbf{k}}\overline{\sigma}}\rangle$ and $\langle
n_{{\mathbf{k}}\overline{\sigma}}\rangle$, respectively. It is this
bifurcation that leads to zeros of the real part of the Green function. The
features in the lower band are created with the momentum-projected operators
$\xi_{k\sigma}=c_{k\sigma}(1-n_{k\overline{\sigma}})$. In the Hubbard model,
the corresponding operators for the lower band are not known exactly. With
these operators, it is easy to see that the physics of the lower Hubbard band
in the HK model is not that of a Fermi liquid. Namely, there are excited
states of the HK model that have no correspondence with those in a FL.
Consider a two-particle excitation. This would be generated by applying the
$\xi^{\dagger}_{k\sigma}$ operator twice. That is, the excited states should
be described by $\xi^{\dagger}_{k\uparrow}\xi^{\dagger}_{k\downarrow}$.
However, this operator is explicitly zero. That is, there are excitations in
the HK model that have no counterpart by acting with the single-particle
operators. At work here is the fact that as long as $U\neq 0$, a Luttinger
surface of zeros is present. There are strictly no pure pole-like excitations.
Hence, there is no FL part of the HK model as long as $U\neq 0$.

We have also shownPhillips  _et al._ (2020) that not only is the Cooper
instability exactly solvable but so is the exact pair susceptibility. The
exact susceptibilityPhillips  _et al._ (2020)

$\displaystyle\chi(i\nu_{n})\equiv\frac{1}{V}\sum_{{\mathbf{k}},{\mathbf{k}}^{\prime}}\int_{0}^{\beta}d\tau\;e^{i\nu_{n}\tau}\langle{Tb_{{\mathbf{k}}}(\tau)b_{{\mathbf{k}}^{\prime}}^{\dagger}}\rangle_{g}$
(6)

can be expressed, in the normal state, in terms of the bare susceptibility

$\displaystyle\chi(i\nu_{n})=\frac{\chi_{0}(i\nu_{n})}{1-g\chi_{0}(i\nu_{n})}$
(7)

which is given by

$\displaystyle\chi_{0}(i\nu_{n})$
$\displaystyle=\chi_{0}^{ll}+\chi_{0}^{uu}+\chi_{0}^{lu}+\chi_{0}^{ul}$ (8)
$\displaystyle\chi_{0}^{ab}$
$\displaystyle=\frac{1}{V}\sum_{{\mathbf{k}}}n_{{\mathbf{k}}\uparrow}^{a}n_{-{\mathbf{k}}\downarrow}^{b}\frac{f(\omega_{{\mathbf{k}}}^{a})+f(\omega_{-{\mathbf{k}}}^{b})-1}{i\nu_{n}-\omega_{{\mathbf{k}}}^{a}-\omega_{-{\mathbf{k}}}^{b}}$
(9)

where for $\omega_{{\mathbf{k}}}^{l}=\xi_{{\mathbf{k}}}$ and
$\omega_{{\mathbf{k}}}^{u}=\xi_{{\mathbf{k}}}+U$,
$n_{{\mathbf{k}}\sigma}^{u}=\langle{n_{{\mathbf{k}}\overline{\sigma}}}\rangle_{0}$
and $n_{{\mathbf{k}}\sigma}^{l}=1-n_{{\mathbf{k}}\sigma}^{u}$, and $f(\omega)$
the Fermi function at temperature $T$, the superscripts $ab$ may represent
$ll$, $uu$, $lu$, or $ul$.

A consequence of the expression of the pair susceptibility is that the
divergence at $\chi_{0}=1/g$ is expected to be a second-order transition to
the bulk superconducting state. As we will see, this divergence is not
coincident with the turn-on of the gap. To calculate the susceptiblity, we
will work with the exact finite temperature occupancy,

$\displaystyle\langle{n_{{\mathbf{k}}\sigma}}\rangle=\frac{1}{2}\langle{n_{{\mathbf{k}}}}\rangle=\frac{e^{-\beta\xi_{{\mathbf{k}}}}+e^{-\beta(2\xi_{{\mathbf{k}}}+U)}}{1+2e^{-\beta\xi_{{\mathbf{k}}}}+e^{-\beta(2\xi_{{\mathbf{k}}}+U)}},$
(10)

so as to give the correct temperature dependence of $\chi(T)$ explicitly.

Figure 1: The temperature dependence of $\chi_{0}$(left) and $\chi$(right) for
a 3D-HK model with superconducting pairing. The energy dispersion was
approximated by a parabola dispersion, and the chemical potential was set to
half filling of the lower Hubbard band, $\mu=0.5W$ where $W$ is the bandwidth
of the lower Hubbard band. We take $W=1$ in the following calculation. The
temperature $T_{c}$ at which $\chi=\chi_{0}/(1-g\chi_{0})$ diverges
($\chi_{0}=1/g$) represents the true superconducting transition onset.

In Fig. 1, we plot the zero-frequency bare susceptibility $\chi_{0}(T)$ as
well as $\chi(T)$. The divergence of $\chi(T=T_{c})$ would imply a diverging
length scale and thus a possible second-order phase transition temperature
$T_{c}$. Above $T_{c}$ there is no extra sigularity in the susceptibility,
which is the same case as in the BCS theory. To reiterate, the onset of the
gap, a mean-field notion, is distinct from the divergence of the
susceptibility.

###  II.1 Mean-field Theory of HK

Given that the susceptibility calculation is exact, we can assert that the
ground state of the HK model with a pairing term is superconducting. To
describe this state, we resort to a mean-field description in the spirit of
BCS. In our previous work, we implemented the mean-field by an appropriately
chosen pairing wave function. To go beyond such a ground-state treatment, we
explicitly diagonalize $H$, obtaining all eigenstates and measuring
observables in grand canonical ensemble, including their full temperature
dependence. The procedure we use has been outlined in the Supplementary
Materials of Ref. Phillips  _et al._ , 2020 and by Zhu, et al.Zhu  _et al._
(2021). In constructing the mean-field, we take advantage that the HK
Hamiltonian does not mix the various k-states. In terms of the pair amplitude,
$\Delta\equiv(g/V)\sum_{{\mathbf{k}}}b_{{\mathbf{k}}}$, we formulate the mean-
field,

$\begin{split}H_{p}&=\frac{g}{V}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}}b_{{\mathbf{k}}}^{\dagger}b_{{\mathbf{k}}^{\prime}}\\\
&=\frac{V}{g}\Delta^{\dagger}\Delta\\\
&=\frac{V}{g}(\delta\Delta^{\dagger}+\overline{\Delta}^{*})(\delta\Delta+\overline{\Delta})\\\
&\approx\sum_{{\mathbf{k}}}(\overline{\Delta}b_{{\mathbf{k}}}^{\dagger}+\overline{\Delta}^{*}b_{{\mathbf{k}}})-\frac{V}{g}|\overline{\Delta}|^{2},\end{split}$
(11)

entirely on the pairing term where we have introduced the average of the
pairing amplitude to be $\overline{\Delta}=\langle\Delta\rangle$. In the last
step, mean-field amounts to ignoring the second-order fluctuation
$O(\delta\Delta^{2})$ term. The mixing between different momentum sectors
averages out because of the momentum-diagonal structure of HK. Thus the mean-
field(MF) HK Hamiltonian can be block diagonalizedPhillips  _et al._ (2020);
Zhu  _et al._ (2021)

$\begin{split}H&\approx\sum_{{\mathbf{k}}\in\rm{HFBZ}}H^{MF}_{{\mathbf{k}}}+\frac{V}{g}|\overline{\Delta}|^{2}\\\
H^{MF}_{{\mathbf{k}}}&=\sum_{s=\pm}\xi_{{\mathbf{k}}}\sum_{\sigma}c_{s{\mathbf{k}}\sigma}^{\dagger}c_{s{\mathbf{k}}\sigma}+Uc_{s{\mathbf{k}}\uparrow}^{\dagger}c_{s{\mathbf{k}}\uparrow}c_{s{\mathbf{k}}\downarrow}^{\dagger}c_{s{\mathbf{k}}\downarrow}\\\
&\qquad-\left(\overline{\Delta}^{*}c_{-s{\mathbf{k}}\downarrow}c_{s{\mathbf{k}}\uparrow}+\overline{\Delta}c_{s{\mathbf{k}}\uparrow}^{\dagger}c_{-s{\mathbf{k}}\downarrow}^{\dagger}\right).\end{split}$
(12)

The summation over k is carried out inside half of the first Brillouin
zone(HFBZ), while the boundary terms with momenta restricted to the edge of
the first Brillouin zone are neglected as they are suppressed by a factor of
$1/N$ and hence vanish in the thermodynamic limit.

The decomposed Hamiltonian lives in a Fock space $F_{{\mathbf{k}}}$ which
contains 4 fermion currents,
$\ket{n_{{\mathbf{k}}\uparrow},n_{{\mathbf{k}}\downarrow},n_{-{\mathbf{k}}\uparrow},n_{-{\mathbf{k}}\downarrow}}$
and spans a 16-dimensional spaceZhu  _et al._ (2021). Due to the fermion
parity conservation of $H_{{\mathbf{k}}}^{MF}$, this Fock space can be
decomposed by parity into $F_{{\mathbf{k}}}=F_{{\mathbf{k}}}^{odd}\oplus
F_{{\mathbf{k}}}^{even}$. The even sector is further block diagonalized into 3
subspaces, $F_{{\mathbf{k}}}^{even}=F_{{\mathbf{k}}}^{PB}\oplus
F_{{\mathbf{k}}}^{ST}\oplus F_{{\mathbf{k}}}^{Mix}$, as shown in Table 1 :

subspace | eigenvalue | degeneracy | basis
---|---|---|---
PB | $E_{{\mathbf{k}}}^{4}\equiv 2\xi_{{\mathbf{k}}}+U$ | 2 | $\ket{1100},\ket{0011}$
ST | $E_{{\mathbf{k}}}^{5}\equiv 2\xi_{{\mathbf{k}}}$ | 3 | $\begin{array}[]{c}\ket{1010},\ket{0101},\\\ (\ket{1001}-\ket{0110})/\sqrt{2}\\\ \end{array}$
Mixing | $\begin{array}[]{c}E_{{\mathbf{k}}}^{1}\equiv E_{1}\\\ E_{{\mathbf{k}}}^{2}\equiv E_{2}\\\ E_{{\mathbf{k}}}^{3}\equiv E_{3}\\\ \end{array}$ | $\begin{array}[]{c}1\\\ 1\\\ 1\end{array}$ | $\begin{array}[]{c}\ket{0000},\\\ (\ket{1001}+\ket{0110})/\sqrt{2},\\\ \ket{1111}\end{array}$
Odd$\times 4$ | $\begin{array}[]{c}E_{{\mathbf{k}}}^{6}\equiv E_{-}\\\ E_{{\mathbf{k}}}^{7}\equiv E_{+}\end{array}$ | $\begin{array}[]{c}4\\\ 4\end{array}$ | $\begin{array}[]{c}\\{\ket{1000},\ket{1110}\\}\\\ \\{\ket{0100},\ket{1101}\\}\\\ \\{\ket{0010},\ket{1011}\\}\\\ \\{\ket{0001},\ket{1011}\\}\\\ \end{array}$
Table 1: Block decomposition and the energy levels of the HKSC mean-field
hamiltonian

  1. 1.

The states in the Pauli Blocking (PB) states, possess an energy level at
$E_{{\mathbf{k}}}^{4}=2\xi_{{\mathbf{k}}}+U$ with degeneracy 2, which have
definite total electron number 2(not participating the superconducting
pairing) and is blocked by the repulsion $U$.

  2. 2.

The states in the Spin Triplet (ST) states, possess an energy level at
$E_{{\mathbf{k}}}^{5}=2\xi_{{\mathbf{k}}}$ with degeneracy 3, which also have
two electrons(not participating the superconducting pairing).

  3. 3.

The 3-dimensional particle number mixing (Mixing) states related by the off-
diagonal superconducting pairing.

The Hamiltonian matrix in the 3-dimensional mixing states subspace is

$\displaystyle\left(\begin{array}[]{ccc}0&-\sqrt{2\overline{\Delta}}&0\\\
-\sqrt{2\overline{\Delta}}&2\xi_{{\mathbf{k}}}&-\sqrt{2\overline{\Delta}}\\\
0&-\sqrt{2\overline{\Delta}}&4\xi_{{\mathbf{k}}}+2U\end{array}\right)$ (16)

with corresponding 3 energy levels ($i=1,2,3$)

$\displaystyle
E_{i}=2\xi_{\mathbf{k}}+\frac{2U}{3}+\frac{4}{\sqrt{3}}E_{\mathbf{k}}^{\text{even
}}\cos\left(\theta_{\mathbf{k}}+\frac{2\pi}{3}i\right)$ (17) $\displaystyle
E_{\mathbf{k}}^{\mathrm{even}}=\sqrt{\left(\xi_{\mathbf{k}}+\frac{U}{2}\right)^{2}+\overline{\Delta}^{2}+\frac{U^{2}}{12}}$
(18)
$\displaystyle\theta_{\mathbf{k}}=\frac{1}{3}\arccos\left[\frac{q_{\mathbf{k}}}{\left(\sqrt{3}E_{\mathbf{k}}^{\text{even
}}\right)^{3}}\right]$ (19) $\displaystyle
q_{\mathbf{k}}=U\left(U^{2}+\frac{9}{2}U\xi_{\mathbf{k}}+\frac{9}{2}\xi_{\mathbf{k}}^{2}-\frac{9}{4}\overline{\Delta}^{2}\right).$
(20)

Similarly, the odd sector can be written into a direct-sum of 4 equivalent
subspace, e.g. {$\ket{1000},\ket{1110}$}. Here the Hamiltonian matrix in each
2-dimensional subspace is

$\displaystyle\left(\begin{array}[]{cc}\xi_{{\mathbf{k}}}&-\overline{\Delta}\\\
-\overline{\Delta}&3\xi_{{\mathbf{k}}}+U\end{array}\right)$ (23)

and the odd sectors share the 2 energy levels,

$\displaystyle E_{\pm}=2\xi_{\mathbf{k}}+\frac{U}{2}\pm
E_{\mathbf{k}}^{\text{odd }}$ (24) $\displaystyle E_{\mathbf{k}}^{\text{odd
}}=\sqrt{\left(\xi_{\mathbf{k}}+\frac{U}{2}\right)^{2}+\overline{\Delta}^{2}}.$
(25)

The energy levels and bases in each subspace are shown in Table 1. Thus we
have all energy levels of the system. The ground state is recognized as the
eigenvector coresponding to $E_{1}$, which is a linear combination of the 3
occupancy states: $\ket{0000}$ from $\Omega_{0}\otimes\Omega_{0}$,
$(\ket{1001}+\ket{0110})/\sqrt{2}$ from $\Omega_{1}\otimes\Omega_{1}$ and
$\ket{1111}$ from $\Omega_{2}\otimes\Omega_{2}$. Our previous variational
treatmentPhillips  _et al._ (2020) concurs with this result.

The work-horse of the mean-field treatment is the self-consistent equation for
the pair amplitude:

$\displaystyle\overline{\Delta}=\frac{g}{V}\sum_{{\mathbf{k}}}\langle{c_{{\mathbf{k}}\uparrow}c_{{\mathbf{-k}}\downarrow}}\rangle.$
(26)

The right-hand side has a dependence on $\overline{\Delta}$ through the
thermal average over the energy levels enumerated in Table 1. We can also
treat the self-consistent equation for $\overline{\Delta}$ as the solution to
the extremum of the free energy,

$\displaystyle F=-k_{B}T\ln Z,$ (27)

where

$\displaystyle Z=\sum_{i,k\in HFBZ}e^{-\beta E_{{\mathbf{k}}}^{i}}.$ (28)

Figure 2: (a) The variation of the free Energy $F$ as a function of the the
pair amplitude $\Delta$ for an HK superconductor at $U/W=1$ and $g=0.3$. The
chemical potential is set at half filling of the lower band ($\mu=0.5$). For
$T>T_{p}$, the global minimum of free energy is $\Delta=0$. At $T=T_{p}$, the
free energy at $\Delta=0$ and $\Delta=\Delta_{p}$ coincides, thus indicating a
first-order transition. When $T_{2}<T<T_{p}$, $\Delta=0$ is still a local
minimum, while the global minimum obtains at finite $\Delta$. For $T=T_{2}$,
the only minimum occurs at finite $\Delta$, and $\Delta=0$ becomes a local
maximum(the inflection point). (b) The dashed purple line and the solid blue
line are the solution to the pair amplitude from the self consistent equation
at different temperatures(the $\Delta-T$ plot) for BCS($U/W=0$) and
HK($U/W=1$). The Green dot shall represent the first-order phase transition
and the red dot represents the second-order phase transition. In both cases,
the axes are scaled according to respective $T_{2}$, and respective zero
temperature pair amplitude $\Delta_{0}$.

Both procedures yield identical results and are tabulated in Fig. 2. Displayed
first is the free energy as $\overline{\Delta}$ is varied. Of first note is
that for $T>T_{p}$, the global minimum of free energy remains zero. At
$T=T_{p}$ ($p$ for pair), there are two degenerate minima, with
$\overline{\Delta}$=0 and $\overline{\Delta}=\Delta_{p}\neq 0$. The technical
definition of $T_{p}$ then is $\partial
F/\partial\Delta|_{\overline{\Delta}=\Delta_{p}}=0$ and $F(\Delta_{p})=F(0)$.
This degeneracy is lifted by lowering the temperature such that
$T_{2}<T<T_{p}$. The solution to
$\partial^{2}F/\partial\Delta^{2}|_{\Delta=0}=0$ defines $T_{2}$ as the
inflection point. As we will see, $T_{2}$ will correspond to the divergence of
the susceptibility.

To corroborate this picture, we solve the self-consistent equations for the
the gap. The dashed curve in Fig. 2(b) corresponds to BCS theory which has a
unique turn-on temperature for the pair amplitude at the red dot, namely
$T=T_{2}$. In HK, however from our evaluation of the free energy, we see
multiple solutions for the gap turn-on distinct from the feature at $T_{2}$.
As shown, the $\Delta-T$ curve for the mean-field HK exhibits a nontrivial
back-folding behavior above $T_{2}$. This multivaluedness confirms that there
are multiple choices of $\Delta$ that make the first derivative of the free
energy with respect to the gap vanish. Only one of these gaps yields the
global minimum of the free energy. The equivalence of the free energy at
$\Delta=0$ and $\Delta_{p}$ implies a first order phase transition at this
temperature $T_{p}$(green dot), which ensures the accuracy of mean-field
treatment since the first order transition does not have any diverging
fluctuations. The degenerate minimum of the free energy at $T=T_{p}$ are shown
in Fig. 2(a) with the green solid curve. There are two mechanisms to lift the
degeneracy: 1) decrease the temperature as shown in Fig. (2a) or 2) decrease
$U$ from the value shown in Fig. (2b). What is not shown is that there is a
critical value of $U$ that is needed to destroy the back-folding in Fig. (2b).
This critical value diminishes as the pairing strength, $g$, decreases. While
it is suggestive at this point that it is Mott physics that leads to the
first-order nature of the superconducting transition, we will confirm this by
a detailed evaluation of the Landau expansion parameters.

Figure 3: Pair amplitude of s-wave and d-wave pairing at the same pairing
strength $g=0.3$ and repulsion strength $U/W=1$.

An analysis of the HK model under d-wave pairing could be carried out
similarly by replacing the pairing term $H_{p}$ in Eq.(3) with a d-wave form

$\displaystyle
H_{d}=\frac{g_{d}}{V}\sum_{{\mathbf{k}},{\mathbf{k}}^{\prime}}\tilde{b}_{{\mathbf{k}}}^{\dagger}\tilde{b}_{{\mathbf{k}}^{\prime}},$
(29)

where
$\tilde{b}_{{\mathbf{k}}}=c_{-{\mathbf{k}}\downarrow}c_{{\mathbf{k}}\uparrow}(\cos
k_{x}-\cos k_{y})$ is the $d$-wave pair creation operator in 2-dimensions.
After performing the same mean-field calculation on
$\tilde{b}_{{\mathbf{k}}}$, we find that the pair amplitude, as shown in Fig.
3, also exhibits a nontrivial back-folding behavior above $T_{2}$, which is
recognized as the key major difference between a BCS superconductor and an HK
superconductor at the mean field level. As there is little difference between
these two cases, in the subsequent analysis only the s-wave case will be
analyzed.

One remaining subtlety is that the calculations in this section are all based
on mean-field theory, which is exact only if the transition is truly first
order. Strictly, they break down at the critical point of a second order
transition. Since the exact pair susceptibility diverges at a temperature
distinct from $T_{p}$, we need to consider two possibilities:

  1. 1.

The phase transition is first-order at $T_{p}$, and the calculation of the
susceptibility cannot be applied to the theory below the transition
temperature at which $T_{p}$($\Delta$ jumps to a finite value instead of
continuously growing from zero when $T=T_{p}$). The mean-field calculation is
exact in this scenario.

  2. 2.

The phase transition is second-order at $T_{c}$. The global minimum of the
free energy at $T_{p}$ shall be excluded due to the divergence of fluctuations
near the second-order phase transition, which destroys the mean-field theory.
Only the non-MF calculation (e.g. susceptibility) is right.

In the first case, the calculation of $T_{c}$ is internally inconsistent as it
presumed that $\Delta_{p}=0$ for $T>T_{c}$. In the second case, we take
advantage of the exact diagonalizability and use the Ginzburg criterion to
estimate the temperature range where the fluctuations diverge (mean-field
theory breaks down).

###  II.2 Ginzburg Criterion

The unusual first-order nature of the phase transition in an HK superconductor
is manifest at the mean-field level. The crucial question is then: Is this
first-order Transition just an artifact? In other words, how valid is the
mean-field theory, or equivalently, what is its range of validity? To address
this, we compute the Ginzburg criterion, presuming that the phase transition
is second-order at $T=T_{c}$.

Recall that the Landau-Ginzburg description of a traditional superconductor is
equivalent to the mean-field treatment of the BCS theory around the transition
pointColeman (2015). We demostrate the same result for HK superconductors. In
an s-wave superconductor which is the case considered here, the order
parameter is given by

$\Psi(x)\propto\langle{c_{\uparrow}(x)c_{\downarrow}(x)}\rangle.$ (30)

In the absence of a magnetic field, the Landau free energy is

$F[\Psi]=\int
d^{d}x\left(a|\Psi(x)|^{2}+\frac{1}{2}b|\Psi(x)|^{4}+K|\nabla\Psi(x)|^{2}\right),$
(31)

where $a$ is directly related to the susceptibility and $K$ is the rigidity.
Close to a second-order critical point $T_{c}$, that is $|\Psi|\ll 1$, we may
neglect the quartic terms. As a result, the free energy,

$\begin{split}F[\Psi,\Psi^{*}]&=\int
d^{d}x\left(a|\Psi(x)|^{2}+K|\nabla\Psi(x)|^{2}\right)\\\
&=\sum_{{\mathbf{k}}}(a+Kk^{2})|\Psi(k)|^{2}\\\
\Psi(x)&=\frac{1}{\sqrt{V}}\sum_{{\mathbf{k}}}\Psi(k)e^{ikx},\end{split}$ (32)

by performing the Fourier transform of the order parameter. The Helmholtz free
energy $A(T)$ is given by

$\begin{split}A(T)&=-k_{B}T\ln Z\\\ &=-k_{B}T\ln\int
D[\Psi,\Psi^{*}]e^{-F[\Psi,\Psi^{*}]/k_{B}T}\\\
&=k_{B}T\sum_{{\mathbf{k}}}\ln\left(\frac{k_{B}T}{a+Kk^{2}}\right).\end{split}$
(33)

We write $a(T)=\alpha t$ with $t=\frac{T-T_{c}}{T_{c}}$. The singular
contribution to the specific heat $C_{V}=-TA^{\prime\prime}(T)$ comes when
differentiating with respect to $T$. The dimensionless heat capacity per unit
cell is

$\begin{split}c&=\frac{C_{V}}{N_{s}k_{B}}\\\
&=\frac{\alpha^{2}\mathsf{a}^{d}}{K^{2}}\int^{\Lambda}\frac{d^{d}k}{(2\pi)^{d}}\frac{1}{(\xi^{-2}+k^{2})^{2}}\\\
&=\frac{\alpha^{2}\mathsf{a}^{d}}{K^{2}}\xi^{4-d}\int^{\Lambda\xi}\frac{d^{d}q}{(2\pi)^{d}}\frac{1}{(1+q^{2})^{2}},\end{split}$
(34)

where $\Lambda\sim\mathsf{a}^{-1}$ is an ultravilot cut-off and $\mathsf{a}$
is the lattice constant(set to unit in the following calculation),
$\xi=(K/a)^{1/2}=(K/\alpha)^{1/2}|t|^{-1/2}$. The fluctuation contribution is
small when

$\frac{\alpha^{2}\mathsf{a}^{d}}{K^{2}}\xi^{4-d}\ll 1,$ (35)

which entails $|t|\gg t_{G}$ where the latter is defined,

$t_{G}=\left(\mathsf{a}^{2}\frac{\alpha}{K}\right)^{\frac{d}{4-d}},$ (36)

as the Ginzburg reduced temperature. This calculation then just hinges on
knowing $\alpha$ and $K$ and we can estimate the temperature region where
mean-field theory is not valid.

To obtain $\alpha$, we calculate the spatial average of the order parameter

$\begin{split}\overline{\Psi}&\equiv\frac{g}{V}\int
d^{d}x\langle{c_{\uparrow}(x)c_{\downarrow}(x)}\rangle\\\
&=\frac{g}{V}\sum_{{\mathbf{k}}}\langle{c_{{\mathbf{k}}\uparrow}c_{{\mathbf{-k}}\downarrow}}\rangle\\\
&=\overline{\Delta}.\end{split}$ (37)

Thus, the Landau free energy parameters can be calculated by differentiating
the free energy density $f$ with respect to $\overline{\Delta}^{2}$

$a=\left.\frac{\partial
f}{\partial\overline{\Delta}^{2}}\right|_{\overline{\Delta}=0},$ (38)

and the parameter $\alpha$ can be calculated by linear fitting the $a(T)$
around the transition temperature as shown in Fig. 4. The key point is that
$\alpha$ does not vary appreciably for the system parameters $U/W$ and $g$.

Figure 4: The dependence of the Landau expansion parameter $\alpha=a/t$ on the
HK interaction strength $U$ at different values of the superconducting pairing
strength $g$. We claim that at $U=0$, $\alpha\approx 2.1$ and at $U>\mu=0.5$,
$\alpha\approx 1.3$ which are around the same scale. The magnitude of $\alpha$
does not vary appreciably for all values of $g$.

We also have the superfluid stiffness defined by

$T[\theta]=\int d^{d}x\frac{1}{2}\rho_{s}(\nabla\theta)^{2},$ (39)

where $\theta$ is the phase of the superfluid
$\Psi(x)=\overline{\Psi}e^{i\theta}$. Thus the relation between $K$ and
$\rho_{s}$ is

$\displaystyle\frac{1}{2}\rho_{s}(\nabla\theta)^{2}=K|\nabla\Psi|^{2}=K\overline{\Psi}^{2}(\nabla\theta)^{2}$
(40) $\displaystyle
K=\frac{\rho_{s}}{2\overline{\Psi}^{2}}=\frac{\rho_{s}}{2\overline{\Delta}^{2}}.$
(41)

From Fig. 4 and Fig. 5, we can read that at $\mu=0.5$, $U>\mu$, the value of
$\alpha$ and $\rho_{s}$ is $\alpha\approx 1.3$,
$K\approx\rho_{s}/2\Delta^{2}\approx 10^{4}$. From Eq. (36), we may estimate
$t_{G}\approx 10^{-11}$. Consequently, the Ginzburg temperature is
sufficiently small to guarantee the validity of the mean-field calculation,
which predicts a first order transition. Hence, the presumption that the phase
transition is second-order shall be ruled out when $T_{p}>T_{2}$. We will
discuss the consequences of this first order transition in more depth.

Figure 5: The rigidity at transition temperature($T\approx T_{c}$, solid line)
and at zero temperature($T\approx 0$, dotted points). The calculation(see
appendix) was performed with $\mu=0.5$, $g=0.2$. The value of $K$ was
calculated with $\overline{\Delta}=10^{-5}$ and $T=T_{c}$(where the
susceptibility diverges).

###  II.3 Phase Diagram

To put this all together, we focus on the phase diagram in the $T-U$ plane
presented in Fig. 6 with a superconducting paring strength $g=0.3$. We observe
that for a small value of $U/W<0.13$, the superconducting transition is
second-order (red line) in which the divergence of the pair susceptibility and
the turn-on of the gap are coincident. In general, the critical value of
$U/W$, denoted by $U_{t}/W$, for the transition to be second order decreases
as $g$ decreases. Hence, it is only the weakly interacting regime in which we
find the traditional result that the superonducting transition is second
order. As pointed out in our previous paper, in no regime (except $U=0$) does
a BCS picture apply. For example, even in the weakly interacting regime, the
transition temperature compared with BCS($U=0$) value increases as $U$
increases implying that the multiband nature of the HK model is at play in
driving the enhancement of superconductivity. The second-order transition line
ends at the tricritical point located at $U_{t}/W=0.13$ and $T_{t}/W=0.0255$.
It is at this point that the local and global minima merge. For larger values
of $U>U_{t}$, the transition becomes a first-order one(green dots). The phase
transition temperature $T_{p}$ saturates when $U/W>0.5$ which coincides the
elimination of the single-occupancy region $\Omega_{1}$ and double occupancy
region $\Omega_{2}$ boundary in the HK model.

We can also compare $T_{2}$ with $T_{c}$. As is evident, they both coincide as
illustrated in Fig. 6. This is significant because $T_{c}$ is computed from a
divergence of the susceptibility while $T_{2}$ follows from a solution to a
mean-field equation.

Figure 6: Phase diagram in the $T-U$ plane of the HK superconductor for
coupling $g=0.3$.The green line corresponds to the first-order transition
temperature $T_{p}$ at which the global minimum of the free energy switches
from $\Delta=0$ to a finite value. The red line corresponds to $T_{2}$ at
which $\Delta=0$ changes from a local minimum into a local maximum. The gray
line is the pair susceptibility diverging temperature $T_{c}$. The gray and
red line coincide, thereby corroborating that the pair-susceptibility
divergence is coincident with the inflection point in the free energy.

For an HK superconductor with first order transition, the temperature $T_{2}$
represents the elimination of the metastable state(local minimum).

In all the computations thus far, either the chemical potential or the filling
has remained fixed. To make contact with the cuprates, a phase diagram of
$T_{c}$ versus filling or chemical potential is needed. This can be done using
the same machinery. Shown in Fig. 7 is a plot of $T_{c}$ versus filling for
two distinct cases: 1) overlapping lower and upper Hubbard bands and 2) no
overlap. This can be engineered simply by changing the value of $U$. In
general electron and hole doping yield the mirror results. When the bands
overlap, a metallic state always ensues, though a non-Fermi liquid onePhillips
_et al._ (2020) and $T_{c}$ is minimized at half-filling though it does not
vanish. A vanishing of $T_{c}$ at the Mott insulating state obtains only for
$U>W$ as is shown in the plots in which $U=1$ and $U=2$. The familiar dome-
shaped phase diagram obtains as dictated by the dome-shaped superfluid
stiffness reported previouslyPhillips  _et al._ (2020). The dome shape here is
a direct consequence of Mott physics.

Figure 7: Phase diagram in the $T_{c}-\langle{n_{\sigma}}\rangle$ plane of the
HK superconductor for coupling $g=0.3$. The blue line ($U=0.5$) represents the
overlapping bands or metallic phase in which $U<W$. The pink and purple
line($U=1,2$) represent the non-overlapping bands (Mott insulator ground state
at half filling, $\langle{n_{\sigma}}\rangle=0.5$).

###  II.4 Landau Expansion Parameters

To complete a theory of a first order transiton, we rewrite the Landau free
energy

$\displaystyle
F[\overline{\Delta}]=a\overline{\Delta}^{2}+b\overline{\Delta}^{4}+c\overline{\Delta}^{6}+O(\overline{\Delta}^{6}).$
(42)

in terms of the 4th order expansion parameter $b$. For a first-order
transition, $b<0$ and the 6th order term be positive to keep the free energy
bounded from below. When the first-order transition obtains, both the value of
the free energy and the first derivative with respect to $\delta$ vanish,

$\displaystyle F[\overline{\Delta}_{p}]$
$\displaystyle=a\overline{\Delta}_{p}^{2}+b\overline{\Delta}_{p}^{4}+c\overline{\Delta}_{p}^{6}=0$
(43) $\displaystyle F^{\prime}[\overline{\Delta}_{p}]$
$\displaystyle=2a\overline{\Delta}_{p}+4b\overline{\Delta}_{p}^{3}+6c\overline{\Delta}_{p}^{5}=0.$
(44)

The solution is $\overline{\Delta}_{p}^{2}=-2a/b$ and $b^{2}=4ac$, thus
proving that $a>0,b=-\sqrt{4ac}<0$ is required for this kind of first order
transition to exist. Recall that for an HK superconductor, extreme
factorizability in momentum sections allows us to write the free energy as

$\displaystyle F[\overline{\Delta}]=\sum_{k\in
HFBZ}F_{{\mathbf{k}}}[\overline{\Delta}].$ (45)

As a result, the free-energy expansion parameters can be decomposed as
$a=\sum_{{\mathbf{k}}}a_{{\mathbf{k}}}$ and
$b=\sum_{{\mathbf{k}}}b_{{\mathbf{k}}}$. To organize our results, we recall
that what makes the HK model a non-Fermi liquid is that at zero-temperature,
as a result of $U$, singly occupied states, $\Omega_{1}$, exist below the
chemical potential which never obtains in a Fermi liquid. This is depicted in
Fig. 8(a). The results for the momentum-resolved Landau expansion coefficients
are shown in panels Fig. 8(b,c).

Figure 8: (a) The occupancy of the HK(U=0.3) and Fermi liquid(FL) at finite
temperature $T/W=0.01$. (b) The 2nd expansion parameter $a_{{\mathbf{k}}}$ of
the Landau Free Energy over the energy levels $\xi_{{\mathbf{k}}}$. (c) The
4th expansion parameter $b_{{\mathbf{k}}}$ of the Landau Free Energy over the
energy levels $\xi_{{\mathbf{k}}}$.

In the BCS case, the negative contribution of $a_{{\mathbf{k}}}$ comes from
the states around the Fermi surface, where the occupancy changes from 2 to
zero. This arises from the sharpness of the Fermi surface. In the HK
superconductor, similar sign changes at the occupancy boundaries from
$\Omega_{2}$ to $\Omega_{1}$ and from $\Omega_{1}$ to $\Omega_{0}$.
$a_{{\mathbf{k}}}$ is always possitive except around the boundary of different
occupancy region. When the positive contribution exactly cancels, the negative
contribution, $a$ vanishes, implying that the free energy changes from a
concave function into a convex one around $\Delta=0$. Most crucial here,
however, is the behaviour of $b_{{\mathbf{k}}}$. For the HK superconductors,
while the distribution of $a_{{\mathbf{k}}}$ follows the BCS case, the value
of $b_{{\mathbf{k}}}$ differs drastically from BCS. In BCS, $b_{{\mathbf{k}}}$
is always positive; thus $b>0$ is true for any temperature. In HK, however,
inside the single occupancy region $\Omega_{1}$, which is not present in any
fermi liquid, $b_{{\mathbf{k}}}$ becomes negative! Together with a surpression
of the positive $b_{{\mathbf{k}}}$ contribution, it is possible that $b<0$ and
a first order transition emerges. Consequently, it is Mottness that drives the
first-order nature of the transition in the HK model. Perhaps this is true in
general.

###  II.5 Condensation Energy

The deviation of an HK superconductor from that of the BCS type is also
manifest from the condensation energy, defined by the energy difference
between the superconducting and normal states,

$\displaystyle E_{\rm cond}=E_{SC}-E_{N}.$ (46)

In a traditional s-wave BCS superconductor, the condensation energy is well-
known to be propotional to the square of the pair amplitude $E_{\rm
cond}^{BCS}=-N(0)\Delta^{2}/2$, where $N(0)$ is the density of states around
the chemical potential. In the HK superconductor, however, this relation holds
up to a modification to the coefficient

$\displaystyle E_{\rm cond}^{HK}=CN(0)\Delta^{2}/2,$ (47)

where $C$ is a pure number. The dependence of $C$ on pairing strength $g$ and
$U$ is plotted in Fig. (9).

Figure 9: The dependence of the condensation energy $E_{\rm cond}$ devided by
$N(0)\Delta^{2}$ at zero temperature in the $g-U$ plane. The chemical
potential is set $\mu/W=0.5$. The bottom red bold line represents the BCS
result $C=-1$. As is evident, HK superconductors have condensation energies
that generically deviate from the BCS result.

Generally, the coefficient $C$ have a smaller value in HK superconductors than
in BCS, demonstrates a stronger tendency to form a superconductor. This result
is surprising but consistent with what we know about the cuprates in which the
pairing temperatures are greatly enhanced despite the local repulsions. In
this case, enhancement arises because of the singly-occupied region. The
minimum value $C\approx-2$ can be achieved at a specific point $U$ smaller
than the chemical potential and quickly goes back to the BCS ratio $C=-1$ as
$U>\mu$, where the double occupancy region $\Omega_{2}$ goes away, and the
boundary between $\Omega_{1}$ and $\Omega_{2}$ is removed. This observation
again emphasizes the significance of the single occupancy region boundary in
the HK superconductor. It may imply a mechanism to enhance the
superconductivity by adding the occupancy boundary, where the Copper pair
formation is enhanced.

###  II.6 Heat capacity

Even when the superconducting gap develops continuously at the BCS transition,
the heat capacity $C(T)$ undergoes a jump discontinuity $\Delta C$ with the
universal ratio $\Delta C/C_{n}(T_{c})=12/7\zeta(3)\approx 1.43$, where
$C_{n}(T_{c})$ is the heat capacity in the normal phase at $T=T_{c}$. The
discontinuity obtains because the BCS Hamiltonian is an effective one that
varies with the temperature through the gap parameter $\Delta(T)$. The same is
true of HK even with a first-order transition as can be seen from the simple
argument. With the free energy given by
$-T\ln\operatorname*{tr}e^{-H/T}=\langle{H}\rangle-TS$, any system with a
temperature-independent Hamiltonian has a heat capacity coefficient
$C/T\equiv\partial_{T}S=\beta^{3}\langle{(H-\langle{H}\rangle)^{2}}\rangle$
that varies continuously with the spectrum of $H$. However, because the gap is
temperature dependent, the heat capacity

$C/T=\beta^{3}\left[\langle{(H-\langle{H}\rangle)^{2}}\rangle-\Delta\frac{\partial\Delta}{\partial{T}}\langle{H}\rangle\right],$
(48)

has an explicit temperature derivative. This term introduces a jump
discontinuity at $T=T_{c}$ since $\Delta(T)\sim\sqrt{T_{c}-T}$ (so
$\Delta\partial_{T}\Delta\sim-1$) on the superconducting side of the
transition whereas $\Delta(T)=0$ (so $\Delta\partial_{T}\Delta=0$) on the
normal side.

In the HKSC where the gap changes discontinuously, the heat capacity still
undergoes such a jump discontinuity, although the size of the discontinuity
$\Delta C$ no longer scales universally with $C_{n}(T_{c})$. Rather as we will
see, it depends on the Mott scale $U$ as well. When the gap changes
discontinuously with temperature, the entropy also changes discontinuously so
that the heat capacity $C=T\partial_{T}S$ is singular at the transition. We
show its behavior in Fig. 10, omitting the singularity. As is evident,
Mottness, as tracked by increasing $U$, enhances the heat capacity jump at
$T_{p}$, the temperature of the first-order transition. Similar enhancements
will be seen as well for the ultrasonic attenuation and the spin-lattice
relaxation rate.

Figure 10:  Heat capacity $C$ near the phase transition at $T=T_{p}$, at pair
coupling $g/W=0.53$ and for representative values of the Mott coupling $U$.
$C$ is normalized by its value $C_{n}(T_{p})$ on the normal side of the phase
transition.

###  II.7 NMR

The key feature of BCS is computability with the mean-field formalism. The
same can be done here as we know all the excited states and hence can
calculate all of the experimental quantities delineated by BCS. As an example,
we compute the spin-lattice relaxation rate which in a BCS superconductor
exhibits a peakZoli (1991) below $T_{c}$. The spins of atomic nuclei relax by
exchanging energy with their environment. As the spins in a superconductor are
in phase, there is an enhancement below $T_{c}$. The relaxation rate $1/T_{1}$
of nuclei in an electronic environment is related to the transverse dynamic
spin susceptibility of the quasiparticles, byColeman (2015)

$\frac{1}{TT_{1}}=\lim_{\omega\rightarrow
0}\frac{2k_{B}}{\gamma_{e}^{2}\hbar^{4}}\sum_{q}|A_{H}(q)|^{2}\frac{\Im{\chi^{s}}(q,\omega)}{\omega}$
(49)

where $\gamma_{e}$ is the electron gyromagnetic ratio, $A_{H}(q)$ is the
hyperfine coupling of the contact interaction with electron spins, and
$\Im{\chi^{s}}$ is the imaginary part of the spin susceptibility.

In a normal metal state, $\frac{\Im{\chi^{s}}(\omega)}{\omega}\sim N(0)^{2}$.
This leads to a linear dependence of the nuclear relaxation rate on
temperature, referred to as _Korringa relaxation_ lawKorringa (1950)

$\frac{1}{TT_{1}}\sim\frac{2k_{B}}{\gamma_{e}^{2}\hbar^{4}}N(0)^{2}\sum_{q}|A_{H}(q)|^{2}={\rm
constant}.$ (50)

Figure 11: The NMR relaxation rate $1/TT_{1}$(normalized by the $1/TT_{1}$
value at the transition temperature ($T=T_{t}$) in normal state) as a function
of temperature at superconducting pairing strength $g=0.53$, with varying
value of $U$. The temperature is scaled by the gap-oppening temperature,
$T_{t}=T_{2}$(second-order transition) for $U=0,U=0.2$, and
$T_{t}=T_{p}$(first order transition) for others. To regulate the divergence
of the delta function when doing momentum summation, we introduced a small
imaginary damping rate $i\epsilon$ to the frequency where $\epsilon/W=0.01$.

In a BCS superconductor, we need to take account of the strongly energy-
dependent quasiparticle density of states $N(E)\rightarrow
N(0)\frac{|E|}{\sqrt{E^{2}-\Delta^{2}}}$. If we further assume a point contact
interaction , $A(q)=A$ is a constant with respect to $q$, the relaxation rate
becomes

$\frac{1}{TT_{1}}\propto\int_{\Delta}^{\infty}dE(\frac{df}{dE})\frac{E^{2}}{E^{2}-\Delta^{2}},$
(51)

which generates the Hebel-Slichter peak right below the transition temperature
$T_{c}$, which is shown in Fig. 11 as $U=0$(blue curve).

By contrast, in the HKSC model, there is no fermionic quasipartical excitation
of the ground state. As a result, we shall calculate the susceptibility
$\chi^{s}(q,\omega)$ from its basic definition. The dynamical susceptibility
in imaginary time is

$\begin{split}\chi^{s}_{ab}(q,i\nu_{n})&=\langle{M_{a}(q)M_{b}(-q)}\rangle\\\
&=\int_{0}^{\beta}\langle{M_{a}(q,\tau)M_{b}(-q,0)}\rangle
e^{i\nu_{n}\tau}.\end{split}$ (52)

For a spin-isotropic system, we have
$\chi^{s}_{ab}(q)=\delta_{ab}\chi^{s}(q)$. Thus, we can calculate the $z$-axis
response to a field appllied along $z$:
$\chi^{s}(q)=\langle{M_{z}(q)M_{z}(-q)}\rangle$. Since

$M_{z}(q)=\sum_{{\mathbf{k}}}(c_{{\mathbf{k-q}},\uparrow}^{\dagger}c_{{\mathbf{k}},\uparrow}-c_{{\mathbf{k-q}},\downarrow}^{\dagger}c_{{\mathbf{k}},\downarrow}),$
(53)

we find that

$\begin{split}\chi^{s}(q,\tau)&=\langle{M_{z}(q,\tau)M_{z}(-q,0)}\rangle=\chi_{r}(q,\tau)-\chi_{a}(q,\tau),\end{split}$
(54)

where

$\begin{split}\chi_{r}(q,\tau)&=\sum_{{\mathbf{k}},\sigma}\langle{c_{{\mathbf{k-q}},\sigma}^{\dagger}(\tau)c_{{\mathbf{k}},\sigma}(\tau)c_{{\mathbf{k}},\sigma}^{\dagger}(0)c_{{\mathbf{k-q}},\sigma}(0)}\rangle\\\
\chi_{a}(q,\tau)&=\sum_{{\mathbf{k}},\sigma}\langle{c_{{\mathbf{k-q}},\sigma}^{\dagger}(\tau)c_{{\mathbf{k}},\sigma}(\tau)c_{{\mathbf{-k+q}},\overline{\sigma}}^{\dagger}(0)c_{{\mathbf{-k}},\overline{\sigma}}(0)}\rangle\end{split}$
(55)

are the regular and anomalous parts of the correlation function. The anomalous
term is non-zero only if the ground state no longer conserves the particle
number as in the superconducting phase. The imaginary part of spin
susceptibility takes the form

$\begin{split}\mathrm{Im}\chi(q,\nu-i0^{+})=\sum_{{\mathbf{k}},\sigma}\sum_{i,j,i^{\prime},j^{\prime}}\left(p_{{\mathbf{k}}}^{j}p_{{\mathbf{k-q}}}^{i^{\prime}}-p_{{\mathbf{k}}}^{i}p_{{\mathbf{k-q}}}^{j^{\prime}}\right)\left[\left|M_{{\mathbf{k}},\sigma}^{ij}\right|^{2}\left|M_{{\mathbf{k-q}},\sigma}^{i^{\prime}j^{\prime}}\right|^{2}+M_{{\mathbf{k}},\sigma}^{ij}M_{-{\mathbf{k}},\overline{\sigma}}^{ji}\overline{M}_{{\mathbf{k-q}},\sigma}^{i^{\prime}j^{\prime}}\overline{M}_{-{\mathbf{k+q}},\overline{\sigma}}^{j^{\prime}i^{\prime}}\right]\delta(\nu-(\omega_{{\mathbf{k}}}^{ji}-\omega_{{\mathbf{k-q}}}^{j^{\prime}i^{\prime}})),\end{split}$
(56)

where $(E_{{\mathbf{k}}}^{i},\ket{\psi_{{\mathbf{k}}}^{i}})_{i=1,...,16}$ is
the eigensystem of $H_{{\mathbf{k}}}$, $p_{{\mathbf{k}}}^{i}=e^{-\beta
E_{{\mathbf{k}}}^{i}}/Z_{{\mathbf{k}}}$ is a Boltzmann weight,
$\omega_{{\mathbf{k}}}^{ij}=E_{{\mathbf{k}}}^{i}-E_{{\mathbf{k}}}^{j}$ is an
excitation energy, and
$M_{{\mathbf{k}}\sigma}^{ij}=\langle{\psi_{{\mathbf{k}}}^{i}|c_{{\mathbf{k}}\sigma}|\psi_{{\mathbf{k}}}^{j}}\rangle$.
Thus we can write the NMR relaxation rate as

$\begin{split}\frac{1}{TT_{1}}&\propto\sum_{q}\left.\frac{\imaginary[\chi(q,\nu-i\delta)]}{\nu}\right|_{\nu\rightarrow
0}\\\
&=\sum_{q}\sum_{{\mathbf{k}},\sigma}\sum_{i,j,i^{\prime},j^{\prime}}(p_{{\mathbf{k}}}^{i}+p_{{\mathbf{k}}}^{j})(p_{{\mathbf{k-q}}}^{i^{\prime}}+p_{{\mathbf{k-q}}}^{j^{\prime}})\\\
&\qquad
f^{\prime}(\omega_{{\mathbf{k}}}^{ij})M_{{\mathbf{k}},{\mathbf{q}},\sigma}^{i,j,i^{\prime},j^{\prime}}\delta(\omega_{{\mathbf{k}}}^{ji}-\omega_{{\mathbf{k-q}}}^{j^{\prime}i^{\prime}}),\end{split}$
(57)

where $M_{{\mathbf{k}},{\mathbf{q}},\sigma}^{i,j,i^{\prime},j^{\prime}}$ is
the term inside the square braket from Eq.(56) and $f(x)=1/({e^{\beta x}+1})$
is the Fermi-distribution function. The NMR relaxation rate of the HKSC model
can be calculated numerically. Thus we use the expression given above together
with the self-consistent pair amplitude to obtain the temperature dependence
of $1/TT_{1}$. Fig. 11 showes the temperature dependence of $1/TT_{1}$ at
$g=0.53$. For $U<0.4$, we have a second-order superconducting transtion and as
$U$ increases, the HS peak shrinks and is completely absent when The first
order transtion takes place at $U\geq 0.4$. At sufficiently low temperature,
$1/TT_{1}$ has an exponential dependence on temperature. Once again, we see
that Mottness is the culprit in leading to significant deviations from the
standard BCS theory. Previous workCavanagh and Powell (2021) has attributed
the absenceOhsugi  _et al._ (1994); Yoshimura  _et al._ (1989); Walstedt  _et
al._ (1994) of the HS peak in the cuprates to spin fluctuations. In the
current work, such fluctuations are absent. Hence, we advocate that Mott
physics alone, the tendency to have single occupancy below the chemical
potential, is sufficient to kill the Hebel-Slichter peak.

##  III Ultrasonic attenuation

Alongside the nuclear spin relaxation rate, another standard observable that
probes the superconducting gap is the attenuation rate $\alpha$ of ultrasonic
phonons transmitted through the sample. The attenuation rate,

$\alpha(q)=\Gamma(q,\nu=\omega_{q}),$ (58)

of phonons with momentum $q$ is given by the on-shell decay rate
$\Gamma(q,\nu=\omega_{q})$ of the phonon propagator, where $\omega_{q}\sim
v_{\mathrm{ph}}q$ is the phonon dispersion at small momenta $q$. Typical
experiments are conducted with phonon frequencies on the order of $f\sim
10^{7}\mathrm{Hz}$ Morse and Bohm (1957). Supposing a speed of sound
$v_{\mathrm{ph}}\sim 10^{4}\mathrm{m/s}$ in the sample, these phonons have
energy $E_{\mathrm{ph}}/k_{B}=\hbar 2\pi f/k_{B}\sim 5\times
10^{-4}\mathrm{K}$ and momentum $q=2\pi f/s\sim 6\times
10^{-7}\text{\AA}^{-1}$. Since $E_{\mathrm{ph}}/k_{B}\ll T_{c}$ and $q\ll
k_{F}\sim\text{\AA}^{-1}$, we focus on the limit

$\alpha\equiv\alpha(q\to 0).$ (59)

At weak electron-phonon coupling, the decay rate

$\Gamma(q,\nu)=\imaginary\chi^{(c)}(q,\nu-i0^{+})$ (60)

is given in turn by the electronic charge susceptibility

$\chi^{(c)}(q,i\nu_{n})\equiv\int_{0}^{\beta}d\tau\;e^{i\nu_{n}\tau}\chi^{(c)}(q,\tau)$
(61)

evaluated at the bosonic Matsubara frequency $\nu_{n}=2\pi n/\beta$. In
imaginary time,

$\displaystyle\quad\chi^{(c)}(q,\tau)$
$\displaystyle\equiv\sum_{{\mathbf{k}},k^{\prime},\sigma,\sigma^{\prime}}\langle{c_{{\mathbf{k-q}},\sigma}^{\dagger}(\tau)c_{{\mathbf{k}},\sigma}(\tau)c_{{\mathbf{k}}^{\prime}+q,\sigma^{\prime}}^{\dagger}(0)c_{{\mathbf{k}}^{\prime},\sigma^{\prime}}(0)}\rangle_{0,c}$
(62) $\displaystyle=\chi_{r}(q,\tau)+\chi_{a}(q,\tau),$ (63)

where $\chi_{r}$ and $\chi_{a}$ are the regular and anomalous correlation
functions introduced in the previous section. The sign between the two terms
is reversed relative to the spin susceptibility. This prevents any coherence
peak in the ultrasonic attenuation rate here, just as it does for BCS
superconductors. At small momenta $q$, the attenuation rate takes the form

$\displaystyle\alpha(q)$
$\displaystyle\sim\omega_{q}\sum_{{\mathbf{k}},\sigma}\sum_{i,j,i^{\prime},j^{\prime}}(p_{{\mathbf{k}}}^{i}+p_{{\mathbf{k}}}^{j})(p_{{\mathbf{k}}}^{i^{\prime}}+p_{{\mathbf{k}}}^{j^{\prime}})$
(64) $\displaystyle\quad\times
f^{\prime}(\omega_{{\mathbf{k}}}^{ji})M_{{\mathbf{k}},0,\sigma}^{i,j,i^{\prime},j^{\prime}}\delta\left(v_{\mathrm{ph}}q-(\omega_{{\mathbf{k}}}^{ji}-\omega_{{\mathbf{k-q}}}^{j^{\prime}i^{\prime}})\right),$

analogous to Eq. (57) for the spin relaxation rate. Its convergence in the
limit $q\to 0$ can be seen as follows for the standard context Tsuneto (1961),
i.e. in $d=3$ dimensions and with an isotropic quadratic dispersion
$\xi_{k}=k^{2}/2m-\mu$. The sum is dominated by the excitation between the
ground state and the lowest-lying odd-parity level, for which
$(i,j)=(i^{\prime},j^{\prime})$ and the delta function resolves to
$\delta(v_{\mathrm{ph}}q-{\mathbf{q}}\cdot\gradient_{\mathbf{k}}\omega_{\mathbf{k}}^{ji})$.
Pulling this back to a form that can be formally integrated over $k$, i.e.
$\delta(k-k_{0})/\absolutevalue{\cdots}$, then extracts a factor of $1/q$,
resulting in an overall $q$-dependence given by the product $\omega_{q}/q\sim
v_{\mathrm{ph}}$.

For general parameters the attenuation rate $\alpha$ must be evaluated
numerically, as shown in Fig. 12. In line with the NMR relaxation rate
calculation, we take an isotropic dispersion that is cut off at some magnitude
of the crystal momentum. Unlike the NMR calculation, however, the probe
momentum $q$ is taken asymptotically to zero instead of being summed over, so
it is necessary in this case to perform the integral directly in momentum
space. For the figure, we have used $v_{\mathrm{ph}}/(\mathsf{a}W)=10^{-2}$
and $q=10^{-7}\pi/\mathsf{a}$. The resulting attenuation rate $\alpha$
decreases monotonically in the superconductor from the normal-phase value
$\alpha_{n}(T_{p})$, changing discontinuously at the phase transition when the
gap $\Delta$ opens discontinuously. At low temperatures it decays
exponentially as $e^{-\Delta/T}$. This allows the gap to be extracted from the
ultrasonic attenuation rate in the HKSC. Once again, we see that as the
overall precipitous fall-off of the ultrasonic attenuation rate is accented as
the Mott parameter $U$ increases. In principle, this trend is experimentally
testable.

Figure 12:  Ultrasonic attenuation rate $\alpha$ in the superconducting phase,
at pair coupling $g/W=0.53$ and for representative values of the Mott coupling
$U$. $\alpha$ is normalized by its value $\alpha_{n}(T_{p})$ on the normal
side of the phase transition at $T=T_{p}$. The curves for $U/W=0.6$ and $0.8$
are difficult to see here because they are covered by the $U/W=1$ curve.

##  IV Final Remarks

Table 2: Summary of the superconducting properties in the HK model and a Fermi liquid(FL). $\chi$ represents the pair susceptibility, $\Delta$ is the pairing gap, H. S. stands for the Hebel-SlichterHebel and Slichter (1957) peak, a key feature of BCS superconductors. property | FL | Mottness (HK)
---|---|---
$\chi$ divergent at | $T_{c}$ | $T_{c}(=T_{2})$
$\Delta$ opens at | $T_{c}$ | $T_{p}(>T_{2})$
$\lim_{g\rightarrow 0}2\Delta_{0}/k_{B}T_{c}$ | 3.53 | $\infty$
$E_{\rm cond}/N(0)\Delta^{2}$ | -1 | $[-2,-1]$
Quasiparticles | Bogoliubons | PHYons111linear combinations of Holon and DoublonPhillips  _et al._ (2020)
Ginzburg reduced $t_{G}$ | $\sim 10^{-9}$ | $\sim 10^{-11}$
NMR $1/TT_{1}$ | H-S peak | no H-S peak
Landau Expansion |  $a=\alpha t$ 222$t=\frac{T-T_{c}}{T_{c}}$ and $\alpha>0$ for both cases$,b>0$ | $a=\alpha t,b<0,c>0$
Ultrasonic attenuation, $\alpha$ |  $\sim e^{-\Delta/T}$ below $T_{c}$ |  $\sim e^{-\Delta/T}$ below $T_{p}<T_{c}$

The primary difficulty in unlocking how superconductivity arises in a doped
Mott insulator is computation with a controlled theory. The HK model enables
such an analysis as it represents a fixed point for quartic Fermionic theories
that break the $\mathbb{Z}_{2}$ symmetry of Fermi liquids. Since this includes
the Hubbard model, it suffices to analyse the tractable HK model. An analogy
with Fermi liquids is relevant here. The relevant physics of a Fermi liquid
follows from the free quadratic theory of Eq. (2) as all short-range
repulsions are irrelevant. That such repulsions are irrelevant follows from
the simple fact that Fermi liquids are local in momentum space. Destruction of
this state, except for pairing, requires an equally local interaction in
momentum space. The HK interaction is just the most relevant interaction in
momentum space that suffices as it maximally breaks the $\mathbb{Z}_{2}$
symmetry of a Fermi liquid. Carrying out a pairing analysis from this starting
point should reveal the key differences with how superconductivity obtains
from doped Mott insulator as opposed to a Fermi liquid.

Table (2) catalogs the differences between superconductivity from a FL with
the Mott counterpart. As determined here, the first key dfference is the
appearance of two energy scales, the pairing temperature, $T_{p}$ and the
temperature at which the susceptiblity diverges. In the absence of Mottness,
only a single scale characterizes superconductivity. Nonetheless, in the HK
model we still find that the mean-field theory is essentially exact as the
Ginzburg reduced temperature is vanishingly small. A key prediction here is
that Mottness makes the underlying transition first order. This can be
confirmed by careful mesurements of the latent heat in the cuprates. Another
key prediction here is that the Hebel-Slichter peak is killed by the strong
correlations of the Mott state. While it had been speculated that
antiferromagnetic correlations diminish the relaxation rateCavanagh and Powell
(2021), what we find here is that in a model that has Mott physics but no
antiferromagnetism, the HS peak does not survive. Experimentally, the best NMR
dataOhsugi  _et al._ (1994); Yoshimura  _et al._ (1989); Walstedt  _et al._
(1994) indicate that on the Cu or O sites, no Hebel-Slichter peak exists. Our
work indicates that the suppression of the HS peak is due entirely the
bifurcation of the spectrum into upper and lower Hubbard bands. Such a
bifurcation prevents the coherence that is typically thought to be the
mechanism behind the HS peak. It is these strong correlations of the Mott
state that lead to a deviation as well from the standard Bogoliubov
quasiparticles and the onset of the composite excitations, PYHonsPhillips  _et
al._ (2020), as the new quasi-excitations above the superconducting ground
state. Such correlations also enhance the condensation energy and lead to a
divergence of $\lim_{T\rightarrow 0}2\Delta_{0}/T_{c}$ in the HK
superconductor. As all of these trends are traceable to the strong
correlations of the Mott state, we conclude that Table (2) should provide the
blueprint for superdconductivity in doped Mott insulators.

Acknowledgements

PWP and JZ thank DMR-2111379 for partial funding of this project. E.W.H. was
supported by the Gordon and Betty Moore Foundation EPiQS Initiative through
the grants GBMF 4305 and GBMF 8691.

## Appendix

###  IV.1 stiffness

The superfluid stiffness is defined by

$F[\theta]=F[0]+\int d^{d}x\frac{1}{2}\rho_{s}(\nabla\theta)^{2},$ (65)

where $F=-\frac{1}{\beta}\log Z$ is the free energy and $\theta$ is the phase
of the superfluid. Consider applying the following twist

$c_{i}^{\dagger}\rightarrow c_{i}^{\dagger}e^{i\phi r_{i,x}}.$ (66)

For charge-$n$e superconductivity, $\nabla\theta=n\phi\widehat{x}$ and we can
calculate the stiffness by

$\begin{split}\rho_{s}&=\frac{1}{n^{2}}\frac{1}{N}\left.\frac{\partial^{2}F[\phi]}{\partial\phi^{2}}\right|_{\phi=0}\\\
&=-\frac{1}{n^{2}}\frac{1}{N}\frac{1}{\beta}\left[\frac{\partial_{\phi}^{2}Z}{Z}-\left(\frac{\partial_{\phi}Z}{Z}\right)^{2}\right]_{\phi=0}.\end{split}$
(67)

Here we have set the lattice constant to 1, and $N$ is the total number of
unit cells. In order to derive the correct stiffness for an HK superconductor,
we work with the Hamiltonian

$H=\sum_{{\mathbf{k}}\sigma}\epsilon_{{\mathbf{k}}}c_{{\mathbf{k}}\sigma}^{\dagger}c_{{\mathbf{k}}\sigma}+\sum_{{\mathbf{k}}}Un_{{\mathbf{k}}\uparrow}n_{{\mathbf{k}}\downarrow}+H_{p},$
(68)

where $H_{p}$ is the general superconducting pairing term

$H_{p}=g\sum_{{\mathbf{q}}}\Delta_{{\mathbf{q}}}^{\dagger}\Delta_{{\mathbf{q}}}$
(69)

and
$\Delta_{{\mathbf{q}}}=\sum_{{\mathbf{k}}}c_{{\mathbf{k}}+{\mathbf{q}}\uparrow}c_{-{\mathbf{k}}\downarrow}$
is the copper pair creation operator with momentum $q$. Under the twist we
introduced in Eq.(66), the fourier transform of the fermion operator becomes

$\begin{split}c_{{\mathbf{k}}}^{\dagger}&=\frac{1}{\sqrt{N}}\sum_{i}c_{i}^{\dagger}e^{i{\mathbf{k}}\cdot{\mathbf{r}}_{i}}\\\
&\rightarrow\frac{1}{\sqrt{N}}\sum_{i}c_{i}^{\dagger}e^{i{\mathbf{k}}\cdot{\mathbf{r}}_{i}+i\phi
r_{i,x}}\\\ &=c_{{\mathbf{k}}+\phi\widehat{e}_{x}}^{\dagger}.\end{split}$ (70)

Thus the Hamiltonian under the twist

$\displaystyle\begin{split}H[\phi]&=\sum_{{\mathbf{k}}\sigma}\epsilon_{{\mathbf{k}}}c_{{\mathbf{k}}+\phi\widehat{e}_{x}\sigma}^{\dagger}c_{{\mathbf{k}}+\phi\widehat{e}_{x}\sigma}\\\
&+\sum_{{\mathbf{k}}}Un_{{\mathbf{k}}\uparrow}n_{{\mathbf{k}}\downarrow}+H_{p}[\phi]\end{split}$
(71) $\displaystyle H_{p}[\phi]$
$\displaystyle=g\sum_{{\mathbf{q}}}\Delta_{{\mathbf{q}}}^{\dagger}[\phi]\Delta_{{\mathbf{q}}}[\phi]$
(72) $\displaystyle\Delta_{{\mathbf{q}}}[\phi]$
$\displaystyle=\sum_{{\mathbf{k}}}c_{{\mathbf{k}}+{\mathbf{q}}+\phi\widehat{e}_{x}\uparrow}c_{-{\mathbf{k}}+\phi\widehat{e}_{x}\downarrow}=\Delta_{{\mathbf{q}}+2\phi\widehat{e}_{x}},$
(73)

where the HK term and the superconducting pairing terms do not change under
the twist. Thus,

$H[\phi]-H[0]=\sum_{{\mathbf{k}}\sigma}(\epsilon_{{\mathbf{k}}}-\epsilon_{{\mathbf{k}}-\phi\widehat{e}_{x}})c_{{\mathbf{k}}\sigma}^{\dagger}c_{{\mathbf{k}}\sigma}.$
(74)

We expand the Hamiltonian to second-order in $\phi$ to obtain

$\displaystyle H[\phi]-H[0]$ $\displaystyle=-\phi
J_{x}+\frac{1}{2}\phi^{2}T_{x}$ (75) $\displaystyle J_{x}$
$\displaystyle=\sum_{{\mathbf{k}}\sigma}(\partial_{k_{x}}\epsilon_{{\mathbf{k}}})c_{{\mathbf{k}}\sigma}^{\dagger}c_{{\mathbf{k}}\sigma}$
(76) $\displaystyle T_{x}$
$\displaystyle=\sum_{{\mathbf{k}}\sigma}(\partial^{2}_{k_{x}}\epsilon_{{\mathbf{k}}})c_{{\mathbf{k}}\sigma}^{\dagger}c_{{\mathbf{k}}\sigma},$
(77)

with $J_{x}$ the total current in the $x$ direction and $T_{x}$ is the kinetic
energy arising from hopping in the $x$ direction. The partition function and
its derivatives are

$\displaystyle\begin{split}Z[\phi]&=\operatorname*{tr}{e^{-\beta H[\phi]}}\\\
&=\operatorname*{tr}{e^{-\beta H[0]}Te^{\int_{0}^{\beta}d\tau\phi
J_{x}(\tau)-\frac{1}{2}\phi^{2}T_{x}(\tau)]}}\end{split}$ (78)
$\displaystyle\begin{split}\left.\frac{\partial_{\phi}Z[\phi]}{Z}\right|_{\phi=0}&=\frac{1}{Z}\operatorname*{tr}{e^{-\beta
H[0]}T\left[\int_{0}^{\beta}d\tau J_{x}\right]}\\\
&=\beta\langle{J_{x}}\rangle=0\end{split}$ (79)
$\displaystyle\begin{split}\left.\frac{\partial^{2}_{\phi}Z[\phi]}{Z}\right|_{\phi=0}&=\frac{1}{Z}\operatorname*{tr}{e^{-\beta
H[0]}T\left[\left(\int_{0}^{\beta}d\tau J_{x}\right)^{2}-\int_{0}^{\beta}d\tau
T_{x}\right]}\\\
&=\beta\left(\int_{0}^{\beta}d\tau\langle{J_{x}(\tau)J_{x}}\rangle\right)-\beta\langle{T_{x}}\rangle.\end{split}$
(80)

Finally, the superfluid stiffness is

$\rho_{s}=\frac{1}{n^{2}N}\left(\langle{T_{x}}\rangle-\int_{0}^{\beta}d\tau\langle{J_{x}(\tau)J_{x}}\rangle\right).$
(81)

The decrease of the superfluid stiffness in the HK modelPhillips  _et al._
(2020) was calculated at $T=0$. For finite temperature, especially near
$T_{c}$, we expect that the ratio $\rho_{s}/\Delta^{2}$ remains finite.

Figure 13: Finite temperature superfluid stiffness. For the sake of
illustrating the ginzburg criterion, we keep the multiple solutions to the
self-consistant gap equation even though only one of them represents the true
global minimum of the free energy. The first order transition happens at
$T_{p}$ where the solid lines change into dotted points for $U=0.2,0.5,0.8$.

Fig. 13 shows the temperature dependence of the superfluid stiffness for
multiple value of $U/W$. Near zero temperature the HK pairing reduce the
stiffness by approximately 2 times. The stiffness close to $T_{c}$, however,
was increased drastically up to 10 times. The increased stiffness guarantees
that the Ginzburg reduced temperature is small in the HK model, and thus
proved the applicability of the mean-field theory of HK superconducting model.

## References

  * Bardeen  _et al._ (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 106, 162 (1957).
  * Corson  _et al._ (1999) J. Corson, R. Mallozzi, J. Orenstein, J. N. Eckstein, and I. Bozovic, Nature 398, 221 (1999).
  * Bollinger  _et al._ (2011) A. T. Bollinger, G. Dubuis, J. Yoon, D. Pavuna, J. Misewich, and I. Božović, Nature 472, 458 (2011).
  * Uemura  _et al._ (1989) Y. J. Uemura, G. M. Luke, B. J. Sternlieb, J. H. Brewer, J. F. Carolan, W. N. Hardy, R. Kadono, J. R. Kempton, R. F. Kiefl, S. R. Kreitzman, P. Mulhern, T. M. Riseman, D. L. Williams, B. X. Yang, S. Uchida, H. Takagi, J. Gopalakrishnan, A. W. Sleight, M. A. Subramanian, C. L. Chien, M. Z. Cieplak, G. Xiao, V. Y. Lee, B. W. Statt, C. E. Stronach, W. J. Kossler, and X. H. Yu, Phys. Rev. Lett. 62, 2317 (1989).
  * Emery and Kivelson (1995) V. J. Emery and S. A. Kivelson, Nature 374, 434 (1995).
  * He  _et al._ (2020) Y. He, S.-D. Chen, Z.-X. Li, D. Zhao, D. Song, Y. Yoshida, H. Eisaki, T. Wu, X.-H. Chen, D.-H. Lu, C. Meingast, T. P. Devereaux, R. J. Birgeneau, M. Hashimoto, D.-H. Lee, and Z.-X. Shen, arXiv e-prints , arXiv:2009.10932 (2020), arXiv:2009.10932 [cond-mat.supr-con] .
  * Norman  _et al._ (2005) M. R. Norman, D. Pines, and C. Kallin, Advances in Physics 54, 715 (2005), https://doi.org/10.1080/00018730500459906 .
  * Miller  _et al._ (2017) T. L. Miller, W. Zhang, H. Eisaki, and A. Lanzara, Phys. Rev. Lett. 118, 097001 (2017).
  * Ghosh  _et al._ (2019) S. S. Ghosh, Y. Quan, and W. E. Pickett, Phys. Rev. B 100, 094521 (2019).
  * Anderson and Haldane (2001) P. W. Anderson and F. D. M. Haldane, Journal of Statistical Physics 103, 425 (2001).
  * Jaeger  _et al._ (1989) H. M. Jaeger, D. B. Haviland, B. G. Orr, and A. M. Goldman, Phys. Rev. B 40, 182 (1989).
  * Yazdani and Kapitulnik (1995) A. Yazdani and A. Kapitulnik, Phys. Rev. Lett. 74, 3037 (1995).
  * Molegraaf  _et al._ (2002) H. J. A. Molegraaf, C. Presura, D. van der Marel, P. H. Kes, and M. Li, Science 295, 2239 (2002).
  * Ferrell and Glover (1958) R. A. Ferrell and R. E. Glover, Phys. Rev. 109, 1398 (1958).
  * Tinkham and Ferrell (1959) M. Tinkham and R. A. Ferrell, Phys. Rev. Lett. 2, 331 (1959).
  * Santander-Syro  _et al._ (2003) A. F. Santander-Syro, R. P. S. M. Lobo, N. Bontemps, Z. Konstantinovic, Z. Z. Li, and H. Raffy, Europhysics Letters (EPL) 62, 568 (2003).
  * Rübhausen  _et al._ (2001) M. Rübhausen, A. Gozar, M. V. Klein, P. Guptasarma, and D. G. Hinks, Phys. Rev. B 63, 224514 (2001).
  * Homes  _et al._ (2004) C. C. Homes, S. V. Dordevic, D. A. Bonn, R. Liang, and W. N. Hardy, Phys. Rev. B 69, 024514 (2004).
  * Phillips  _et al._ (2020) P. W. Phillips, L. Yeo, and E. W. Huang, Nature Physics 16, 1175 (2020).
  * Hatsugai and Kohmoto (1992) Y. Hatsugai and M. Kohmoto, Journal of the Physical Society of Japan 61, 2056 (1992), https://doi.org/10.1143/JPSJ.61.2056 .
  * Huang  _et al._ (2022) E. W. Huang, G. L. Nave, and P. W. Phillips, Nature Physics (2022), 10.1038/s41567-022-01529-8.
  * Sellin and Babaev (2016) K. A. H. Sellin and E. Babaev, Phys. Rev. B 93, 054524 (2016).
  * Huang and Lin (2020) W.-M. Huang and H.-H. Lin, Scientific Reports 10, 7439 (2020).
  * Padilha and Continentino (2009) I. T. Padilha and M. Continentino, Physica B: Condensed Matter 404, 2920 (2009).
  * Sarasua and Continentino (2002) L. G. Sarasua and M. A. Continentino, Phys. Rev. B 65, 184503 (2002).
  * Hebel and Slichter (1957) L. C. Hebel and C. P. Slichter, Phys. Rev. 107, 901 (1957).
  * Ohsugi  _et al._ (1994) S. Ohsugi, Y. Kitaoka, K. Ishida, G.-q. Zheng, and K. Asayama, Journal of the Physical Society of Japan 63, 700 (1994).
  * Yoshimura  _et al._ (1989) K. Yoshimura, T. Imai, T. Shimizu, Y. Ueda, K. Kosuge, and H. Yasuoka, Journal of the Physical Society of Japan 58, 3057 (1989).
  * Walstedt  _et al._ (1994) R. Walstedt, B. Shastry, and S. Cheong, Physical review letters 72, 3610 (1994).
  * Zoli (1991) M. Zoli, Journal of the Physical Society of Japan 60, 3837 (1991), https://doi.org/10.1143/JPSJ.60.3837 .
  * Cavanagh and Powell (2021) D. C. Cavanagh and B. J. Powell, Phys. Rev. Research 3, 013241 (2021).
  * Setty (2021) C. Setty, Phys. Rev. B 103, 014501 (2021).
  * Zhu  _et al._ (2021) H.-S. Zhu, Z. Li, Q. Han, and Z. D. Wang, Phys. Rev. B 103, 024514 (2021).
  * Coleman (2015) P. Coleman,  _Introduction to many-body physics_ (Cambridge University Press, 2015).
  * Korringa (1950) J. Korringa, Physica 16, 601 (1950).
  * Morse and Bohm (1957) R. W. Morse and H. V. Bohm, Physical Review 108, 1094 (1957).
  * Tsuneto (1961) T. Tsuneto, Physical Review 121, 402 (1961).