# Pair-density-wave and $s\pm\mathrm{i}d$ superconductivity
in a strongly coupled, lightly doped Kondo insulator

Fangze Liu  Department of Physics, Stanford University, Stanford, California
94305, USA  Stanford Institute for Materials and Energy Sciences, SLAC
National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025,
USA  Zhaoyu Han<EMAIL_ADDRESS>Department of Physics, Stanford
University, Stanford, California 94305, USA

###### Abstract

We investigate the large Kondo coupling limit of the Kondo-Heisenberg model on
one- and two-dimensional lattices. Focusing on the possible superconducting
states when slightly doping the Kondo insulator state, we identify different
pairing modes to be most stable in different parameter regimes. Possibilities
include uniform $s$-wave, pair-density-wave with momentum $\pi$ (in both one
and two dimensions) and uniform $s\pm\mathrm{i}d_{x^{2}-y^{2}}$-wave (in two
dimensions). We attribute these exotic pairing states to the presence of
various pair-hopping terms with a “wrong” sign in the effective model, a
mechanism that is likely universal for inducing pairing states with spatially
modulated pair wavefunctions.

The Kondo-Heisenberg model is a paradigmatic model of strongly correlated
electronic systems, attracting interest as a model of various materials
Coleman (2015); Tsunetsugu  _et al._ (1997); Gulacsi (2006); Zhang  _et al._
(2020); Zhang and Vishwanath (2022); Gall  _et al._ (2021); Sompet  _et al._
(2022); Song and Bernevig (2022); Chou and Das Sarma (2023); Hu  _et al._
(2023); Shankar  _et al._ (2023) and as a prototypical model for various
exotic physical phenomena including quantum criticality Doniach (1977); Si and
Steglich (2010); Paul  _et al._ (2007); Si  _et al._ (2014); Komijani and
Coleman (2019), fractionalization Senthil  _et al._ (2003, 2004), odd
frequency pairing Coleman  _et al._ (1994); Tsvelik (2019) and pair density
wave Zachar (2001); Berg  _et al._ (2010); Chen  _et al._ (2023); Tsvelik
(2019); May-Mann  _et al._ (2020). Although it has received considerable
theoretical investigations, the majority of research efforts have adopted
bosonization techniques in one dimension Tsunetsugu  _et al._ (1997); Gulacsi
(2006); Gulácsi* (2004); Zachar and Tsvelik (2001); Zachar (2001); Emery  _et
al._ (1999) or large-$N$ techniques Read and Newns (1983a); Coleman (1984);
Read and Sachdev (1991) in two or higher dimensions Coleman (2015); Coleman
and Andrei (1989); Wang  _et al._ (2020); Bernhard and Lacroix (2015); Zhu
_et al._ (2008); Iglesias  _et al._ (1997); Flint and Coleman (2010a); Read
and Newns (1983b); Auerbach and Levin (1986); Coleman (1987); Saremi and Lee
(2007); Zhong  _et al._ (2015); Read  _et al._ (1984); Liu  _et al._ (2014);
Flint and Coleman (2010b), and, with a few exceptions Lacroix (1985); Sigrist
_et al._ (1992); Tsunetsugu  _et al._ (1997); Gulacsi (2006); Bastide and
Lacroix (1987); Coleman and Nevidomskyy (2010), focused on relatively weak
coupling regimes (i.e. the Kondo coupling strength is not strong compared to
the bandwidth of the itinerant electrons).

In pursuit of an in-depth understanding of this important model, this work
focuses on the large Kondo coupling regime, where the (antiferromagnetic)
Kondo coupling $J_{\text{K}}$ is much greater than the electron hopping
amplitude $|t|$ and the local moment (antiferromagnetic) Heisenberg coupling
$J_{\text{H}}$. In this limit, a mapping to the infinite $U$ Hubbard model
Lacroix (1985) and a strong coupling expansion can be justified, based on
which we explicitly derive the low-energy effective model. This model is
similar to the $t-J$ model derived from the strong coupling limit of the
Hubbard model but features various additional pair-hopping terms. Based on
this effective model, we investigate the superconducting (SC) phase diagram at
a filling fraction slightly away from one electron per unit cell and at a low
temperature, by means of numerical mean-field (MF) theory. For each set of
parameters ($t/J_{\text{K}}$ and $J_{\text{H}}/J_{\text{K}}$), we solve the
pairing wavefunction self-consistently and compute the free energy at every
Cooper pair momentum $\bm{q}$, allowing us to determine the optimal pairing
momentum, and, if the pairing is at certain high-symmetry momentum, the
pairing symmetry. The phase diagrams are summarized for the 2D and 1D cases in
Figs. 1& 2, respectively. Remarkably, in the 2D scenario, a
$s\pm\mathrm{i}d_{x^{2}-y^{2}}$ pairing state featuring broken time-reversal
symmetry is found to be stable within a regime of the phase diagram. More
interestingly, within a similar parameter regime for both cases, a pair-
density-wave (PDW) Agterberg  _et al._ (2020) with momentum $\pi$ is found to
be the most stable state, in possible agreement with the 1D result obtained by
a bosonization method Zachar (2001) and a numerical density-matrix-
renormalization-group study Berg  _et al._ (2010). We explain the observed
exotic pairing phases from a strong-pairing perspective based on the
observation that the leading pair-hopping terms in the effective model have a
“wrong” sign, which is likely a general mechanism for such exotic pairing
momenta and/or symmetries (see, e.g. Refs. Han  _et al._ (2020); Wu  _et al._
(2023), for similar examples).

Model and Method.  In this paper, we study the Kondo-Heisenberg model:

$\displaystyle\hat{H}=$ $\displaystyle-t\sum_{\langle
ij\rangle,\sigma}(\hat{c}^{\dagger}_{i\sigma}\hat{c}_{j\sigma}+\text{h.c.})+J_{\text{K}}\sum_{i}\hat{\bm{s}}_{i}\cdot\hat{\bm{S}}_{i}$
$\displaystyle\ \ +J_{\text{H}}\sum_{\langle
ij\rangle}\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{j}$ (1)

where $\hat{c}_{i\sigma}$ annihilates a spin-$\sigma$, itinerant electron on
site-$i$, and $\hat{\bm{s}}_{i}$ and $\hat{\bm{S}}_{i}$ respectively represent
the electron spin and the local moment on site $i$. While this model is
definable on any lattice, for concreteness we will focus on the 1D chain and
the 2D square lattice. Due to the bipartite nature of the lattices, there is a
particle-hole symmetry generated by
$\hat{c}_{i\sigma}\rightarrow(-1)^{i}\hat{c}^{\dagger}_{i\sigma}$, allowing us
to concentrate on the hole-doped side, where
$n\equiv\frac{1}{N}\sum_{i}\langle\hat{c}^{\dagger}_{i\sigma}\hat{c}_{i\sigma}\rangle<1$
($N$ is the system size). Furthermore, since the sign of $t$ can be trivially
altered by a gauge transformation
$\hat{c}_{i\sigma}\rightarrow(-1)^{i}\hat{c}_{i\sigma}$, we assume $t>0$
without loss of generality.

In this work, we consider the large $J_{\text{K}}$ limit by regarding
$t/J_{\text{K}}$ and $J_{\text{H}}/J_{\text{K}}$ as small parameters. To the
zeroth order of the analysis and when electron filling $n<1$, each site has
three possible states:
$(|\Uparrow\downarrow\rangle-|\Downarrow\uparrow\rangle)/\sqrt{2}$ (Kondo
singlet) or $|\varnothing\updownarrow\rangle$, where $\Updownarrow$ represents
the spin of the electron whereas $\updownarrow$ represents the local moment,
and $\varnothing$ indicates the absence of any itinerant electron. Different
tensor-product combinations of these states form the low-energy Hilbert space,
$\mathcal{H}_{\text{eff}}$, which is separated from all the other states by an
energy gap $\sim J_{\text{K}}$. To describe the low-energy physics within
$\mathcal{H}_{\text{eff}}$, we define a set of fermionic operators
$\hat{h}_{i\sigma}$ to effectively describe the holes doped into the system.
The mapping between these hole operators and the operators in
$\mathcal{H}_{\text{eff}}$ can be locally established as Lacroix (1985)

$\displaystyle\hat{h}_{i\updownarrow}\leftrightarrow\frac{1}{\sqrt{2}}(|\Uparrow\downarrow\rangle-|\Downarrow\uparrow\rangle)\langle\varnothing\updownarrow|,$
(2)

thus $\hat{h}_{i\sigma}$ annihilates a spin-$\sigma$, charge-$-1$ object
relative to the ‘vacuum’ of $\mathcal{H}_{\text{eff}}$, which refers to the
strong coupling Kondo insulator state with a Kondo singlet on every site. Note
that, to faithfully map between in the physical Hilbert space
$\mathcal{H}_{\text{eff}}$ and the Fock space of the hole operators, it needs
to be further recognized that two holes cannot simultaneously occupy the same
site.

We then perform a perturbation expansion to derive a low-energy effective
Hamiltonian. Leveraging the above equivalence mapping of Hilbert space, we
express the result in terms of the hole operators, including all terms to the
zeroth and the first order in powers of $1/J_{\text{K}}$:

$\displaystyle H_{\text{eff}}=$
$\displaystyle\hat{P}(\hat{H}_{t}+\hat{H}_{\tau}+\hat{H}_{V})\hat{P}$ (3)
$\displaystyle\hat{H}_{t}=$ $\displaystyle t_{1}\sum_{\langle
ij\rangle\sigma}\left(\hat{h}_{i\sigma}^{\dagger}\hat{h}_{j\sigma}+\text{h.c.}\right)$
$\displaystyle+t_{2}\sum_{\langle
ijk\rangle\sigma}\left(\hat{h}_{i\sigma}^{\dagger}\hat{h}_{k\sigma}+\text{h.c.}\right)$
(4) $\displaystyle\hat{H}_{V}=$ $\displaystyle J_{\text{H}}\sum_{\langle
ij\rangle}\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{j}+V\sum_{\langle
ij\rangle}\hat{n}_{i}^{h}\hat{n}_{j}^{h}$
$\displaystyle-J^{\prime}\sum_{\langle
ijk\rangle}\hat{\bm{S}}_{i}\cdot\hat{\bm{S}}_{k}(1-\hat{n}_{j}^{h})$ (5)
$\displaystyle\hat{H}_{\tau}=$ $\displaystyle\sum_{\langle
ijk\rangle}\Bigg{[}t_{1}^{\prime}\sum_{\sigma}\hat{h}^{\dagger}_{i\sigma}\hat{n}^{h}_{k}\hat{h}_{j\sigma}+\tau_{1}\hat{\xi}^{\dagger}_{ik}\hat{\xi}_{kj}+(i\leftrightarrow
k)$ $\displaystyle\ \ \ \ \ \ \
+t_{2}^{\prime}\sum_{\sigma}\hat{h}^{\dagger}_{i\sigma}\hat{n}^{h}_{j}\hat{h}_{k\sigma}+\tau_{2}\hat{\xi}^{\dagger}_{ij}\hat{\xi}_{jk}+\text{h.c.}\Bigg{]}$
(6)

where $\langle ijk\rangle$ represents a triplet of sites in which site $j$ is
a nearest-neighbor of two distinct sites $i$ and $k$,
$\hat{n}^{h}_{i}\equiv\sum_{\sigma}\hat{h}^{\dagger}_{i\sigma}\hat{h}_{i\sigma}$
is the hole density on site-$i$, and $\hat{P}$ is a projector enforcing the
Hilbert constraint, i.e. excluding the states with double occupation of holes
on any site. The effective parameters are
$t_{1}=\frac{t}{2}-\frac{3tJ_{\text{H}}}{4J_{\text{K}}}$,
$t_{2}=\frac{t^{2}}{6J_{\text{K}}}$,
$V=(\frac{5t^{2}}{6J_{\text{K}}}+\frac{9J_{\text{H}}^{2}}{32J_{\text{K}}})$,
$J^{\prime}=\frac{J_{\text{H}}^{2}}{2J_{\text{K}}}$,
$t_{1}^{\prime}=-\frac{tJ_{\text{H}}}{8J_{\text{K}}}$,
$\tau_{1}=\frac{tJ_{\text{H}}}{2J_{\text{K}}}$,
$t_{2}^{\prime}=\frac{t^{2}}{12J_{\text{K}}}$, and
$\tau_{2}=\frac{t^{2}}{2J_{\text{K}}}$. For convenience, we have defined
$\hat{\xi}_{ij}\equiv(\hat{h}_{i\uparrow}\hat{h}_{j\downarrow}+\hat{h}_{j\uparrow}\hat{h}_{i\downarrow})/\sqrt{2}$,
the singlet annihilation operator on sites $i$ and $j$. We note that a similar
mapping and expansion have been done for the Kondo lattice model without
Heisenberg coupling Tsunetsugu  _et al._ (1997); Gulacsi (2006); Bastide and
Lacroix (1987), and our results agree with the existing literature upon
setting $J_{\text{H}}=0$.

To mitigate the complexities of this Hamiltonian, we invoke an exact rewriting
(for any $i,j$)

$\displaystyle\hat{P}\hat{h}^{\dagger}_{i\sigma}\hat{h}_{j\sigma}\hat{P}=(1-\hat{n}^{h}_{i\bar{\sigma}})\hat{h}^{\dagger}_{i\sigma}\hat{h}_{j\sigma}(1-\hat{n}^{h}_{j\bar{\sigma}})$
(7)

to equivalently implement the projection. Then, we consider the dilute hole
limit, i.e. $n^{h}\equiv 1-n\ll 1$. In this limit, the expectation values of
all pairs of fermion operators, i.e. $\hat{h}^{\dagger}\hat{h}$, are bounded
by $n^{h}$. Therefore, we neglect the terms consisting of more than four
fermion operators, as they are of order $\mathcal{O}[(n^{h})^{3}]$ and thus
insignificant relative to other terms.

After these manipulations of the Hamiltonian, writing in momentum space, we
obtain a standard interacting Hamiltonian for fermions:

$\displaystyle H_{\text{eff}}\approx$
$\displaystyle\sum_{\bm{k}\sigma}\epsilon_{\bm{k}}\hat{h}^{\dagger}_{\bm{k}\sigma}\hat{h}_{\bm{k}\sigma}$
$\displaystyle+\frac{1}{N}\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime};\\\
\sigma,\sigma^{\prime};\bm{q}\end{subarray}}\Gamma^{\sigma,\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime};\bm{q}}\hat{h}^{\dagger}_{\frac{\bm{q}}{2}-\bm{k}\sigma}\hat{h}^{\dagger}_{\frac{\bm{q}}{2}+\bm{k}\sigma^{\prime}}\hat{h}_{\frac{\bm{q}}{2}+\bm{k}^{\prime}\sigma^{\prime}}\hat{h}_{\frac{\bm{q}}{2}-\bm{k}^{\prime}\sigma}$
(8)

where the expressions of the bare dispersion $\epsilon_{\bm{k}}$ and
interacting vertex
$\Gamma^{\sigma,\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime};\bm{q}}$ are
explicitly given in Supplemental Materials (SM) 111See [url] for the explicit
expressions of the effective model, the detailed discussions about the mean
field calculation and the results. .

Finally, we perform MF analysis for the possible SC states in the system. Due
to the spin rotation symmetry, we can focus on the sector with
$\sigma^{\prime}=\bar{\sigma}$, since the other two triplet states are
degenerate with the one with $S^{z}=0$. Then for each possible Cooper pair
momentum $\bm{q}$, we adopt the Bogoliubov-de-Gennes MF ansatz (note that in
the equation below, $\bm{q}$ is no longer a dummy variable)

$\displaystyle H^{(\bm{q})}_{\text{MF}}=$
$\displaystyle\sum_{\bm{k}\sigma}\epsilon_{\bm{k}}\hat{h}^{\dagger}_{\bm{k}\sigma}\hat{h}_{\bm{k}\sigma}$
$\displaystyle+\frac{2}{N}\sum_{\bm{k},\bm{k}^{\prime}}\Gamma^{\uparrow\downarrow}_{\bm{k},\bm{k}^{\prime};\bm{q}}\left[\hat{h}^{\dagger}_{\frac{\bm{q}}{2}-\bm{k}\downarrow}\hat{h}^{\dagger}_{\frac{\bm{q}}{2}+\bm{k}\uparrow}\phi_{\bm{k}^{\prime}}^{(\bm{q})}+\text{h.c.}\right]$
(9)

where

$\displaystyle\phi^{(\bm{q})}_{\bm{k}^{\prime}}\equiv\langle\hat{h}_{\frac{\bm{q}}{2}+\bm{k}^{\prime}\uparrow}\hat{h}_{\frac{\bm{q}}{2}+\bm{k}^{\prime}\downarrow}\rangle_{\text{MF}}$
(10)

gives the MF self-consistency equation that we will solve by numerical
iteration.

Figure 1:  Condensation energy in 2D square lattice. (a) Condensation energy
difference (per site) between $\bm{q}=(\pi,0)$ and $\bm{q}=(0,0)$, i.e.
$(F_{\text{H}}^{(\pi,0)}-F_{\text{H}}^{(0,0)})/N$. The black dashed lines
serve as visual guides to demarcate distinct phase regions. (b) Examples of
condensation energies at several high-symmetry $\bm{q}$’s along the vertical
dotted line cut in (a). (c) The full $\bm{q}$ dependence for a few
representative choices of parameters for different phases marked in (a).
Simulations are done with a $N=12\times 12$ lattice for (a) and (b), and a
$N=24\times 24$ lattice for (c).  Figure 2:  Condensation energy in 1D
lattice. (a) Condensation energy difference (per site) between $q=\pi$ and
$q=0$, i.e. $(F_{\text{H}}^{(\pi)}-F_{\text{H}}^{(0)})/N$. The black dashed
lines serve as visual guides to demarcate distinct phase regions. (b) Examples
of condensation energies at several high-symmetry $\bm{q}$’s along the
vertical dotted line cut in (a). (c) The full $\bm{q}$ dependence for a few
representative choices of parameters for different phases marked in (a).
Simulations are done with a $N=128$ chain for (a) and (b), and a $N=256$ chain
for (c).

Numerical results.  For each set of parameters
$\\{t/J_{\text{K}},J_{\text{H}}/J_{\text{K}},n^{h}\\}$, we perform MF
calculation on various different $\bm{q}$ values. For each specific $\bm{q}$,
we solve the MF equation Eq. 10 and obtain a solution with the lowest
Helmholtz free energy, $F^{(\bm{q})}_{\text{H}}$. By comparing the free
energies for different $\bm{q}$ values, the optimal SC states can then be
determined. The condensation energy can be further obtained by comparing the
free energy $F^{(\bm{q})}_{\text{H}}$ with that of a system without SC order,
denoted as $F^{0}_{\text{H}}$.

For all the computations presented in this study, we set $J_{\text{K}}=10$ as
the large energy scale and explore the system’s behavior at $T=1/20$, the
lowest temperature at which we can attain well-converged solutions. For these
calculations, we properly choose the chemical potential to ensure a hole
density of $n^{h}=1/8$. We systematically explore a range of relatively small
values for $J_{\text{H}}\in[0,5]$ and $t\in[0,3]$ are studied. The primary
results are presented in the main text, while more detailed data can be
accessed in the SM Note (1).

We first investigate the 2D square lattice that is of most interest. The
results of the condensation energy density are summarized in Fig. 1. In Fig.
1a&b, it is evident that over a broad region at large $J_{\text{H}}$, a PDW
state with Cooper pair momentum $\bm{q}=(\pi,0)$ (or $(0,\pi)$) is
energetically more favorable, and we have verified that there is no other
competing $\bm{q}$ within the entire Brillouin zone. To provide a better
illustration, we select three representative sets of parameters and plot their
condensation energy density as a function of $\bm{q}$ in Fig. 1c. For
$(t,J_{\text{H}})=(1.3,4.5)$, it is clear that the energy minimum is located
at $(\pi,0)$ (or $(0,\pi)$). It should also be noted that the three curves
have notable qualitative distinctions. Clearly, two of them (with
$(t,J_{\text{H}})=(1.3,4.5)$ and $(0.5,2.0)$) have a small curvature around
the minimum and a narrow bandwidth relative to the condensation energy,
whereas the other point (with $(t,J_{\text{H}})=(1.3,1.0)$) has the opposite
features. This observation suggests that the $J_{\text{H}}\gtrsim 2t$ region
is in a strong pairing regime that can be more suitably described by a Bose-
Einstein Condensation (BEC) of preformed pairs, whereas the
$J_{\text{H}}\lesssim 2t$ region is a weak pairing region and can be
effectively described by the Bardeen-Cooper-Schrieffer (BCS) theory.

Figure 3:  The magnitude of pairing order parameters, as defined in Eqs.
11&12, for the uniform pairing states ($\bm{q}=(0,0)$) in 2D. Simulations are
done on a $12\times 12$ lattice. Note that as in Fig. 1a, distinct phase
regions are delineated by black dashed lines; within the dashed-out region,
the uniform pairing state is less favorable than the PDW state. Figure 4:  The
magnitude (represented by circle radius) and phase (indicated by color) of the
normalized pair field
$\phi_{\bm{r}}^{(\bm{q})}/\sqrt{\sum_{\bm{r}}|\phi_{\bm{r}}^{(\bm{q})}|^{2}}$
(i.e., the wavefunction of the Cooper pairs) as a function of $\bm{r}$ (the
relative coordinate between the two electrons in a cooper pair), at several
representative parameter points marked in Fig. 1a. The
$\phi^{(\bm{q})}_{\bm{r}=(0,0)}$ is manually set to zero in these plots to
respect the Hilbert space constraint. The parameters for these representative
points are as follows: (left) $t=1.3,J_{\text{H}}=1,\bm{q}=(0,0)$; (middle)
$t=0.5,J_{\text{H}}=2,\bm{q}=(0,0)$; (right)
$t=1.3,J_{\text{H}}=4.5,\bm{q}=(\pi,0)$. Simulations are done on a $24\times
24$ 2D square lattice.

To further investigate the pairing symmetries of the uniform pairing
($\bm{q}=(0,0)$) states occupying most parts of the phase diagram in Fig. 1a,
we compute several order parameters defined as:

$\displaystyle
O^{\ell}\equiv\frac{1}{N}\sum_{\bm{k}}f^{\ell}_{\bm{k}}\phi^{(0,0)}_{\bm{k}}$
(11)

where $\ell=s,d_{xy},d_{x^{2}-y^{2}},(p_{x},p_{y})$ are the irreducible
representations of the $D_{4}$ group, and $f^{\ell}$ are the corresponding
form factors. For the pairing symmetries that are non-zero in our case, we
take $f^{s}_{\bm{k}}\equiv\cos{k_{x}}+\cos{k_{y}}$, and
$f^{d_{x^{2}-y^{2}}}_{\bm{k}}\equiv\cos{k_{x}}-\cos{k_{y}}$. To detect time-
reversal symmetry breaking, we compute

$\displaystyle O_{\mathcal{T}}\equiv
1-\frac{|\sum_{\bm{k}}(\phi^{(0,0)}_{\bm{k}})^{2}|}{\sum_{\bm{k}}|\phi^{(0,0)}_{\bm{k}}|^{2}}$
(12)

The amplitudes of these order parameters are plotted in Fig. 3. It is probably
not surprising to see that the BCS uniform pairing state is a pure $s$-wave
state. However, interestingly, we find the BEC uniform phase has coexisting
$s$ and $d_{x^{2}-y^{2}}$ pairing components, and the time-reversal symmetry
is also spontaneously broken, suggesting an exotic
$s\pm\mathrm{i}d_{x^{2}-y^{2}}$ pairing. This finding gains further support
through the direct visualization of the pairing wavefunctions in real space in
Fig. 4, where it can be directly seen that the relative phase between the pair
fields on the nearest neighbor bonds in $x$ and $y$ directions is $\pi/2$. It
is also remarkable that in the dashed-out regime where the uniform pairing
state gives way to the PDW state, the uniform pairing state itself crossovers
from $s+\mathrm{i}d_{x^{2}-y^{2}}$ to a $d_{x^{2}-y^{2}}$-wave state, and has
competitive energy compared to the PDW state (Fig. 1c).

Although mean-field theories are generically less reliable in 1D due to strong
fluctuations, we nonetheless performed the same analysis for the 1D chain
case, with the aim of facilitating comparison with existing results. The
outcomes, as depicted in Fig. 2, closely resemble the findings in the 2D
scenario, and the differences compared to the 2D case are 1) the PDW state
with $\bm{q}=\pi$ is favorable in an even broader regime, and 2) the BEC
uniform pairing state is no longer exotic. It is encouraging to note that, a
density-matrix-renormalization-group study has found PDW to be the ground
state at $J_{\text{H}}/J_{\text{K}}=1$, $t/J_{\text{K}}=1/2$ Berg  _et al._
(2010); May-Mann  _et al._ (2020), a point that is possibly connected to the
PDW regime in Fig. 2a after extrapolating the phase boundary.

A possible mechanism of the exotic SC states.  As seen in Fig. 4, we find that
the interesting PDW and $s\pm\mathrm{i}d_{x^{2}-y^{2}}$ states have dominant
pairing amplitude on the nearest neighbor bonds. On the other hand, from Fig.
1c it can be seen that the energy gain associated with the pair formation
$\sim|F^{(\bm{q})}_{\text{H}}/(n^{h}N)|$ (or the single-particle gap presented
in SM Note (1)), is much higher than the phase stiffness
$\sim\nabla^{2}_{\bm{q}}F^{(\bm{q})}_{\text{H}}/N$ at the optimal $\bm{q}$, so
it can be concluded that these pairs are pre-formed before phase coherence
develops. This can be intuitively understood by the presence of a strong
$J_{\text{H}}$ which stabilizes such a local singlet pairing at a relatively
high energy scale. This observation motivates us to take a perspective
starting from these preformed “bond dimers” by considering an effective dimer
theory (subject to hard-core constraints that are relatively unimportant in
the dilute limit due to the low collision probability):

$\displaystyle\hat{H}_{\mathrm{dimer}}=-\sum_{\langle ij\rangle,\langle
mn\rangle}\left(\tau_{ij,mn}\hat{\xi}^{\dagger}_{ij}\hat{\xi}_{mn}+\textrm{h.c.}\right),$
(13)

where $\tau_{ij,mn}$ is the effective pair hopping amplitude between bond
$\langle ij\rangle$ and bond $\langle mn\rangle$. From the form of the
effective Hamiltonian in Eq. 3, the leading terms that can contribute to the
boson hopping matrix are the $t_{2}$ terms in Eq. 4 and $t_{2}^{\prime}$,
$\tau_{2}$ terms in Eq. 6, which can move a bond dimer to another bond with a
shared site. The crucial thing that allows the exotic pairing states to arise
in this system is that these dominant terms contribute negatively to the
hopping matrix $\tau_{ij,mn}$, circumventing the limitation of the Perron-
Frobenius theorem that always gives rise to a uniform $s$-wave pairing state
and is applicable when all matrix entries are non-negative. Actually, these
leading contributions yield an exactly flat boson band at low energy, which
opens room for the higher-order perturbations in the boson hopping matrix to
lift the degeneracy and lead to an exotic pairing momentum and/or symmetry.
This picture for PDW based on the “wrong” signs of certain pair-hopping terms
seems to be a general, strong coupling mechanism generalizable and applicable
to other systems Han  _et al._ (2020); Wu  _et al._ (2023).

Acknowledgement.  We thank helpful discussions with Steven Kivelson and
Srinivas Raghu. This work is funded by the Department of Energy, Office of
Basic Energy Sciences, Division of Materials Sciences and Engineering, under
contract DE-AC02-76SF00515 at Stanford.

## References

  * Coleman (2015) P. Coleman, arXiv preprint arXiv:1509.05769 (2015).
  * Tsunetsugu  _et al._ (1997) H. Tsunetsugu, M. Sigrist, and K. Ueda, Rev. Mod. Phys. 69, 809 (1997).
  * Gulacsi (2006) M. Gulacsi, Philosophical Magazine 86, 1907 (2006).
  * Zhang  _et al._ (2020) G.-M. Zhang, Y.-f. Yang, and F.-C. Zhang, Phys. Rev. B 101, 020501 (2020).
  * Zhang and Vishwanath (2022) Y.-H. Zhang and A. Vishwanath, Phys. Rev. B 106, 045103 (2022).
  * Gall  _et al._ (2021) M. Gall, N. Wurz, J. Samland, C. F. Chan, and M. Köhl, Nature 589, 40 (2021).
  * Sompet  _et al._ (2022) P. Sompet, S. Hirthe, D. Bourgund, T. Chalopin, J. Bibo, J. Koepsell, P. Bojović, R. Verresen, F. Pollmann, G. Salomon, _et al._ , Nature 606, 484 (2022).
  * Song and Bernevig (2022) Z.-D. Song and B. A. Bernevig, Phys. Rev. Lett. 129, 047601 (2022).
  * Chou and Das Sarma (2023) Y.-Z. Chou and S. Das Sarma, Phys. Rev. Lett. 131, 026501 (2023).
  * Hu  _et al._ (2023) H. Hu, B. A. Bernevig, and A. M. Tsvelik, Phys. Rev. Lett. 131, 026502 (2023).
  * Shankar  _et al._ (2023) A. S. Shankar, D. O. Oriekhov, A. K. Mitchell, and L. Fritz, Phys. Rev. B 107, 245102 (2023).
  * Doniach (1977) S. Doniach, Physica B+C 91, 231 (1977).
  * Si and Steglich (2010) Q. Si and F. Steglich, Science 329, 1161 (2010), https://www.science.org/doi/pdf/10.1126/science.1191195 .
  * Paul  _et al._ (2007) I. Paul, C. Pépin, and M. R. Norman, Phys. Rev. Lett. 98, 026402 (2007).
  * Si  _et al._ (2014) Q. Si, J. H. Pixley, E. Nica, S. J. Yamamoto, P. Goswami, R. Yu, and S. Kirchner, Journal of the Physical Society of Japan 83, 061005 (2014).
  * Komijani and Coleman (2019) Y. Komijani and P. Coleman, Phys. Rev. Lett. 122, 217001 (2019).
  * Senthil  _et al._ (2003) T. Senthil, S. Sachdev, and M. Vojta, Phys. Rev. Lett. 90, 216403 (2003).
  * Senthil  _et al._ (2004) T. Senthil, M. Vojta, and S. Sachdev, Phys. Rev. B 69, 035111 (2004).
  * Coleman  _et al._ (1994) P. Coleman, E. Miranda, and A. Tsvelik, Phys. Rev. B 49, 8955 (1994).
  * Tsvelik (2019) A. M. Tsvelik, Proceedings of the National Academy of Sciences 116, 12729 (2019), https://www.pnas.org/doi/pdf/10.1073/pnas.1902928116 .
  * Zachar (2001) O. Zachar, Phys. Rev. B 63, 205104 (2001).
  * Berg  _et al._ (2010) E. Berg, E. Fradkin, and S. A. Kivelson, Phys. Rev. Lett. 105, 146403 (2010).
  * Chen  _et al._ (2023) J. Chen, J. Wang, and Y.-f. Yang, arXiv preprint arXiv:2308.11414 (2023).
  * May-Mann  _et al._ (2020) J. May-Mann, R. Levy, R. Soto-Garrido, G. Y. Cho, B. K. Clark, and E. Fradkin, Phys. Rev. B 101, 165133 (2020).
  * Gulácsi* (2004) M. Gulácsi*, Advances in physics 53, 769 (2004).
  * Zachar and Tsvelik (2001) O. Zachar and A. M. Tsvelik, Phys. Rev. B 64, 033103 (2001).
  * Emery  _et al._ (1999) V. J. Emery, S. A. Kivelson, and O. Zachar, Phys. Rev. B 59, 15641 (1999).
  * Read and Newns (1983a) N. Read and D. Newns, Journal of Physics C: Solid State Physics 16, 3273 (1983a).
  * Coleman (1984) P. Coleman, Phys. Rev. B 29, 3035 (1984).
  * Read and Sachdev (1991) N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991).
  * Coleman and Andrei (1989) P. Coleman and N. Andrei, Journal of Physics: Condensed Matter 1, 4057 (1989).
  * Wang  _et al._ (2020) J. Wang, Y.-Y. Chang, C.-Y. Mou, S. Kirchner, and C.-H. Chung, Phys. Rev. B 102, 115133 (2020).
  * Bernhard and Lacroix (2015) B. H. Bernhard and C. Lacroix, Phys. Rev. B 92, 094401 (2015).
  * Zhu  _et al._ (2008) J.-X. Zhu, I. Martin, and A. R. Bishop, Phys. Rev. Lett. 100, 236403 (2008).
  * Iglesias  _et al._ (1997) J. R. Iglesias, C. Lacroix, and B. Coqblin, Phys. Rev. B 56, 11820 (1997).
  * Flint and Coleman (2010a) R. Flint and P. Coleman, Phys. Rev. Lett. 105, 246404 (2010a).
  * Read and Newns (1983b) N. Read and D. Newns, Journal of Physics C: Solid State Physics 16, L1055 (1983b).
  * Auerbach and Levin (1986) A. Auerbach and K. Levin, Phys. Rev. Lett. 57, 877 (1986).
  * Coleman (1987) P. Coleman, Phys. Rev. B 35, 5072 (1987).
  * Saremi and Lee (2007) S. Saremi and P. A. Lee, Phys. Rev. B 75, 165110 (2007).
  * Zhong  _et al._ (2015) Y. Zhong, L. Zhang, H.-T. Lu, and H.-G. Luo, The European Physical Journal B 88, 238 (2015).
  * Read  _et al._ (1984) N. Read, D. M. Newns, and S. Doniach, Phys. Rev. B 30, 3841 (1984).
  * Liu  _et al._ (2014) Y. Liu, G.-M. Zhang, and L. Yu, Chinese Physics Letters 31, 087102 (2014).
  * Flint and Coleman (2010b) R. Flint and P. Coleman, Physical review letters 105, 246404 (2010b).
  * Lacroix (1985) C. Lacroix, Solid state communications 54, 991 (1985).
  * Sigrist  _et al._ (1992) M. Sigrist, H. Tsunetsugu, K. Ueda, and T. M. Rice, Phys. Rev. B 46, 13838 (1992).
  * Bastide and Lacroix (1987) C. Bastide and C. Lacroix, Europhysics Letters 4, 935 (1987).
  * Coleman and Nevidomskyy (2010) P. Coleman and A. H. Nevidomskyy, Journal of Low Temperature Physics 161, 182 (2010).
  * Agterberg  _et al._ (2020) D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * Han  _et al._ (2020) Z. Han, S. A. Kivelson, and H. Yao, Phys. Rev. Lett. 125, 167001 (2020).
  * Wu  _et al._ (2023) Y.-M. Wu, P. A. Nosov, A. A. Patel, and S. Raghu, Phys. Rev. Lett. 130, 026001 (2023).
  * Note (1) See [url] for the explicit expressions of the effective model, the detailed discussions about the mean field calculation and the results.